Index: /trunk/psLib/src/math/psFunctions.c =================================================================== --- /trunk/psLib/src/math/psFunctions.c (revision 4567) +++ /trunk/psLib/src/math/psFunctions.c (revision 4568) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-15 23:19:41 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-16 00:06:32 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1818,88 +1818,123 @@ must exist n+1 points in "knots". -XXX: Ensure that domain[i+1] != domain[i] +XXX: Is this really needed anymore? -XXX: What should be the defualty type for knots be? psF32 is assumed. +XXX: Ensure that knots[i+1] != knots[i] + +XXX: What should be the default type for knots be? psF32 is assumed. *****************************************************************************/ -psSpline1D *psSpline1DAlloc( int numSplines, - int order, - float min, - float max) -{ +psSpline1D *psSpline1DAlloc(int numSplines, + int order, + float min, + float max) +{ + printf("HEY: psSpline1DAlloc()\n"); PS_ASSERT_INT_NONNEGATIVE(numSplines, NULL); PS_ASSERT_INT_NONNEGATIVE(order, NULL); PS_ASSERT_FLOAT_NON_EQUAL(max, min, NULL); - psSpline1D *tmp = NULL; - psS32 i; - psF32 tmpDomain; - psF32 width; - - tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - tmp->n = numSplines; - - tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); - } - - // This should be set by the psVectorFitSpline1D() - tmp->p_psDeriv2 = NULL; - - tmp->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); - width = (max - min) / ((psF32) numSplines); - - tmp->knots->data.F32[0] = min; - tmpDomain = min+width; - for (i=1;iknots->data.F32[i] = tmpDomain; - tmpDomain+= width; - } - tmp->knots->data.F32[numSplines] = max; - - psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); - return(tmp); -} - + psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + tmpSpline->n = numSplines; + + // + // XXX: We might have to allocate single or double polynomials depending on the type + // of the psVector bounds. For now, all knots and spline polynomials are 32-bit. + // + tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (psS32 i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + } + + // This will be computed by psVectorFitSpline1D() + tmpSpline->p_psDeriv2 = NULL; + + // + // XXX:Ensure that the knots are distinct, and monotonic. + // + tmpSpline->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); + if (tmpSpline->knots == NULL) { + printf("BAD\n"); + } else { + printf("GOOD\n"); + } + psF32 width = (max - min) / ((psF32) numSplines); + tmpSpline->knots->data.F32[0] = min; + for (psS32 i=1;iknots->data.F32[i] = min + (width * (psF32) i); + } + tmpSpline->knots->data.F32[numSplines] = max; + + psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); + return(tmpSpline); +} /***************************************************************************** -XXX: What should be the defualty type for knots be? psF32 is assumed. +XXX: Is there a psLib function for this? + *****************************************************************************/ +psVector *PsVectorDup(psVector *in) +{ + psVector *out = psVectorAlloc(in->n, in->type.type); + + if (in->type.type == PS_TYPE_F32) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F32[i] = in->data.F32[i]; + } + } else if (in->type.type == PS_TYPE_F64) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F64[i] = in->data.F64[i]; + } + } else { + printf("XXX: Generate an error here.\n"); + return(NULL); + } + return(out); +} + +/***************************************************************************** +XXX: What should be the default type for knots, spline polys? psF32 is assumed. *****************************************************************************/ psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, int order) { + printf("HEY: psSpline1DAllocGeneric()\n"); PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); PS_ASSERT_VECTOR_NON_EMPTY(bounds, NULL); + PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL); PS_ASSERT_INT_NONNEGATIVE(order, NULL); - psSpline1D *tmpSpline = NULL; - unsigned int i; - unsigned int numSplines; - - tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - - numSplines = bounds->n - 1; + psSpline1D *tmpSpline = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + psS32 numSplines = bounds->n - 1; tmpSpline->n = numSplines; + // + // XXX: We might have to allocate single or double polynomials depending on the type + // of the psVector bounds. For now, all knots and spline polynomials are 32-bit. + // tmpSpline->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); } - // This should be set by the psVectorFitSpline1D() + // This will be computed by psVectorFitSpline1D() tmpSpline->p_psDeriv2 = NULL; - tmpSpline->knots = psVectorAlloc(bounds->n, PS_TYPE_F32); - - for (i=0;in;i++) { - tmpSpline->knots->data.F32[i] = bounds->data.F32[i]; - if (i<(bounds->n-1)) { - if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { - psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n"); - } - } - } - - psMemSetDeallocator(tmpSpline,(psFreeFunc)spline1DFree); + // + // Ensure that all knots are distinct. + // XXX:Ensure that the knots are monotonic. + // + for (psS32 i=0;in-1;i++) { + if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { + psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n"); + return(NULL); + } + } + tmpSpline->knots = PsVectorDup((psVector *) bounds); + if (tmpSpline->knots == NULL) { + printf("BAD\n"); + } else { + printf("GOOD\n"); + } + + psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); return(tmpSpline); } Index: /trunk/psLib/src/math/psFunctions.h =================================================================== --- /trunk/psLib/src/math/psFunctions.h (revision 4567) +++ /trunk/psLib/src/math/psFunctions.h (revision 4568) @@ -12,6 +12,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-12 19:12:01 $ + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-16 00:06:32 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -434,9 +434,9 @@ typedef struct { - unsigned int n; ///< The number of spline polynomials + int n; ///< The number of spline pieces psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials + psVector *knots; ///< The boundaries between each spline piece. Size is n+1. psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. - psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. - psVector *knots; ///< The boundaries between each spline piece. Size is n+1. + psF32 *p_psDomains; ///< The boundaries between each spline piece. Size is n+1. } psSpline1D; Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 4567) +++ /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, Index: /trunk/psLib/src/sys/psErrorCodes.c =================================================================== --- /trunk/psLib/src/sys/psErrorCodes.c (revision 4567) +++ /trunk/psLib/src/sys/psErrorCodes.c (revision 4568) @@ -1,5 +1,4 @@ /** @file psErrorCodes.c - * - * @brief Contains the error codes for the error classes +* @brief Contains the error codes for the error classes * * @ingroup ErrorHandling @@ -7,6 +6,9 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-15 02:33:54 $ + * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-16 00:06:33 $ +======= + * @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-16 00:06:33 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii