Index: /trunk/psLib/src/dataManip/psFunctions.c =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.c (revision 1858) +++ /trunk/psLib/src/dataManip/psFunctions.c (revision 1859) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-22 01:30:21 $ + * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-23 04:56:40 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -467,4 +467,6 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { if (myPoly->mask[loop_x] == 0) { + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10, + "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]); polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -1687,5 +1689,5 @@ /***************************************************************************** -VectorBinDisectF32(): This is a private function which takes as input a +p_psVectorBinDisectF32(): This is a private function which takes as input a vector of floating point data as well as a single floating point values. The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for @@ -1696,7 +1698,7 @@ XXX: Macro this for a few different types. *****************************************************************************/ -int VectorBinDisectF32(float *bins, - int numBins, - float x) +int p_psVectorBinDisectF32(float *bins, + int numBins, + float x) { int min; @@ -1704,11 +1706,11 @@ int mid; - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Calling VectorBinDisectF32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Calling p_psVectorBinDisectF32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", + "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1720,11 +1722,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); if (x == bins[mid]) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Exiting VectorBinDisectF32(): bin %d\n", mid); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid); return(mid); } else if (x < bins[mid]) { @@ -1736,15 +1738,15 @@ } - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Exiting VectorBinDisectF32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min); return(min); } /***************************************************************************** -VectorBinDisectS32(): integer version of above. +p_psVectorBinDisectS32(): integer version of above. *****************************************************************************/ -int VectorBinDisectS32(int *bins, - int numBins, - int x) +int p_psVectorBinDisectS32(int *bins, + int numBins, + int x) { int min; @@ -1752,11 +1754,11 @@ int mid; - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Calling VectorBinDisectS32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Calling p_psVectorBinDisectS32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", + "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1768,11 +1770,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); if (x == bins[mid]) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Exiting VectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); return(min); } else if (x < bins[mid]) { @@ -1784,11 +1786,11 @@ } - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Exiting VectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); return(min); } /***************************************************************************** -p_psVectorBinDisect(): A wrapper to the above VectorBinDisect(). +p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). *****************************************************************************/ int p_psVectorBinDisect(psVector *bins, @@ -1801,7 +1803,7 @@ if (x->type.type == PS_TYPE_S32) { - return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); + return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); } else if (x->type.type == PS_TYPE_F32) { - return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); + return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); } else { // XXX: Generate error message. @@ -1879,5 +1881,5 @@ "---- p_ps1DInterpolateF32() begin ----\n"); - binNum = VectorBinDisectF32(domain, n, x); + binNum = p_psVectorBinDisectF32(domain, n, x); if (0 == numIntPoints%2) { @@ -1953,5 +1955,5 @@ n = spline->n; - binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x); + binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x); if (binNum == -1) { psLogMsg(__func__, PS_LOG_WARN, Index: /trunk/psLib/src/dataManip/psFunctions.h =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.h (revision 1858) +++ /trunk/psLib/src/dataManip/psFunctions.h (revision 1859) @@ -12,6 +12,6 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-09-22 01:30:21 $ +* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-09-23 04:56:41 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -407,4 +407,12 @@ psVector *data); +int p_psVectorBinDisectF32(float *bins, + int numBins, + float x); + +int p_psVectorBinDisectS32(int *bins, + int numBins, + int x); + int p_psVectorBinDisect(psVector *bins, psScalar *x); @@ -415,4 +423,9 @@ psScalar *x); +float p_psNRSpline1DEval(psSpline1D *spline, + const psVector* restrict x, + const psVector* restrict y, + float X); + /* \} */// End of MathGroup Functions Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 1858) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 1859) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-22 02:42:49 $ + * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-23 04:56:41 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -43,12 +43,9 @@ #include "psMinimize.h" #include "psMatrix.h" -#include "psTrace.h" +#include "psConstants.h" /*****************************************************************************/ /* DEFINE STATEMENTS */ /*****************************************************************************/ - -#define MAX_LMM_ITERATIONS 100 -#define MAX_MINIMIZE_ITERATIONS 100 /** Preprocessor macro to generate error on a NULL 1DPolynomial */ @@ -96,5 +93,4 @@ /* GLOBAL VARIABLES */ /*****************************************************************************/ - static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL; //static psMinimizePowellFunc PowellFunc = NULL; @@ -116,6 +112,10 @@ input parameter "x" between 0 and input parameter polyOrder. The result is returned as a psVector sums. + +XXX: change name *****************************************************************************/ -void p_psBuildSums1D(double x, int polyOrder, psVector* sums) +void p_psBuildSums1D(double x, + int polyOrder, + psVector* sums) { int i = 0; @@ -148,13 +148,12 @@ /***************************************************************************** - CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a -tabulated function at n points, this routine calcualtes the second +tabulated function at n points, this routine calculates the second derivatives of the interpolating cubic splines at those n points. The first and second derivatives at the endpoints, undefined in the SDR, are -here defined to be 0.0. - -XXX: This algorithm is very similar to that in Numerical Recipes. +here defined to be 0.0. They can be modified via ypo and yp1. + +XXX: This algorithm is derived from the Numerical Recipes. *****************************************************************************/ float *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) @@ -174,4 +173,5 @@ float *X = (float *) & (x->data.F32[0]); float *Y = (float *) & (y->data.F32[0]); + float qn; if (x == NULL) { @@ -183,10 +183,8 @@ } - // XXX: The first and second derivatives at the endpoints, undefined in - // the SDR, are here defined to be 0.0. - u[0]= 0.0; - u[n-1]= 0.0; - derivs2[0] = 0.0; - derivs2[n-1] = 0.0; + // XXX: The second derivatives at the endpoints, undefined in the SDR, + // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV. + derivs2[0] = -0.5; + u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - LEFT_SPLINE_DERIV); for (i=1;i<=(n-2);i++) { @@ -194,6 +192,6 @@ p = sig * derivs2[i-1] + 2.0; derivs2[i] = (sig - 1.0) / p; - u[i] = (Y[i+1] - Y[i])/(X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]); - u[i] = (6.0 * u[i] / (X[i+1] - X[i-1]) - sig * u[i-1]) / p; + u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1])); + u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p; psTrace(".psLib.dataManip.CalculateSecondDerivs", 6, @@ -205,9 +203,15 @@ } + qn = 0.5; + u[n-1] = (3.0/(X[n-1]-X[n-2])) * (RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2])); + derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0); + for (k=(n-2);k>=0;k--) { derivs2[k] = derivs2[k] * derivs2[k+1] + u[k]; + psTrace(".psLib.dataManip.CalculateSecondDerivs", 6, "derivs2[%d] is %f\n", k, derivs2[k]); } + if (mustFreeX == true) { psFree(X); @@ -224,6 +228,4 @@ /* FUNCTION IMPLEMENTATION - PUBLIC */ /*****************************************************************************/ - - /***************************************************************************** @@ -238,8 +240,11 @@ (3) ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 + (4) ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0 -Where - H = x[1]-x[0], A = (x[1]-x)/H, B = (x-x[0])/H +Where: + H = x[1]-x[0] + A = (x[1]-x)/H + B = (x-x[0])/H The bulk of the code in this routine is the expansion of the above equation -into a polynomial in terms of x. This gets pretty complicated. +into a polynomial in terms of x, and then saving the coefficients of the +powers of x in the spline polynomials. This gets pretty complicated. XXX: usage of yErr is not specified in IfA documentation. @@ -262,4 +267,7 @@ } mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y); + for (i=0;in;i++) + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); for (i=0;ispline[i])->coeff[0] = y->data.F32[i] * (x->data.F32[i+1]/H); + (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H); // From (2) - ((mySpline->spline[i])->coeff[0])-= (y->data.F32[i+1] * x->data.F32[i]/H); + ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H); // From (3) tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); @@ -289,17 +297,17 @@ // // From (1) - (mySpline->spline[i])->coeff[1] = - ((mySpline->p_psDeriv2)[i]) / H; + (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H; // From (2) - (mySpline->spline[i])->coeff[1]-= ((mySpline->p_psDeriv2)[i+1]) / H; + ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H); // From (3) - tmp = - (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); + tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); tmp+= (1.0 / H); tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0; - (mySpline->spline[i])->coeff[1]+= tmp; + ((mySpline->spline[i])->coeff[1])+= tmp; // From (4) - tmp = (x->data.F32[i] + x->data.F32[i]) / (H * H * H); + tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H); tmp-= (1.0 / H); tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0; - (mySpline->spline[i])->coeff[1]+= tmp; + ((mySpline->spline[i])->coeff[1])+= tmp; // @@ -307,7 +315,7 @@ // // From (3) - (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / 6.0; + (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H); // From (4) - (mySpline->spline[i])->coeff[2]-= ((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / 6.0; + ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H)); // @@ -317,5 +325,5 @@ (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H); // From (4) - (mySpline->spline[i])->coeff[3]+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); + ((mySpline->spline[i])->coeff[3])+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, @@ -334,4 +342,48 @@ return(mySpline); } + +/****************************************************************************** +p_psNRSpline1DEval(): This routine does NR-style evaluation of cubic splines. +It takes advantage of the 2nd derivatives of the cubic splines, which are +stored in the psSPline1D data structure, and computes the interpolated value +directly, without computing (or using) the interpolating cubic spline +polynomial. + +This routine is here mostly for a sanity check on the psLib function +evalSpline() which computes the interpolated value based on the cubic spline +polynomials which are stored in psSpline1D. + *****************************************************************************/ +float p_psNRSpline1DEval(psSpline1D *spline, + const psVector* restrict x, + const psVector* restrict y, + float X) +{ + int n; + int klo; + int khi; + float H; + float A; + float B; + float C; + float D; + float Y; + + n = spline->n; + klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X); + khi = klo + 1; + H = (spline->domains)[khi] - (spline->domains)[klo]; + A = ((spline->domains)[khi] - X) / H; + B = (X - (spline->domains)[klo]) / H; + C = ((A*A*A)-A) * (H*H/6.0); + D = ((B*B*B)-B) * (H*H/6.0); + + Y = (A * y->data.F32[klo]) + + (B * y->data.F32[khi]) + + (C * (spline->p_psDeriv2)[klo]) + + (D * (spline->p_psDeriv2)[khi]); + + return(Y); +} + /****************************************************************************** @@ -687,10 +739,9 @@ } - -#define STEP_SIZE 0.10 /****************************************************************************** - This routine takes as input an arbitrary function, and the parameter to - vary, and the line along which it must vary. This function produces as - output a bracket [a, b, c] such that +p_psDetermineBracket(): This routine takes as input an arbitrary function, +and the parameter to vary, and the line along which it must vary. This +function produces as output a bracket [a, b, c] such that + f(param + b * line) is less than f(param + a * line) and f(param + c * line). @@ -700,4 +751,6 @@ smaller/larger than b. Repeat this process until a local minimum is found. + + *****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, @@ -719,5 +772,5 @@ float new_cDir = 0.0; psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); - float stepSize = STEP_SIZE; + float stepSize = DETERMINE_BRACKET_STEP_SIZE; psVector *tmp = NULL; int i = 0; @@ -869,15 +922,4 @@ return(NULL); } - - - - - - - - - - - /****************************************************************************** @@ -1061,17 +1103,4 @@ - - - - - - - - - - - - - /****************************************************************************** This routine must minimize a possibly multi-dimensional function. The Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 1858) +++ /trunk/psLib/src/math/psMinimize.c (revision 1859) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-22 02:42:49 $ + * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-23 04:56:41 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -43,12 +43,9 @@ #include "psMinimize.h" #include "psMatrix.h" -#include "psTrace.h" +#include "psConstants.h" /*****************************************************************************/ /* DEFINE STATEMENTS */ /*****************************************************************************/ - -#define MAX_LMM_ITERATIONS 100 -#define MAX_MINIMIZE_ITERATIONS 100 /** Preprocessor macro to generate error on a NULL 1DPolynomial */ @@ -96,5 +93,4 @@ /* GLOBAL VARIABLES */ /*****************************************************************************/ - static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL; //static psMinimizePowellFunc PowellFunc = NULL; @@ -116,6 +112,10 @@ input parameter "x" between 0 and input parameter polyOrder. The result is returned as a psVector sums. + +XXX: change name *****************************************************************************/ -void p_psBuildSums1D(double x, int polyOrder, psVector* sums) +void p_psBuildSums1D(double x, + int polyOrder, + psVector* sums) { int i = 0; @@ -148,13 +148,12 @@ /***************************************************************************** - CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a -tabulated function at n points, this routine calcualtes the second +tabulated function at n points, this routine calculates the second derivatives of the interpolating cubic splines at those n points. The first and second derivatives at the endpoints, undefined in the SDR, are -here defined to be 0.0. - -XXX: This algorithm is very similar to that in Numerical Recipes. +here defined to be 0.0. They can be modified via ypo and yp1. + +XXX: This algorithm is derived from the Numerical Recipes. *****************************************************************************/ float *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) @@ -174,4 +173,5 @@ float *X = (float *) & (x->data.F32[0]); float *Y = (float *) & (y->data.F32[0]); + float qn; if (x == NULL) { @@ -183,10 +183,8 @@ } - // XXX: The first and second derivatives at the endpoints, undefined in - // the SDR, are here defined to be 0.0. - u[0]= 0.0; - u[n-1]= 0.0; - derivs2[0] = 0.0; - derivs2[n-1] = 0.0; + // XXX: The second derivatives at the endpoints, undefined in the SDR, + // are set in psConstants.h: LEFT_SPLINE_DERIV, RIGHT_SPLINE_DERIV. + derivs2[0] = -0.5; + u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - LEFT_SPLINE_DERIV); for (i=1;i<=(n-2);i++) { @@ -194,6 +192,6 @@ p = sig * derivs2[i-1] + 2.0; derivs2[i] = (sig - 1.0) / p; - u[i] = (Y[i+1] - Y[i])/(X[i+1]-X[i]) - (Y[i]-Y[i-1])/(X[i]-X[i-1]); - u[i] = (6.0 * u[i] / (X[i+1] - X[i-1]) - sig * u[i-1]) / p; + u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1])); + u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p; psTrace(".psLib.dataManip.CalculateSecondDerivs", 6, @@ -205,9 +203,15 @@ } + qn = 0.5; + u[n-1] = (3.0/(X[n-1]-X[n-2])) * (RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2])); + derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0); + for (k=(n-2);k>=0;k--) { derivs2[k] = derivs2[k] * derivs2[k+1] + u[k]; + psTrace(".psLib.dataManip.CalculateSecondDerivs", 6, "derivs2[%d] is %f\n", k, derivs2[k]); } + if (mustFreeX == true) { psFree(X); @@ -224,6 +228,4 @@ /* FUNCTION IMPLEMENTATION - PUBLIC */ /*****************************************************************************/ - - /***************************************************************************** @@ -238,8 +240,11 @@ (3) ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 + (4) ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0 -Where - H = x[1]-x[0], A = (x[1]-x)/H, B = (x-x[0])/H +Where: + H = x[1]-x[0] + A = (x[1]-x)/H + B = (x-x[0])/H The bulk of the code in this routine is the expansion of the above equation -into a polynomial in terms of x. This gets pretty complicated. +into a polynomial in terms of x, and then saving the coefficients of the +powers of x in the spline polynomials. This gets pretty complicated. XXX: usage of yErr is not specified in IfA documentation. @@ -262,4 +267,7 @@ } mySpline->p_psDeriv2 = CalculateSecondDerivs(x, y); + for (i=0;in;i++) + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); for (i=0;ispline[i])->coeff[0] = y->data.F32[i] * (x->data.F32[i+1]/H); + (mySpline->spline[i])->coeff[0] = (y->data.F32[i] * x->data.F32[i+1]/H); // From (2) - ((mySpline->spline[i])->coeff[0])-= (y->data.F32[i+1] * x->data.F32[i]/H); + ((mySpline->spline[i])->coeff[0])-= ((y->data.F32[i+1] * x->data.F32[i])/H); // From (3) tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); @@ -289,17 +297,17 @@ // // From (1) - (mySpline->spline[i])->coeff[1] = - ((mySpline->p_psDeriv2)[i]) / H; + (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H; // From (2) - (mySpline->spline[i])->coeff[1]-= ((mySpline->p_psDeriv2)[i+1]) / H; + ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H); // From (3) - tmp = - (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); + tmp = -3.0 * (x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); tmp+= (1.0 / H); tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0; - (mySpline->spline[i])->coeff[1]+= tmp; + ((mySpline->spline[i])->coeff[1])+= tmp; // From (4) - tmp = (x->data.F32[i] + x->data.F32[i]) / (H * H * H); + tmp = 3.0 * (x->data.F32[i] * x->data.F32[i]) / (H * H * H); tmp-= (1.0 / H); tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0; - (mySpline->spline[i])->coeff[1]+= tmp; + ((mySpline->spline[i])->coeff[1])+= tmp; // @@ -307,7 +315,7 @@ // // From (3) - (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / 6.0; + (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x->data.F32[i+1] / (6.0 * H); // From (4) - (mySpline->spline[i])->coeff[2]-= ((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / 6.0; + ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x->data.F32[i] / (6.0 * H)); // @@ -317,5 +325,5 @@ (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H); // From (4) - (mySpline->spline[i])->coeff[3]+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); + ((mySpline->spline[i])->coeff[3])+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, @@ -334,4 +342,48 @@ return(mySpline); } + +/****************************************************************************** +p_psNRSpline1DEval(): This routine does NR-style evaluation of cubic splines. +It takes advantage of the 2nd derivatives of the cubic splines, which are +stored in the psSPline1D data structure, and computes the interpolated value +directly, without computing (or using) the interpolating cubic spline +polynomial. + +This routine is here mostly for a sanity check on the psLib function +evalSpline() which computes the interpolated value based on the cubic spline +polynomials which are stored in psSpline1D. + *****************************************************************************/ +float p_psNRSpline1DEval(psSpline1D *spline, + const psVector* restrict x, + const psVector* restrict y, + float X) +{ + int n; + int klo; + int khi; + float H; + float A; + float B; + float C; + float D; + float Y; + + n = spline->n; + klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X); + khi = klo + 1; + H = (spline->domains)[khi] - (spline->domains)[klo]; + A = ((spline->domains)[khi] - X) / H; + B = (X - (spline->domains)[klo]) / H; + C = ((A*A*A)-A) * (H*H/6.0); + D = ((B*B*B)-B) * (H*H/6.0); + + Y = (A * y->data.F32[klo]) + + (B * y->data.F32[khi]) + + (C * (spline->p_psDeriv2)[klo]) + + (D * (spline->p_psDeriv2)[khi]); + + return(Y); +} + /****************************************************************************** @@ -687,10 +739,9 @@ } - -#define STEP_SIZE 0.10 /****************************************************************************** - This routine takes as input an arbitrary function, and the parameter to - vary, and the line along which it must vary. This function produces as - output a bracket [a, b, c] such that +p_psDetermineBracket(): This routine takes as input an arbitrary function, +and the parameter to vary, and the line along which it must vary. This +function produces as output a bracket [a, b, c] such that + f(param + b * line) is less than f(param + a * line) and f(param + c * line). @@ -700,4 +751,6 @@ smaller/larger than b. Repeat this process until a local minimum is found. + + *****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, @@ -719,5 +772,5 @@ float new_cDir = 0.0; psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); - float stepSize = STEP_SIZE; + float stepSize = DETERMINE_BRACKET_STEP_SIZE; psVector *tmp = NULL; int i = 0; @@ -869,15 +922,4 @@ return(NULL); } - - - - - - - - - - - /****************************************************************************** @@ -1061,17 +1103,4 @@ - - - - - - - - - - - - - /****************************************************************************** This routine must minimize a possibly multi-dimensional function. The Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 1858) +++ /trunk/psLib/src/math/psPolynomial.c (revision 1859) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-22 01:30:21 $ + * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-23 04:56:40 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -467,4 +467,6 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { if (myPoly->mask[loop_x] == 0) { + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10, + "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]); polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -1687,5 +1689,5 @@ /***************************************************************************** -VectorBinDisectF32(): This is a private function which takes as input a +p_psVectorBinDisectF32(): This is a private function which takes as input a vector of floating point data as well as a single floating point values. The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for @@ -1696,7 +1698,7 @@ XXX: Macro this for a few different types. *****************************************************************************/ -int VectorBinDisectF32(float *bins, - int numBins, - float x) +int p_psVectorBinDisectF32(float *bins, + int numBins, + float x) { int min; @@ -1704,11 +1706,11 @@ int mid; - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Calling VectorBinDisectF32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Calling p_psVectorBinDisectF32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", + "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1720,11 +1722,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); if (x == bins[mid]) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Exiting VectorBinDisectF32(): bin %d\n", mid); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid); return(mid); } else if (x < bins[mid]) { @@ -1736,15 +1738,15 @@ } - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Exiting VectorBinDisectF32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min); return(min); } /***************************************************************************** -VectorBinDisectS32(): integer version of above. +p_psVectorBinDisectS32(): integer version of above. *****************************************************************************/ -int VectorBinDisectS32(int *bins, - int numBins, - int x) +int p_psVectorBinDisectS32(int *bins, + int numBins, + int x) { int min; @@ -1752,11 +1754,11 @@ int mid; - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Calling VectorBinDisectS32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Calling p_psVectorBinDisectS32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", + "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1768,11 +1770,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); if (x == bins[mid]) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Exiting VectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); return(min); } else if (x < bins[mid]) { @@ -1784,11 +1786,11 @@ } - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Exiting VectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); return(min); } /***************************************************************************** -p_psVectorBinDisect(): A wrapper to the above VectorBinDisect(). +p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). *****************************************************************************/ int p_psVectorBinDisect(psVector *bins, @@ -1801,7 +1803,7 @@ if (x->type.type == PS_TYPE_S32) { - return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); + return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); } else if (x->type.type == PS_TYPE_F32) { - return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); + return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); } else { // XXX: Generate error message. @@ -1879,5 +1881,5 @@ "---- p_ps1DInterpolateF32() begin ----\n"); - binNum = VectorBinDisectF32(domain, n, x); + binNum = p_psVectorBinDisectF32(domain, n, x); if (0 == numIntPoints%2) { @@ -1953,5 +1955,5 @@ n = spline->n; - binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x); + binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x); if (binNum == -1) { psLogMsg(__func__, PS_LOG_WARN, Index: /trunk/psLib/src/math/psPolynomial.h =================================================================== --- /trunk/psLib/src/math/psPolynomial.h (revision 1858) +++ /trunk/psLib/src/math/psPolynomial.h (revision 1859) @@ -12,6 +12,6 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-09-22 01:30:21 $ +* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-09-23 04:56:41 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -407,4 +407,12 @@ psVector *data); +int p_psVectorBinDisectF32(float *bins, + int numBins, + float x); + +int p_psVectorBinDisectS32(int *bins, + int numBins, + int x); + int p_psVectorBinDisect(psVector *bins, psScalar *x); @@ -415,4 +423,9 @@ psScalar *x); +float p_psNRSpline1DEval(psSpline1D *spline, + const psVector* restrict x, + const psVector* restrict y, + float X); + /* \} */// End of MathGroup Functions Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 1858) +++ /trunk/psLib/src/math/psSpline.c (revision 1859) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-22 01:30:21 $ + * @version $Revision: 1.39 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-23 04:56:40 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -467,4 +467,6 @@ for (loop_x = 0; loop_x < myPoly->n; loop_x++) { if (myPoly->mask[loop_x] == 0) { + psTrace(".psLib.dataManip.psFunctions.p_psOrdPolynomial1DEval", 10, + "polysum+= sum*coeff [%f+= (%f * %f)\n", polySum, xSum, myPoly->coeff[loop_x]); polySum += xSum * myPoly->coeff[loop_x]; xSum *= x; @@ -1687,5 +1689,5 @@ /***************************************************************************** -VectorBinDisectF32(): This is a private function which takes as input a +p_psVectorBinDisectF32(): This is a private function which takes as input a vector of floating point data as well as a single floating point values. The input vector values are assumed to be non-decreasing (v[i-1] <= v[j] for @@ -1696,7 +1698,7 @@ XXX: Macro this for a few different types. *****************************************************************************/ -int VectorBinDisectF32(float *bins, - int numBins, - float x) +int p_psVectorBinDisectF32(float *bins, + int numBins, + float x) { int min; @@ -1704,11 +1706,11 @@ int mid; - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Calling VectorBinDisectF32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Calling p_psVectorBinDisectF32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "VectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", + "p_psVectorBinDisectF32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1720,11 +1722,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); if (x == bins[mid]) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Exiting VectorBinDisectF32(): bin %d\n", mid); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Exiting p_psVectorBinDisectF32(): bin %d\n", mid); return(mid); } else if (x < bins[mid]) { @@ -1736,15 +1738,15 @@ } - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectF32", 4, - "---- Exiting VectorBinDisectF32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectF32", 4, + "---- Exiting p_psVectorBinDisectF32(): bin %d\n", min); return(min); } /***************************************************************************** -VectorBinDisectS32(): integer version of above. +p_psVectorBinDisectS32(): integer version of above. *****************************************************************************/ -int VectorBinDisectS32(int *bins, - int numBins, - int x) +int p_psVectorBinDisectS32(int *bins, + int numBins, + int x) { int min; @@ -1752,11 +1754,11 @@ int mid; - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Calling VectorBinDisectS32(%f)\n", x); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Calling p_psVectorBinDisectS32(%f)\n", x); if ((x < bins[0]) || (x > bins[numBins-1])) { psLogMsg(__func__, PS_LOG_WARN, - "VectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", + "p_psVectorBinDisectS32(): ordinate %f is outside vector range (%f - %f).", x, bins[0], bins[numBins-1]); return(-1); @@ -1768,11 +1770,11 @@ while (min != max) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", min, mid, max, x, bins[mid]); if (x == bins[mid]) { - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Exiting VectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); return(min); } else if (x < bins[mid]) { @@ -1784,11 +1786,11 @@ } - psTrace(".psLib.dataManip.psFunctions.VectorBinDisectS32", 4, - "---- Exiting VectorBinDisectS32(): bin %d\n", min); + psTrace(".psLib.dataManip.psFunctions.p_psVectorBinDisectS32", 4, + "---- Exiting p_psVectorBinDisectS32(): bin %d\n", min); return(min); } /***************************************************************************** -p_psVectorBinDisect(): A wrapper to the above VectorBinDisect(). +p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). *****************************************************************************/ int p_psVectorBinDisect(psVector *bins, @@ -1801,7 +1803,7 @@ if (x->type.type == PS_TYPE_S32) { - return(VectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); + return(p_psVectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); } else if (x->type.type == PS_TYPE_F32) { - return(VectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); + return(p_psVectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); } else { // XXX: Generate error message. @@ -1879,5 +1881,5 @@ "---- p_ps1DInterpolateF32() begin ----\n"); - binNum = VectorBinDisectF32(domain, n, x); + binNum = p_psVectorBinDisectF32(domain, n, x); if (0 == numIntPoints%2) { @@ -1953,5 +1955,5 @@ n = spline->n; - binNum = VectorBinDisectF32(spline->domains, (spline->n)+1, x); + binNum = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, x); if (binNum == -1) { psLogMsg(__func__, PS_LOG_WARN, Index: /trunk/psLib/src/math/psSpline.h =================================================================== --- /trunk/psLib/src/math/psSpline.h (revision 1858) +++ /trunk/psLib/src/math/psSpline.h (revision 1859) @@ -12,6 +12,6 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-09-22 01:30:21 $ +* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-09-23 04:56:41 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -407,4 +407,12 @@ psVector *data); +int p_psVectorBinDisectF32(float *bins, + int numBins, + float x); + +int p_psVectorBinDisectS32(int *bins, + int numBins, + int x); + int p_psVectorBinDisect(psVector *bins, psScalar *x); @@ -415,4 +423,9 @@ psScalar *x); +float p_psNRSpline1DEval(psSpline1D *spline, + const psVector* restrict x, + const psVector* restrict y, + float X); + /* \} */// End of MathGroup Functions