Changeset 2228

Timestamp:

Oct 28, 2004, 12:47:00 PM (22 years ago)

Author:

gusciora

Message:

Modifying the minimization algorithms.

Location:

trunk/psLib/src

Files:

• dataManip/psConstants.h (modified) (2 diffs)
• dataManip/psMinimize.c (modified) (29 diffs)
• dataManip/psMinimize.h (modified) (3 diffs)
• dataManip/psStats.c (modified) (12 diffs)
• math/psConstants.h (modified) (2 diffs)
• math/psMinimize.c (modified) (29 diffs)
• math/psMinimize.h (modified) (3 diffs)
• math/psStats.c (modified) (12 diffs)
• sys/psTrace.h (modified) (2 diffs)
• sysUtils/psTrace.h (modified) (2 diffs)

trunk/psLib/src/dataManip/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 00:22:53 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -208 +208 @@
 p_psMemSetPersistent(NAME->data.V, true); \
 
+#define PS_VECTOR_DECLARE_ALLOC_STATIC(NAME, SIZE, TYPE) \
+static psVector *NAME = NULL; \
+if (NAME == NULL) { \
+    NAME = psVectorAlloc(SIZE, TYPE); \
+    p_psMemSetPersistent(NAME, true); \
+} \
 
 

trunk/psLib/src/dataManip/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.72 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 23:31:43 $
+ *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -408 +408 @@
 }
 
+
+/******************************************************************************
+XXX: We assume unnormalized gaussians.
+XXX: Currently, yErr is ignored.
+ *****************************************************************************/
+psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
+                                  const psVector *params,
+                                  const psArray *coords)
+{
+    psTrace(".psLib.dataManip.psMinimize", 4,
+            "---- psMinimizeLMChi2Gauss1D() begin ----\n");
+
+    PS_PTR_CHECK_NULL(coords, NULL);
+    PS_PTR_CHECK_NULL(params, NULL);
+    float x;
+    int i;
+    float mean = params->data.F32[0];
+    float stdev = params->data.F32[1];
+    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
+
+    psTrace(".psLib.dataManip.psMinimize", 6,
+            "(mean, stdev) is (%f, %f)\n", mean, stdev);
+
+    if (deriv == NULL) {
+        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
+    }
+
+    for (i=0;i<coords->n;i++) {
+        x = ((psVector *) (coords->data[i]))->data.F32[0];
+        out->data.F32[i] = psGaussian(x, mean, stdev, false);
+        //        psTrace(".psLib.dataManip.psMinimize", 6,
+        //            "out->data.F32[%d]: (x, m, s) is (%f, %f, %f) is %f\n", i, x, mean, stdev, out->data.F32[i]);
+    }
+
+    for (i=0;i<coords->n;i++) {
+        x = ((psVector *) (coords->data[i]))->data.F32[0];
+        float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
+        deriv->data.F32[i][0] = tmp / (stdev * stdev);
+
+        tmp = (x - mean) * (x - mean) *
+              psGaussian(x, mean, stdev, 0);
+        deriv->data.F32[i][1] = tmp / (stdev * stdev * stdev);
+    }
+
+    psTrace(".psLib.dataManip.psMinimize", 4,
+            "---- psMinimizeLMChi2Gauss1D() end ----\n");
+    return(out);
+}
+
+
 /******************************************************************************
 psMinimizeLMChi2():  This routine will take an procedure which calculates
@@ -429 +479 @@
 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,
+>     will have to convert all F32 input vectors to F64 regardless.  So,
      the F64 port might be.
  *****************************************************************************/
@@ -457 +507 @@
     }
 
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+    psTrace(".psLib.dataManip.psMinimize", 4,
             "---- psMinimizeLMChi2() begin ----\n");
     psS32 numData = y->n;
@@ -490 +540 @@
     float newChi2 = 0.0;
     float lamda = 0.00005;
-
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+    lamda = 0.05;
+
+    psTrace(".psLib.dataManip.psMinimize", 6,
             "min->maxIter is %d\n", min->maxIter);
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+    psTrace(".psLib.dataManip.psMinimize", 6,
             "min->tol is %f\n", min->tol);
 
@@ -504 +555 @@
 
     while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "------------------------------------------------------\n");
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
 
@@ -515 +566 @@
         currValueVec = func(deriv, params, x);
         for (n=0;n<numData;n++) {
-            currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
+            currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) *
+                       (currValueVec->data.F32[n] - y->data.F32[n]);
+            psTrace(".psLib.dataManip.psMinimize", 6,
+                    "data[%d], chi2 calculation+= (%f * %f)\n", n,
+                    currValueVec->data.F32[n], y->data.F32[n]);
         }
 
         for (p=0;p<numParams;p++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
         }
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+        psTrace(".psLib.dataManip.psMinimize", 6,
                 "Current chi-squared is (%f)\n", currChi2);
 
@@ -548 +603 @@
             // XXX: multiple by -1 here?
             (beta->data.F64[p])*= -1.0;
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
         }
@@ -580 +635 @@
         }
         for (j=0;j<numParams;j++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+            psTrace(".psLib.dataManip.psMinimize", 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");
+                psTrace(".psLib.dataManip.psMinimize", 6, "%f ", A->data.F64[j][k]);
+            }
+            psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n");
         }
 
@@ -597 +652 @@
         //
         for (i=0;i<numParams;i++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
@@ -607 +662 @@
         }
 
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+        psTrace(".psLib.dataManip.psMinimize", 6,
                 "Calling func() with new parameters:\n");
         for (i=0;i<numParams;i++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
         }
@@ -621 +676 @@
         newValueVec = func(deriv, newParams, x);
         for (n=0;n<numData;n++) {
-            newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
+            newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) *
+                      (newValueVec->data.F32[n] - y->data.F32[n]);
+
         }
         psFree(newValueVec);
 
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
 
@@ -640 +697 @@
             }
             lamda*= 0.1;
+            psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n");
         } else {
             lamda*= 10.0;
-        }
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n");
+        }
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "lamda is %f\n", lamda);
         min->iter++;
@@ -662 +721 @@
     }
 
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+    psTrace(".psLib.dataManip.psMinimize", 4,
             "---- psMinimizeLMChi2() end (false) ----\n");
     return(false);
@@ -942 +1001 @@
 }
 
+#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
+for (int i=0;i<BASE->n;i++) { \
+    if (BASEMASK->data.U8[i] == 0) { \
+        OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \
+    } else { \
+        OUT->data.F32[i] = BASE->data.F32[i]; \
+    } \
+} \
+
+
 /******************************************************************************
 p_psDetermineBracket():  This routine takes as input an arbitrary function,
@@ -947 +1016 @@
 function produces as output a bracket [a, b, c] such that
  
-        f(param + b * line) is less than f(param + a * line) and
-                                         f(param + c * line).
- 
-Algorithm: XXX completely ad hoc: start with the user-supplied starting
+f(param + b * line) is less than f(param + a * line) and
+f(param + c * line).
+ 
+Algorithm: 
+XXX completely ad hoc: 
+start with the user-supplied starting
 parameter and call that b.  Calculate a/c as a fractional amount
 smaller/larger than b.  Repeat this process until a local minimum is found.
  
-XXX: new algorithm: start at x=0, expand in one direction until the function
+XXX: 
+new algorithm: 
+start at x=0, expand in one direction until the function
 decreases.  Then you have two points in the bracket.  Keep going until it
 increases, or x is too large.  If thst does not work, expand in the other
 direction.
  
-XXX: This is F32 only
- *****************************************************************************/
+XXX: 
+This is F32 only.
+ 
+XXX: 
+output bracket vector should be an input as well.
+*****************************************************************************/
 psVector *p_psDetermineBracket(psVector *params,
                                psVector *line,
@@ -966 +1043 @@
                                const psArray *coords,
                                psMinimizePowellFunc func)
+{
+    float a = 0.0;
+    float b = 0.0;
+    float c = 0.0;
+    float fa = 0.0;
+    float fb = 0.0;
+    float fc = 0.0;
+    psS32 iter = 100;
+    float aDir = 0.0;
+    float cDir = 0.0;
+    float new_aDir = 0.0;
+    float new_cDir = 0.0;
+    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
+    float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
+    psVector *tmp = NULL;
+    psS32 i = 0;
+    psS32 null = 0;
+
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+            "---- p_psDetermineBracket() begin ----\n");
+
+    for (i=0;i<params->n;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            if (line->data.F32[i] >= FLT_EPSILON) {
+                null = 1;
+            }
+        }
+    }
+
+    if (null == 0) {
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 2,
+                "p_psDetermineBracket() called with zero line vector.\n");
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                "---- p_psDetermineBracket() end (NULL) ----\n");
+        psFree(bracket);
+        return(NULL);
+    }
+
+    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
+
+    b = 0;
+    a = -stepSize;
+    c = stepSize;
+
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
+    fa = func(tmp, coords);
+
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
+    fb = func(tmp, coords);
+
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
+    fc = func(tmp, coords);
+
+    if (fa < fb) {
+        aDir = -1;
+    } else {
+        aDir = 1;
+    }
+
+    if (fc < fb) {
+        cDir = -1;
+    } else {
+        cDir = 1;
+    }
+
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+            "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc);
+
+    while (iter > 0) {
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                "psDetermineBracket(): iteration %d\n", iter);
+        if ((fb < fa) && (fb < fc)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+        stepSize*= (1.0 + stepSize);
+        a =- stepSize;
+        c =+ stepSize;
+
+        for (i=0;i<params->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
+            } else {
+                tmp->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fa = func(tmp, coords);
+
+        for (i=0;i<params->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
+            } else {
+                tmp->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fc = func(tmp, coords);
+
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
+
+        if (fa < fb) {
+            new_aDir = -1;
+        } else {
+            new_aDir = 1;
+        }
+
+        if (fc < fb) {
+            new_cDir = -1;
+        } else {
+            new_cDir = 1;
+        }
+        if ((new_aDir == 1) && (aDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+
+        if ((new_cDir == 1) && (cDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+        aDir = new_aDir;
+        cDir = new_cDir;
+        iter--;
+    }
+    psFree(tmp);
+    psFree(bracket);
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+            "---- p_psDetermineBracket() end (NULL) ----\n");
+    return(NULL);
+}
+
+psVector *p_psDetermineBracket2(psVector *params,
+                                psVector *line,
+                                const psVector *paramMask,
+                                const psArray *coords,
+                                psMinimizePowellFunc func)
 {
     float a = 0.0;
@@ -1365 +1593 @@
     float fqp = 0.0;
     float diff = 0.0;
+    int iterationNumber = 0;
 
     psTrace(".psLib.dataManip.psMinimizePowell", 4,
@@ -1407 +1636 @@
     }
 
-    min->iter = 0;
-    while (min->iter < min->maxIter) {
-        min->iter++;
+    while (iterationNumber < min->maxIter) {
+        iterationNumber++;
         psTrace(".psLib.dataManip.psMinimizePowell", 6,
-                "psMinimizePowell() iteration %d\n", min->iter);
+                "psMinimizePowell() iteration %d\n", iterationNumber);
 
         // 3: For each dimension in params, move Q only in the vector v[i] to
@@ -1418 +1646 @@
         baseFuncVal = func(Q, coords);
         currFuncVal = baseFuncVal;
+        psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                "Current function value is %f\n", currFuncVal);
+
         diff = 0.0;
         biggestDiff = 0;
@@ -1423 +1654 @@
         for (i=0;i<numDims;i++) {
             if (myParamMask->data.U8[i] == 0) {
-                dummyMin.maxIter = 100;
+                dummyMin.maxIter = 20;
                 dummyMin.tol = 0.01;
                 mul = p_psLineMin(&dummyMin, Q, v[i], myParamMask, coords, func);
@@ -1433 +1664 @@
             }
         }
-
+        psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                "New function value is %f\n", currFuncVal);
         // 4: Set the vector u = Q - P
         for (i=0;i<numDims;i++) {
@@ -1450 +1682 @@
 
         // 5: Move Q only in the direction u, and minimize the function.
-        psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
         for (i=0;i<numDims;i++) {
             psTrace(".psLib.dataManip.psMinimizePowell", 6,
@@ -1470 +1701 @@
                 psFree(myParamMask);
             }
+            min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
                     "---- psMinimizePowell() end (true) ----\n");
@@ -1513 +1745 @@
             psFree(pQP);
             psFree(Q);
+            if (paramMask == NULL) {
+                psFree(myParamMask);
+            }
+            min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
                     "---- psMinimizePowell() end (true) ----\n");
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             return(true);
         }
@@ -1532 +1765 @@
     psFree(pQP);
     psFree(Q);
+    min->iter = iterationNumber;
     psTrace(".psLib.dataManip.psMinimizePowell", 4,
             "---- psMinimizePowell() end (false) ----\n");
@@ -1549 +1783 @@
                        const psArray *coords)
 {
+    psTrace(".psLib.dataManip.myPowellChi2Func", 4,
+            "---- myPowellChi2Func() begin ----\n");
+    PS_VECTOR_CHECK_NULL(params, NAN);
+    PS_VECTOR_CHECK_EMPTY(params, NAN);
+    PS_VECTOR_CHECK_NULL(myValue, NAN);
+    PS_VECTOR_CHECK_EMPTY(myValue, NAN);
+    PS_PTR_CHECK_NULL(coords, NAN);
+
     float chi2 = 0.0;
     float d;
@@ -1555 +1797 @@
 
     tmp = Chi2PowellFunc(params, coords);
-    for (i=0;i<coords->n;i++) {
-        d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i];
-        chi2+= d * d;
+    if (myError == NULL) {
+        for (i=0;i<coords->n;i++) {
+            d = (tmp->data.F32[i] - myValue->data.F32[i]);
+            chi2+= d * d;
+        }
+    } else {
+        for (i=0;i<coords->n;i++) {
+            d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i];
+            chi2+= d * d;
+        }
     }
     psFree(tmp);
+    psTrace(".psLib.dataManip.myPowellChi2Func", 4,
+            "---- myPowellChi2Func() end (%f) ----\n", chi2);
     return(chi2);
 }
@@ -1584 +1835 @@
     myValue = (psVector *) value;
     myError = (psVector *) error;
+
     Chi2PowellFunc = func;
 

trunk/psLib/src/dataManip/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 23:31:43 $
+ *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -81 +81 @@
                                   const psArray *x);
 
-/*
-psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
-                                  const psVector *params,
-                                  const psArray *coords)
-{
-    psVector *x = (psVector *) coords->data[0];
-    psVector *y = (psVector *) coords->data[1];
-    psVector *yErr = NULL;
-    int i;
-    int j;
-    float normalization = params->data.F32[0];
-    float center = params->data.F32[1];
-    float stdev = params->data.F32[2];
-    psVector *out = psVectorAlloc(x->n, PS_TYPE_F32);
- 
- 
-    if (deriv == NULL) {
-        deriv = psImageAlloc(params->n, x->n, PS_TYPE_F32);
-    }
- 
-    if (coords->n == 2) {
-        yErr = (psVector *) coords->data[2];
-    }
- 
-    for (i=0;i<x->n;i++) {
-        out->data.F32[i] = psGaussian(x->data.F32[i], center, stdev, 1);
-    }
- 
-    for (i=0;i<x->n;i++) {
-        for (j=0;j<params->n;j++) {
-            deriv->data.F32[i][j] = the derivative of the gaussian w.r.t. the j-th parameter at the coord x[i]
-        }
-    }
-    return(out);
-}
-*/
-
 psBool psMinimizeLMChi2(psMinimization *min,
                         psImage *covar,
@@ -137 +100 @@
                         psMinimizePowellFunc func);
 
+psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
+                                  const psVector *params,
+                                  const psArray *coords);
+
 typedef
 psVector *(*psMinimizeChi2PowellFunc)(const psVector *params,

trunk/psLib/src/dataManip/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 00:46:04 $
-n *
+ *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
+ *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -123 +123 @@
     default:
         return false;
-    }
-}
-
-/*****************************************************************************
-p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
-function that simply prints a vector to STDOUT.  Used primarily for
-debugging.
- *****************************************************************************/
-
-void p_psVectorPrint(psVector* myVector, psVector* maskVector, psU32 maskVal, psStats* stats)
-{
-    psS32 i = 0;                  // Loop index variable.
-
-    for (i = 0; i < myVector->n; i++) {
-        if (maskVector != NULL)
-            printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
-        else
-            printf("Element %d is %f\n", i, myVector->data.F32[i]);
     }
 }
@@ -428 +410 @@
     // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    unsortedVector->n = unsortedVector->nalloc;
-
     sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    sortedVector->n = sortedVector->nalloc;
 
     // Determine if we must only use data points within a min/max range.
@@ -579 +558 @@
     // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    unsortedVector->n = unsortedVector->nalloc;
     sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    sortedVector->n = sortedVector->nalloc;
 
     // Determine if we must only use data points within a min/max range.
@@ -817 +794 @@
 p_psNormalizeVectorF32(myData): this is a private function which normalizes the
 elements of a vector to a range between 0.0 and 1.0.
+ 
+XXX: how about an arbitrary p_psNormalizeVectorF32(myData, min, max)?
  *****************************************************************************/
 void p_psNormalizeVectorF32(psVector* myData)
@@ -847 +826 @@
 elements of a vector to a range between -1.0 and 1.0.
 XXX: 0-1 or -1:1?
+ 
+XXX: how about an arbitrary p_psNormalizeVectorF64(myData, min, max)?
  *****************************************************************************/
 void p_psNormalizeVectorF64(psVector* myData)
@@ -876 +857 @@
 }
 
+// XXX: how about an arbitrary p_psNormalizeVector(myData, min, max)?
 void p_psNormalizeVector(psVector* myData)
 {
@@ -885 +867 @@
         psError(__func__, "Unalowable data type.\n");
     }
-}
-
-/*****************************************************************************
-p_psGaussian(myData, myParams): an internal function, used by robust stats,
-intended to compute the Gaussian with a specified sigma/mean, at the
-specified data point.
- *****************************************************************************/
-float p_psGaussian(const psVector* restrict myData,
-                   const psVector* restrict myParams)
-{
-    float x = myData->data.F32[0];
-    float mean = myParams->data.F32[0];
-    float stdev = myParams->data.F32[1];
-    float tmp = exp(-((x - mean) * (x - mean)) / (2.0 * stdev * stdev));
-
-    tmp /= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
-
-    // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
-    return (tmp);
-}
-
-/*****************************************************************************
-p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which
-calculates the specified partial derivative of the above Gaussian function.
- *****************************************************************************/
-float p_psGaussianDeriv(const psVector* restrict myData,
-                        const psVector* restrict myParams, psS32 whichParam)
-{
-    float x = myData->data.F32[0];
-    float mean = myParams->data.F32[0];
-    float stdev = myParams->data.F32[1];
-    float tmp = 0.0;
-
-    if (whichParam == 0) {
-        // Return the derivative w.r.t. the mean.
-        tmp = (x - mean) * p_psGaussian(myData, myParams);
-        tmp /= (stdev * stdev);
-    } else if (whichParam == 1) {
-        // Return the derivative w.r.t. the stdev.
-        tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams);
-        tmp /= (stdev * stdev * stdev);
-    }
-    printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp);
-
-    return (tmp);
-}
-
-/*****************************************************************************
-p_psQuadratic(myData, myParams): an internal function, used by robust stats,
-intended to compute a quadratic, with the specified parameters, at the
-specified data point.
- *****************************************************************************/
-float p_psQuadratic(const psVector* restrict myParams, const psVector* restrict myCoords)
-{
-    float x = myCoords->data.F32[0];
-    float A = myParams->data.F32[0];
-    float B = myParams->data.F32[1];
-    float C = myParams->data.F32[2];
-    float tmp = 0.0;
-
-    tmp = (A * x * x) + (B * x) + C;
-    return (tmp);
-}
-
-/*****************************************************************************
-p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which
-calculates the specified partial derivative of the above quadratic function.
- *****************************************************************************/
-float p_psQuadraticDeriv(const psVector* restrict myParams,
-                         const psVector* restrict myCoords, psS32 whichParamDeriv)
-{
-    float x = myCoords->data.F32[0];
-    float tmp = 0.0;
-
-    if (whichParamDeriv == 0) {
-        tmp = x * x;
-    } else if (whichParamDeriv == 1) {
-        tmp = x;
-    } else if (whichParamDeriv == 2) {
-        tmp = 1.0;
-    }
-
-    return (tmp);
 }
 
@@ -1022 +921 @@
 and i+1 is used for x[i]).  It then determines for what value x does that
 quadratic f(x) = yVal (the input parameter).
- 
-// XXX: Use static variables.
 *****************************************************************************/
 float fitQuadraticSearchForYThenReturnX(psVector *xVec,
@@ -1030 +927 @@
                                         float yVal)
 {
-    //    static psVector* x = NULL;
-    //    static psVector* y = NULL;
-    //    static psVector* yErr = NULL;
-    psVector* x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* y = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* yErr = psVectorAlloc(3, PS_TYPE_F64);
-
-    /*
-        if (x == NULL) {
-            x = psVectorAlloc(3, PS_TYPE_F64);
-            p_psMemSetPersistent(x, true);
-        }
-        if (y == NULL) {
-            y = psVectorAlloc(3, PS_TYPE_F64);
-            p_psMemSetPersistent(y, true);
-        }
-        if (yErr == NULL) {
-            yErr = psVectorAlloc(3, PS_TYPE_F64);
-            p_psMemSetPersistent(yErr, true);
-        }
-    */
-    psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
+    PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);
+    PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);
+    PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);
+    static psPolynomial1D* myPoly = NULL;
+
+    if (myPoly == NULL) {
+        myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
+        p_psMemSetPersistent(myPoly, true);
+    }
     float tmpFloat;
 
@@ -1090 +974 @@
     }
 
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myPoly);
     return(tmpFloat);
 }
@@ -1475 +1355 @@
 
 /******************************************************************************
-psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure
+psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure
 with the specifed bounds.
  

trunk/psLib/src/math/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 00:22:53 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -208 +208 @@
 p_psMemSetPersistent(NAME->data.V, true); \
 
+#define PS_VECTOR_DECLARE_ALLOC_STATIC(NAME, SIZE, TYPE) \
+static psVector *NAME = NULL; \
+if (NAME == NULL) { \
+    NAME = psVectorAlloc(SIZE, TYPE); \
+    p_psMemSetPersistent(NAME, true); \
+} \
 
 

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.72 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 23:31:43 $
+ *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -408 +408 @@
 }
 
+
+/******************************************************************************
+XXX: We assume unnormalized gaussians.
+XXX: Currently, yErr is ignored.
+ *****************************************************************************/
+psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
+                                  const psVector *params,
+                                  const psArray *coords)
+{
+    psTrace(".psLib.dataManip.psMinimize", 4,
+            "---- psMinimizeLMChi2Gauss1D() begin ----\n");
+
+    PS_PTR_CHECK_NULL(coords, NULL);
+    PS_PTR_CHECK_NULL(params, NULL);
+    float x;
+    int i;
+    float mean = params->data.F32[0];
+    float stdev = params->data.F32[1];
+    psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32);
+
+    psTrace(".psLib.dataManip.psMinimize", 6,
+            "(mean, stdev) is (%f, %f)\n", mean, stdev);
+
+    if (deriv == NULL) {
+        deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32);
+    }
+
+    for (i=0;i<coords->n;i++) {
+        x = ((psVector *) (coords->data[i]))->data.F32[0];
+        out->data.F32[i] = psGaussian(x, mean, stdev, false);
+        //        psTrace(".psLib.dataManip.psMinimize", 6,
+        //            "out->data.F32[%d]: (x, m, s) is (%f, %f, %f) is %f\n", i, x, mean, stdev, out->data.F32[i]);
+    }
+
+    for (i=0;i<coords->n;i++) {
+        x = ((psVector *) (coords->data[i]))->data.F32[0];
+        float tmp = (x - mean) * psGaussian(x, mean, stdev, false);
+        deriv->data.F32[i][0] = tmp / (stdev * stdev);
+
+        tmp = (x - mean) * (x - mean) *
+              psGaussian(x, mean, stdev, 0);
+        deriv->data.F32[i][1] = tmp / (stdev * stdev * stdev);
+    }
+
+    psTrace(".psLib.dataManip.psMinimize", 4,
+            "---- psMinimizeLMChi2Gauss1D() end ----\n");
+    return(out);
+}
+
+
 /******************************************************************************
 psMinimizeLMChi2():  This routine will take an procedure which calculates
@@ -429 +479 @@
 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,
+>     will have to convert all F32 input vectors to F64 regardless.  So,
      the F64 port might be.
  *****************************************************************************/
@@ -457 +507 @@
     }
 
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+    psTrace(".psLib.dataManip.psMinimize", 4,
             "---- psMinimizeLMChi2() begin ----\n");
     psS32 numData = y->n;
@@ -490 +540 @@
     float newChi2 = 0.0;
     float lamda = 0.00005;
-
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+    lamda = 0.05;
+
+    psTrace(".psLib.dataManip.psMinimize", 6,
             "min->maxIter is %d\n", min->maxIter);
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+    psTrace(".psLib.dataManip.psMinimize", 6,
             "min->tol is %f\n", min->tol);
 
@@ -504 +555 @@
 
     while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "------------------------------------------------------\n");
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
 
@@ -515 +566 @@
         currValueVec = func(deriv, params, x);
         for (n=0;n<numData;n++) {
-            currChi2+= (currValueVec->data.F32[n] * currValueVec->data.F32[n]);
+            currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) *
+                       (currValueVec->data.F32[n] - y->data.F32[n]);
+            psTrace(".psLib.dataManip.psMinimize", 6,
+                    "data[%d], chi2 calculation+= (%f * %f)\n", n,
+                    currValueVec->data.F32[n], y->data.F32[n]);
         }
 
         for (p=0;p<numParams;p++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
         }
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+        psTrace(".psLib.dataManip.psMinimize", 6,
                 "Current chi-squared is (%f)\n", currChi2);
 
@@ -548 +603 @@
             // XXX: multiple by -1 here?
             (beta->data.F64[p])*= -1.0;
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
         }
@@ -580 +635 @@
         }
         for (j=0;j<numParams;j++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n");
+            psTrace(".psLib.dataManip.psMinimize", 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");
+                psTrace(".psLib.dataManip.psMinimize", 6, "%f ", A->data.F64[j][k]);
+            }
+            psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n");
         }
 
@@ -597 +652 @@
         //
         for (i=0;i<numParams;i++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
@@ -607 +662 @@
         }
 
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+        psTrace(".psLib.dataManip.psMinimize", 6,
                 "Calling func() with new parameters:\n");
         for (i=0;i<numParams;i++) {
-            psTrace(".psLib.dataManip.psMinimizeLMChi2", 6,
+            psTrace(".psLib.dataManip.psMinimize", 6,
                     "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
         }
@@ -621 +676 @@
         newValueVec = func(deriv, newParams, x);
         for (n=0;n<numData;n++) {
-            newChi2+= (newValueVec->data.F32[n] * newValueVec->data.F32[n]);
+            newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) *
+                      (newValueVec->data.F32[n] - y->data.F32[n]);
+
         }
         psFree(newValueVec);
 
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
 
@@ -640 +697 @@
             }
             lamda*= 0.1;
+            psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n");
         } else {
             lamda*= 10.0;
-        }
-        psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+            psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n");
+        }
+        psTrace(".psLib.dataManip.psMinimize", 4,
                 "lamda is %f\n", lamda);
         min->iter++;
@@ -662 +721 @@
     }
 
-    psTrace(".psLib.dataManip.psMinimizeLMChi2", 4,
+    psTrace(".psLib.dataManip.psMinimize", 4,
             "---- psMinimizeLMChi2() end (false) ----\n");
     return(false);
@@ -942 +1001 @@
 }
 
+#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
+for (int i=0;i<BASE->n;i++) { \
+    if (BASEMASK->data.U8[i] == 0) { \
+        OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \
+    } else { \
+        OUT->data.F32[i] = BASE->data.F32[i]; \
+    } \
+} \
+
+
 /******************************************************************************
 p_psDetermineBracket():  This routine takes as input an arbitrary function,
@@ -947 +1016 @@
 function produces as output a bracket [a, b, c] such that
  
-        f(param + b * line) is less than f(param + a * line) and
-                                         f(param + c * line).
- 
-Algorithm: XXX completely ad hoc: start with the user-supplied starting
+f(param + b * line) is less than f(param + a * line) and
+f(param + c * line).
+ 
+Algorithm: 
+XXX completely ad hoc: 
+start with the user-supplied starting
 parameter and call that b.  Calculate a/c as a fractional amount
 smaller/larger than b.  Repeat this process until a local minimum is found.
  
-XXX: new algorithm: start at x=0, expand in one direction until the function
+XXX: 
+new algorithm: 
+start at x=0, expand in one direction until the function
 decreases.  Then you have two points in the bracket.  Keep going until it
 increases, or x is too large.  If thst does not work, expand in the other
 direction.
  
-XXX: This is F32 only
- *****************************************************************************/
+XXX: 
+This is F32 only.
+ 
+XXX: 
+output bracket vector should be an input as well.
+*****************************************************************************/
 psVector *p_psDetermineBracket(psVector *params,
                                psVector *line,
@@ -966 +1043 @@
                                const psArray *coords,
                                psMinimizePowellFunc func)
+{
+    float a = 0.0;
+    float b = 0.0;
+    float c = 0.0;
+    float fa = 0.0;
+    float fb = 0.0;
+    float fc = 0.0;
+    psS32 iter = 100;
+    float aDir = 0.0;
+    float cDir = 0.0;
+    float new_aDir = 0.0;
+    float new_cDir = 0.0;
+    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
+    float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE;
+    psVector *tmp = NULL;
+    psS32 i = 0;
+    psS32 null = 0;
+
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+            "---- p_psDetermineBracket() begin ----\n");
+
+    for (i=0;i<params->n;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            if (line->data.F32[i] >= FLT_EPSILON) {
+                null = 1;
+            }
+        }
+    }
+
+    if (null == 0) {
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 2,
+                "p_psDetermineBracket() called with zero line vector.\n");
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                "---- p_psDetermineBracket() end (NULL) ----\n");
+        psFree(bracket);
+        return(NULL);
+    }
+
+    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
+
+    b = 0;
+    a = -stepSize;
+    c = stepSize;
+
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
+    fa = func(tmp, coords);
+
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
+    fb = func(tmp, coords);
+
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
+    fc = func(tmp, coords);
+
+    if (fa < fb) {
+        aDir = -1;
+    } else {
+        aDir = 1;
+    }
+
+    if (fc < fb) {
+        cDir = -1;
+    } else {
+        cDir = 1;
+    }
+
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+            "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc);
+
+    while (iter > 0) {
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                "psDetermineBracket(): iteration %d\n", iter);
+        if ((fb < fa) && (fb < fc)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+        stepSize*= (1.0 + stepSize);
+        a =- stepSize;
+        c =+ stepSize;
+
+        for (i=0;i<params->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
+            } else {
+                tmp->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fa = func(tmp, coords);
+
+        for (i=0;i<params->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
+            } else {
+                tmp->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fc = func(tmp, coords);
+
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
+
+        if (fa < fb) {
+            new_aDir = -1;
+        } else {
+            new_aDir = 1;
+        }
+
+        if (fc < fb) {
+            new_cDir = -1;
+        } else {
+            new_cDir = 1;
+        }
+        if ((new_aDir == 1) && (aDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+
+        if ((new_cDir == 1) && (cDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+        aDir = new_aDir;
+        cDir = new_cDir;
+        iter--;
+    }
+    psFree(tmp);
+    psFree(bracket);
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+            "---- p_psDetermineBracket() end (NULL) ----\n");
+    return(NULL);
+}
+
+psVector *p_psDetermineBracket2(psVector *params,
+                                psVector *line,
+                                const psVector *paramMask,
+                                const psArray *coords,
+                                psMinimizePowellFunc func)
 {
     float a = 0.0;
@@ -1365 +1593 @@
     float fqp = 0.0;
     float diff = 0.0;
+    int iterationNumber = 0;
 
     psTrace(".psLib.dataManip.psMinimizePowell", 4,
@@ -1407 +1636 @@
     }
 
-    min->iter = 0;
-    while (min->iter < min->maxIter) {
-        min->iter++;
+    while (iterationNumber < min->maxIter) {
+        iterationNumber++;
         psTrace(".psLib.dataManip.psMinimizePowell", 6,
-                "psMinimizePowell() iteration %d\n", min->iter);
+                "psMinimizePowell() iteration %d\n", iterationNumber);
 
         // 3: For each dimension in params, move Q only in the vector v[i] to
@@ -1418 +1646 @@
         baseFuncVal = func(Q, coords);
         currFuncVal = baseFuncVal;
+        psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                "Current function value is %f\n", currFuncVal);
+
         diff = 0.0;
         biggestDiff = 0;
@@ -1423 +1654 @@
         for (i=0;i<numDims;i++) {
             if (myParamMask->data.U8[i] == 0) {
-                dummyMin.maxIter = 100;
+                dummyMin.maxIter = 20;
                 dummyMin.tol = 0.01;
                 mul = p_psLineMin(&dummyMin, Q, v[i], myParamMask, coords, func);
@@ -1433 +1664 @@
             }
         }
-
+        psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                "New function value is %f\n", currFuncVal);
         // 4: Set the vector u = Q - P
         for (i=0;i<numDims;i++) {
@@ -1450 +1682 @@
 
         // 5: Move Q only in the direction u, and minimize the function.
-        psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
         for (i=0;i<numDims;i++) {
             psTrace(".psLib.dataManip.psMinimizePowell", 6,
@@ -1470 +1701 @@
                 psFree(myParamMask);
             }
+            min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
                     "---- psMinimizePowell() end (true) ----\n");
@@ -1513 +1745 @@
             psFree(pQP);
             psFree(Q);
+            if (paramMask == NULL) {
+                psFree(myParamMask);
+            }
+            min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
                     "---- psMinimizePowell() end (true) ----\n");
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             return(true);
         }
@@ -1532 +1765 @@
     psFree(pQP);
     psFree(Q);
+    min->iter = iterationNumber;
     psTrace(".psLib.dataManip.psMinimizePowell", 4,
             "---- psMinimizePowell() end (false) ----\n");
@@ -1549 +1783 @@
                        const psArray *coords)
 {
+    psTrace(".psLib.dataManip.myPowellChi2Func", 4,
+            "---- myPowellChi2Func() begin ----\n");
+    PS_VECTOR_CHECK_NULL(params, NAN);
+    PS_VECTOR_CHECK_EMPTY(params, NAN);
+    PS_VECTOR_CHECK_NULL(myValue, NAN);
+    PS_VECTOR_CHECK_EMPTY(myValue, NAN);
+    PS_PTR_CHECK_NULL(coords, NAN);
+
     float chi2 = 0.0;
     float d;
@@ -1555 +1797 @@
 
     tmp = Chi2PowellFunc(params, coords);
-    for (i=0;i<coords->n;i++) {
-        d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i];
-        chi2+= d * d;
+    if (myError == NULL) {
+        for (i=0;i<coords->n;i++) {
+            d = (tmp->data.F32[i] - myValue->data.F32[i]);
+            chi2+= d * d;
+        }
+    } else {
+        for (i=0;i<coords->n;i++) {
+            d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i];
+            chi2+= d * d;
+        }
     }
     psFree(tmp);
+    psTrace(".psLib.dataManip.myPowellChi2Func", 4,
+            "---- myPowellChi2Func() end (%f) ----\n", chi2);
     return(chi2);
 }
@@ -1584 +1835 @@
     myValue = (psVector *) value;
     myError = (psVector *) error;
+
     Chi2PowellFunc = func;
 

trunk/psLib/src/math/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 23:31:43 $
+ *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -81 +81 @@
                                   const psArray *x);
 
-/*
-psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
-                                  const psVector *params,
-                                  const psArray *coords)
-{
-    psVector *x = (psVector *) coords->data[0];
-    psVector *y = (psVector *) coords->data[1];
-    psVector *yErr = NULL;
-    int i;
-    int j;
-    float normalization = params->data.F32[0];
-    float center = params->data.F32[1];
-    float stdev = params->data.F32[2];
-    psVector *out = psVectorAlloc(x->n, PS_TYPE_F32);
- 
- 
-    if (deriv == NULL) {
-        deriv = psImageAlloc(params->n, x->n, PS_TYPE_F32);
-    }
- 
-    if (coords->n == 2) {
-        yErr = (psVector *) coords->data[2];
-    }
- 
-    for (i=0;i<x->n;i++) {
-        out->data.F32[i] = psGaussian(x->data.F32[i], center, stdev, 1);
-    }
- 
-    for (i=0;i<x->n;i++) {
-        for (j=0;j<params->n;j++) {
-            deriv->data.F32[i][j] = the derivative of the gaussian w.r.t. the j-th parameter at the coord x[i]
-        }
-    }
-    return(out);
-}
-*/
-
 psBool psMinimizeLMChi2(psMinimization *min,
                         psImage *covar,
@@ -137 +100 @@
                         psMinimizePowellFunc func);
 
+psVector *psMinimizeLMChi2Gauss1D(psImage *deriv,
+                                  const psVector *params,
+                                  const psArray *coords);
+
 typedef
 psVector *(*psMinimizeChi2PowellFunc)(const psVector *params,

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 00:46:04 $
-n *
+ *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:46:57 $
+ *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
  */
@@ -123 +123 @@
     default:
         return false;
-    }
-}
-
-/*****************************************************************************
-p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
-function that simply prints a vector to STDOUT.  Used primarily for
-debugging.
- *****************************************************************************/
-
-void p_psVectorPrint(psVector* myVector, psVector* maskVector, psU32 maskVal, psStats* stats)
-{
-    psS32 i = 0;                  // Loop index variable.
-
-    for (i = 0; i < myVector->n; i++) {
-        if (maskVector != NULL)
-            printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
-        else
-            printf("Element %d is %f\n", i, myVector->data.F32[i]);
     }
 }
@@ -428 +410 @@
     // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    unsortedVector->n = unsortedVector->nalloc;
-
     sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    sortedVector->n = sortedVector->nalloc;
 
     // Determine if we must only use data points within a min/max range.
@@ -579 +558 @@
     // Allocate temporary vectors for the data.
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    unsortedVector->n = unsortedVector->nalloc;
     sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
-    sortedVector->n = sortedVector->nalloc;
 
     // Determine if we must only use data points within a min/max range.
@@ -817 +794 @@
 p_psNormalizeVectorF32(myData): this is a private function which normalizes the
 elements of a vector to a range between 0.0 and 1.0.
+ 
+XXX: how about an arbitrary p_psNormalizeVectorF32(myData, min, max)?
  *****************************************************************************/
 void p_psNormalizeVectorF32(psVector* myData)
@@ -847 +826 @@
 elements of a vector to a range between -1.0 and 1.0.
 XXX: 0-1 or -1:1?
+ 
+XXX: how about an arbitrary p_psNormalizeVectorF64(myData, min, max)?
  *****************************************************************************/
 void p_psNormalizeVectorF64(psVector* myData)
@@ -876 +857 @@
 }
 
+// XXX: how about an arbitrary p_psNormalizeVector(myData, min, max)?
 void p_psNormalizeVector(psVector* myData)
 {
@@ -885 +867 @@
         psError(__func__, "Unalowable data type.\n");
     }
-}
-
-/*****************************************************************************
-p_psGaussian(myData, myParams): an internal function, used by robust stats,
-intended to compute the Gaussian with a specified sigma/mean, at the
-specified data point.
- *****************************************************************************/
-float p_psGaussian(const psVector* restrict myData,
-                   const psVector* restrict myParams)
-{
-    float x = myData->data.F32[0];
-    float mean = myParams->data.F32[0];
-    float stdev = myParams->data.F32[1];
-    float tmp = exp(-((x - mean) * (x - mean)) / (2.0 * stdev * stdev));
-
-    tmp /= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
-
-    // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
-    return (tmp);
-}
-
-/*****************************************************************************
-p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which
-calculates the specified partial derivative of the above Gaussian function.
- *****************************************************************************/
-float p_psGaussianDeriv(const psVector* restrict myData,
-                        const psVector* restrict myParams, psS32 whichParam)
-{
-    float x = myData->data.F32[0];
-    float mean = myParams->data.F32[0];
-    float stdev = myParams->data.F32[1];
-    float tmp = 0.0;
-
-    if (whichParam == 0) {
-        // Return the derivative w.r.t. the mean.
-        tmp = (x - mean) * p_psGaussian(myData, myParams);
-        tmp /= (stdev * stdev);
-    } else if (whichParam == 1) {
-        // Return the derivative w.r.t. the stdev.
-        tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams);
-        tmp /= (stdev * stdev * stdev);
-    }
-    printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp);
-
-    return (tmp);
-}
-
-/*****************************************************************************
-p_psQuadratic(myData, myParams): an internal function, used by robust stats,
-intended to compute a quadratic, with the specified parameters, at the
-specified data point.
- *****************************************************************************/
-float p_psQuadratic(const psVector* restrict myParams, const psVector* restrict myCoords)
-{
-    float x = myCoords->data.F32[0];
-    float A = myParams->data.F32[0];
-    float B = myParams->data.F32[1];
-    float C = myParams->data.F32[2];
-    float tmp = 0.0;
-
-    tmp = (A * x * x) + (B * x) + C;
-    return (tmp);
-}
-
-/*****************************************************************************
-p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which
-calculates the specified partial derivative of the above quadratic function.
- *****************************************************************************/
-float p_psQuadraticDeriv(const psVector* restrict myParams,
-                         const psVector* restrict myCoords, psS32 whichParamDeriv)
-{
-    float x = myCoords->data.F32[0];
-    float tmp = 0.0;
-
-    if (whichParamDeriv == 0) {
-        tmp = x * x;
-    } else if (whichParamDeriv == 1) {
-        tmp = x;
-    } else if (whichParamDeriv == 2) {
-        tmp = 1.0;
-    }
-
-    return (tmp);
 }
 
@@ -1022 +921 @@
 and i+1 is used for x[i]).  It then determines for what value x does that
 quadratic f(x) = yVal (the input parameter).
- 
-// XXX: Use static variables.
 *****************************************************************************/
 float fitQuadraticSearchForYThenReturnX(psVector *xVec,
@@ -1030 +927 @@
                                         float yVal)
 {
-    //    static psVector* x = NULL;
-    //    static psVector* y = NULL;
-    //    static psVector* yErr = NULL;
-    psVector* x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* y = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* yErr = psVectorAlloc(3, PS_TYPE_F64);
-
-    /*
-        if (x == NULL) {
-            x = psVectorAlloc(3, PS_TYPE_F64);
-            p_psMemSetPersistent(x, true);
-        }
-        if (y == NULL) {
-            y = psVectorAlloc(3, PS_TYPE_F64);
-            p_psMemSetPersistent(y, true);
-        }
-        if (yErr == NULL) {
-            yErr = psVectorAlloc(3, PS_TYPE_F64);
-            p_psMemSetPersistent(yErr, true);
-        }
-    */
-    psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
+    PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);
+    PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);
+    PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);
+    static psPolynomial1D* myPoly = NULL;
+
+    if (myPoly == NULL) {
+        myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
+        p_psMemSetPersistent(myPoly, true);
+    }
     float tmpFloat;
 
@@ -1090 +974 @@
     }
 
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
-    psFree(myPoly);
     return(tmpFloat);
 }
@@ -1475 +1355 @@
 
 /******************************************************************************
-psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure
+psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure
 with the specifed bounds.
  

trunk/psLib/src/sys/psTrace.h

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 23:31:44 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:47:00 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -18 +18 @@
 
 #define PS_UNKNOWN_TRACE_LEVEL -9999   // we don't know this name's level
-#define PS_DEFAULT_TRACE_LEVEL 0
+#define PS_DEFAULT_TRACE_LEVEL -1
 
 /** \addtogroup LogTrace

trunk/psLib/src/sysUtils/psTrace.h

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-27 23:31:44 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-28 22:47:00 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -18 +18 @@
 
 #define PS_UNKNOWN_TRACE_LEVEL -9999   // we don't know this name's level
-#define PS_DEFAULT_TRACE_LEVEL 0
+#define PS_DEFAULT_TRACE_LEVEL -1
 
 /** \addtogroup LogTrace