Changeset 1953 for trunk/psLib/src/dataManip/psMinimize.c

Timestamp:

Oct 4, 2004, 3:43:58 PM (22 years ago)

Author:

gusciora

Message:

Added F32 types for the VectorFit routine.

File:

• trunk/psLib/src/dataManip/psMinimize.c (modified) (13 diffs)

trunk/psLib/src/dataManip/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.51 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-10-03 23:35:47 $
+ *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-05 01:43:58 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -24 +24 @@
 /* INCLUDE FILES                                                             */
 /*****************************************************************************/
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <stdarg.h>
-#include <float.h>
-#include <math.h>
-
-#include <gsl/gsl_multifit_nlin.h>
-#include <gsl/gsl_multimin.h>
-#include <gsl/gsl_rng.h>
-#include <gsl/gsl_randist.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_blas.h>
-
-#include "psMemory.h"
-#include "psVector.h"
-#include "psImage.h"
-#include "psTrace.h"
-#include "psLogMsg.h"
-#include "psError.h"
-#include "psAbort.h"
-#include "psFunctions.h"
 #include "psMinimize.h"
-#include "psMatrix.h"
-#include "psConstants.h"
-
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
@@ -655 +630 @@
 }
 
+#define GEN_STATIC_RECYCLED_VECTOR(NAME, SIZE, TYPE) \
+static psVector *NAME = NULL; \
+NAME = psVectorRecycle(NAME, SIZE, TYPE); \
+p_psMemSetPersistent(NAME, true); \
+p_psMemSetPersistent(NAME->data.V, true); \
+
 
 /******************************************************************************
@@ -661 +642 @@
 coefficients of that polynomial.
  
-XXX: yErr is currently ignored.
- 
-XXX: must add type F32 (currently F64 only).
- 
-XXX: Use private name?
- *****************************************************************************/
+XXX: 
+yErr is currently ignored.
+ 
+XXX: 
+Use private name?
+*****************************************************************************/
 psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
         const psVector* restrict x,
@@ -675 +656 @@
     int k;
     int n = x->n;
-    psVector *f = psVectorAlloc(n, PS_TYPE_F64);
     double fac;
     double sum;
-    // XXX: Use static memory here.
-    psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
+    GEN_STATIC_RECYCLED_VECTOR(f, n, PS_TYPE_F64);
     psScalar *fScalar;
+    psScalar tmpScalar;
+    tmpScalar.type.type = PS_TYPE_F64;
 
     // XXX: These assignments appear too simple to warrant code and
@@ -690 +671 @@
     double bpa = 0.5 * (max+min);  // 0
 
-    // XXX: Eliminate this later by generating a F64 version of the
-    // LaGrange interpolation routines.
-    PS_VECTOR_F64_TO_F32(x, x32);
-    PS_VECTOR_F64_TO_F32(y, y32);
-
     // In this loop, we first calculate the values of X for which the
     // Chebyshev polynomials are zero (see NR, section 5.4).  Then we
@@ -707 +683 @@
         double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));
         double X = (Y + bma + bpa) - 1.0;
-        tmpScalar->data.F32 = (float) X;
+        tmpScalar.data.F64 = X;
 
         // We interpolate against are tabluated x,y vectors to determine the
         // function value at X.
-        fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32,
-                                        3, tmpScalar);
-
-        f->data.F64[i] = (double) fScalar->data.F32;
+        fScalar = p_psVectorInterpolate((psVector *) x,
+                                        (psVector *) y,
+                                        3,
+                                        &tmpScalar);
+
+        f->data.F64[i] = fScalar->data.F64;
         psFree(fScalar);
 
@@ -724 +702 @@
     // We have the values for f() at the zero points, we now calculate the
     // coefficients of the Chebyshev polynomial: NR 5.8.7.
+
     fac = 2.0/((float) n);
-
     // XXX: is this loop bound correct?
     for (j=0;j<myPoly->n;j++) {
@@ -737 +715 @@
     }
 
-    // XXX: Must free memory.
-    psFree(f);
-    psFree(x32);
-    psFree(y32);
-    psFree(tmpScalar);
-
     return(myPoly);
 }
@@ -750 +722 @@
 polynomial of degree myPoly to the data points (x, y) and return the
 coefficients of that polynomial.
- 
-XXX: must add type F32 (currently F64 only).
  
 XXX: Use private name?
@@ -875 +845 @@
 polynomial.
  
-XXX: must add type F32 (currently F64 only).
+XXX: type F32 is done via vector conversion only.
  *****************************************************************************/
 psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
@@ -882 +852 @@
                                         const psVector* restrict yErr)
 {
-    bool mustFreeMyYErr = false;
-    bool mustFreemyX = false;
     int i;
     psPolynomial1D *tmpPoly;
-    psVector *myX = NULL;
-    psVector *myYErr = NULL;
-
     PS_CHECK_NULL_1DPOLY_ACTION(myPoly, 0);
     PS_CHECK_NULL_VECTOR_ACTION(y, 0);
     PS_CHECK_EMPTY_VECTOR_ACTION(y, 0);
+    psVector *x64 = NULL;
+    psVector *y64 = NULL;
+    psVector *yErr64 = NULL;
+    static psVector *x64Static = NULL;
+    static psVector *y64Static = NULL;
+    static psVector *yErr64Static = NULL;
+
+    if (y->type.type == PS_TYPE_F64) {
+        y64 = (psVector *) y;
+    } else if (y->type.type == PS_TYPE_F32) {
+        y64Static = psVectorRecycle(y64Static, y->n, PS_TYPE_F64);
+        p_psMemSetPersistent(y64Static, true);
+        p_psMemSetPersistent(y64Static->data.V, true);
+        y64 = y64Static;
+    } else {
+        // XXX: psError() bad type.
+        psAbort(__func__, "Bad type for y64 (%d)", y->type.type);
+    }
+
 
     // If yErr==NULL, set all errors equal.
     if (yErr == NULL) {
-        myYErr = psVectorAlloc(y->n, y->type.type);
-        mustFreeMyYErr = true;
-
-        if (y->type.type == PS_TYPE_F32) {
-            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<myYErr->n;i++) {
-                myYErr->data.F64[i] = 1.0;
-            }
-        }
+        yErr64Static = psVectorRecycle(yErr64Static, y->n, PS_TYPE_F64);
+        p_psMemSetPersistent(yErr64Static, true);
+        p_psMemSetPersistent(yErr64Static->data.V, true);
+
+        for (i=0;i<yErr64Static->n;i++) {
+            yErr64Static->data.F64[i] = 1.0;
+        }
+        yErr64 = yErr64Static;
     } else {
-        myYErr = (psVector *) yErr;
-    }
-
-    // If x==NULL, create an myX vector with x values set to (0:n), and if
+        if (yErr->type.type == PS_TYPE_F64) {
+            yErr64 = (psVector *) yErr;
+        } else if (yErr->type.type == PS_TYPE_F32) {
+            yErr64Static = psVectorRecycle(yErr64Static, yErr->n, PS_TYPE_F64);
+            p_psMemSetPersistent(yErr64Static, true);
+            p_psMemSetPersistent(yErr64Static->data.V, true);
+        } else {
+            // XXX: psError() bad type.
+            psAbort(__func__, "Bad type for yErr64");
+        }
+    }
+
+    // If x==NULL, create an x64 vector with x values set to (0:n), and if
     // this is a Chebyshev polynomial, we must scale to (-1:1).
-
     // XXX: Verify that this is the correct action.
     if (x == NULL) {
-        myX = psVectorAlloc(y->n, y->type.type);
-        mustFreemyX = true;
-
-        if (y->type.type == PS_TYPE_F32) {
-            if (myPoly->type == PS_POLYNOMIAL_ORD) {
-                for (i=0;i<myX->n;i++) {
-                    myX->data.F32[i] = (float) i;
-                }
-            } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-                float min = 0.0;
-                float max = (float) (y->n - 1);
-
-                for (i=0;i<myX->n;i++) {
-                    myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) /
-                                       (0.5 * (max - min));
-                }
-            }
-        } else if (y->type.type == PS_TYPE_F64) {
-            if (myPoly->type == PS_POLYNOMIAL_ORD) {
-                for (i=0;i<myX->n;i++) {
-                    myX->data.F64[i] = (float) i;
-                }
-            } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-                double min = 0.0;
-                double max = (double) (y->n - 1);
-
-                for (i=0;i<myX->n;i++) {
-                    myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) /
-                                       (0.5 * (max - min));
-                }
+        x64Static = psVectorRecycle(x64Static, y->n, PS_TYPE_F64);
+        p_psMemSetPersistent(x64Static, true);
+        p_psMemSetPersistent(x64Static->data.V, true);
+        x64 = x64Static;
+
+        if (myPoly->type == PS_POLYNOMIAL_ORD) {
+            for (i=0;i<x64->n;i++) {
+                x64->data.F64[i] = (float) i;
+            }
+        } else if (myPoly->type == PS_POLYNOMIAL_CHEB) {
+            double min = 0.0;
+            double max = (double) (y->n - 1);
+
+            for (i=0;i<x64->n;i++) {
+                x64->data.F64[i] = (((double) i) - 0.5 * (min + max)) /
+                                   (0.5 * (max - min));
             }
         }
     } else {
-        myX = (psVector *) x;
-    }
-
-    PS_CHECK_VECTOR_SIZE_EQUAL(y, myX);
-    PS_CHECK_VECTOR_SIZE_EQUAL(y, myYErr);
+        if (x->type.type == PS_TYPE_F64) {
+            x64 = (psVector *) x;
+        } else if (x->type.type == PS_TYPE_F32) {
+            x64Static = psVectorRecycle(x64Static, x->n, PS_TYPE_F64);
+            p_psMemSetPersistent(x64Static, true);
+            p_psMemSetPersistent(x64Static->data.V, true);
+            x64 = x64Static;
+        } else {
+            // XXX: psError() bad type.
+            psAbort(__func__, "Bad type for x64");
+        }
+    }
+
+    PS_CHECK_VECTOR_SIZE_EQUAL(y64, x64);
+    PS_CHECK_VECTOR_SIZE_EQUAL(y64, yErr64);
 
     // Call the appropriate vector fitting routine.
     if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, myX, y, myYErr);
+        tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64);
     } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
-        tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, myX, y, myYErr);
+        tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64);
     } else {
         // XXX: psErrorMsg()
@@ -966 +950 @@
 
     // Free any allocated memory.
-    if (mustFreeMyYErr == true) {
-        psFree(myYErr);
-    }
-    if (mustFreemyX == true) {
-        psFree(myX);
-    }
+    /*
+        if ((yErr == NULL) || (yErr->type.type != PS_TYPE_F64)) {
+            psFree(yErr64);
+        }
+        if ((x == NULL) || (x->type.type != PS_TYPE_F64)) {
+            psFree(x64);
+        }
+        if ((y == NULL) || (y->type.type != PS_TYPE_F64)) {
+            psFree(y64);
+        }
+    */
 
     return(myPoly);