Index: trunk/psLib/src/math/Makefile.am =================================================================== --- trunk/psLib/src/math/Makefile.am (revision 4446) +++ trunk/psLib/src/math/Makefile.am (revision 4540) @@ -1,38 +1,27 @@ -#Makefile for dataManip functions of psLib +#Makefile for math functions of psLib # -INCLUDES = \ - -I$(top_srcdir)/src/astro \ - -I$(top_srcdir)/src/db \ - -I$(top_srcdir)/src/fft \ - -I$(top_srcdir)/src/fits \ - -I$(top_srcdir)/src/imageops \ - -I$(top_srcdir)/src/math \ - -I$(top_srcdir)/src/mathtypes \ - -I$(top_srcdir)/src/sys \ - -I$(top_srcdir)/src/types \ - -I$(top_srcdir)/src/xml \ - $(all_includes) +noinst_LTLIBRARIES = libpslibmath.la -noinst_LTLIBRARIES = libpslibdataManip.la +libpslibmath_la_SOURCES = \ + psBinaryOp.c \ + psCompare.c \ + psMatrix.c \ + psMinimize.c \ + psRandom.c \ + psFunctions.c \ + psStats.c \ + psUnaryOp.c -libpslibdataManip_la_SOURCES = psUnaryOp.c psBinaryOp.c psStats.c \ - psFunctions.c psMatrix.c psVectorFFT.c psMinimize.c psRandom.c - -BUILT_SOURCES = psDataManipErrors.h -EXTRA_DIST = psDataManipErrors.dat psDataManipErrors.h dataManip.i - -psDataManipErrors.h: psDataManipErrors.dat - $(top_srcdir)/src/psParseErrorCodes --data=$? $@ +EXTRA_DIST = math.i pslibincludedir = $(includedir) pslibinclude_HEADERS = \ + psBinaryOp.h \ + psCompare.h \ psConstants.h \ - psStats.h \ + psMatrix.h \ + psMinimize.h \ + psRandom.h \ psFunctions.h \ - psMatrix.h \ - psBinaryOp.h \ - psUnaryOp.h \ - psVectorFFT.h \ - psMinimize.h \ - psRandom.h - + psStats.h \ + psUnaryOp.h Index: trunk/psLib/src/math/math.i =================================================================== --- trunk/psLib/src/math/math.i (revision 4540) +++ trunk/psLib/src/math/math.i (revision 4540) @@ -0,0 +1,11 @@ +/* math headers */ + +%include "psBinaryOp.h" +%include "psCompare.h" +%include "psConstants.h" +%include "psMatrix.h" +%include "psMinimize.h" +%include "psFunctions.h" +%include "psRandom.h" +%include "psStats.h" +%include "psUnaryOp.h" Index: trunk/psLib/src/math/psBinaryOp.c =================================================================== --- trunk/psLib/src/math/psBinaryOp.c (revision 4446) +++ trunk/psLib/src/math/psBinaryOp.c (revision 4540) @@ -30,6 +30,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-22 02:05:41 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -50,5 +50,5 @@ #include "psLogMsg.h" #include "psConstants.h" -#include "psDataManipErrors.h" +#include "psErrorText.h" /***************************************************************************** @@ -426,5 +426,5 @@ } -psPtr psBinaryOp(psPtr out, const psPtr in1, const char *op, const psPtr in2) +psMathType* psBinaryOp(psPtr out, const psPtr in1, const char *op, const psPtr in2) { Index: trunk/psLib/src/math/psBinaryOp.h =================================================================== --- trunk/psLib/src/math/psBinaryOp.h (revision 4446) +++ trunk/psLib/src/math/psBinaryOp.h (revision 4540) @@ -30,6 +30,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-22 02:05:41 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -55,5 +55,5 @@ * @return psType* : Pointer to either psImage or psVector. */ -psType* psBinaryOp( +psMathType* psBinaryOp( psPtr out, ///< Output type, either psImage or psVector. const psPtr in1, ///< First input, either psImage or psVector. Index: trunk/psLib/src/math/psFunctions.c =================================================================== --- trunk/psLib/src/math/psFunctions.c (revision 4540) +++ trunk/psLib/src/math/psFunctions.c (revision 4540) @@ -0,0 +1,2201 @@ +/** @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.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ + * + * 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 +#include + +#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 dPolynomial1DFree(psDPolynomial1D* poly); +static void dPolynomial2DFree(psDPolynomial2D* poly); +static void dPolynomial3DFree(psDPolynomial3D* poly); +static void dPolynomial4DFree(psDPolynomial4D* 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); +} + +static void dPolynomial1DFree(psDPolynomial1D* poly) +{ + psFree(poly->coeff); + psFree(poly->coeffErr); + psFree(poly->mask); +} + +static void dPolynomial2DFree(psDPolynomial2D* poly) +{ + for (unsigned int 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 dPolynomial3DFree(psDPolynomial3D* 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 dPolynomial4DFree(psDPolynomial4D* 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; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + + psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, + "---- Calling ordPolynomial1DEval(%f)\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 %f\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 [%f+= (%f * %f)\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_FLOAT_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; + psF32 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_F32); + if(poly->mask[n-1] == 0) { + d->data.F32[n-1] = poly->coeff[n-1]; + } else { + d->data.F32[n-1] = 0.0; + } + + d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]); + if(poly->mask[n-2] == 0) { + d->data.F32[n-2] += poly->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]); + if(poly->mask[i] == 0) { + d->data.F32[i] += poly->coeff[i]; + } + } + + tmp = (x * d->data.F32[1]) - + (d->data.F32[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; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + psF32 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_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_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; + psF32 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_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); + psS32 loop_x = 0; + psS32 loop_y = 0; + psS32 loop_z = 0; + psS32 i = 0; + psF32 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_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_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; + psF32 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;in; loop_x++) { + if (poly->mask[loop_x] == 0) { + 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? +static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* poly) +{ + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + psVector *d; + psS32 n; + psS32 i; + psF64 tmp; + + n = poly->n; + 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); +} + +static psF64 dOrdPolynomial2DEval(psF64 x, + psF64 y, + const psDPolynomial2D* poly) +{ + 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 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* poly) +{ + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + 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 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly) +{ + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_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 dChebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly) +{ + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); + PS_ASSERT_FLOAT_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. This routine makes use of the + GSL routines for generating both uniformly distributed numbers and the + Gaussian distribution as well. + +XXX: There is no way to seed the random generator. + *****************************************************************************/ +psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) +{ + PS_ASSERT_INT_NONNEGATIVE(Npts, NULL); + + psVector* gauss = NULL; + const gsl_rng_type *T = NULL; + gsl_rng *r = NULL; + psS32 i = 0; + + + gauss = psVectorAlloc(Npts, PS_TYPE_F32); + gauss->n = Npts; + gsl_rng_env_setup(); + T = gsl_rng_default; + r = gsl_rng_alloc(T); + + for (i = 0; i < Npts; i++) { + gauss->data.F32[i] = mean + gsl_ran_gaussian(r, sigma); + } + + // XXX: Should I free r, T as well? This is a memory leak. + return(gauss); +} + +/***************************************************************************** + This routine must allocate memory for the polynomial structures. + *****************************************************************************/ +psPolynomial1D* psPolynomial1DAlloc(int n, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(n, NULL); + + int i = 0; + psPolynomial1D* newPoly = NULL; + + newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); + psMemSetDeallocator(newPoly, (psFreeFunc) polynomial1DFree); + + newPoly->type = type; + newPoly->n = n; + newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32)); + newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32)); + 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( int nX, int nY, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(nX, NULL); + PS_ASSERT_INT_POSITIVE(nY, NULL); + + int x = 0; + 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 = (psF32 **)psAlloc(nX * sizeof(psF32 *)); + newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *)); + newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); + newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); + 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( int nX, int nY, int nZ, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(nX, NULL); + PS_ASSERT_INT_POSITIVE(nY, NULL); + PS_ASSERT_INT_POSITIVE(nZ, NULL); + + psS32 x = 0; + psS32 y = 0; + psS32 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 = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + 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( int nX, int nY, int nZ, 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); + + psS32 x = 0; + psS32 y = 0; + psS32 z = 0; + psS32 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 = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); + newPoly->coeffErr = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); + newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); + newPoly->coeffErr[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); + newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); + newPoly->coeffErr[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); + newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *)); + for (z = 0; z < nZ; z++) { + newPoly->coeff[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); + newPoly->coeffErr[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); + 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); +} + +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_F32, NULL); + + psVector *tmp; + + tmp = psVectorAlloc(x->n, PS_TYPE_F32); + for (psS32 i=0;in;i++) { + tmp->data.F32[i] = psPolynomial1DEval(poly, x->data.F32[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_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, 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_F32); + + // Evaluate the polynomial at the specified points + for (psS32 i=0; idata.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[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_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, 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_F32); + + // Evaluate polynomial + for (psS32 i = 0; i < vecLen; i++) { + tmp->data.F32[i] = psPolynomial3DEval(poly, + x->data.F32[i], + y->data.F32[i], + z->data.F32[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_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(t, NULL); + PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, 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_F32); + + // Evaluate polynomial + for (psS32 i = 0; i < vecLen; i++) { + tmp->data.F32[i] = psPolynomial4DEval(poly, + x->data.F32[i], + y->data.F32[i], + z->data.F32[i], + t->data.F32[i]); + } + + // Return output vector + return(tmp); +} + + +psDPolynomial1D* psDPolynomial1DAlloc( int n, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(n, NULL); + + unsigned int i = 0; + psDPolynomial1D* newPoly = NULL; + + newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); + psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial1DFree); + + newPoly->type = type; + newPoly->n = n; + newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64)); + newPoly->coeffErr = (psF64 *)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); +} + +psDPolynomial2D* psDPolynomial2DAlloc( int nX, int nY, + psPolynomialType type) +{ + PS_ASSERT_INT_POSITIVE(nX, NULL); + PS_ASSERT_INT_POSITIVE(nY, NULL); + + unsigned int x = 0; + unsigned int y = 0; + psDPolynomial2D* newPoly = NULL; + + newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); + psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial2DFree); + + newPoly->type = type; + newPoly->nX = nX; + newPoly->nY = nY; + + newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *)); + newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *)); + newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); + newPoly->coeffErr[x] = (psF64 *)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); +} + +psDPolynomial3D* psDPolynomial3DAlloc( int nX, int nY, 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; + psDPolynomial3D* newPoly = NULL; + + newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); + psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial3DFree); + + newPoly->type = type; + newPoly->nX = nX; + newPoly->nY = nY; + newPoly->nZ = nZ; + + newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); + newPoly->coeffErr[x][y] = (psF64 *)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); +} + +psDPolynomial4D* psDPolynomial4DAlloc( int nX, int nY, int nZ, 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; + psDPolynomial4D* newPoly = NULL; + + newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); + psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial4DFree); + + newPoly->type = type; + newPoly->nX = nX; + newPoly->nY = nY; + newPoly->nZ = nZ; + newPoly->nT = nT; + + newPoly->coeff = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); + newPoly->coeffErr = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); + newPoly->mask = (char ****)psAlloc(nX * sizeof(char ***)); + for (x = 0; x < nX; x++) { + newPoly->coeff[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); + newPoly->coeffErr[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); + newPoly->mask[x] = (char ***)psAlloc(nY * sizeof(char **)); + for (y = 0; y < nY; y++) { + newPoly->coeff[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); + newPoly->coeffErr[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); + newPoly->mask[x][y] = (char **)psAlloc(nZ * sizeof(char *)); + for (z = 0; z < nZ; z++) { + newPoly->coeff[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); + newPoly->coeffErr[x][y][z] = (psF64 *)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); +} + + +psF64 psDPolynomial1DEval(const psDPolynomial1D* poly, psF64 x) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial1DEval(x, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial1DEval(x, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psDPolynomial1DEvalVector(const psDPolynomial1D *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] = psDPolynomial1DEval(poly, + x->data.F64[i]); + } + + return(tmp); +} + + +psF64 psDPolynomial2DEval(const psDPolynomial2D* poly, + psF64 x, + psF64 y) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial2DEval(x, y, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial2DEval(x, y, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psDPolynomial2DEvalVector(const psDPolynomial2D *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 output vector length from minimum length of input vectors + if (y->n < vecLen) { + vecLen = y->n; + } + + // Allocate output vector + tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + + // Evaluate the polynomial + for (psS32 i = 0; i < vecLen; i++) { + tmp->data.F64[i] = psDPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); + } + + // Return output vector + return(tmp); +} + + +psF64 psDPolynomial3DEval(const psDPolynomial3D* poly, + psF64 x, + psF64 y, + psF64 z) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial3DEval(x, y, z, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial3DEval(x, y, z, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psDPolynomial3DEvalVector(const psDPolynomial3D *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 size of output vector from min of 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] = psDPolynomial3DEval(poly, + x->data.F64[i], + y->data.F64[i], + z->data.F64[i]); + } + + // Return output vector + return(tmp); +} + +psF64 psDPolynomial4DEval(const psDPolynomial4D* poly, + psF64 x, + psF64 y, + psF64 z, + psF64 t) +{ + PS_ASSERT_POLY_NON_NULL(poly, NAN); + + if (poly->type == PS_POLYNOMIAL_ORD) { + return(dOrdPolynomial4DEval(x,y,z,t, poly)); + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + return(dChebPolynomial4DEval(x,y,z,t, poly)); + } else { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, + PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, + poly->type); + } + return(NAN); +} + +psVector *psDPolynomial4DEvalVector(const psDPolynomial4D *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 the 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 the polynomial + for (psS32 i = 0; i < vecLen; i++) { + tmp->data.F64[i] = psDPolynomial4DEval(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: Ensure that domain[i+1] != domain[i] + +XXX: What should be the defualty type for knots be? psF32 is assumed. + *****************************************************************************/ +psSpline1D *psSpline1DAlloc( int numSplines, + 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 *tmp = NULL; + psS32 i; + psF32 tmpDomain; + psF32 width; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + tmp->n = numSplines; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + } + + // This should be set by the psVectorFitSpline1D() + tmp->p_psDeriv2 = NULL; + + tmp->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); + width = (max - min) / ((psF32) numSplines); + + tmp->knots->data.F32[0] = min; + tmpDomain = min+width; + for (i=1;iknots->data.F32[i] = tmpDomain; + tmpDomain+= width; + } + tmp->knots->data.F32[numSplines] = max; + + psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); + return(tmp); +} + + +/***************************************************************************** +XXX: What should be the defualty type for knots be? psF32 is assumed. + *****************************************************************************/ +psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, + int order) +{ + PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); + PS_ASSERT_VECTOR_NON_EMPTY(bounds, NULL); + PS_ASSERT_INT_NONNEGATIVE(order, NULL); + + psSpline1D *tmp = NULL; + unsigned int i; + unsigned int numSplines; + + tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); + + numSplines = bounds->n - 1; + tmp->n = numSplines; + + tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); + for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); + } + + // This should be set by the psVectorFitSpline1D() + tmp->p_psDeriv2 = NULL; + + tmp->knots = psVectorAlloc(bounds->n, PS_TYPE_F32); + + for (i=0;in;i++) { + tmp->knots->data.F32[i] = bounds->data.F32[i]; + if (i<(bounds->n-1)) { + if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { + psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n"); + } + } + } + + psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); + return(tmp); +} + +/***************************************************************************** +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(). + *****************************************************************************/ +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: trunk/psLib/src/math/psFunctions.h =================================================================== --- trunk/psLib/src/math/psFunctions.h (revision 4540) +++ trunk/psLib/src/math/psFunctions.h (revision 4540) @@ -0,0 +1,519 @@ +/** @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.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ + * + * 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 + psElemType ctype; ///< Polynomial precision + int n; ///< Number of terms + psF32 *coeff; ///< Coefficients + psF32 *coeffErr; ///< Error in coefficients + char *mask; ///< Coefficient mask +} +psPolynomial1D; + +/** Two-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + psElemType ctype; ///< Polynomial precision + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + psF32 **coeff; ///< Coefficients + psF32 **coeffErr; ///< Error in coefficients + char **mask; ///< Coefficients mask +} +psPolynomial2D; + +/** Three-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + psElemType ctype; ///< Polynomial precision + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + int nZ; ///< Number of terms in z + psF32 ***coeff; ///< Coefficients + psF32 ***coeffErr; ///< Error in coefficients + char ***mask; ///< Coefficients mask +} +psPolynomial3D; + +/** Four-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + psElemType ctype; ///< Polynomial precision + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + int nZ; ///< Number of terms in z + int nT; ///< Number of terms in t + psF32 ****coeff; ///< Coefficients + psF32 ****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( + int n, ///< Number of terms + psPolynomialType type ///< Polynomial Type +); + +/** Allocates a 2-D polynomial structure + * + * @return psPolynomial2D* new 2-D polynomial struct + */ +psPolynomial2D* psPolynomial2DAlloc( + int nX, ///< Number of terms in x + 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( + int nX, ///< Number of terms in x + int nY, ///< Number of terms in y + 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( + int nX, ///< Number of terms in x + int nY, ///< Number of terms in y + int nZ, ///< Number of terms in z + int nT, ///< Number of terms in t + psPolynomialType type ///< Polynomial Type +); + +/** 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 +); + +/*****************************************************************************/ + +/* Double-precision polynomials, mainly for use in astrometry */ + +/** Double-precision one-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int n; ///< Number of terms + psF64 *coeff; ///< Coefficients + psF64 *coeffErr; ///< Error in coefficients + char *mask; ///< Coefficient mask +} +psDPolynomial1D; + +/** Double-precision two-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + psF64 **coeff; ///< Coefficients + psF64 **coeffErr; ///< Error in coefficients + char **mask; ///< Coefficients mask +} +psDPolynomial2D; + +/** Double-precision three-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int nX; ///< Number of terms in x + int nY; ///< Number of terms in y + int nZ; ///< Number of terms in z + psF64 ***coeff; ///< Coefficients + psF64 ***coeffErr; ///< Error in coefficients + char ***mask; ///< Coefficient mask +} +psDPolynomial3D; + +/** Double-precision four-dimensional polynomial */ +typedef struct +{ + psPolynomialType type; ///< Polynomial type + int nX; ///< Number of terms in w + int nY; ///< Number of terms in x + int nZ; ///< Number of terms in y + int nT; ///< Number of terms in z + psF64 ****coeff; ///< Coefficients + psF64 ****coeffErr; ///< Error in coefficients + char ****mask; ///< Coefficients mask +} +psDPolynomial4D; + +/** Allocates a double-precision 1-D polynomial structure with n terms + * + * @return psPolynomial1D* new double-precision 1-D polynomial struct + */ +psDPolynomial1D* psDPolynomial1DAlloc( + int n, ///< Number of terms + psPolynomialType type ///< Polynomial Type +); + +/** Allocates a double-precision 2-D polynomial structure + * + * @return psPolynomial2D* new double-precision 2-D polynomial struct + */ +psDPolynomial2D* psDPolynomial2DAlloc( + int nX, ///< Number of terms in x + int nY, ///< Number of terms in y + psPolynomialType type ///< Polynomial Type +); + +/** Allocates a double-precision 3-D polynomial structure + * + * @return psPolynomial3D* new double-precision 3-D polynomial struct + */ +psDPolynomial3D* psDPolynomial3DAlloc( + int nX, ///< Number of terms in x + int nY, ///< Number of terms in y + int nZ, ///< Number of terms in z + psPolynomialType type ///< Polynomial Type +); + +/** Allocates a double-precision 4-D polynomial structure + * + * @return psPolynomial4D* new double-precision 4-D polynomial struct + */ +psDPolynomial4D* psDPolynomial4DAlloc( + int nX, ///< Number of terms in w + int nY, ///< Number of terms in x + int nZ, ///< Number of terms in y + int nT, ///< Number of terms in z + psPolynomialType type ///< Polynomial Type +); + +/** Evaluates a double-precision 1-D polynomial at specific coordinates. + * + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial1DEval( + const psDPolynomial1D* poly, ///< Coefficients for the polynomial + psF64 x ///< Value at which to evaluate +); + +/** Evaluates a double-precision 2-D polynomial at specific coordinates. + * + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial2DEval( + const psDPolynomial2D* poly, ///< Coefficients for the polynomial + psF64 x, ///< Value x at which to evaluate + psF64 y ///< Value y at which to evaluate +); + +/** Evaluates a double-precision 3-D polynomial at specific coordinates. + * + * @return psF64 result of polynomial at given location + */ +psF64 psDPolynomial3DEval( + const psDPolynomial3D* poly, ///< Coefficients for the polynomial + psF64 x, ///< Value x at which to evaluate + psF64 y, ///< Value y at which to evaluate + psF64 z ///< Value z at which to evaluate +); + +/** Evaluates a double-precision 4-D polynomial at specific coordinates. + * + * @return psF64 result of polynomial at given location + */ +psF64 psDPolynomial4DEval( + const psDPolynomial4D* poly, ///< Coefficients for the polynomial + psF64 x, ///< Value w at which to evaluate + psF64 y, ///< Value x at which to evaluate + psF64 z, ///< Value y at which to evaluate + psF64 t ///< Value z at which to evaluate +); + +/** Evaluates a double-precision 1-D polynomial at specific sets of coordinates. + * + * @return psVector* results of polynomial at given locations + */ +psVector *psDPolynomial1DEvalVector( + const psDPolynomial1D *poly, ///< Coefficients for the polynomial + const psVector *x ///< x locations at which to evaluate +); + +/** Evaluates a double-precision 2-D polynomial at specific sets of coordinates. + * + * @return psVector* results of polynomial at given locations + */ +psVector *psDPolynomial2DEvalVector( + const psDPolynomial2D *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 double-precision 3-D polynomial at specific sets of coordinates. + * + * @return psVector* results of polynomial at given locations + */ +psVector *psDPolynomial3DEvalVector( + const psDPolynomial3D *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 double-precision 4-D polynomial at specific sets of coordinates. + * + * @return psVector* results of polynomial at given locations + */ +psVector *psDPolynomial4DEvalVector( + const psDPolynomial4D *poly, ///< Coefficients for the polynomial + const psVector *x, ///< w locations at which to evaluate + const psVector *y, ///< x locations at which to evaluate + const psVector *z, ///< y locations at which to evaluate + const psVector *t ///< z locations at which to evaluate +); + +/** One-Dimensional Spline */ +typedef struct +{ + unsigned int n; ///< The number of spline polynomials + psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials + psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. + psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. + psVector *knots; ///< The boundaries between each spline piece. Size is n+1. +} +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( + int n, ///< Number of spline polynomials + 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 + 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 + Index: trunk/psLib/src/math/psMatrix.c =================================================================== --- trunk/psLib/src/math/psMatrix.c (revision 4446) +++ trunk/psLib/src/math/psMatrix.c (revision 4540) @@ -21,6 +21,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-25 01:15:01 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -43,5 +43,5 @@ #include "psMatrix.h" #include "psConstants.h" -#include "psDataManipErrors.h" +#include "psErrorText.h" Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 4446) +++ trunk/psLib/src/math/psMinimize.c (revision 4540) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.123 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-20 22:42:30 $ + * @version $Revision: 1.124 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -1357,11 +1357,11 @@ XXX: type F32 is done via vector conversion only. *****************************************************************************/ -psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly, +psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* poly, const psVector* x, const psVector* y, const psVector* yErr) { - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL); + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_INT_NONNEGATIVE(poly->n, NULL); PS_ASSERT_VECTOR_NON_NULL(y, NULL); PS_ASSERT_VECTOR_NON_EMPTY(y, NULL); @@ -1389,5 +1389,5 @@ if (x == NULL) { PS_VECTOR_GEN_X_INDEX_STATIC_F64(x64Static, y->n); - if (myPoly->type == PS_POLYNOMIAL_CHEB) { + if (poly->type == PS_POLYNOMIAL_CHEB) { p_psNormalizeVectorRangeF64(x64Static, -1.0, 1.0); } @@ -1396,5 +1396,5 @@ PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); PS_VECTOR_CONVERT_F32_TO_F64_STATIC(x, x64, x64Static); - if (myPoly->type == PS_POLYNOMIAL_CHEB) { + if (poly->type == PS_POLYNOMIAL_CHEB) { p_psNormalizeVectorRangeF64(x64, -1.0, 1.0); } @@ -1405,8 +1405,8 @@ // Call the appropriate vector fitting routine. psPolynomial1D *rc = NULL; - if (myPoly->type == PS_POLYNOMIAL_CHEB) { - rc = vectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64); - } else if (myPoly->type == PS_POLYNOMIAL_ORD) { - rc = vectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64); + if (poly->type == PS_POLYNOMIAL_CHEB) { + rc = vectorFitPolynomial1DCheby(poly, x64, y64, yErr64); + } else if (poly->type == PS_POLYNOMIAL_ORD) { + rc = vectorFitPolynomial1DOrd(poly, x64, y64, yErr64); } else { psError(PS_ERR_BAD_PARAMETER_VALUE, true, @@ -1419,5 +1419,5 @@ } - return(myPoly); + return(poly); } @@ -1432,6 +1432,6 @@ psMinimization *min = psAlloc(sizeof(psMinimization)); - min->maxIter = maxIter; - min->tol = tol; + *(int*)&min->maxIter = maxIter; + *(float*)&min->tol = tol; min->value = 0.0; min->iter = 0; @@ -1993,6 +1993,6 @@ for (i=0;idata.U8[i] == 0) { - dummyMin.maxIter = PS_MINIMIZE_POWELL_LINEMIN_MAX_ITERATIONS; - dummyMin.tol = PS_MINIMIZE_POWELL_LINEMIN_ERROR_TOLERANCE; + *(int*)&dummyMin.maxIter = PS_MINIMIZE_POWELL_LINEMIN_MAX_ITERATIONS; + *(float*)&dummyMin.tol = PS_MINIMIZE_POWELL_LINEMIN_ERROR_TOLERANCE; mul = p_psLineMin(&dummyMin, Q, Index: trunk/psLib/src/math/psMinimize.h =================================================================== --- trunk/psLib/src/math/psMinimize.h (revision 4446) +++ trunk/psLib/src/math/psMinimize.h (revision 4540) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-21 03:01:37 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -38,19 +38,19 @@ /** A data structure for minimization routines. - * + * * Contains numerical analysis parameters/values */ typedef struct { - psS32 maxIter; ///< Convergence limit - psF32 tol; ///< Error Tolerance - psF32 value; ///< Value of function at minimum - psS32 iter; ///< Number of iterations to date - psF32 lastDelta; ///< The last difference for the fit + const int maxIter; ///< Convergence limit + const float tol; ///< Error Tolerance + float value; ///< Value of function at minimum + int iter; ///< Number of iterations to date + float lastDelta; ///< The last difference for the fit } psMinimization; /** Allocates a psMinimization structure. - * + * * @return psMinimization* : a new psMinimization struct */ @@ -74,5 +74,5 @@ */ psPolynomial1D* psVectorFitPolynomial1D( - psPolynomial1D* myPoly, ///< Polynomial to fit + psPolynomial1D* poly, ///< Polynomial to fit const psVector* x, ///< Ordinates (or NULL to just use the indices) const psVector* y, ///< Coordinates @@ -81,8 +81,8 @@ /** Derive a one-dimensional spline fit. - * - * Given a psSpline1D data structure and a set of x,y vectors, this routine + * + * Given a psSpline1D data structure and a set of x,y vectors, this routine * generates the linear splines which satisfy those data points. - * + * * @return psSpline1D*: the calculated one-dimensional splines */ @@ -96,17 +96,17 @@ /** Specifies the format of a user-defined function that the general Levenberg- * Marquardt minimizer routine will accept. - * - * @return psF64: the single float value of the function given the parameters, + * + * @return float: the single float value of the function given the parameters, * positions, and derivatives. */ typedef -psF64 (*psMinimizeLMChi2Func)( +float (*psMinimizeLMChi2Func)( psVector *deriv, ///< derivatives of the function - psVector *params, ///< the parameters used to evaluate the function - psVector *x ///< positions for evaluation + const psVector *params, ///< the parameters used to evaluate the function + const psVector *x ///< positions for evaluation ); /** Minimizes a specified function based on the Levenberg-Marquardt method. - * + * * @return bool: True if successful. */ @@ -123,5 +123,5 @@ /** Use specified alpha, beta, params to generate a new guess for Alpha, Beta, Params - * + * * @return psBool: True if successful. */ @@ -138,5 +138,5 @@ /** Function used to set parameters for generating "best guess" in minimizing Chi-Squared value. - * + * * @return psF64: Chi-squared value for new guess */ @@ -153,10 +153,10 @@ /** Specifies the format of a user-defined function that the general Powell * minimizer routine will accept. - * - * @return psF32: the single float value of the function given the parameters + * + * @return float: the single float value of the function given the parameters * and coordinate vectors. */ typedef -psF32 (*psMinimizePowellFunc)( +float (*psMinimizePowellFunc)( const psVector *params, ///< Parameters used to evaluate the function const psArray *coords ///< Coordinates at which to evaluate @@ -164,5 +164,5 @@ /** Minimizes a specified function based on the Powell method. - * + * * @return bool: True if successful. */ @@ -177,5 +177,5 @@ /** Calculates the one-dimensional Gaussian in a format acceptable to the Levenberg- * Marquardt minimizer routine. - * + * * @return psVector*: Calculated values */ @@ -188,5 +188,5 @@ /** Calculates the one-dimensional Gaussian in a format acceptable to the Powell * chi-squared minimizer routine. - * + * * @return psVector*: Calculated values */ @@ -198,5 +198,5 @@ /** Specifies the format of a user-defined function that the general Powell chi- * squared minimizer routine will accept. - * + * * @return psVector*: Calculated values given the parameters and coordinates. */ @@ -208,5 +208,5 @@ /** Minimizes a specified function based on the Powell chi-squared method. - * + * * @return bool: True is successful. */ @@ -222,5 +222,5 @@ /** Gauss-Jordan numerical solver. - * + * * @return bool: True if successful. */ Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 4446) +++ (revision ) @@ -1,2201 +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.114 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-29 03:23:19 $ - * - * 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 -#include - -#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 "psDataManipErrors.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 dPolynomial1DFree(psDPolynomial1D* poly); -static void dPolynomial2DFree(psDPolynomial2D* poly); -static void dPolynomial3DFree(psDPolynomial3D* poly); -static void dPolynomial4DFree(psDPolynomial4D* 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); -} - -static void dPolynomial1DFree(psDPolynomial1D* poly) -{ - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void dPolynomial2DFree(psDPolynomial2D* poly) -{ - for (unsigned int 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 dPolynomial3DFree(psDPolynomial3D* 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 dPolynomial4DFree(psDPolynomial4D* 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 psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* poly) -{ - psS32 loop_x = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - - psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, - "---- Calling ordPolynomial1DEval(%f)\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 %f\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 [%f+= (%f * %f)\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 psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* poly) -{ - PS_ASSERT_FLOAT_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; - psF32 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_F32); - if(poly->mask[n-1] == 0) { - d->data.F32[n-1] = poly->coeff[n-1]; - } else { - d->data.F32[n-1] = 0.0; - } - - d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]); - if(poly->mask[n-2] == 0) { - d->data.F32[n-2] += poly->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]); - if(poly->mask[i] == 0) { - d->data.F32[i] += poly->coeff[i]; - } - } - - tmp = (x * d->data.F32[1]) - - (d->data.F32[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 psF32 ordPolynomial2DEval(psF32 x, - psF32 y, - const psPolynomial2D* poly) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - psS32 loop_x = 0; - psS32 loop_y = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - psF32 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 psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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; - psF32 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 psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - psF32 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 psF32 chebPolynomial4DEval(psF32 x, psF32 y, psF32 z, psF32 t, const psPolynomial4D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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; - psF32 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;in; loop_x++) { - if (poly->mask[loop_x] == 0) { - 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? -static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - psVector *d; - psS32 n; - psS32 i; - psF64 tmp; - - n = poly->n; - 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); -} - -static psF64 dOrdPolynomial2DEval(psF64 x, - psF64 y, - const psDPolynomial2D* poly) -{ - 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 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - 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 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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 dChebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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. This routine makes use of the - GSL routines for generating both uniformly distributed numbers and the - Gaussian distribution as well. - -XXX: There is no way to seed the random generator. - *****************************************************************************/ -psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) -{ - PS_ASSERT_INT_NONNEGATIVE(Npts, NULL); - - psVector* gauss = NULL; - const gsl_rng_type *T = NULL; - gsl_rng *r = NULL; - psS32 i = 0; - - - gauss = psVectorAlloc(Npts, PS_TYPE_F32); - gauss->n = Npts; - gsl_rng_env_setup(); - T = gsl_rng_default; - r = gsl_rng_alloc(T); - - for (i = 0; i < Npts; i++) { - gauss->data.F32[i] = mean + gsl_ran_gaussian(r, sigma); - } - - // XXX: Should I free r, T as well? This is a memory leak. - return(gauss); -} - -/***************************************************************************** - This routine must allocate memory for the polynomial structures. - *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(unsigned int n, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(n, NULL); - - int i = 0; - psPolynomial1D* newPoly = NULL; - - newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); - psMemSetDeallocator(newPoly, (psFreeFunc) polynomial1DFree); - - newPoly->type = type; - newPoly->n = n; - newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32)); - newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32)); - newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); - 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); - - int x = 0; - 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 = (psF32 **)psAlloc(nX * sizeof(psF32 *)); - newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *)); - newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); - newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); - newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); - } - 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); - - psS32 x = 0; - psS32 y = 0; - psS32 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 = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); - newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); - newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); - newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); - newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); - newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); - newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); - } - } - 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); - - psS32 x = 0; - psS32 y = 0; - psS32 z = 0; - psS32 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 = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); - newPoly->coeffErr = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); - newPoly->mask = (psU8 ****)psAlloc(nX * sizeof(psU8 ***)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); - newPoly->coeffErr[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); - newPoly->mask[x] = (psU8 ***)psAlloc(nY * sizeof(psU8 **)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); - newPoly->coeffErr[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); - newPoly->mask[x][y] = (psU8 **)psAlloc(nZ * sizeof(psU8 *)); - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); - newPoly->coeffErr[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); - newPoly->mask[x][y][z] = (psU8 *)psAlloc(nT * sizeof(psU8)); - } - } - } - 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); -} - -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_F32, NULL); - - psVector *tmp; - - tmp = psVectorAlloc(x->n, PS_TYPE_F32); - for (psS32 i=0;in;i++) { - tmp->data.F32[i] = psPolynomial1DEval(poly, x->data.F32[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_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, 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_F32); - - // Evaluate the polynomial at the specified points - for (psS32 i=0; idata.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[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_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, 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_F32); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F32[i] = psPolynomial3DEval(poly, - x->data.F32[i], - y->data.F32[i], - z->data.F32[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_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, 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_F32); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F32[i] = psPolynomial4DEval(poly, - x->data.F32[i], - y->data.F32[i], - z->data.F32[i], - t->data.F32[i]); - } - - // Return output vector - return(tmp); -} - - -psDPolynomial1D* psDPolynomial1DAlloc(unsigned int n, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(n, NULL); - - unsigned int i = 0; - psDPolynomial1D* newPoly = NULL; - - newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial1DFree); - - newPoly->type = type; - newPoly->n = n; - newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64)); - newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64)); - newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); - for (i = 0; i < n; i++) { - newPoly->coeff[i] = 0.0; - newPoly->coeffErr[i] = 0.0; - newPoly->mask[i] = 0; - } - - return(newPoly); -} - -psDPolynomial2D* psDPolynomial2DAlloc(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; - psDPolynomial2D* newPoly = NULL; - - newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial2DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - - newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *)); - newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *)); - newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); - newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); - newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); - } - 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); -} - -psDPolynomial3D* psDPolynomial3DAlloc(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; - psDPolynomial3D* newPoly = NULL; - - newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial3DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - - newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); - } - } - 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); -} - -psDPolynomial4D* psDPolynomial4DAlloc(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; - psDPolynomial4D* newPoly = NULL; - - newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial4DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - newPoly->nT = nT; - - newPoly->coeff = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); - newPoly->coeffErr = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); - newPoly->mask = (psU8 ****)psAlloc(nX * sizeof(psU8 ***)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); - newPoly->coeffErr[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); - newPoly->mask[x] = (psU8 ***)psAlloc(nY * sizeof(psU8 **)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); - newPoly->coeffErr[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); - newPoly->mask[x][y] = (psU8 **)psAlloc(nZ * sizeof(psU8 *)); - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); - newPoly->coeffErr[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); - newPoly->mask[x][y][z] = (psU8 *)psAlloc(nT * sizeof(psU8)); - } - } - } - 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); -} - - -psF64 psDPolynomial1DEval(const psDPolynomial1D* poly, psF64 x) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial1DEval(x, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial1DEval(x, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial1DEvalVector(const psDPolynomial1D *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] = psDPolynomial1DEval(poly, - x->data.F64[i]); - } - - return(tmp); -} - - -psF64 psDPolynomial2DEval(const psDPolynomial2D* poly, - psF64 x, - psF64 y) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial2DEval(x, y, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial2DEval(x, y, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial2DEvalVector(const psDPolynomial2D *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 output vector length from minimum length of input vectors - if (y->n < vecLen) { - vecLen = y->n; - } - - // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate the polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); - } - - // Return output vector - return(tmp); -} - - -psF64 psDPolynomial3DEval(const psDPolynomial3D* poly, - psF64 x, - psF64 y, - psF64 z) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial3DEval(x, y, z, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial3DEval(x, y, z, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial3DEvalVector(const psDPolynomial3D *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 size of output vector from min of 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] = psDPolynomial3DEval(poly, - x->data.F64[i], - y->data.F64[i], - z->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -psF64 psDPolynomial4DEval(const psDPolynomial4D* poly, - psF64 x, - psF64 y, - psF64 z, - psF64 t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial4DEval(x,y,z,t, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial4DEval(x,y,z,t, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial4DEvalVector(const psDPolynomial4D *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 the 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 the polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial4DEval(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: Ensure that domain[i+1] != domain[i] - -XXX: What should be the defualty type for knots be? psF32 is assumed. - *****************************************************************************/ -psSpline1D *psSpline1DAlloc(int numSplines, - 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 *tmp = NULL; - psS32 i; - psF32 tmpDomain; - psF32 width; - - tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - tmp->n = numSplines; - - tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); - } - - // This should be set by the psVectorFitSpline1D() - tmp->p_psDeriv2 = NULL; - - tmp->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); - width = (max - min) / ((psF32) numSplines); - - tmp->knots->data.F32[0] = min; - tmpDomain = min+width; - for (i=1;iknots->data.F32[i] = tmpDomain; - tmpDomain+= width; - } - tmp->knots->data.F32[numSplines] = max; - - psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); - return(tmp); -} - - -/***************************************************************************** -XXX: What should be the defualty type for knots be? psF32 is assumed. - *****************************************************************************/ -psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, - int order) -{ - PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(bounds, NULL); - PS_ASSERT_INT_NONNEGATIVE(order, NULL); - - psSpline1D *tmp = NULL; - psS32 i; - psS32 numSplines; - - tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - - numSplines = bounds->n - 1; - tmp->n = numSplines; - - tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); - } - - // This should be set by the psVectorFitSpline1D() - tmp->p_psDeriv2 = NULL; - - tmp->knots = psVectorAlloc(bounds->n, PS_TYPE_F32); - - for (i=0;in;i++) { - tmp->knots->data.F32[i] = bounds->data.F32[i]; - if (i<(bounds->n-1)) { - if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { - psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n"); - } - } - } - - psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); - return(tmp); -} - -/***************************************************************************** -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(). - *****************************************************************************/ -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, - psS32 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); - - psS32 binNum; - psS32 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); - - psS32 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: trunk/psLib/src/math/psPolynomial.h =================================================================== --- trunk/psLib/src/math/psPolynomial.h (revision 4446) +++ (revision ) @@ -1,515 +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.51 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-29 03:23:19 $ - * - * 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 - psF32 *coeff; ///< Coefficients - psF32 *coeffErr; ///< Error in coefficients - psU8 *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 - psF32 **coeff; ///< Coefficients - psF32 **coeffErr; ///< Error in coefficients - psU8 **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 - psF32 ***coeff; ///< Coefficients - psF32 ***coeffErr; ///< Error in coefficients - psU8 ***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 - psF32 ****coeff; ///< Coefficients - psF32 ****coeffErr; ///< Error in coefficients - psU8 ****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 -); - -/** 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 -); - -/*****************************************************************************/ - -/* Double-precision polynomials, mainly for use in astrometry */ - -/** Double-precision one-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int n; ///< Number of terms - psF64 *coeff; ///< Coefficients - psF64 *coeffErr; ///< Error in coefficients - psU8 *mask; ///< Coefficient mask -} -psDPolynomial1D; - -/** Double-precision 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 - psU8 **mask; ///< Coefficients mask -} -psDPolynomial2D; - -/** Double-precision 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 - psU8 ***mask; ///< Coefficient mask -} -psDPolynomial3D; - -/** Double-precision four-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in w - unsigned int nY; ///< Number of terms in x - unsigned int nZ; ///< Number of terms in y - unsigned int nT; ///< Number of terms in z - psF64 ****coeff; ///< Coefficients - psF64 ****coeffErr; ///< Error in coefficients - psU8 ****mask; ///< Coefficients mask -} -psDPolynomial4D; - -/** Allocates a double-precision 1-D polynomial structure with n terms - * - * @return psPolynomial1D* new double-precision 1-D polynomial struct - */ -psDPolynomial1D* psDPolynomial1DAlloc( - unsigned int n, ///< Number of terms - psPolynomialType type ///< Polynomial Type -); - -/** Allocates a double-precision 2-D polynomial structure - * - * @return psPolynomial2D* new double-precision 2-D polynomial struct - */ -psDPolynomial2D* psDPolynomial2DAlloc( - unsigned int nX, ///< Number of terms in x - unsigned int nY, ///< Number of terms in y - psPolynomialType type ///< Polynomial Type -); - -/** Allocates a double-precision 3-D polynomial structure - * - * @return psPolynomial3D* new double-precision 3-D polynomial struct - */ -psDPolynomial3D* psDPolynomial3DAlloc( - 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 double-precision 4-D polynomial structure - * - * @return psPolynomial4D* new double-precision 4-D polynomial struct - */ -psDPolynomial4D* psDPolynomial4DAlloc( - unsigned int nX, ///< Number of terms in w - unsigned int nY, ///< Number of terms in x - unsigned int nZ, ///< Number of terms in y - unsigned int nT, ///< Number of terms in z - psPolynomialType type ///< Polynomial Type -); - -/** Evaluates a double-precision 1-D polynomial at specific coordinates. - * - * @return psF32 result of polynomial at given location - */ -psF64 psDPolynomial1DEval( - const psDPolynomial1D* poly, ///< Coefficients for the polynomial - psF64 x ///< Value at which to evaluate -); - -/** Evaluates a double-precision 2-D polynomial at specific coordinates. - * - * @return psF32 result of polynomial at given location - */ -psF64 psDPolynomial2DEval( - const psDPolynomial2D* poly, ///< Coefficients for the polynomial - psF64 x, ///< Value x at which to evaluate - psF64 y ///< Value y at which to evaluate -); - -/** Evaluates a double-precision 3-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psDPolynomial3DEval( - const psDPolynomial3D* poly, ///< Coefficients for the polynomial - psF64 x, ///< Value x at which to evaluate - psF64 y, ///< Value y at which to evaluate - psF64 z ///< Value z at which to evaluate -); - -/** Evaluates a double-precision 4-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psDPolynomial4DEval( - const psDPolynomial4D* poly, ///< Coefficients for the polynomial - psF64 x, ///< Value w at which to evaluate - psF64 y, ///< Value x at which to evaluate - psF64 z, ///< Value y at which to evaluate - psF64 t ///< Value z at which to evaluate -); - -/** Evaluates a double-precision 1-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial1DEvalVector( - const psDPolynomial1D *poly, ///< Coefficients for the polynomial - const psVector *x ///< x locations at which to evaluate -); - -/** Evaluates a double-precision 2-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial2DEvalVector( - const psDPolynomial2D *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 double-precision 3-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial3DEvalVector( - const psDPolynomial3D *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 double-precision 4-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial4DEvalVector( - const psDPolynomial4D *poly, ///< Coefficients for the polynomial - const psVector *x, ///< w locations at which to evaluate - const psVector *y, ///< x locations at which to evaluate - const psVector *z, ///< y locations at which to evaluate - const psVector *t ///< z locations at which to evaluate -); - -/** One-Dimensional Spline */ -typedef struct -{ - psS32 n; ///< The number of spline polynomials - psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials - psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. - psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. - psVector *knots; ///< The boundaries between each spline piece. Size is n+1. -} -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( - int n, ///< Number of spline polynomials - 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 - 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 - psS32 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 - Index: trunk/psLib/src/math/psRandom.c =================================================================== --- trunk/psLib/src/math/psRandom.c (revision 4446) +++ trunk/psLib/src/math/psRandom.c (revision 4540) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-06-21 03:01:37 $ +* @version $Revision: 1.10 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -31,5 +31,5 @@ #include "psLogMsg.h" #include "psConstants.h" -#include "psDataManipErrors.h" +#include "psErrorText.h" psU64 p_psRandomGetSystemSeed() Index: trunk/psLib/src/math/psSpline.c =================================================================== --- trunk/psLib/src/math/psSpline.c (revision 4446) +++ (revision ) @@ -1,2201 +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.114 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-29 03:23:19 $ - * - * 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 -#include - -#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 "psDataManipErrors.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 dPolynomial1DFree(psDPolynomial1D* poly); -static void dPolynomial2DFree(psDPolynomial2D* poly); -static void dPolynomial3DFree(psDPolynomial3D* poly); -static void dPolynomial4DFree(psDPolynomial4D* 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); -} - -static void dPolynomial1DFree(psDPolynomial1D* poly) -{ - psFree(poly->coeff); - psFree(poly->coeffErr); - psFree(poly->mask); -} - -static void dPolynomial2DFree(psDPolynomial2D* poly) -{ - for (unsigned int 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 dPolynomial3DFree(psDPolynomial3D* 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 dPolynomial4DFree(psDPolynomial4D* 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 psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* poly) -{ - psS32 loop_x = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - - psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, - "---- Calling ordPolynomial1DEval(%f)\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 %f\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 [%f+= (%f * %f)\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 psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* poly) -{ - PS_ASSERT_FLOAT_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; - psF32 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_F32); - if(poly->mask[n-1] == 0) { - d->data.F32[n-1] = poly->coeff[n-1]; - } else { - d->data.F32[n-1] = 0.0; - } - - d->data.F32[n-2] = (2.0 * x * d->data.F32[n-1]); - if(poly->mask[n-2] == 0) { - d->data.F32[n-2] += poly->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]); - if(poly->mask[i] == 0) { - d->data.F32[i] += poly->coeff[i]; - } - } - - tmp = (x * d->data.F32[1]) - - (d->data.F32[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 psF32 ordPolynomial2DEval(psF32 x, - psF32 y, - const psPolynomial2D* poly) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - psS32 loop_x = 0; - psS32 loop_y = 0; - psF32 polySum = 0.0; - psF32 xSum = 1.0; - psF32 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 psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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; - psF32 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 psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - psS32 loop_x = 0; - psS32 loop_y = 0; - psS32 loop_z = 0; - psS32 i = 0; - psF32 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 psF32 chebPolynomial4DEval(psF32 x, psF32 y, psF32 z, psF32 t, const psPolynomial4D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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; - psF32 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;in; loop_x++) { - if (poly->mask[loop_x] == 0) { - 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? -static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - psVector *d; - psS32 n; - psS32 i; - psF64 tmp; - - n = poly->n; - 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); -} - -static psF64 dOrdPolynomial2DEval(psF64 x, - psF64 y, - const psDPolynomial2D* poly) -{ - 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 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - 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 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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 dChebPolynomial4DEval(psF64 x, psF64 y, psF64 z, psF64 t, const psDPolynomial4D* poly) -{ - PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_WITHIN_RANGE(z, -1.0, 1.0, 0.0); - PS_ASSERT_FLOAT_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. This routine makes use of the - GSL routines for generating both uniformly distributed numbers and the - Gaussian distribution as well. - -XXX: There is no way to seed the random generator. - *****************************************************************************/ -psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) -{ - PS_ASSERT_INT_NONNEGATIVE(Npts, NULL); - - psVector* gauss = NULL; - const gsl_rng_type *T = NULL; - gsl_rng *r = NULL; - psS32 i = 0; - - - gauss = psVectorAlloc(Npts, PS_TYPE_F32); - gauss->n = Npts; - gsl_rng_env_setup(); - T = gsl_rng_default; - r = gsl_rng_alloc(T); - - for (i = 0; i < Npts; i++) { - gauss->data.F32[i] = mean + gsl_ran_gaussian(r, sigma); - } - - // XXX: Should I free r, T as well? This is a memory leak. - return(gauss); -} - -/***************************************************************************** - This routine must allocate memory for the polynomial structures. - *****************************************************************************/ -psPolynomial1D* psPolynomial1DAlloc(unsigned int n, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(n, NULL); - - int i = 0; - psPolynomial1D* newPoly = NULL; - - newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); - psMemSetDeallocator(newPoly, (psFreeFunc) polynomial1DFree); - - newPoly->type = type; - newPoly->n = n; - newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32)); - newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32)); - newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); - 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); - - int x = 0; - 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 = (psF32 **)psAlloc(nX * sizeof(psF32 *)); - newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *)); - newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); - newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); - newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); - } - 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); - - psS32 x = 0; - psS32 y = 0; - psS32 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 = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); - newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); - newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); - newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); - newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); - newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); - newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); - } - } - 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); - - psS32 x = 0; - psS32 y = 0; - psS32 z = 0; - psS32 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 = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); - newPoly->coeffErr = (psF32 ****)psAlloc(nX * sizeof(psF32 ***)); - newPoly->mask = (psU8 ****)psAlloc(nX * sizeof(psU8 ***)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); - newPoly->coeffErr[x] = (psF32 ***)psAlloc(nY * sizeof(psF32 **)); - newPoly->mask[x] = (psU8 ***)psAlloc(nY * sizeof(psU8 **)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); - newPoly->coeffErr[x][y] = (psF32 **)psAlloc(nZ * sizeof(psF32 *)); - newPoly->mask[x][y] = (psU8 **)psAlloc(nZ * sizeof(psU8 *)); - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); - newPoly->coeffErr[x][y][z] = (psF32 *)psAlloc(nT * sizeof(psF32)); - newPoly->mask[x][y][z] = (psU8 *)psAlloc(nT * sizeof(psU8)); - } - } - } - 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); -} - -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_F32, NULL); - - psVector *tmp; - - tmp = psVectorAlloc(x->n, PS_TYPE_F32); - for (psS32 i=0;in;i++) { - tmp->data.F32[i] = psPolynomial1DEval(poly, x->data.F32[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_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, 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_F32); - - // Evaluate the polynomial at the specified points - for (psS32 i=0; idata.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[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_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, 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_F32); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F32[i] = psPolynomial3DEval(poly, - x->data.F32[i], - y->data.F32[i], - z->data.F32[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_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, 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_F32); - - // Evaluate polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F32[i] = psPolynomial4DEval(poly, - x->data.F32[i], - y->data.F32[i], - z->data.F32[i], - t->data.F32[i]); - } - - // Return output vector - return(tmp); -} - - -psDPolynomial1D* psDPolynomial1DAlloc(unsigned int n, - psPolynomialType type) -{ - PS_ASSERT_INT_POSITIVE(n, NULL); - - unsigned int i = 0; - psDPolynomial1D* newPoly = NULL; - - newPoly = (psDPolynomial1D* ) psAlloc(sizeof(psDPolynomial1D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial1DFree); - - newPoly->type = type; - newPoly->n = n; - newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64)); - newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64)); - newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); - for (i = 0; i < n; i++) { - newPoly->coeff[i] = 0.0; - newPoly->coeffErr[i] = 0.0; - newPoly->mask[i] = 0; - } - - return(newPoly); -} - -psDPolynomial2D* psDPolynomial2DAlloc(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; - psDPolynomial2D* newPoly = NULL; - - newPoly = (psDPolynomial2D* ) psAlloc(sizeof(psDPolynomial2D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial2DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - - newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *)); - newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *)); - newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); - newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); - newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); - } - 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); -} - -psDPolynomial3D* psDPolynomial3DAlloc(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; - psDPolynomial3D* newPoly = NULL; - - newPoly = (psDPolynomial3D* ) psAlloc(sizeof(psDPolynomial3D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial3DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - - newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); - newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); - newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); - newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); - } - } - 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); -} - -psDPolynomial4D* psDPolynomial4DAlloc(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; - psDPolynomial4D* newPoly = NULL; - - newPoly = (psDPolynomial4D* ) psAlloc(sizeof(psDPolynomial4D)); - psMemSetDeallocator(newPoly, (psFreeFunc) dPolynomial4DFree); - - newPoly->type = type; - newPoly->nX = nX; - newPoly->nY = nY; - newPoly->nZ = nZ; - newPoly->nT = nT; - - newPoly->coeff = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); - newPoly->coeffErr = (psF64 ****)psAlloc(nX * sizeof(psF64 ***)); - newPoly->mask = (psU8 ****)psAlloc(nX * sizeof(psU8 ***)); - for (x = 0; x < nX; x++) { - newPoly->coeff[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); - newPoly->coeffErr[x] = (psF64 ***)psAlloc(nY * sizeof(psF64 **)); - newPoly->mask[x] = (psU8 ***)psAlloc(nY * sizeof(psU8 **)); - for (y = 0; y < nY; y++) { - newPoly->coeff[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); - newPoly->coeffErr[x][y] = (psF64 **)psAlloc(nZ * sizeof(psF64 *)); - newPoly->mask[x][y] = (psU8 **)psAlloc(nZ * sizeof(psU8 *)); - for (z = 0; z < nZ; z++) { - newPoly->coeff[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); - newPoly->coeffErr[x][y][z] = (psF64 *)psAlloc(nT * sizeof(psF64)); - newPoly->mask[x][y][z] = (psU8 *)psAlloc(nT * sizeof(psU8)); - } - } - } - 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); -} - - -psF64 psDPolynomial1DEval(const psDPolynomial1D* poly, psF64 x) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial1DEval(x, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial1DEval(x, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial1DEvalVector(const psDPolynomial1D *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] = psDPolynomial1DEval(poly, - x->data.F64[i]); - } - - return(tmp); -} - - -psF64 psDPolynomial2DEval(const psDPolynomial2D* poly, - psF64 x, - psF64 y) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial2DEval(x, y, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial2DEval(x, y, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial2DEvalVector(const psDPolynomial2D *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 output vector length from minimum length of input vectors - if (y->n < vecLen) { - vecLen = y->n; - } - - // Allocate output vector - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate the polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]); - } - - // Return output vector - return(tmp); -} - - -psF64 psDPolynomial3DEval(const psDPolynomial3D* poly, - psF64 x, - psF64 y, - psF64 z) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial3DEval(x, y, z, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial3DEval(x, y, z, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial3DEvalVector(const psDPolynomial3D *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 size of output vector from min of 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] = psDPolynomial3DEval(poly, - x->data.F64[i], - y->data.F64[i], - z->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -psF64 psDPolynomial4DEval(const psDPolynomial4D* poly, - psF64 x, - psF64 y, - psF64 z, - psF64 t) -{ - PS_ASSERT_POLY_NON_NULL(poly, NAN); - - if (poly->type == PS_POLYNOMIAL_ORD) { - return(dOrdPolynomial4DEval(x,y,z,t, poly)); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - return(dChebPolynomial4DEval(x,y,z,t, poly)); - } else { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psFunctions_INVALID_POLYNOMIAL_TYPE, - poly->type); - } - return(NAN); -} - -psVector *psDPolynomial4DEvalVector(const psDPolynomial4D *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 the 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 the polynomial - for (psS32 i = 0; i < vecLen; i++) { - tmp->data.F64[i] = psDPolynomial4DEval(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: Ensure that domain[i+1] != domain[i] - -XXX: What should be the defualty type for knots be? psF32 is assumed. - *****************************************************************************/ -psSpline1D *psSpline1DAlloc(int numSplines, - 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 *tmp = NULL; - psS32 i; - psF32 tmpDomain; - psF32 width; - - tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - tmp->n = numSplines; - - tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); - } - - // This should be set by the psVectorFitSpline1D() - tmp->p_psDeriv2 = NULL; - - tmp->knots = psVectorAlloc(numSplines+1, PS_TYPE_F32); - width = (max - min) / ((psF32) numSplines); - - tmp->knots->data.F32[0] = min; - tmpDomain = min+width; - for (i=1;iknots->data.F32[i] = tmpDomain; - tmpDomain+= width; - } - tmp->knots->data.F32[numSplines] = max; - - psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); - return(tmp); -} - - -/***************************************************************************** -XXX: What should be the defualty type for knots be? psF32 is assumed. - *****************************************************************************/ -psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, - int order) -{ - PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(bounds, NULL); - PS_ASSERT_INT_NONNEGATIVE(order, NULL); - - psSpline1D *tmp = NULL; - psS32 i; - psS32 numSplines; - - tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); - - numSplines = bounds->n - 1; - tmp->n = numSplines; - - tmp->spline = (psPolynomial1D **) psAlloc(numSplines * sizeof(psPolynomial1D *)); - for (i=0;ispline)[i] = psPolynomial1DAlloc(order+1, PS_POLYNOMIAL_ORD); - } - - // This should be set by the psVectorFitSpline1D() - tmp->p_psDeriv2 = NULL; - - tmp->knots = psVectorAlloc(bounds->n, PS_TYPE_F32); - - for (i=0;in;i++) { - tmp->knots->data.F32[i] = bounds->data.F32[i]; - if (i<(bounds->n-1)) { - if (FLT_EPSILON >= fabs(bounds->data.F32[i+1]-bounds->data.F32[i])) { - psError(PS_ERR_UNKNOWN, true, "data points must be distinct\n"); - } - } - } - - psMemSetDeallocator(tmp,(psFreeFunc)spline1DFree); - return(tmp); -} - -/***************************************************************************** -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(). - *****************************************************************************/ -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, - psS32 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); - - psS32 binNum; - psS32 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); - - psS32 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: trunk/psLib/src/math/psSpline.h =================================================================== --- trunk/psLib/src/math/psSpline.h (revision 4446) +++ (revision ) @@ -1,515 +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.51 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-29 03:23:19 $ - * - * 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 - psF32 *coeff; ///< Coefficients - psF32 *coeffErr; ///< Error in coefficients - psU8 *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 - psF32 **coeff; ///< Coefficients - psF32 **coeffErr; ///< Error in coefficients - psU8 **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 - psF32 ***coeff; ///< Coefficients - psF32 ***coeffErr; ///< Error in coefficients - psU8 ***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 - psF32 ****coeff; ///< Coefficients - psF32 ****coeffErr; ///< Error in coefficients - psU8 ****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 -); - -/** 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 -); - -/*****************************************************************************/ - -/* Double-precision polynomials, mainly for use in astrometry */ - -/** Double-precision one-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int n; ///< Number of terms - psF64 *coeff; ///< Coefficients - psF64 *coeffErr; ///< Error in coefficients - psU8 *mask; ///< Coefficient mask -} -psDPolynomial1D; - -/** Double-precision 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 - psU8 **mask; ///< Coefficients mask -} -psDPolynomial2D; - -/** Double-precision 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 - psU8 ***mask; ///< Coefficient mask -} -psDPolynomial3D; - -/** Double-precision four-dimensional polynomial */ -typedef struct -{ - psPolynomialType type; ///< Polynomial type - unsigned int nX; ///< Number of terms in w - unsigned int nY; ///< Number of terms in x - unsigned int nZ; ///< Number of terms in y - unsigned int nT; ///< Number of terms in z - psF64 ****coeff; ///< Coefficients - psF64 ****coeffErr; ///< Error in coefficients - psU8 ****mask; ///< Coefficients mask -} -psDPolynomial4D; - -/** Allocates a double-precision 1-D polynomial structure with n terms - * - * @return psPolynomial1D* new double-precision 1-D polynomial struct - */ -psDPolynomial1D* psDPolynomial1DAlloc( - unsigned int n, ///< Number of terms - psPolynomialType type ///< Polynomial Type -); - -/** Allocates a double-precision 2-D polynomial structure - * - * @return psPolynomial2D* new double-precision 2-D polynomial struct - */ -psDPolynomial2D* psDPolynomial2DAlloc( - unsigned int nX, ///< Number of terms in x - unsigned int nY, ///< Number of terms in y - psPolynomialType type ///< Polynomial Type -); - -/** Allocates a double-precision 3-D polynomial structure - * - * @return psPolynomial3D* new double-precision 3-D polynomial struct - */ -psDPolynomial3D* psDPolynomial3DAlloc( - 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 double-precision 4-D polynomial structure - * - * @return psPolynomial4D* new double-precision 4-D polynomial struct - */ -psDPolynomial4D* psDPolynomial4DAlloc( - unsigned int nX, ///< Number of terms in w - unsigned int nY, ///< Number of terms in x - unsigned int nZ, ///< Number of terms in y - unsigned int nT, ///< Number of terms in z - psPolynomialType type ///< Polynomial Type -); - -/** Evaluates a double-precision 1-D polynomial at specific coordinates. - * - * @return psF32 result of polynomial at given location - */ -psF64 psDPolynomial1DEval( - const psDPolynomial1D* poly, ///< Coefficients for the polynomial - psF64 x ///< Value at which to evaluate -); - -/** Evaluates a double-precision 2-D polynomial at specific coordinates. - * - * @return psF32 result of polynomial at given location - */ -psF64 psDPolynomial2DEval( - const psDPolynomial2D* poly, ///< Coefficients for the polynomial - psF64 x, ///< Value x at which to evaluate - psF64 y ///< Value y at which to evaluate -); - -/** Evaluates a double-precision 3-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psDPolynomial3DEval( - const psDPolynomial3D* poly, ///< Coefficients for the polynomial - psF64 x, ///< Value x at which to evaluate - psF64 y, ///< Value y at which to evaluate - psF64 z ///< Value z at which to evaluate -); - -/** Evaluates a double-precision 4-D polynomial at specific coordinates. - * - * @return psF64 result of polynomial at given location - */ -psF64 psDPolynomial4DEval( - const psDPolynomial4D* poly, ///< Coefficients for the polynomial - psF64 x, ///< Value w at which to evaluate - psF64 y, ///< Value x at which to evaluate - psF64 z, ///< Value y at which to evaluate - psF64 t ///< Value z at which to evaluate -); - -/** Evaluates a double-precision 1-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial1DEvalVector( - const psDPolynomial1D *poly, ///< Coefficients for the polynomial - const psVector *x ///< x locations at which to evaluate -); - -/** Evaluates a double-precision 2-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial2DEvalVector( - const psDPolynomial2D *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 double-precision 3-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial3DEvalVector( - const psDPolynomial3D *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 double-precision 4-D polynomial at specific sets of coordinates. - * - * @return psVector* results of polynomial at given locations - */ -psVector *psDPolynomial4DEvalVector( - const psDPolynomial4D *poly, ///< Coefficients for the polynomial - const psVector *x, ///< w locations at which to evaluate - const psVector *y, ///< x locations at which to evaluate - const psVector *z, ///< y locations at which to evaluate - const psVector *t ///< z locations at which to evaluate -); - -/** One-Dimensional Spline */ -typedef struct -{ - psS32 n; ///< The number of spline polynomials - psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials - psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. - psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. - psVector *knots; ///< The boundaries between each spline piece. Size is n+1. -} -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( - int n, ///< Number of spline polynomials - 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 - 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 - psS32 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 - Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 4446) +++ trunk/psLib/src/math/psStats.c (revision 4540) @@ -14,6 +14,6 @@ * stats->binsize * - * @version $Revision: 1.136 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-28 20:17:52 $ + * @version $Revision: 1.137 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -41,5 +41,5 @@ #include "psConstants.h" -#include "psDataManipErrors.h" +#include "psErrorText.h" /*****************************************************************************/ @@ -77,5 +77,6 @@ psBool p_psGetStatValue(const psStats* stats, psF64 *value) { - switch (stats->options & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE)) { + // switch (stats->options & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE)) { + switch (stats->options & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) { case PS_STAT_SAMPLE_MEAN: *value = stats->sampleMean; @@ -658,5 +659,5 @@ x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMin = p_psVectorBinDisect(robustHistogram->bounds, &x); + jMin = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); if (jMin < 0) { psError(PS_ERR_UNEXPECTED_NULL, @@ -673,5 +674,5 @@ x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMax = p_psVectorBinDisect(robustHistogram->bounds, &x); + jMax = p_psVectorBinDisect( *(psVector* *)&robustHistogram->bounds, &x); if (jMax < 0) { psError(PS_ERR_UNEXPECTED_NULL, @@ -1705,5 +1706,5 @@ //Determine id that is okay. stats->robustMedian = fitQuadraticSearchForYThenReturnX( - robustHistogram->bounds, + *(psVector* *)&robustHistogram->bounds, cumulativeRobustSums, medianBinNum, @@ -1721,5 +1722,5 @@ //Determine id that is okay. stats->robustLQ = fitQuadraticSearchForYThenReturnX( - robustHistogram->bounds, + *(psVector* *)&robustHistogram->bounds, cumulativeRobustSums, LQBinNum, @@ -1728,5 +1729,5 @@ //Determine id that is okay. stats->robustUQ = fitQuadraticSearchForYThenReturnX( - robustHistogram->bounds, + *(psVector* *)&robustHistogram->bounds, cumulativeRobustSums, UQBinNum, @@ -1812,5 +1813,5 @@ psMemSetDeallocator(newHist, (psFreeFunc) histogramFree); newHist->bounds = psVectorAlloc(n + 1, PS_TYPE_F32); - newHist->bounds->n = newHist->bounds->nalloc; + *(int *)&newHist->bounds->n = newHist->bounds->nalloc; // Calculate the bounds for each bin. @@ -1859,5 +1860,5 @@ psMemSetDeallocator(newHist, (psFreeFunc) histogramFree); newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32); - newHist->bounds->n = newHist->bounds->nalloc; + *(int *)&newHist->bounds->n = newHist->bounds->nalloc; for (i = 0; i < bounds->n; i++) { newHist->bounds->data.F32[i] = bounds->data.F32[i]; @@ -2067,5 +2068,5 @@ // correct bin number requires a bit more work. tmpScalar.data.F32 = inF32->data.F32[i]; - binNum = p_psVectorBinDisect(out->bounds, &tmpScalar); + binNum = p_psVectorBinDisect( *(psVector* *)&out->bounds, &tmpScalar); if (binNum < 0) { psLogMsg(__func__, PS_LOG_WARN, Index: trunk/psLib/src/math/psStats.h =================================================================== --- trunk/psLib/src/math/psStats.h (revision 4446) +++ trunk/psLib/src/math/psStats.h (revision 4540) @@ -14,6 +14,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-28 20:17:52 $ + * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -52,5 +52,5 @@ PS_STAT_USE_RANGE = 0x002000, ///< Range PS_STAT_USE_BINSIZE = 0x004000, ///< Binsize - PS_STAT_ROBUST_FOR_SAMPLE = 0x008000 ///< Robust for sample + // PS_STAT_ROBUST_FOR_SAMPLE = 0x008000 ///< Robust for sample } psStatsOptions; @@ -60,26 +60,26 @@ typedef struct { - psF64 sampleMean; ///< formal mean of sample - psF64 sampleMedian; ///< formal median of sample - psF64 sampleStdev; ///< standard deviation of sample - psF64 sampleUQ; ///< upper quartile of sample - psF64 sampleLQ; ///< lower quartile of sample - psF64 robustMean; ///< robust mean of array - psF64 robustMedian; ///< robust median of array - psF64 robustMode; ///< Robust mode of array - psF64 robustStdev; ///< robust standard deviation of array - psF64 robustUQ; ///< robust upper quartile - psF64 robustLQ; ///< robust lower quartile - psS32 robustN50; ///< Number of points in Gaussian fit. XXX: This is currently never set. - psS32 robustNfit; ///< The number of points between the quartiles. - psF64 clippedMean; ///< Nsigma clipped mean - psF64 clippedStdev; ///< standard deviation after clipping - psS32 clippedNvalues; ///< Number of data points used for clipped mean: This value is never used. - psF64 clipSigma; ///< Nsigma used for clipping; user input - psS32 clipIter; ///< Number of clipping iterations; user input - psF64 min; ///< minimum data value in array - psF64 max; ///< maximum data value in array - psF64 binsize; ///< binsize for robust fit (input/ouput) - psStatsOptions options; ///< bitmask of calculated values + double sampleMean; ///< formal mean of sample + double sampleMedian; ///< formal median of sample + double sampleStdev; ///< standard deviation of sample + double sampleUQ; ///< upper quartile of sample + double sampleLQ; ///< lower quartile of sample + double robustMean; ///< robust mean of array + double robustMedian; ///< robust median of array + double robustMode; ///< Robust mode of array + double robustStdev; ///< robust standard deviation of array + double robustUQ; ///< robust upper quartile + double robustLQ; ///< robust lower quartile + int robustN50; ///< Number of points in Gaussian fit. XXX: This is currently never set. + int robustNfit; ///< The number of points between the quartiles. + double clippedMean; ///< Nsigma clipped mean + double clippedStdev; ///< standard deviation after clipping + int clippedNvalues; ///< Number of data points used for clipped mean: This value is never used. + double clipSigma; ///< Nsigma used for clipping; user input + int clipIter; ///< Number of clipping iterations; user input + double min; ///< minimum data value in array + double max; ///< maximum data value in array + double binsize; ///< binsize for robust fit (input/ouput) + psStatsOptions options; ///< bitmask of calculated values } psStats; @@ -90,9 +90,9 @@ */ psStats* psVectorStats( - psStats* stats, ///< stats structure defines stats to be calculated and how - const psVector* in, ///< Vector to be analysed. - const psVector* errors, ///< Errors. - const psVector* mask, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL - psMaskType maskVal ///< Only mask elements with one of these bits set in maskVector + psStats* stats, ///< stats structure defines stats to be calculated and how + const psVector* in, ///< Vector to be analysed. + const psVector* errors, ///< Errors. + const psVector* mask, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL + psMaskType maskVal ///< Only mask elements with one of these bits set in maskVector ); @@ -121,9 +121,9 @@ typedef struct { - psVector* bounds; ///< Bounds for the bins (type F32) + const psVector* bounds; ///< Bounds for the bins (type F32) psVector* nums; ///< Number in each of the bins (INT) - psS32 minNum; ///< Number below the minimum - psS32 maxNum; ///< Number above the maximum - psBool uniform; ///< Is it a uniform distribution? + int minNum; ///< Number below the minimum + int maxNum; ///< Number above the maximum + bool uniform; ///< Is it a uniform distribution? } psHistogram; Index: trunk/psLib/src/math/psUnaryOp.c =================================================================== --- trunk/psLib/src/math/psUnaryOp.c (revision 4446) +++ trunk/psLib/src/math/psUnaryOp.c (revision 4540) @@ -30,6 +30,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-28 20:17:52 $ + * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -51,5 +51,5 @@ #include "psLogMsg.h" #include "psConstants.h" -#include "psDataManipErrors.h" +#include "psErrorText.h" /***************************************************************************** @@ -304,5 +304,5 @@ } -psPtr psUnaryOp(psPtr out, const psPtr in, const char *op) +psMathType* psUnaryOp(psPtr out, const psPtr in, const char *op) { #define psUnaryOp_EXIT { \ Index: trunk/psLib/src/math/psUnaryOp.h =================================================================== --- trunk/psLib/src/math/psUnaryOp.h (revision 4446) +++ trunk/psLib/src/math/psUnaryOp.h (revision 4540) @@ -30,6 +30,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-06-28 20:17:52 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-07-12 19:12:01 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -59,5 +59,5 @@ * @return psType* : Pointer to either psImage or psVector. */ -psType* psUnaryOp( +psMathType* psUnaryOp( psPtr out, ///< Output type, either psImage or psVector. const psPtr in, ///< Input, either psImage or psVector.