Changeset 2250

Timestamp:

Nov 1, 2004, 1:57:08 PM (22 years ago)

Author:

gusciora

Message:

Fixed the robust stats routines for stdev/mean. Fixed minimization routines
for fitting gaussians. CVS:

---

Enter Log. Lines beginning with `CVS:' are removed automatically CVS: CVS:
Committing in . CVS: CVS: Modified Files: CVS:
src/dataManip/psMinimize.c src/dataManip/psMinimize.h CVS:
src/dataManip/psStats.c src/dataManip/makedir/psMinimize.d CVS:
src/dataManip/makedir/psStats.d CVS: test/dataManip/tst_psMinimize07.c

---

Location:

trunk/psLib

Files:

• src/dataManip/psMinimize.c (modified) (25 diffs)
• src/dataManip/psMinimize.h (modified) (2 diffs)
• src/dataManip/psStats.c (modified) (4 diffs)
• src/math/psMinimize.c (modified) (25 diffs)
• src/math/psMinimize.h (modified) (2 diffs)
• src/math/psStats.c (modified) (4 diffs)
• test/dataManip/tst_psMinimize07.c (modified) (6 diffs)
• test/dataManip/tst_psStats07.c (modified) (2 diffs)

trunk/psLib/src/dataManip/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-01 18:57:05 $
+ *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-01 23:57:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -572 +572 @@
                        (currValueVec->data.F32[n] - y->data.F32[n]);
             psTrace(".psLib.dataManip.psMinimize", 6,
-                    "data[%d], chi2 calculation+= (%f * %f)\n", n,
+                    "data[%d], chi2 calculation+= (%f - %f)^2\n", n,
                     currValueVec->data.F32[n], y->data.F32[n]);
         }
@@ -1004 +1004 @@
 }
 
+// This macro takes as input the vector BASE and adds a multiple of the vector
+// LINE to it.  We assume BASEMASK is non-null.
 #define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
 for (int i=0;i<BASE->n;i++) { \
@@ -1036 +1038 @@
 and the parameter to vary, and the line along which it must vary.  This
 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).
+f(param + b * line) < f(param + a * line)
+f(param + b * line) < f(param + c * line)
+a < b < c
  
 Algorithm:
@@ -1192 +1194 @@
 
 
+#define RETURN_FINAL_BRACKET(d) \
+if (a < c) { \
+    bracket->data.F32[0] = a; \
+    bracket->data.F32[1] = b; \
+    bracket->data.F32[2] = c; \
+} else { \
+    bracket->data.F32[0] = c; \
+    bracket->data.F32[1] = b; \
+    bracket->data.F32[2] = a; \
+} \
+psTrace(".psLib.dataManip.p_psDetermineBracket", 4, \
+        "---- p_psDetermineBracket() end ----\n"); \
+psTrace(".psLib.dataManip.p_psDetermineBracket", 4, "Final bracket (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); \
+return(bracket); \
+
 #define PS_DETERMINE_BRACKET_MAX_ITERATIONS 100
 psVector *p_psDetermineBracket2(psVector *params,
@@ -1206 +1223 @@
     float fc = 0.0;
     psS32 iter = 0;
-    psVector *tmp = NULL;
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
     psBool boolLineIsNull = true;
     float prevMin = 0.0;
     int countMin = 0;
-
-    psTrace(".PS_SEG.PS_PWD.PS_FILE.__func__", 6,
-            "---- __func__() begin COOL");
-
-    psTrace(".psLib.dataManip.__func__", 4,
-            "---- p_psDetermineBracket() COOL ----\n");
 
     psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
@@ -1230 +1241 @@
     }
 
-    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
-    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
-    a = -0.5;
-    b = 0.0;
-    c = 0.5;
-
-    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
-    fb = func(tmp, coords);
-
+    // We determine in what x-direction does the function decrease.
+    a = 0.0;
+    fa = func(params, coords);
+    b = 0.5;
+    iter = 0;
     do {
-        a*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
-        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
-        fa = func(tmp, coords);
-    } while (fabs(fb - fa) < FLT_EPSILON);
-
-    do {
-        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
-        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
-        fc = func(tmp, coords);
-    } while (fabs(fc - fb) < FLT_EPSILON);
-
-    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
-            "initial (fa, Fb, Fc) is (%f, %f, %f)\n", fa, fb, fc);
-
-    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);
-    }
-
-    // The object of this code is to set the initial point a to 0.0, and
-    // ensure that (fb < fa).
-    if (fc < fb) {
-        a = b;
-        b = c;
-        fa = fb;
-        fb = fc;
-    } else if (fa < fb) {
-        float aSave = a;
-        float faSave = fa;
-
+        b*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
+        fb = func(tmp, coords);
+    } while ((fabs(fb - fa) < FLT_EPSILON) && (iter++ < 100));
+
+    if (fb > fa) {
         a = b;
         fa = fb;
-
-        b = aSave;
-        fb = faSave;
+        b = 0.0;
+        fb = func(params, coords);
     }
     c = b;
 
-    // At this point we have (a, b) and we know that (fa > fb).  Initially, c=b;
+    // At this point we have (a, b) and we know that (fa >= fb).  Initially, c=b;
     // We keep stretching b out further from "a" until (fc > previous fc).  If
     // that happens, then we have our bracket.
+    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
     iter = 0;
     while (iter < PS_DETERMINE_BRACKET_MAX_ITERATIONS) {
         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
                 "psDetermineBracket(): iterationA %d\n", iter);
-        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
+        c+= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE) * (c - a);
 
         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
@@ -1296 +1274 @@
 
         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
-                "IterA(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
+                "Iteration(%d) (bracket): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
 
         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);
+            RETURN_FINAL_BRACKET();
         } else {
             b = c;
@@ -1312 +1284 @@
 
         // This code maintains a count of how many times the minimum fc has
-        // statyed the same.  If it gets too high, we exit this loop.
+        // stayed the same.  If it gets too high, we exit this loop.
         if (fc == prevMin) {
             countMin++;
@@ -1320 +1292 @@
         prevMin = fc;
         if (countMin == 10) {
-            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);
+            RETURN_FINAL_BRACKET();
         }
 
@@ -1332 +1298 @@
     }
 
-    // Okay, let's try the other direction.
-
-    psFree(tmp);
     psFree(bracket);
     psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
-            "---- p_psDetermineBracket() end (NULL) ----\n");
+            "---- p_psDetermineBracket() end (NULL) (BAD) ----\n");
     return(NULL);
 }
@@ -1351 +1314 @@
 function.
 XXX: This is F32 only
+XXX: Since this is an internal function, many of the parameter checks are
+     redundant.
  *****************************************************************************/
+#define PS_LINEMIN_MAX_ITERATIONS 30
 float p_psLineMin(psMinimization *min,
                   psVector *params,
@@ -1359 +1325 @@
                   psMinimizePowellFunc func)
 {
+    PS_PTR_CHECK_NULL(min, NAN);
+    PS_VECTOR_CHECK_NULL(params, NAN);
+    PS_VECTOR_CHECK_EMPTY(params, NAN);
+    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
+    PS_VECTOR_CHECK_NULL(line, NAN);
+    PS_VECTOR_CHECK_EMPTY(line, NAN);
+    PS_VECTOR_CHECK_TYPE(line, PS_TYPE_F32, NAN);
+    PS_VECTOR_CHECK_NULL(paramMask, NAN);
+    PS_VECTOR_CHECK_EMPTY(paramMask, NAN);
+    PS_VECTOR_CHECK_TYPE(paramMask, PS_TYPE_U8, NAN);
+    PS_PTR_CHECK_NULL(coords, NAN);
+    PS_PTR_CHECK_NULL(func, NAN);
     psVector *bracket;
     float a = 0.0;
@@ -1369 +1347 @@
     float fn = 0.0;
     float mul = 0.0;
-    psVector *tmpa = NULL;
-    psVector *tmpb = NULL;
-    psVector *tmpc = NULL;
-    psVector *tmpn = NULL;
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpc, params->n, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpn, params->n, PS_TYPE_F32);
     psS32 i = 0;
     psS32 boolLineIsNull = true;
-
-    psTrace(".psLib.dataManip.p_psLineMin", 4,
-            "---- p_psLineMin() begin ----\n");
-    psTrace(".psLib.dataManip.p_psLineMin", 6,
-            "min->maxIter is %d\n", min->maxIter);
-    psTrace(".psLib.dataManip.p_psLineMin", 6,
-            "min->tol is %f\n", min->tol);
-
-    if (paramMask != NULL) {
-        for (i=0;i<params->n;i++) {
-            if (paramMask->data.U8[i] == 0) {
-                if (line->data.F32[i] >= FLT_EPSILON) {
-                    boolLineIsNull = false;
-                    psTrace(".psLib.dataManip.p_psLineMin", 4,
-                            "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
-                }
-            }
-        }
-    } else {
-        for (i=0;i<params->n;i++) {
-            if (line->data.F32[i] >= FLT_EPSILON) {
-                boolLineIsNull = false;
-                psTrace(".psLib.dataManip.p_psLineMin", 4,
-                        "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
-            }
-        }
-    }
+    psS32 numIterations = 0;
+
+    psTrace(".psLib.dataManip.p_psLineMin", 4, "---- p_psLineMin() begin ----\n");
+    PS_VECTOR_F32_CHECK_ZERO_VECTOR(line, boolLineIsNull);
 
     if (boolLineIsNull == true) {
+        min->value = func(params, coords);
         psTrace(".psLib.dataManip.p_psLineMin", 2,
-                "p_psLineMin() called with zero line vector.\n");
+                "p_psLineMin() called with zero line vector.  Return 0.0.  Function valus is %f\n", min->value);
         return(0.0);
     }
-
-    tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
-    tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
-    tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
-    tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
 
     for (i=0;i<params->n;i++) {
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "params, paramMask, line [%d] is (%f %d %f)\n", i,
+                "(params, paramMask, line)[%d] is (%f %d %f)\n", i,
                 params->data.F32[i],
                 paramMask->data.U8[i],
@@ -1424 +1375 @@
     bracket = p_psDetermineBracket2(params, line, paramMask, coords, func);
     if (bracket == NULL) {
-        psFree(tmpa);
-        psFree(tmpb);
-        psFree(tmpc);
-        psFree(tmpn);
         psError(__func__, "(1) Could not bracket minimum.");
         return(NAN);
     }
-
-    min->iter = 0;
-    while (min->iter < min->maxIter) {
-        min->iter++;
+    numIterations = 0;
+    while (numIterations < PS_LINEMIN_MAX_ITERATIONS) {
+        numIterations++;
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "p_psLineMin(): iteration %d\n", min->iter);
+                "p_psLineMin(): iteration %d\n", numIterations);
 
         a = bracket->data.F32[0];
         b = bracket->data.F32[1];
         c = bracket->data.F32[2];
-
-        for (i=0;i<params->n;i++) {
-            if (paramMask->data.U8[i] == 0) {
-                tmpa->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
-                tmpb->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
-                tmpc->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
-            } else {
-                tmpa->data.F32[i] = params->data.F32[i];
-                tmpb->data.F32[i] = params->data.F32[i];
-                tmpc->data.F32[i] = params->data.F32[i];
-            }
-        }
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpa, a);
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpb, b);
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpc, c);
         fa = func(tmpa, coords);
         fb = func(tmpb, coords);
         fc = func(tmpc, coords);
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter,
-                a, b, c, fa, fb, fc);
+                "LineMin: f(%f %f %f) is (%f %f %f)\n", a, b, c, fa, fb, fc);
 
         // We determine which is the biggest segment in [a,b,c] then split
@@ -1495 +1431 @@
         }
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
+                "LineMin: new bracket is (%f %f %f)\n", bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
 
         mul = bracket->data.F32[1];
         if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) {
-            for (i=0;i<params->n;i++) {
-                if (paramMask->data.U8[i] == 0) {
-                    params->data.F32[i]+= mul * line->data.F32[i];
-                }
-            }
+            PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
             min->value = func(params, coords);
             psFree(bracket);
-            psFree(tmpa);
-            psFree(tmpb);
-            psFree(tmpc);
-            psFree(tmpn);
             psTrace(".psLib.dataManip.p_psLineMin", 4,
-                    "---- p_psLineMin() end.a (%f) ----\n", mul);
+                    "---- p_psLineMin() end.a (%f) (%f) ----\n", mul, min->value);
             return(mul);
         }
     }
 
+    mul = bracket->data.F32[1];
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
+    min->value = func(params, coords);
+    psTrace(".psLib.dataManip.p_psLineMin", 4,
+            "---- p_psLineMin() end.b (%f) %f ----\n", mul, min->value);
+
     psFree(bracket);
-    psFree(tmpa);
-    psFree(tmpb);
-    psFree(tmpc);
-    psFree(tmpn);
-
-    psTrace(".psLib.dataManip.p_psLineMin", 4,
-            "---- p_psLineMin() end.b (0.0) ----\n");
-    return(0.0);
+    return(mul);
 }
 
@@ -1560 +1487 @@
 
     psS32 numDims = params->n;
+    PS_VECTOR_GEN_STATIC_RECYCLED(pQP, numDims, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(u, numDims, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(Q, numDims, PS_TYPE_F32);
     psS32 i = 0;
     psS32 j = 0;
     psVector **v = NULL;
-    psVector *pQP = NULL;
-    psVector *Q = NULL;
-    psVector *u = NULL;
     psVector *myParamMask = NULL;  // XXX: Use a recylced vector here.
     psMinimization dummyMin;
@@ -1586 +1513 @@
             "min->tol is %f\n", min->tol);
 
-    pQP = psVectorAlloc(numDims, PS_TYPE_F32);
-    u = psVectorAlloc(numDims, PS_TYPE_F32);
-    Q = psVectorAlloc(numDims, PS_TYPE_F32);
-
     if (paramMask == NULL) {
-        myParamMask = psVectorAlloc(params->n, PS_TYPE_U8);
+        myParamMask = psVectorRecycle(myParamMask, params->n, PS_TYPE_U8);
+        p_psMemSetPersistent(myParamMask, true);
+        p_psMemSetPersistent(myParamMask->data.V, true);
         for (i=0;i<myParamMask->n;i++) {
             myParamMask->data.U8[i] = 0;
@@ -1643 +1568 @@
                 if (isnan(mul)) {
                     psError(__func__, "Could not perform line minimization (1).\n");
-                    psFree(pQP);
-                    psFree(u);
-                    psFree(Q);
-                    if (paramMask == NULL) {
-                        psFree(myParamMask);
-                    }
-                    for (i=0;i<numDims;i++) {
-                        psFree(v[i]);
+                    for (j=0;j<numDims;j++) {
+                        psFree(v[j]);
                     }
                     psFree(v);
                     return(false);
                 }
+                psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                        "LineMin along dimension %d has multiple %f\n", i, mul);
+
                 if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
                     biggestDiff = fabs(dummyMin.value - currFuncVal);
@@ -1664 +1586 @@
         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++) {
@@ -1670 +1593 @@
 
                 psTrace(".psLib.dataManip.psMinimizePowell", 6,
-                        "u=Q-P (%f = %f - %f)\n", u->data.F32[i],
+                        "u[i]=Q[i]-P[i] (%f = %f - %f)\n", u->data.F32[i],
                         Q->data.F32[i],
                         params->data.F32[i]);
@@ -1688 +1611 @@
         if (isnan(mul)) {
             psError(__func__, "Could not perform line minimization. (2)\n");
-            psFree(pQP);
-            psFree(u);
-            psFree(Q);
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             for (i=0;i<numDims;i++) {
                 psFree(v[i]);
@@ -1707 +1624 @@
             }
             psFree(v);
-            psFree(u);
-            psFree(pQP);
-            psFree(Q);
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
-                    "---- psMinimizePowell() end (true) ----\n");
+                    "---- psMinimizePowell() end (1)(true) ----\n");
             return(true);
         }
@@ -1754 +1665 @@
             }
             psFree(v);
-            psFree(u);
-            psFree(pQP);
-            psFree(Q);
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
-                    "---- psMinimizePowell() end (true) ----\n");
+                    "---- psMinimizePowell() end (2) (true) ----\n");
             return(true);
         }
     }
 
-    if (paramMask == NULL) {
-        psFree(myParamMask);
-    }
     for (i=0;i<numDims;i++) {
         psFree(v[i]);
     }
     psFree(v);
-    psFree(u);
-    psFree(pQP);
-    psFree(Q);
     min->iter = iterationNumber;
     psTrace(".psLib.dataManip.psMinimizePowell", 4,
-            "---- psMinimizePowell() end (false) ----\n");
+            "---- psMinimizePowell() end (0) (false) ----\n");
     return(false);
+}
+
+
+/******************************************************************************
+XXX: We assume unnormalized gaussians.
+XXX: Currently, yErr is ignored.
+ *****************************************************************************/
+psVector *psMinimizePowellChi2Gauss1D(const psVector *params,
+                                      const psArray *coords)
+{
+    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);
+
+    for (i=0;i<coords->n;i++) {
+        x = ((psVector *) (coords->data[i]))->data.F32[0];
+        out->data.F32[i] = psGaussian(x, mean, stdev, false);
+    }
+
+    return(out);
 }
 
@@ -1822 +1746 @@
     psFree(tmp);
     psTrace(".psLib.dataManip.myPowellChi2Func", 4,
-            "---- myPowellChi2Func() end (%f) ----\n", chi2);
+            "---- myPowellChi2Func() end (chi2 is %f) ----\n", chi2);
     return(chi2);
 }

trunk/psLib/src/dataManip/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 22:46:57 $
+ *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-01 23:57:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -104 +104 @@
                                   const psArray *coords);
 
+psVector *psMinimizePowellChi2Gauss1D(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.80 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 22:46:57 $
+ *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-01 23:57:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -790 +790 @@
     psFree(tmpMask);
 }
+
+void p_psNormalizeVectorF32_0(psVector* myData)
+{
+    float min = (float)HUGE;
+    float max = (float)-HUGE;
+    psS32 i = 0;
+
+    for (i = 0; i < myData->n; i++) {
+        if (myData->data.F32[i] < min) {
+            min = myData->data.F32[i];
+        }
+        if (myData->data.F32[i] > max) {
+            max = myData->data.F32[i];
+        }
+    }
+
+    //  float range = max - min;
+    //    for (i = 0; i < myData->n; i++) {
+    //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
+    //    }
+    for (i = 0; i < myData->n; i++) {
+        myData->data.F32[i] = (myData->data.F32[i] - min) / (max - min);
+    }
+}
+
+
 
 /*****************************************************************************
@@ -1162 +1188 @@
     myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
 
+
+    // Using the above (myMean, myStdev) as initial estimates, we fit a
+    // Gaussian to the robustHistogramVector.
+    psMinimization *min = psMinimizationAlloc(100, 0.1);
+    psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
+    psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
+    psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
+
+    p_psNormalizeVectorF32_0(robustHistogramVector);
+    for (i=0;i<robustHistogramVector->n;i++) {
+        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
+        ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
+        y->data.F32[i] = robustHistogramVector->data.F32[i];
+    }
+
+    myParams->data.F32[0] = myMean;
+    myParams->data.F32[1] = myStdev;
+    psMinimizeLMChi2(min,
+                     NULL,
+                     myParams,
+                     NULL,
+                     myCoords,
+                     y,
+                     NULL,
+                     (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
+    psFree(min);
+    psFree(myParams);
+    psFree(myCoords);
+    psFree(y);
     /**************************************************************************
     Set the appropriate members in the output stats struct.
     **************************************************************************/
     if (stats->options & PS_STAT_ROBUST_MEAN) {
-        stats->robustMean = myMean;
+        if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
+            printf("WARNING: the fitted Gaussian has more than 10% error for the mean.\n");
+            printf("Using the calculated mean (calc, fit) is (%f, %f)\n", myMean, myParams->data.F32[0]);
+            stats->robustMean = myMean;
+        } else {
+            stats->robustMean = myParams->data.F32[0];
+        }
     }
 
@@ -1174 +1235 @@
 
     if (stats->options & PS_STAT_ROBUST_STDEV) {
-        stats->robustStdev = myStdev;
+        if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
+            printf("WARNING: the fitted Gaussian has more than 10% error for the stdev.\n");
+            printf("Using the calculated stdev (calc, fit) is (%f, %f)\n", myStdev, myParams->data.F32[1]);
+            stats->robustStdev = myStdev;
+        } else {
+            stats->robustStdev = myParams->data.F32[1];
+        }
     }
 

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.79 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-01 18:57:05 $
+ *  @version $Revision: 1.80 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-01 23:57:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -572 +572 @@
                        (currValueVec->data.F32[n] - y->data.F32[n]);
             psTrace(".psLib.dataManip.psMinimize", 6,
-                    "data[%d], chi2 calculation+= (%f * %f)\n", n,
+                    "data[%d], chi2 calculation+= (%f - %f)^2\n", n,
                     currValueVec->data.F32[n], y->data.F32[n]);
         }
@@ -1004 +1004 @@
 }
 
+// This macro takes as input the vector BASE and adds a multiple of the vector
+// LINE to it.  We assume BASEMASK is non-null.
 #define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \
 for (int i=0;i<BASE->n;i++) { \
@@ -1036 +1038 @@
 and the parameter to vary, and the line along which it must vary.  This
 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).
+f(param + b * line) < f(param + a * line)
+f(param + b * line) < f(param + c * line)
+a < b < c
  
 Algorithm:
@@ -1192 +1194 @@
 
 
+#define RETURN_FINAL_BRACKET(d) \
+if (a < c) { \
+    bracket->data.F32[0] = a; \
+    bracket->data.F32[1] = b; \
+    bracket->data.F32[2] = c; \
+} else { \
+    bracket->data.F32[0] = c; \
+    bracket->data.F32[1] = b; \
+    bracket->data.F32[2] = a; \
+} \
+psTrace(".psLib.dataManip.p_psDetermineBracket", 4, \
+        "---- p_psDetermineBracket() end ----\n"); \
+psTrace(".psLib.dataManip.p_psDetermineBracket", 4, "Final bracket (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); \
+return(bracket); \
+
 #define PS_DETERMINE_BRACKET_MAX_ITERATIONS 100
 psVector *p_psDetermineBracket2(psVector *params,
@@ -1206 +1223 @@
     float fc = 0.0;
     psS32 iter = 0;
-    psVector *tmp = NULL;
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32);
     psBool boolLineIsNull = true;
     float prevMin = 0.0;
     int countMin = 0;
-
-    psTrace(".PS_SEG.PS_PWD.PS_FILE.__func__", 6,
-            "---- __func__() begin COOL");
-
-    psTrace(".psLib.dataManip.__func__", 4,
-            "---- p_psDetermineBracket() COOL ----\n");
 
     psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
@@ -1230 +1241 @@
     }
 
-    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
-    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
-    a = -0.5;
-    b = 0.0;
-    c = 0.5;
-
-    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
-    fb = func(tmp, coords);
-
+    // We determine in what x-direction does the function decrease.
+    a = 0.0;
+    fa = func(params, coords);
+    b = 0.5;
+    iter = 0;
     do {
-        a*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
-        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a);
-        fa = func(tmp, coords);
-    } while (fabs(fb - fa) < FLT_EPSILON);
-
-    do {
-        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
-        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
-        fc = func(tmp, coords);
-    } while (fabs(fc - fb) < FLT_EPSILON);
-
-    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
-            "initial (fa, Fb, Fc) is (%f, %f, %f)\n", fa, fb, fc);
-
-    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);
-    }
-
-    // The object of this code is to set the initial point a to 0.0, and
-    // ensure that (fb < fa).
-    if (fc < fb) {
-        a = b;
-        b = c;
-        fa = fb;
-        fb = fc;
-    } else if (fa < fb) {
-        float aSave = a;
-        float faSave = fa;
-
+        b*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b);
+        fb = func(tmp, coords);
+    } while ((fabs(fb - fa) < FLT_EPSILON) && (iter++ < 100));
+
+    if (fb > fa) {
         a = b;
         fa = fb;
-
-        b = aSave;
-        fb = faSave;
+        b = 0.0;
+        fb = func(params, coords);
     }
     c = b;
 
-    // At this point we have (a, b) and we know that (fa > fb).  Initially, c=b;
+    // At this point we have (a, b) and we know that (fa >= fb).  Initially, c=b;
     // We keep stretching b out further from "a" until (fc > previous fc).  If
     // that happens, then we have our bracket.
+    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
     iter = 0;
     while (iter < PS_DETERMINE_BRACKET_MAX_ITERATIONS) {
         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
                 "psDetermineBracket(): iterationA %d\n", iter);
-        c*= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE);
+        c+= (1.0 + PS_DETERMINE_BRACKET_STEP_SIZE) * (c - a);
 
         PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c);
@@ -1296 +1274 @@
 
         psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
-                "IterA(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
+                "Iteration(%d) (bracket): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
 
         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);
+            RETURN_FINAL_BRACKET();
         } else {
             b = c;
@@ -1312 +1284 @@
 
         // This code maintains a count of how many times the minimum fc has
-        // statyed the same.  If it gets too high, we exit this loop.
+        // stayed the same.  If it gets too high, we exit this loop.
         if (fc == prevMin) {
             countMin++;
@@ -1320 +1292 @@
         prevMin = fc;
         if (countMin == 10) {
-            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);
+            RETURN_FINAL_BRACKET();
         }
 
@@ -1332 +1298 @@
     }
 
-    // Okay, let's try the other direction.
-
-    psFree(tmp);
     psFree(bracket);
     psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
-            "---- p_psDetermineBracket() end (NULL) ----\n");
+            "---- p_psDetermineBracket() end (NULL) (BAD) ----\n");
     return(NULL);
 }
@@ -1351 +1314 @@
 function.
 XXX: This is F32 only
+XXX: Since this is an internal function, many of the parameter checks are
+     redundant.
  *****************************************************************************/
+#define PS_LINEMIN_MAX_ITERATIONS 30
 float p_psLineMin(psMinimization *min,
                   psVector *params,
@@ -1359 +1325 @@
                   psMinimizePowellFunc func)
 {
+    PS_PTR_CHECK_NULL(min, NAN);
+    PS_VECTOR_CHECK_NULL(params, NAN);
+    PS_VECTOR_CHECK_EMPTY(params, NAN);
+    PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN);
+    PS_VECTOR_CHECK_NULL(line, NAN);
+    PS_VECTOR_CHECK_EMPTY(line, NAN);
+    PS_VECTOR_CHECK_TYPE(line, PS_TYPE_F32, NAN);
+    PS_VECTOR_CHECK_NULL(paramMask, NAN);
+    PS_VECTOR_CHECK_EMPTY(paramMask, NAN);
+    PS_VECTOR_CHECK_TYPE(paramMask, PS_TYPE_U8, NAN);
+    PS_PTR_CHECK_NULL(coords, NAN);
+    PS_PTR_CHECK_NULL(func, NAN);
     psVector *bracket;
     float a = 0.0;
@@ -1369 +1347 @@
     float fn = 0.0;
     float mul = 0.0;
-    psVector *tmpa = NULL;
-    psVector *tmpb = NULL;
-    psVector *tmpc = NULL;
-    psVector *tmpn = NULL;
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpc, params->n, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(tmpn, params->n, PS_TYPE_F32);
     psS32 i = 0;
     psS32 boolLineIsNull = true;
-
-    psTrace(".psLib.dataManip.p_psLineMin", 4,
-            "---- p_psLineMin() begin ----\n");
-    psTrace(".psLib.dataManip.p_psLineMin", 6,
-            "min->maxIter is %d\n", min->maxIter);
-    psTrace(".psLib.dataManip.p_psLineMin", 6,
-            "min->tol is %f\n", min->tol);
-
-    if (paramMask != NULL) {
-        for (i=0;i<params->n;i++) {
-            if (paramMask->data.U8[i] == 0) {
-                if (line->data.F32[i] >= FLT_EPSILON) {
-                    boolLineIsNull = false;
-                    psTrace(".psLib.dataManip.p_psLineMin", 4,
-                            "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
-                }
-            }
-        }
-    } else {
-        for (i=0;i<params->n;i++) {
-            if (line->data.F32[i] >= FLT_EPSILON) {
-                boolLineIsNull = false;
-                psTrace(".psLib.dataManip.p_psLineMin", 4,
-                        "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
-            }
-        }
-    }
+    psS32 numIterations = 0;
+
+    psTrace(".psLib.dataManip.p_psLineMin", 4, "---- p_psLineMin() begin ----\n");
+    PS_VECTOR_F32_CHECK_ZERO_VECTOR(line, boolLineIsNull);
 
     if (boolLineIsNull == true) {
+        min->value = func(params, coords);
         psTrace(".psLib.dataManip.p_psLineMin", 2,
-                "p_psLineMin() called with zero line vector.\n");
+                "p_psLineMin() called with zero line vector.  Return 0.0.  Function valus is %f\n", min->value);
         return(0.0);
     }
-
-    tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
-    tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
-    tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
-    tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
 
     for (i=0;i<params->n;i++) {
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "params, paramMask, line [%d] is (%f %d %f)\n", i,
+                "(params, paramMask, line)[%d] is (%f %d %f)\n", i,
                 params->data.F32[i],
                 paramMask->data.U8[i],
@@ -1424 +1375 @@
     bracket = p_psDetermineBracket2(params, line, paramMask, coords, func);
     if (bracket == NULL) {
-        psFree(tmpa);
-        psFree(tmpb);
-        psFree(tmpc);
-        psFree(tmpn);
         psError(__func__, "(1) Could not bracket minimum.");
         return(NAN);
     }
-
-    min->iter = 0;
-    while (min->iter < min->maxIter) {
-        min->iter++;
+    numIterations = 0;
+    while (numIterations < PS_LINEMIN_MAX_ITERATIONS) {
+        numIterations++;
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "p_psLineMin(): iteration %d\n", min->iter);
+                "p_psLineMin(): iteration %d\n", numIterations);
 
         a = bracket->data.F32[0];
         b = bracket->data.F32[1];
         c = bracket->data.F32[2];
-
-        for (i=0;i<params->n;i++) {
-            if (paramMask->data.U8[i] == 0) {
-                tmpa->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
-                tmpb->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
-                tmpc->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
-            } else {
-                tmpa->data.F32[i] = params->data.F32[i];
-                tmpb->data.F32[i] = params->data.F32[i];
-                tmpc->data.F32[i] = params->data.F32[i];
-            }
-        }
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpa, a);
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpb, b);
+        PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmpc, c);
         fa = func(tmpa, coords);
         fb = func(tmpb, coords);
         fc = func(tmpc, coords);
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter,
-                a, b, c, fa, fb, fc);
+                "LineMin: f(%f %f %f) is (%f %f %f)\n", a, b, c, fa, fb, fc);
 
         // We determine which is the biggest segment in [a,b,c] then split
@@ -1495 +1431 @@
         }
         psTrace(".psLib.dataManip.p_psLineMin", 6,
-                "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
+                "LineMin: new bracket is (%f %f %f)\n", bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
 
         mul = bracket->data.F32[1];
         if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) {
-            for (i=0;i<params->n;i++) {
-                if (paramMask->data.U8[i] == 0) {
-                    params->data.F32[i]+= mul * line->data.F32[i];
-                }
-            }
+            PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
             min->value = func(params, coords);
             psFree(bracket);
-            psFree(tmpa);
-            psFree(tmpb);
-            psFree(tmpc);
-            psFree(tmpn);
             psTrace(".psLib.dataManip.p_psLineMin", 4,
-                    "---- p_psLineMin() end.a (%f) ----\n", mul);
+                    "---- p_psLineMin() end.a (%f) (%f) ----\n", mul, min->value);
             return(mul);
         }
     }
 
+    mul = bracket->data.F32[1];
+    PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
+    min->value = func(params, coords);
+    psTrace(".psLib.dataManip.p_psLineMin", 4,
+            "---- p_psLineMin() end.b (%f) %f ----\n", mul, min->value);
+
     psFree(bracket);
-    psFree(tmpa);
-    psFree(tmpb);
-    psFree(tmpc);
-    psFree(tmpn);
-
-    psTrace(".psLib.dataManip.p_psLineMin", 4,
-            "---- p_psLineMin() end.b (0.0) ----\n");
-    return(0.0);
+    return(mul);
 }
 
@@ -1560 +1487 @@
 
     psS32 numDims = params->n;
+    PS_VECTOR_GEN_STATIC_RECYCLED(pQP, numDims, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(u, numDims, PS_TYPE_F32);
+    PS_VECTOR_GEN_STATIC_RECYCLED(Q, numDims, PS_TYPE_F32);
     psS32 i = 0;
     psS32 j = 0;
     psVector **v = NULL;
-    psVector *pQP = NULL;
-    psVector *Q = NULL;
-    psVector *u = NULL;
     psVector *myParamMask = NULL;  // XXX: Use a recylced vector here.
     psMinimization dummyMin;
@@ -1586 +1513 @@
             "min->tol is %f\n", min->tol);
 
-    pQP = psVectorAlloc(numDims, PS_TYPE_F32);
-    u = psVectorAlloc(numDims, PS_TYPE_F32);
-    Q = psVectorAlloc(numDims, PS_TYPE_F32);
-
     if (paramMask == NULL) {
-        myParamMask = psVectorAlloc(params->n, PS_TYPE_U8);
+        myParamMask = psVectorRecycle(myParamMask, params->n, PS_TYPE_U8);
+        p_psMemSetPersistent(myParamMask, true);
+        p_psMemSetPersistent(myParamMask->data.V, true);
         for (i=0;i<myParamMask->n;i++) {
             myParamMask->data.U8[i] = 0;
@@ -1643 +1568 @@
                 if (isnan(mul)) {
                     psError(__func__, "Could not perform line minimization (1).\n");
-                    psFree(pQP);
-                    psFree(u);
-                    psFree(Q);
-                    if (paramMask == NULL) {
-                        psFree(myParamMask);
-                    }
-                    for (i=0;i<numDims;i++) {
-                        psFree(v[i]);
+                    for (j=0;j<numDims;j++) {
+                        psFree(v[j]);
                     }
                     psFree(v);
                     return(false);
                 }
+                psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                        "LineMin along dimension %d has multiple %f\n", i, mul);
+
                 if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
                     biggestDiff = fabs(dummyMin.value - currFuncVal);
@@ -1664 +1586 @@
         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++) {
@@ -1670 +1593 @@
 
                 psTrace(".psLib.dataManip.psMinimizePowell", 6,
-                        "u=Q-P (%f = %f - %f)\n", u->data.F32[i],
+                        "u[i]=Q[i]-P[i] (%f = %f - %f)\n", u->data.F32[i],
                         Q->data.F32[i],
                         params->data.F32[i]);
@@ -1688 +1611 @@
         if (isnan(mul)) {
             psError(__func__, "Could not perform line minimization. (2)\n");
-            psFree(pQP);
-            psFree(u);
-            psFree(Q);
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             for (i=0;i<numDims;i++) {
                 psFree(v[i]);
@@ -1707 +1624 @@
             }
             psFree(v);
-            psFree(u);
-            psFree(pQP);
-            psFree(Q);
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
-                    "---- psMinimizePowell() end (true) ----\n");
+                    "---- psMinimizePowell() end (1)(true) ----\n");
             return(true);
         }
@@ -1754 +1665 @@
             }
             psFree(v);
-            psFree(u);
-            psFree(pQP);
-            psFree(Q);
-            if (paramMask == NULL) {
-                psFree(myParamMask);
-            }
             min->iter = iterationNumber;
             psTrace(".psLib.dataManip.psMinimizePowell", 4,
-                    "---- psMinimizePowell() end (true) ----\n");
+                    "---- psMinimizePowell() end (2) (true) ----\n");
             return(true);
         }
     }
 
-    if (paramMask == NULL) {
-        psFree(myParamMask);
-    }
     for (i=0;i<numDims;i++) {
         psFree(v[i]);
     }
     psFree(v);
-    psFree(u);
-    psFree(pQP);
-    psFree(Q);
     min->iter = iterationNumber;
     psTrace(".psLib.dataManip.psMinimizePowell", 4,
-            "---- psMinimizePowell() end (false) ----\n");
+            "---- psMinimizePowell() end (0) (false) ----\n");
     return(false);
+}
+
+
+/******************************************************************************
+XXX: We assume unnormalized gaussians.
+XXX: Currently, yErr is ignored.
+ *****************************************************************************/
+psVector *psMinimizePowellChi2Gauss1D(const psVector *params,
+                                      const psArray *coords)
+{
+    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);
+
+    for (i=0;i<coords->n;i++) {
+        x = ((psVector *) (coords->data[i]))->data.F32[0];
+        out->data.F32[i] = psGaussian(x, mean, stdev, false);
+    }
+
+    return(out);
 }
 
@@ -1822 +1746 @@
     psFree(tmp);
     psTrace(".psLib.dataManip.myPowellChi2Func", 4,
-            "---- myPowellChi2Func() end (%f) ----\n", chi2);
+            "---- myPowellChi2Func() end (chi2 is %f) ----\n", chi2);
     return(chi2);
 }

trunk/psLib/src/math/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 22:46:57 $
+ *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-01 23:57:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -104 +104 @@
                                   const psArray *coords);
 
+psVector *psMinimizePowellChi2Gauss1D(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.80 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-28 22:46:57 $
+ *  @version $Revision: 1.81 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-11-01 23:57:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -790 +790 @@
     psFree(tmpMask);
 }
+
+void p_psNormalizeVectorF32_0(psVector* myData)
+{
+    float min = (float)HUGE;
+    float max = (float)-HUGE;
+    psS32 i = 0;
+
+    for (i = 0; i < myData->n; i++) {
+        if (myData->data.F32[i] < min) {
+            min = myData->data.F32[i];
+        }
+        if (myData->data.F32[i] > max) {
+            max = myData->data.F32[i];
+        }
+    }
+
+    //  float range = max - min;
+    //    for (i = 0; i < myData->n; i++) {
+    //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
+    //    }
+    for (i = 0; i < myData->n; i++) {
+        myData->data.F32[i] = (myData->data.F32[i] - min) / (max - min);
+    }
+}
+
+
 
 /*****************************************************************************
@@ -1162 +1188 @@
     myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
 
+
+    // Using the above (myMean, myStdev) as initial estimates, we fit a
+    // Gaussian to the robustHistogramVector.
+    psMinimization *min = psMinimizationAlloc(100, 0.1);
+    psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
+    psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
+    psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
+
+    p_psNormalizeVectorF32_0(robustHistogramVector);
+    for (i=0;i<robustHistogramVector->n;i++) {
+        myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
+        ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
+        y->data.F32[i] = robustHistogramVector->data.F32[i];
+    }
+
+    myParams->data.F32[0] = myMean;
+    myParams->data.F32[1] = myStdev;
+    psMinimizeLMChi2(min,
+                     NULL,
+                     myParams,
+                     NULL,
+                     myCoords,
+                     y,
+                     NULL,
+                     (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
+    psFree(min);
+    psFree(myParams);
+    psFree(myCoords);
+    psFree(y);
     /**************************************************************************
     Set the appropriate members in the output stats struct.
     **************************************************************************/
     if (stats->options & PS_STAT_ROBUST_MEAN) {
-        stats->robustMean = myMean;
+        if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
+            printf("WARNING: the fitted Gaussian has more than 10% error for the mean.\n");
+            printf("Using the calculated mean (calc, fit) is (%f, %f)\n", myMean, myParams->data.F32[0]);
+            stats->robustMean = myMean;
+        } else {
+            stats->robustMean = myParams->data.F32[0];
+        }
     }
 
@@ -1174 +1235 @@
 
     if (stats->options & PS_STAT_ROBUST_STDEV) {
-        stats->robustStdev = myStdev;
+        if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
+            printf("WARNING: the fitted Gaussian has more than 10% error for the stdev.\n");
+            printf("Using the calculated stdev (calc, fit) is (%f, %f)\n", myStdev, myParams->data.F32[1]);
+            stats->robustStdev = myStdev;
+        } else {
+            stats->robustStdev = myParams->data.F32[1];
+        }
     }
 

trunk/psLib/test/dataManip/tst_psMinimize07.c

@@ -10 +10 @@
 #define N 24
 #define NUM_PARAMS 2
-#define MEAN 8.5
-#define STDEV 1.5
+#define MEAN 14.5
+#define STDEV 1.2
 
 float expectedParm[NUM_PARAMS];
 psS32 testStatus = true;
 
-psVector *genGaussianVector()
+// This routine generates a vector of length "n" with the specified mean and
+// stdev.
+psVector *genGaussianVector(int n, float mean, float stdev)
 {
-    psVector *y = psVectorAlloc(N, PS_TYPE_F32);
+    psVector *y = psVectorAlloc(n, PS_TYPE_F32);
 
-    for (psS32 i=0;i<N;i++) {
-        y->data.F32[i] = psGaussian((float) i, MEAN, STDEV, false);
+    for (psS32 i=0;i<n;i++) {
+        y->data.F32[i] = psGaussian((float) i, mean, stdev, false);
     }
     return(y);
 }
 
-psVector *myFunc(psVector *params,
-                 psArray *coords)
-{
-    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);
-
-    for (i=0;i<coords->n;i++) {
-        x = ((psVector *) (coords->data[i]))->data.F32[0];
-        out->data.F32[i] = psGaussian(x, mean, stdev, false);
-    }
-
-    return(out);
-}
-
+// This routine tries to fit a Gaussian to a set of data points via the
+// psMinimizeLMChi2() and psMinimizeLMChi2Gauss1D() functions.
 psS32 t00()
 {
@@ -56 +40 @@
     psMinimization *min;
 
-    psTraceSetLevel(".psLib.dataManip.psMinimize", 6);
-    psTraceSetLevel(".psLib.dataManip.p_psDetermineBracket", 0);
+    psTraceSetLevel(".psLib.dataManip.psMinimize", 0);
     /**************************************************************************
      *************************************************************************/
@@ -63 +46 @@
     myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32);
     myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
-    myParams->data.F32[0] = (float) (N/2);
-    myParams->data.F32[1] = 0.5;
+    myParams->data.F32[0] = MEAN * 0.7;
+    myParams->data.F32[1] = STDEV * 0.7;
     myCoords = psArrayAlloc(N);
-    psVector *y = genGaussianVector();
+    psVector *y = genGaussianVector(N, MEAN, STDEV);
 
     for (i=0;i<N;i++) {
         myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
         ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
+        printf("Initial data %d: (%f, %f)\n", i, ((psVector *) (myCoords->data[i]))->data.F32[0], y->data.F32[i]);
     }
 
@@ -105 +89 @@
 }
 
+// This routine tries to fit a Gaussian to a set of data points via the
+// psMinimizeChi2Powell() function.
 psS32 t01()
 {
@@ -115 +101 @@
     psMinimization *min;
 
-    psTraceSetLevel(".psLib.dataManip.psMinimizePowell", 10);
-    psTraceSetLevel(".psLib.dataManip.myPowellChi2Func", 10);
+    psTraceSetLevel(".psLib.dataManip.psMinimize", 0);
+    psTraceSetLevel(".psLib.dataManip.psMinimizeLMChi2Gauss1D", 0);
+
+    psTraceSetLevel(".psLib.dataManip.psMinimizePowell", 0);
+    psTraceSetLevel(".psLib.dataManip.myPowellChi2Func", 0);
+    psTraceSetLevel(".psLib.dataManip.p_psLineMin", 0);
+    psTraceSetLevel(".psLib.dataManip.p_psDetermineBracket", 0);
     psTraceSetLevel(".psLib.dataManip.psFunctions.psGaussian", 0);
     psTraceSetLevel(".psLib.dataManip.psFunctions", 0);
-    psTraceSetLevel(".", 10);
-    psTracePrintLevels();
+    psTraceSetLevel(".", 0);
+    //    psTracePrintLevels();
     /**************************************************************************
      *************************************************************************/
@@ -126 +117 @@
     myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32);
     myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
-    myParams->data.F32[0] = (float) (N/2);
-    myParams->data.F32[0] = 2.0 * MEAN;
-    myParams->data.F32[1] = 0.5;
-    myParams->data.F32[1] = 1.5 * STDEV;
+    myParams->data.F32[0] = MEAN * 0.7;
+    myParams->data.F32[1] = STDEV * 0.7;
     myCoords = psArrayAlloc(N);
-    psVector *y = genGaussianVector();
+    psVector *y = genGaussianVector(N, MEAN, STDEV);
 
     for (i=0;i<N;i++) {
         myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
         ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i;
+        printf("Initial data %d: (%f, %f)\n", i, ((psVector *) (myCoords->data[i]))->data.F32[0], y->data.F32[i]);
     }
 
+    //    psMinimizeChi2Powell(min, myParams, NULL, myCoords, y, NULL,
+    //                         (psMinimizeChi2PowellFunc) myFunc);
     psMinimizeChi2Powell(min, myParams, NULL, myCoords, y, NULL,
-                         (psMinimizeChi2PowellFunc) myFunc);
+                         (psMinimizeChi2PowellFunc) psMinimizePowellChi2Gauss1D);
 
     printf("\nThe chi-squared is %f\n", min->value);

trunk/psLib/test/dataManip/tst_psStats07.c

@@ -13 +13 @@
 #define MEAN 32.0
 #define STDEV 2.0
-#define ERROR_TOLERANCE 0.10
+#define ERROR_TOLERANCE 0.15
 
 psS32 t00()
@@ -38 +38 @@
     psS32 realNfitNoMask = N / 4;
 
-    psTraceSetLevel(".psLib.dataManip.psStats", 10);
+    psTraceSetLevel(".psLib.dataManip", 0);
+    psTraceSetLevel(".psLib.dataManip.psMinimize", 0);
+    psTraceSetLevel(".psLib.dataManip.psMinimizeLMChi2Gauss1D", 0);
+    psTraceSetLevel(".psLib.dataManip.psFunctions.psGaussian", 0);
+    psTraceSetLevel(".psLib.dataManip.psFunctions", 0);
 
     /*************************************************************************/