From sources-request at octave dot org Mon Aug 22 13:03:53 2005
Subject: ppval
From: "Ray Muzic" <muzic at uhrad dot com>
To: octave sources mailing list <sources at octave dot org>
Date: Fri, 19 Aug 2005 16:11:21 -0400

[Perhaps someone would like to adapt this and replace ppval.m in
octave-forge?  --jwe]


John,

I see your wish list for octave includes ppval.
Below is a good start.  It is a MATLAB mex source I wrote to
implement ppval within MATLAB.  It is much faster than the 
non-compiled ppval.m MATLAB version.

Perhaps you can adapt my c-code to work with Octave.

You may use it per GPL.

/* ppval.c: mex implementation of pieceise polynomial evaluation */
/* to build/compile: mex ppval.c                                 */
/* Copyright (c) 2002 Raymond F. Muzic, Jr.                        */
/* ray muzic, 20020410 original version                          */

#include "mex.h"
#include "matrix.h"

#include <sys/timeb.h>
#include <time.h>

//static double timeAll = 0;
//static double timeFind = 0;

/* helper function to findidx.  bisect searches using bisection */
int bisect(double v, int l, int u, double *a)
{
   int m;
   if ((u - l) <= 1) return l;
   m = (l+u)/2;
   if (v >= a[m]) return bisect(v, m, u, a);
   return bisect(v, l, m, a);
}


/* findidx reutrns index i of a such a[i] <= v < a[i+1] 
   with i constrained l <= i <= u
 */
int findidx(double v, int l, int u, double *a)
{
   if (a[l] >= v) return l;
   if (a[u] < v) return u;
   return bisect(v, l, u, a);
}


void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
   double lx;  /* local x coordinate                         */
   
   double *pc;     /* polynomial coefficients                    */
   int pieces;     /* number of intervals = number of breaks - 1 */
   int dim;        /* dimension                                  */
   int order;      /* polynomial order + 1                       */
   double *breaks; /* breaks                                     */
   double *x;      /* independent variable                       */
   char form[3];
   int i, j, k, l, n;
   double d;
   double *out;
   int numInD;
   int numOutD;
   int outDim[100];  /* assume number of dimensions <= 100 */
   int *inDim;

//   struct _timeb timeStart, timeFinish;   /* for determination of run
time */
//   struct _timeb timeFindStart, timeFindFinish;   /* for determination
of run time */
   
//   _ftime( &timeStart );

   if (nrhs == 0) { mexPrintf("ppval (mex): mex version of ppval.  See
ppval.m for usage.\nppval (mex) was created by Ray Muzic
rfm2 at po dot cwru dot edu dot   User assumes all risk.\n"); return; }
   if (nrhs != 2) mexErrMsgTxt("ppval (mex): Two inputs required:
ppval(pp, x)");
   if (mxGetNumberOfElements(mxGetField(prhs[0], 0, "form")) != 2)
mexErrMsgTxt("ppval (mex): First argument must be 'pp' form");
   mxGetString(mxGetField(prhs[0], 0, "form"), form, 3);
   if (strcmp(form, "pp")) mexErrMsgTxt("ppval (mex): First argument
must be 'pp' form");


   /* pp struct */
   breaks = mxGetPr(mxGetField(prhs[0], 0, "breaks"));
   pc = mxGetPr(mxGetField(prhs[0], 0, "coefs"));
   pieces = (int )(*mxGetPr(mxGetField(prhs[0], 0, "pieces")));
   order = (int )(*mxGetPr(mxGetField(prhs[0], 0, "order")));
   dim = (int )(*mxGetPr(mxGetField(prhs[0], 0, "dim")));

   /* independent variable */ 
   x = mxGetPr(prhs[1]);
   
   /* create output variable */
   /*  be a little fancier than line below in order to support  */
   /*  dim > 1 and prhs[1] as multi-dimensional array           */
   /*  plhs[0] = mxCreateDoubleMatrix(dim, n, 1, mxREAL);       */
   numInD = mxGetNumberOfDimensions(prhs[1]);
   inDim = mxGetDimensions(prhs[1]);
   
   /* do not change inDim!  It is NOT a copy of dimension info, it is a
reference to it */
   for (i = 0; i < numInD; i++) outDim[i] = inDim[i];
   outDim[0] *= dim;
   numOutD = numInD;

   plhs[0] = mxCreateNumericArray(numOutD, outDim, mxDOUBLE_CLASS,
mxREAL);
   n = mxGetNumberOfElements(prhs[1]);
   out = mxGetPr(plhs[0]);

   for (i = 0; i < n; i++) {
//      _ftime( &timeFindStart );
      j = findidx(x[i], 0, pieces-1, breaks); /* breaks[j] <= x[i] <
breaks[j+1] */
//      _ftime( &timeFindFinish );
      lx = x[i] - breaks[j];     /* local coordinate */
      for (k = 0; k < dim; k++) {
         for (l = 1, d = pc[dim*j+k]; l < order; l++) d = d * lx +
pc[dim*j+k+l*(dim*pieces)];
         *out++ = d;
      }
   }

//    _ftime( &timeFinish );


//    timeAll += 	     ((double)timeFinish.time +
(double)timeFinish.millitm/1000.) - 
//	     ((double)timeStart.time + (double)timeStart.millitm/1000.);
//    timeFind +=        ((double)timeFindFinish.time +
(double)timeFindFinish.millitm/1000.) - 
//	     ((double)timeFindStart.time +
(double)timeFindStart.millitm/1000.);
//    mexPrintf("whole: %f,   find: %f\n", timeAll, timeFind);

}


Raymond F. Muzic, Jr., Ph.D.
Associate Professor, Radiology, Oncology & Biomedical Engineering
Case Western Reserve University

Email                           Snail Mail
muzic at ieee dot org                  Nuclear Medicine
muzic at uhrad dot com                 Univ. Hospitals Cleveland 
raymond dot muzic at case dot edu          11100 Euclid Avenue 
Voice: 216/844-3543, Lab: 7350  Cleveland, OH 44106
FAX:   216/844-1048

You got to be careful if you don't know where you're going, 
because you might not get there. - Yogi Berra
  

I have learned to use the word "impossible" with the 
greatest caution. - Wernher von Braun  

Research is what I'm doing when I don't know what I'm doing. - Wernher
von Braun 



 
CONFIDENTIALITY NOTICE: This e-mail message, including any attachments,
is for the sole use of the intended recipient(s) and may contain
confidential and privileged information. Any unauthorized review, use,
disclosure or distribution is prohibited. If you are not the intended
recipient, please contact the sender by reply e-mail and destroy all
copies of the original message.