Index: unk/psLib/src/math/psFunctions.c =================================================================== --- /trunk/psLib/src/math/psFunctions.c (revision 4963) +++ (revision ) @@ -1,1532 +1,0 @@ -/** @file psFunctions.c -* -* @brief Contains basic function allocation, deallocation, and evaluation -* routines. -* -* This file will hold the functions for allocated, freeing, and evaluating -* polynomials. It also contains a Gaussian functions. -* -* @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-07 21:35:50 $ -* -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii -* -* XXX: Should the "coeffErr[]" be used as well? Bug ???. Ignore coeffErr -* -* XXX: In the various polyAlloc(n) functions, n is really the order of the -* polynomial plus 1. To create a 2nd-order polynomial, n == 3. -*/ -/*****************************************************************************/ -/* INCLUDE FILES */ -/*****************************************************************************/ -#include -#include -#include -#include - -#include "psRandom.h" -#include "psMemory.h" -#include "psVector.h" -#include "psScalar.h" -#include "psTrace.h" -#include "psError.h" -#include "psLogMsg.h" -#include "psFunctions.h" -#include "psConstants.h" -#include "psErrorText.h" - -/*****************************************************************************/ -/* DEFINE STATEMENTS */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* TYPE DEFINITIONS */ -/*****************************************************************************/ -static void polynomial1DFree(psPolynomial1D* poly); -static void polynomial2DFree(psPolynomial2D* poly); -static void polynomial3DFree(psPolynomial3D* poly); -static void polynomial4DFree(psPolynomial4D* poly); -static void spline1DFree(psSpline1D *tmpSpline); -static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); -static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); - -/*****************************************************************************/ -/* GLOBAL VARIABLES */ -/*****************************************************************************/ - -// None - -/*****************************************************************************/ -/* FILE STATIC VARIABLES */ -/*****************************************************************************/ - -// None - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - LOCAL */ -/*****************************************************************************/ - -static void spline1DFree(psSpline1D *tmpSpline) -{ - psS32 i; - - if (tmpSpline == NULL) { - return; - } - - if (tmpSpline->spline != NULL) { - for (i=0;in;i++) { - psFree((tmpSpline->spline)[i]); - } - psFree(tmpSpline->spline); - } - - if (tmpSpline->p_psDeriv2 != NULL) { - psFree(tmpSpline->p_psDeriv2); - } - psFree(tmpSpline->knots); - - return; -} - -static void polynomial1DFree(psPolynomial1D* poly) -{ - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void polynomial2DFree(psPolynomial2D* poly) -{ - unsigned int x = 0; - - for (x = 0; x < poly->nX; x++) { - psFree(poly->coeff[x]); - psFree(poly->coeffErr[x]); - psFree(poly->mask[x]); - } - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void polynomial3DFree(psPolynomial3D* poly) -{ - unsigned int x = 0; - unsigned int y = 0; - - for (x = 0; x < poly->nX; x++) { - for (y = 0; y < poly->nY; y++) { - psFree(poly->coeff[x][y]); - psFree(poly->coeffErr[x][y]); - psFree(poly->mask[x][y]); - } - psFree(poly->coeff[x]); - psFree(poly->coeffErr[x]); - psFree(poly->mask[x]); - } - - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void polynomial4DFree(psPolynomial4D* poly) -{ - unsigned int x = 0; - unsigned int y = 0; - unsigned int z = 0; - - for (x = 0; x < poly->nX; x++) { - for (y = 0; y < poly->nY; y++) { - for (z = 0; z < poly->nZ; z++) { - psFree(poly->coeff[x][y][z]); - psFree(poly->coeffErr[x][y][z]); - psFree(poly->mask[x][y][z]); - } - psFree(poly->coeff[x][y]); - psFree(poly->coeffErr[x][y]); - psFree(poly->mask[x][y]); - } - psFree(poly->coeff[x]); - psFree(poly->coeffErr[x]); - psFree(poly->mask[x]); - } - - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -/***************************************************************************** -createChebyshevPolys(n): this routine takes as input the required order n, -and returns as output as a pointer to an array of n psPolynomial1D -structures, corresponding to the first n Chebyshev polynomials. - -XXX: The output should be static since the Chebyshev polynomials might be -used frequently and the data structure created here does not contain the -outer coefficients of the Chebyshev polynomials. - *****************************************************************************/ -static psPolynomial1D **createChebyshevPolys(psS32 maxChebyPoly) -{ - PS_ASSERT_INT_NONNEGATIVE(maxChebyPoly, NULL); - - psPolynomial1D **chebPolys = NULL; - - chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *)); - for (psS32 i = 0; i < maxChebyPoly; i++) { - chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD); - } - - // Create the Chebyshev polynomials. - // Polynomial i has i-th order. - chebPolys[0]->coeff[0] = 1; - - // XXX: Bug 296 - if (maxChebyPoly > 1) { - chebPolys[1]->coeff[1] = 1; - - for (psS32 i = 2; i < maxChebyPoly; i++) { - for (psS32 j = 0; j < chebPolys[i - 1]->n; j++) { - chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j]; - } - for (psS32 j = 0; j < chebPolys[i - 2]->n; j++) { - chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; - } - } - } else { - // XXX: Code this. - printf("WARNING: %d-order chebyshev polynomials not correctly implemented.\n", maxChebyPoly); - } - - return (chebPolys); -} - -/***************************************************************************** - Polynomial coefficients will be accessed in [w][x][y][z] fashion. - *****************************************************************************/ -static psF64 ordPolynomial1DEval(psF64 x, const psPolynomial1D* poly) -{ - psS32 loop_x = 0; - psF64 polySum = 0.0; - psF64 xSum = 1.0; - - psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, - "---- Calling ordPolynomial1DEval(%lf)\n", x); - psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, - "Polynomial order is %d\n", poly->n); - for (loop_x = 0; loop_x < poly->n; loop_x++) { - psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, - "Polynomial coeff[%d] is %lf\n", loop_x, poly->coeff[loop_x]); - } - - for (loop_x = 0; loop_x < poly->n; loop_x++) { - if (poly->mask[loop_x] == 0) { - psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 10, - "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]); - polySum += xSum * poly->coeff[loop_x]; - } - xSum *= x; - } - - return(polySum); -} - -// XXX: You can do this without having to psAlloc() vector d. -// XXX: How does the mask vector effect Crenshaw's formula? -// XXX: We assume that x is scaled between -1.0 and 1.0; -static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly) -{ - PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - // XXX: Create a macro for this in psConstants.h - if (poly->n < 1) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: Chebyshev polynomial as order %d.", poly->n); - return(NAN); - } - psVector *d; - psS32 n = poly->n; - psS32 i; - psF64 tmp = 0.0; - - // Special case where the Chebyshev poly is constant. - if (n == 1) { - if (poly->mask[0] == 0) { - tmp += poly->coeff[0]; - } - return(tmp); - } - - // Special case where the Chebyshev poly is linear. - if (n == 2) { - if (poly->mask[0] == 0) { - tmp+= poly->coeff[0]; - } - if (poly->mask[1] == 0) { - tmp+= poly->coeff[1] * x; - } - return(tmp); - } - - // General case where the Chebyshev poly has 2 or more terms. - d = psVectorAlloc(n, PS_TYPE_F64); - if(poly->mask[n-1] == 0) { - d->data.F64[n-1] = poly->coeff[n-1]; - } else { - d->data.F64[n-1] = 0.0; - } - - d->data.F64[n-2] = (2.0 * x * d->data.F64[n-1]); - if(poly->mask[n-2] == 0) { - d->data.F64[n-2] += poly->coeff[n-2]; - } - - for (i=n-3;i>=1;i--) { - d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - - (d->data.F64[i+2]); - if(poly->mask[i] == 0) { - d->data.F64[i] += poly->coeff[i]; - } - } - - tmp = (x * d->data.F64[1]) - - (d->data.F64[2]); - if(poly->mask[0] == 0) { - tmp += (0.5 * poly->coeff[0]); - } - psFree(d); - return(tmp); - - /* This is old code that does not use Clenshaw's formula. Get rid of it. - - psS32 n; - psS32 i; - psF32 tmp; - psPolynomial1D **chebPolys = NULL; - - n = poly->n; - chebPolys = createChebyshevPolys(n); - - tmp = 0.0; - for (i=0;in;i++) { - tmp+= (poly->coeff[i] * psPolynomial1DEval(x, chebPolys[i])); - } - tmp-= (poly->coeff[0]/2.0); - - - return(tmp); - */ -} - -static psF64 ordPolynomial2DEval(psF64 x, - psF64 y, - const psPolynomial2D* poly) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - psS32 loop_x = 0; - psS32 loop_y = 0; - psF64 polySum = 0.0; - psF64 xSum = 1.0; - psF64 ySum = 1.0; - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - if (poly->mask[loop_x][loop_y] == 0) { - polySum += ySum * poly->coeff[loop_x][loop_y]; - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -static psF64 chebPolynomial2DEval(psF64 x, psF64 y, const psPolynomial2D* poly) -{ - PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 i = 0; - psF64 polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - chebPolys = createChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - if (poly->mask[loop_x][loop_y] == 0) { - polySum += poly->coeff[loop_x][loop_y] * - psPolynomial1DEval(chebPolys[loop_x], x) * - psPolynomial1DEval(chebPolys[loop_y], y); - } - } - } - for (i=0;inX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - if (poly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += zSum * poly->coeff[loop_x][loop_y][loop_z]; - } - zSum *= z; - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -static psF64 chebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psPolynomial3D* poly) -{ - PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - psF64 polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - if (poly->nZ > maxChebyPoly) { - maxChebyPoly = poly->nZ; - } - chebPolys = createChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - if (poly->mask[loop_x][loop_y][loop_z] == 0) { - polySum += poly->coeff[loop_x][loop_y][loop_z] * - psPolynomial1DEval(chebPolys[loop_x], x) * - psPolynomial1DEval(chebPolys[loop_y], y) * - psPolynomial1DEval(chebPolys[loop_z], z); - } - } - } - } - - for (i=0;inX; loop_x++) { - ySum = xSum; - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - zSum = ySum; - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - tSum = zSum; - for (loop_t = 0; loop_t < poly->nT; loop_t++) { - if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { - polySum += tSum * poly->coeff[loop_x][loop_y][loop_z][loop_t]; - } - tSum *= t; - } - zSum *= z; - } - ySum *= y; - } - xSum *= x; - } - - return(polySum); -} - -static psF64 chebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psPolynomial4D* poly) -{ - PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_DOUBLE_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - PS_ASSERT_DOUBLE_WITHIN_RANGE(t, -1.0, 1.0, 0.0); - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 loop_t = 0; - psS32 i = 0; - psF64 polySum = 0.0; - psPolynomial1D* *chebPolys = NULL; - psS32 maxChebyPoly = 0; - - // Determine how many Chebyshev polynomials - // are needed, then create them. - maxChebyPoly = poly->nX; - if (poly->nY > maxChebyPoly) { - maxChebyPoly = poly->nY; - } - if (poly->nZ > maxChebyPoly) { - maxChebyPoly = poly->nZ; - } - if (poly->nT > maxChebyPoly) { - maxChebyPoly = poly->nT; - } - chebPolys = createChebyshevPolys(maxChebyPoly); - - for (loop_x = 0; loop_x < poly->nX; loop_x++) { - for (loop_y = 0; loop_y < poly->nY; loop_y++) { - for (loop_z = 0; loop_z < poly->nZ; loop_z++) { - for (loop_t = 0; loop_t < poly->nT; loop_t++) { - if (poly->mask[loop_x][loop_y][loop_z][loop_t] == 0) { - polySum += poly->coeff[loop_x][loop_y][loop_z][loop_t] * - psPolynomial1DEval(chebPolys[loop_x], x) * - psPolynomial1DEval(chebPolys[loop_y], y) * - psPolynomial1DEval(chebPolys[loop_z], z) * - psPolynomial1DEval(chebPolys[loop_t], t); - } - } - } - } - } - - for (i=0;idata.TYPE, true); \ - \ - psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 4, \ - "---- fullInterpolate1D##TYPE() begin (%d-order at x=%f) (%d data points)----\n", n-1, x, n); \ - \ - for (i=0;idata.TYPE[i] = range[i]; \ - psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 6, \ - "p->data.TYPE[%d] is %f\n", i, p->data.TYPE[i]); \ - \ - } \ - \ - /* From NR, during each iteration of the m loop, we are computing the \ - p_{i ... i+m} terms. \ - */ \ - for (m=1;mdata.TYPE[i] = (((x-domain[i+m]) * p->data.TYPE[i]) + \ - ((domain[i]-x) * p->data.TYPE[i+1])) / \ - (domain[i] - domain[i+m]); \ - /*printf("((%f-%f * %f) + (%f-%f * %f)) / (%f - %f)\n", x, domain[i+m], p->data.TYPE[i], domain[i], x, p->data.TYPE[i+1], domain[i], domain[i+m]); \ - */ \ - psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 6, \ - "p->data.TYPE[%d] is %f\n", i, p->data.TYPE[i]); \ - } \ - } \ - psTrace(".psLib.dataManip.psFunctions.fullInterpolate1D##TYPE", 4, \ - "---- fullInterpolate1D##TYPE() end ----\n"); \ - \ - return(p->data.TYPE[0]); \ -} \ - -/* -FUNC_MACRO_FULL_INTERPOLATE_1D(U8) -FUNC_MACRO_FULL_INTERPOLATE_1D(U16) -FUNC_MACRO_FULL_INTERPOLATE_1D(U32) -FUNC_MACRO_FULL_INTERPOLATE_1D(U64) -FUNC_MACRO_FULL_INTERPOLATE_1D(S8) -FUNC_MACRO_FULL_INTERPOLATE_1D(S16) -FUNC_MACRO_FULL_INTERPOLATE_1D(S32) -FUNC_MACRO_FULL_INTERPOLATE_1D(S64) -FUNC_MACRO_FULL_INTERPOLATE_1D(F64) -*/ -FUNC_MACRO_FULL_INTERPOLATE_1D(F32) - - -/***************************************************************************** -interpolate1DF32(): this is the base 1-D flat memory routine to perform -LaGrange interpolation. - *****************************************************************************/ -static psF32 interpolate1DF32(psF32 *domain, - psF32 *range, - psS32 n, - psS32 order, - psF32 x) -{ - PS_ASSERT_PTR_NON_NULL(domain, NAN) - PS_ASSERT_PTR_NON_NULL(range, NAN) - // XXX: Check valid values for n, order, and x? - - psS32 binNum; - psS32 numIntPoints = order+1; - psS32 origin; - - psTrace(".psLib.dataManip.psFunctions.interpolate1DF32", 4, - "---- interpolate1DF32() begin ----\n"); - - binNum = vectorBinDisectF32(domain, n, x); - - if (0 == numIntPoints%2) { - origin = binNum - ((numIntPoints/2) - 1); - } else { - origin = binNum - (numIntPoints/2); - if ((x-domain[binNum]) > (domain[binNum+1]-x)) { - // x is closer to binNum+1. - origin = 1 + (binNum - (numIntPoints/2)); - } - } - if (origin < 0) { - origin = 0; - } - if ((origin + numIntPoints) > n) { - origin = n - numIntPoints; - } - - psTrace(".psLib.dataManip.psFunctions.interpolate1DF32", 4, - "---- interpolate1DF32() end ----\n"); - return(fullInterpolate1DF32(&domain[origin], &range[origin], order+1, x)); -} - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - PUBLIC */ -/*****************************************************************************/ - -/***************************************************************************** - Evaluate a non-normalized Gaussian with the given mean and sigma at the - given coordianate. Note that this is not a Gaussian deviate. The - evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] - *****************************************************************************/ -float psGaussian(float x, float mean, float sigma, bool normal) -{ - psF32 tmp = 1.0; - - psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, - "---- psGaussian() begin ----\n"); - - if (normal == true) { - tmp = 1.0 / sqrtf(2.0 * M_PI * (sigma * sigma)); - } - - psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, - "---- psGaussian() end ----\n"); - return(tmp * exp(-((x - mean) * (x - mean)) / (2.0 * sigma * sigma))); -} - -/***************************************************************************** - p_psGaussianDev() - This private routine (formerly a psLib API routine) creates a psVector of the - specified size and type F32 and fills it with a random Gaussian distribution - of numbers with the specified mean and sigma. - -XXX: It's possible to have a different seed everutime. However, for now, -for testability, we use a common seed. - *****************************************************************************/ -#define PS_XXX_GAUSSIAN_SEED 1995 -psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) -{ - PS_ASSERT_INT_NONNEGATIVE(Npts, NULL); - - // psRandom *r = psRandomAlloc(PS_RANDOM_TAUS, p_psRandomGetSystemSeed()); - psRandom *r = psRandomAlloc(PS_RANDOM_TAUS, PS_XXX_GAUSSIAN_SEED); - psVector* gauss = psVectorAlloc(Npts, PS_TYPE_F32); - for (psS32 i = 0; i < Npts; i++) { - gauss->data.F32[i] = mean + p_psRandomGaussian(r, sigma); - } - psFree(r); - - return(gauss); -} - -/***************************************************************************** - This routine must allocate memory for the polynomial structures. - *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(unsigned int n, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(n, NULL); - - unsigned int i = 0; - psPolynomial1D* newPoly = NULL; - - newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); - psMemSetDeallocator(newPoly, (psFreeFunc) polynomial1DFree); - - newPoly->type = type; - newPoly->n = n; - newPoly->coeff = psAlloc(n * sizeof(psF64)); - newPoly->coeffErr = psAlloc(n * sizeof(psF64)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); - for (i = 0; i < n; i++) { - newPoly->coeff[i] = 0.0; - newPoly->coeffErr[i] = 0.0; - newPoly->mask[i] = 0; - } - - return(newPoly); -} - -psPolynomial2D* psPolynomial2DAlloc( unsigned int nX, - unsigned int nY, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(nX, NULL); - PS_ASSERT_INT_POSITIVE(nY, NULL); - - unsigned int x = 0; - unsigned int y = 0; - psPolynomial2D* newPoly = NULL; - - newPoly = (psPolynomial2D* ) psAlloc(sizeof(psPolynomial2D)); - psMemSetDeallocator(newPoly, (psFreeFunc) polynomial2DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - - newPoly->coeff = psAlloc(nX * sizeof(psF64 *)); - newPoly->coeffErr = psAlloc(nX * sizeof(psF64 *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = psAlloc(nY * sizeof(psF64)); - newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = 0.0; - newPoly->coeffErr[x][y] = 0.0; - newPoly->mask[x][y] = 0; - } - } - - return(newPoly); -} - -psPolynomial3D* psPolynomial3DAlloc( unsigned int nX, - unsigned int nY, - unsigned int nZ, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(nX, NULL); - PS_ASSERT_INT_POSITIVE(nY, NULL); - PS_ASSERT_INT_POSITIVE(nZ, NULL); - - unsigned int x = 0; - unsigned int y = 0; - unsigned int z = 0; - psPolynomial3D* newPoly = NULL; - - newPoly = (psPolynomial3D* ) psAlloc(sizeof(psPolynomial3D)); - psMemSetDeallocator(newPoly, (psFreeFunc) polynomial3DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - - newPoly->coeff = psAlloc(nX * sizeof(psF64 **)); - newPoly->coeffErr = psAlloc(nX * sizeof(psF64 **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 *)); - newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64)); - newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); - } - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = 0.0; - newPoly->coeffErr[x][y][z] = 0.0; - newPoly->mask[x][y][z] = 0; - } - } - } - - return(newPoly); -} - -psPolynomial4D* psPolynomial4DAlloc( unsigned int nX, - unsigned int nY, - unsigned int nZ, - unsigned int nT, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(nX, NULL); - PS_ASSERT_INT_POSITIVE(nY, NULL); - PS_ASSERT_INT_POSITIVE(nZ, NULL); - PS_ASSERT_INT_POSITIVE(nT, NULL); - - unsigned int x = 0; - unsigned int y = 0; - unsigned int z = 0; - unsigned int t = 0; - psPolynomial4D* newPoly = NULL; - - newPoly = (psPolynomial4D* ) psAlloc(sizeof(psPolynomial4D)); - psMemSetDeallocator(newPoly, (psFreeFunc) polynomial4DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - newPoly->nT = nT; - - newPoly->coeff = psAlloc(nX * sizeof(psF64 ***)); - newPoly->coeffErr = psAlloc(nX * sizeof(psF64 ***)); - newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = psAlloc(nY * sizeof(psF64 **)); - newPoly->coeffErr[x] = psAlloc(nY * sizeof(psF64 **)); - newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = psAlloc(nZ * sizeof(psF64 *)); - newPoly->coeffErr[x][y] = psAlloc(nZ * sizeof(psF64 *)); - newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *)); - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = psAlloc(nT * sizeof(psF64)); - newPoly->coeffErr[x][y][z] = psAlloc(nT * sizeof(psF64)); - newPoly->mask[x][y][z] = (char *)psAlloc(nT * sizeof(char)); - } - } - } - for (x = 0; x < nX; x++) { - for (y = 0; y < nY; y++) { - for (z = 0; z < nZ; z++) { - for (t = 0; t < nT; t++) { - newPoly->coeff[x][y][z][t] = 0.0; - newPoly->coeffErr[x][y][z][t] = 0.0; - newPoly->mask[x][y][z][t] = 0; - } - } - } - } - - return(newPoly); -} - - -bool psMemCheckPolynomial1D(psPtr ptr) -{ - return ( psMemGetDeallocator(ptr) == (psFreeFunc)polynomial1DFree ); -} - -bool psMemCheckPolynomial2D(psPtr ptr) -{ - return ( psMemGetDeallocator(ptr) == (psFreeFunc)polynomial2DFree ); -} - -bool psMemCheckPolynomial3D(psPtr ptr) -{ - return ( psMemGetDeallocator(ptr) == (psFreeFunc)polynomial3DFree ); -} - -bool psMemCheckPolynomial4D(psPtr ptr) -{ - return ( psMemGetDeallocator(ptr) == (psFreeFunc)polynomial4DFree ); -} - -bool psMemCheckSpline1D(psPtr ptr) -{ - return ( psMemGetDeallocator(ptr) == (psFreeFunc)spline1DFree ); -} - - -psF64 psPolynomial1DEval(const psPolynomial1D* poly, - psF64 x) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(ordPolynomial1DEval(x, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(chebPolynomial1DEval(x, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psPolynomial1DEvalVector(const psPolynomial1D *poly, - const psVector *x) -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - - psVector *tmp; - - tmp = psVectorAlloc(x->n, PS_TYPE_F64); - for (psS32 i=0;in;i++) { - tmp->data.F64[i] = psPolynomial1DEval(poly, x->data.F64[i]); - } - - return(tmp); -} - -psF64 psPolynomial2DEval(const psPolynomial2D* poly, - psF64 x, - psF64 y) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(ordPolynomial2DEval(x, y, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(chebPolynomial2DEval(x, y, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psPolynomial2DEvalVector(const psPolynomial2D *poly, - const psVector *x, - const psVector *y) - -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine the length of the output vector to by the minimum of the x,y vectors - if (y->n < vecLen) { - vecLen = y->n; - } - - // Create output vector to return - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate the polynomial at the specified points - for (psS32 i=0; idata.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -psF64 psPolynomial3DEval(const psPolynomial3D* poly, - psF64 x, - psF64 y, - psF64 z) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(ordPolynomial3DEval(x, y, z, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(chebPolynomial3DEval(x, y, z, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psPolynomial3DEvalVector(const psPolynomial3D *poly, - const psVector *x, - const psVector *y, - const psVector *z) - -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine the length of output vector from min of the input vectors - if (y->n < vecLen) { - vecLen = y->n; - } - if (z->n < vecLen) { - vecLen = z->n; - } - - // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psPolynomial3DEval(poly, - x->data.F64[i], - y->data.F64[i], - z->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -psF64 psPolynomial4DEval(const psPolynomial4D* poly, - psF64 x, - psF64 y, - psF64 z, - psF64 t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(ordPolynomial4DEval(x,y,z,t, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(chebPolynomial4DEval(x,y,z,t, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psPolynomial4DEvalVector(const psPolynomial4D *poly, - const psVector *x, - const psVector *y, - const psVector *z, - const psVector *t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine output vector size from min of input vectors - if (z->n < vecLen) { - vecLen = z->n; - } - if (y->n < vecLen) { - vecLen = y->n; - } - if (t->n < vecLen) { - vecLen = t->n; - } - - // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psPolynomial4DEval(poly, - x->data.F64[i], - y->data.F64[i], - z->data.F64[i], - t->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -//typedef struct { -// psS32 n; -// psPolynomial1D **spline; -// psF32 *p_psDeriv2; -// psVector *knots; -//} psSpline1D; - -/***************************************************************************** - NOTE: "n" specifies the number of spline polynomials. Therefore, there - must exist n+1 points in "knots". - -XXX: Is this really needed anymore? - -XXX: Ensure that knots[i+1] != knots[i] - -XXX: What should be the default type for knots be? psF32 is assumed. - *****************************************************************************/ -psSpline1D *psSpline1DAlloc(unsigned int numSplines, - unsigned int order, - float min, - float max) -{ - PS_ASSERT_INT_NONNEGATIVE(numSplines, NULL); - PS_ASSERT_INT_NONNEGATIVE(order, NULL); - PS_ASSERT_FLOAT_NON_EQUAL(max, min, NULL); - - 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 (unsigned int i=0; i < numSplines; i++) { - (tmpSpline->spline)[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); - psF32 width = (max - min) / ((psF32) numSplines); - tmpSpline->knots->data.F32[0] = min; - for (unsigned int i=1; i < numSplines; i++) { - tmpSpline->knots->data.F32[i] = min + (width * (psF32) i); - } - tmpSpline->knots->data.F32[numSplines] = max; - - psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); - return(tmpSpline); -} - -/***************************************************************************** -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, - unsigned int order) -{ - 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 = (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 (unsigned int i=0; i < numSplines; i++) { - (tmpSpline->spline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); - } - - // This will be computed by psVectorFitSpline1D() - tmpSpline->p_psDeriv2 = NULL; - - // - // Ensure that all knots are distinct. - // XXX:Ensure that the knots are monotonic. - // - for (unsigned int i=0; i < bounds->n-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 ([%d] %f %f)\n", i, bounds->data.F32[i], bounds->data.F32[i+1]); - return(NULL); - } - } - tmpSpline->knots = PsVectorDup((psVector *) bounds); - - psMemSetDeallocator(tmpSpline, (psFreeFunc)spline1DFree); - return(tmpSpline); -} - -/***************************************************************************** -vectorBinDisectF32(): This is a macro for a private function which takes as -input a vector an array of data as well as a single value for that data. The -input vector values are assumed to be non-decreasing (v[i-1] <= v[i] for all -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 (-2 or -1). - *****************************************************************************/ -#define FUNC_MACRO_VECTOR_BIN_DISECT(TYPE) \ -static psS32 vectorBinDisect##TYPE(ps##TYPE *bins, \ - psS32 numBins, \ - ps##TYPE x) \ -{ \ - psS32 min; \ - psS32 max; \ - psS32 mid; \ - \ - psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ - "---- Calling vectorBinDisect##TYPE(%f)\n", x); \ - \ - if (x < bins[0]) { \ - psLogMsg(__func__, PS_LOG_WARN, \ - "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ - return(-2); \ - } \ - \ - if (x > bins[numBins-1]) { \ - psLogMsg(__func__, PS_LOG_WARN, \ - "vectorBinDisect%s(): ordinate %f is outside vector range (%f - %f).", \ - #TYPE, x, bins[0], bins[numBins-1]); \ - return(-1); \ - } \ - \ - min = 0; \ - max = numBins-2; \ - mid = ((max+1)-min)/2; \ - \ - while (min != max) { \ - psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ - "(min, mid, max) is (%d, %d, %d): (x, bins) is (%f, %f)\n", \ - min, mid, max, x, bins[mid]); \ - \ - if (x == bins[mid]) { \ - psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ - "---- Exiting vectorBinDisect##TYPE(): bin %d\n", mid); \ - return(mid); \ - } else if (x < bins[mid]) { \ - max = mid-1; \ - } else { \ - min = mid; \ - } \ - mid = ((max+1)+min)/2; \ - } \ - \ - psTrace(".psLib.dataManip.psFunctions.vectorBinDisect##TYPE", 4, \ - "---- Exiting vectorBinDisect##TYPE(): bin %d\n", min); \ - return(min); \ -} \ - -FUNC_MACRO_VECTOR_BIN_DISECT(S8) -FUNC_MACRO_VECTOR_BIN_DISECT(S16) -FUNC_MACRO_VECTOR_BIN_DISECT(S32) -FUNC_MACRO_VECTOR_BIN_DISECT(S64) -FUNC_MACRO_VECTOR_BIN_DISECT(U8) -FUNC_MACRO_VECTOR_BIN_DISECT(U16) -FUNC_MACRO_VECTOR_BIN_DISECT(U32) -FUNC_MACRO_VECTOR_BIN_DISECT(U64) -FUNC_MACRO_VECTOR_BIN_DISECT(F32) -FUNC_MACRO_VECTOR_BIN_DISECT(F64) - -/***************************************************************************** -p_psVectorBinDisect(): A wrapper to the above p_psVectorBinDisect(). - -XXX: Assert that the psVector and psScalar have the same type. - *****************************************************************************/ -psS32 p_psVectorBinDisect(psVector *bins, - psScalar *x) -{ - PS_ASSERT_VECTOR_NON_NULL(bins, -4); - PS_ASSERT_VECTOR_NON_EMPTY(bins, -4); - PS_ASSERT_PTR_NON_NULL(x, -6); - PS_ASSERT_PTR_TYPE_EQUAL(x, bins, -3); - char* strType; - - switch (x->type.type) { - case PS_TYPE_U8: - return(vectorBinDisectU8(bins->data.U8, bins->n, x->data.U8)); - case PS_TYPE_U16: - return(vectorBinDisectU16(bins->data.U16, bins->n, x->data.U16)); - case PS_TYPE_U32: - return(vectorBinDisectU32(bins->data.U32, bins->n, x->data.U32)); - case PS_TYPE_U64: - return(vectorBinDisectU64(bins->data.U64, bins->n, x->data.U64)); - case PS_TYPE_S8: - return(vectorBinDisectS8(bins->data.S8, bins->n, x->data.S8)); - case PS_TYPE_S16: - return(vectorBinDisectS16(bins->data.S16, bins->n, x->data.S16)); - case PS_TYPE_S32: - return(vectorBinDisectS32(bins->data.S32, bins->n, x->data.S32)); - case PS_TYPE_S64: - return(vectorBinDisectS64(bins->data.S64, bins->n, x->data.S64)); - case PS_TYPE_F32: - return(vectorBinDisectF32(bins->data.F32, bins->n, x->data.F32)); - case PS_TYPE_F64: - return(vectorBinDisectF64(bins->data.F64, bins->n, x->data.F64)); - case PS_TYPE_C32: - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, - strType); - return 0; - case PS_TYPE_C64: - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, - strType); - return 0; - case PS_TYPE_BOOL: - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, - strType); - return 0; - } - return(-3); -} - -/***************************************************************************** -p_psVectorInterpolate(): This routine will take as input psVectors domain and -range, and the x value, assumed to lie with the domain vector. It produces -as output the LaGrange interpolated value of a polynomial of the specified -order around the point x. - -XXX: This stuff does not currently work with a mask. - -XXX: add another psScalar argument for the result. - -XXX: The VectorCopy routines seg fault when I declare range32 as static. - *****************************************************************************/ -psScalar *p_psVectorInterpolate(psVector *domain, - psVector *range, - int order, - psScalar *x) -{ - PS_ASSERT_VECTOR_NON_NULL(domain, NULL); - PS_ASSERT_VECTOR_NON_NULL(range, NULL); - PS_ASSERT_PTR_NON_NULL(x, NULL); - PS_ASSERT_INT_NONNEGATIVE(order, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(domain, range, NULL); - PS_ASSERT_PTR_TYPE_EQUAL(domain, range, NULL); - PS_ASSERT_PTR_TYPE_EQUAL(domain, x, NULL); - - psVector *range32 = NULL; - psVector *domain32 = NULL; - psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, - "---- p_psVectorInterpolate() begin ----\n"); - - if (order > (domain->n - 1)) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - PS_ERRORTEXT_psFunctions_NOT_ENOUGH_DATAPOINTS, - order); - return(NULL); - } - - if (x->type.type == PS_TYPE_F32) { - psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, - "---- p_psVectorInterpolate() end ----\n"); - return(psScalarAlloc(interpolate1DF32(domain->data.F32, - range->data.F32, - domain->n, - order, - x->data.F32), PS_TYPE_F32)); - } else if (x->type.type == PS_TYPE_F64) { - // XXX: use recycled vectors here. - range32 = psVectorCopy(range32, range, PS_TYPE_F32); - domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); - - psScalar *tmpScalar = psScalarAlloc((psF64) - interpolate1DF32(domain32->data.F32, - range32->data.F32, - domain32->n, - order, - (psF32) x->data.F64), PS_TYPE_F64); - psFree(range32); - psFree(domain32); - - psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, - "---- p_psVectorInterpolate() end ----\n"); - // XXX: Convert data type to F64? - return(tmpScalar); - - } else { - char* strType; - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, - strType); - } - - psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, - "return(NULL)\n"); - psTrace(".psLib.dataManip.psFunctions.p_psVectorInterpolate", 4, - "---- p_psVectorInterpolate() end ----\n"); - - return(NULL); -} - - -/***************************************************************************** -psSpline1DEval(): this routine takes an existing spline of arbitrary order -and an independent x value. Each determines which spline that x corresponds -to by doing a bracket disection on the knots of the spline data structure -(vectorBinDisectF32()). Then it evaluates the spline at that x location -by a call to the 1D polynomial functions. - -XXX: The spline eval functions require input and output to be F32. however - the spline fit functions require F32 and F64. - -XXX: This only works if spline0>knots if psF32. Must add support for psU32 and -psF64. - *****************************************************************************/ -float psSpline1DEval( - const psSpline1D *spline, - float x -) -{ - PS_ASSERT_PTR_NON_NULL(spline, NAN); - PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN); - PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NAN); - - unsigned int binNum; - unsigned int n; - - n = spline->n; - //XXX binNum = vectorBinDisectF32(spline->domains, (spline->n)+1, x); - binNum = vectorBinDisectF32(spline->knots->data.F32, (spline->n)+1, x); - if (binNum < 0) { - psLogMsg(__func__, PS_LOG_WARN, - "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f).", - x, spline->knots->data.F32[0], - spline->knots->data.F32[n-1]); - - if (x < spline->knots->data.F32[0]) { - return(psPolynomial1DEval(spline->spline[0], - x)); - } else if (x > spline->knots->data.F32[n-1]) { - return(psPolynomial1DEval(spline->spline[n-1], - x)); - } - } - - return(psPolynomial1DEval(spline->spline[binNum], - x)); -} - -// XXX: The spline eval functions require input and output to be F32. -// however the spline fit functions require F32 and F64. -psVector *psSpline1DEvalVector( - const psSpline1D *spline, - const psVector *x -) -{ - PS_ASSERT_PTR_NON_NULL(spline, NULL); - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); - PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL); - - unsigned int i; - psVector *tmpVector; - - tmpVector = psVectorAlloc(x->n, PS_TYPE_F32); - if (x->type.type == PS_TYPE_F32) { - for (i=0;in;i++) { - tmpVector->data.F32[i] = psSpline1DEval( - spline, - x->data.F32[i] - ); - } - } else if (x->type.type == PS_TYPE_F64) { - for (i=0;in;i++) { - tmpVector->data.F32[i] = psSpline1DEval( - spline, - (psF32) x->data.F64[i] - ); - } - } else { - char* strType; - PS_TYPE_NAME(strType,x->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, - PS_ERRORTEXT_psFunctions_TYPE_NOT_SUPPORTED, - strType); - return(NULL); - } - - return(tmpVector); -} Index: unk/psLib/src/math/psFunctions.h =================================================================== --- /trunk/psLib/src/math/psFunctions.h (revision 4963) +++ (revision ) @@ -1,385 +1,0 @@ -/** @file psFunctions.h - * @brief Standard Mathematical Functions. - * @ingroup Stats - * - * This file will hold the prototypes for procedures which allocate, free, - * and evaluate various polynomials. Those polynomial structures are also - * defined here. - * - * @ingroup Stats - * - * @author Someone at IfA - * @author GLG, MHPCC - * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-08-31 22:28:35 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ - -#ifndef PS_FUNCTIONS_H -#define PS_FUNCTIONS_H - -#include -#include -#include -#include - -#include "psVector.h" -#include "psScalar.h" - -/** \addtogroup Stats - * \{ - */ - -/** Evaluate a non-normalized Gaussian with the given mean and sigma at the - * given coordianate. - * - * Note that this is not a Gaussian deviate. The evaluated Gaussian is: - * \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] - * - * @return float value on the gaussian curve given the input parameters - */ -float psGaussian( - float x, ///< Value at which to evaluate - float mean, ///< Mean for the Gaussian - float sigma, ///< Standard deviation for the Gaussian - bool normal ///< Indicates whether result should be normalized -); - -/** Produce a vector of random numbers from a Gaussian distribution with - * the specified mean and sigma - * - * @return psVector* vector of random numbers - * - */ -psVector* p_psGaussianDev( - psF32 mean, ///< The mean of the Gaussian - psF32 sigma, ///< The sigma of the Gaussian - psS32 Npts ///< The size of the vector -); - -/** Polynomial Type. - * - * Enumeration for Polynomial types. - */ -typedef enum { - PS_POLYNOMIAL_ORD, ///< Ordinary Polynomial - PS_POLYNOMIAL_CHEB ///< Chebyshev Polynomial -} -psPolynomialType; - -/** One-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int n; ///< Number of terms - psF64 *coeff; ///< Coefficients - psF64 *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask -} -psPolynomial1D; - -/** Two-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in x - unsigned int nY; ///< Number of terms in y - psF64 **coeff; ///< Coefficients - psF64 **coeffErr; ///< Error in coefficients - char **mask; ///< Coefficients mask -} -psPolynomial2D; - -/** Three-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in x - unsigned int nY; ///< Number of terms in y - unsigned int nZ; ///< Number of terms in z - psF64 ***coeff; ///< Coefficients - psF64 ***coeffErr; ///< Error in coefficients - char ***mask; ///< Coefficients mask -} -psPolynomial3D; - -/** Four-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in x - unsigned int nY; ///< Number of terms in y - unsigned int nZ; ///< Number of terms in z - unsigned int nT; ///< Number of terms in t - psF64 ****coeff; ///< Coefficients - psF64 ****coeffErr; ///< Error in coefficients - char ****mask; ///< Coefficients mask -} -psPolynomial4D; - - -/** Allocates a psPolynomial1D structure with n terms - * - * @return psPolynomial1D* new 1-D polynomial struct - */ -psPolynomial1D* psPolynomial1DAlloc( - unsigned int n, ///< Number of terms - psPolynomialType type ///< Polynomial Type -); - -/** Allocates a 2-D polynomial structure - * - * @return psPolynomial2D* new 2-D polynomial struct - */ -psPolynomial2D* psPolynomial2DAlloc( - unsigned int nX, ///< Number of terms in x - unsigned int nY, ///< Number of terms in y - psPolynomialType type ///< Polynomial Type -); - -/** Allocates a 3-D polynomial structure - * - * @return psPolynomial3D* new 3-D polynomial struct - */ -psPolynomial3D* psPolynomial3DAlloc( - unsigned int nX, ///< Number of terms in x - unsigned int nY, ///< Number of terms in y - unsigned int nZ, ///< Number of terms in z - psPolynomialType type ///< Polynomial Type -); - -/** Allocates a 4-D polynomial structure - * - * @return psPolynomial4D* new 4-D polynomial struct - */ -psPolynomial4D* psPolynomial4DAlloc( - unsigned int nX, ///< Number of terms in x - unsigned int nY, ///< Number of terms in y - unsigned int nZ, ///< Number of terms in z - unsigned int nT, ///< Number of terms in t - psPolynomialType type ///< Polynomial Type -); - -/** Checks the type of a particular pointer. - * - * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. - * - * @return bool: True if the pointer matches a psPolynomial1D structure, false otherwise. - */ -bool psMemCheckPolynomial1D( - psPtr ptr ///< the pointer whose type to check -); - -/** Checks the type of a particular pointer. - * - * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. - * - * @return bool: True if the pointer matches a psPolynomial2D structure, false otherwise. - */ -bool psMemCheckPolynomial2D( - psPtr ptr ///< the pointer whose type to check -); - -/** Checks the type of a particular pointer. - * - * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. - * - * @return bool: True if the pointer matches a psPolynomial3D structure, false otherwise. - */ -bool psMemCheckPolynomial3D( - psPtr ptr ///< the pointer whose type to check -); - -/** Checks the type of a particular pointer. - * - * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. - * - * @return bool: True if the pointer matches a psPolynomial4D structure, false otherwise. - */ -bool psMemCheckPolynomial4D( - psPtr ptr ///< the pointer whose type to check -); - -/** Checks the type of a particular pointer. - * - * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. - * - * @return bool: True if the pointer matches a psSpline1D structure, false otherwise. - */ -bool psMemCheckSpline1D( - psPtr ptr ///< the pointer whose type to check -); - -/** Evaluates a 1-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psPolynomial1DEval( - const psPolynomial1D* poly, ///< Coefficients for the polynomial - psF64 x ///< location at which to evaluate -); - -/** Evaluates a 2-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psPolynomial2DEval( - const psPolynomial2D* poly, ///< Coefficients for the polynomial - psF64 x, ///< x location at which to evaluate - psF64 y ///< y location at which to evaluate -); - -/** Evaluates a 3-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psPolynomial3DEval( - const psPolynomial3D* poly, ///< Coefficients for the polynomial - psF64 x, ///< x location at which to evaluate - psF64 y, ///< y location at which to evaluate - psF64 z ///< z location at which to evaluate -); - -/** Evaluates a 4-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psPolynomial4DEval( - const psPolynomial4D* poly, ///< Coefficients for the polynomial - psF64 x, ///< x location at which to evaluate - psF64 y, ///< y location at which to evaluate - psF64 z, ///< z location at which to evaluate - psF64 t ///< t location at which to evaluate -); - -/** Evaluates a 1-D polynomial at specific sets of coordinates - * - * @return psVector* results of polynomials at given locations - */ -psVector *psPolynomial1DEvalVector( - const psPolynomial1D *poly, ///< Coefficients for the polynomial - const psVector *x ///< x locations at which to evaluate -); - -/** Evaluates a 2-D polynomial at specific sets of coordinates - * - * @return psVector* results of polynomial at given locations - */ -psVector *psPolynomial2DEvalVector( - const psPolynomial2D *poly, ///< Coefficients for the polynomial - const psVector *x, ///< x locations at which to evaluate - const psVector *y ///< y locations at which to evaluate -); - -/** Evaluates a 3-D polynomial at specific sets of coordinates - * - * @return psVector* results of polynomial at given locations - */ -psVector *psPolynomial3DEvalVector( - const psPolynomial3D *poly, ///< Coefficients for the polynomial - const psVector *x, ///< x locations at which to evaluate - const psVector *y, ///< y locations at which to evaluate - const psVector *z ///< z locations at which to evaluate -); - -/** Evaluates a 4-D polynomial at specific sets of coordinates - * - * @return psVector* results of polynomial at given locations - */ -psVector *psPolynomial4DEvalVector( - const psPolynomial4D *poly, ///< Coefficients for the polynomial - const psVector *x, ///< x locations at which to evaluate - const psVector *y, ///< y locations at which to evaluate - const psVector *z, ///< z locations at which to evaluate - const psVector *t ///< t locations at which to evaluate -); - -/** One-Dimensional Spline */ -typedef struct -{ - unsigned 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 *p_psDomains; ///< The boundaries between each spline piece. Size is n+1. -} -psSpline1D; - -/** Allocates a psSpline1D structure - * - * Allocator for psSpline1D where the bounds are implicitly specified through specifying - * min and max values along with the number of splines. - * - * @return psSpline1D* new 1-D spline struct - */ -psSpline1D *psSpline1DAlloc( - unsigned int n, ///< Number of spline polynomials - unsigned int order, ///< Order of spline polynomials - float min, ///< Lower boundary value of spline polynomials - float max ///< Upper boundary value of spline polynomials -); - -/** Allocates a psSpline1D structure - * - * Allocator for psSpline1D where the bounds are explicitly specified. - * - * @return psSpline1D* new 1-D spline struct - */ -psSpline1D *psSpline1DAllocGeneric( - const psVector *bounds, ///< Bounds for spline polynomials - unsigned int order ///< Order of spline polynomials -); - -/** Evaluates 1-D spline polynomials at a specific coordinate. - * - * @return float result of spline polynomials evaluated at given location - */ -float psSpline1DEval( - const psSpline1D *spline, ///< Coefficients for spline polynomials - float x ///< location at which to evaluate -); - -/** Evaluates 1-D spline polynomials at a set of specific coordinates. - * - * @return psVector* results of spline polynomials evaluated at given locations - */ -psVector *psSpline1DEvalVector( - const psSpline1D *spline, ///< Coefficients of spline polynomials - const psVector *x ///< locations at which to evaluate -); - -/** Performs a binary disection on a given vector. - * Searches through an array of data for a specified value. - * - * @return psS32 corresponding index number of specified value - */ -psS32 p_psVectorBinDisect( - psVector *bins, ///< Array of non-decreasing values - psScalar *x ///< Target value to find -); - -/** Interpolates a series of data points for evaluation at a specific coordinate. Uses a - * Lagrange interpolation method. - * - * @return psScalar* Lagrange interpolation value at given location - */ -psScalar *p_psVectorInterpolate( - psVector *domain, ///< Domain (x coords) for interpolation - psVector *range, ///< Range (y coords) for interpolation - int order, ///< Order of interpolation function - psScalar *x ///< Location at which to evaluate -); - -#if 0 -psF32 p_psNRSpline1DEval(psSpline1D *spline, - const psVector* x, - const psVector* y, - psF32 X); -#endif // #if 0 - -/** \} */ // End of MathGroup Functions - -#endif // #ifndef PS_FUNCTIONS_H -