Changeset 14320 for trunk/psLib/src/math/psMinimizeLMM.c

Timestamp:

Jul 19, 2007, 2:22:24 PM (19 years ago)

Author:

magnier

Message:

dLinear, the linearized change in the chisq, was being calculated
AFTER the parameters were set to their limits. this is inconsistent
with the assumptions, and resulted in impossible negative values for
dLinear, which in turn caused problems for the fitting convergence.

I moved the call to dLinear into GuessABP, and modified the APIs a bit
to deal with this.

File:

• trunk/psLib/src/math/psMinimizeLMM.c (modified) (12 diffs)

trunk/psLib/src/math/psMinimizeLMM.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.31 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-06-26 19:27:04 $
+ *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-20 00:22:24 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -63 +63 @@
     const psVector *paramMask,
     psMinimizeLMLimitFunc checkLimits,
-    psF32 lambda)
+    psF32 lambda,
+    psF32 *dLinear)
 {
     PS_ASSERT_VECTOR_TYPE(Alpha,     PS_TYPE_F32,  false);
@@ -135 +136 @@
     # endif
 
+    // measure linear model prediction
+    // (we must do this before truncating Beta below)
+    if (dLinear) {
+        *dLinear = psMinLM_dLinear(Beta, beta, paramMask, lambda);
+    }
+
     // apply Beta to get new Params values
     for (int j = 0; j < params->n; j++) {
@@ -196 +203 @@
     }
 
-    bool rcBool = psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, 0.0);
+    bool rcBool = psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL);
     if (rcBool == false) {
         psTrace ("psLib.math", 5, "psMinLM_GuessABP() returned FALSE.\n");
@@ -223 +230 @@
     const psVector *Beta,
     const psVector *beta,
+    const psVector *paramMask,
     psF32 lambda)
 {
@@ -230 +238 @@
     psF32 *B = Beta->data.F32;
     psF32 *b = beta->data.F32;
+
+    float dh = 0.0, sh = 0.0, Sh = 0.0;
     for (int i = 0; i < beta->n; i++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[i])) continue;
+        dh = lambda*B[i] + b[i];
+        sh = 0.5*B[i]*dh;
+        Sh += sh;
+        // fprintf (stderr, "%f, %f  * %f %f %f * :  %f : %f + %f\n", B[i], b[i], dh, sh, Sh, lambda*PS_SQR(B[i]) + B[i]*b[i], lambda*PS_SQR(B[i]), B[i]*b[i]);
         dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i];
     }
@@ -391 +406 @@
     psF32 Chisq = 0.0;
     psF32 lambda = 0.001;
+    psF32 dLinear = 0.0;
 
     // the user provides the error or NULL.  we need to convert
@@ -422 +438 @@
 
         // set a new guess for Alpha, Beta, Params
-        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)) {
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
             min->iter ++;
             lambda *= 10.0;
@@ -428 +444 @@
         }
 
-        // measure linear model prediction
-        psF32 dLinear = psMinLM_dLinear(Beta, beta, lambda);
-
         // dump some useful info if trace is defined
         if (psTraceGetLevel("psLib.math") >= 6) {
-            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)");
-            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)");
+            p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)");
         }
         if (psTraceGetLevel("psLib.math") >= 5) {
@@ -454 +468 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
-                Chisq, min->lastDelta, rho);
+        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);
 
         // dump some useful info if trace is defined
@@ -470 +484 @@
             beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
             params = psVectorCopy(params, Params, PS_TYPE_F32);
-            lambda *= 0.1;
+            lambda *= 0.25;
         } else {
             lambda *= 10.0;
@@ -480 +494 @@
     // construct & return the covariance matrix (if requested)
     if (covar != NULL) {
-        if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0)) {
+        if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
             psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
         }