Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 2221) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 2228) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.72 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-27 23:31:43 $ + * @version $Revision: 1.73 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -408,4 +408,54 @@ } + +/****************************************************************************** +XXX: We assume unnormalized gaussians. +XXX: Currently, yErr is ignored. + *****************************************************************************/ +psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, + const psVector *params, + const psArray *coords) +{ + psTrace(".psLib.dataManip.psMinimize", 4, + "---- psMinimizeLMChi2Gauss1D() begin ----\n"); + + PS_PTR_CHECK_NULL(coords, NULL); + PS_PTR_CHECK_NULL(params, NULL); + float x; + int i; + float mean = params->data.F32[0]; + float stdev = params->data.F32[1]; + psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); + + psTrace(".psLib.dataManip.psMinimize", 6, + "(mean, stdev) is (%f, %f)\n", mean, stdev); + + if (deriv == NULL) { + deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); + } + + for (i=0;in;i++) { + x = ((psVector *) (coords->data[i]))->data.F32[0]; + out->data.F32[i] = psGaussian(x, mean, stdev, false); + // psTrace(".psLib.dataManip.psMinimize", 6, + // "out->data.F32[%d]: (x, m, s) is (%f, %f, %f) is %f\n", i, x, mean, stdev, out->data.F32[i]); + } + + for (i=0;in;i++) { + x = ((psVector *) (coords->data[i]))->data.F32[0]; + float tmp = (x - mean) * psGaussian(x, mean, stdev, false); + deriv->data.F32[i][0] = tmp / (stdev * stdev); + + tmp = (x - mean) * (x - mean) * + psGaussian(x, mean, stdev, 0); + deriv->data.F32[i][1] = tmp / (stdev * stdev * stdev); + } + + psTrace(".psLib.dataManip.psMinimize", 4, + "---- psMinimizeLMChi2Gauss1D() end ----\n"); + return(out); +} + + /****************************************************************************** psMinimizeLMChi2(): This routine will take an procedure which calculates @@ -429,5 +479,5 @@ 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, +> will have to convert all F32 input vectors to F64 regardless. So, the F64 port might be. *****************************************************************************/ @@ -457,5 +507,5 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2() begin ----\n"); psS32 numData = y->n; @@ -490,8 +540,9 @@ float newChi2 = 0.0; float lamda = 0.00005; - - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + lamda = 0.05; + + psTrace(".psLib.dataManip.psMinimize", 6, "min->maxIter is %d\n", min->maxIter); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "min->tol is %f\n", min->tol); @@ -504,7 +555,7 @@ while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "------------------------------------------------------\n"); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "Iteration %d. Delta is %f\n", min->iter, min->lastDelta); @@ -515,12 +566,16 @@ currValueVec = func(deriv, params, x); for (n=0;ndata.F32[n] * currValueVec->data.F32[n]); + currChi2+= (currValueVec->data.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)\n", n, + currValueVec->data.F32[n], y->data.F32[n]); } for (p=0;pdata.F32[%d] is %f.\n", p, params->data.F32[p]); } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "Current chi-squared is (%f)\n", currChi2); @@ -548,5 +603,5 @@ // XXX: multiple by -1 here? (beta->data.F64[p])*= -1.0; - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]); } @@ -580,9 +635,9 @@ } for (j=0;jdata.F64[j][k]); - } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n"); + psTrace(".psLib.dataManip.psMinimize", 6, "%f ", A->data.F64[j][k]); + } + psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n"); } @@ -597,5 +652,5 @@ // for (i=0;idata.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]); if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { @@ -607,8 +662,8 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + 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]); } @@ -621,9 +676,11 @@ newValueVec = func(deriv, newParams, x); for (n=0;ndata.F32[n] * newValueVec->data.F32[n]); + newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) * + (newValueVec->data.F32[n] - y->data.F32[n]); + } psFree(newValueVec); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2); @@ -640,8 +697,10 @@ } lamda*= 0.1; + psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n"); } else { lamda*= 10.0; - } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n"); + } + psTrace(".psLib.dataManip.psMinimize", 4, "lamda is %f\n", lamda); min->iter++; @@ -662,5 +721,5 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2() end (false) ----\n"); return(false); @@ -942,4 +1001,14 @@ } +#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ +for (int i=0;in;i++) { \ + if (BASEMASK->data.U8[i] == 0) { \ + OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ + } else { \ + OUT->data.F32[i] = BASE->data.F32[i]; \ + } \ +} \ + + /****************************************************************************** p_psDetermineBracket(): This routine takes as input an arbitrary function, @@ -947,18 +1016,26 @@ 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). - -Algorithm: XXX completely ad hoc: start with the user-supplied starting +f(param + b * line) is less than f(param + a * line) and +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. -XXX: new algorithm: start at x=0, expand in one direction until the function +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. -XXX: This is F32 only - *****************************************************************************/ +XXX: +This is F32 only. + +XXX: +output bracket vector should be an input as well. +*****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, psVector *line, @@ -966,4 +1043,155 @@ const psArray *coords, 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; + psS32 iter = 100; + float aDir = 0.0; + float cDir = 0.0; + float new_aDir = 0.0; + float new_cDir = 0.0; + psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); + float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; + psVector *tmp = NULL; + psS32 i = 0; + psS32 null = 0; + + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() begin ----\n"); + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + if (line->data.F32[i] >= FLT_EPSILON) { + null = 1; + } + } + } + + if (null == 0) { + psTrace(".psLib.dataManip.p_psDetermineBracket", 2, + "p_psDetermineBracket() called with zero line vector.\n"); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end (NULL) ----\n"); + psFree(bracket); + return(NULL); + } + + tmp = psVectorAlloc(params->n, PS_TYPE_F32); + + b = 0; + a = -stepSize; + c = stepSize; + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a); + fa = func(tmp, coords); + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b); + fb = func(tmp, coords); + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c); + fc = func(tmp, coords); + + if (fa < fb) { + aDir = -1; + } else { + aDir = 1; + } + + if (fc < fb) { + cDir = -1; + } else { + cDir = 1; + } + + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); + + while (iter > 0) { + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "psDetermineBracket(): iteration %d\n", iter); + if ((fb < fa) && (fb < fc)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + stepSize*= (1.0 + stepSize); + a =- stepSize; + c =+ stepSize; + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]); + } else { + tmp->data.F32[i] = params->data.F32[i]; + } + } + fa = func(tmp, coords); + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]); + } else { + tmp->data.F32[i] = params->data.F32[i]; + } + } + fc = func(tmp, coords); + + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc); + + if (fa < fb) { + new_aDir = -1; + } else { + new_aDir = 1; + } + + if (fc < fb) { + new_cDir = -1; + } else { + new_cDir = 1; + } + if ((new_aDir == 1) && (aDir == -1)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + + if ((new_cDir == 1) && (cDir == -1)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + aDir = new_aDir; + cDir = new_cDir; + iter--; + } + psFree(tmp); + psFree(bracket); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end (NULL) ----\n"); + return(NULL); +} + +psVector *p_psDetermineBracket2(psVector *params, + psVector *line, + const psVector *paramMask, + const psArray *coords, + psMinimizePowellFunc func) { float a = 0.0; @@ -1365,4 +1593,5 @@ float fqp = 0.0; float diff = 0.0; + int iterationNumber = 0; psTrace(".psLib.dataManip.psMinimizePowell", 4, @@ -1407,9 +1636,8 @@ } - min->iter = 0; - while (min->iter < min->maxIter) { - min->iter++; + while (iterationNumber < min->maxIter) { + iterationNumber++; psTrace(".psLib.dataManip.psMinimizePowell", 6, - "psMinimizePowell() iteration %d\n", min->iter); + "psMinimizePowell() iteration %d\n", iterationNumber); // 3: For each dimension in params, move Q only in the vector v[i] to @@ -1418,4 +1646,7 @@ baseFuncVal = func(Q, coords); currFuncVal = baseFuncVal; + psTrace(".psLib.dataManip.psMinimizePowell", 6, + "Current function value is %f\n", currFuncVal); + diff = 0.0; biggestDiff = 0; @@ -1423,5 +1654,5 @@ for (i=0;idata.U8[i] == 0) { - dummyMin.maxIter = 100; + dummyMin.maxIter = 20; dummyMin.tol = 0.01; mul = p_psLineMin(&dummyMin, Q, v[i], myParamMask, coords, func); @@ -1433,5 +1664,6 @@ } } - + psTrace(".psLib.dataManip.psMinimizePowell", 6, + "New function value is %f\n", currFuncVal); // 4: Set the vector u = Q - P for (i=0;iiter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (true) ----\n"); @@ -1513,9 +1745,10 @@ psFree(pQP); psFree(Q); + if (paramMask == NULL) { + psFree(myParamMask); + } + min->iter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (true) ----\n"); - if (paramMask == NULL) { - psFree(myParamMask); - } return(true); } @@ -1532,4 +1765,5 @@ psFree(pQP); psFree(Q); + min->iter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (false) ----\n"); @@ -1549,4 +1783,12 @@ const psArray *coords) { + psTrace(".psLib.dataManip.myPowellChi2Func", 4, + "---- myPowellChi2Func() begin ----\n"); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_EMPTY(params, NAN); + PS_VECTOR_CHECK_NULL(myValue, NAN); + PS_VECTOR_CHECK_EMPTY(myValue, NAN); + PS_PTR_CHECK_NULL(coords, NAN); + float chi2 = 0.0; float d; @@ -1555,9 +1797,18 @@ tmp = Chi2PowellFunc(params, coords); - for (i=0;in;i++) { - d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i]; - chi2+= d * d; + if (myError == NULL) { + for (i=0;in;i++) { + d = (tmp->data.F32[i] - myValue->data.F32[i]); + chi2+= d * d; + } + } else { + for (i=0;in;i++) { + d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i]; + chi2+= d * d; + } } psFree(tmp); + psTrace(".psLib.dataManip.myPowellChi2Func", 4, + "---- myPowellChi2Func() end (%f) ----\n", chi2); return(chi2); } @@ -1584,4 +1835,5 @@ myValue = (psVector *) value; myError = (psVector *) error; + Chi2PowellFunc = func;