Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 2741) +++ trunk/psLib/src/math/psMinimize.c (revision 2788) @@ -9,14 +9,10 @@ * @author GLG, MHPCC * - * @version $Revision: 1.95 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-17 00:18:31 $ + * @version $Revision: 1.96 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * 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 redefine the LMM functions. * */ @@ -64,10 +60,10 @@ XXX: Use a static vector. *****************************************************************************/ -void psBuildSums1D(double x, +void psBuildSums1D(psF64 x, psS32 polyOrder, psVector* sums) { psS32 i = 0; - double xSum = 0.0; + psF64 xSum = 0.0; if (sums == NULL) { @@ -100,5 +96,5 @@ XXX: do an F64 version? *****************************************************************************/ -float *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) +psF32 *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) const psVector* restrict y) ///< Coordinates { @@ -108,12 +104,12 @@ psS32 i; psS32 k; - float sig; - float p; + psF32 sig; + psF32 p; psS32 n = y->n; - float *u = (float *) psAlloc(n * sizeof(float)); - float *derivs2 = (float *) psAlloc(n * sizeof(float)); - float *X = (float *) & (x->data.F32[0]); - float *Y = (float *) & (y->data.F32[0]); - float qn; + psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *X = (psF32 *) & (x->data.F32[0]); + psF32 *Y = (psF32 *) & (y->data.F32[0]); + psF32 qn; // XXX: The second derivatives at the endpoints, undefined in the SDR, @@ -168,8 +164,8 @@ *****************************************************************************/ /* -float p_psNRSpline1DEval(psSpline1D *spline, +psF32 p_psNRSpline1DEval(psSpline1D *spline, const psVector* restrict x, const psVector* restrict y, - float X) + psF32 X) { PS_PTR_CHECK_NULL(spline, NAN); @@ -185,10 +181,10 @@ psS32 klo; psS32 khi; - float H; - float A; - float B; - float C; - float D; - float Y; + psF32 H; + psF32 A; + psF32 B; + psF32 C; + psF32 D; + psF32 Y; n = spline->n; @@ -258,8 +254,8 @@ "---- psVectorFitSpline1D() begin ----\n"); psS32 numSplines = (y->n)-1; - float tmp; - float H; + psF32 tmp; + psF32 H; psS32 i; - float slope; + psF32 slope; psVector *x32 = NULL; psVector *y32 = NULL; @@ -421,8 +417,8 @@ psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2Gauss1D() begin ----\n"); - float x; - int i; - float mean = params->data.F32[0]; - float stdev = params->data.F32[1]; + psF32 x; + psS32 i; + psF32 mean = params->data.F32[0]; + psF32 stdev = params->data.F32[1]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -443,5 +439,5 @@ for (i=0;in;i++) { x = ((psVector *) (coords->data[i]))->data.F32[0]; - float tmp = (x - mean) * psGaussian(x, mean, stdev, false); + psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false); deriv->data.F32[i][0] = tmp / (stdev * stdev); tmp = (x - mean) * (x - mean) * @@ -491,10 +487,10 @@ PS_PTR_CHECK_NULL(params, NULL); - double normalization = params->data.F32[0]; - double x0 = params->data.F32[1]; - double y0 = params->data.F32[2]; - double sigmaX = params->data.F32[3]; - double sigmaY = params->data.F32[4]; - double theta = params->data.F32[5]; + psF64 normalization = params->data.F32[0]; + psF64 x0 = params->data.F32[1]; + psF64 y0 = params->data.F32[2]; + psF64 sigmaX = params->data.F32[3]; + psF64 sigmaY = params->data.F32[4]; + psF64 theta = params->data.F32[5]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -502,5 +498,5 @@ deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); } else { - // XXX: Check size of derivative + PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL); } @@ -508,14 +504,14 @@ "---- psMinimizeLMChi2Gauss2D() begin ----\n"); - for (int i=0;in;i++) { - double x = ((psVector *) coords->data[i])->data.F32[0]; - double y = ((psVector *) coords->data[i])->data.F32[0]; - - double u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); - double v = (x-x0)*cos(theta) + (y-y0)*sin(theta); - - double flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + - v*v/(2.0 * sigmaY * sigmaY)))/ - (2.0 * M_PI * sigmaX * sigmaY); + for (psS32 i=0;in;i++) { + psF64 x = ((psVector *) coords->data[i])->data.F32[0]; + psF64 y = ((psVector *) coords->data[i])->data.F32[0]; + + psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); + psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta); + + psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + + v*v/(2.0 * sigmaY * sigmaY)))/ + (2.0 * M_PI * sigmaX * sigmaY); out->data.F32[i] = flux; @@ -584,6 +580,4 @@ PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL); } - - // XXX: Generate code that modifies covar matrix if not-NULL. if (covar != NULL) { PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL); @@ -610,17 +604,10 @@ psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); - - // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *)); - // for (i=0;idata.F64[j][k] = - (double) ((1.0 + lamda) * alpha->data.F32[j][k]); + (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]); } else { - A->data.F64[j][k] = (double) alpha->data.F32[j][k]; + A->data.F64[j][k] = (psF64) alpha->data.F32[j][k]; } } @@ -741,5 +728,5 @@ } else { newParams->data.F32[i] = params->data.F32[i] - - (float) paramDeltasF64->data.F64[i]; + (psF32) paramDeltasF64->data.F64[i]; } } @@ -777,5 +764,5 @@ // We already masked params. for (i=0;idata.F32[i] = (float) newParams->data.F32[i]; + params->data.F32[i] = (psF32) newParams->data.F32[i]; } lamda*= 0.1; @@ -826,6 +813,6 @@ psS32 k; psS32 n = x->n; - double fac; - double sum; + psF64 fac; + psF64 sum; PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64); psScalar *fScalar; @@ -836,8 +823,8 @@ // variable declarations. I retain them here to maintain coherence // with the NR code. - double min = -1.0; - double max = 1.0; - double bma = 0.5 * (max-min); // 1 - double bpa = 0.5 * (max+min); // 0 + psF64 min = -1.0; + psF64 max = 1.0; + psF64 bma = 0.5 * (max-min); // 1 + psF64 bpa = 0.5 * (max+min); // 0 // In this loop, we first calculate the values of X for which the @@ -851,6 +838,6 @@ for (psS32 i=0;in;j++) { @@ -879,5 +866,5 @@ for (k=0;kdata.F64[k] * - cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n)); + cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n)); } @@ -1075,5 +1062,5 @@ *****************************************************************************/ psMinimization *psMinimizationAlloc(psS32 maxIter, - float tol) + psF32 tol) { PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL); @@ -1092,5 +1079,5 @@ // LINE to it. We assume BASEMASK is non-null. #define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if (BASEMASK->data.U8[i] == 0) { \ OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ @@ -1102,5 +1089,5 @@ #define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \ BOOL_VAR = true; \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \ BOOL_VAR = false; \ @@ -1111,5 +1098,5 @@ #define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \ BOOL_VAR = true; \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \ BOOL_VAR = false; \ @@ -1152,17 +1139,17 @@ psMinimizePowellFunc func) { - float a = 0.0; - float b = 0.0; - float c = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; psS32 iter = 100; - float aDir = 0.0; - float cDir = 0.0; - float new_aDir = 0.0; - float new_cDir = 0.0; + psF32 aDir = 0.0; + psF32 cDir = 0.0; + psF32 new_aDir = 0.0; + psF32 new_cDir = 0.0; psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); - float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; + psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; psVector *tmp = NULL; psBool boolLineIsNull = true; @@ -1301,15 +1288,15 @@ psMinimizePowellFunc func) { - float a = 0.0; - float b = 0.0; - float c = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; psS32 iter = 0; PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32); psBool boolLineIsNull = true; - float prevMin = 0.0; - int countMin = 0; + psF32 prevMin = 0.0; + psS32 countMin = 0; psTrace(".psLib.dataManip.p_psDetermineBracket", 4, @@ -1403,5 +1390,5 @@ *****************************************************************************/ #define PS_LINEMIN_MAX_ITERATIONS 30 -float p_psLineMin(psMinimization *min, +psF32 p_psLineMin(psMinimization *min, psVector *params, psVector *line, @@ -1423,13 +1410,13 @@ PS_PTR_CHECK_NULL(func, NAN); psVector *bracket; - float a = 0.0; - float b = 0.0; - float c = 0.0; - float n = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; - float fn = 0.0; - float mul = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 n = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; + psF32 fn = 0.0; + psF32 mul = 0.0; PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32); PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32); @@ -1544,5 +1531,5 @@ This routine must minimize a possibly multi-dimensional function. The function to be minimized "func" is: - float func(psVector *params, psArray *coords) + psF32 func(psVector *params, psArray *coords) The "params" are the parameters of the function which are varied. The data points at which the function is varied are in the argument "coords" which is @@ -1578,10 +1565,10 @@ psVector *myParamMask = NULL; psMinimization dummyMin; - float mul = 0.0; - float baseFuncVal = 0.0; - float currFuncVal = 0.0; + psF32 mul = 0.0; + psF32 baseFuncVal = 0.0; + psF32 currFuncVal = 0.0; psS32 biggestIter = 0; - float biggestDiff = 0.0; - int iterationNumber = 0; + psF32 biggestDiff = 0.0; + psS32 iterationNumber = 0; psTrace(".psLib.dataManip.psMinimizePowell", 4, @@ -1711,9 +1698,9 @@ } } - float fqp = func(pQP, coords); - float term1 = (baseFuncVal - currFuncVal) - biggestDiff; + psF32 fqp = func(pQP, coords); + psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff; term1*= term1; term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp); - float term2 = baseFuncVal - fqp; + psF32 term2 = baseFuncVal - fqp; term2*= term2 * biggestDiff; if (term1 < term2) { @@ -1760,8 +1747,8 @@ PS_PTR_CHECK_NULL(params, NULL); - float x; - int i; - float mean = params->data.F32[0]; - float stdev = params->data.F32[1]; + psF32 x; + psS32 i; + psF32 mean = params->data.F32[0]; + psF32 stdev = params->data.F32[1]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -1783,5 +1770,5 @@ XXX: This is F32 only *****************************************************************************/ -float myPowellChi2Func(const psVector *params, +psF32 myPowellChi2Func(const psVector *params, const psArray *coords) { @@ -1794,6 +1781,6 @@ PS_PTR_CHECK_NULL(coords, NAN); - float chi2 = 0.0; - float d; + psF32 chi2 = 0.0; + psF32 d; psS32 i; psVector *tmp;