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

Timestamp:

Sep 28, 2004, 1:27:37 PM (22 years ago)

Author:

gusciora

Message:

Chebyshev polynomila fitting routines.

File:

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

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-27 23:41:42 $
+ *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-28 23:27:37 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -656 +656 @@
 
 
-/*****************************************************************************
-CreateChebyshevPolys(n): this routine takes as input the required order n,
-and returns as output as a pointer to an array of n psPolynomial1D
-structures, corresponding to the first n Chebyshev polynomials.
- 
-XXX: The output should be static since the Chebyshev polynomials might be
-used frequently and the data structure created here does not contain the
-outer coefficients of the Chebyshev polynomials.
+/******************************************************************************
+p_psVectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
+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).
+ 
+XXX: Use private name?
  *****************************************************************************/
-static psPolynomial1D **CreateChebyshevPolys(int maxChebyPoly)
-{
-    psPolynomial1D **chebPolys = NULL;
-    int i = 0;
-    int j = 0;
-
-    chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *));
-    for (i = 0; i < maxChebyPoly; i++) {
-        chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD);
-    }
-
-    // Create the Chebyshev polynomials.
-    // Polynomial i has i-th order.
-    chebPolys[0]->coeff[0] = 1;
-    chebPolys[1]->coeff[1] = 1;
-    for (i = 2; i < maxChebyPoly; i++) {
-        for (j = 0; j < chebPolys[i - 1]->n; j++) {
-            chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j];
-        }
-        for (j = 0; j < chebPolys[i - 2]->n; j++) {
-            chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
-        }
-    }
-
-    return (chebPolys);
-}
-
-
-
-psPolynomial1D* p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
+psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly,
         const psVector* restrict x,
         const psVector* restrict y,
@@ -702 +675 @@
     int k;
     int n = x->n;
-    psPolynomial1D **chebPolys = CreateChebyshevPolys(myPoly->n);
-    ;
-    float fac;
-    float sum;
-    /*
-        fac = 2.0/((float) n);
-        for (j=0;j<n;j++) {
-            sum = 0.0;
-            for (k=0;k<n;k++) {
-                sum+= y->data.F64[k] *
-                      cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n)); 
-    //            sum+= y->data.F64[k] *
-    //                  cos(M_PI * ((float) j) * (-0.5 + ((float) k)) / ((float) n)) *
-    //                  cos(M_PI * (-0.5 + ((float) k)) / ((float) n)); 
-     
-            }
-            myPoly->coeff[j] = fac * sum;
-     
-        }
-        return(myPoly);
-    */
-
+    psVector *f = psVectorAlloc(n, PS_TYPE_F64);
+    double fac;
+    double sum;
+    // XXX: Use static memory here.
+    psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
+    psScalar *fScalar;
+
+
+    // XXX: These assignments appear too simple to warrant code and
+    // variable declarations.  I retain them here to maintain coherence
+    // with the NR code.
+    double min = -1.0;
+    double max = 1.0;
+    double bma = 0.5 * (max-min);  // 1
+    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
+    // calculate the value of the function we are fitting the Chebyshev
+    // polynomials to at those values of X.  This is a bit tricky since
+    // we don't know that function.  So, we instead do 3-order LaGrange
+    // interpolation at the point X for the psVectors x,y for which we
+    // are fitting this ChebyShev polynomial to.
+
+    for (int i=0;i<n;i++) {
+        // NR 5.8.4
+        double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));
+        double X = (Y + bma + bpa) - 1.0;
+        tmpScalar->data.F32 = (float) 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;
+
+        psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 6,
+                "(x, X, y, f(X)) is (%f, %f, %f, %f)\n",
+                x->data.F64[i], X, y->data.F64[i], f->data.F64[i]);
+    }
+
+    // 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);
-    for (j=0;j<myPoly->n;j++) {
+    for (j=0;j<n;j++) {
         sum = 0.0;
-        for (k=1;k<n;k++) {
-            sum+= (y->data.F64[k] *
-                   psPolynomial1DEval((float) x->data.F64[k], chebPolys[j]));
+        for (k=0;k<n;k++) {
+            sum+= f->data.F64[k] *
+                  cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n));
         }
         myPoly->coeff[j] = fac * sum;
     }
+
+    // XXX: Must free memory.
+    psFree(x32);
+    psFree(y32);
+    psFree(tmpScalar);
+    psFree(fScalar);
     return(myPoly);
 }
 
 /******************************************************************************
-psVectorFitPolynomial1D():  This routine must fit a polynomial of degree
-myPoly to the data points (x, y) and return the coefficients of that
-polynomial, as well as the error for each data point (yErr).
+p_psVectorFitPolynomial1DOrd():  This routine will fit an ordinary
+polynomial of degree myPoly to the data points (x, y) and return the
+coefficients of that polynomial.
  
 XXX: yErr is currently ignored.
@@ -745 +750 @@
 XXX: must add type F32 (currently F64 only).
  
-XXX: Must do: if x==NULL, use the index vector (???).
- 
-XXX: Must do: if yErr==NULL, set all errors equal.
- 
-XXX: type is currently ignored.  Must implement this for Chebyshev
-polynomials.
+XXX: Use private name?
  *****************************************************************************/
-psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
-                                        const psVector* restrict x,
-                                        const psVector* restrict y,
-                                        const psVector* restrict yErr)
+psPolynomial1D* p_psVectorFitPolynomial1DOrd(psPolynomial1D* myPoly,
+        const psVector* restrict x,
+        const psVector* restrict y,
+        const psVector* restrict yErr)
 {
-    // XXX: Create a p_psVectorFitPolynomial1DOrd() function.  Here, we
-    // should only be calling that or the Cheby function.
-
-    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
-        return(p_psVectorFitPolynomial1DCheby(myPoly, x, y, yErr));
-    }
-
     int polyOrder = myPoly->n;
     psImage* A = NULL;
@@ -783 +776 @@
     // }
 
+    // XXX: Some of these are redundant.
     PS_CHECK_NULL_1DPOLY(myPoly);
     PS_CHECK_NULL_VECTOR(y);
     PS_CHECK_EMPTY_VECTOR(y);
-
-    // XXX: Verify that this is the correct action.
-    if (x == NULL) {
-        x = psVectorAlloc(y->n, PS_TYPE_F32);
-        for (i=0;i<x->n;i++) {
-            x->data.F32[i] = (float) i;
-        }
-    }
-
     PS_CHECK_EMPTY_VECTOR(x);
     PS_CHECK_NULL_VECTOR(yErr);
@@ -863 +848 @@
     return (myPoly);
 }
+
+/******************************************************************************
+psVectorFitPolynomial1D():  This routine must fit a 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).
+ *****************************************************************************/
+psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
+                                        const psVector* restrict x,
+                                        const psVector* restrict y,
+                                        const psVector* restrict yErr)
+{
+    bool mustFreeMyYErr = false;
+    bool mustFreemyX = false;
+    int i;
+    psPolynomial1D *tmpPoly;
+    psVector *myX = NULL;
+    psVector *myYErr = NULL;
+
+    PS_CHECK_NULL_1DPOLY(myPoly);
+    PS_CHECK_NULL_VECTOR(y);
+    PS_CHECK_EMPTY_VECTOR(y);
+
+    // 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<yErr->n;i++) {
+                myYErr->data.F32[i] = 1.0;
+            }
+        } else if (y->type.type == PS_TYPE_F64) {
+            for (i=0;i<yErr->n;i++) {
+                myYErr->data.F64[i] = 1.0;
+            }
+        }
+    } else {
+        myYErr = (psVector *) yErr;
+    }
+
+    // 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).
+
+    // 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<yErr->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<yErr->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<yErr->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<yErr->n;i++) {
+                    myX->data.F64[i] = (((float) 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);
+
+    // Call the appropriate vector fitting routine.
+    if (myPoly->type == PS_POLYNOMIAL_CHEB) {
+        tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, myX, y, myYErr);
+    } else if (myPoly->type == PS_POLYNOMIAL_ORD) {
+        tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, myX, y, myYErr);
+    } else {
+        // XXX: psErrorMsg()
+        return(NULL);
+    }
+
+    // Free any allocated memory.
+    if (mustFreeMyYErr == true) {
+        psFree(myYErr);
+    }
+    if (mustFreemyX == true) {
+        psFree(myX);
+    }
+
+    return(myPoly);
+}
+