Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 2106) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 2197) @@ -9,6 +9,6 @@ * @author GLF, MHPCC * - * @version $Revision: 1.58 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 01:44:48 $ + * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -636,4 +636,261 @@ psFree(currValueVec); psFree(newValueVec); + + if ((min->iter < min->maxIter) || + (min->lastDelta <= min->tol)) { + return(true); + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "---- psMinimizeLMChi2() end (false) ----\n"); + return(false); +} + +/****************************************************************************** +psMyMinimizeLMChi2(): 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. + *****************************************************************************/ +bool psMyMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMyMinimizeLMChi2Func func) +{ + PS_CHECK_NULL_PTR_RETURN_NULL(min); + PS_CHECK_NULL_VECTOR_RETURN_NULL(params); + PS_CHECK_EMPTY_VECTOR_RETURN_NULL(params); + PS_CHECK_NULL_PTR_RETURN_NULL(x); + PS_CHECK_NULL_VECTOR_RETURN_NULL(y); + PS_CHECK_EMPTY_VECTOR_RETURN_NULL(y); + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y); + if (paramMask != NULL) { + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(params, paramMask); + } + if (yErr != NULL) { + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(y, yErr); + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "---- psMinimizeLMChi2() begin ----\n"); + int numData = y->n; + int numParams = params->n; + int i; + int j; + int k; + int l; + int n; + int p; + psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64); + + 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); + + // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *)); + // for (i=0;imaxIter is %d\n", min->maxIter); + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "min->tol is %f\n", min->tol); + + for (p=0;pdata.F32[p] = params->data.F32[p]; + } + + min->lastDelta = HUGE; + min->iter = 0; + + while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "------------------------------------------------------\n"); + psTrace(".psLib.dataManip.psMinimizeLMChi2", 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); + for (n=0;ndata.F32[n] * currValueVec->data.F32[n]); + } + + for (p=0;pdata.F32[%d] is %f.\n", p, params->data.F32[p]); + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "Current chi-squared is (%f)\n", currChi2); + + // + // Mask elements of the derivative for each data point. + // + for (p=0;pdata.U8[p] != 0)) { + for (n=0;ndata.F32[n][p] = 0.0; + } + } + } + + // + // Calculate the BETA vector. + // + for (p=0;pdata.F64[p] = 0.0; + for (n=0;ndata.F64[p])+= + (y->data.F32[n] - currValueVec->data.F32[n]) * + deriv->data.F32[n][p]; + } + // XXX: multiple by -1 here? + (beta->data.F64[p])*= -1.0; + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]); + } + psFree(currValueVec); + + // + // Calculate the ALPHA matrix. + // + 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] = + (double) ((1.0 + lamda) * alpha->data.F32[j][k]); + } else { + A->data.F64[j][k] = (double) alpha->data.F32[j][k]; + } + } + } + for (j=0;jdata.F64[j][k]); + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n"); + } + + // + // Solve A * alpha = Beta + // + 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] - + (float) paramDeltasF64->data.F64[i]; + } + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 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] * newValueVec->data.F32[n]); + } + psFree(newValueVec); + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 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; + min->value = newChi2; + + // We already masked params. + for (i=0;idata.F32[i] = (float) newParams->data.F32[i]; + } + lamda*= 0.1; + } else { + lamda*= 10.0; + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 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) ||