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

Timestamp:

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

Author:

eugene

Message:

add new psMinimizeLMChi2 function with new convergence criterion; add threaded version of psImageSmoothNoMask

File:

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

trunk/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;
@@ -360 +364 @@
 }
 
+/******************************************************************************
+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;
+}
 
 /******************************************************************************
@@ -370 +398 @@
 XXX Make an F64 version?
   *****************************************************************************/
-bool psMinimizeLMChi2(
+bool psMinimizeLMChi2_Old(
     psMinimization *min,
     psImage *covar,
@@ -507 +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
@@ -580 +608 @@
 }
 
+/******************************************************************************
+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");
+    }
+
+    int nFitParams = Beta->n;
+    psImage *alpha   = psImageAlloc(nFitParams, nFitParams, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(nFitParams, 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;
+    psF32 nu = 2.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 - nFitParams;
+
+    // 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)");
+    }
+
+    bool done = (min->iter >= min->maxIter);
+    while (!done) {
+        psTrace("psLib.math", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+
+        if (min->chisqConvergence) {
+          psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        } else {
+          psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*nFitParams, min->maxTol*nFitParams);
+        }
+
+        // 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;
+            // ALT? lambda *= 2.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 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("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
+        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;
+            // ALT lambda *= 2.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: %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
+        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)
+
+        // 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->value = Chisq;
+            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
+            beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
+            params = psVectorCopy(params, Params, PS_TYPE_F32);
+        } 
+        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);
+        
+        // 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 / nDOF;
+        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*nFitParams);
+        }
+    }
+    psTrace("psLib.math", 5, "chisq: %f, last delta: %f, rParSigma: %f, Niter: %d\n", min->value, min->lastDelta, min->rParSigma, 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->chisqConvergence) {
+      if (min->lastDelta <= min->maxTol) {
+        psTrace("psLib.math", 6, "---- end (true) ----\n");
+        return(true);
+      }
+    } else {
+      if (min->rParSigma <= min->maxTol*nFitParams) {
+        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) {
 
@@ -625 +966 @@
     min->iter = 0;
     min->lastDelta = NAN;
+    min->rParSigma = NAN;
     min->maxChisqDOF = NAN;
 
+    // we default to the old algorithm for convergence
+    min->chisqConvergence = true;
+    min->gainFactorMode = 0;
+    
     return(min);
 }