Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 1859) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 1860) @@ -9,10 +9,14 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-23 04:56:41 $ + * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-23 05:43:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * * XXX: must follow coding name standards on local functions. + * + * XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via + * LMM) is not addressed here. We are waiting for subsequent SDRs + * which will redfeine the LMM functions. * */ @@ -38,4 +42,5 @@ #include "psImage.h" #include "psTrace.h" +#include "psLogMsg.h" #include "psError.h" #include "psAbort.h" @@ -94,5 +99,4 @@ /*****************************************************************************/ static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL; -//static psMinimizePowellFunc PowellFunc = NULL; static psVector *myValue; static psVector *myError; @@ -121,4 +125,7 @@ int i = 0; double xSum = 0.0; + if (sums == NULL) { + sums = psVectorAlloc(polyOrder, PS_TYPE_F64); + } xSum = 1.0; @@ -139,5 +146,4 @@ tmp = psVectorAlloc(x, PS_TYPE_F32); - for (i = 0; i < x; i++) { tmp->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0; @@ -224,122 +230,4 @@ } - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - PUBLIC */ -/*****************************************************************************/ - -/***************************************************************************** -psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y -vectors, this routine generates the cublic splines which satisfy those data -points. - -The formula for calculating the spline polynomials is derived from Numerical -Recipes in C. The basic idea is that the polynomial is - (1) y = (A * y[0]) + - (2) (B * y[1]) + - (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 -The bulk of the code in this routine is the expansion of the above equation -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. - *****************************************************************************/ - -psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, ///< The spline which will be generated. - const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) - const psVector* restrict y, ///< Coordinates - const psVector* restrict yErr) ///< Errors in coordinates, or NULL -{ - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() begin ----\n"); - int numSplines = (y->n)-1; - float tmp; - float H; - int i; - - if (mySpline == NULL) { - //XXX psErrorMsg() - } - 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;idata.F32[i+1] - x->data.F32[i]; - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H); - // - // ******** Calulate 0-order term ******** - // - // From (1) - (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); - // From (3) - tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); - tmp-= (x->data.F32[i+1] / H); - tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - // From (4) - tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H); - tmp+= (x->data.F32[i] / H); - tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - - // - // ******** Calulate 1-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H; - // From (2) - ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H); - // From (3) - 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; - // From (4) - 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; - - // - // ******** Calulate 2-order term ******** - // - // From (3) - (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 * H)); - - // - // ******** Calulate 3-order term ******** - // - // From (3) - (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); - - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); - - } - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() end ----\n"); - return(mySpline); -} /****************************************************************************** @@ -387,4 +275,175 @@ +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +/***************************************************************************** +psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y +vectors, this routine generates the linear or cublic splines which satisfy +those data points. + +The formula for calculating the spline polynomials is derived from Numerical +Recipes in C. The basic idea is that the polynomial is + (1) y = (A * y[0]) + + (2) (B * y[1]) + + (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 +The bulk of the code in this routine is the expansion of the above equation +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. +XXX: Must do: if yErr==NULL, set all errors equal. + +XXX: Is the x argument redundant? What do we do if the x argument is + supplied, but does not equal the domains specified in mySpline? + +XXX: can psSpline be NULL? + +XXX: Implemented in F32 only, must add F64. + *****************************************************************************/ +psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, ///< The spline which will be generated. + const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) + const psVector* restrict y, ///< Coordinates + const psVector* restrict yErr) ///< Errors in coordinates, or NULL +{ + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() begin ----\n"); + int numSplines = (y->n)-1; + float tmp; + float H; + int i; + float slope; + + if (mySpline == NULL) { + //XXX psErrorMsg() + } + if (y == NULL) { + //XXX psErrorMsg() + } + + if (y->n != (1 + mySpline->n)) { + psLogMsg(__func__, PS_LOG_WARN, + "data size / spline size mismatch (%d %d)\n", + y->n, mySpline->n); + // XXX: psErrorMsg() + } + + // If these are linear splines, which means their polynomials will have + // two coefficients, then we do the simple calculation. + if (2 == (mySpline->spline[0])->n) { + for (i=0;in;i++) { + slope = (y->data.F32[i+1] - y->data.F32[i]) / + (mySpline->domains[i+1] - mySpline->domains[i]); + (mySpline->spline[i])->coeff[0] = y->data.F32[i] - + (slope * mySpline->domains[i]); + + (mySpline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- mySpline %d coeffs are (%f, %f)\n", i, + (mySpline->spline[i])->coeff[0], + (mySpline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Exiting psSpline1DGen()\n"); + return((psSpline1D *) mySpline); + } + + // Check if these are cubic splines (n==4). If not, psError. + if (4 != (mySpline->spline[0])->n) { + psLogMsg(__func__, PS_LOG_WARN, + "Don't know how to generate %d-order splines.", + (mySpline->spline[0])->n-1); + // XXX: psErrorMsg() + return(NULL); + } + + // If we get here, then we know these are cubic splines. We first + // generate the second derivatives at each data point. + 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]); + + // We generate the coefficients of the spline polynomials. I can't + // concisely explain how this code works. See above function comments + // and Numerical Recipes in C. + for (i=0;idata.F32[i+1] - x->data.F32[i]; + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H); + // + // ******** Calculate 0-order term ******** + // + // From (1) + (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); + // From (3) + tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); + tmp-= (x->data.F32[i+1] / H); + tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + // From (4) + tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H); + tmp+= (x->data.F32[i] / H); + tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + + // + // ******** Calculate 1-order term ******** + // + // From (1) + (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H; + // From (2) + ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H); + // From (3) + 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; + // From (4) + 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; + + // + // ******** Calculate 2-order term ******** + // + // From (3) + (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 * H)); + + // + // ******** Calculate 3-order term ******** + // + // From (3) + (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); + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); + + } + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() end ----\n"); + return(mySpline); +} + /****************************************************************************** psMinimizeLMChi2(): This routine will take an procedure which calculates an @@ -398,4 +457,16 @@ A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; + +XXX: This is very different from what is specified in the SDR. Must +coordinate with IfA on new SDR. + +XXX: Do vector/image recycles. + +XXX: probably yErr will be part of the SDR. + +XXX: This must work for both F32 and F64. F32 is currently implemented. + Note: since the LUD routines are only implemented in F64, then we + will have to convert all F32 input vectors to F64 regardless. So, + the F64 port might be. *****************************************************************************/ bool psMinimizeLMChi2(psMinimization *min, @@ -450,5 +521,4 @@ min->lastDelta = HUGE; - min->lastDelta = 12345.0; min->iter = 0; @@ -627,5 +697,11 @@ polynomial, as well as the error for each data point (yErr). -XXX: NOTE: yErr is currently ignored. +XXX: yErr is currently ignored. + +XXX: must add type F32 (currently F64 only). + +XXX: Must do: if x==NULL, use the index vector (???). + +XXX: Must do: if yErr==NULL, set all errors equal. *****************************************************************************/ psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly, @@ -726,4 +802,7 @@ } + +/****************************************************************************** + *****************************************************************************/ psMinimization *psMinimizationAlloc(int maxIter, float tol) @@ -747,10 +826,12 @@ f(param + c * line). - Algorithm: XXX completely ad hoc: start with the user-supplied starting - parameter and call that b. Calculate a/c as a fractional amount - smaller/larger than b. Repeat this process until a local minimum is - found. - - +Algorithm: XXX completely ad hoc: start with the user-supplied starting +parameter and call that b. Calculate a/c as a fractional amount +smaller/larger than b. Repeat this process until a local minimum is found. + +XXX: new algorithm: start at x=0, expand in one direction until the function +decreases. Then you have two points in the bracket. Keep going until it +increases, or x is too large. If thst does not work, expand in the other +direction. *****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, @@ -926,4 +1007,6 @@ This routine must minimize a possibly multi-dimensional function along a vector defined by line. + +XXX: Use a p_psName(). *****************************************************************************/ float psLineMin(psMinimization *min, @@ -1113,5 +1196,7 @@ XXX: We do not use Brent's method. -XXX: We do not use the +XXX: The SDR is silent about data types. F32 is implemented here. + +XXX: Define constants for default dummy number of iterations and tolerance. *****************************************************************************/ bool psMinimizePowell(psMinimization *min, Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 1859) +++ /trunk/psLib/src/math/psMinimize.c (revision 1860) @@ -9,10 +9,14 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-23 04:56:41 $ + * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-23 05:43:25 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * * XXX: must follow coding name standards on local functions. + * + * XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via + * LMM) is not addressed here. We are waiting for subsequent SDRs + * which will redfeine the LMM functions. * */ @@ -38,4 +42,5 @@ #include "psImage.h" #include "psTrace.h" +#include "psLogMsg.h" #include "psError.h" #include "psAbort.h" @@ -94,5 +99,4 @@ /*****************************************************************************/ static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL; -//static psMinimizePowellFunc PowellFunc = NULL; static psVector *myValue; static psVector *myError; @@ -121,4 +125,7 @@ int i = 0; double xSum = 0.0; + if (sums == NULL) { + sums = psVectorAlloc(polyOrder, PS_TYPE_F64); + } xSum = 1.0; @@ -139,5 +146,4 @@ tmp = psVectorAlloc(x, PS_TYPE_F32); - for (i = 0; i < x; i++) { tmp->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0; @@ -224,122 +230,4 @@ } - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - PUBLIC */ -/*****************************************************************************/ - -/***************************************************************************** -psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y -vectors, this routine generates the cublic splines which satisfy those data -points. - -The formula for calculating the spline polynomials is derived from Numerical -Recipes in C. The basic idea is that the polynomial is - (1) y = (A * y[0]) + - (2) (B * y[1]) + - (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 -The bulk of the code in this routine is the expansion of the above equation -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. - *****************************************************************************/ - -psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, ///< The spline which will be generated. - const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) - const psVector* restrict y, ///< Coordinates - const psVector* restrict yErr) ///< Errors in coordinates, or NULL -{ - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() begin ----\n"); - int numSplines = (y->n)-1; - float tmp; - float H; - int i; - - if (mySpline == NULL) { - //XXX psErrorMsg() - } - 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;idata.F32[i+1] - x->data.F32[i]; - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H); - // - // ******** Calulate 0-order term ******** - // - // From (1) - (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); - // From (3) - tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); - tmp-= (x->data.F32[i+1] / H); - tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - // From (4) - tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H); - tmp+= (x->data.F32[i] / H); - tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - - // - // ******** Calulate 1-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H; - // From (2) - ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H); - // From (3) - 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; - // From (4) - 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; - - // - // ******** Calulate 2-order term ******** - // - // From (3) - (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 * H)); - - // - // ******** Calulate 3-order term ******** - // - // From (3) - (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); - - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, - "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); - - } - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() end ----\n"); - return(mySpline); -} /****************************************************************************** @@ -387,4 +275,175 @@ +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +/***************************************************************************** +psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y +vectors, this routine generates the linear or cublic splines which satisfy +those data points. + +The formula for calculating the spline polynomials is derived from Numerical +Recipes in C. The basic idea is that the polynomial is + (1) y = (A * y[0]) + + (2) (B * y[1]) + + (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 +The bulk of the code in this routine is the expansion of the above equation +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. +XXX: Must do: if yErr==NULL, set all errors equal. + +XXX: Is the x argument redundant? What do we do if the x argument is + supplied, but does not equal the domains specified in mySpline? + +XXX: can psSpline be NULL? + +XXX: Implemented in F32 only, must add F64. + *****************************************************************************/ +psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, ///< The spline which will be generated. + const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) + const psVector* restrict y, ///< Coordinates + const psVector* restrict yErr) ///< Errors in coordinates, or NULL +{ + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() begin ----\n"); + int numSplines = (y->n)-1; + float tmp; + float H; + int i; + float slope; + + if (mySpline == NULL) { + //XXX psErrorMsg() + } + if (y == NULL) { + //XXX psErrorMsg() + } + + if (y->n != (1 + mySpline->n)) { + psLogMsg(__func__, PS_LOG_WARN, + "data size / spline size mismatch (%d %d)\n", + y->n, mySpline->n); + // XXX: psErrorMsg() + } + + // If these are linear splines, which means their polynomials will have + // two coefficients, then we do the simple calculation. + if (2 == (mySpline->spline[0])->n) { + for (i=0;in;i++) { + slope = (y->data.F32[i+1] - y->data.F32[i]) / + (mySpline->domains[i+1] - mySpline->domains[i]); + (mySpline->spline[i])->coeff[0] = y->data.F32[i] - + (slope * mySpline->domains[i]); + + (mySpline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- mySpline %d coeffs are (%f, %f)\n", i, + (mySpline->spline[i])->coeff[0], + (mySpline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psFunctions.psSpline1DGen", 4, + "---- Exiting psSpline1DGen()\n"); + return((psSpline1D *) mySpline); + } + + // Check if these are cubic splines (n==4). If not, psError. + if (4 != (mySpline->spline[0])->n) { + psLogMsg(__func__, PS_LOG_WARN, + "Don't know how to generate %d-order splines.", + (mySpline->spline[0])->n-1); + // XXX: psErrorMsg() + return(NULL); + } + + // If we get here, then we know these are cubic splines. We first + // generate the second derivatives at each data point. + 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]); + + // We generate the coefficients of the spline polynomials. I can't + // concisely explain how this code works. See above function comments + // and Numerical Recipes in C. + for (i=0;idata.F32[i+1] - x->data.F32[i]; + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "X data (%f - %f) (%f)\n", x->data.F32[i], x->data.F32[i+1], H); + // + // ******** Calculate 0-order term ******** + // + // From (1) + (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); + // From (3) + tmp = (x->data.F32[i+1] * x->data.F32[i+1] * x->data.F32[i+1]) / (H * H * H); + tmp-= (x->data.F32[i+1] / H); + tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + // From (4) + tmp = -(x->data.F32[i] * x->data.F32[i] * x->data.F32[i]) / (H * H * H); + tmp+= (x->data.F32[i] / H); + tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + + // + // ******** Calculate 1-order term ******** + // + // From (1) + (mySpline->spline[i])->coeff[1] = -(y->data.F32[i]) / H; + // From (2) + ((mySpline->spline[i])->coeff[1])+= (y->data.F32[i+1] / H); + // From (3) + 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; + // From (4) + 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; + + // + // ******** Calculate 2-order term ******** + // + // From (3) + (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 * H)); + + // + // ******** Calculate 3-order term ******** + // + // From (3) + (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); + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); + + } + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() end ----\n"); + return(mySpline); +} + /****************************************************************************** psMinimizeLMChi2(): This routine will take an procedure which calculates an @@ -398,4 +457,16 @@ A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; + +XXX: This is very different from what is specified in the SDR. Must +coordinate with IfA on new SDR. + +XXX: Do vector/image recycles. + +XXX: probably yErr will be part of the SDR. + +XXX: This must work for both F32 and F64. F32 is currently implemented. + Note: since the LUD routines are only implemented in F64, then we + will have to convert all F32 input vectors to F64 regardless. So, + the F64 port might be. *****************************************************************************/ bool psMinimizeLMChi2(psMinimization *min, @@ -450,5 +521,4 @@ min->lastDelta = HUGE; - min->lastDelta = 12345.0; min->iter = 0; @@ -627,5 +697,11 @@ polynomial, as well as the error for each data point (yErr). -XXX: NOTE: yErr is currently ignored. +XXX: yErr is currently ignored. + +XXX: must add type F32 (currently F64 only). + +XXX: Must do: if x==NULL, use the index vector (???). + +XXX: Must do: if yErr==NULL, set all errors equal. *****************************************************************************/ psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly, @@ -726,4 +802,7 @@ } + +/****************************************************************************** + *****************************************************************************/ psMinimization *psMinimizationAlloc(int maxIter, float tol) @@ -747,10 +826,12 @@ f(param + c * line). - Algorithm: XXX completely ad hoc: start with the user-supplied starting - parameter and call that b. Calculate a/c as a fractional amount - smaller/larger than b. Repeat this process until a local minimum is - found. - - +Algorithm: XXX completely ad hoc: start with the user-supplied starting +parameter and call that b. Calculate a/c as a fractional amount +smaller/larger than b. Repeat this process until a local minimum is found. + +XXX: new algorithm: start at x=0, expand in one direction until the function +decreases. Then you have two points in the bracket. Keep going until it +increases, or x is too large. If thst does not work, expand in the other +direction. *****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, @@ -926,4 +1007,6 @@ This routine must minimize a possibly multi-dimensional function along a vector defined by line. + +XXX: Use a p_psName(). *****************************************************************************/ float psLineMin(psMinimization *min, @@ -1113,5 +1196,7 @@ XXX: We do not use Brent's method. -XXX: We do not use the +XXX: The SDR is silent about data types. F32 is implemented here. + +XXX: Define constants for default dummy number of iterations and tolerance. *****************************************************************************/ bool psMinimizePowell(psMinimization *min,