Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 1848) +++ 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