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

Timestamp:

Sep 27, 2004, 1:41:42 PM (22 years ago)

Author:

gusciora

Message:

...

File:

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

trunk/psLib/src/dataManip/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-24 20:08:22 $
+ *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-27 23:41:42 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -53 +53 @@
 /* DEFINE STATEMENTS                                                         */
 /*****************************************************************************/
-
-/** Preprocessor macro to generate error on a NULL 1DPolynomial */
-#define PS_CHECK_NULL_1DPOLY(NAME)                                                          \
-if (NAME == NULL || NAME->coeff == NULL) {                                                         \
-    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
-}
-
-/** Preprocessor macro to generate error on a NULL vector */
-#define PS_CHECK_NULL_VECTOR(NAME)                                                          \
-if (NAME == NULL || NAME->data.V == NULL) {                                                         \
-    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
-}
-
-/** Preprocessor macro to generate error for zero length vector */
-#define PS_CHECK_EMPTY_VECTOR(NAME)                                                          \
-if (NAME->n < 1) {                                                                                  \
-    psError(__func__,"Invalid operation: %s has zero n value.", #NAME);                             \
-}
-
-/** Preprocessor macro to generate error on differing size vectors */
-#define PS_CHECK_VECTOR_SIZE_EQUAL(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); \
-}
-
-/** Preprocessor macro to generate error on a NULL image */
-#define PS_CHECK_NULL_IMAGE(NAME)                                                           \
-if (NAME == NULL || NAME->data.V == NULL) {                                                         \
-    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
-}
-
-/** Preprocessor macro to generate error for zero length rows or columns */
-#define PS_CHECK_EMPTY_IMAGE(NAME)                                                           \
-if (NAME->numCols < 1 || NAME->numRows < 1) {                                                       \
-    psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME,              \
-            NAME->numCols, NAME->numRows);                                                          \
-}
 
 /*****************************************************************************/
@@ -692 +655 @@
 }
 
+
+/*****************************************************************************
+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.
+ *****************************************************************************/
+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,
+        const psVector* restrict x,
+        const psVector* restrict y,
+        const psVector* restrict yErr)
+{
+    int j;
+    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);
+    */
+
+    fac = 2.0/((float) n);
+    for (j=0;j<myPoly->n;j++) {
+        sum = 0.0;
+        for (k=1;k<n;k++) {
+            sum+= (y->data.F64[k] *
+                   psPolynomial1DEval((float) x->data.F64[k], chebPolys[j]));
+        }
+        myPoly->coeff[j] = fac * sum;
+    }
+    return(myPoly);
+}
+
 /******************************************************************************
 psVectorFitPolynomial1D():  This routine must fit a polynomial of degree
@@ -704 +748 @@
  
 XXX: Must do: if yErr==NULL, set all errors equal.
+ 
+XXX: type is currently ignored.  Must implement this for Chebyshev
+polynomials.
  *****************************************************************************/
 psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
@@ -710 +757 @@
                                         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;