Changeset 29060 for branches/sc_branches/trunkTest/psLib/src/math/psMinimizeLMM.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• psLib/src/math/psMinimizeLMM.c (modified) (9 diffs)

branches/sc_branches/trunkTest/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);