Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 2247) +++ trunk/psLib/src/math/psMinimize.c (revision 2250) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.79 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-01 18:57:05 $ + * @version $Revision: 1.80 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-01 23:57:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -572,5 +572,5 @@ (currValueVec->data.F32[n] - y->data.F32[n]); psTrace(".psLib.dataManip.psMinimize", 6, - "data[%d], chi2 calculation+= (%f * %f)\n", n, + "data[%d], chi2 calculation+= (%f - %f)^2\n", n, currValueVec->data.F32[n], y->data.F32[n]); } @@ -1004,4 +1004,6 @@ } +// This macro takes as input the vector BASE and adds a multiple of the vector +// 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++) { \ @@ -1036,7 +1038,7 @@ and the parameter to vary, and the line along which it must vary. This 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). +f(param + b * line) < f(param + a * line) +f(param + b * line) < f(param + c * line) +a < b < c Algorithm: @@ -1192,4 +1194,19 @@ +#define RETURN_FINAL_BRACKET(d) \ +if (a < c) { \ + bracket->data.F32[0] = a; \ + bracket->data.F32[1] = b; \ + bracket->data.F32[2] = c; \ +} else { \ + bracket->data.F32[0] = c; \ + bracket->data.F32[1] = b; \ + bracket->data.F32[2] = a; \ +} \ +psTrace(".psLib.dataManip.p_psDetermineBracket", 4, \ + "---- p_psDetermineBracket() end ----\n"); \ +psTrace(".psLib.dataManip.p_psDetermineBracket", 4, "Final bracket (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); \ +return(bracket); \ + #define PS_DETERMINE_BRACKET_MAX_ITERATIONS 100 psVector *p_psDetermineBracket2(psVector *params, @@ -1206,14 +1223,8 @@ float fc = 0.0; psS32 iter = 0; - psVector *tmp = NULL; + PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32); psBool boolLineIsNull = true; float prevMin = 0.0; int countMin = 0; - - psTrace(".PS_SEG.PS_PWD.PS_FILE.__func__", 6, - "---- __func__() begin COOL"); - - psTrace(".psLib.dataManip.__func__", 4, - "---- p_psDetermineBracket() COOL ----\n"); psTrace(".psLib.dataManip.p_psDetermineBracket", 4, @@ -1230,65 +1241,32 @@ } - psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); - tmp = psVectorAlloc(params->n, PS_TYPE_F32); - a = -0.5; - b = 0.0; - c = 0.5; - - PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b); - fb = func(tmp, coords); - + // We determine in what x-direction does the function decrease. + a = 0.0; + fa = func(params, coords); + b = 0.5; + iter = 0; do { - a*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE); - PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a); - fa = func(tmp, coords); - } while (fabs(fb - fa) < FLT_EPSILON); - - do { - c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE); - PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c); - fc = func(tmp, coords); - } while (fabs(fc - fb) < FLT_EPSILON); - - psTrace(".psLib.dataManip.p_psDetermineBracket", 6, - "initial (fa, Fb, Fc) is (%f, %f, %f)\n", fa, fb, fc); - - 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); - } - - // The object of this code is to set the initial point a to 0.0, and - // ensure that (fb < fa). - if (fc < fb) { - a = b; - b = c; - fa = fb; - fb = fc; - } else if (fa < fb) { - float aSave = a; - float faSave = fa; - + b*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE); + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b); + fb = func(tmp, coords); + } while ((fabs(fb - fa) < FLT_EPSILON) && (iter++ < 100)); + + if (fb > fa) { a = b; fa = fb; - - b = aSave; - fb = faSave; + b = 0.0; + fb = func(params, coords); } c = b; - // At this point we have (a, b) and we know that (fa > fb). Initially, c=b; + // At this point we have (a, b) and we know that (fa >= fb). Initially, c=b; // We keep stretching b out further from "a" until (fc > previous fc). If // that happens, then we have our bracket. + psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); iter = 0; while (iter < PS_DETERMINE_BRACKET_MAX_ITERATIONS) { psTrace(".psLib.dataManip.p_psDetermineBracket", 6, "psDetermineBracket(): iterationA %d\n", iter); - c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE); + c+= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE) * (c - a); PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c); @@ -1296,14 +1274,8 @@ psTrace(".psLib.dataManip.p_psDetermineBracket", 6, - "IterA(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc); + "Iteration(%d) (bracket): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc); 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); + RETURN_FINAL_BRACKET(); } else { b = c; @@ -1312,5 +1284,5 @@ // This code maintains a count of how many times the minimum fc has - // statyed the same. If it gets too high, we exit this loop. + // stayed the same. If it gets too high, we exit this loop. if (fc == prevMin) { countMin++; @@ -1320,11 +1292,5 @@ prevMin = fc; if (countMin == 10) { - 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); + RETURN_FINAL_BRACKET(); } @@ -1332,10 +1298,7 @@ } - // Okay, let's try the other direction. - - psFree(tmp); psFree(bracket); psTrace(".psLib.dataManip.p_psDetermineBracket", 4, - "---- p_psDetermineBracket() end (NULL) ----\n"); + "---- p_psDetermineBracket() end (NULL) (BAD) ----\n"); return(NULL); } @@ -1351,5 +1314,8 @@ function. XXX: This is F32 only +XXX: Since this is an internal function, many of the parameter checks are + redundant. *****************************************************************************/ +#define PS_LINEMIN_MAX_ITERATIONS 30 float p_psLineMin(psMinimization *min, psVector *params, @@ -1359,4 +1325,16 @@ psMinimizePowellFunc func) { + PS_PTR_CHECK_NULL(min, NAN); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_EMPTY(params, NAN); + PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_NULL(line, NAN); + PS_VECTOR_CHECK_EMPTY(line, NAN); + PS_VECTOR_CHECK_TYPE(line, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_NULL(paramMask, NAN); + PS_VECTOR_CHECK_EMPTY(paramMask, NAN); + PS_VECTOR_CHECK_TYPE(paramMask, PS_TYPE_U8, NAN); + PS_PTR_CHECK_NULL(coords, NAN); + PS_PTR_CHECK_NULL(func, NAN); psVector *bracket; float a = 0.0; @@ -1369,52 +1347,25 @@ float fn = 0.0; float mul = 0.0; - psVector *tmpa = NULL; - psVector *tmpb = NULL; - psVector *tmpc = NULL; - psVector *tmpn = NULL; + PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32); + PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32); + PS_VECTOR_GEN_STATIC_RECYCLED(tmpc, params->n, PS_TYPE_F32); + PS_VECTOR_GEN_STATIC_RECYCLED(tmpn, params->n, PS_TYPE_F32); psS32 i = 0; psS32 boolLineIsNull = true; - - psTrace(".psLib.dataManip.p_psLineMin", 4, - "---- p_psLineMin() begin ----\n"); - psTrace(".psLib.dataManip.p_psLineMin", 6, - "min->maxIter is %d\n", min->maxIter); - psTrace(".psLib.dataManip.p_psLineMin", 6, - "min->tol is %f\n", min->tol); - - if (paramMask != NULL) { - for (i=0;in;i++) { - if (paramMask->data.U8[i] == 0) { - if (line->data.F32[i] >= FLT_EPSILON) { - boolLineIsNull = false; - psTrace(".psLib.dataManip.p_psLineMin", 4, - "line->data.F32[%d] is %f\n", i, line->data.F32[i]); - } - } - } - } else { - for (i=0;in;i++) { - if (line->data.F32[i] >= FLT_EPSILON) { - boolLineIsNull = false; - psTrace(".psLib.dataManip.p_psLineMin", 4, - "line->data.F32[%d] is %f\n", i, line->data.F32[i]); - } - } - } + psS32 numIterations = 0; + + psTrace(".psLib.dataManip.p_psLineMin", 4, "---- p_psLineMin() begin ----\n"); + PS_VECTOR_F32_CHECK_ZERO_VECTOR(line, boolLineIsNull); if (boolLineIsNull == true) { + min->value = func(params, coords); psTrace(".psLib.dataManip.p_psLineMin", 2, - "p_psLineMin() called with zero line vector.\n"); + "p_psLineMin() called with zero line vector. Return 0.0. Function valus is %f\n", min->value); return(0.0); } - - tmpa = psVectorAlloc(params->n, PS_TYPE_F32); - tmpb = psVectorAlloc(params->n, PS_TYPE_F32); - tmpc = psVectorAlloc(params->n, PS_TYPE_F32); - tmpn = psVectorAlloc(params->n, PS_TYPE_F32); for (i=0;in;i++) { psTrace(".psLib.dataManip.p_psLineMin", 6, - "params, paramMask, line [%d] is (%f %d %f)\n", i, + "(params, paramMask, line)[%d] is (%f %d %f)\n", i, params->data.F32[i], paramMask->data.U8[i], @@ -1424,39 +1375,24 @@ bracket = p_psDetermineBracket2(params, line, paramMask, coords, func); if (bracket == NULL) { - psFree(tmpa); - psFree(tmpb); - psFree(tmpc); - psFree(tmpn); psError(__func__, "(1) Could not bracket minimum."); return(NAN); } - - min->iter = 0; - while (min->iter < min->maxIter) { - min->iter++; + numIterations = 0; + while (numIterations < PS_LINEMIN_MAX_ITERATIONS) { + numIterations++; psTrace(".psLib.dataManip.p_psLineMin", 6, - "p_psLineMin(): iteration %d\n", min->iter); + "p_psLineMin(): iteration %d\n", numIterations); a = bracket->data.F32[0]; b = bracket->data.F32[1]; c = bracket->data.F32[2]; - - for (i=0;in;i++) { - if (paramMask->data.U8[i] == 0) { - tmpa->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]); - tmpb->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]); - tmpc->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]); - } else { - tmpa->data.F32[i] = params->data.F32[i]; - tmpb->data.F32[i] = params->data.F32[i]; - tmpc->data.F32[i] = params->data.F32[i]; - } - } + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpa, a); + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpb, b); + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpc, c); fa = func(tmpa, coords); fb = func(tmpb, coords); fc = func(tmpc, coords); psTrace(".psLib.dataManip.p_psLineMin", 6, - "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, - a, b, c, fa, fb, fc); + "LineMin: f(%f %f %f) is (%f %f %f)\n", a, b, c, fa, fb, fc); // We determine which is the biggest segment in [a,b,c] then split @@ -1495,34 +1431,25 @@ } psTrace(".psLib.dataManip.p_psLineMin", 6, - "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]); + "LineMin: new bracket is (%f %f %f)\n", bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]); mul = bracket->data.F32[1]; if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) { - for (i=0;in;i++) { - if (paramMask->data.U8[i] == 0) { - params->data.F32[i]+= mul * line->data.F32[i]; - } - } + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul); min->value = func(params, coords); psFree(bracket); - psFree(tmpa); - psFree(tmpb); - psFree(tmpc); - psFree(tmpn); psTrace(".psLib.dataManip.p_psLineMin", 4, - "---- p_psLineMin() end.a (%f) ----\n", mul); + "---- p_psLineMin() end.a (%f) (%f) ----\n", mul, min->value); return(mul); } } + mul = bracket->data.F32[1]; + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul); + min->value = func(params, coords); + psTrace(".psLib.dataManip.p_psLineMin", 4, + "---- p_psLineMin() end.b (%f) %f ----\n", mul, min->value); + psFree(bracket); - psFree(tmpa); - psFree(tmpb); - psFree(tmpc); - psFree(tmpn); - - psTrace(".psLib.dataManip.p_psLineMin", 4, - "---- p_psLineMin() end.b (0.0) ----\n"); - return(0.0); + return(mul); } @@ -1560,10 +1487,10 @@ psS32 numDims = params->n; + PS_VECTOR_GEN_STATIC_RECYCLED(pQP, numDims, PS_TYPE_F32); + PS_VECTOR_GEN_STATIC_RECYCLED(u, numDims, PS_TYPE_F32); + PS_VECTOR_GEN_STATIC_RECYCLED(Q, numDims, PS_TYPE_F32); psS32 i = 0; psS32 j = 0; psVector **v = NULL; - psVector *pQP = NULL; - psVector *Q = NULL; - psVector *u = NULL; psVector *myParamMask = NULL; // XXX: Use a recylced vector here. psMinimization dummyMin; @@ -1586,10 +1513,8 @@ "min->tol is %f\n", min->tol); - pQP = psVectorAlloc(numDims, PS_TYPE_F32); - u = psVectorAlloc(numDims, PS_TYPE_F32); - Q = psVectorAlloc(numDims, PS_TYPE_F32); - if (paramMask == NULL) { - myParamMask = psVectorAlloc(params->n, PS_TYPE_U8); + myParamMask = psVectorRecycle(myParamMask, params->n, PS_TYPE_U8); + p_psMemSetPersistent(myParamMask, true); + p_psMemSetPersistent(myParamMask->data.V, true); for (i=0;in;i++) { myParamMask->data.U8[i] = 0; @@ -1643,16 +1568,13 @@ if (isnan(mul)) { psError(__func__, "Could not perform line minimization (1).\n"); - psFree(pQP); - psFree(u); - psFree(Q); - if (paramMask == NULL) { - psFree(myParamMask); - } - for (i=0;i biggestDiff) { biggestDiff = fabs(dummyMin.value - currFuncVal); @@ -1664,4 +1586,5 @@ psTrace(".psLib.dataManip.psMinimizePowell", 6, "New function value is %f\n", currFuncVal); + // 4: Set the vector u = Q - P for (i=0;idata.F32[i], + "u[i]=Q[i]-P[i] (%f = %f - %f)\n", u->data.F32[i], Q->data.F32[i], params->data.F32[i]); @@ -1688,10 +1611,4 @@ if (isnan(mul)) { psError(__func__, "Could not perform line minimization. (2)\n"); - psFree(pQP); - psFree(u); - psFree(Q); - if (paramMask == NULL) { - psFree(myParamMask); - } for (i=0;iiter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, - "---- psMinimizePowell() end (true) ----\n"); + "---- psMinimizePowell() end (1)(true) ----\n"); return(true); } @@ -1754,31 +1665,44 @@ } psFree(v); - psFree(u); - psFree(pQP); - psFree(Q); - if (paramMask == NULL) { - psFree(myParamMask); - } min->iter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, - "---- psMinimizePowell() end (true) ----\n"); + "---- psMinimizePowell() end (2) (true) ----\n"); return(true); } } - if (paramMask == NULL) { - psFree(myParamMask); - } for (i=0;iiter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, - "---- psMinimizePowell() end (false) ----\n"); + "---- psMinimizePowell() end (0) (false) ----\n"); return(false); +} + + +/****************************************************************************** +XXX: We assume unnormalized gaussians. +XXX: Currently, yErr is ignored. + *****************************************************************************/ +psVector *psMinimizePowellChi2Gauss1D(const psVector *params, + const psArray *coords) +{ + 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); + + for (i=0;in;i++) { + x = ((psVector *) (coords->data[i]))->data.F32[0]; + out->data.F32[i] = psGaussian(x, mean, stdev, false); + } + + return(out); } @@ -1822,5 +1746,5 @@ psFree(tmp); psTrace(".psLib.dataManip.myPowellChi2Func", 4, - "---- myPowellChi2Func() end (%f) ----\n", chi2); + "---- myPowellChi2Func() end (chi2 is %f) ----\n", chi2); return(chi2); }