Changeset 28655

Timestamp:

Jul 11, 2010, 3:18:21 PM (16 years ago)

Author:

eugene

Message:

add maxTol and maxChisqDOF to psMin to allow for better control over fit success and convergence

Location:

branches/eam_branches/ipp-20100621/psLib/src/math

Files:

• psMinimizeLMM.c (modified) (9 diffs)
• psMinimizeLMM.h (modified) (3 diffs)
• psMinimizePowell.c (modified) (4 diffs)
• psStats.c (modified) (1 diff)

branches/eam_branches/ipp-20100621/psLib/src/math/psMinimizeLMM.c

@@ -441 +441 @@
     }
 
+    // number of degrees of freedom for this fit
+    int nDOF = dy->n - params->n;
+
     // calculate initial alpha and beta, set chisq (min->value)
     min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
@@ -456 +459 @@
     }
 
-    // iterate until the tolerance is reached, or give up
-    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
+    // iterate until: (a) nIter = min->iter or (b) (chisq / ndof) < maxChisq and deltaChisq < minTol (but don't stop unless Chisq is finite)
+    bool done = (min->iter >= min->maxIter);
+    while (!done) {
         psTrace("psLib.math", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psLib.math", 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
+        psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
 
         // set a new guess for Alpha, Beta, Params
         if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
             min->iter ++;
+            if (min->iter >=  min->maxIter) break;
             lambda *= 10.0;
             continue;
@@ -482 +487 @@
         if (isnan(Chisq)) {
             min->iter ++;
+            if (min->iter >= min->maxIter) break;
             lambda *= 10.0;
             continue;
@@ -492 +498 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f\n", min->value,
-                Chisq, min->lastDelta, dLinear, rho, lambda);
-
-        psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f\n", min->value,
-                Chisq, min->lastDelta, dLinear, rho, lambda);
+        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF);
+
+        psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda);
 
         // dump some useful info if trace is defined
@@ -506 +510 @@
         /* rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho) */
         if (rho >= -1e-6) {
-            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
+            min->lastDelta = (min->value - Chisq) / nDOF;
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
@@ -516 +520 @@
         }
         min->iter++;
+
+        done = (min->iter >= min->maxIter);
+        
+        // check for convergence:
+        float chisqDOF = Chisq / nDOF;
+        if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) {
+            done |= (min->lastDelta < min->minTol);
+        }
     }
     psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
@@ -549 +561 @@
         psFree(dy);
     }
-    if (min->iter == min->maxIter) {
-        psTrace("psLib.math", 3, "---- end (false) ----\n");
-        return(false);
-    }
-    psTrace("psLib.math", 3, "---- end (true) ----\n");
-    return(true);
+
+    // if the last improvement was at least as good as maxTol, accept the fit:
+    if (min->lastDelta <= min->maxTol) {
+        psTrace("psLib.math", 6, "---- end (true) ----\n");
+        return(true);
+    }
+    psTrace("psLib.math", 6, "---- end (false) ----\n");
+    return(false);
 }
 
@@ -588 +602 @@
 }
 
-psMinimization *psMinimizationAlloc(int maxIter,
-                                    float tol)
+psMinimization *psMinimizationAlloc(int maxIter, float minTol, float maxTol)
 {
     PS_ASSERT_INT_NONNEGATIVE(maxIter, NULL);
@@ -595 +608 @@
     psMinimization *min = psAlloc(sizeof(psMinimization));
     psMemSetDeallocator(min, (psFreeFunc)minimizationFree);
+
     P_PSMINIMIZATION_SET_MAXITER(min,maxIter);
-    P_PSMINIMIZATION_SET_TOL(min,tol);
+    P_PSMINIMIZATION_SET_MIN_TOL(min,minTol);
+    P_PSMINIMIZATION_SET_MAX_TOL(min,maxTol);
+
     min->value = 0.0;
     min->iter = 0;
     min->lastDelta = NAN;
+    min->maxChisqDOF = NAN;
 
     return(min);

branches/eam_branches/ipp-20100621/psLib/src/math/psMinimizeLMM.h

@@ -33 +33 @@
 #define PS_MAX_LMM_ITERATIONS 100
 #define PS_MAX_MINIMIZE_ITERATIONS 100
-#define P_PSMINIMIZATION_SET_MAXITER(m,val) *(int*)&m->maxIter = val
-        #define P_PSMINIMIZATION_SET_TOL(m,val) *(float*)&m->tol = val
+#define P_PSMINIMIZATION_SET_MAXITER(m,val) { *(int*)&m->maxIter  = val; }
+#define P_PSMINIMIZATION_SET_MIN_TOL(m,val) { *(float*)&m->minTol = val; }
+#define P_PSMINIMIZATION_SET_MAX_TOL(m,val) { *(float*)&m->maxTol = val; }
 
                 typedef enum {
@@ -75 +76 @@
 typedef struct
 {
-    const int maxIter;                 ///< Convergence limit
-    const float tol;                   ///< Error Tolerance
+    const int maxIter;                  ///< Convergence limit
+    const float minTol;                 ///< Convergence Tolerance (stop if we reach this value)
+    const float maxTol;                 ///< Max Tolerance (accept fit if last improvement was this good)
     float value;                       ///< Value of function at minimum
     int iter;                          ///< Number of iterations to date
     float lastDelta;                   ///< The last difference for the fit
+    float maxChisqDOF;                 ///< for Chisq minimization, require that we reach here before checking tolerance
 }
 psMinimization;
@@ -102 +105 @@
 psMinimization *psMinimizationAlloc(
     int maxIter,                       ///< Number of minimization iterations to perform.
-    float tol                          ///< Requested error tolerance
+    float minTol,                      ///< stop if tolerance is less than this
+    float maxTol                       ///< accept fit if tolerance is less than this
 ) PS_ATTR_MALLOC;
 

branches/eam_branches/ipp-20100621/psLib/src/math/psMinimizePowell.c

@@ -457 +457 @@
 
         mul = bracket->data.F32[1];
-        if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) {
+        if ((fabs(a-b) < min->minTol) && (fabs(b-c) < min->minTol)) {
             PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
             min->value = func(params, coords);
@@ -524 +524 @@
     psTrace("psLib.math", 4, "---- psMinimizePowell() begin ----\n");
     psTrace("psLib.math", 6, "min->maxIter is %d\n", min->maxIter);
-    psTrace("psLib.math", 6, "min->tol is %f\n", min->tol);
+    psTrace("psLib.math", 6, "min->minTol is %f\n", min->minTol);
 
     if (paramMask == NULL) {
@@ -573 +573 @@
         for (i=0;i<numDims;i++) {
             if (myParamMask->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) {
+
                 P_PSMINIMIZATION_SET_MAXITER((&dummyMin),PS_MINIMIZE_POWELL_LINEMIN_MAX_ITERATIONS);
-                *(float*)&dummyMin.tol = PS_MINIMIZE_POWELL_LINEMIN_ERROR_TOLERANCE;
+                P_PSMINIMIZATION_SET_MIN_TOL((&dummyMin),PS_MINIMIZE_POWELL_LINEMIN_ERROR_TOLERANCE);
+
                 mul = LineMin(&dummyMin, Q, ((psVector *) v->data[i]),
                               myParamMask, coords, func);
@@ -677 +679 @@
 
         // 8: Go to step 3 until the change is less than some tolerance.
-        if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
+        if (fabs(baseFuncVal - currFuncVal) <= min->minTol) {
             psFree(v);
             psFree(pQP);

branches/eam_branches/ipp-20100621/psLib/src/math/psStats.c

@@ -1211 +1211 @@
 
         // Fit a Gaussian to the data.
-        psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
+        psMinimization *minimizer = psMinimizationAlloc(100, 0.01, 1.0); // The minimizer information
         psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
         // Initial guess for the mean (index 0) and var (index 1).