Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 4540) +++ trunk/psLib/src/math/psMinimize.c (revision 4568) @@ -9,6 +9,6 @@ * @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,4 +166,6 @@ XXX: spline->knots must be psF32 + +XXXX: Remove this for next code shipment. *****************************************************************************/ /* @@ -262,5 +264,4 @@ PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL); } - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1D() begin ----\n"); @@ -425,4 +426,213 @@ } + + 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;in;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;in;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;idata.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,