Changeset 35623

Timestamp:

Jun 3, 2013, 2:05:58 PM (13 years ago)

Author:

eugene

Message:

use Madson option for gain factor; build-time choice of _Old vs _Alt for LMM

Location:

branches/eam_branches/ipp-20130509/psLib

Files:

• src/math/psMinimizeLMM.c (modified) (11 diffs)
• src/math/psMinimizeLMM.h (modified) (1 diff)
• test/math/tap_psMinimizeLMM_Trail.c (modified) (1 diff)

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

@@ -365 +365 @@
 
 /******************************************************************************
+psMinimizeLMChi2():  wrapper to call either _Old or _Alt
+  *****************************************************************************/
+bool psMinimizeLMChi2(
+    psMinimization *min,
+    psImage *covar,
+    psVector *params,
+    psMinConstraint *constraint,
+    const psArray *x,
+    const psVector *y,
+    const psVector *yWt,
+    psMinimizeLMChi2Func func)
+{
+    bool status = psMinimizeLMChi2_Alt(
+        min,
+        covar,
+        params,
+        constraint,
+        x,
+        y,
+        yWt,
+        func);
+    return status;
+}
+
+/******************************************************************************
 psMinimizeLMChi2():  This routine will take an procedure which calculates an
 arbitrary function and it's derivative and minimize the chi-squared match
@@ -373 +398 @@
 XXX Make an F64 version?
   *****************************************************************************/
-bool psMinimizeLMChi2(
+bool psMinimizeLMChi2_Old(
     psMinimization *min,
     psImage *covar,
@@ -510 +535 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        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);
+        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, 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: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda);
 
         // dump some useful info if trace is defined
@@ -662 +687 @@
     psF32 Chisq = 0.0;
     psF32 lambda = 0.001;
+    psF32 nu = 2.0;
     psF32 dLinear = 0.0;
 
@@ -706 +732 @@
             min->iter ++;
             if (min->iter >=  min->maxIter) break;
-            lambda *= 10.0;
+            // lambda *= 10.0;
+            lambda *= 2.0;
             continue;
         }
@@ -732 +759 @@
         // 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];
+        // 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++;
         }
         rParSigma = sqrt(rParSigma);
         psTrace("psLib.math", 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
+        fprintf (stderr, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
 
         // calculate Chisq for new guess, update Alpha & Beta
@@ -743 +779 @@
             min->iter ++;
             if (min->iter >= min->maxIter) break;
-            lambda *= 10.0;
+            // lambda *= 10.0;
+            lambda *= 2.0;
             continue;
         }
@@ -753 +790 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        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);
+        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, 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: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda);
 
         // dump some useful info if trace is defined
@@ -767 +804 @@
 
         // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho)
+
+        // 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;
@@ -777 +817 @@
             lambda *= 10.0;
         }
+# endif
+
+        // rho is positive if the new chisq is smaller
+        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);
+        } 
+
+# 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;
+        }
+# endif
         min->iter++;
 
@@ -831 +898 @@
 
     // 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) {
         psTrace("psLib.math", 6, "---- end (true) ----\n");

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

@@ -134 +134 @@
 );
 
+/** Minimizes a specified function based on the Levenberg-Marquardt method.
+ *
+ *  @return bool:   True if successful.
+ */
+bool psMinimizeLMChi2_Old(
+    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
+);
+
 /** Minimizes a specified function based on the Levenberg-Marquardt method
     Uses alternative convergence criterion.

branches/eam_branches/ipp-20130509/psLib/test/math/tap_psMinimizeLMM_Trail.c

@@ -501 +501 @@
 
 # if (1)        
-        psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
+        // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
         bool tmpBool = psMinimizeLMChi2_Alt(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc);
 # else