Changeset 2197 for trunk/psLib/src/math/psMinimize.c

Timestamp:

Oct 26, 2004, 11:25:59 AM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/math/psMinimize.c (modified) (2 diffs)

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author GLF, MHPCC
  *
- *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-14 01:44:48 $
+ *  @version $Revision: 1.59 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -636 +636 @@
     psFree(currValueVec);
     psFree(newValueVec);
+
+    if ((min->iter < min->maxIter) ||
+            (min->lastDelta <= min->tol)) {
+        return(true);
+    }
+
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            "---- psMinimizeLMChi2() end (false) ----\n");
+    return(false);
+}
+
+/******************************************************************************
+psMyMinimizeLMChi2():  This routine will take an procedure which calculates
+an arbitrary function and it's derivative and minimize the chi-squared match
+between that function at the specified coords and the specified value at
+those coords.
+ 
+XXX: Do this:
+ After checking that all entries in the paramMask are 1 or 0, when
+ forming the A matrix from alpha, try this:
+ 
+     A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
+ 
+XXX: This is very different from what is specified in the SDR.  Must
+coordinate with IfA on new SDR.
+ 
+XXX: Do vector/image recycles.
+ 
+XXX: probably yErr will be part of the SDR.
+ 
+XXX: This must work for both F32 and F64.  F32 is currently implemented.
+     Note: since the LUD routines are only implemented in F64, then we
+     will have to convert all F32 input vectors to F64 regardless.  So,
+     the F64 port might be.
+ *****************************************************************************/
+bool psMyMinimizeLMChi2(psMinimization *min,
+                        psImage *covar,
+                        psVector *params,
+                        const psVector *paramMask,
+                        const psArray *x,
+                        const psVector *y,
+                        const psVector *yErr,
+                        psMyMinimizeLMChi2Func func)
+{
+    PS_CHECK_NULL_PTR_RETURN_NULL(min);
+    PS_CHECK_NULL_VECTOR_RETURN_NULL(params);
+    PS_CHECK_EMPTY_VECTOR_RETURN_NULL(params);
+    PS_CHECK_NULL_PTR_RETURN_NULL(x);
+    PS_CHECK_NULL_VECTOR_RETURN_NULL(y);
+    PS_CHECK_EMPTY_VECTOR_RETURN_NULL(y);
+    PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y);
+    if (paramMask != NULL) {
+        PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(params, paramMask);
+    }
+    if (yErr != NULL) {
+        PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(y, yErr);
+    }
+
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            "---- psMinimizeLMChi2() begin ----\n");
+    int numData = y->n;
+    int numParams = params->n;
+    int i;
+    int j;
+    int k;
+    int l;
+    int n;
+    int p;
+    psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64);
+
+    psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
+    psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);
+    psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);
+
+    psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);
+    psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+    psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
+
+    //    psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));
+    //    for (i=0;i<numData;i++) {
+    //        deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);
+    //    }
+    psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32);
+
+    psVector *currValueVec;
+    psVector *newValueVec;
+
+    float currChi2 = 0.0;
+    float newChi2 = 0.0;
+    float lamda = 0.00005;
+
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            "min->maxIter is %d\n", min->maxIter);
+    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            "min->tol is %f\n", min->tol);
+
+    for (p=0;p<numParams;p++) {
+        origParams->data.F32[p] = params->data.F32[p];
+    }
+
+    min->lastDelta = HUGE;
+    min->iter = 0;
+
+    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "------------------------------------------------------\n");
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
+
+        //
+        // Calculate the current values and chi-squared of the function.
+        //
+        currChi2 = 0.0;
+        currValueVec = func(deriv, params, x);
+        for (n=0;n<numData;n++) {
+            currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
+        }
+
+        for (p=0;p<numParams;p++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                "Current chi-squared is (%f)\n", currChi2);
+
+        //
+        // Mask elements of the derivative for each data point.
+        //
+        for (p=0;p<numParams;p++) {
+            if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
+                for (n=0;n<numData;n++) {
+                    deriv->data.F32[n][p] = 0.0;
+                }
+            }
+        }
+
+        //
+        // Calculate the BETA vector.
+        //
+        for (p=0;p<numParams;p++) {
+            beta->data.F64[p] = 0.0;
+            for (n=0;n<numData;n++) {
+                (beta->data.F64[p])+=
+                    (y->data.F32[n] - currValueVec->data.F32[n]) *
+                    deriv->data.F32[n][p];
+            }
+            // XXX: multiple by -1 here?
+            (beta->data.F64[p])*= -1.0;
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
+        }
+        psFree(currValueVec);
+
+        //
+        // Calculate the ALPHA matrix.
+        //
+        for (k=0;k<numParams;k++) {
+            for (l=0;l<numParams;l++) {
+                alpha->data.F32[k][l] = 0.0;
+                for (n=0;n<numData;n++) {
+                    alpha->data.F32[k][l]+= deriv->data.F32[n][k] *
+                                            deriv->data.F32[n][l];
+                }
+            }
+        }
+
+        //
+        // Calculate the matrix A.
+        //
+        for (j=0;j<numParams;j++) {
+            for (k=0;k<numParams;k++) {
+                if (j == k) {
+                    A->data.F64[j][k] =
+                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
+                } else {
+                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
+                }
+            }
+        }
+        for (j=0;j<numParams;j++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+            for (k=0;k<numParams;k++) {
+                psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "%f ", A->data.F64[j][k]);
+            }
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+        }
+
+        //
+        // Solve A * alpha = Beta
+        //
+        aOut = psMatrixLUD(aOut, perm, A);
+        paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
+
+        //
+        // Mask any masked parameters.
+        //
+        for (i=0;i<numParams;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
+            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
+                newParams->data.F32[i] = origParams->data.F32[i];
+            } else {
+                newParams->data.F32[i] = params->data.F32[i] -
+                                         (float) paramDeltasF64->data.F64[i];
+            }
+        }
+
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                "Calling func() with new parameters:\n");
+        for (i=0;i<numParams;i++) {
+            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+                    "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
+        }
+
+
+        //
+        // Calculate new function values.
+        //
+        newChi2 = 0.0;
+        newValueVec = func(deriv, newParams, x);
+        for (n=0;n<numData;n++) {
+            newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
+        }
+        psFree(newValueVec);
+
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
+
+        //
+        // If the new chi-squared is lower, then keep it.
+        //
+        if (currChi2 > newChi2) {
+            min->lastDelta = currChi2 - newChi2;
+            min->value = newChi2;
+
+            // We already masked params.
+            for (i=0;i<numParams;i++) {
+                params->data.F32[i] = (float) newParams->data.F32[i];
+            }
+            lamda*= 0.1;
+        } else {
+            lamda*= 10.0;
+        }
+        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+                "lamda is %f\n", lamda);
+        min->iter++;
+    }
+    psFree(beta);
+    psFree(perm);
+    psFree(paramDeltasF64);
+    psFree(origParams);
+    psFree(newParams);
+    psFree(alpha);
+    psFree(A);
+    psFree(aOut);
+    psFree(deriv);
 
     if ((min->iter < min->maxIter) ||