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

Timestamp:

Jul 15, 2005, 2:06:33 PM (21 years ago)

Author:

gusciora

Message:

...

File:

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

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.124 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-07-12 19:12:01 $
+ *  @version $Revision: 1.125 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-07-16 00:06:32 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -166 +166 @@
  
 XXX: spline->knots must be psF32
+ 
+XXXX: Remove this for next code shipment.
  *****************************************************************************/
 /*
@@ -262 +264 @@
         PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL);
     }
-
     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
             "---- psVectorFitSpline1D() begin ----\n");
@@ -425 +426 @@
 
     }
+
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+            "---- psVectorFitSpline1D() end ----\n");
+    return(mySpline);
+}
+
+
+
+
+
+
+
+
+
+/*****************************************************************************
+psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y
+ 
+xF32 and yF32 are internal psVectors which are used to hold the psF32 versions
+of the input data, if necessary.  xPtr and yPtr are pointers to either xF32 or
+the x argument.  All computation is done on xPtr and yPtr.  xF32 and yF32 will
+simply be psFree() at the end.
+ 
+XXX: nKnots makes no sense.  This number is always equal to the size of the x
+an y vectors.
+ 
+XXX: How do we specify the spline order?  For now, order=3.
+ *****************************************************************************/
+#define PS_XXX_SPLINE_ORDER 3
+psSpline1D *psVectorFitSpline1DNEW(const psVector* x,        ///< Ordinates (or NULL to just use the indices)
+                                   const psVector* y,        ///< Coordinates
+                                   int nKnots)
+{
+    psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n");
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
+    //    PS_ASSERT_INT_EQUAL(y->n, nKnots);
+
+    //
+    // The following code ensures that xPtr points to a psF32 version of the
+    // ordinate data.
+    //
+    psVector *xF32 = NULL;
+    psVector *xPtr = NULL;
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
+        if (PS_TYPE_F64 == x->type.type) {
+            xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+            for (psS32 i = 0 ; i < x->n ; i++) {
+                xF32->data.F32[i] = (psF32) x->data.F64[i];
+            }
+            xPtr = xF32;
+        } else if (PS_TYPE_F32 == x->type.type) {
+            xPtr = (psVector *) x;
+        } else {
+            printf("XXX: Gen Error message: x is wrong type.\n");
+            return(NULL);
+        }
+    } else {
+        // If x==NULL, create an x32 vector with x values set to (0:n).
+        xF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+        for (psS32 i = 0 ; i < x->n ; i++) {
+            xF32->data.F32[i] = (psF32) i;
+        }
+        xPtr = xF32;
+    }
+
+    //
+    // If y is of type psF64, then create a new vector yF32 and convert the
+    // y elements.  Regardless of y's type, we create a yPtr which will be
+    // used in the remainder of this function.
+    //
+    psVector *yF32 = NULL;
+    psVector *yPtr = NULL;
+    if (PS_TYPE_F64 == y->type.type) {
+        yF32 = psVectorAlloc(y->n, PS_TYPE_F32);
+        for (psS32 i = 0 ; i < y->n ; i++) {
+            yF32->data.F32[i] = (psF32) y->data.F64[i];
+        }
+        yPtr = yF32;
+    } else {
+        yPtr = (psVector *) y;
+    }
+
+    psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER);
+
+    psS32 numSplines = nKnots - 1;
+    psF32 tmp;
+    psF32 H;
+    psS32 i;
+    psF32 slope;
+    // XXX: get rid of x32 and y32 (this is from old code)
+    psVector *x32 = xPtr;
+    psVector *y32 = yPtr;
+
+
+    // If these are linear splines, which means their polynomials will have
+    // two coefficients, then we do the simple calculation.
+    if (1 == PS_XXX_SPLINE_ORDER) {
+        for (i=0;i<mySpline->n;i++) {
+            slope = (y32->data.F32[i+1] - y32->data.F32[i]) /
+                    (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]);
+            (mySpline->spline[i])->coeff[0] = y32->data.F32[i] -
+                                              (slope * mySpline->knots->data.F32[i]);
+
+            (mySpline->spline[i])->coeff[1] = slope;
+            psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                    "---- mySpline %d coeffs are (%f, %f)\n", i,
+                    (mySpline->spline[i])->coeff[0],
+                    (mySpline->spline[i])->coeff[1]);
+        }
+        psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4,
+                "---- Exiting psVectorFitSpline1D()()\n");
+        return((psSpline1D *) mySpline);
+    }
+
+    //
+    // Check if these are cubic splines (n==4).  If not, psError.
+    //
+    if (3 != PS_XXX_SPLINE_ORDER) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER);
+        return(NULL);
+    }
+
+    //
+    // If we get here, then we know these are cubic splines.  We first
+    // generate the second derivatives at each data point.
+    //
+    mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32);
+    for (i=0;i<y32->n;i++)
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]);
+
+    //
+    // We generate the coefficients of the spline polynomials.  I can't
+    // concisely explain how this code works.  See above function comments
+    // and Numerical Recipes in C.
+    //
+    for (i=0;i<numSplines;i++) {
+        H = x32->data.F32[i+1] - x32->data.F32[i];
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,
+                "x data (%f - %f) (%f)\n",
+                x32->data.F32[i],
+                x32->data.F32[i+1], H);
+        //
+        // ******** Calculate 0-order term ********
+        //
+        // From (1)
+        (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H);
+        // From (2)
+        ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H);
+        // From (3)
+        tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp-= (x32->data.F32[i+1] / H);
+        tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+        // From (4)
+        tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp+= (x32->data.F32[i] / H);
+        tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[0])+= tmp;
+
+        //
+        // ******** Calculate 1-order term ********
+        //
+        // From (1)
+        (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H;
+        // From (2)
+        ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H);
+        // From (3)
+        tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H);
+        tmp+= (1.0 / H);
+        tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
+        // From (4)
+        tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H);
+        tmp-= (1.0 / H);
+        tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0;
+        ((mySpline->spline[i])->coeff[1])+= tmp;
+
+        //
+        // ******** Calculate 2-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H);
+        // From (4)
+        ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H));
+
+        //
+        // ******** Calculate 3-order term ********
+        //
+        // From (3)
+        (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H);
+        // From (4)
+        ((mySpline->spline[i])->coeff[3])+=  ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H);
+
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]);
+        psTrace(".psLib.dataManip.psVectorFitSpline1D", 6,
+                "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]);
+
+    }
+
+    psFree(xF32);
+    psFree(yF32);
 
     psTrace(".psLib.dataManip.psVectorFitSpline1D", 4,