Changeset 2197

Timestamp:

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

Author:

gusciora

Message:

...

Location:

trunk/psLib

Files:

• src/astro/psCoord.c (modified) (4 diffs)
• src/astronomy/psCoord.c (modified) (4 diffs)
• src/dataManip/psConstants.h (modified) (6 diffs)
• src/dataManip/psFunctions.c (modified) (1 diff)
• src/dataManip/psMinimize.c (modified) (2 diffs)
• src/dataManip/psMinimize.h (modified) (3 diffs)
• src/dataManip/psStats.c (modified) (29 diffs)
• src/math/psConstants.h (modified) (6 diffs)
• src/math/psMinimize.c (modified) (2 diffs)
• src/math/psMinimize.h (modified) (3 diffs)
• src/math/psPolynomial.c (modified) (1 diff)
• src/math/psSpline.c (modified) (1 diff)
• src/math/psStats.c (modified) (29 diffs)
• test/dataManip/tst_psMinimize06.c (modified) (5 diffs)
• test/dataManip/tst_psStats07.c (modified) (4 diffs)

trunk/psLib/src/astro/psCoord.c

@@ -10 +10 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-10-14 01:17:50 $
+*  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-10-26 21:24:42 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -58 +58 @@
 /*****************************************************************************/
 
-static float p_psArg(float x, float y);
-
-/******************************************************************************
-XXX: Verify this arc tan function.
-XXX: macro
-XXX: I don't think this arg() is correct.
- *****************************************************************************/
-
-float p_psArg(float x,
-              float y)
-{
-    if (y!= 0.0) {
-        return(atan(y/x));
-    } else {
-        return(0.5 * M_PI);
-    }
-    return (0.0);
-    /*
-        if (x > 0) {
-            return (atan(y / x));
-        } else if ((x == 0) && (y == 0)) {
-            return (0.5 * M_PI);
-        } else if ((x == 0) && (y == 0)) {
-            return (-0.5 * M_PI);
-        } else if ((x == 0) && (y == 0)) {
-            return (M_PI + atan(y / x));
-        } else if ((x == 0) && (y == 0)) {
-            return (-M_PI + atan(y / x));
-        }
-        psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y);
-    */
-}
-
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
@@ -292 +259 @@
     if (projection->type == PS_PROJ_TAN) {
         R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
-        tmp->d = p_psArg(-coord->y, coord->x);
+        tmp->d = atan2(-coord->y, coord->x);
         tmp->r = atan(180.0 / (R * M_PI));
 
     } else if (projection->type == PS_PROJ_SIN) {
         R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
-        tmp->d = p_psArg(-coord->y, coord->x);
+        tmp->d = atan2(-coord->y, coord->x);
         tmp->r = acos((R * M_PI) / 180.0);
 
@@ -314 +281 @@
         chu2 *= chu2;
         chu = sqrt(1.0 - chu1 - chu2);
-        tmp->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
+        tmp->d = 2.0 * atan2((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
         tmp->r = asin((coord->y * chu * M_PI) / 180.0);
 

trunk/psLib/src/astronomy/psCoord.c

@@ -10 +10 @@
 *  @author George Gusciora, MHPCC
 *
-*  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-10-14 01:17:50 $
+*  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-10-26 21:24:42 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -58 +58 @@
 /*****************************************************************************/
 
-static float p_psArg(float x, float y);
-
-/******************************************************************************
-XXX: Verify this arc tan function.
-XXX: macro
-XXX: I don't think this arg() is correct.
- *****************************************************************************/
-
-float p_psArg(float x,
-              float y)
-{
-    if (y!= 0.0) {
-        return(atan(y/x));
-    } else {
-        return(0.5 * M_PI);
-    }
-    return (0.0);
-    /*
-        if (x > 0) {
-            return (atan(y / x));
-        } else if ((x == 0) && (y == 0)) {
-            return (0.5 * M_PI);
-        } else if ((x == 0) && (y == 0)) {
-            return (-0.5 * M_PI);
-        } else if ((x == 0) && (y == 0)) {
-            return (M_PI + atan(y / x));
-        } else if ((x == 0) && (y == 0)) {
-            return (-M_PI + atan(y / x));
-        }
-        psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y);
-    */
-}
-
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
@@ -292 +259 @@
     if (projection->type == PS_PROJ_TAN) {
         R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
-        tmp->d = p_psArg(-coord->y, coord->x);
+        tmp->d = atan2(-coord->y, coord->x);
         tmp->r = atan(180.0 / (R * M_PI));
 
     } else if (projection->type == PS_PROJ_SIN) {
         R = sqrt((coord->x * coord->x) + (coord->y * coord->y));
-        tmp->d = p_psArg(-coord->y, coord->x);
+        tmp->d = atan2(-coord->y, coord->x);
         tmp->r = acos((R * M_PI) / 180.0);
 
@@ -314 +281 @@
         chu2 *= chu2;
         chu = sqrt(1.0 - chu1 - chu2);
-        tmp->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
+        tmp->d = 2.0 * atan2((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0);
         tmp->r = asin((coord->y * chu * M_PI) / 180.0);
 

trunk/psLib/src/dataManip/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-14 23:43:12 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -107 +107 @@
 } \
 
+/** Preprocessor macro to generate error on a NULL vector */
+#define PS_CHECK_NULL_VECTOR_RETURN_F0(NAME) \
+if (NAME == NULL || NAME->data.V == NULL) { \
+    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
+    return(0.0); \
+} \
+
 /** Preprocessor macro to generate error on a NULL poniter */
 #define PS_CHECK_NULL_PTR_RETURN_NULL(NAME) \
@@ -113 +120 @@
     return(NULL); \
 }
+
+/** Preprocessor macro to generate error on a NULL poniter */
+#define PS_CHECK_NULL_PTR_RETURN_F0(NAME) \
+if (NAME == NULL) { \
+    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
+    return(0.0); \
+}
+
 /** Preprocessor macro to generate error on a NULL poniter */
 #define PS_CHECK_NULL_PTR_RETURN_ZERO(NAME) \
@@ -127 +142 @@
 } \
 
+/** Preprocessor macro to generate error for zero length vector */
+#define PS_CHECK_EMPTY_VECTOR_RETURN_F0(NAME) \
+if (NAME->n < 1) { \
+    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
+    return(0.0); \
+} \
+
 /** Preprocessor macro to generate error on differing size vectors */
 #define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(VEC1, VEC2) \
@@ -132 +154 @@
     psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     return(NULL); \
+}
+
+/** Preprocessor macro to generate error on differing size vectors */
+#define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_F0(VEC1, VEC2) \
+if (VEC1->n != VEC2->n) { \
+    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
+    return(0.0); \
 }
 
@@ -199 +228 @@
 
 #define PS_MAX(A, B) \
-((A > B) ? A : B) \
+(((A) > (B)) ? (A) : (B)) \
 
 #define PS_MIN(A, B) \
-((A < B) ? A : B) \
-
+(((A) < (B)) ? (A) : (B)) \
+

trunk/psLib/src/dataManip/psFunctions.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-14 01:44:48 $
+ *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psLib/src/dataManip/psMinimize.c

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

trunk/psLib/src/dataManip/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-06 23:38:49 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:43 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -84 +84 @@
 
 
+/*
+   Bug 203:
+ 
+typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params,
+const psVector *x);
+ 
+This function takes the current guess for the parameters for which we are
+trying to get the best values (params), and a single vector (x) of
+conditions from the array of x values fed into the minimiser.
+ 
+    bool psMinimizeLMChi2(psMinimization *min,
+                          psImage *covar,
+                          psVector *params,
+                          const psVector *paramMask,
+                          const psArray *x,
+                          const psVector *y,
+                          const psVector *yErr,
+                          psMinimizeLMChi2Func func);
+ 
+This takes the minimization specs (min), returns the covariance matrix
+(covar), takes the best guess of initial parameters (params), the parameter
+mask (paramMask), and takes multiple vectors of conditions in an array (x),
+the corresponding measured values (y) and errors (yErr) and the function to
+fit (func).
+ 
+For example, for GRB afterglows, I have flux as a function of time and
+frequency.  So I stuff into the "psArray *x" all my time and frequency values
+(so a whole heap of vectors of size 2), I have the measured values in "y" and
+errors in "yErr".  Then each of the time-frequency pairs are passed to my model
+function with the current parameters, and the model function returns the flux
+for that time-frequency pair, and the derivative with respect to each of the
+parameters.
+ 
+This seems reasonable.  The only thing we could change would be to have the
+function be defined:
+ 
+typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,
+                                          const psVector *params,
+                                          const psArray *x);
+ 
+So it would return the model value for each of the measurements at once, and
+return for each the derivatives (so that it returns a matrix).
+ 
+What do you think?
+ 
+ 
+ 
+I'm not sure I understand how LM chi-squared minimization will work with the
+following function that you define:
+ 
+    typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,
+                                              const psVector *params,
+                                              const psArray *x);
+ 
+The chi-squared minimization algorithm, as defined in NR, requires that
+function to be minimized be evaluated at each data point, and that all
+derivatives, with respect to each parameter, be calculated at each data
+point.  In the above, can I assume that
+ 
+    x is an array of psVectors, with each vector corresponding to a single
+ data point.
+ 
+    The returned value has the same length as x.  It contains the value of
+ the function at each data point in x.
+ 
+    deriv: an n-by-p matrix where "n" is the number of data points, and "p"
+ is the number of parameters.  The [i][j] element of this matrix
+ holds the derivative of the function at the i-th data point with
+ respect to the j-th parameter.
+ 
+ 
+ 
+*/
+
+
 // XXX: What if any of these arguments are NULL?
+
+
 typedef
 float (*psMinimizeLMChi2Func)(psVector *deriv,
@@ -98 +175 @@
                       psMinimizeLMChi2Func func);
 
+
+typedef
+psVector* (*psMyMinimizeLMChi2Func)(psImage *deriv,
+                                    const psVector *params,
+                                    const psArray *x);
+
+bool psMyMinimizeLMChi2(psMinimization *min,
+                        psImage *covar,
+                        psVector *params,
+                        const psVector *paramMask,
+                        const psArray *x,
+                        const psVector *y,
+                        const psVector *yErr,
+                        psMyMinimizeLMChi2Func func);
+
 typedef
 float (*psMinimizePowellFunc)(const psVector *params,

trunk/psLib/src/dataManip/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.67 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-13 02:40:13 $
+ *  @version $Revision: 1.68 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:43 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -40 +40 @@
 /*****************************************************************************/
 // will use robust statistical methods.
-#define GAUSS_WIDTH 20                     // The width of the Gaussian or boxcar smoothing.
+#define GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
 #define CLIPPED_NUM_ITER_LB 1
 #define CLIPPED_NUM_ITER_UB 10
@@ -57 +57 @@
 /* TYPE DEFINITIONS                                                          */
 /*****************************************************************************/
+psVector* p_psConvertToF32(psVector* in);
 
 /*****************************************************************************/
@@ -423 +424 @@
     float rangeMax = 0.0;       // Exclude date above this
 
-    // Determine if the number of data points exceed a threshold which will
-    // cause to generate robust stats, as opposed to exact stats.
-    // XXX: This code is no longer used.  We now calculate the median exactly
-    // regardless of the vector size.
-    /*
-     * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-     * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-     * 
-     * // Store the robust quartiles into the sample quartile members. stats->sampleMedian =
-     * stats2->robustMedian;
-     * 
-     * // Free temporary data buffers. psFree(stats2);
-     * 
-     * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options |
-     * PS_STAT_ROBUST_FOR_SAMPLE;
-     * 
-     * return; } */
-
     // Determine how many data points fit inside this min/max range
     // and are not masked, IF the maskVector is not NULL>
@@ -507 +490 @@
 
 /******************************************************************************
-    This routine smoothes the data in the input robustHistogram with a
-    Gaussian of width sigma.
- 
-    XXX: Must consult with IfA on the proper values for the Gaussian
-    smoothing.  Currently, the Gaussian coefficients are such that only
-    the center coefficient is non-zero.  This is because "e" is being
-    raised to a very large negative number for all points except the
-    center point.
- 
-    XXX: use a static variable for gaussianCoefs[] and compute them once.
- *****************************************************************************/
-psVector* p_psVectorsmoothHistGaussian(psHistogram* robustHistogram,
+p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
+robustHistogram with a Gaussian of width sigma.
+ *****************************************************************************/
+#define GAUSS_WIDTH_OTHER 5.0
+psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
                                        float sigma)
 {
     int i = 0;                  // Loop index variable
     int j = 0;                  // Loop index variable
-    float denom = 0.0;          // Temporary variable
-    float expo = 0.0;           // Temporary variable
-    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
-    psVector* smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32);
-
-    for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) {
-        if (fabs(sigma) >= FLT_EPSILON) {
-            // If sigma does not equal zero, then we use Gaussian smoothing.
-            #ifdef  DARWIN
-            denom = (float)sqrt(2.0 * M_PI * sigma * sigma);
-            #else
-
-            denom = sqrtf(2.0 * M_PI * sigma * sigma);
-            #endif
-
-            expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH));
-            expo /= (2.0 * sigma * sigma);
-            gaussianCoefs[i] = exp(expo / denom);
-
-            //printf("p_psVectorsmoothHistGaussian(): gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
-            // NOTE: Gaussian smoothing just isn't working with low sigma
-            // values.  The problem is that the Gaussian coefficients are
-            // all zero, except for the middle coefficient, which is exactly
-            // one.  Therefore, I'm using boxcar smoothing.
-            gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH));
-            // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
-        } else {
-            /* If sigma equals zero (all pixels have the same value) the above code will divide by zero.
-             * Therefore, we don't need to smooth the data. */
-            for (i = 0; i < robustHistogram->nums->n; i++) {
-                smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i];
-            }
-            return (smooth);
-        }
-    }
-
-    // Perform the actual smoothing.
-    for (i = 0; i < robustHistogram->nums->n; i++) {
+    float iMid;
+    float jMid;
+    int numBins = robustHistogram->nums->n;
+    int numBounds = robustHistogram->bounds->n;
+    psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32);
+    int jMin = 0;
+    int jMax = 0;
+    float firstBound = robustHistogram->bounds->data.F32[0];
+    float lastBound = robustHistogram->bounds->data.F32[numBounds-1];
+    psScalar x;
+
+    x.type.type = PS_TYPE_F32;
+    for (i = 0; i < numBins; i++) {
+        iMid = (robustHistogram->bounds->data.F32[i] +
+                robustHistogram->bounds->data.F32[i+1]) / 2.0;
+
+        // YYY: The p_psVectorBinDisect() routine does much of the work of
+        // the following conditionals, however, it also reports a warning
+        // message.  I don't want the warning message so I reproduce the
+        // conditionals here.  Maybe p_psVectorBinDisect() should not produce
+        // warnings?
+
+        x.data.F32 = iMid - (GAUSS_WIDTH * sigma);
+        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+            jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
+        } else if (x.data.F32 <= firstBound) {
+            jMin = 0;
+        } else if (x.data.F32 >= lastBound) {
+            jMin = robustHistogram->bounds->n - 1;
+
+        }
+
+        x.data.F32 = iMid + (GAUSS_WIDTH * sigma);
+        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+            jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
+        } else if (x.data.F32 <= firstBound) {
+            jMax = 0;
+        } else if (x.data.F32 >= lastBound) {
+            jMax = robustHistogram->bounds->n - 1;
+
+        }
+
         smooth->data.F32[i] = 0.0;
-        for (j = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) {
-            if (((j + i) >= 0) && ((j + i) < smooth->n)) {
-                smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] *
-                                        (float)robustHistogram->nums->data.S32[j + i]);
-            }
-        }
-    }
-    return (smooth);
+        for (j = jMin ; j <= jMax ; j++) {
+            jMid = (robustHistogram->bounds->data.F32[j] +
+                    robustHistogram->bounds->data.F32[j+1]) / 2.0;
+            smooth->data.F32[i] +=
+                ((float) robustHistogram->nums->data.U32[j]) *
+                psGaussian(jMid, iMid, sigma, true);
+        }
+    }
+
+    return(smooth);
 }
 
@@ -913 +895 @@
 specified data point.
  *****************************************************************************/
-float p_psGaussian(const psVector* restrict myData, const psVector* restrict myParams)
+float p_psGaussian(const psVector* restrict myData,
+                   const psVector* restrict myParams)
 {
     float x = myData->data.F32[0];
@@ -930 +913 @@
 calculates the specified partial derivative of the above Gaussian function.
  *****************************************************************************/
-float p_psGaussianDeriv(const psVector* restrict myData, const psVector* restrict myParams, int whichParam)
+float p_psGaussianDeriv(const psVector* restrict myData,
+                        const psVector* restrict myParams, int whichParam)
 {
     float x = myData->data.F32[0];
@@ -1036 +1020 @@
 
 /******************************************************************************
-p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given
-the specified histogram, with the specified pre-calculated cumulativeSums,
-this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum",
-and then finds the point x such that f(x) == fitFloat.
+fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): This routine
+takes psVectors of x/y pairs as input, and fits a quadratic to the 3 points
+surrounding element binNum in the vectors (the midpoint between element i
+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: This function is currently not being used.
- *****************************************************************************/
-float p_psFitQuadratic(psHistogram* histogram,
-                       psVector* cumulativeSums,
-                       int binNum,
-                       float fitFloat)
-{
+// XXX: Use static variables.
+*****************************************************************************/
+float fitQuadraticSearchForYThenReturnX(psVector *xVec,
+                                        psVector *yVec,
+                                        int binNum,
+                                        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);
-
-    if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) {
-        x->data.F64[0] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]);
-        x->data.F64[1] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]);
-        x->data.F64[2] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]);
-
-        y->data.F64[0] = cumulativeSums->data.F32[binNum - 1];
-        y->data.F64[1] = cumulativeSums->data.F32[binNum];
-        y->data.F64[2] = cumulativeSums->data.F32[binNum + 1];
-
-        if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) {
-            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
-        }
-
+    float tmpFloat;
+
+    if ((binNum > 0) && (binNum < (yVec->n - 2))) {
+        x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
+        x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
+        x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
+        y->data.F64[0] = yVec->data.F32[binNum - 1];
+        y->data.F64[1] = yVec->data.F32[binNum];
+        y->data.F64[2] = yVec->data.F32[binNum + 1];
+
+        // Ensure that yVal is within the range of the bins we are using.
+        if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) {
+            printf("((%f), %f, %f)\n", yVal, y->data.F64[0], y->data.F64[2]);
+            psError(__func__, "yVal not within y-range\n");
+        }
         yErr->data.F64[0] = 1.0;
         yErr->data.F64[1] = 1.0;
         yErr->data.F64[2] = 1.0;
 
+        // Determine the coefficients of the polynomial.
         myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-        return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat));
+        // Call p_ps1DPolyMedian(), which does a binary search on the
+        // polynomial, looking for the value x such that f(x) = yVal
+        tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);
     } else {
-        return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum]));
+        if (binNum == 0) {
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        } else if (binNum == (xVec->n - 1)) {
+            // XXX: Is this right?
+            tmpFloat = xVec->data.F32[binNum];
+        } else if (binNum == (xVec->n - 2)) {
+            // XXX: Is this right?
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        }
     }
 
@@ -1083 +1097 @@
     psFree(yErr);
     psFree(myPoly);
-    return (0.0);
-}
+    return(tmpFloat);
+}
+
+#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
+(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
 
 /******************************************************************************
@@ -1111 +1128 @@
 *****************************************************************************/
 void p_psVectorRobustStats(const psVector* restrict myVector,
-                           const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
+                           const psVector* restrict maskVector,
+                           unsigned int maskVal,
+                           psStats* stats)
 {
     psHistogram* robustHistogram = NULL;
@@ -1118 +1137 @@
     int LQBinNum = -1;          // Bin num for lower quartile
     int UQBinNum = -1;          // Bin num for upper quartile
+    int medianBinNum = -1;
     int i = 0;                  // Loop index variable
-    int maxBinNum = 0;
-    float maxBinCount = 0.0;
+    int modeBinNum = 0;
+    float modeBinCount = 0.0;
     float dL = 0.0;
     int numBins = 0;
-    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);
-
-    // psImage* domain;
-    // psVector* errors;
-    // psVector* data;
-    // psVector* initialGuess;
-    // psVector* theParams;
-    // float chiSq=0.0;
-    // float max = -HUGE;
-    // float max_pos;
     float myMean = 0.0;
     float myStdev = 0.0;
@@ -1139 +1149 @@
     float sumSquares = 0.0;
     float sumDiffs = 0.0;
-    psVector* x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* y = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* yErr = psVectorAlloc(3, PS_TYPE_F64);
-    psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
-    psVector* cumulativeRobustSumsFullRange = NULL;
-    psVector* cumulativeRobustSumsDlRange = NULL;
+    psVector* cumulativeRobustSums = NULL;
     float sumRobust = 0.0;
     float sumN50 = 0.0;
     float sumNfit = 0.0;
-    float cumulativeMedian = 0.0;
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F32;
+    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);
 
     // Compute the initial bin size of the robust histogram.  This is done
@@ -1188 +1195 @@
         p_psVectorMax(myVector, maskVector, maskVal, stats);
     }
+
     // Create the histogram structure.  NOTE: we can not specify the bin size
     // precisely since the argument to psHistogramAlloc() is the number of
@@ -1198 +1206 @@
     psVectorHistogram(robustHistogram, myVector, maskVector, maskVal);
 
-    // Smooth the histogram.
-    // XXX: is that the right stdev?
-    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram,
+    // Smooth the histogram, Gaussian-style.
+    robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram,
                             tmpStats->clippedStdev / 4.0f);
 
@@ -1208 +1215 @@
 
     /**************************************************************************
-    Determine the lower/upper quartiles.
+    Determine the median/lower/upper quartile bin numbers.
+
+    We define a vector called "cumulativeRobustSums" where the value at
+    index position i is equal to the sum of bins 0:i.  This will be used in
+    determining the median and lower/upper quartiles.
     **************************************************************************/
-    // We define a vector called "cumulativeRobustSums..." where the value at
-    // index position i is equal to the sum of bins 0:i.  This will be used
-    // now and later in determining the lower/upper quartiles.
-    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0];
+    cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
+    cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0];
     for (i = 1; i < robustHistogramVector->n; i++) {
-        cumulativeRobustSumsFullRange->data.F32[i] =
-            cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
-    }
-    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
-
-    // Determine the bin number containing the lower quartile point.
-    LQBinNum = -1;
-    for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) {
-        if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) {
-            LQBinNum = i;
-            break;
-        }
-    }
-
-    // Determine the bin number containing the upper quartile point.
-    UQBinNum = -1;
-    for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) {
-        if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) {
-            UQBinNum = i;
-            break;
-        }
-    }
-
-    if ((LQBinNum == -1) || (UQBinNum == -1)) {
-        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
-    }
-
+        cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] +
+                                            robustHistogramVector->data.F32[i];
+    }
+    sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];
+
+    tmpScalar.data.F32 = sumRobust / 4.0;
+    LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+    tmpScalar.data.F32 = 3.0 * sumRobust / 4.0;
+    UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+    tmpScalar.data.F32 = sumRobust / 2.0;
+    medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+
+    if ((LQBinNum < 0) || (UQBinNum < 0)) {
+        psError(__func__, "Could not determine the robust lower/upper quartile bin numbers.");
+    }
+    if (medianBinNum < 0) {
+        psError(__func__, "Could not determine the robust lower/upper quartile bin numbers.");
+    }
     /**************************************************************************
     Determine the mode in the range LQ:UQ.
     **************************************************************************/
     // Determine the bin with the peak value in the range LQ to UQ.
-    maxBinNum = LQBinNum;
-    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
-    sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum];
+    modeBinNum = LQBinNum;
+    modeBinCount = robustHistogramVector->data.F32[LQBinNum];
+    sumN50 = (float)robustHistogram->nums->data.U32[LQBinNum];
     for (i = LQBinNum + 1; i <= UQBinNum; i++) {
-        if (robustHistogramVector->data.F32[i] > maxBinCount) {
-            maxBinNum = i;
-            maxBinCount = robustHistogramVector->data.F32[i];
-        }
-        sumN50 += (float)robustHistogram->nums->data.S32[i];
-    }
-
-    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
+        if (robustHistogramVector->data.F32[i] > modeBinCount) {
+            modeBinNum = i;
+            modeBinCount = robustHistogramVector->data.F32[i];
+        }
+        sumN50 += (float)robustHistogram->nums->data.U32[i];
+    }
+
     dL = (UQBinNum - LQBinNum) / 4;
-
-    printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum);
-
     /**************************************************************************
     Determine the mean/stdev for the bins in the range mode-dL to mode+dL
     **************************************************************************/
-    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    for (i = 0; i < robustHistogramVector->n; i++) {
-        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
-    }
-    sumNfit = 0.0;
-    cumulativeMedian = 0.0;
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
-        if ((0 <= i) && (i < robustHistogramVector->n)) {
-            cumulativeRobustSumsDlRange->data.F32[i] =
-                cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
-            cumulativeMedian += robustHistogramVector->data.F32[i];
-            sumNfit += (float)robustHistogram->nums->data.S32[i];
-        }
-    }
-
     // Calculate the mean of the smoothed robust histogram in the range
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     // that bin (this is a non-exact approximation).
+    sumNfit = 0.0;
     myMean = 0.0;
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
+    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
         if ((0 <= i) && (i < robustHistogramVector->n)) {
-            myMean += (robustHistogramVector->data.F32[i]) * 0.5 *
-                      (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]);
+            myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i);
             countFloat += robustHistogramVector->data.F32[i];
         }
+
+        sumNfit += (float)robustHistogram->nums->data.U32[i];
     }
     myMean /= countFloat;
@@ -1297 +1279 @@
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     // that bin.
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
+    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
         if ((0 <= i) && (i < robustHistogramVector->n)) {
-            diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] +
-                           robustHistogram->bounds->data.F32[i])) - myMean;
+            diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean;
             sumSquares += diff * diff * robustHistogramVector->data.F32[i];
             sumDiffs += diff * robustHistogramVector->data.F32[i];
@@ -1315 +1296 @@
 
     if (stats->options & PS_STAT_ROBUST_MODE) {
-        stats->robustMode = 0.5 *
-                            (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]);
+        stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum);
     }
 
@@ -1322 +1302 @@
         stats->robustStdev = myStdev;
     }
-    // To determine the median (and later, the lower/upper quartile), we fit
-    // a quadratic to the three bins surrounding the bin containing the median.
-    // The quadratic y=f(x) with x being the midpoint of each bin, and y being
-    // the cumulative number of data points in all bins up to, and including,
-    // this bin.  We then solve the quadratic for
-
-    // XXX: Create a "p_psQuadraticFitAndSolveForX() function.
-
+
+    // To determine the median, we fit a quadratic y=f(x) to the three bins
+    // surrounding the bin containing the median (x is the midpoint of each
+    // bin and y is the value of each bin).  Then we figure out what value
+    // of x corresponds to f(x) being the median (half of all points).
     if (stats->options & PS_STAT_ROBUST_MEDIAN) {
-        if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum - 1] +
-                              robustHistogram->bounds->data.F32[maxBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum] +
-                              robustHistogram->bounds->data.F32[maxBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum + 1] +
-                              robustHistogram->bounds->data.F32[maxBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum];
-            y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1];
-
-            // Ensure that cumulativeMedian/2 is actually within the range of the bins
-            // we are using.
-            cumulativeMedian *= 0.5;
-            if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) {
-                printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]);
-                psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n");
-            }
-            // XXX: yErr is not currently used by psVectorFitPolynomial1D().  We
-            // may have to set this meaningfully later.
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            // Determine the coefficients of the polynomial.
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            // Call p_ps1DPolyMedian(), which does a binary search on the
-            // polynomial, looking for the value x such that
-            // f(x) = cumulativeMedian.
-            stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian);
-        } else {
-            // If the mode is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] +
-                                         robustHistogram->bounds->data.F32[maxBinNum]);
-        }
-    }
-    // The lower/upper quartile calculations are very similar to the median
-    // calculations.  We fit a quadratic to the array containing the
-    // cumulative data points in each bin, then search for x such that
-    // f(x) equals the lower/upper quartile exactly.
-    //
+        // Take a psVector.  Fit a polynomial y = f(x) to the 3 data elements
+        // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0.
+        stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                              cumulativeRobustSums,
+                              medianBinNum,
+                              sumRobust/2.0);
+    }
+
+    // To determine the quartiles, we fit a quadratic y=f(x) to the three bins
+    // surrounding the bin containing LQ/UQ (x is the midpoint of each
+    // bin and y is the value of each bin).  Then we figure out what value
+    // of x corresponds to f(x) being the LQ/UQ.
     if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-        countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
-
-        if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum - 1] +
-                              robustHistogram->bounds->data.F32[LQBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum] +
-                              robustHistogram->bounds->data.F32[LQBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum + 1] +
-                              robustHistogram->bounds->data.F32[LQBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum];
-            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1];
-
-            if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n");
-            }
-
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0);
-
-        } else {
-            // If the LQ is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] +
-                                     robustHistogram->bounds->data.F32[LQBinNum]);
-        }
-
-        if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum - 1] +
-                              robustHistogram->bounds->data.F32[UQBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum] +
-                              robustHistogram->bounds->data.F32[UQBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum + 1] +
-                              robustHistogram->bounds->data.F32[UQBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum];
-            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1];
-
-            if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) &&
-                    ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n");
-            }
-
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            stats->robustUQ = p_ps1DPolyMedian(myPoly,
-                                               x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0);
-        } else {
-            // If the UQ is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] +
-                                     robustHistogram->bounds->data.F32[UQBinNum]);
-        }
-    }
+        countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];
+
+        stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                          cumulativeRobustSums,
+                          LQBinNum,
+                          countFloat/4.0);
+        stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                          cumulativeRobustSums,
+                          UQBinNum,
+                          3.0 * countFloat/4.0);
+    }
+    // XXX: I think sumNfit == sumN50 here.
     stats->robustNfit = sumNfit;
     stats->robustN50 = sumN50;
 
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
     psFree(tmpStats);
     psFree(robustHistogram);
-    psFree(myPoly);
-    psFree(cumulativeRobustSumsFullRange);
-    psFree(cumulativeRobustSumsDlRange);
+    psFree(robustHistogramVector);
+    psFree(cumulativeRobustSums);
 }
 
@@ -1468 +1350 @@
     // I am commenting this code out, and replacing it with code which
     // calculates the mean directly on the robustHistogram.
-    //
-    // XXX: Replace this with the LaGrange interpolation stuff.
  
     domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
@@ -1479 +1359 @@
     max_pos = -1;
     for (i=0;i<robustHistogramVector->n;i++) {
-        domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] +
-                                        robustHistogram->bounds->data.F32[i]);
+        domain->data.F32[i][0] =  PS_BIN_MIDPOINT(robustHistogram, i);bounds->data.F32[modeBinNum + 1]);
         data->data.F32[i] = (float) robustHistogramVector->data.F32[i];
         errors->data.F32[i] = 1.0;
@@ -1695 +1574 @@
     psScalar tmpScalar;
     tmpScalar.type.type = PS_TYPE_F32;
+    psVector* inF32;
+    int mustFreeTmp = 1;
 
     // NOTE: Verify that this is the correct action.
@@ -1706 +1587 @@
 
     if (out->nums->type.type != PS_TYPE_U32) {
-        psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member.");
+        psAbort(__func__, "Only data type PS_TYPE_U32 for out->nums member.");
     }
     // NOTE: Verify that this is the correct action.
@@ -1713 +1594 @@
     }
 
-    if (in->type.type != PS_TYPE_F32) {
-        psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type);
-    }
-
     if (mask != NULL) {
         if (in->n != mask->n) {
-            psAbort(__func__, "Vector data and vector mask are of different sizes.");
+            psError(__func__, "Vector data and vector mask are of different sizes.");
         }
         if (mask->type.type != PS_TYPE_U8) {
-            psAbort(__func__, "Vector mask must be type PS_TYPE_U8");
-        }
+            psError(__func__, "Vector mask must be type PS_TYPE_U8");
+        }
+    }
+
+    inF32 = p_psConvertToF32((psVector *) in);
+    if (inF32 == NULL) {
+        inF32 = (psVector *) in;
+        mustFreeTmp = 0;
     }
 
@@ -1730 +1613 @@
 
     numBins = out->nums->n;
-    for (i = 0; i < in->n; i++) {
+    for (i = 0; i < inF32->n; i++) {
         // Check if this pixel is masked, and if so, skip it.
         if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
             // Check if this pixel is below the minimum value, and if so
             // count it, then skip it.
-            if (in->data.F32[i] < out->bounds->data.F32[0]) {
+            if (inF32->data.F32[i] < out->bounds->data.F32[0]) {
                 out->minNum++;
                 // Check if this pixel is above the maximum value, and if so
                 // count it, then skip it.
-            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
+            } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) {
                 out->maxNum++;
             } else {
@@ -1746 +1629 @@
                 if (out->uniform == true) {
                     binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
-                    binNum = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize);
+                    binNum = (int)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
 
                     // NOTE: This next if-statement really shouldn't be necessary.
@@ -1755 +1638 @@
                     }
 
-                    (out->nums->data.S32[binNum])++;
+                    (out->nums->data.U32[binNum])++;
 
                     // If this is a non-uniform histogram, determining the correct
                     // bin number requires a bit more work.
                 } else {
-                    tmpScalar.data.F32 = in->data.F32[i];
-                    tmp = p_psVectorBinDisect(out->bounds,
-                                              &tmpScalar);
-                    if (tmp == -1) {
-                        //                      XXX: Generate warning message
+                    tmpScalar.data.F32 = inF32->data.F32[i];
+                    tmp = p_psVectorBinDisect(out->bounds, &tmpScalar);
+                    if (tmp < 0) {
+                        // XXX: Generate warning message
                     } else {
-                        (out->nums->data.S32[tmp])++;
+                        (out->nums->data.U32[tmp])++;
                     }
                 }
@@ -1772 +1654 @@
         }
     }
+
+    if (mustFreeTmp == 1) {
+        psFree(inF32);
+    }
     return (out);
 }
 
 /******************************************************************************
-p_psConvertToF32(myVector, maskVector, maskVal, stats): this is the cheap
-way to support a variety of vector data types: we simply convert the input
-vector to F32 at the beginning, and write all of our functions in F32.  If
-the vast majority of all vector stat operations are F32 (or any other single
-type), then this is probably the best way to go.  Otherwise, when the
-algorithms stablize, we will then macro everything and put type support in
-the various stat functions.
+p_psConvertToF32(in): this is the cheap way to support a variety of vector
+data types: we simply convert the input vector to F32 at the beginning, and
+write all of our functions in F32.  If the vast majority of all vector stat
+operations are F32 (or any other single type), then this is probably the
+best way to go.  Otherwise, when the algorithms stablize, we will then macro
+everything and put type support in the various stat functions.
  
 XXX: Should the default data type be F64?  Since we are buying Athlons...
  *****************************************************************************/
-psVector* p_psConvertToF32(psStats* stats, psVector* in, psVector* mask, unsigned int maskVal)
+psVector* p_psConvertToF32(psVector* in)
 {
     int i = 0;
@@ -1800 +1685 @@
         for (i = 0; i < in->n; i++) {
             tmp->data.F32[i] = (float)in->data.U16[i];
+        }
+    } else if (in->type.type == PS_TYPE_U8) {
+        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
+        for (i = 0; i < in->n; i++) {
+            tmp->data.F32[i] = (float)in->data.U8[i];
         }
     } else if (in->type.type == PS_TYPE_F64) {
@@ -1839 +1729 @@
     int mustFreeTmp = 1;
 
-    // NOTE: Verify that this is the correct action.
     if (in == NULL) {
-        return (stats);
-    }
-    if (stats == NULL) {
-        return (NULL);
-    }
-
-    inF32 = p_psConvertToF32(stats, in, mask, maskVal);
+        psError(__func__, "in is NULL\n");
+        return(stats);
+    }
+    PS_CHECK_NULL_PTR_RETURN_NULL(stats);
+
+    inF32 = p_psConvertToF32((psVector *) in);
     if (inF32 == NULL) {
         inF32 = in;
         mustFreeTmp = 0;
     }
-    // XXX: Should we abort if (stats->min == stats->max)?
     if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
-        psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
+        psError(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
+        return(stats);
     }
 

trunk/psLib/src/math/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-14 23:43:12 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -107 +107 @@
 } \
 
+/** Preprocessor macro to generate error on a NULL vector */
+#define PS_CHECK_NULL_VECTOR_RETURN_F0(NAME) \
+if (NAME == NULL || NAME->data.V == NULL) { \
+    psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \
+    return(0.0); \
+} \
+
 /** Preprocessor macro to generate error on a NULL poniter */
 #define PS_CHECK_NULL_PTR_RETURN_NULL(NAME) \
@@ -113 +120 @@
     return(NULL); \
 }
+
+/** Preprocessor macro to generate error on a NULL poniter */
+#define PS_CHECK_NULL_PTR_RETURN_F0(NAME) \
+if (NAME == NULL) { \
+    psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \
+    return(0.0); \
+}
+
 /** Preprocessor macro to generate error on a NULL poniter */
 #define PS_CHECK_NULL_PTR_RETURN_ZERO(NAME) \
@@ -127 +142 @@
 } \
 
+/** Preprocessor macro to generate error for zero length vector */
+#define PS_CHECK_EMPTY_VECTOR_RETURN_F0(NAME) \
+if (NAME->n < 1) { \
+    psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \
+    return(0.0); \
+} \
+
 /** Preprocessor macro to generate error on differing size vectors */
 #define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(VEC1, VEC2) \
@@ -132 +154 @@
     psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     return(NULL); \
+}
+
+/** Preprocessor macro to generate error on differing size vectors */
+#define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_F0(VEC1, VEC2) \
+if (VEC1->n != VEC2->n) { \
+    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
+    return(0.0); \
 }
 
@@ -199 +228 @@
 
 #define PS_MAX(A, B) \
-((A > B) ? A : B) \
+(((A) > (B)) ? (A) : (B)) \
 
 #define PS_MIN(A, B) \
-((A < B) ? A : B) \
-
+(((A) < (B)) ? (A) : (B)) \
+

trunk/psLib/src/math/psMinimize.c

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

trunk/psLib/src/math/psMinimize.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-06 23:38:49 $
+ *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:43 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -84 +84 @@
 
 
+/*
+   Bug 203:
+ 
+typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params,
+const psVector *x);
+ 
+This function takes the current guess for the parameters for which we are
+trying to get the best values (params), and a single vector (x) of
+conditions from the array of x values fed into the minimiser.
+ 
+    bool psMinimizeLMChi2(psMinimization *min,
+                          psImage *covar,
+                          psVector *params,
+                          const psVector *paramMask,
+                          const psArray *x,
+                          const psVector *y,
+                          const psVector *yErr,
+                          psMinimizeLMChi2Func func);
+ 
+This takes the minimization specs (min), returns the covariance matrix
+(covar), takes the best guess of initial parameters (params), the parameter
+mask (paramMask), and takes multiple vectors of conditions in an array (x),
+the corresponding measured values (y) and errors (yErr) and the function to
+fit (func).
+ 
+For example, for GRB afterglows, I have flux as a function of time and
+frequency.  So I stuff into the "psArray *x" all my time and frequency values
+(so a whole heap of vectors of size 2), I have the measured values in "y" and
+errors in "yErr".  Then each of the time-frequency pairs are passed to my model
+function with the current parameters, and the model function returns the flux
+for that time-frequency pair, and the derivative with respect to each of the
+parameters.
+ 
+This seems reasonable.  The only thing we could change would be to have the
+function be defined:
+ 
+typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,
+                                          const psVector *params,
+                                          const psArray *x);
+ 
+So it would return the model value for each of the measurements at once, and
+return for each the derivatives (so that it returns a matrix).
+ 
+What do you think?
+ 
+ 
+ 
+I'm not sure I understand how LM chi-squared minimization will work with the
+following function that you define:
+ 
+    typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv,
+                                              const psVector *params,
+                                              const psArray *x);
+ 
+The chi-squared minimization algorithm, as defined in NR, requires that
+function to be minimized be evaluated at each data point, and that all
+derivatives, with respect to each parameter, be calculated at each data
+point.  In the above, can I assume that
+ 
+    x is an array of psVectors, with each vector corresponding to a single
+ data point.
+ 
+    The returned value has the same length as x.  It contains the value of
+ the function at each data point in x.
+ 
+    deriv: an n-by-p matrix where "n" is the number of data points, and "p"
+ is the number of parameters.  The [i][j] element of this matrix
+ holds the derivative of the function at the i-th data point with
+ respect to the j-th parameter.
+ 
+ 
+ 
+*/
+
+
 // XXX: What if any of these arguments are NULL?
+
+
 typedef
 float (*psMinimizeLMChi2Func)(psVector *deriv,
@@ -98 +175 @@
                       psMinimizeLMChi2Func func);
 
+
+typedef
+psVector* (*psMyMinimizeLMChi2Func)(psImage *deriv,
+                                    const psVector *params,
+                                    const psArray *x);
+
+bool psMyMinimizeLMChi2(psMinimization *min,
+                        psImage *covar,
+                        psVector *params,
+                        const psVector *paramMask,
+                        const psArray *x,
+                        const psVector *y,
+                        const psVector *yErr,
+                        psMyMinimizeLMChi2Func func);
+
 typedef
 float (*psMinimizePowellFunc)(const psVector *params,

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-14 01:44:48 $
+ *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-14 01:44:48 $
+ *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.67 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-13 02:40:13 $
+ *  @version $Revision: 1.68 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-26 21:24:43 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -40 +40 @@
 /*****************************************************************************/
 // will use robust statistical methods.
-#define GAUSS_WIDTH 20                     // The width of the Gaussian or boxcar smoothing.
+#define GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
 #define CLIPPED_NUM_ITER_LB 1
 #define CLIPPED_NUM_ITER_UB 10
@@ -57 +57 @@
 /* TYPE DEFINITIONS                                                          */
 /*****************************************************************************/
+psVector* p_psConvertToF32(psVector* in);
 
 /*****************************************************************************/
@@ -423 +424 @@
     float rangeMax = 0.0;       // Exclude date above this
 
-    // Determine if the number of data points exceed a threshold which will
-    // cause to generate robust stats, as opposed to exact stats.
-    // XXX: This code is no longer used.  We now calculate the median exactly
-    // regardless of the vector size.
-    /*
-     * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
-     * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
-     * 
-     * // Store the robust quartiles into the sample quartile members. stats->sampleMedian =
-     * stats2->robustMedian;
-     * 
-     * // Free temporary data buffers. psFree(stats2);
-     * 
-     * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options |
-     * PS_STAT_ROBUST_FOR_SAMPLE;
-     * 
-     * return; } */
-
     // Determine how many data points fit inside this min/max range
     // and are not masked, IF the maskVector is not NULL>
@@ -507 +490 @@
 
 /******************************************************************************
-    This routine smoothes the data in the input robustHistogram with a
-    Gaussian of width sigma.
- 
-    XXX: Must consult with IfA on the proper values for the Gaussian
-    smoothing.  Currently, the Gaussian coefficients are such that only
-    the center coefficient is non-zero.  This is because "e" is being
-    raised to a very large negative number for all points except the
-    center point.
- 
-    XXX: use a static variable for gaussianCoefs[] and compute them once.
- *****************************************************************************/
-psVector* p_psVectorsmoothHistGaussian(psHistogram* robustHistogram,
+p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
+robustHistogram with a Gaussian of width sigma.
+ *****************************************************************************/
+#define GAUSS_WIDTH_OTHER 5.0
+psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram,
                                        float sigma)
 {
     int i = 0;                  // Loop index variable
     int j = 0;                  // Loop index variable
-    float denom = 0.0;          // Temporary variable
-    float expo = 0.0;           // Temporary variable
-    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
-    psVector* smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32);
-
-    for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) {
-        if (fabs(sigma) >= FLT_EPSILON) {
-            // If sigma does not equal zero, then we use Gaussian smoothing.
-            #ifdef  DARWIN
-            denom = (float)sqrt(2.0 * M_PI * sigma * sigma);
-            #else
-
-            denom = sqrtf(2.0 * M_PI * sigma * sigma);
-            #endif
-
-            expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH));
-            expo /= (2.0 * sigma * sigma);
-            gaussianCoefs[i] = exp(expo / denom);
-
-            //printf("p_psVectorsmoothHistGaussian(): gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
-            // NOTE: Gaussian smoothing just isn't working with low sigma
-            // values.  The problem is that the Gaussian coefficients are
-            // all zero, except for the middle coefficient, which is exactly
-            // one.  Therefore, I'm using boxcar smoothing.
-            gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH));
-            // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
-        } else {
-            /* If sigma equals zero (all pixels have the same value) the above code will divide by zero.
-             * Therefore, we don't need to smooth the data. */
-            for (i = 0; i < robustHistogram->nums->n; i++) {
-                smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i];
-            }
-            return (smooth);
-        }
-    }
-
-    // Perform the actual smoothing.
-    for (i = 0; i < robustHistogram->nums->n; i++) {
+    float iMid;
+    float jMid;
+    int numBins = robustHistogram->nums->n;
+    int numBounds = robustHistogram->bounds->n;
+    psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32);
+    int jMin = 0;
+    int jMax = 0;
+    float firstBound = robustHistogram->bounds->data.F32[0];
+    float lastBound = robustHistogram->bounds->data.F32[numBounds-1];
+    psScalar x;
+
+    x.type.type = PS_TYPE_F32;
+    for (i = 0; i < numBins; i++) {
+        iMid = (robustHistogram->bounds->data.F32[i] +
+                robustHistogram->bounds->data.F32[i+1]) / 2.0;
+
+        // YYY: The p_psVectorBinDisect() routine does much of the work of
+        // the following conditionals, however, it also reports a warning
+        // message.  I don't want the warning message so I reproduce the
+        // conditionals here.  Maybe p_psVectorBinDisect() should not produce
+        // warnings?
+
+        x.data.F32 = iMid - (GAUSS_WIDTH * sigma);
+        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+            jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
+        } else if (x.data.F32 <= firstBound) {
+            jMin = 0;
+        } else if (x.data.F32 >= lastBound) {
+            jMin = robustHistogram->bounds->n - 1;
+
+        }
+
+        x.data.F32 = iMid + (GAUSS_WIDTH * sigma);
+        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
+            jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
+        } else if (x.data.F32 <= firstBound) {
+            jMax = 0;
+        } else if (x.data.F32 >= lastBound) {
+            jMax = robustHistogram->bounds->n - 1;
+
+        }
+
         smooth->data.F32[i] = 0.0;
-        for (j = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) {
-            if (((j + i) >= 0) && ((j + i) < smooth->n)) {
-                smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] *
-                                        (float)robustHistogram->nums->data.S32[j + i]);
-            }
-        }
-    }
-    return (smooth);
+        for (j = jMin ; j <= jMax ; j++) {
+            jMid = (robustHistogram->bounds->data.F32[j] +
+                    robustHistogram->bounds->data.F32[j+1]) / 2.0;
+            smooth->data.F32[i] +=
+                ((float) robustHistogram->nums->data.U32[j]) *
+                psGaussian(jMid, iMid, sigma, true);
+        }
+    }
+
+    return(smooth);
 }
 
@@ -913 +895 @@
 specified data point.
  *****************************************************************************/
-float p_psGaussian(const psVector* restrict myData, const psVector* restrict myParams)
+float p_psGaussian(const psVector* restrict myData,
+                   const psVector* restrict myParams)
 {
     float x = myData->data.F32[0];
@@ -930 +913 @@
 calculates the specified partial derivative of the above Gaussian function.
  *****************************************************************************/
-float p_psGaussianDeriv(const psVector* restrict myData, const psVector* restrict myParams, int whichParam)
+float p_psGaussianDeriv(const psVector* restrict myData,
+                        const psVector* restrict myParams, int whichParam)
 {
     float x = myData->data.F32[0];
@@ -1036 +1020 @@
 
 /******************************************************************************
-p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given
-the specified histogram, with the specified pre-calculated cumulativeSums,
-this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum",
-and then finds the point x such that f(x) == fitFloat.
+fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): This routine
+takes psVectors of x/y pairs as input, and fits a quadratic to the 3 points
+surrounding element binNum in the vectors (the midpoint between element i
+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: This function is currently not being used.
- *****************************************************************************/
-float p_psFitQuadratic(psHistogram* histogram,
-                       psVector* cumulativeSums,
-                       int binNum,
-                       float fitFloat)
-{
+// XXX: Use static variables.
+*****************************************************************************/
+float fitQuadraticSearchForYThenReturnX(psVector *xVec,
+                                        psVector *yVec,
+                                        int binNum,
+                                        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);
-
-    if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) {
-        x->data.F64[0] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]);
-        x->data.F64[1] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]);
-        x->data.F64[2] = (double)0.5 *
-                         (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]);
-
-        y->data.F64[0] = cumulativeSums->data.F32[binNum - 1];
-        y->data.F64[1] = cumulativeSums->data.F32[binNum];
-        y->data.F64[2] = cumulativeSums->data.F32[binNum + 1];
-
-        if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) {
-            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
-        }
-
+    float tmpFloat;
+
+    if ((binNum > 0) && (binNum < (yVec->n - 2))) {
+        x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
+        x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
+        x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
+        y->data.F64[0] = yVec->data.F32[binNum - 1];
+        y->data.F64[1] = yVec->data.F32[binNum];
+        y->data.F64[2] = yVec->data.F32[binNum + 1];
+
+        // Ensure that yVal is within the range of the bins we are using.
+        if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) {
+            printf("((%f), %f, %f)\n", yVal, y->data.F64[0], y->data.F64[2]);
+            psError(__func__, "yVal not within y-range\n");
+        }
         yErr->data.F64[0] = 1.0;
         yErr->data.F64[1] = 1.0;
         yErr->data.F64[2] = 1.0;
 
+        // Determine the coefficients of the polynomial.
         myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-        return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat));
+        // Call p_ps1DPolyMedian(), which does a binary search on the
+        // polynomial, looking for the value x such that f(x) = yVal
+        tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);
     } else {
-        return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum]));
+        if (binNum == 0) {
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        } else if (binNum == (xVec->n - 1)) {
+            // XXX: Is this right?
+            tmpFloat = xVec->data.F32[binNum];
+        } else if (binNum == (xVec->n - 2)) {
+            // XXX: Is this right?
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
+                              xVec->data.F32[binNum + 1]);
+        }
     }
 
@@ -1083 +1097 @@
     psFree(yErr);
     psFree(myPoly);
-    return (0.0);
-}
+    return(tmpFloat);
+}
+
+#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
+(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
 
 /******************************************************************************
@@ -1111 +1128 @@
 *****************************************************************************/
 void p_psVectorRobustStats(const psVector* restrict myVector,
-                           const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
+                           const psVector* restrict maskVector,
+                           unsigned int maskVal,
+                           psStats* stats)
 {
     psHistogram* robustHistogram = NULL;
@@ -1118 +1137 @@
     int LQBinNum = -1;          // Bin num for lower quartile
     int UQBinNum = -1;          // Bin num for upper quartile
+    int medianBinNum = -1;
     int i = 0;                  // Loop index variable
-    int maxBinNum = 0;
-    float maxBinCount = 0.0;
+    int modeBinNum = 0;
+    float modeBinCount = 0.0;
     float dL = 0.0;
     int numBins = 0;
-    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);
-
-    // psImage* domain;
-    // psVector* errors;
-    // psVector* data;
-    // psVector* initialGuess;
-    // psVector* theParams;
-    // float chiSq=0.0;
-    // float max = -HUGE;
-    // float max_pos;
     float myMean = 0.0;
     float myStdev = 0.0;
@@ -1139 +1149 @@
     float sumSquares = 0.0;
     float sumDiffs = 0.0;
-    psVector* x = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* y = psVectorAlloc(3, PS_TYPE_F64);
-    psVector* yErr = psVectorAlloc(3, PS_TYPE_F64);
-    psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
-    psVector* cumulativeRobustSumsFullRange = NULL;
-    psVector* cumulativeRobustSumsDlRange = NULL;
+    psVector* cumulativeRobustSums = NULL;
     float sumRobust = 0.0;
     float sumN50 = 0.0;
     float sumNfit = 0.0;
-    float cumulativeMedian = 0.0;
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F32;
+    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);
 
     // Compute the initial bin size of the robust histogram.  This is done
@@ -1188 +1195 @@
         p_psVectorMax(myVector, maskVector, maskVal, stats);
     }
+
     // Create the histogram structure.  NOTE: we can not specify the bin size
     // precisely since the argument to psHistogramAlloc() is the number of
@@ -1198 +1206 @@
     psVectorHistogram(robustHistogram, myVector, maskVector, maskVal);
 
-    // Smooth the histogram.
-    // XXX: is that the right stdev?
-    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram,
+    // Smooth the histogram, Gaussian-style.
+    robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram,
                             tmpStats->clippedStdev / 4.0f);
 
@@ -1208 +1215 @@
 
     /**************************************************************************
-    Determine the lower/upper quartiles.
+    Determine the median/lower/upper quartile bin numbers.
+
+    We define a vector called "cumulativeRobustSums" where the value at
+    index position i is equal to the sum of bins 0:i.  This will be used in
+    determining the median and lower/upper quartiles.
     **************************************************************************/
-    // We define a vector called "cumulativeRobustSums..." where the value at
-    // index position i is equal to the sum of bins 0:i.  This will be used
-    // now and later in determining the lower/upper quartiles.
-    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0];
+    cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
+    cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0];
     for (i = 1; i < robustHistogramVector->n; i++) {
-        cumulativeRobustSumsFullRange->data.F32[i] =
-            cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
-    }
-    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
-
-    // Determine the bin number containing the lower quartile point.
-    LQBinNum = -1;
-    for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) {
-        if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) {
-            LQBinNum = i;
-            break;
-        }
-    }
-
-    // Determine the bin number containing the upper quartile point.
-    UQBinNum = -1;
-    for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) {
-        if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) {
-            UQBinNum = i;
-            break;
-        }
-    }
-
-    if ((LQBinNum == -1) || (UQBinNum == -1)) {
-        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
-    }
-
+        cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] +
+                                            robustHistogramVector->data.F32[i];
+    }
+    sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];
+
+    tmpScalar.data.F32 = sumRobust / 4.0;
+    LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+    tmpScalar.data.F32 = 3.0 * sumRobust / 4.0;
+    UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+    tmpScalar.data.F32 = sumRobust / 2.0;
+    medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
+
+    if ((LQBinNum < 0) || (UQBinNum < 0)) {
+        psError(__func__, "Could not determine the robust lower/upper quartile bin numbers.");
+    }
+    if (medianBinNum < 0) {
+        psError(__func__, "Could not determine the robust lower/upper quartile bin numbers.");
+    }
     /**************************************************************************
     Determine the mode in the range LQ:UQ.
     **************************************************************************/
     // Determine the bin with the peak value in the range LQ to UQ.
-    maxBinNum = LQBinNum;
-    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
-    sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum];
+    modeBinNum = LQBinNum;
+    modeBinCount = robustHistogramVector->data.F32[LQBinNum];
+    sumN50 = (float)robustHistogram->nums->data.U32[LQBinNum];
     for (i = LQBinNum + 1; i <= UQBinNum; i++) {
-        if (robustHistogramVector->data.F32[i] > maxBinCount) {
-            maxBinNum = i;
-            maxBinCount = robustHistogramVector->data.F32[i];
-        }
-        sumN50 += (float)robustHistogram->nums->data.S32[i];
-    }
-
-    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
+        if (robustHistogramVector->data.F32[i] > modeBinCount) {
+            modeBinNum = i;
+            modeBinCount = robustHistogramVector->data.F32[i];
+        }
+        sumN50 += (float)robustHistogram->nums->data.U32[i];
+    }
+
     dL = (UQBinNum - LQBinNum) / 4;
-
-    printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum);
-
     /**************************************************************************
     Determine the mean/stdev for the bins in the range mode-dL to mode+dL
     **************************************************************************/
-    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
-    for (i = 0; i < robustHistogramVector->n; i++) {
-        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
-    }
-    sumNfit = 0.0;
-    cumulativeMedian = 0.0;
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
-        if ((0 <= i) && (i < robustHistogramVector->n)) {
-            cumulativeRobustSumsDlRange->data.F32[i] =
-                cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
-            cumulativeMedian += robustHistogramVector->data.F32[i];
-            sumNfit += (float)robustHistogram->nums->data.S32[i];
-        }
-    }
-
     // Calculate the mean of the smoothed robust histogram in the range
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     // that bin (this is a non-exact approximation).
+    sumNfit = 0.0;
     myMean = 0.0;
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
+    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
         if ((0 <= i) && (i < robustHistogramVector->n)) {
-            myMean += (robustHistogramVector->data.F32[i]) * 0.5 *
-                      (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]);
+            myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i);
             countFloat += robustHistogramVector->data.F32[i];
         }
+
+        sumNfit += (float)robustHistogram->nums->data.U32[i];
     }
     myMean /= countFloat;
@@ -1297 +1279 @@
     // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
     // that bin.
-    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
+    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
         if ((0 <= i) && (i < robustHistogramVector->n)) {
-            diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] +
-                           robustHistogram->bounds->data.F32[i])) - myMean;
+            diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean;
             sumSquares += diff * diff * robustHistogramVector->data.F32[i];
             sumDiffs += diff * robustHistogramVector->data.F32[i];
@@ -1315 +1296 @@
 
     if (stats->options & PS_STAT_ROBUST_MODE) {
-        stats->robustMode = 0.5 *
-                            (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]);
+        stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum);
     }
 
@@ -1322 +1302 @@
         stats->robustStdev = myStdev;
     }
-    // To determine the median (and later, the lower/upper quartile), we fit
-    // a quadratic to the three bins surrounding the bin containing the median.
-    // The quadratic y=f(x) with x being the midpoint of each bin, and y being
-    // the cumulative number of data points in all bins up to, and including,
-    // this bin.  We then solve the quadratic for
-
-    // XXX: Create a "p_psQuadraticFitAndSolveForX() function.
-
+
+    // To determine the median, we fit a quadratic y=f(x) to the three bins
+    // surrounding the bin containing the median (x is the midpoint of each
+    // bin and y is the value of each bin).  Then we figure out what value
+    // of x corresponds to f(x) being the median (half of all points).
     if (stats->options & PS_STAT_ROBUST_MEDIAN) {
-        if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum - 1] +
-                              robustHistogram->bounds->data.F32[maxBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum] +
-                              robustHistogram->bounds->data.F32[maxBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[maxBinNum + 1] +
-                              robustHistogram->bounds->data.F32[maxBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum];
-            y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1];
-
-            // Ensure that cumulativeMedian/2 is actually within the range of the bins
-            // we are using.
-            cumulativeMedian *= 0.5;
-            if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) {
-                printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]);
-                psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n");
-            }
-            // XXX: yErr is not currently used by psVectorFitPolynomial1D().  We
-            // may have to set this meaningfully later.
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            // Determine the coefficients of the polynomial.
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            // Call p_ps1DPolyMedian(), which does a binary search on the
-            // polynomial, looking for the value x such that
-            // f(x) = cumulativeMedian.
-            stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian);
-        } else {
-            // If the mode is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] +
-                                         robustHistogram->bounds->data.F32[maxBinNum]);
-        }
-    }
-    // The lower/upper quartile calculations are very similar to the median
-    // calculations.  We fit a quadratic to the array containing the
-    // cumulative data points in each bin, then search for x such that
-    // f(x) equals the lower/upper quartile exactly.
-    //
+        // Take a psVector.  Fit a polynomial y = f(x) to the 3 data elements
+        // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0.
+        stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                              cumulativeRobustSums,
+                              medianBinNum,
+                              sumRobust/2.0);
+    }
+
+    // To determine the quartiles, we fit a quadratic y=f(x) to the three bins
+    // surrounding the bin containing LQ/UQ (x is the midpoint of each
+    // bin and y is the value of each bin).  Then we figure out what value
+    // of x corresponds to f(x) being the LQ/UQ.
     if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-        countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];
-
-        if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum - 1] +
-                              robustHistogram->bounds->data.F32[LQBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum] +
-                              robustHistogram->bounds->data.F32[LQBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[LQBinNum + 1] +
-                              robustHistogram->bounds->data.F32[LQBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum];
-            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1];
-
-            if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n");
-            }
-
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0);
-
-        } else {
-            // If the LQ is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] +
-                                     robustHistogram->bounds->data.F32[LQBinNum]);
-        }
-
-        if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) {
-            x->data.F64[0] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum - 1] +
-                              robustHistogram->bounds->data.F32[UQBinNum]);
-            x->data.F64[1] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum] +
-                              robustHistogram->bounds->data.F32[UQBinNum + 1]);
-            x->data.F64[2] = (double)0.5 *
-                             (robustHistogram->bounds->data.F32[UQBinNum + 1] +
-                              robustHistogram->bounds->data.F32[UQBinNum + 2]);
-
-            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1];
-            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum];
-            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1];
-
-            if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) &&
-                    ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) {
-                psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n");
-            }
-
-            yErr->data.F64[0] = 1.0;
-            yErr->data.F64[1] = 1.0;
-            yErr->data.F64[2] = 1.0;
-
-            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
-            stats->robustUQ = p_ps1DPolyMedian(myPoly,
-                                               x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0);
-        } else {
-            // If the UQ is the first/last histogram bin, then simply use
-            // the midpoint of that bin.
-            stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] +
-                                     robustHistogram->bounds->data.F32[UQBinNum]);
-        }
-    }
+        countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];
+
+        stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                          cumulativeRobustSums,
+                          LQBinNum,
+                          countFloat/4.0);
+        stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds,
+                          cumulativeRobustSums,
+                          UQBinNum,
+                          3.0 * countFloat/4.0);
+    }
+    // XXX: I think sumNfit == sumN50 here.
     stats->robustNfit = sumNfit;
     stats->robustN50 = sumN50;
 
-    psFree(x);
-    psFree(y);
-    psFree(yErr);
     psFree(tmpStats);
     psFree(robustHistogram);
-    psFree(myPoly);
-    psFree(cumulativeRobustSumsFullRange);
-    psFree(cumulativeRobustSumsDlRange);
+    psFree(robustHistogramVector);
+    psFree(cumulativeRobustSums);
 }
 
@@ -1468 +1350 @@
     // I am commenting this code out, and replacing it with code which
     // calculates the mean directly on the robustHistogram.
-    //
-    // XXX: Replace this with the LaGrange interpolation stuff.
  
     domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
@@ -1479 +1359 @@
     max_pos = -1;
     for (i=0;i<robustHistogramVector->n;i++) {
-        domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] +
-                                        robustHistogram->bounds->data.F32[i]);
+        domain->data.F32[i][0] =  PS_BIN_MIDPOINT(robustHistogram, i);bounds->data.F32[modeBinNum + 1]);
         data->data.F32[i] = (float) robustHistogramVector->data.F32[i];
         errors->data.F32[i] = 1.0;
@@ -1695 +1574 @@
     psScalar tmpScalar;
     tmpScalar.type.type = PS_TYPE_F32;
+    psVector* inF32;
+    int mustFreeTmp = 1;
 
     // NOTE: Verify that this is the correct action.
@@ -1706 +1587 @@
 
     if (out->nums->type.type != PS_TYPE_U32) {
-        psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member.");
+        psAbort(__func__, "Only data type PS_TYPE_U32 for out->nums member.");
     }
     // NOTE: Verify that this is the correct action.
@@ -1713 +1594 @@
     }
 
-    if (in->type.type != PS_TYPE_F32) {
-        psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type);
-    }
-
     if (mask != NULL) {
         if (in->n != mask->n) {
-            psAbort(__func__, "Vector data and vector mask are of different sizes.");
+            psError(__func__, "Vector data and vector mask are of different sizes.");
         }
         if (mask->type.type != PS_TYPE_U8) {
-            psAbort(__func__, "Vector mask must be type PS_TYPE_U8");
-        }
+            psError(__func__, "Vector mask must be type PS_TYPE_U8");
+        }
+    }
+
+    inF32 = p_psConvertToF32((psVector *) in);
+    if (inF32 == NULL) {
+        inF32 = (psVector *) in;
+        mustFreeTmp = 0;
     }
 
@@ -1730 +1613 @@
 
     numBins = out->nums->n;
-    for (i = 0; i < in->n; i++) {
+    for (i = 0; i < inF32->n; i++) {
         // Check if this pixel is masked, and if so, skip it.
         if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
             // Check if this pixel is below the minimum value, and if so
             // count it, then skip it.
-            if (in->data.F32[i] < out->bounds->data.F32[0]) {
+            if (inF32->data.F32[i] < out->bounds->data.F32[0]) {
                 out->minNum++;
                 // Check if this pixel is above the maximum value, and if so
                 // count it, then skip it.
-            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
+            } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) {
                 out->maxNum++;
             } else {
@@ -1746 +1629 @@
                 if (out->uniform == true) {
                     binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
-                    binNum = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize);
+                    binNum = (int)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
 
                     // NOTE: This next if-statement really shouldn't be necessary.
@@ -1755 +1638 @@
                     }
 
-                    (out->nums->data.S32[binNum])++;
+                    (out->nums->data.U32[binNum])++;
 
                     // If this is a non-uniform histogram, determining the correct
                     // bin number requires a bit more work.
                 } else {
-                    tmpScalar.data.F32 = in->data.F32[i];
-                    tmp = p_psVectorBinDisect(out->bounds,
-                                              &tmpScalar);
-                    if (tmp == -1) {
-                        //                      XXX: Generate warning message
+                    tmpScalar.data.F32 = inF32->data.F32[i];
+                    tmp = p_psVectorBinDisect(out->bounds, &tmpScalar);
+                    if (tmp < 0) {
+                        // XXX: Generate warning message
                     } else {
-                        (out->nums->data.S32[tmp])++;
+                        (out->nums->data.U32[tmp])++;
                     }
                 }
@@ -1772 +1654 @@
         }
     }
+
+    if (mustFreeTmp == 1) {
+        psFree(inF32);
+    }
     return (out);
 }
 
 /******************************************************************************
-p_psConvertToF32(myVector, maskVector, maskVal, stats): this is the cheap
-way to support a variety of vector data types: we simply convert the input
-vector to F32 at the beginning, and write all of our functions in F32.  If
-the vast majority of all vector stat operations are F32 (or any other single
-type), then this is probably the best way to go.  Otherwise, when the
-algorithms stablize, we will then macro everything and put type support in
-the various stat functions.
+p_psConvertToF32(in): this is the cheap way to support a variety of vector
+data types: we simply convert the input vector to F32 at the beginning, and
+write all of our functions in F32.  If the vast majority of all vector stat
+operations are F32 (or any other single type), then this is probably the
+best way to go.  Otherwise, when the algorithms stablize, we will then macro
+everything and put type support in the various stat functions.
  
 XXX: Should the default data type be F64?  Since we are buying Athlons...
  *****************************************************************************/
-psVector* p_psConvertToF32(psStats* stats, psVector* in, psVector* mask, unsigned int maskVal)
+psVector* p_psConvertToF32(psVector* in)
 {
     int i = 0;
@@ -1800 +1685 @@
         for (i = 0; i < in->n; i++) {
             tmp->data.F32[i] = (float)in->data.U16[i];
+        }
+    } else if (in->type.type == PS_TYPE_U8) {
+        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
+        for (i = 0; i < in->n; i++) {
+            tmp->data.F32[i] = (float)in->data.U8[i];
         }
     } else if (in->type.type == PS_TYPE_F64) {
@@ -1839 +1729 @@
     int mustFreeTmp = 1;
 
-    // NOTE: Verify that this is the correct action.
     if (in == NULL) {
-        return (stats);
-    }
-    if (stats == NULL) {
-        return (NULL);
-    }
-
-    inF32 = p_psConvertToF32(stats, in, mask, maskVal);
+        psError(__func__, "in is NULL\n");
+        return(stats);
+    }
+    PS_CHECK_NULL_PTR_RETURN_NULL(stats);
+
+    inF32 = p_psConvertToF32((psVector *) in);
     if (inF32 == NULL) {
         inF32 = in;
         mustFreeTmp = 0;
     }
-    // XXX: Should we abort if (stats->min == stats->max)?
     if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
-        psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
+        psError(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
+        return(stats);
     }
 

trunk/psLib/test/dataManip/tst_psMinimize06.c

@@ -74 +74 @@
 }
 
-
-int main()
+psVector *myFunc2(psImage *myDeriv,
+                  psVector *myParams,
+                  psArray *myCoords)
+{
+    psVector *sum = psVectorAlloc(myCoords->n, PS_TYPE_F32);
+    int i;
+    int j;
+
+
+    if (myDeriv == NULL) {
+        myDeriv = psImageAlloc(myParams->n, myCoords->n, PS_TYPE_F32);
+        psError(__func__, "myDeriv is NULL.\n");
+    }
+
+    for (i=0;i<N;i++) {
+        sum->data.F32[i] = MIN_VALUE;
+        for (j=0;j<NUM_PARAMS;j++) {
+            sum->data.F32[i]+= (myParams->data.F32[j] - expectedParm[j]) *
+                               (myParams->data.F32[j] - expectedParm[j]);
+
+            myDeriv->data.F32[i][j] = (2.0 * myParams->data.F32[j]) -
+                                      (2.0 * expectedParm[j]);
+        }
+    }
+
+    return(sum);
+}
+
+int t01()
+{
+    int currentId = psMemGetId();
+    int memLeaks = 0;
+    int i = 0;
+    psArray *myCoords;
+    psVector *myParams;
+    psImage *myCovar;
+    psMinimization *min;
+    psVector *y;
+
+    psTraceSetLevel(".psLib", 0);
+    /**************************************************************************
+     *************************************************************************/
+    min = psMinimizationAlloc(NUM_ITERATIONS, ERR_TOL);
+    myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32);
+    myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32);
+    myCoords = psArrayAlloc(N);
+    y = psVectorAlloc(N, PS_TYPE_F32);
+
+    for (i=0;i<N;i++) {
+        myCoords->data[i] = (psPTR *) psVectorAlloc(2, PS_TYPE_F32);
+        ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) (i+10);
+        ((psVector *) (myCoords->data[i]))->data.F32[1] = (float) (i+3);
+        y->data.F32[i] = (float) i;
+    }
+    for (i=0;i<NUM_PARAMS;i++) {
+        expectedParm[i] = 2.42 + (float) (2 * i);
+        myParams->data.F32[i] = (float) i;
+        myParams->data.F32[i] = expectedParm[i] * 1.3;
+        myParams->data.F32[i] = (float) (5 + i);
+        myParams->data.F32[i] = 0.0;
+    }
+
+    psMyMinimizeLMChi2(min,
+                       myCovar,
+                       myParams,
+                       NULL,
+                       myCoords,
+                       y,
+                       NULL,
+                       (psMyMinimizeLMChi2Func) myFunc2);
+
+    printf("\nThe chi-squared is %f\n", min->value);
+    for (i=0;i<NUM_PARAMS;i++) {
+        printf("Parameter %d at the minimum is %f (expected: %f)\n", i,
+               myParams->data.F32[i], expectedParm[i]);
+    }
+
+    psFree(min);
+    psFree(myCovar);
+    psFree(myParams);
+    psFree(myCoords);
+    psFree(y);
+
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,NULL);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+    return (!testStatus);
+}
+
+int t02()
 {
     int currentId = psMemGetId();
@@ -85 +175 @@
     psImage *myCovar;
     psMinimization *min;
-    psVector *x = psVectorAlloc(N, PS_TYPE_F32);
     psVector *y = psVectorAlloc(N, PS_TYPE_F32);
 
     psTraceSetLevel(".psLib", 0);
+    t02();
     /**************************************************************************
      *************************************************************************/
@@ -101 +191 @@
         ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) (i+10);
         ((psVector *) (myCoords->data[i]))->data.F32[1] = (float) (i+3);
-        x->data.F32[i] = (float) i;
         y->data.F32[i] = (float) i;
     }
@@ -131 +220 @@
     psFree(myParamMask);
     psFree(min);
-    psFree(x);
     psFree(y);
     psFree(myCovar);
@@ -142 +230 @@
     return (!testStatus);
 }
+
+int main()
+{
+    t01();
+    //    t02();
+}

trunk/psLib/test/dataManip/tst_psStats07.c

@@ -10 +10 @@
 #include <math.h>
 
-#define N 200
+#define N 90
 #define MEAN 32.0
 #define STDEV 2.0
 #define ERROR_TOLERANCE 0.10
 
-int main()
+int t00()
 {
     psStats * myStats = NULL;
@@ -32 +32 @@
     float realMedianNoMask = MEAN;
     float realModeNoMask = MEAN;
-    float realStdevNoMask = STDEV * 0.33;
+    float realStdevNoMask = STDEV * 0.20;
     float realLQNoMask = MEAN - ( 0.6 * STDEV );
     float realUQNoMask = MEAN + ( 0.6 * STDEV );
     int realN50NoMask = N / 4;
     int realNfitNoMask = N / 4;
-    float realMeanWithMask = MEAN;
-    float realMedianWithMask = MEAN;
-    float realModeWithMask = MEAN;
-    float realStdevWithMask = STDEV;
-    float realLQWithMask = MEAN;
-    float realUQWithMask = MEAN;
-    int realN50WithMask = N / 4;
-    int realNfitWithMask = N / 4;
-
+
+    psTraceSetLevel(".psLib.dataManip.psStats", 10);
 
     /*************************************************************************/
@@ -232 +225 @@
                  "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask",
                  testStatus );
-
-
-    return ( 0 );
+    return(testStatus);
+
+    /*************************************************************************/
+    /*  Deallocate data structures                                           */
+    /*************************************************************************/
+    printPositiveTestHeader( stdout,
+                             "psStats functions",
+                             "psStats(): deallocating memory" );
+
+    psFree( myStats );
+    psFree( myVector );
+    psFree( maskVector );
+
+    psMemCheckCorruption( 1 );
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if ( 0 != memLeaks ) {
+        psAbort( __func__, "Memory Leaks! (%d leaks)", memLeaks );
+    }
+
+    printFooter( stdout,
+                 "psVector functions",
+                 "psStats(): deallocating memory",
+                 testStatus );
+
+    return ( !globalTestStatus );
+}
+
+
+int t01()
+{
+    psStats * myStats = NULL;
+    int testStatus = true;
+    int globalTestStatus = true;
+    int i = 0;
+    psVector *myVector = NULL;
+    psVector *maskVector = NULL;
+    // NOTE: These values were calculated by running the function on the data.
+    // A: They must be changed if we adjust the number of data points.
+    // B: We don't really know that they are correct.
+    int count = 0;
+    int currentId = psMemGetId();
+    int memLeaks = 0;
+    float realMeanWithMask = MEAN;
+    float realMedianWithMask = MEAN;
+    float realModeWithMask = MEAN;
+    float realStdevWithMask = STDEV * 0.20;
+    float realLQWithMask = MEAN;
+    float realUQWithMask = MEAN;
+    int realN50WithMask = N / 4;
+    int realNfitWithMask = N / 4;
+
+    psTraceSetLevel(".psLib.dataManip.psStats", 10);
+
+    /*************************************************************************/
+    /*  Allocate and initialize data structures                              */
+    /*************************************************************************/
+    myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN |
+                            PS_STAT_ROBUST_MEDIAN |
+                            PS_STAT_ROBUST_MODE |
+                            PS_STAT_ROBUST_STDEV |
+                            PS_STAT_ROBUST_QUARTILE );
+
+    maskVector = psVectorAlloc( N, PS_TYPE_U8 );
+    maskVector->n = N;
+    myVector = psGaussianDev( MEAN, STDEV, N );
+    // Set the mask vector and calculate the expected maximum.
+    for ( i = 0;i < N;i++ ) {
+        if ( i < ( N / 2 ) ) {
+            maskVector->data.U8[ i ] = 0;
+            count++;
+        } else {
+            maskVector->data.U8[ i ] = 1;
+        }
+    }
 
     /*************************************************************************/
@@ -422 +486 @@
     return ( !globalTestStatus );
 }
+
+int main()
+{
+    t00();
+    t01();
+}