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

Timestamp:

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

Author:

gusciora

Message:

Modified the psVectorfit() tests and code.

File:

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

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: