Changeset 1945

Timestamp:

Oct 3, 2004, 1:35:47 PM (22 years ago)

Author:

gusciora

Message:

Modified the psVectorfit() tests and code.

Location:

trunk/psLib

Files:

• src/dataManip/psConstants.h (modified) (4 diffs)
• src/dataManip/psFunctions.c (modified) (6 diffs)
• src/dataManip/psMinimize.c (modified) (30 diffs)
• src/math/psConstants.h (modified) (4 diffs)
• src/math/psMinimize.c (modified) (30 diffs)
• src/math/psPolynomial.c (modified) (6 diffs)
• src/math/psSpline.c (modified) (6 diffs)
• test/dataManip/Makefile (modified) (2 diffs)
• test/dataManip/tst_psMinimize04.c (modified) (6 diffs)
• test/dataManip/tst_psMinimize04b.c (modified) (6 diffs)

trunk/psLib/src/dataManip/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-28 23:27:37 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-03 23:35:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -30 +30 @@
 }
 
+/** Preprocessor macro to generate error on a NULL vector */
+#define PS_CHECK_NULL_VECTOR_ACTION(NAME, ACTION) \
+if (NAME == NULL || NAME->data.V == NULL) { \
+    if (ACTION == 0) { \
+        psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
+        return(NULL); \
+    } \
+    if (ACTION == 1) { \
+        psAbort(__func__, "Invalid operation: %s or its data is NULL.", #NAME); \
+    } \
+}
+
 /** Preprocessor macro to generate error on a NULL poniter */
 #define PS_CHECK_NULL_PTR(NAME) \
@@ -40 +52 @@
 if (NAME->n < 1) { \
     psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
+}
+
+/** Preprocessor macro to generate error for zero length vector */
+#define PS_CHECK_EMPTY_VECTOR_ACTION(NAME, ACTION) \
+if (NAME->n < 1) { \
+    if (ACTION == 0) { \
+        psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
+        return(NULL); \
+    } \
+    if (ACTION == 1) { \
+        psAbort(__func__,"Invalid operation: %s has zero n value.", #NAME); \
+    } \
 }
 
@@ -67 +91 @@
 }
 
+/** Preprocessor macro to generate error on a NULL 1DPolynomial */
+#define PS_CHECK_NULL_1DPOLY_ACTION(NAME, ACTION)                                                          \
+if (NAME == NULL || NAME->coeff == NULL) {                                                         \
+    if (ACTION == 0) { \
+        psError(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
+        return(NULL); \
+    } \
+    if (ACTION == 1) { \
+        psAbort(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
+    } \
+}
+
 #define PS_PRINT_VECTOR(NAME) \
 for (int my_i=0;my_i<NAME->n;my_i++) { \

trunk/psLib/src/dataManip/psFunctions.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-28 23:27:37 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-03 23:35:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -486 +486 @@
 float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
 {
-    /*
     psVector *d;
     int n;
@@ -508 +507 @@
     psFree(d);
     return(tmp);
-    */
+    /*
 
     int n;
@@ -527 +526 @@
 
     return(tmp);
+    */
 }
 
@@ -1869 +1869 @@
     p_psMemSetPersistent(p, true);
     p_psMemSetPersistent(p->data.F32, true);
+    /*
+        psVector *p = psVectorAlloc(n, PS_TYPE_F32);
+        float tmp;
+    */
 
     psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
@@ -1901 +1905 @@
             "---- p_ps1DFullInterpolateF32() end ----\n");
 
+    /*
+        tmp = p->data.F32[0];
+        psFree(p);
+        return(tmp);
+    */
     return(p->data.F32[0]);
 }

trunk/psLib/src/dataManip/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-28 23:27:37 $
+ *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-03 23:35:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -682 +682 @@
     psScalar *fScalar;
 
-
     // XXX: These assignments appear too simple to warrant code and
     // variable declarations.  I retain them here to maintain coherence
@@ -714 +713 @@
         fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32,
                                         3, tmpScalar);
+
         f->data.F64[i] = (double) fScalar->data.F32;
+        psFree(fScalar);
 
         psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 6,
@@ -724 +725 @@
     // coefficients of the Chebyshev polynomial: NR 5.8.7.
     fac = 2.0/((float) n);
-    for (j=0;j<n;j++) {
+
+    // XXX: is this loop bound correct?
+    for (j=0;j<myPoly->n;j++) {
         sum = 0.0;
         for (k=0;k<n;k++) {
@@ -730 +733 @@
                   cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n));
         }
+
         myPoly->coeff[j] = fac * sum;
     }
 
     // XXX: Must free memory.
+    psFree(f);
     psFree(x32);
     psFree(y32);
     psFree(tmpScalar);
-    psFree(fScalar);
+
     return(myPoly);
 }
@@ -745 +750 @@
 polynomial of degree myPoly to the data points (x, y) and return the
 coefficients of that polynomial.
- 
-XXX: yErr is currently ignored.
  
 XXX: must add type F32 (currently F64 only).
@@ -795 +798 @@
 
     // Initialize data structures.
-    for (i = 0; i < (polyOrder); i++) {
+    for (i = 0; i < polyOrder; i++) {
         B->data.F64[i] = 0.0;
         coeffs->data.F64[i] = 0.0;
         outPerm->data.F64[i] = 0.0;
-        for (j = 0; j < (polyOrder); j++) {
+        for (j = 0; j < polyOrder; j++) {
             A->data.F64[i][j] = 0.0;
             ALUD->data.F64[i][j] = 0.0;
@@ -809 +812 @@
 
     // Build the B and A data structs.
-    for (i = 0; i < X->n; i++) {
-        p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
-
-        for (k = 0; k < (polyOrder); k++) {
-            B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k];
-        }
-
-        for (k = 0; k < (polyOrder); k++) {
-            for (j = 0; j < (polyOrder); j++) {
-                A->data.F64[k][j] += xSums->data.F64[k + j];
+    if (yErr == NULL) {
+        for (i = 0; i < X->n; i++) {
+            p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
+
+            for (k = 0; k < polyOrder; k++) {
+                B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k];
+            }
+
+            for (k = 0; k < polyOrder; k++) {
+                for (j = 0; j < polyOrder; j++) {
+                    A->data.F64[k][j] += xSums->data.F64[k + j];
+                }
+            }
+        }
+    } else {
+        for (i = 0; i < X->n; i++) {
+            p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
+
+            for (k = 0; k < polyOrder; k++) {
+                B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k] /
+                                  yErr->data.F64[i];
+            }
+
+            for (k = 0; k < polyOrder; k++) {
+                for (j = 0; j < polyOrder; j++) {
+                    A->data.F64[k][j] += xSums->data.F64[k + j] /
+                                         yErr->data.F64[i];
+                }
             }
         }
@@ -826 +847 @@
     coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
 
-    for (k = 0; k < (polyOrder); k++) {
+    for (k = 0; k < polyOrder; k++) {
         myPoly->coeff[k] = coeffs->data.F64[k];
         // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]);
@@ -853 +874 @@
 myPoly to the data points (x, y) and return the coefficients of that
 polynomial.
- 
-XXX: yErr is currently ignored.
  
 XXX: must add type F32 (currently F64 only).
@@ -870 +889 @@
     psVector *myYErr = NULL;
 
-    PS_CHECK_NULL_1DPOLY(myPoly);
-    PS_CHECK_NULL_VECTOR(y);
-    PS_CHECK_EMPTY_VECTOR(y);
+    PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
+    PS_CHECK_NULL_VECTOR_ACTION(y, 0);
+    PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
 
     // If yErr==NULL, set all errors equal.
@@ -880 +899 @@
 
         if (y->type.type == PS_TYPE_F32) {
-            for (i=0;i<yErr->n;i++) {
+            for (i=0;i<myYErr->n;i++) {
                 myYErr->data.F32[i] = 1.0;
             }
         } else if (y->type.type == PS_TYPE_F64) {
-            for (i=0;i<yErr->n;i++) {
+            for (i=0;i<myYErr->n;i++) {
                 myYErr->data.F64[i] = 1.0;
             }
@@ -893 +912 @@
 
     // If x==NULL, create an myX vector with x values set to (0:n), and if
-    // this is a CHebyshev polynomial, we must scale to (-1:1).
+    // this is a Chebyshev polynomial, we must scale to (-1:1).
 
     // XXX: Verify that this is the correct action.
@@ -902 +921 @@
         if (y->type.type == PS_TYPE_F32) {
             if (myPoly->type == PS_POLYNOMIAL_ORD) {
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F32[i] = (float) i;
                 }
@@ -909 +928 @@
                 float max = (float) (y->n - 1);
 
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) /
                                        (0.5 * (max - min));
@@ -916 +935 @@
         } else if (y->type.type == PS_TYPE_F64) {
             if (myPoly->type == PS_POLYNOMIAL_ORD) {
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F64[i] = (float) i;
                 }
@@ -923 +942 @@
                 double max = (double) (y->n - 1);
 
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) /
                                        (0.5 * (max - min));
@@ -1161 +1180 @@
 
 /******************************************************************************
-    This routine must minimize a possibly multi-dimensional function
-    along a vector defined by line.
- 
-XXX: Use a p_psName().
+This routine takes as input a possibly multi-dimensional function, along
+with an initial guess at the parameters of that function and vector "line"
+of the same size as the parameter vector.  It will minimize the function
+along that vector anr returns the offset along that vector at which the
+minimum is determined.
+ 
+XXX: This routine is not very efficient in terms of total evaluations of the
+function.
  *****************************************************************************/
-float psLineMin(psMinimization *min,
-                psVector *params,
-                psVector *line,
-                const psVector *paramMask,
-                const psArray *coords,
-                psMinimizePowellFunc func)
+float p_psLineMin(psMinimization *min,
+                  psVector *params,
+                  psVector *line,
+                  const psVector *paramMask,
+                  const psArray *coords,
+                  psMinimizePowellFunc func)
 {
     psVector *bracket;
@@ -1190 +1213 @@
     int null = 0;
 
-    psTrace(".psLib.dataManip.psLineMin", 4,
-            "---- psLineMin() begin ----\n");
-    psTrace(".psLib.dataManip.psLineMin", 6,
+    psTrace(".psLib.dataManip.p_psLineMin", 4,
+            "---- p_psLineMin() begin ----\n");
+    psTrace(".psLib.dataManip.p_psLineMin", 6,
             "min->maxIter is %d\n", min->maxIter);
-    psTrace(".psLib.dataManip.psLineMin", 6,
+    psTrace(".psLib.dataManip.p_psLineMin", 6,
             "min->tol is %f\n", min->tol);
 
@@ -1202 +1225 @@
                 if (line->data.F32[i] >= FLT_EPSILON) {
                     null = 1;
-                    psTrace(".psLib.dataManip.psLineMin", 4,
+                    psTrace(".psLib.dataManip.p_psLineMin", 4,
                             "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
                 }
@@ -1211 +1234 @@
             if (line->data.F32[i] >= FLT_EPSILON) {
                 null = 1;
-                psTrace(".psLib.dataManip.psLineMin", 4,
+                psTrace(".psLib.dataManip.p_psLineMin", 4,
                         "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
             }
@@ -1218 +1241 @@
 
     if (null == 0) {
-        psTrace(".psLib.dataManip.psLineMin", 2,
-                "psLineMin() called with zero line vector.\n");
+        psTrace(".psLib.dataManip.p_psLineMin", 2,
+                "p_psLineMin() called with zero line vector.\n");
         return(0.0);
     }
@@ -1229 +1252 @@
 
     for (i=0;i<params->n;i++) {
-        psTrace(".psLib.dataManip.psLineMin", 6,
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
                 "params, paramMask, line [%d] is (%f %d %f)\n", i,
                 params->data.F32[i],
@@ -1248 +1271 @@
     while (min->iter < min->maxIter) {
         min->iter++;
-        psTrace(".psLib.dataManip.psLineMin", 6,
-                "psLineMin(): iteration %d\n", min->iter);
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
+                "p_psLineMin(): iteration %d\n", min->iter);
 
         a = bracket->data.F32[0];
@@ -1269 +1292 @@
         fb = func(tmpb, coords);
         fc = func(tmpc, coords);
-        psTrace(".psLib.dataManip.psLineMin", 6,
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
                 "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter,
                 a, b, c, fa, fb, fc);
@@ -1307 +1330 @@
             }
         }
-        psTrace(".psLib.dataManip.psLineMin", 6,
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
                 "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
 
         mul = bracket->data.F32[1];
-        if ((fabs(a-b) < min->tol) &&
-                (fabs(b-c) < min->tol)) {
+        if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) {
             for (i=0;i<params->n;i++) {
                 if (paramMask->data.U8[i] == 0) {
@@ -1324 +1346 @@
             psFree(tmpc);
             psFree(tmpn);
-            psTrace(".psLib.dataManip.psLineMin", 4,
-                    "---- psLineMin() end.a (%f) ----\n", mul);
+            psTrace(".psLib.dataManip.p_psLineMin", 4,
+                    "---- p_psLineMin() end.a (%f) ----\n", mul);
             return(mul);
         }
@@ -1336 +1358 @@
     psFree(tmpn);
 
-    psTrace(".psLib.dataManip.psLineMin", 4,
-            "---- psLineMin() end.b (0.0) ----\n");
+    psTrace(".psLib.dataManip.p_psLineMin", 4,
+            "---- p_psLineMin() end.b (0.0) ----\n");
     return(0.0);
 }
@@ -1429 +1451 @@
                 dummyMin.maxIter = 100;
                 dummyMin.tol = 0.01;
-                mul = psLineMin(&dummyMin, Q, v[i], paramMask, coords, func);
+                mul = p_psLineMin(&dummyMin, Q, v[i], paramMask, coords, func);
                 if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
                     biggestDiff = fabs(dummyMin.value - currFuncVal);
@@ -1460 +1482 @@
         }
 
-        mul = psLineMin(min, params, u, paramMask, coords, func);
+        mul = p_psLineMin(min, params, u, paramMask, coords, func);
 
         // 6:

trunk/psLib/src/math/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-28 23:27:37 $
+ *  @version $Revision: 1.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-03 23:35:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -30 +30 @@
 }
 
+/** Preprocessor macro to generate error on a NULL vector */
+#define PS_CHECK_NULL_VECTOR_ACTION(NAME, ACTION) \
+if (NAME == NULL || NAME->data.V == NULL) { \
+    if (ACTION == 0) { \
+        psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \
+        return(NULL); \
+    } \
+    if (ACTION == 1) { \
+        psAbort(__func__, "Invalid operation: %s or its data is NULL.", #NAME); \
+    } \
+}
+
 /** Preprocessor macro to generate error on a NULL poniter */
 #define PS_CHECK_NULL_PTR(NAME) \
@@ -40 +52 @@
 if (NAME->n < 1) { \
     psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
+}
+
+/** Preprocessor macro to generate error for zero length vector */
+#define PS_CHECK_EMPTY_VECTOR_ACTION(NAME, ACTION) \
+if (NAME->n < 1) { \
+    if (ACTION == 0) { \
+        psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \
+        return(NULL); \
+    } \
+    if (ACTION == 1) { \
+        psAbort(__func__,"Invalid operation: %s has zero n value.", #NAME); \
+    } \
 }
 
@@ -67 +91 @@
 }
 
+/** Preprocessor macro to generate error on a NULL 1DPolynomial */
+#define PS_CHECK_NULL_1DPOLY_ACTION(NAME, ACTION)                                                          \
+if (NAME == NULL || NAME->coeff == NULL) {                                                         \
+    if (ACTION == 0) { \
+        psError(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
+        return(NULL); \
+    } \
+    if (ACTION == 1) { \
+        psAbort(__func__,"Invalid operation: %s or its coeffs is NULL.", #NAME);                          \
+    } \
+}
+
 #define PS_PRINT_VECTOR(NAME) \
 for (int my_i=0;my_i<NAME->n;my_i++) { \

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-28 23:27:37 $
+ *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-03 23:35:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -682 +682 @@
     psScalar *fScalar;
 
-
     // XXX: These assignments appear too simple to warrant code and
     // variable declarations.  I retain them here to maintain coherence
@@ -714 +713 @@
         fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32,
                                         3, tmpScalar);
+
         f->data.F64[i] = (double) fScalar->data.F32;
+        psFree(fScalar);
 
         psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 6,
@@ -724 +725 @@
     // coefficients of the Chebyshev polynomial: NR 5.8.7.
     fac = 2.0/((float) n);
-    for (j=0;j<n;j++) {
+
+    // XXX: is this loop bound correct?
+    for (j=0;j<myPoly->n;j++) {
         sum = 0.0;
         for (k=0;k<n;k++) {
@@ -730 +733 @@
                   cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n));
         }
+
         myPoly->coeff[j] = fac * sum;
     }
 
     // XXX: Must free memory.
+    psFree(f);
     psFree(x32);
     psFree(y32);
     psFree(tmpScalar);
-    psFree(fScalar);
+
     return(myPoly);
 }
@@ -745 +750 @@
 polynomial of degree myPoly to the data points (x, y) and return the
 coefficients of that polynomial.
- 
-XXX: yErr is currently ignored.
  
 XXX: must add type F32 (currently F64 only).
@@ -795 +798 @@
 
     // Initialize data structures.
-    for (i = 0; i < (polyOrder); i++) {
+    for (i = 0; i < polyOrder; i++) {
         B->data.F64[i] = 0.0;
         coeffs->data.F64[i] = 0.0;
         outPerm->data.F64[i] = 0.0;
-        for (j = 0; j < (polyOrder); j++) {
+        for (j = 0; j < polyOrder; j++) {
             A->data.F64[i][j] = 0.0;
             ALUD->data.F64[i][j] = 0.0;
@@ -809 +812 @@
 
     // Build the B and A data structs.
-    for (i = 0; i < X->n; i++) {
-        p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
-
-        for (k = 0; k < (polyOrder); k++) {
-            B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k];
-        }
-
-        for (k = 0; k < (polyOrder); k++) {
-            for (j = 0; j < (polyOrder); j++) {
-                A->data.F64[k][j] += xSums->data.F64[k + j];
+    if (yErr == NULL) {
+        for (i = 0; i < X->n; i++) {
+            p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
+
+            for (k = 0; k < polyOrder; k++) {
+                B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k];
+            }
+
+            for (k = 0; k < polyOrder; k++) {
+                for (j = 0; j < polyOrder; j++) {
+                    A->data.F64[k][j] += xSums->data.F64[k + j];
+                }
+            }
+        }
+    } else {
+        for (i = 0; i < X->n; i++) {
+            p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums);
+
+            for (k = 0; k < polyOrder; k++) {
+                B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k] /
+                                  yErr->data.F64[i];
+            }
+
+            for (k = 0; k < polyOrder; k++) {
+                for (j = 0; j < polyOrder; j++) {
+                    A->data.F64[k][j] += xSums->data.F64[k + j] /
+                                         yErr->data.F64[i];
+                }
             }
         }
@@ -826 +847 @@
     coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
 
-    for (k = 0; k < (polyOrder); k++) {
+    for (k = 0; k < polyOrder; k++) {
         myPoly->coeff[k] = coeffs->data.F64[k];
         // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]);
@@ -853 +874 @@
 myPoly to the data points (x, y) and return the coefficients of that
 polynomial.
- 
-XXX: yErr is currently ignored.
  
 XXX: must add type F32 (currently F64 only).
@@ -870 +889 @@
     psVector *myYErr = NULL;
 
-    PS_CHECK_NULL_1DPOLY(myPoly);
-    PS_CHECK_NULL_VECTOR(y);
-    PS_CHECK_EMPTY_VECTOR(y);
+    PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
+    PS_CHECK_NULL_VECTOR_ACTION(y, 0);
+    PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
 
     // If yErr==NULL, set all errors equal.
@@ -880 +899 @@
 
         if (y->type.type == PS_TYPE_F32) {
-            for (i=0;i<yErr->n;i++) {
+            for (i=0;i<myYErr->n;i++) {
                 myYErr->data.F32[i] = 1.0;
             }
         } else if (y->type.type == PS_TYPE_F64) {
-            for (i=0;i<yErr->n;i++) {
+            for (i=0;i<myYErr->n;i++) {
                 myYErr->data.F64[i] = 1.0;
             }
@@ -893 +912 @@
 
     // If x==NULL, create an myX vector with x values set to (0:n), and if
-    // this is a CHebyshev polynomial, we must scale to (-1:1).
+    // this is a Chebyshev polynomial, we must scale to (-1:1).
 
     // XXX: Verify that this is the correct action.
@@ -902 +921 @@
         if (y->type.type == PS_TYPE_F32) {
             if (myPoly->type == PS_POLYNOMIAL_ORD) {
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F32[i] = (float) i;
                 }
@@ -909 +928 @@
                 float max = (float) (y->n - 1);
 
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) /
                                        (0.5 * (max - min));
@@ -916 +935 @@
         } else if (y->type.type == PS_TYPE_F64) {
             if (myPoly->type == PS_POLYNOMIAL_ORD) {
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F64[i] = (float) i;
                 }
@@ -923 +942 @@
                 double max = (double) (y->n - 1);
 
-                for (i=0;i<yErr->n;i++) {
+                for (i=0;i<myX->n;i++) {
                     myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) /
                                        (0.5 * (max - min));
@@ -1161 +1180 @@
 
 /******************************************************************************
-    This routine must minimize a possibly multi-dimensional function
-    along a vector defined by line.
- 
-XXX: Use a p_psName().
+This routine takes as input a possibly multi-dimensional function, along
+with an initial guess at the parameters of that function and vector "line"
+of the same size as the parameter vector.  It will minimize the function
+along that vector anr returns the offset along that vector at which the
+minimum is determined.
+ 
+XXX: This routine is not very efficient in terms of total evaluations of the
+function.
  *****************************************************************************/
-float psLineMin(psMinimization *min,
-                psVector *params,
-                psVector *line,
-                const psVector *paramMask,
-                const psArray *coords,
-                psMinimizePowellFunc func)
+float p_psLineMin(psMinimization *min,
+                  psVector *params,
+                  psVector *line,
+                  const psVector *paramMask,
+                  const psArray *coords,
+                  psMinimizePowellFunc func)
 {
     psVector *bracket;
@@ -1190 +1213 @@
     int null = 0;
 
-    psTrace(".psLib.dataManip.psLineMin", 4,
-            "---- psLineMin() begin ----\n");
-    psTrace(".psLib.dataManip.psLineMin", 6,
+    psTrace(".psLib.dataManip.p_psLineMin", 4,
+            "---- p_psLineMin() begin ----\n");
+    psTrace(".psLib.dataManip.p_psLineMin", 6,
             "min->maxIter is %d\n", min->maxIter);
-    psTrace(".psLib.dataManip.psLineMin", 6,
+    psTrace(".psLib.dataManip.p_psLineMin", 6,
             "min->tol is %f\n", min->tol);
 
@@ -1202 +1225 @@
                 if (line->data.F32[i] >= FLT_EPSILON) {
                     null = 1;
-                    psTrace(".psLib.dataManip.psLineMin", 4,
+                    psTrace(".psLib.dataManip.p_psLineMin", 4,
                             "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
                 }
@@ -1211 +1234 @@
             if (line->data.F32[i] >= FLT_EPSILON) {
                 null = 1;
-                psTrace(".psLib.dataManip.psLineMin", 4,
+                psTrace(".psLib.dataManip.p_psLineMin", 4,
                         "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
             }
@@ -1218 +1241 @@
 
     if (null == 0) {
-        psTrace(".psLib.dataManip.psLineMin", 2,
-                "psLineMin() called with zero line vector.\n");
+        psTrace(".psLib.dataManip.p_psLineMin", 2,
+                "p_psLineMin() called with zero line vector.\n");
         return(0.0);
     }
@@ -1229 +1252 @@
 
     for (i=0;i<params->n;i++) {
-        psTrace(".psLib.dataManip.psLineMin", 6,
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
                 "params, paramMask, line [%d] is (%f %d %f)\n", i,
                 params->data.F32[i],
@@ -1248 +1271 @@
     while (min->iter < min->maxIter) {
         min->iter++;
-        psTrace(".psLib.dataManip.psLineMin", 6,
-                "psLineMin(): iteration %d\n", min->iter);
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
+                "p_psLineMin(): iteration %d\n", min->iter);
 
         a = bracket->data.F32[0];
@@ -1269 +1292 @@
         fb = func(tmpb, coords);
         fc = func(tmpc, coords);
-        psTrace(".psLib.dataManip.psLineMin", 6,
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
                 "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter,
                 a, b, c, fa, fb, fc);
@@ -1307 +1330 @@
             }
         }
-        psTrace(".psLib.dataManip.psLineMin", 6,
+        psTrace(".psLib.dataManip.p_psLineMin", 6,
                 "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
 
         mul = bracket->data.F32[1];
-        if ((fabs(a-b) < min->tol) &&
-                (fabs(b-c) < min->tol)) {
+        if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) {
             for (i=0;i<params->n;i++) {
                 if (paramMask->data.U8[i] == 0) {
@@ -1324 +1346 @@
             psFree(tmpc);
             psFree(tmpn);
-            psTrace(".psLib.dataManip.psLineMin", 4,
-                    "---- psLineMin() end.a (%f) ----\n", mul);
+            psTrace(".psLib.dataManip.p_psLineMin", 4,
+                    "---- p_psLineMin() end.a (%f) ----\n", mul);
             return(mul);
         }
@@ -1336 +1358 @@
     psFree(tmpn);
 
-    psTrace(".psLib.dataManip.psLineMin", 4,
-            "---- psLineMin() end.b (0.0) ----\n");
+    psTrace(".psLib.dataManip.p_psLineMin", 4,
+            "---- p_psLineMin() end.b (0.0) ----\n");
     return(0.0);
 }
@@ -1429 +1451 @@
                 dummyMin.maxIter = 100;
                 dummyMin.tol = 0.01;
-                mul = psLineMin(&dummyMin, Q, v[i], paramMask, coords, func);
+                mul = p_psLineMin(&dummyMin, Q, v[i], paramMask, coords, func);
                 if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
                     biggestDiff = fabs(dummyMin.value - currFuncVal);
@@ -1460 +1482 @@
         }
 
-        mul = psLineMin(min, params, u, paramMask, coords, func);
+        mul = p_psLineMin(min, params, u, paramMask, coords, func);
 
         // 6:

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-28 23:27:37 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-03 23:35:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -486 +486 @@
 float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
 {
-    /*
     psVector *d;
     int n;
@@ -508 +507 @@
     psFree(d);
     return(tmp);
-    */
+    /*
 
     int n;
@@ -527 +526 @@
 
     return(tmp);
+    */
 }
 
@@ -1869 +1869 @@
     p_psMemSetPersistent(p, true);
     p_psMemSetPersistent(p->data.F32, true);
+    /*
+        psVector *p = psVectorAlloc(n, PS_TYPE_F32);
+        float tmp;
+    */
 
     psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
@@ -1901 +1905 @@
             "---- p_ps1DFullInterpolateF32() end ----\n");
 
+    /*
+        tmp = p->data.F32[0];
+        psFree(p);
+        return(tmp);
+    */
     return(p->data.F32[0]);
 }

trunk/psLib/src/math/psSpline.c

@@ -7 +7 @@
  *  polynomials.  It also contains a Gaussian functions.
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-28 23:27:37 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-03 23:35:47 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -486 +486 @@
 float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly)
 {
-    /*
     psVector *d;
     int n;
@@ -508 +507 @@
     psFree(d);
     return(tmp);
-    */
+    /*
 
     int n;
@@ -527 +526 @@
 
     return(tmp);
+    */
 }
 
@@ -1869 +1869 @@
     p_psMemSetPersistent(p, true);
     p_psMemSetPersistent(p->data.F32, true);
+    /*
+        psVector *p = psVectorAlloc(n, PS_TYPE_F32);
+        float tmp;
+    */
 
     psTrace(".psLib.dataManip.psFunctions.p_ps1DFullInterpolateF32", 4,
@@ -1901 +1905 @@
             "---- p_ps1DFullInterpolateF32() end ----\n");
 
+    /*
+        tmp = p->data.F32[0];
+        psFree(p);
+        return(tmp);
+    */
     return(p->data.F32[0]);
 }

trunk/psLib/test/dataManip/Makefile

@@ -3 +3 @@
 ##  Makefile:   test/sysUtils
 ##
-##  $Revision: 1.47 $  $Name: not supported by cvs2svn $
-##  $Date: 2004-09-28 00:30:32 $
+##  $Revision: 1.48 $  $Name: not supported by cvs2svn $
+##  $Date: 2004-10-03 23:35:47 $
 ##
 ##  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -41 +41 @@
 tst_psMatrixVectorArithmetic03 \
 tst_psMinimize04 \
+tst_psMinimize04b \
 tst_psMinimize05 \
 tst_psMinimize06 \

trunk/psLib/test/dataManip/tst_psMinimize04.c

@@ -1 +1 @@
 /*****************************************************************************
-    This routine must ensure that the psVectorFitPolynomial1D works correctly.
+This routine must ensure that the psVectorFitPolynomial1D works correctly.
+We create a vectors of data points (x and y), and populate them with the
+values from an arbitrary function setData().  We then call
+psVectorFitPolynomial1D() with a regular polynomial data structure.  We then
+evaluate the polynomial with the coefficients generated above and determine
+if they are within an error tolerance of the expected values.
+ 
+    t00(): all input vectors are non-NULL.
+    t01(): yErr is NULL.
+    t02(): x, yErr is NULL.
+    t03(): x, y, yErr is NULL.
  *****************************************************************************/
 #include <stdio.h>
@@ -13 +23 @@
 #define NUM_DATA 10
 #define POLY_ORDER 5
-
-double setData(double A,
-               double B,
-               double C,
-               double x)
+#define A 3.0
+#define B 2.0
+#define C 3.0
+#define ERROR_TOLERANCE 0.10
+#define YERR 10.0
+double setData(double x)
 {
     return(A + (B * x) + (C * x * x));
 }
 
-int main()
+int t00()
 {
     psPolynomial1D *myPoly = NULL;
@@ -32 +43 @@
     int testStatus = true;
     int memLeaks = 0;
+    double expectData;
+    double actualData;
 
     myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_ORD);
@@ -40 +53 @@
     for (i=0;i<NUM_DATA;i++) {
         x->data.F64[i] = (double) i;
-        y->data.F64[i] = setData(3.0, 2.0, 3.0, x->data.F64[i]);
-        yErr->data.F64[i] = 0.1;
+        y->data.F64[i] = setData(x->data.F64[i]);
+        yErr->data.F64[i] = YERR;
         printf("Original data %d: (%.1f %.1f)\n", i, x->data.F64[i], y->data.F64[i]);
     }
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): equal errors in yErr");
 
     psVectorFitPolynomial1D(myPoly, x, y, yErr);
@@ -52 +69 @@
 
     for (i=0;i<NUM_DATA;i++) {
-        printf("Fitted data %d: (%.1f %.1f)\n", i, x->data.F64[i],
-               setData(myPoly->coeff[0], myPoly->coeff[1],
-                       myPoly->coeff[2], x->data.F64[i]));
-    }
-
-
-    psMemCheckCorruption(1);
-    printFooter(stdout,
-                "psMinimize functions",
-                "psMinimize(): no masks",
-                testStatus);
+        expectData = setData(x->data.F64[i]);
+        actualData = psPolynomial1DEval(x->data.F64[i], myPoly);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+            testStatus = false;
+        } else {
+            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+        }
+    }
 
     psMemCheckCorruption(1);
@@ -69 +86 @@
     psFree(y);
     psFree(yErr);
-
-    psMemCheckCorruption(1);
-    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
-    if (0 != memLeaks) {
-        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
-    }
-
-    return (!testStatus);
-}
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): equal errors in yErr",
+                testStatus);
+
+    return (!testStatus);
+}
+
+int t01()
+{
+    psPolynomial1D *myPoly = NULL;
+    psVector *x = NULL;
+    psVector *y = NULL;
+    int i = 0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+    double expectData;
+    double actualData;
+
+    myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_ORD);
+    x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+
+    for (i=0;i<NUM_DATA;i++) {
+        x->data.F64[i] = (double) i;
+        y->data.F64[i] = setData(x->data.F64[i]);
+        printf("Original data %d: (%.1f %.1f)\n", i, x->data.F64[i], y->data.F64[i]);
+    }
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): yErr is NULL");
+
+    psVectorFitPolynomial1D(myPoly, x, y, NULL);
+
+    for (i=0;i<POLY_ORDER+1;i++) {
+        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i]);
+    }
+
+    for (i=0;i<NUM_DATA;i++) {
+        expectData = setData(x->data.F64[i]);
+        actualData = psPolynomial1DEval(x->data.F64[i], myPoly);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+            testStatus = false;
+        } else {
+            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+        }
+    }
+
+    psMemCheckCorruption(1);
+    psFree(myPoly);
+    psFree(x);
+    psFree(y);
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): yErr is NULL",
+                testStatus);
+
+    return (!testStatus);
+}
+
+int t02()
+{
+    psPolynomial1D *myPoly = NULL;
+    psVector *y = NULL;
+    int i = 0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+    double expectData;
+    double actualData;
+
+    myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_ORD);
+    y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+
+    for (i=0;i<NUM_DATA;i++) {
+        y->data.F64[i] = setData((double) i);
+        printf("Original data %d: (%.1f)\n", i, y->data.F64[i]);
+    }
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): x, yErr is NULL");
+
+    psVectorFitPolynomial1D(myPoly, NULL, y, NULL);
+
+    for (i=0;i<POLY_ORDER+1;i++) {
+        printf("Polynomial coefficient %d is %0.1f\n", i, myPoly->coeff[i]);
+    }
+
+    for (i=0;i<NUM_DATA;i++) {
+        expectData = setData((double) i);
+        actualData = psPolynomial1DEval((double) i, myPoly);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, (double) i, actualData, expectData);
+            testStatus = false;
+        } else {
+            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, (double) i, actualData, expectData);
+        }
+    }
+
+    psMemCheckCorruption(1);
+    psFree(myPoly);
+    psFree(y);
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): x, yErr is NULL",
+                testStatus);
+
+    return (!testStatus);
+}
+
+int t03()
+{
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): all inputs are NULL");
+
+    psVectorFitPolynomial1D(NULL, NULL, NULL, NULL);
+
+
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): all inputs are NULL",
+                testStatus);
+
+    return (!testStatus);
+}
+
+
+int main()
+{
+    t00();
+    t01();
+    t02();
+    t03();
+}

trunk/psLib/test/dataManip/tst_psMinimize04b.c

@@ -11 +11 @@
 #include "psMinimize.h"
 #include <math.h>
-#define NUM_DATA 15
+#define NUM_DATA 20
 #define POLY_ORDER 10
 #define A 2.0
@@ -18 +18 @@
 #define D 4.0
 #define E 5.0
+#define ERROR_TOLERANCE 0.10
+#define IGNORE (ERROR_TOLERANCE * NUM_DATA)
 
 double setData(double x)
@@ -24 +26 @@
 }
 
-int main()
+int t00()
 {
     psPolynomial1D *myPoly = NULL;
@@ -31 +33 @@
     psVector *yErr = NULL;
     int i = 0;
-    //    int currentId = psMemGetId();
-    int testStatus = true;
-    int memLeaks = 0;
-    float sum=0.0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
 
     myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_CHEB);
@@ -47 +48 @@
     }
     p_psNormalizeVectorF64(x);
-    for (i=0;i<NUM_DATA;i++) {
-        //        printf("Original data %d: (%.1f %.1f)\n", i, x->data.F64[i], y->data.F64[i]);
-    }
-
-    sum=0.0;
-    for (i=0;i<NUM_DATA;i++)
-        sum+=y->data.F64[i];
-    printf("c0 is %f\n", 2.0 * sum/((float) NUM_DATA));
-
-    sum=0.0;
-    for (i=0;i<NUM_DATA;i++)
-        sum+=y->data.F64[i] * x->data.F64[i];
-    printf("c1 is %f\n", 2.0 * sum/((float) NUM_DATA));
-    sum=0.0;
-
-    for (i=0;i<NUM_DATA;i++)
-        sum+=y->data.F64[i] * ((2.0 * x->data.F64[i] * x->data.F64[i]) - 1.0);
-    printf("c2 is %f\n", 2.0 * sum/((float) NUM_DATA));
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): CHEB, equal errors in yErr");
 
     psVectorFitPolynomial1D(myPoly, x, y, yErr);
@@ -72 +59 @@
     }
 
-    for (i=0;i<NUM_DATA;i++) {
-        printf("Fitted data %d: (%.1f %.1f) should be %.1f\n", i, x->data.F64[i],
-               psPolynomial1DEval((float) x->data.F64[i], myPoly), y->data.F64[i]);
-    }
-
-
-    //    psMemCheckCorruption(1);
-    printFooter(stdout,
-                "psMinimize functions",
-                "psMinimize(): no masks",
-                testStatus);
-
-    //    psMemCheckCorruption(1);
+    // We don't test the first or last few data items.
+    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
+        double expectData = y->data.F64[i];
+        double actualData = psPolynomial1DEval(x->data.F64[i], myPoly);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+            testStatus = false;
+        } else {
+            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+        }
+    }
+
+    psMemCheckCorruption(1);
     psFree(myPoly);
     psFree(x);
     psFree(y);
     psFree(yErr);
-
-    //    psMemCheckCorruption(1);
-    //    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
-    if (0 != memLeaks) {
-        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
-    }
-
-    return (!testStatus);
-}
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): CHEB, equal errors in yErr",
+                testStatus);
+    return (!testStatus);
+}
+
+
+int t01()
+{
+    psPolynomial1D *myPoly = NULL;
+    psVector *x = NULL;
+    psVector *y = NULL;
+    int i = 0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+
+    myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_CHEB);
+    x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+
+    for (i=0;i<NUM_DATA;i++) {
+        x->data.F64[i] = (double) i;
+        y->data.F64[i] = setData(x->data.F64[i]);
+    }
+    p_psNormalizeVectorF64(x);
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): CHEB, yErr is NULL");
+
+    psVectorFitPolynomial1D(myPoly, x, y, NULL);
+
+    // We don't test the first or last few data items.
+    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
+        double expectData = y->data.F64[i];
+        double actualData = psPolynomial1DEval(x->data.F64[i], myPoly);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+            testStatus = false;
+        } else {
+            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+        }
+    }
+
+    psMemCheckCorruption(1);
+    psFree(myPoly);
+    psFree(x);
+    psFree(y);
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): CHEB, yErr is NULL",
+                testStatus);
+    return (!testStatus);
+}
+
+
+int t02()
+{
+    psPolynomial1D *myPoly = NULL;
+    psVector *x = NULL;
+    psVector *y = NULL;
+    int i = 0;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+
+    myPoly = psPolynomial1DAlloc(POLY_ORDER+1, PS_POLYNOMIAL_CHEB);
+    x = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+    y = psVectorAlloc(NUM_DATA, PS_TYPE_F64);
+
+    for (i=0;i<NUM_DATA;i++) {
+        x->data.F64[i] = (double) i;
+        y->data.F64[i] = setData(x->data.F64[i]);
+    }
+    p_psNormalizeVectorF64(x);
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): CHEB, x, yErr is NULL");
+
+    psVectorFitPolynomial1D(myPoly, NULL, y, NULL);
+
+    // We don't test the first or last few data items.
+    for (i=IGNORE;i<NUM_DATA-IGNORE;i++) {
+        double expectData = y->data.F64[i];
+        double actualData = psPolynomial1DEval(x->data.F64[i], myPoly);
+        if (fabs(actualData-expectData) > fabs(ERROR_TOLERANCE * expectData)) {
+            printf("ERROR: Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+            testStatus = false;
+        } else {
+            printf("Fitted data %d: (%.1f %.1f), expected was (%.1f)\n",
+                   i, x->data.F64[i], actualData, expectData);
+        }
+    }
+
+    psMemCheckCorruption(1);
+    psFree(myPoly);
+    psFree(x);
+    psFree(y);
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): CHEB, x, yErr is NULL",
+                testStatus);
+    return (!testStatus);
+}
+
+int t03()
+{
+    psPolynomial1D *myPoly = NULL;
+    int currentId = psMemGetId();
+    int testStatus = true;
+    int memLeaks = 0;
+
+    printPositiveTestHeader(stdout,
+                            "psMinimize functions",
+                            "psVectorFitPolynomial1D(): CHEB, yErr is NULL");
+
+    myPoly = psVectorFitPolynomial1D(NULL, NULL, NULL, NULL);
+    if (myPoly != NULL) {
+        printf("ERROR: psVectorFitPolynomial1D() returned a non-NULL polynomial.\n");
+        testStatus = false;
+    }
+
+    psMemCheckCorruption(1);
+    memLeaks = psMemCheckLeaks(currentId,NULL,stderr);
+    if (0 != memLeaks) {
+        psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks);
+    }
+    printFooter(stdout,
+                "psMinimize functions",
+                "psVectorFitPolynomial1D(): CHEB, yErr is NULL",
+                testStatus);
+    return (!testStatus);
+}
+
+int main()
+{
+    t00();
+    t01();
+    t02();
+    t03();
+}