Changeset 35768 for trunk/psModules/src/objects/pmPCM_MinimizeChisq.c

Timestamp:

Jul 3, 2013, 2:37:22 PM (13 years ago)

Author:

eugene

Message:

deprecate KiiOpen,KiiClose (now KapaOpen,etc); major rework of psEllipse translations : use common functions pmModelAxesToParams and pmModelParamsToAxes ; use new convergence method in pmPCM_MinimizeChisq; add convergence crerition options to psMinimization; threaded versions of pmPSFtryFitEXT and pmPSFtryFitPSF

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/pmPCM_MinimizeChisq.c (modified) (6 diffs)

trunk/psModules/src/objects/pmPCM_MinimizeChisq.c

@@ -81 +81 @@
         psAbort ("programming error: no unmasked parameters to be fit\n");
     }
+    psAssert (pcm->nPar == Beta->n, "did we set the masked parameters correctly??");
 
     // allocate internal arrays (current vs Guess)
-    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
-    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
+    psImage *alpha   = psImageAlloc(pcm->nPar, pcm->nPar, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(pcm->nPar, PS_TYPE_F32);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
 
@@ -90 +91 @@
     psF32 lambda = 0.001;
     psF32 dLinear = 0.0;
+    psF32 nu = 2.0;
 
 # if (USE_FFT && PRE_CONVOLVE)
@@ -120 +122 @@
     bool done = (min->iter >= min->maxIter);
     while (!done) {
-        psTrace("psModules.objects", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psModules.objects", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        psTrace(FACILITY, 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+
+        if (min->chisqConvergence) {
+            psTrace(FACILITY, 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        } else {
+            psTrace(FACILITY, 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*pcm->nPar, min->maxTol*pcm->nPar);
+        }
 
         // set a new guess for Alpha, Beta, Params
@@ -140 +147 @@
             p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
         }
+
+        // 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, while params and Params have all 
+
+        // dParSigma = Alpha[i][i] : error (squared) on parameter i
+        // dParDelta = Params->data.F32[i] - params->data.F32[i]     : change on parameter i
+        float rParSigma = 0.0;
+        for (int j = 0, J = 0; j < Params->n; j++) {
+            if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
+                continue;
+            }
+            rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J];
+            J++;
+        }
+        rParSigma = sqrt(rParSigma);
+        psTrace(FACILITY, 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
@@ -164 +194 @@
         }
 
-        /* if (Chisq < min->value) {  */
+        // change in chisq/nDOF since last minimum
+        min->lastDelta = (min->value - Chisq) / pcm->nDOF;
+
+        // 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
+        // rho is positive if the new chisq is smaller
         if (rho >= -1e-6) {
-            min->lastDelta = (min->value - Chisq) / (source->pixels->numCols*source->pixels->numRows - params->n);
             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;
 
             // save the new convolved model image
             psFree (source->modelFlux);
             source->modelFlux = pmPCMdataSaveImage(pcm);
-        } else {
-            lambda *= 10.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;
+        }
         min->iter++;
 
+        // ending conditions:
+        // 1) hard limit : too many iterations
         done = (min->iter >= min->maxIter);
         
-        // check for convergence:
+        // 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 / pcm->nDOF;
-        if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) {
+        if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) {
+            continue;
+        }
+
+        // delta-chisq or rParSigma ?
+        if (min->chisqConvergence) {
             done |= (min->lastDelta < min->minTol);
+        } else {
+            done |= (rParSigma < min->minTol*pcm->nPar);
         }
     }
@@ -220 +300 @@
 
     // if the last improvement was at least as good as maxTol, accept the fit:
-    if (min->lastDelta <= min->maxTol) {
-        psTrace(FACILITY, 6, "---- end (true) ----\n");
-        return(true);
+    if (min->chisqConvergence) {
+        if (min->lastDelta <= min->maxTol) {
+            psTrace(FACILITY, 6, "---- end (true) ----\n");
+            return(true);
+        }
+    } else {
+        if (min->rParSigma <= min->maxTol*pcm->nPar) {
+            psTrace(FACILITY, 6, "---- end (true) ----\n");
+            return(true);
+        }
     }
     psTrace(FACILITY, 6, "---- end (false) ----\n");