Changeset 35617

Timestamp:

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

Author:

eugene

Message:

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

Location:

branches/eam_branches/ipp-20130509/psLib

Files:

• src/math/psMinimizeLMM.c (modified) (4 diffs)
• src/math/psMinimizeLMM.h (modified) (1 diff)
• test/math/Makefile.am (modified) (1 diff)
• test/pstap/src/pstap.h (modified) (1 diff)

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);

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

@@ -134 +134 @@
 );
 
+/** Minimizes a specified function based on the Levenberg-Marquardt method
+    Uses alternative convergence criterion.
+ *
+ *  @return bool:   True if successful.
+ */
+bool psMinimizeLMChi2_Alt(
+    psMinimization *min,               ///< Minimization specification
+    psImage *covar,                    ///< Covariance matrix
+    psVector *params,                  ///< "Best Guess" for the parameters that minimize func
+    psMinConstraint *constraint, ///< Constraints on the parameters
+    const psArray *x,                  ///< Measurement ordinates of multiple vectors
+    const psVector *y,                 ///< Measurement coordinates
+    const psVector *yWt,               ///< Errors in the measurement coordinates
+    psMinimizeLMChi2Func func          ///< Specified function
+);
+
 bool psMinimizeGaussNewtonDelta (
     psVector *delta,

branches/eam_branches/ipp-20130509/psLib/test/math/Makefile.am

@@ -48 +48 @@
         tap_psMatrixVectorArithmetic04 \
         tap_psMinimizePowell \
+        tap_psMinimizeLMM \
+        tap_psMinimizeLMM_Alt \
+        tap_psMinimizeLMM_Trail \
         tap_psSpline1D \
         tap_psPolynomialMD \

branches/eam_branches/ipp-20130509/psLib/test/pstap/src/pstap.h

@@ -144 +144 @@
     } \
 }
+
+# define ok_float      is_float     
+# define ok_float_tol  is_float_tol 
+# define ok_double     is_double    
+# define ok_double_tol is_double_tol
+# define ok_str        is_str         
+# define ok_strn       is_strn      
+# define ok_int        is_int         
+# define ok_long       is_long      
+# define ok_bool       is_bool