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

Timestamp:

May 31, 2013, 6:19:08 AM (13 years ago)

Author:

eugene

Message:

add a new version of psMinimizeLMM using a better convergence criterion -- delta parameters not delta chisq

File:

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

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

@@ -360 +360 @@
 }
 
-
 /******************************************************************************
 psMinimizeLMChi2():  This routine will take an procedure which calculates an
@@ -580 +579 @@
 }
 
+/******************************************************************************
+psMinimizeLMChi2(): This routine takes a function-pointer (func) which calculates an arbitrary
+function and it's derivatives and minimizes the chi-squared match between that function at the
+specified points and the specified value at those points.
+
+The original version of this function used a convergence criterion based on the change in
+chisq.  this has problems since it depends on the choice of points used to measure the fit.
+(consider a gaussian on a background : it 100 pixels are used -- and some or most contribute to
+the chisq -- and the delta-chisq is 10%, then the same change in model fit will yield a
+delta-chisq of 1% if 1000 pixels are used (all but the 100 measuring the background)).
+
+This implementation uses changes to the parameters and stops if they are no longer significant.
+
+This requires F32 input data; all internal calls use F32.
+  *****************************************************************************/
+bool psMinimizeLMChi2_Alt(
+    psMinimization *min,
+    psImage *covar,
+    psVector *params,
+    psMinConstraint *constraint,
+    const psArray *x,
+    const psVector *y,
+    const psVector *yWt,
+    psMinimizeLMChi2Func func)
+{
+    psTrace("psLib.math", 3, "---- begin ----\n");
+    PS_ASSERT_PTR_NON_NULL(min, false);
+    PS_ASSERT_VECTOR_NON_NULL(params, false);
+    PS_ASSERT_VECTOR_NON_EMPTY(params, false);
+    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
+    psVector *paramMask = NULL;
+    if (constraint != NULL) {
+        paramMask = constraint->paramMask;
+        if (paramMask != NULL) {
+            PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
+            PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false);
+        }
+    }
+    PS_ASSERT_PTR_NON_NULL(x, false);
+    for (psS32 i = 0 ; i < x->n ; i++) {
+        psVector *coord = (psVector *) (x->data[i]);
+        PS_ASSERT_VECTOR_NON_NULL(coord, false);
+        PS_ASSERT_VECTOR_TYPE(coord, PS_TYPE_F32, false);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_NON_EMPTY(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, false);
+    if (yWt != NULL) {
+        PS_ASSERT_VECTOR_TYPE(yWt, PS_TYPE_F32, false);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, yWt, false);
+    }
+    PS_ASSERT_PTR_NON_NULL(func, false);
+
+    psMinimizeLMLimitFunc checkLimits = NULL;
+    if (constraint) {
+        checkLimits = constraint->checkLimits;
+    }
+
+    // this function has test and current values for several things
+    // the current best value is in lower case
+    // the next guess value is in upper case
+
+    // allocate internal arrays (current vs Guess)
+    psImage *Alpha = NULL;
+    psVector *Beta = NULL;
+
+    // Alpha & Beta only contain elements to represent the unmasked parameters
+    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
+        psAbort ("programming error: no unmasked parameters to be fit\n");
+    }
+
+    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
+    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
+
+    psVector *dy     = NULL;
+    psF32 Chisq = 0.0;
+    psF32 lambda = 0.001;
+    psF32 dLinear = 0.0;
+
+    // the user provides the error or NULL.  we need to convert
+    // to appropriate weights
+    if (yWt != NULL) {
+        dy = (psVector *) yWt;
+    } else {
+        dy = psVectorAlloc(y->n, PS_TYPE_F32);
+        psVectorInit(dy, 1.0);
+    }
+
+    // number of degrees of freedom for this fit
+    int nDOF = dy->n - params->n;
+
+    // calculate initial alpha and beta, set chisq (min->value)
+    min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
+    if (isnan(min->value)) {
+        min->iter = min->maxIter;
+        psFree(alpha);
+        psFree(Alpha);
+        psFree(beta);
+        psFree(Beta);
+        psFree(Params);
+        return(false);
+    }
+    // dump some useful info if trace is defined
+    if (psTraceGetLevel("psLib.math") >= 6) {
+        p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
+        p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
+    }
+    if (psTraceGetLevel("psLib.math") >= 5) {
+        p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
+    }
+
+    // iterate until: (a) nIter = min->iter or (b) (chisq / ndof) < maxChisq and deltaChisq < minTol (but don't stop unless Chisq is finite)
+    bool done = (min->iter >= min->maxIter);
+    while (!done) {
+        psTrace("psLib.math", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+        psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+
+        // set a new guess for Alpha, Beta, Params
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
+            min->iter ++;
+            if (min->iter >=  min->maxIter) break;
+            lambda *= 10.0;
+            continue;
+        }
+
+        // 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_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)");
+        }
+        if (psTraceGetLevel("psLib.math") >= 5) {
+            p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
+        }
+
+        // calculate Chisq for new guess, update Alpha & Beta
+        Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func);
+        if (isnan(Chisq)) {
+            min->iter ++;
+            if (min->iter >= min->maxIter) break;
+            lambda *= 10.0;
+            continue;
+        }
+
+        // convergence criterion:
+        // compare the delta (min->value - Chisq) with the
+        // expected delta from the linear model (dLinear)
+        // accept new guess if it is an improvement (rho > 0), or else increase lambda
+        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);
+
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel("psLib.math") >= 6) {
+            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)");
+            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)");
+        }
+
+        // change in chisq/nDOF since last minimum
+        min->lastDelta = (min->value - Chisq) / nDOF;
+
+        /* rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho) */
+        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);
+            lambda *= 0.25;
+        } else {
+            lambda *= 10.0;
+        }
+        min->iter++;
+
+        // ending conditions:
+        // 1) hard limit : too many iterations
+        done = (min->iter >= min->maxIter);
+        
+        // 2) require chisqDOF < maxChisqDOF (if maxChisqDOF is not NAN)
+        float chisqDOF = Chisq / nDOF;
+        if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) {
+            continue;
+        }
+
+        // 3) require deltaChi > 1e-6
+        if (min->lastDelta < 1e-5) {
+            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);
+    }
+    psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+
+    // construct & return the covariance matrix (if requested)
+    if (covar != NULL) {
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
+            psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
+        }
+        // set covar values which are not masked
+        psImageInit (covar, 0.0);
+        for (int j = 0, J = 0; j < params->n; j++) {
+            if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
+                covar->data.F32[j][j] = 1.0;
+                continue;
+            }
+            for (int k = 0, K = 0; k < params->n; k++) {
+                if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue;
+                covar->data.F32[j][k] = Alpha->data.F32[J][K];
+                K++;
+            }
+            J++;
+        }
+    }
+
+    // free the internal temporary data
+    psFree(alpha);
+    psFree(Alpha);
+    psFree(beta);
+    psFree(Beta);
+    psFree(Params);
+    if (yWt == NULL) {
+        psFree(dy);
+    }
+
+    // if the last improvement was at least as good as maxTol, accept the fit:
+    if (min->lastDelta <= min->maxTol) {
+        psTrace("psLib.math", 6, "---- end (true) ----\n");
+        return(true);
+    }
+    psTrace("psLib.math", 6, "---- end (false) ----\n");
+    return(false);
+}
+
 bool psMinLM_AllocAB (psImage **Alpha, psVector **Beta, const psVector *params, const psVector *paramMask) {