Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 4729) +++ /trunk/psLib/src/math/psMinimize.c (revision 4730) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.127 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-19 02:55:54 $ + * @version $Revision: 1.128 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-08-08 21:42:07 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -54,4 +54,59 @@ /* FUNCTION IMPLEMENTATION - LOCAL */ /*****************************************************************************/ + +// measure the distance to the minimum assuming a linear model +bool psMinimizeGaussNewtonDelta (psVector *delta, + const psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func) +{ + + // allocate internal arrays (current vs Guess) + psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); + psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *dy = psVectorAlloc (y->n, PS_TYPE_F32); + + // the user provides the error or NULL. we need to convert + // to appropriate weights + if (yErr != NULL) { + for (int i = 0; i < dy->n; i++) { + dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]); + } + } else { + for (int i = 0; i < dy->n; i++) { + dy->data.F32[i] = 1.0; + } + } + + p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func); + p_psMinLM_GuessABP (Alpha, delta, Params, alpha, beta, params, paramMask, 0.0); + + psFree (alpha); + psFree (Alpha); + psFree (beta); + psFree (Params); + psFree (dy); + return (true); +} + + +// measure linear model prediction +psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda) +{ + + /* get linear model prediction */ + psF64 dLinear = 0; + psF64 *B = Beta->data.F64; + psF64 *b = beta->data.F64; + for (int i = 0; i < beta->n; i++) { + dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i]; + } + return (0.5*dLinear); +} /****************************************************************************** @@ -642,4 +697,6 @@ + + /****************************************************************************** XXX: We assume unnormalized gaussians. @@ -769,15 +826,27 @@ } +//XXX: What's this for? psF64 p_psImageGetElementF64(psImage *a, int i, int j); - +/****************************************************************************** +psMinimizeLMChi2(): This routine will take an procedure which calculates +an arbitrary function and it's derivative and minimize the chi-squared match +between that function at the specified coords and the specified value at +those coords. + // XXX EAM this is my re-implementation of MinLM -bool psMinimizeLMChi2(psMinimization *min, - psImage *covar, - psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yErr, - psMinimizeLMChi2Func func) + +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. + *****************************************************************************/ +psBool psMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func) { PS_ASSERT_PTR_NON_NULL(min, NULL); @@ -799,11 +868,12 @@ psVector *beta = psVectorAlloc (params->n, PS_TYPE_F64); psVector *Beta = psVectorAlloc (params->n, PS_TYPE_F64); - psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64); + psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32); psVector *dy = NULL; - psF64 chisq = 0.0; psF64 Chisq = 0.0; psF64 lambda = 0.001; - // the initial guess on params is provided by the user + // XXX EAM: why is this needed here? the value is not used, and the memory + // is allocated above. However, if I drop it, I get weird answers or + // crashes. Params = psVectorCopy (Params, params, PS_TYPE_F32); @@ -822,36 +892,61 @@ // calculate initial alpha and beta, set chisq (min->value) - min->value = p_psMinLM_SetABX (alpha, beta, params, x, y, dy, func); + min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func); # ifndef PS_NO_TRACE // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) { + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { p_psImagePrint (psTraceGetDestination(), alpha, "alpha guess"); p_psVectorPrint (psTraceGetDestination(), beta, "beta guess"); p_psVectorPrint (psTraceGetDestination(), params, "params guess"); } + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { + // XXX: Where is this? + // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); + } # endif /* PS_NO_TRACE */ - + min->lastDelta = min->tol + 1.0; // iterate until the tolerance is reached, or give up while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { - // set a new guess for Alpha, Beta, Params - p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, lambda); + p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, lambda); + + // measure linear model prediction + psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda); # ifndef PS_NO_TRACE // dump some useful info if trace is defined - if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) { + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); } + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { + // XXX: Where is this? + // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); + } # endif /* PS_NO_TRACE */ // calculate Chisq for new guess, update Alpha & Beta - Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, x, y, dy, func); - psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta); - + Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func); + + // XXX EAM alternate convergence criterion: + // compare the delta (min->value - Chisq) with the + // expected delta from the linear model (dLinear) // accept new guess (if improvement), or increase lambda - if (Chisq < min->value) { + psF64 rho = (min->value - Chisq) / dLinear; + + psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho); + # ifndef PS_NO_TRACE + // dump some useful info if trace is defined + if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) { + p_psImagePrint (psTraceGetDestination(), Alpha, "alpha guess"); + p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess"); + p_psVectorPrint (psTraceGetDestination(), Params, "params guess"); + } + # endif /* PS_NO_TRACE */ + + /* if (Chisq < min->value) { */ + if (rho > 0.0) { min->lastDelta = (min->value - Chisq) / (dy->n - params->n); min->value = Chisq; @@ -864,6 +959,13 @@ } min->iter ++; - } - psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta); + + // printf("CONDITIONS: (%f > %f) && (%d < %d)\n", min->lastDelta, min->tol, min->iter, min->maxIter); + } + psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter); + + // construct & return the covariance matrix (if requested) + if (covar != NULL) { + p_psMinLM_GuessABP (covar, Beta, Params, alpha, beta, params, paramMask, 0.0); + } // free the internal temporary data @@ -874,4 +976,8 @@ psFree (Params); psFree (dy); + + if (min->iter == min->maxIter) { + return (false); + } return (true); } @@ -882,5 +988,6 @@ psF64 p_psMinLM_SetABX (psImage *alpha, psVector *beta, - psVector *params, + const psVector *params, + const psVector *paramMask, const psArray *x, const psVector *y, @@ -898,7 +1005,7 @@ for (int j = 0; j < params->n; j++) { for (int k = 0; k < params->n; k++) { - alpha->data.F64[j][k] = 0.0; - } - beta->data.F64[j] = 0.0; + alpha->data.F64[j][k] = 0; + } + beta->data.F64[j] = 0; } chisq = 0.0; @@ -912,6 +1019,10 @@ for (int j = 0; j < params->n; j++) { + if ((paramMask != NULL) && (paramMask->data.U8[j])) + continue; weight = deriv->data.F32[j] * dy->data.F32[i]; for (int k = 0; k <= j; k++) { + if ((paramMask != NULL) && (paramMask->data.U8[k])) + continue; alpha->data.F64[j][k] += weight * deriv->data.F32[k]; } @@ -926,4 +1037,15 @@ } } + + // fill in pivots if we apply a mask + if (paramMask != NULL) { + for (int j = 0; j < params->n; j++) { + if (paramMask->data.U8[j]) { + alpha->data.F64[j][j] = 1; + beta->data.F64[j] = 1; + } + } + } + psFree (deriv); return (chisq); @@ -934,7 +1056,8 @@ psVector *Beta, psVector *Params, - psImage *alpha, - psVector *beta, - psVector *params, + const psImage *alpha, + const psVector *beta, + const psVector *params, + const psVector *paramMask, psF64 lambda) { @@ -948,9 +1071,11 @@ // LU decomposition version - psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using LUD version"); + psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n"); // set new guess values (creates matrix A) A = psImageCopy (NULL, alpha, PS_TYPE_F64); for (int j = 0; j < params->n; j++) { + if ((paramMask != NULL) && (paramMask->data.U8[j])) + continue; A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); } @@ -964,5 +1089,5 @@ # else // gauss-jordan version - psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using Gauss-J version"); + psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version"); // set new guess values (creates matrix A) @@ -970,4 +1095,6 @@ Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64); for (int j = 0; j < params->n; j++) { + if ((paramMask != NULL) && (paramMask->data.U8[j])) + continue; Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); } @@ -976,6 +1103,8 @@ # endif - // apply beta to get new params values + // apply Beta to get new Params values for (int j = 0; j < params->n; j++) { + if ((paramMask != NULL) && (paramMask->data.U8[j])) + continue; Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; } @@ -1095,282 +1224,4 @@ psFree (indxc); return (false); -} - -/****************************************************************************** -psMinimizeLMChi2(): This routine will take an procedure which calculates -an arbitrary function and it's derivative and minimize the chi-squared match -between that function at the specified coords and the specified value at -those coords. - -XXX: Do this: - After checking that all entries in the paramMask are 1 or 0, when - forming the A matrix from alpha, try this: - - 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. - -XXX: Must update the covar matrix. - *****************************************************************************/ -psBool psMinimizeLMChi2Old(psMinimization *min, - psImage *covar, - psVector *params, - const psVector *paramMask, - const psArray *x, - const psVector *y, - const psVector *yErr, - psMinimizeLMChi2Func func) -{ - PS_ASSERT_PTR_NON_NULL(min, NULL); - PS_ASSERT_VECTOR_NON_NULL(params, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(params, NULL); - PS_ASSERT_PTR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL); - PS_ASSERT_PTR_NON_NULL(func, NULL); - - if (paramMask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, NULL); - } - if (yErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, NULL); - } - if (covar != NULL) { - PS_ASSERT_IMAGE_SIZE(covar, params->n, params->n, NULL); - } - - psTrace(".psLib.dataManip.psMinimize", 4, - "---- psMinimizeLMChi2() begin ----\n"); - psS32 numData = y->n; - psS32 numParams = params->n; - psS32 i; - psS32 j; - psS32 k; - psS32 l; - psS32 n; - psS32 p; - psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64); - psVector *perm = NULL; - - psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64); - psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32); - psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32); - - psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32); - psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); - psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); - psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32); - psVector *currValueVec = NULL; - psVector *newValueVec = NULL; - psF32 currChi2 = 0.0; - psF32 newChi2 = 0.0; - psF32 lamda = 0.00005; // XXX EAM : this starting value is VERY small (lamda is mis-spelt) - lamda = 0.05; // XXX EAM : this starting value is quite large (lamda is mis-spelt) - - psTrace(".psLib.dataManip.psMinimize", 6, - "min->maxIter is %d\n", min->maxIter); - psTrace(".psLib.dataManip.psMinimize", 6, - "min->tol is %f\n", min->tol); - - for (p=0;pdata.F32[p] = params->data.F32[p]; - } - - min->lastDelta = PS_MAX_F32; - min->iter = 0; - - while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { - psTrace(".psLib.dataManip.psMinimize", 4, - "------------------------------------------------------\n"); - psTrace(".psLib.dataManip.psMinimize", 4, - "Iteration %d. Delta is %f\n", min->iter, min->lastDelta); - - // - // Calculate the current values and chi-squared of the function. - // - currChi2 = 0.0; - // currValueVec = func(deriv, params, x); - - // XXX EAM: use BinaryOp ? - // t1 = BinaryOp (NULL, currValueVec, "-", y); - // t1 = BinaryOp (t1, t1, "*", t1); - - // XXX EAM: this ignores yErr - for (n=0;ndata.F32[n] - y->data.F32[n]) * - (currValueVec->data.F32[n] - y->data.F32[n]); - psTrace(".psLib.dataManip.psMinimize", 6, - "data[%d], chi2 calculation+= (%f - %f)^2\n", n, - currValueVec->data.F32[n], y->data.F32[n]); - } - - // XXX EAM: this is just for tracing - for (p=0;pdata.F32[%d] is %f.\n", p, params->data.F32[p]); - } - psTrace(".psLib.dataManip.psMinimize", 6, - "Current chi-squared is (%f)\n", currChi2); - - // - // Mask elements of the derivative for each data point. - // XXX EAM : is this necessary? probably not... - for (p=0;pdata.U8[p] != 0)) { - for (n=0;ndata.F32[n][p] = 0.0; - } - } - } - - // - // Calculate the BETA vector. - // XXX EAM: I think this is wrong - for (p=0;pdata.U8[p] != 0)) { - continue; - } - beta->data.F64[p] = 0.0; - for (n=0;ndata.F64[p])+= - (y->data.F32[n] - currValueVec->data.F32[n]) * - deriv->data.F32[n][p]; - } - // XXX: multiply by -1 here? - (beta->data.F64[p])*= -1.0; - psTrace(".psLib.dataManip.psMinimize", 6, - "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]); - } - psFree(currValueVec); - - // - // Calculate the ALPHA matrix. - // XXX EAM: also wrong? (missing yErr) - for (k=0;kdata.F32[k][l] = 0.0; - for (n=0;ndata.F32[k][l]+= deriv->data.F32[n][k] * - deriv->data.F32[n][l]; - } - } - } - - // - // Calculate the matrix A. - // - for (j=0;jdata.F64[j][k] = - (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]); - } else { - A->data.F64[j][k] = (psF64) alpha->data.F32[j][k]; - } - } - } - for (j=0;jdata.F64[j][k]); - } - psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n"); - } - - // - // Solve A * alpha = Beta - // - // XXX: How do we know if these functions were successful? - // - aOut = psMatrixLUD(aOut, &perm, A); - paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm); - - // - // Mask any masked parameters. - // - for (i=0;idata.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]); - if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { - newParams->data.F32[i] = origParams->data.F32[i]; - } else { - newParams->data.F32[i] = params->data.F32[i] - - (psF32) paramDeltasF64->data.F64[i]; - } - } - - psTrace(".psLib.dataManip.psMinimize", 6, - "Calling func() with new parameters:\n"); - for (i=0;idata.F32[%d] is %f.\n", i, newParams->data.F32[i]); - } - - - // - // Calculate new function values. - // - newChi2 = 0.0; - // newValueVec = func(deriv, newParams, x); - for (n=0;ndata.F32[n] - y->data.F32[n]) * - (newValueVec->data.F32[n] - y->data.F32[n]); - - } - psFree(newValueVec); - - psTrace(".psLib.dataManip.psMinimize", 4, - "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2); - - // - // If the new chi-squared is lower, then keep it. - // - if (currChi2 > newChi2) { - min->lastDelta = (currChi2 - newChi2)/currChi2; - min->value = newChi2; - - // We already masked params. - for (i=0;idata.F32[i] = (psF32) newParams->data.F32[i]; - } - lamda*= 0.1; - psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n"); - } else { - lamda*= 10.0; - psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n"); - } - psTrace(".psLib.dataManip.psMinimize", 4, - "lamda is %f\n", lamda); - min->iter++; - } - psFree(beta); - psFree(perm); - psFree(paramDeltasF64); - psFree(origParams); - psFree(newParams); - psFree(alpha); - psFree(A); - psFree(aOut); - psFree(deriv); - - if ((min->iter < min->maxIter) || - (min->lastDelta <= min->tol)) { - return(true); - } - - psTrace(".psLib.dataManip.psMinimize", 4, - "---- psMinimizeLMChi2() end (false) ----\n"); - return(false); } Index: /trunk/psLib/src/math/psMinimize.h =================================================================== --- /trunk/psLib/src/math/psMinimize.h (revision 4729) +++ /trunk/psLib/src/math/psMinimize.h (revision 4730) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-16 00:42:28 $ + * @version $Revision: 1.51 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-08-08 21:42:07 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -135,7 +135,8 @@ psVector *Beta, ///< New Beta guess psVector *Params, ///< New Params guess - psImage *alpha, ///< Old Alpha guess - psVector *beta, ///< Old Beta guess - psVector *params, ///< Old Params guess + const psImage *alpha, ///< Old Alpha guess + const psVector *beta, ///< Old Beta guess + const psVector *params, ///< Old Params guess + const psVector *paramMask, ///< Param Mask psF64 lambda ///< Factor used in update //XXX Got Better Desc? @@ -149,5 +150,6 @@ psImage *alpha, ///< alpha guess psVector *beta, ///< beta guess - psVector *params, ///< params guess + const psVector *params, ///< params guess + const psVector *paramMask, ///< param mask const psArray *x, ///< Measurement ordinates const psVector *y, ///< Measurement coordinates @@ -155,4 +157,5 @@ psMinimizeLMChi2Func func ///< Specified function ); + /** Specifies the format of a user-defined function that the general Powell