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

Timestamp:

Jun 3, 2013, 8:40:02 AM (13 years ago)

Author:

eugene

Message:

new convergence criterion for LMM (in _Alt only for now)

File:

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

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

@@ -325 +325 @@
         chisq += PS_SQR(delta) * dy->data.F32[i];
 
+        // XXX remove this later:
+        // psVector *tmp = x->data[i];
+        // fprintf (stderr, "%f %f  %f %f  %f\n", tmp->data.F32[0], tmp->data.F32[1], y->data.F32[i], dy->data.F32[i], ymodel);
+
         if (isnan(dy->data.F32[i])) goto escape;
         if (isnan(delta)) goto escape;
@@ -716 +720 @@
         }
 
+        // calculate the parameter change (rParDelta) and error radius (rParSigma)
+        //    rParDelta : radius of parameter change;
+        //    rParSigma : radius of parameter error 
+        
+        // note that (before SetABX) Alpha[i][i] is the covariance matrix and
+        // Beta is the actual parameter change for this pass
+
+        // note that Alpha & Beta only represent unmasked parameters
+
+        // dParSigma = Alpha[i][i] : error (squared) on parameter i
+        // dParDelta = Beta[i]^2   : change (squared) in 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];
+        }
+        rParSigma = sqrt(rParSigma);
+        psTrace("psLib.math", 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
+
         // calculate Chisq for new guess, update Alpha & Beta
         Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func);
@@ -744 +766 @@
         min->lastDelta = (min->value - Chisq) / nDOF;
 
-        /* rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho) */
+        // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho)
+        // XXX : Madsen gives suggestion for better use of rho
         if (rho >= -1e-6) {
             min->value = Chisq;
@@ -766 +789 @@
         }
 
-        // 3) require deltaChi > 1e-6
-        if (min->lastDelta < 1e-5) {
+        // 3) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change)
+        if (min->lastDelta < -1e-6) {
             continue;
         }
 
         // 4) require rParDelta < \eta * rParSigma
-        //    rParDelta : radius of parameter change;
-        //    rParSigma : radius of parameter error 
-        
-        // note that alpha & beta only represent unmasked parameters
-
-        // R(Par_sigma) = \sum (alpha[i][i])
-        float rParSigma = 0.0;
-        for (int j = 0; j < beta->n; j++) {
-            rParSigma += Alpha->data.F32[j][j];
-        }
-        rParSigma = sqrt(rParSigma);
-
-        // note that beta or Beta is the actual parameter change for this pass : 
-        // new param = old param - beta.  I need to be sure I'm using the new 
-        // or the old one and also the one before or after the matrix eqn Ax = B is solved.
-        float rParDelta = 0.0;
-        for (int j = 0; j < beta->n; j++) {
-            rParDelta += PS_SQR(lastBeta->data.F32[j] - Beta->data.F32[j]);
-        }
-        rParDelta = sqrt(rParDelta);
-
-        psTrace("psLib.math", 5, "rParDelta : %f, rParSigma: %f, Niter: %d\n", rParDelta, rParSigma, min->iter);
+        done |= (rParSigma < 0.01);
     }
     psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);