Index: /trunk/psLib/src/dataManip/psFunctions.c =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.c (revision 1717) +++ /trunk/psLib/src/dataManip/psFunctions.c (revision 1718) @@ -1,3 +1,2 @@ - /** @file psFunctions.c * @@ -8,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-11 20:07:45 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-08 06:00:49 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -32,4 +31,5 @@ #include "psAbort.h" #include "psFunctions.h" +#include "psLogMsg.h" #include @@ -42,5 +42,43 @@ /*****************************************************************************/ -// None +#define PS_CONVERT_VECTOR_F32(IN, OUT) \ +if (IN->type.type != PS_TYPE_F32) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \ + \ + if (IN->type.type == PS_TYPE_F64) { \ + for (int i=0;in;i++) { \ + OUT->data.F32[i] = (float) (IN->data.F64[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + +#define PS_CONVERT_VECTOR_F64(IN, OUT) \ +if (IN->type.type != PS_TYPE_F64) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \ + \ + if (IN->type.type == PS_TYPE_F32) { \ + for (int i=0;in;i++) { \ + OUT->data.F64[i] = (float) (IN->data.F32[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + + /*****************************************************************************/ @@ -107,5 +145,5 @@ } - return (tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); + return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); } @@ -138,11 +176,18 @@ // NOTE: Should I free r as well? - return (gauss); -} + return(gauss); +} + +//XXX: remove this +//typedef enum { +// PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial +// PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial +//} psPolynomialType; /***************************************************************************** This routine must allocate memory for the polynomial structures. *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(int n) +psPolynomial1D* psPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -151,4 +196,6 @@ newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (float *)psAlloc(n * sizeof(float)); @@ -161,8 +208,9 @@ } - return (newPoly); -} - -psPolynomial2D* psPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -172,4 +220,6 @@ newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -191,8 +241,9 @@ } - return (newPoly); -} - -psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -203,4 +254,6 @@ newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -230,8 +283,9 @@ } - return (newPoly); -} - -psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -243,4 +297,6 @@ newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -278,5 +334,5 @@ } - return (newPoly); + return(newPoly); } @@ -359,5 +415,5 @@ XXX: Should the "coeffErr[]" should be used as well? *****************************************************************************/ -float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly) { int loop_x = 0; @@ -371,5 +427,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -383,8 +439,67 @@ } - return (polySum); -} - -float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) + return(polySum); +} + +// XXX: You can do this without having to psAlloc() vector d. +float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + psVector *d; + int n; + int i; + float tmp; + + n = myPoly->n; + d = psVectorAlloc(n, PS_TYPE_F32); + d->data.F32[n-1] = myPoly->coeff[n-1]; + d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2]; + for (i=n-3;i>=1;i--) { + d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - + (d->data.F32[i+2]) + + (myPoly->coeff[i]); + } + + tmp = (x * d->data.F32[1]) - + (d->data.F32[2]) + + (0.5 * myPoly->coeff[0]); + + psFree(d); + return(tmp); +} + +float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial1DEvalVector(const psVector *x, + const psPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F32(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmp->data.F32[i] = psPolynomial1DEval(x->data.F32[i], myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + return(tmp); +} + + +float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) { int loop_x = 0; @@ -403,8 +518,61 @@ } - return (polySum); -} - -float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) + return(polySum); +} + +float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + return(0.0); +} + +float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + + +psVector *psPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial2DEval(myX->data.F32[i], + myY->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + return(tmp); +} + + + +float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) { int loop_x = 0; @@ -429,8 +597,73 @@ } - return (polySum); -} - -float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) + return(polySum); +} + +float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + return(0.0); +} + +float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial3DEval(myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + + +float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) { int loop_w = 0; @@ -461,8 +694,83 @@ } - return (polySum); -} - -psDPolynomial1D* psDPolynomial1DAlloc(int n) + return(polySum); +} + +float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + return(0.0); +} + +float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(w, myW); + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial4DEval(myW->data.F32[i], + myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F32) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + +psDPolynomial1D* psDPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -471,4 +779,6 @@ newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (double *)psAlloc(n * sizeof(double)); @@ -481,8 +791,9 @@ } - return (newPoly); -} - -psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -492,4 +803,6 @@ newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -511,8 +824,9 @@ } - return (newPoly); -} - -psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -523,4 +837,6 @@ newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -550,8 +866,9 @@ } - return (newPoly); -} - -psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -563,4 +880,6 @@ newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -598,5 +917,5 @@ } - return (newPoly); + return(newPoly); } @@ -673,5 +992,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly) { int loop_x = 0; @@ -681,5 +1000,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -689,8 +1008,47 @@ } - return (polySum); -} - -double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + return(0.0); +} + +double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial1DEvalVector(const psVector *x, + const psDPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F64(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F64); + for (i=0;in;i++) { + tmp->data.F64[i] = psDPolynomial1DEval(x->data.F64[i], myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + return(tmp); +} + + + +double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) { int loop_x = 0; @@ -709,8 +1067,60 @@ } - return (polySum); -} - -double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + return(0.0); +} + +double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psDPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial2DEval(myX->data.F64[i], + myY->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + return(tmp); +} + + + +double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) { int loop_x = 0; @@ -735,8 +1145,75 @@ } - return (polySum); -} - -double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + return(0.0); +} + +double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial3DEval(myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + + + + + +double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) { int loop_w = 0; @@ -767,4 +1244,222 @@ } - return (polySum); -} + return(polySum); +} + +double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + return(0.0); +} + +double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(w, myW); + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial4DEval(myW->data.F64[i], + myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F64) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + +//typedef struct { +// int n; +// psPolynomial1D **spline; +// float *domains; +//} psSpline1D; + +/***************************************************************************** + + NOTE: "n" specifies the number of spline polynomials. Therefore, there + must exist n+1 points in "domains". + *****************************************************************************/ +psSpline1D *psSpline1DAlloc(int numSplines, + int order, + float min, + float max) +{ + psSpline1D *tmp = NULL; + int i; + float tmpDomain; + float width; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + tmp->n = numSplines; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); + width = (max - min) / ((float) numSplines); + tmpDomain = min; + for (i=0;idomains)[i] = tmpDomain; + tmpDomain+= width; + } + + return(tmp); +} + +psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, + int order) +{ + psSpline1D *tmp = NULL; + int i; + int numSplines; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + numSplines = bounds->n - 1; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); + for (i=0;in;i++) { + (tmp->domains)[i] = bounds->data.F32[i]; + } + + return(tmp); +} + +/***************************************************************************** +VectorBinDisect(): This is a private function which takes as input a vector +of floating point data as well as a single floating point values. The input +vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i). +This routine does a binary disection of the vector and returns "i" such +that (v[i] <= x <= v[i+1). If x lies outside the range of v[], +then this routine prints a warning message and returns -1. + *****************************************************************************/ +int VectorBinDisect(float *bins, + int numBins, + float x) +{ + int min; + int max; + int mid; + + if ((x < bins[0]) || + (x > bins[numBins-1])) { + psLogMsg(__func__, PS_LOG_WARN, + "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", + x, bins[0], bins[numBins-1]); + return(-1); + } + + min = 0; + max = numBins-1; + mid = (max-min)/2; + + while (min != max) { + mid = (max-min)/2; + + if (x < bins[mid]) { + max = mid; + } else { + min = mid; + } + } + + return(min); +} + +float psSpline1DEval(const psSpline1D *spline, + float x) +{ + int binNum; + int n; + + n = spline->n; + binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); + if (binNum == -1) { + psLogMsg(__func__, PS_LOG_WARN, + "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", + x, (spline->domains)[0], + (spline->domains)[n-1]); + + if (x < (spline->domains)[0]) { + return(psPolynomial1DEval(x, spline->spline[0])); + } else if (x > (spline->domains)[n-1]) { + return(psPolynomial1DEval(x, spline->spline[n-1])); + } + } + + return(psPolynomial1DEval(x, spline->spline[binNum])); +} + +psVector *psSpline1DEvalVector(const psVector *x, + const psSpline1D *spline) +{ + int i; + psVector *tmpVector; + + tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]); + } + + return(tmpVector); +} + Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 1717) +++ /trunk/psLib/src/math/psPolynomial.c (revision 1718) @@ -1,3 +1,2 @@ - /** @file psFunctions.c * @@ -8,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-11 20:07:45 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-08 06:00:49 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -32,4 +31,5 @@ #include "psAbort.h" #include "psFunctions.h" +#include "psLogMsg.h" #include @@ -42,5 +42,43 @@ /*****************************************************************************/ -// None +#define PS_CONVERT_VECTOR_F32(IN, OUT) \ +if (IN->type.type != PS_TYPE_F32) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \ + \ + if (IN->type.type == PS_TYPE_F64) { \ + for (int i=0;in;i++) { \ + OUT->data.F32[i] = (float) (IN->data.F64[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + +#define PS_CONVERT_VECTOR_F64(IN, OUT) \ +if (IN->type.type != PS_TYPE_F64) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \ + \ + if (IN->type.type == PS_TYPE_F32) { \ + for (int i=0;in;i++) { \ + OUT->data.F64[i] = (float) (IN->data.F32[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + + /*****************************************************************************/ @@ -107,5 +145,5 @@ } - return (tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); + return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); } @@ -138,11 +176,18 @@ // NOTE: Should I free r as well? - return (gauss); -} + return(gauss); +} + +//XXX: remove this +//typedef enum { +// PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial +// PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial +//} psPolynomialType; /***************************************************************************** This routine must allocate memory for the polynomial structures. *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(int n) +psPolynomial1D* psPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -151,4 +196,6 @@ newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (float *)psAlloc(n * sizeof(float)); @@ -161,8 +208,9 @@ } - return (newPoly); -} - -psPolynomial2D* psPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -172,4 +220,6 @@ newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -191,8 +241,9 @@ } - return (newPoly); -} - -psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -203,4 +254,6 @@ newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -230,8 +283,9 @@ } - return (newPoly); -} - -psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -243,4 +297,6 @@ newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -278,5 +334,5 @@ } - return (newPoly); + return(newPoly); } @@ -359,5 +415,5 @@ XXX: Should the "coeffErr[]" should be used as well? *****************************************************************************/ -float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly) { int loop_x = 0; @@ -371,5 +427,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -383,8 +439,67 @@ } - return (polySum); -} - -float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) + return(polySum); +} + +// XXX: You can do this without having to psAlloc() vector d. +float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + psVector *d; + int n; + int i; + float tmp; + + n = myPoly->n; + d = psVectorAlloc(n, PS_TYPE_F32); + d->data.F32[n-1] = myPoly->coeff[n-1]; + d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2]; + for (i=n-3;i>=1;i--) { + d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - + (d->data.F32[i+2]) + + (myPoly->coeff[i]); + } + + tmp = (x * d->data.F32[1]) - + (d->data.F32[2]) + + (0.5 * myPoly->coeff[0]); + + psFree(d); + return(tmp); +} + +float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial1DEvalVector(const psVector *x, + const psPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F32(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmp->data.F32[i] = psPolynomial1DEval(x->data.F32[i], myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + return(tmp); +} + + +float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) { int loop_x = 0; @@ -403,8 +518,61 @@ } - return (polySum); -} - -float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) + return(polySum); +} + +float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + return(0.0); +} + +float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + + +psVector *psPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial2DEval(myX->data.F32[i], + myY->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + return(tmp); +} + + + +float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) { int loop_x = 0; @@ -429,8 +597,73 @@ } - return (polySum); -} - -float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) + return(polySum); +} + +float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + return(0.0); +} + +float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial3DEval(myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + + +float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) { int loop_w = 0; @@ -461,8 +694,83 @@ } - return (polySum); -} - -psDPolynomial1D* psDPolynomial1DAlloc(int n) + return(polySum); +} + +float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + return(0.0); +} + +float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(w, myW); + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial4DEval(myW->data.F32[i], + myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F32) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + +psDPolynomial1D* psDPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -471,4 +779,6 @@ newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (double *)psAlloc(n * sizeof(double)); @@ -481,8 +791,9 @@ } - return (newPoly); -} - -psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -492,4 +803,6 @@ newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -511,8 +824,9 @@ } - return (newPoly); -} - -psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -523,4 +837,6 @@ newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -550,8 +866,9 @@ } - return (newPoly); -} - -psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -563,4 +880,6 @@ newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -598,5 +917,5 @@ } - return (newPoly); + return(newPoly); } @@ -673,5 +992,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly) { int loop_x = 0; @@ -681,5 +1000,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -689,8 +1008,47 @@ } - return (polySum); -} - -double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + return(0.0); +} + +double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial1DEvalVector(const psVector *x, + const psDPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F64(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F64); + for (i=0;in;i++) { + tmp->data.F64[i] = psDPolynomial1DEval(x->data.F64[i], myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + return(tmp); +} + + + +double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) { int loop_x = 0; @@ -709,8 +1067,60 @@ } - return (polySum); -} - -double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + return(0.0); +} + +double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psDPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial2DEval(myX->data.F64[i], + myY->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + return(tmp); +} + + + +double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) { int loop_x = 0; @@ -735,8 +1145,75 @@ } - return (polySum); -} - -double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + return(0.0); +} + +double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial3DEval(myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + + + + + +double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) { int loop_w = 0; @@ -767,4 +1244,222 @@ } - return (polySum); -} + return(polySum); +} + +double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + return(0.0); +} + +double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(w, myW); + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial4DEval(myW->data.F64[i], + myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F64) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + +//typedef struct { +// int n; +// psPolynomial1D **spline; +// float *domains; +//} psSpline1D; + +/***************************************************************************** + + NOTE: "n" specifies the number of spline polynomials. Therefore, there + must exist n+1 points in "domains". + *****************************************************************************/ +psSpline1D *psSpline1DAlloc(int numSplines, + int order, + float min, + float max) +{ + psSpline1D *tmp = NULL; + int i; + float tmpDomain; + float width; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + tmp->n = numSplines; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); + width = (max - min) / ((float) numSplines); + tmpDomain = min; + for (i=0;idomains)[i] = tmpDomain; + tmpDomain+= width; + } + + return(tmp); +} + +psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, + int order) +{ + psSpline1D *tmp = NULL; + int i; + int numSplines; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + numSplines = bounds->n - 1; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); + for (i=0;in;i++) { + (tmp->domains)[i] = bounds->data.F32[i]; + } + + return(tmp); +} + +/***************************************************************************** +VectorBinDisect(): This is a private function which takes as input a vector +of floating point data as well as a single floating point values. The input +vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i). +This routine does a binary disection of the vector and returns "i" such +that (v[i] <= x <= v[i+1). If x lies outside the range of v[], +then this routine prints a warning message and returns -1. + *****************************************************************************/ +int VectorBinDisect(float *bins, + int numBins, + float x) +{ + int min; + int max; + int mid; + + if ((x < bins[0]) || + (x > bins[numBins-1])) { + psLogMsg(__func__, PS_LOG_WARN, + "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", + x, bins[0], bins[numBins-1]); + return(-1); + } + + min = 0; + max = numBins-1; + mid = (max-min)/2; + + while (min != max) { + mid = (max-min)/2; + + if (x < bins[mid]) { + max = mid; + } else { + min = mid; + } + } + + return(min); +} + +float psSpline1DEval(const psSpline1D *spline, + float x) +{ + int binNum; + int n; + + n = spline->n; + binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); + if (binNum == -1) { + psLogMsg(__func__, PS_LOG_WARN, + "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", + x, (spline->domains)[0], + (spline->domains)[n-1]); + + if (x < (spline->domains)[0]) { + return(psPolynomial1DEval(x, spline->spline[0])); + } else if (x > (spline->domains)[n-1]) { + return(psPolynomial1DEval(x, spline->spline[n-1])); + } + } + + return(psPolynomial1DEval(x, spline->spline[binNum])); +} + +psVector *psSpline1DEvalVector(const psVector *x, + const psSpline1D *spline) +{ + int i; + psVector *tmpVector; + + tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]); + } + + return(tmpVector); +} + Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 1717) +++ /trunk/psLib/src/math/psSpline.c (revision 1718) @@ -1,3 +1,2 @@ - /** @file psFunctions.c * @@ -8,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-11 20:07:45 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-08 06:00:49 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -32,4 +31,5 @@ #include "psAbort.h" #include "psFunctions.h" +#include "psLogMsg.h" #include @@ -42,5 +42,43 @@ /*****************************************************************************/ -// None +#define PS_CONVERT_VECTOR_F32(IN, OUT) \ +if (IN->type.type != PS_TYPE_F32) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F32); \ + \ + if (IN->type.type == PS_TYPE_F64) { \ + for (int i=0;in;i++) { \ + OUT->data.F32[i] = (float) (IN->data.F64[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + +#define PS_CONVERT_VECTOR_F64(IN, OUT) \ +if (IN->type.type != PS_TYPE_F64) { \ + psLogMsg(__func__, PS_LOG_WARN, \ + "Input vector has incorrect type (%d).\n", \ + IN->type); \ + OUT = psVectorAlloc(IN->n, PS_TYPE_F64); \ + \ + if (IN->type.type == PS_TYPE_F32) { \ + for (int i=0;in;i++) { \ + OUT->data.F64[i] = (float) (IN->data.F32[i]); \ + } \ + }\ + else { \ + psAbort(__func__, "Wrong type.\n"); \ + } \ +} else { \ + OUT = (psVector *) IN; \ +} + + /*****************************************************************************/ @@ -107,5 +145,5 @@ } - return (tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); + return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); } @@ -138,11 +176,18 @@ // NOTE: Should I free r as well? - return (gauss); -} + return(gauss); +} + +//XXX: remove this +//typedef enum { +// PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial +// PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial +//} psPolynomialType; /***************************************************************************** This routine must allocate memory for the polynomial structures. *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(int n) +psPolynomial1D* psPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -151,4 +196,6 @@ newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (float *)psAlloc(n * sizeof(float)); @@ -161,8 +208,9 @@ } - return (newPoly); -} - -psPolynomial2D* psPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psPolynomial2D* psPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -172,4 +220,6 @@ newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -191,8 +241,9 @@ } - return (newPoly); -} - -psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial3D* psPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -203,4 +254,6 @@ newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -230,8 +283,9 @@ } - return (newPoly); -} - -psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psPolynomial4D* psPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -243,4 +297,6 @@ newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) polynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -278,5 +334,5 @@ } - return (newPoly); + return(newPoly); } @@ -359,5 +415,5 @@ XXX: Should the "coeffErr[]" should be used as well? *****************************************************************************/ -float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +float p_psOrgPolynomial1DEval(float x, const psPolynomial1D* myPoly) { int loop_x = 0; @@ -371,5 +427,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -383,8 +439,67 @@ } - return (polySum); -} - -float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) + return(polySum); +} + +// XXX: You can do this without having to psAlloc() vector d. +float p_psChebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + psVector *d; + int n; + int i; + float tmp; + + n = myPoly->n; + d = psVectorAlloc(n, PS_TYPE_F32); + d->data.F32[n-1] = myPoly->coeff[n-1]; + d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]) + myPoly->coeff[n-2]; + for (i=n-3;i>=1;i--) { + d->data.F32[i] = (2.0 * x * d->data.F32[i+1]) - + (d->data.F32[i+2]) + + (myPoly->coeff[i]); + } + + tmp = (x * d->data.F32[1]) - + (d->data.F32[2]) + + (0.5 * myPoly->coeff[0]); + + psFree(d); + return(tmp); +} + +float psPolynomial1DEval(float x, const psPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial1DEvalVector(const psVector *x, + const psPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F32(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmp->data.F32[i] = psPolynomial1DEval(x->data.F32[i], myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + return(tmp); +} + + +float p_psOrgPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) { int loop_x = 0; @@ -403,8 +518,61 @@ } - return (polySum); -} - -float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) + return(polySum); +} + +float p_psChebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + return(0.0); +} + +float psPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + + +psVector *psPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial2DEval(myX->data.F32[i], + myY->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + return(tmp); +} + + + +float p_psOrgPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) { int loop_x = 0; @@ -429,8 +597,73 @@ } - return (polySum); -} - -float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) + return(polySum); +} + +float p_psChebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + return(0.0); +} + +float psPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial3DEval(myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + + +float p_psOrgPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) { int loop_w = 0; @@ -461,8 +694,83 @@ } - return (polySum); -} - -psDPolynomial1D* psDPolynomial1DAlloc(int n) + return(polySum); +} + +float p_psChebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + return(0.0); +} + +float psPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F32(w, myW); + PS_CONVERT_VECTOR_F32(x, myX); + PS_CONVERT_VECTOR_F32(y, myY); + PS_CONVERT_VECTOR_F32(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F32); + for (i=0;idata.F32[i] = psPolynomial4DEval(myW->data.F32[i], + myX->data.F32[i], + myY->data.F32[i], + myZ->data.F32[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F32) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F32) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F32) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F32) { + psFree(myZ); + } + return(tmp); +} + + + + + +psDPolynomial1D* psDPolynomial1DAlloc(int n, + psPolynomialType type) { int i = 0; @@ -471,4 +779,6 @@ newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial1DFree); + + newPoly->type = type; newPoly->n = n; newPoly->coeff = (double *)psAlloc(n * sizeof(double)); @@ -481,8 +791,9 @@ } - return (newPoly); -} - -psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY) + return(newPoly); +} + +psDPolynomial2D* psDPolynomial2DAlloc(int nX, int nY, + psPolynomialType type) { int x = 0; @@ -492,4 +803,6 @@ newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial2DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -511,8 +824,9 @@ } - return (newPoly); -} - -psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial3D* psDPolynomial3DAlloc(int nX, int nY, int nZ, + psPolynomialType type) { int x = 0; @@ -523,4 +837,6 @@ newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial3DFree); + + newPoly->type = type; newPoly->nX = nX; newPoly->nY = nY; @@ -550,8 +866,9 @@ } - return (newPoly); -} - -psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ) + return(newPoly); +} + +psDPolynomial4D* psDPolynomial4DAlloc(int nW, int nX, int nY, int nZ, + psPolynomialType type) { int w = 0; @@ -563,4 +880,6 @@ newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); p_psMemSetDeallocator(newPoly, (psFreeFcn) dPolynomial4DFree); + + newPoly->type = type; newPoly->nW = nW; newPoly->nX = nX; @@ -598,5 +917,5 @@ } - return (newPoly); + return(newPoly); } @@ -673,5 +992,5 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +double p_psDOrgPolynomial1DEval(double x, const psDPolynomial1D* myPoly) { int loop_x = 0; @@ -681,5 +1000,5 @@ // NOTE: Do we want to flag this case? if (myPoly->n == 0) { - return (1.0); + return(1.0); } @@ -689,8 +1008,47 @@ } - return (polySum); -} - -double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + return(0.0); +} + +double psDPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial1DEval(x, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial1DEval(x, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial1DEvalVector(const psVector *x, + const psDPolynomial1D *myPoly) +{ + psVector *tmp; + psVector *myX; + int i; + + PS_CONVERT_VECTOR_F64(x, myX); + + tmp = psVectorAlloc(x->n, PS_TYPE_F64); + for (i=0;in;i++) { + tmp->data.F64[i] = psDPolynomial1DEval(x->data.F64[i], myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + return(tmp); +} + + + +double p_psDOrgPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) { int loop_x = 0; @@ -709,8 +1067,60 @@ } - return (polySum); -} - -double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + return(0.0); +} + +double psDPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial2DEval(x, y, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial2DEval(x, y, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial2DEvalVector(const psVector *x, + const psVector *y, + const psDPolynomial2D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial2DEval(myX->data.F64[i], + myY->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + return(tmp); +} + + + +double p_psDOrgPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) { int loop_x = 0; @@ -735,8 +1145,75 @@ } - return (polySum); -} - -double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) + return(polySum); +} + +double p_psDChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + return(0.0); +} + +double psDPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial3DEval(x, y, z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial3DEval(x, y, z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial3DEvalVector(const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial3D *myPoly) +{ + psVector *tmp; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial3DEval(myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + + + + + +double p_psDOrgPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) { int loop_w = 0; @@ -767,4 +1244,222 @@ } - return (polySum); -} + return(polySum); +} + +double p_psDChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + return(0.0); +} + +double psDPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +{ + if (myPoly->type == PS_POLYNOMIAL_ORD) { + return(p_psDOrgPolynomial4DEval(w,x,y,z, myPoly)); + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + return(p_psDChebPolynomial4DEval(w,x,y,z, myPoly)); + } else { + psAbort(__func__, "Unknown polynomial type 0x%x\n", myPoly->type); + } + return(0.0); +} + +psVector *psDPolynomial4DEvalVector(const psVector *w, + const psVector *x, + const psVector *y, + const psVector *z, + const psDPolynomial4D *myPoly) +{ + psVector *tmp; + psVector *myW; + psVector *myX; + psVector *myY; + psVector *myZ; + int i; + int vecLen=x->n; + + PS_CONVERT_VECTOR_F64(w, myW); + PS_CONVERT_VECTOR_F64(x, myX); + PS_CONVERT_VECTOR_F64(y, myY); + PS_CONVERT_VECTOR_F64(z, myZ); + vecLen=w->n; + if (y->n < vecLen) { + vecLen = y->n; + } + if (x->n < vecLen) { + vecLen = x->n; + } + if (z->n < vecLen) { + vecLen = z->n; + } + + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + for (i=0;idata.F64[i] = psDPolynomial4DEval(myW->data.F64[i], + myX->data.F64[i], + myY->data.F64[i], + myZ->data.F64[i], + myPoly); + } + + if (w->type.type != PS_TYPE_F64) { + psFree(myW); + } + if (x->type.type != PS_TYPE_F64) { + psFree(myX); + } + if (y->type.type != PS_TYPE_F64) { + psFree(myY); + } + if (z->type.type != PS_TYPE_F64) { + psFree(myZ); + } + return(tmp); +} + + + + +//typedef struct { +// int n; +// psPolynomial1D **spline; +// float *domains; +//} psSpline1D; + +/***************************************************************************** + + NOTE: "n" specifies the number of spline polynomials. Therefore, there + must exist n+1 points in "domains". + *****************************************************************************/ +psSpline1D *psSpline1DAlloc(int numSplines, + int order, + float min, + float max) +{ + psSpline1D *tmp = NULL; + int i; + float tmpDomain; + float width; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + tmp->n = numSplines; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); + width = (max - min) / ((float) numSplines); + tmpDomain = min; + for (i=0;idomains)[i] = tmpDomain; + tmpDomain+= width; + } + + return(tmp); +} + +psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, + int order) +{ + psSpline1D *tmp = NULL; + int i; + int numSplines; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + numSplines = bounds->n - 1; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order, PS_POLYNOMIAL_ORD); + } + + tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); + for (i=0;in;i++) { + (tmp->domains)[i] = bounds->data.F32[i]; + } + + return(tmp); +} + +/***************************************************************************** +VectorBinDisect(): This is a private function which takes as input a vector +of floating point data as well as a single floating point values. The input +vector values are assumed to be non-decreasing (v[i-1] <= v[j] for all j>=i). +This routine does a binary disection of the vector and returns "i" such +that (v[i] <= x <= v[i+1). If x lies outside the range of v[], +then this routine prints a warning message and returns -1. + *****************************************************************************/ +int VectorBinDisect(float *bins, + int numBins, + float x) +{ + int min; + int max; + int mid; + + if ((x < bins[0]) || + (x > bins[numBins-1])) { + psLogMsg(__func__, PS_LOG_WARN, + "VectorBinDisect(): ordinate %f is outside vector range (%f - %f).", + x, bins[0], bins[numBins-1]); + return(-1); + } + + min = 0; + max = numBins-1; + mid = (max-min)/2; + + while (min != max) { + mid = (max-min)/2; + + if (x < bins[mid]) { + max = mid; + } else { + min = mid; + } + } + + return(min); +} + +float psSpline1DEval(const psSpline1D *spline, + float x) +{ + int binNum; + int n; + + n = spline->n; + binNum = VectorBinDisect(spline->domains, (spline->n)+1, x); + if (binNum == -1) { + psLogMsg(__func__, PS_LOG_WARN, + "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", + x, (spline->domains)[0], + (spline->domains)[n-1]); + + if (x < (spline->domains)[0]) { + return(psPolynomial1DEval(x, spline->spline[0])); + } else if (x > (spline->domains)[n-1]) { + return(psPolynomial1DEval(x, spline->spline[n-1])); + } + } + + return(psPolynomial1DEval(x, spline->spline[binNum])); +} + +psVector *psSpline1DEvalVector(const psVector *x, + const psSpline1D *spline) +{ + int i; + psVector *tmpVector; + + tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); + for (i=0;in;i++) { + tmpVector->data.F32[i] = psSpline1DEval(spline, x->data.F32[i]); + } + + return(tmpVector); +} + Index: /trunk/psLib/src/sys/psTrace.c =================================================================== --- /trunk/psLib/src/sys/psTrace.c (revision 1717) +++ /trunk/psLib/src/sys/psTrace.c (revision 1718) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-08 00:09:06 $ + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-08 05:59:41 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -27,4 +27,12 @@ not the same thing as a node's "level", which corresponds to the trace level of that node. + +I think the following is the correct behavior, but not sure: + PS_UNKNOWN_TRACE_LEVEL: We never set the level of a component to this + value. This value is only used when psTraceGetLevel is called with + a bad component name, or if the component root is undefined, I think. + + PS_DEFAULT_TRACE_LEVEL: This should only be used when adding the + intermediate components of a long name. Ie. the "B" in .A.B.C *****************************************************************************/ @@ -40,8 +48,5 @@ #include "psString.h" #include "psError.h" -#include "psLogMsg.h" - -#include "psSysUtilsErrors.h" -#define ERRORNAME_PREFIX PS_ERRORNAME_DOMAIN ".psTrace." +#include "psAbort.h" static p_psComponent* cRoot = NULL; // The root of the trace component @@ -98,4 +103,6 @@ /***************************************************************************** Set all trace levels to zero. + +XXX: Currently, no function calls this routine. *****************************************************************************/ void p_psTraceReset(p_psComponent* currentNode) @@ -110,7 +117,6 @@ for (i = 0; i < currentNode->n; i++) { if (NULL == currentNode->subcomp[i]) { - psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, - PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, - i, currentNode->name); + psError(__func__, + "Sub-component %d of node %s in the trace tree is NULL.\n", i, currentNode->name); } else { p_psTraceReset(currentNode->subcomp[i]); @@ -125,5 +131,5 @@ void psTraceReset() { - psFree(cRoot); + componentFree(cRoot); cRoot = NULL; } @@ -135,5 +141,5 @@ to ANSI-C. *****************************************************************************/ -static bool componentAdd(const char *addNodeName, int level) +static void componentAdd(const char *addNodeName, int level) { int i = 0; // Loop index variable. @@ -146,7 +152,5 @@ // XXX: Verify that this is the correct behavior. if (strcmp("", addNodeName) == 0) { - psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_NULL,true, - PS_ERRORTEXT_psTrace_ADD_NULL_COMPONENT); - return false; + psAbort(__func__, "Failed to add null component to trace tree.\n"); } @@ -154,12 +158,10 @@ if (strcmp(".", addNodeName) == 0) { cRoot->level = level; - return true; + return; } if (addNodeName[0] != '.') { - psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTrace_MALFORMED_COMPONENT_NAME, - addNodeName); - return false; + printf("ERROR: failed to add %s to the root component tree.\n", addNodeName); + exit(1); } @@ -186,16 +188,21 @@ currentNode->subcomp = psRealloc(currentNode->subcomp, (currentNode->n + 1) * sizeof(p_psComponent* )); - currentNode->n++; + currentNode->n = (currentNode->n) + 1; if (pname == NULL) { - currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, level); + // This is the final component to add. + currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, level); } else { - currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, currentNode->level); - } - currentNode = currentNode->subcomp[currentNode->n - 1]; - } - } - - return true; + // We are adding an intermediate component. The trace level + // is not defined. An undefined trace level inherits the + // trace level of it's parent. However, we do not set that + // specifically here since that would inheritance to be a + // static, one-time, type of behavior. + + currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, PS_DEFAULT_TRACE_LEVEL); + } + currentNode = currentNode->subcomp[(currentNode->n) - 1]; + } + } } @@ -211,6 +218,6 @@ zero *****************************************************************************/ -bool psTraceSetLevel(const char *comp, // component of interest - int level) // desired trace level +int psTraceSetLevel(const char *comp, // component of interest + int level) // desired trace level { // If the root component tree does not exist, then initialize it. @@ -219,11 +226,8 @@ } // Add the new component to the component tree. - if ( !componentAdd(comp, level) ) { - psErrorMsg(ERRORNAME_PREFIX "psTraceSetLevel", PS_ERR_UNKNOWN, false, - PS_ERRORTEXT_psTrace_FAILED_TO_ADD_COMPONENT,level,comp); - return false; - } + componentAdd(comp, level); + // return 0 on success. - return true; + return 0; } @@ -249,4 +253,5 @@ p_psComponent* currentNode = cRoot; int i = 0; + int defaultLevel = 0; if (NULL == currentNode) { @@ -262,4 +267,5 @@ } + defaultLevel = cRoot->level; strcpy(name, aname); pname = &name[1]; @@ -268,22 +274,31 @@ for (i = 0; i < currentNode->n; i++) { if (NULL == currentNode->subcomp[i]) { - psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, - PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, - i, currentNode->name); + psError(__func__, + "Sub-component %d of node %s in trace tree is NULL.\n", i, currentNode->name); } if (strcmp(currentNode->subcomp[i]->name, firstComponent) == 0) { currentNode = currentNode->subcomp[i]; + // For level inheritance purpose, we save the level of this + // component if it is not DEFAULT. + if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { + defaultLevel = currentNode->level; + } + // Determine if thisis the last component: if (pname == NULL) { - return (currentNode->level); + if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { + return (currentNode->level); + } else { + return(PS_DEFAULT_TRACE_LEVEL); + } } } } } - return (PS_UNKNOWN_TRACE_LEVEL); -} - -/***************************************************************************** - psGetTraceLevel() + return(defaultLevel); +} + +/***************************************************************************** + psTraceLevelGet() Return a trace level of "name" in the root component tree. If the exact string of components in "name" does not exist in the root @@ -303,4 +318,37 @@ // Search the component root tree, determine the trace level. return (doGetTraceLevel(name)); +} + +/***************************************************************************** + doPrintTraceLevelsOld() + This function recursively searches the component tree supplied by the + parameter "comp" and prints the name and level of each component. + Inputs: + comp: a node in the component tree. + level: the level of that node + Outputs: + none + Returns: + null + *****************************************************************************/ +static void doPrintTraceLevelsOld(const p_psComponent* comp, int depth) +{ + int i = 0; + + if (comp->name[0] == '\0') { + printf("%*s%-*s %d\n", depth, "", 20 - depth, + "(root)", (comp->level == PS_DEFAULT_TRACE_LEVEL) ? 0 : comp->level); + } else { + if (comp->level == PS_DEFAULT_TRACE_LEVEL) { + printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); + } else { + printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, + comp->level); + } + } + + for (i = 0; i < comp->n; i++) { + doPrintTraceLevelsOld(comp->subcomp[i], depth + 1); + } } @@ -317,23 +365,31 @@ null *****************************************************************************/ -static void doPrintTraceLevels(const p_psComponent* comp, int depth) +static void doPrintTraceLevels(const p_psComponent* comp, + int depth, + int defLevel) { int i = 0; if (comp->name[0] == '\0') { - printf("%*s%-*s %d\n", depth, "", 20 - depth, - "(root)", (comp->level == PS_UNKNOWN_TRACE_LEVEL) ? 0 : comp->level); + psAbort(__func__, "component name is NULL\n"); } else { - if (comp->level == PS_UNKNOWN_TRACE_LEVEL) { - printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); + if (comp->level == PS_DEFAULT_TRACE_LEVEL) { + printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, + defLevel); } else { - printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level); + printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, + comp->level); } } for (i = 0; i < comp->n; i++) { - doPrintTraceLevels(comp->subcomp[i], depth + 1); - } -} + if (comp->level == PS_DEFAULT_TRACE_LEVEL) { + doPrintTraceLevels(comp->subcomp[i], depth + 1, defLevel); + } else { + doPrintTraceLevels(comp->subcomp[i], depth + 1, comp->level); + } + } +} + /***************************************************************************** @@ -353,5 +409,5 @@ } - doPrintTraceLevels(cRoot, 0); + doPrintTraceLevels(cRoot, 0, PS_DEFAULT_TRACE_LEVEL); } @@ -379,8 +435,5 @@ if (NULL == comp) { - psErrorMsg(ERRORNAME_PREFIX "psTrace", PS_ERR_BAD_PARAMETER_NULL, true, - PS_ERRORTEXT_psTrace_NULL_TRACETREE, - __func__); - return; + psError(__func__, "p_psTrace() called on a NULL trace level tree\n"); } // Only display this message if it's trace level is less than the level Index: /trunk/psLib/src/sysUtils/psTrace.c =================================================================== --- /trunk/psLib/src/sysUtils/psTrace.c (revision 1717) +++ /trunk/psLib/src/sysUtils/psTrace.c (revision 1718) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-08 00:09:06 $ + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-08 05:59:41 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -27,4 +27,12 @@ not the same thing as a node's "level", which corresponds to the trace level of that node. + +I think the following is the correct behavior, but not sure: + PS_UNKNOWN_TRACE_LEVEL: We never set the level of a component to this + value. This value is only used when psTraceGetLevel is called with + a bad component name, or if the component root is undefined, I think. + + PS_DEFAULT_TRACE_LEVEL: This should only be used when adding the + intermediate components of a long name. Ie. the "B" in .A.B.C *****************************************************************************/ @@ -40,8 +48,5 @@ #include "psString.h" #include "psError.h" -#include "psLogMsg.h" - -#include "psSysUtilsErrors.h" -#define ERRORNAME_PREFIX PS_ERRORNAME_DOMAIN ".psTrace." +#include "psAbort.h" static p_psComponent* cRoot = NULL; // The root of the trace component @@ -98,4 +103,6 @@ /***************************************************************************** Set all trace levels to zero. + +XXX: Currently, no function calls this routine. *****************************************************************************/ void p_psTraceReset(p_psComponent* currentNode) @@ -110,7 +117,6 @@ for (i = 0; i < currentNode->n; i++) { if (NULL == currentNode->subcomp[i]) { - psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, - PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, - i, currentNode->name); + psError(__func__, + "Sub-component %d of node %s in the trace tree is NULL.\n", i, currentNode->name); } else { p_psTraceReset(currentNode->subcomp[i]); @@ -125,5 +131,5 @@ void psTraceReset() { - psFree(cRoot); + componentFree(cRoot); cRoot = NULL; } @@ -135,5 +141,5 @@ to ANSI-C. *****************************************************************************/ -static bool componentAdd(const char *addNodeName, int level) +static void componentAdd(const char *addNodeName, int level) { int i = 0; // Loop index variable. @@ -146,7 +152,5 @@ // XXX: Verify that this is the correct behavior. if (strcmp("", addNodeName) == 0) { - psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_NULL,true, - PS_ERRORTEXT_psTrace_ADD_NULL_COMPONENT); - return false; + psAbort(__func__, "Failed to add null component to trace tree.\n"); } @@ -154,12 +158,10 @@ if (strcmp(".", addNodeName) == 0) { cRoot->level = level; - return true; + return; } if (addNodeName[0] != '.') { - psErrorMsg(ERRORNAME_PREFIX "componentAdd",PS_ERR_BAD_PARAMETER_VALUE,true, - PS_ERRORTEXT_psTrace_MALFORMED_COMPONENT_NAME, - addNodeName); - return false; + printf("ERROR: failed to add %s to the root component tree.\n", addNodeName); + exit(1); } @@ -186,16 +188,21 @@ currentNode->subcomp = psRealloc(currentNode->subcomp, (currentNode->n + 1) * sizeof(p_psComponent* )); - currentNode->n++; + currentNode->n = (currentNode->n) + 1; if (pname == NULL) { - currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, level); + // This is the final component to add. + currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, level); } else { - currentNode->subcomp[currentNode->n - 1] = componentAlloc(firstComponent, currentNode->level); - } - currentNode = currentNode->subcomp[currentNode->n - 1]; - } - } - - return true; + // We are adding an intermediate component. The trace level + // is not defined. An undefined trace level inherits the + // trace level of it's parent. However, we do not set that + // specifically here since that would inheritance to be a + // static, one-time, type of behavior. + + currentNode->subcomp[(currentNode->n) - 1] = componentAlloc(firstComponent, PS_DEFAULT_TRACE_LEVEL); + } + currentNode = currentNode->subcomp[(currentNode->n) - 1]; + } + } } @@ -211,6 +218,6 @@ zero *****************************************************************************/ -bool psTraceSetLevel(const char *comp, // component of interest - int level) // desired trace level +int psTraceSetLevel(const char *comp, // component of interest + int level) // desired trace level { // If the root component tree does not exist, then initialize it. @@ -219,11 +226,8 @@ } // Add the new component to the component tree. - if ( !componentAdd(comp, level) ) { - psErrorMsg(ERRORNAME_PREFIX "psTraceSetLevel", PS_ERR_UNKNOWN, false, - PS_ERRORTEXT_psTrace_FAILED_TO_ADD_COMPONENT,level,comp); - return false; - } + componentAdd(comp, level); + // return 0 on success. - return true; + return 0; } @@ -249,4 +253,5 @@ p_psComponent* currentNode = cRoot; int i = 0; + int defaultLevel = 0; if (NULL == currentNode) { @@ -262,4 +267,5 @@ } + defaultLevel = cRoot->level; strcpy(name, aname); pname = &name[1]; @@ -268,22 +274,31 @@ for (i = 0; i < currentNode->n; i++) { if (NULL == currentNode->subcomp[i]) { - psLogMsg(ERRORNAME_PREFIX "p_psTraceReset", PS_LOG_WARN, - PS_ERRORTEXT_psTrace_NULL_SUBCOMPONENT, - i, currentNode->name); + psError(__func__, + "Sub-component %d of node %s in trace tree is NULL.\n", i, currentNode->name); } if (strcmp(currentNode->subcomp[i]->name, firstComponent) == 0) { currentNode = currentNode->subcomp[i]; + // For level inheritance purpose, we save the level of this + // component if it is not DEFAULT. + if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { + defaultLevel = currentNode->level; + } + // Determine if thisis the last component: if (pname == NULL) { - return (currentNode->level); + if (currentNode->level != PS_DEFAULT_TRACE_LEVEL) { + return (currentNode->level); + } else { + return(PS_DEFAULT_TRACE_LEVEL); + } } } } } - return (PS_UNKNOWN_TRACE_LEVEL); -} - -/***************************************************************************** - psGetTraceLevel() + return(defaultLevel); +} + +/***************************************************************************** + psTraceLevelGet() Return a trace level of "name" in the root component tree. If the exact string of components in "name" does not exist in the root @@ -303,4 +318,37 @@ // Search the component root tree, determine the trace level. return (doGetTraceLevel(name)); +} + +/***************************************************************************** + doPrintTraceLevelsOld() + This function recursively searches the component tree supplied by the + parameter "comp" and prints the name and level of each component. + Inputs: + comp: a node in the component tree. + level: the level of that node + Outputs: + none + Returns: + null + *****************************************************************************/ +static void doPrintTraceLevelsOld(const p_psComponent* comp, int depth) +{ + int i = 0; + + if (comp->name[0] == '\0') { + printf("%*s%-*s %d\n", depth, "", 20 - depth, + "(root)", (comp->level == PS_DEFAULT_TRACE_LEVEL) ? 0 : comp->level); + } else { + if (comp->level == PS_DEFAULT_TRACE_LEVEL) { + printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); + } else { + printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, + comp->level); + } + } + + for (i = 0; i < comp->n; i++) { + doPrintTraceLevelsOld(comp->subcomp[i], depth + 1); + } } @@ -317,23 +365,31 @@ null *****************************************************************************/ -static void doPrintTraceLevels(const p_psComponent* comp, int depth) +static void doPrintTraceLevels(const p_psComponent* comp, + int depth, + int defLevel) { int i = 0; if (comp->name[0] == '\0') { - printf("%*s%-*s %d\n", depth, "", 20 - depth, - "(root)", (comp->level == PS_UNKNOWN_TRACE_LEVEL) ? 0 : comp->level); + psAbort(__func__, "component name is NULL\n"); } else { - if (comp->level == PS_UNKNOWN_TRACE_LEVEL) { - printf("%*s%-*s %s\n", depth, "", 20 - depth, comp->name, "."); + if (comp->level == PS_DEFAULT_TRACE_LEVEL) { + printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, + defLevel); } else { - printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, comp->level); + printf("%*s%-*s %d\n", depth, "", 20 - depth, comp->name, + comp->level); } } for (i = 0; i < comp->n; i++) { - doPrintTraceLevels(comp->subcomp[i], depth + 1); - } -} + if (comp->level == PS_DEFAULT_TRACE_LEVEL) { + doPrintTraceLevels(comp->subcomp[i], depth + 1, defLevel); + } else { + doPrintTraceLevels(comp->subcomp[i], depth + 1, comp->level); + } + } +} + /***************************************************************************** @@ -353,5 +409,5 @@ } - doPrintTraceLevels(cRoot, 0); + doPrintTraceLevels(cRoot, 0, PS_DEFAULT_TRACE_LEVEL); } @@ -379,8 +435,5 @@ if (NULL == comp) { - psErrorMsg(ERRORNAME_PREFIX "psTrace", PS_ERR_BAD_PARAMETER_NULL, true, - PS_ERRORTEXT_psTrace_NULL_TRACETREE, - __func__); - return; + psError(__func__, "p_psTrace() called on a NULL trace level tree\n"); } // Only display this message if it's trace level is less than the level