Changeset 35637 for branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.c

Timestamp:

Jun 9, 2013, 3:27:54 PM (13 years ago)

Author:

eugene

Message:

choice of convergence method via user selections

File:

• branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.c (modified) (13 diffs)

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

@@ -680 +680 @@
     }
 
-    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
-    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
+    int nFitParams = Beta->n;
+    psImage *alpha   = psImageAlloc(nFitParams, nFitParams, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(nFitParams, PS_TYPE_F32);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
 
@@ -687 +688 @@
     psF32 Chisq = 0.0;
     psF32 lambda = 0.001;
+    psF32 dLinear = 0.0;
     psF32 nu = 2.0;
-    psF32 dLinear = 0.0;
 
     // the user provides the error or NULL.  we need to convert
@@ -700 +701 @@
 
     // number of degrees of freedom for this fit
-    int nDOF = dy->n - params->n;
+    int nDOF = dy->n - nFitParams;
 
     // calculate initial alpha and beta, set chisq (min->value)
@@ -722 +723 @@
     }
 
-    // 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.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+
+        if (min->chisqConvergence) {
+          psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        } else {
+          psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*nFitParams, min->maxTol*nFitParams);
+        }
 
         // set a new guess for Alpha, Beta, Params
@@ -732 +737 @@
             min->iter ++;
             if (min->iter >=  min->maxIter) break;
-            // lambda *= 10.0;
-            lambda *= 2.0;
+            lambda *= 10.0;
+            // ALT? lambda *= 2.0;
             continue;
         }
@@ -754 +759 @@
         // Beta is the actual parameter change for this pass
 
-        // note that Alpha & Beta only represent unmasked parameters
+        // note that Alpha & Beta only represent unmasked parameters, while params and Params have all 
 
         // dParSigma = Alpha[i][i] : error (squared) on parameter i
-        // dParDelta = Beta[i]^2   : change (squared) in parameter i
+        // dParDelta = Params->data.F32[i] - params->data.F32[i]     : change on parameter i
         float rParSigma = 0.0;
-        // for (int j = 0; j < Beta->n; j++) {
-        //     rParSigma += PS_SQR(Beta->data.F32[j]) / Alpha->data.F32[j][j];
-        // }
         for (int j = 0, J = 0; j < Params->n; j++) {
             if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
                 continue;
             }
-            // fprintf (stderr, "del: %f  - %f  vs %f  : %f\n", Params->data.F32[j], params->data.F32[j], sqrt(Alpha->data.F32[J][J]), PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J]);
             rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J];
             J++;
@@ -779 +780 @@
             min->iter ++;
             if (min->iter >= min->maxIter) break;
-            // lambda *= 10.0;
-            lambda *= 2.0;
+            lambda *= 10.0;
+            // ALT lambda *= 2.0;
             continue;
         }
@@ -807 +808 @@
         // XXX the old version of lambda changes:
         // XXX : Madsen gives suggestion for better use of rho
-# if (0)
-        if (rho >= -1e-6) {
-            min->value = Chisq;
-            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
-            beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
-            params = psVectorCopy(params, Params, PS_TYPE_F32);
-            lambda *= 0.25;
-        } else {
-            lambda *= 10.0;
-        }
-# endif
-
         // rho is positive if the new chisq is smaller
         if (rho >= -1e-6) {
@@ -826 +815 @@
             params = psVectorCopy(params, Params, PS_TYPE_F32);
         } 
-
-# if (0)
-        // adjust the gain ratio (lambda) based on rho
-        if (rho < 0.25) {
-            lambda *= 2.0;
-        } 
-        if (rho > 0.75) {
-            lambda *= 0.333;
-        }
-# else
-        if (rho > 0.0) {
-          lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
-          nu = 2.0;
-        } else {
-          lambda *= nu;
-          nu *= 2.0;
+        switch (min->gainFactorMode) {
+          case 0:
+            if (rho >= -1e-6) {
+              lambda *= 0.25;
+            } else {
+              lambda *= 10.0;
+            }
+            break;
+
+          case 1:
+            // adjust the gain ratio (lambda) based on rho
+            if (rho < 0.25) {
+              lambda *= 2.0;
+            } 
+            if (rho > 0.75) {
+              lambda *= 0.333;
+            }
+            break;
+
+          case 2:
+            if (rho > 0.0) {
+              lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
+              nu = 2.0;
+            } else {
+              lambda *= nu;
+              nu *= 2.0;
+            }
+            break;
         }
-# endif
         min->iter++;
 
@@ -850 +850 @@
         done = (min->iter >= min->maxIter);
         
+        // 2) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change)
+        if (min->lastDelta < -1e-6) {
+            continue;
+        }
+
+        // save this value in case we stop iterating
+        min->rParSigma = rParSigma;
+
         // 2) require chisqDOF < maxChisqDOF (if maxChisqDOF is not NAN)
+        // keep iterating regardless of rParSigma in this case
         float chisqDOF = Chisq / nDOF;
         if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) {
@@ -856 +865 @@
         }
 
-        // 3) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change)
-        if (min->lastDelta < -1e-6) {
-            continue;
+        // delta-chisq or rParSigma ?
+        if (min->chisqConvergence) {
+          done |= (min->lastDelta < min->minTol);
+        } else {
+          done |= (rParSigma < min->minTol*nFitParams);
         }
-
-        // 4) require rParDelta < \eta * rParSigma
-        done |= (rParSigma < 0.01);
-    }
-    psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+    }
+    psTrace("psLib.math", 5, "chisq: %f, last delta: %f, rParSigma: %f, Niter: %d\n", min->value, min->lastDelta, min->rParSigma, min->iter);
 
     // construct & return the covariance matrix (if requested)
@@ -898 +906 @@
 
     // if the last improvement was at least as good as maxTol, accept the fit:
-    // XXX : maybe re-work this to use rParSigma & maxTol
-    if (min->lastDelta <= min->maxTol) {
+    if (min->chisqConvergence) {
+      if (min->lastDelta <= min->maxTol) {
         psTrace("psLib.math", 6, "---- end (true) ----\n");
         return(true);
+      }
+    } else {
+      if (min->rParSigma <= min->maxTol*nFitParams) {
+        psTrace("psLib.math", 6, "---- end (true) ----\n");
+        return(true);
+      }
     }
     psTrace("psLib.math", 6, "---- end (false) ----\n");
@@ -952 +966 @@
     min->iter = 0;
     min->lastDelta = NAN;
+    min->rParSigma = NAN;
     min->maxChisqDOF = NAN;
 
+    // we default to the old algorithm for convergence
+    min->chisqConvergence = true;
+    min->gainFactorMode = 0;
+    
     return(min);
 }