Index: /trunk/psLib/src/dataManip/psConstants.h =================================================================== --- /trunk/psLib/src/dataManip/psConstants.h (revision 2787) +++ /trunk/psLib/src/dataManip/psConstants.h (revision 2788) @@ -6,8 +6,14 @@ * @author GLG, MHPCC * - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:42:07 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * XXX: Add parenthesis around all arguments so that these macros can be + * called with complex expressions. + * + * + * */ Index: /trunk/psLib/src/dataManip/psFunctions.c =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.c (revision 2787) +++ /trunk/psLib/src/dataManip/psFunctions.c (revision 2788) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.74 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:42:07 $ + * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -55,5 +55,5 @@ static void dPolynomial4DFree(psDPolynomial4D* myPoly); static void spline1DFree(psSpline1D *tmpSpline); -static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x); +static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); @@ -277,9 +277,9 @@ *****************************************************************************/ -static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly) +static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) { psS32 loop_x = 0; - float polySum = 0.0; - float xSum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, @@ -307,5 +307,5 @@ // XXX: How does the mask vector effect Crenshaw's formula? // XXX: We assume that x is scaled between -1.0 and 1.0; -static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) { // PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0); @@ -314,5 +314,5 @@ psS32 n; psS32 i; - float tmp; + psF32 tmp; n = myPoly->n; @@ -337,5 +337,5 @@ psS32 n; psS32 i; - float tmp; + psF32 tmp; psPolynomial1D **chebPolys = NULL; @@ -355,5 +355,5 @@ } -static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -361,7 +361,7 @@ psS32 loop_x = 0; psS32 loop_y = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -379,5 +379,5 @@ } -static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -386,5 +386,5 @@ psS32 loop_y = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -414,13 +414,13 @@ } -static float ordPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; + psF32 zSum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -442,5 +442,5 @@ } -static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) { psS32 loop_x = 0; @@ -448,5 +448,5 @@ psS32 loop_z = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -483,5 +483,5 @@ } -static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -489,9 +489,9 @@ psS32 loop_y = 0; psS32 loop_z = 0; - float polySum = 0.0; - float wSum = 1.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; + psF32 polySum = 0.0; + psF32 wSum = 1.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; + psF32 zSum = 1.0; for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { @@ -517,5 +517,5 @@ } -static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -524,5 +524,5 @@ psS32 loop_z = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -568,9 +568,9 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) { psS32 loop_x = 0; - double polySum = 0.0; - double xSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; for (loop_x = 0; loop_x < myPoly->n; loop_x++) { @@ -586,10 +586,10 @@ // XXX: You can do this without having to psAlloc() vector d. // XXX: How does the mask vector effect Crenshaw's formula? -static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) { psVector *d; psS32 n; psS32 i; - double tmp; + psF64 tmp; n = myPoly->n; @@ -611,11 +611,11 @@ } -static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -633,10 +633,10 @@ } -static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -667,13 +667,13 @@ } -static double dOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -695,5 +695,5 @@ } -static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) { psS32 loop_x = 0; @@ -701,5 +701,5 @@ psS32 loop_z = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -736,5 +736,5 @@ } -static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -742,9 +742,9 @@ psS32 loop_y = 0; psS32 loop_z = 0; - double polySum = 0.0; - double wSum = 1.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; + psF64 polySum = 0.0; + psF64 wSum = 1.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { @@ -770,5 +770,5 @@ } -static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -777,5 +777,5 @@ psS32 loop_z = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -828,5 +828,5 @@ *****************************************************************************/ #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ -static float fullInterpolate1D##TYPE(ps##TYPE *domain, \ +static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ ps##TYPE *range, \ psS32 n, \ @@ -896,9 +896,9 @@ LaGrange interpolation. *****************************************************************************/ -static float interpolate1DF32(float *domain, - float *range, +static psF32 interpolate1DF32(psF32 *domain, + psF32 *range, psS32 n, psS32 order, - float x) + psF32 x) { psS32 binNum; @@ -941,7 +941,7 @@ evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] *****************************************************************************/ -float psGaussian(float x, float mean, float sigma, psBool normal) -{ - float tmp = 1.0; +psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal) +{ + psF32 tmp = 1.0; psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, @@ -967,5 +967,5 @@ XXX: There is no way to seed the random generator. *****************************************************************************/ -psVector* p_psGaussianDev(float mean, float sigma, psS32 Npts) +psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) { PS_INT_CHECK_NON_NEGATIVE(Npts, NULL); @@ -1007,7 +1007,7 @@ newPoly->type = type; newPoly->n = n; - newPoly->coeff = (float *)psAlloc(n * sizeof(float)); - newPoly->coeffErr = (float *)psAlloc(n * sizeof(float)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); + 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; @@ -1036,11 +1036,11 @@ newPoly->nY = nY; - newPoly->coeff = (float **)psAlloc(nX * sizeof(float *)); - newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + 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] = (float *)psAlloc(nY * sizeof(float)); - newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + 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++) { @@ -1075,15 +1075,15 @@ newPoly->nZ = nZ; - newPoly->coeff = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + 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] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + 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] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + 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)); } } @@ -1124,19 +1124,19 @@ newPoly->nZ = nZ; - newPoly->coeff = (float ****)psAlloc(nW * sizeof(float ***)); - newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***)); - newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); + newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); + newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); } } @@ -1157,5 +1157,5 @@ } -float psPolynomial1DEval(const psPolynomial1D* myPoly, float x) +psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1195,5 +1195,5 @@ } -float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, float y) +psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1250,5 +1250,5 @@ } -float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, float z) +psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1316,5 +1316,5 @@ } -float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, float z) +psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1406,7 +1406,7 @@ newPoly->type = type; newPoly->n = n; - newPoly->coeff = (double *)psAlloc(n * sizeof(double)); - newPoly->coeffErr = (double *)psAlloc(n * sizeof(double)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); + 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; @@ -1435,11 +1435,11 @@ newPoly->nY = nY; - newPoly->coeff = (double **)psAlloc(nX * sizeof(double *)); - newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + 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] = (double *)psAlloc(nY * sizeof(double)); - newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + 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++) { @@ -1474,15 +1474,15 @@ newPoly->nZ = nZ; - newPoly->coeff = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + 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] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + 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] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + 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)); } } @@ -1523,19 +1523,19 @@ newPoly->nZ = nZ; - newPoly->coeff = (double ****)psAlloc(nW * sizeof(double ***)); - newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***)); - newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); + newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); + newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); + newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); + newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); } } @@ -1557,5 +1557,5 @@ -double psDPolynomial1DEval(const psDPolynomial1D* myPoly, double x) +psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1598,7 +1598,7 @@ -double psDPolynomial2DEval(const psDPolynomial2D* myPoly, - double x, - double y) +psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly, + psF64 x, + psF64 y) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1655,8 +1655,8 @@ -double psDPolynomial3DEval(const psDPolynomial3D* myPoly, - double x, - double y, - double z) +psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly, + psF64 x, + psF64 y, + psF64 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1724,9 +1724,9 @@ } -double psDPolynomial4DEval(const psDPolynomial4D* myPoly, - double w, - double x, - double y, - double z) +psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly, + psF64 w, + psF64 x, + psF64 y, + psF64 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1810,6 +1810,6 @@ // psS32 n; // psPolynomial1D **spline; -// float *p_psDeriv2; -// float *domains; +// psF32 *p_psDeriv2; +// psF32 *domains; //} psSpline1D; @@ -1822,6 +1822,6 @@ psSpline1D *psSpline1DAlloc(psS32 numSplines, psS32 order, - float min, - float max) + psF32 min, + psF32 max) { PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL); @@ -1831,6 +1831,6 @@ psSpline1D *tmp = NULL; psS32 i; - float tmpDomain; - float width; + psF32 tmpDomain; + psF32 width; tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); @@ -1845,6 +1845,6 @@ tmp->p_psDeriv2 = NULL; - tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); - width = (max - min) / ((float) numSplines); + tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32)); + width = (max - min) / ((psF32) numSplines); (tmp->domains)[0] = min; @@ -1884,5 +1884,5 @@ } - tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); + tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32)); for (i=0;in;i++) { @@ -1977,5 +1977,5 @@ { PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3); - char* strType; + psS8* strType; switch (x->type.type) { @@ -2074,10 +2074,10 @@ domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); - psScalar *tmpScalar = psScalarAlloc((double) + psScalar *tmpScalar = psScalarAlloc((psF64) interpolate1DF32(domain32->data.F32, range32->data.F32, domain32->n, order, - (float) x->data.F64), PS_TYPE_F64); + (psF32) x->data.F64), PS_TYPE_F64); psFree(range32); psFree(domain32); @@ -2089,5 +2089,5 @@ } else { - char* strType; + psS8* strType; PS_TYPE_NAME(strType,x->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, @@ -2115,6 +2115,6 @@ the spline fit functions require F32 and F64. *****************************************************************************/ -float psSpline1DEval( - float x, +psF32 psSpline1DEval( + psF32 x, const psSpline1D *spline ) @@ -2172,10 +2172,10 @@ for (i=0;in;i++) { tmpVector->data.F32[i] = psSpline1DEval( - (float) x->data.F64[i], + (psF32) x->data.F64[i], spline ); } } else { - char* strType; + psS8* strType; PS_TYPE_NAME(strType,x->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, Index: /trunk/psLib/src/dataManip/psFunctions.h =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.h (revision 2787) +++ /trunk/psLib/src/dataManip/psFunctions.h (revision 2788) @@ -12,6 +12,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-12-02 21:12:51 $ +* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -39,10 +39,10 @@ * \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 stddev, ///< Standard deviation for the Gaussian + * @return psF32 value on the gaussian curve given the input parameters + */ +psF32 psGaussian( + psF32 x, ///< Value at which to evaluate + psF32 mean, ///< Mean for the Gaussian + psF32 stddev, ///< Standard deviation for the Gaussian psBool normal ///< Indicates whether result should be normalized ); @@ -55,6 +55,6 @@ */ psVector* p_psGaussianDev( - float mean, ///< The mean of the Gaussian - float sigma, ///< The sigma of the Gaussian + psF32 mean, ///< The mean of the Gaussian + psF32 sigma, ///< The sigma of the Gaussian psS32 Npts ///< The size of the vector ); @@ -70,7 +70,7 @@ psPolynomialType type; ///< Polynomial type psS32 n; ///< Number of terms - float *coeff; ///< Coefficients - float *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + psF32 *coeff; ///< Coefficients + psF32 *coeffErr; ///< Error in coefficients + psU8 *mask; ///< Coefficient mask } psPolynomial1D; @@ -82,7 +82,7 @@ psS32 nX; ///< Number of terms in x psS32 nY; ///< Number of terms in y - float **coeff; ///< Coefficients - float **coeffErr; ///< Error in coefficients - char **mask; ///< Coefficients mask + psF32 **coeff; ///< Coefficients + psF32 **coeffErr; ///< Error in coefficients + psU8 **mask; ///< Coefficients mask } psPolynomial2D; @@ -95,7 +95,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - float ***coeff; ///< Coefficients - float ***coeffErr; ///< Error in coefficients - char ***mask; ///< Coefficients mask + psF32 ***coeff; ///< Coefficients + psF32 ***coeffErr; ///< Error in coefficients + psU8 ***mask; ///< Coefficients mask } psPolynomial3D; @@ -109,7 +109,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - float ****coeff; ///< Coefficients - float ****coeffErr; ///< Error in coefficients - char ****mask; ///< Coefficients mask + psF32 ****coeff; ///< Coefficients + psF32 ****coeffErr; ///< Error in coefficients + psU8 ****mask; ///< Coefficients mask } psPolynomial4D; @@ -160,42 +160,42 @@ /** Evaluates a 1-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial1DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial1DEval( const psPolynomial1D* myPoly, ///< Coefficients for the polynomial - float x ///< location at which to evaluate + psF32 x ///< location at which to evaluate ); /** Evaluates a 2-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial2DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial2DEval( const psPolynomial2D* myPoly, ///< Coefficients for the polynomial - float x, ///< x location at which to evaluate - float y ///< y location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y ///< y location at which to evaluate ); /** Evaluates a 3-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial3DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial3DEval( const psPolynomial3D* myPoly, ///< Coefficients for the polynomial - float x, ///< x location at which to evaluate - float y, ///< y location at which to evaluate - float z ///< z location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y, ///< y location at which to evaluate + psF32 z ///< z location at which to evaluate ); /** Evaluates a 4-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial4DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial4DEval( const psPolynomial4D* myPoly, ///< Coefficients for the polynomial - float w, ///< w location at which to evaluate - float x, ///< x location at which to evaluate - float y, ///< y location at which to evaluate - float z ///< z location at which to evaluate + psF32 w, ///< w location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y, ///< y location at which to evaluate + psF32 z ///< z location at which to evaluate ); @@ -235,7 +235,7 @@ psPolynomialType type; ///< Polynomial type psS32 n; ///< Number of terms - double *coeff; ///< Coefficients - double *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + psF64 *coeff; ///< Coefficients + psF64 *coeffErr; ///< Error in coefficients + psU8 *mask; ///< Coefficient mask } psDPolynomial1D; @@ -247,7 +247,7 @@ psS32 nX; ///< Number of terms in x psS32 nY; ///< Number of terms in y - double **coeff; ///< Coefficients - double **coeffErr; ///< Error in coefficients - char **mask; ///< Coefficients mask + psF64 **coeff; ///< Coefficients + psF64 **coeffErr; ///< Error in coefficients + psU8 **mask; ///< Coefficients mask } psDPolynomial2D; @@ -260,7 +260,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - double ***coeff; ///< Coefficients - double ***coeffErr; ///< Error in coefficients - char ***mask; ///< Coefficient mask + psF64 ***coeff; ///< Coefficients + psF64 ***coeffErr; ///< Error in coefficients + psU8 ***mask; ///< Coefficient mask } psDPolynomial3D; @@ -274,7 +274,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - double ****coeff; ///< Coefficients - double ****coeffErr; ///< Error in coefficients - char ****mask; ///< Coefficients mask + psF64 ****coeff; ///< Coefficients + psF64 ****coeffErr; ///< Error in coefficients + psU8 ****mask; ///< Coefficients mask } psDPolynomial4D; @@ -324,42 +324,42 @@ /** Evaluates a double-precision 1-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -double psDPolynomial1DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial1DEval( const psDPolynomial1D* myPoly, ///< Coefficients for the polynomial - double x ///< Value at which to evaluate + psF64 x ///< Value at which to evaluate ); /** Evaluates a double-precision 2-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -double psDPolynomial2DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial2DEval( const psDPolynomial2D* myPoly, ///< Coefficients for the polynomial - double x, ///< Value x at which to evaluate - double y ///< Value y at which to evaluate + 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 float result of polynomial at given location - */ -double psDPolynomial3DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial3DEval( const psDPolynomial3D* myPoly, ///< Coefficients for the polynomial - double x, ///< Value x at which to evaluate - double y, ///< Value y at which to evaluate - double z ///< Value z at which to evaluate + 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 float result of polynomial at given location - */ -double psDPolynomial4DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial4DEval( const psDPolynomial4D* myPoly, ///< Coefficients for the polynomial - double w, ///< Value w at which to evaluate - double x, ///< Value x at which to evaluate - double y, ///< Value y at which to evaluate - double z ///< Value z at which to evaluate + psF64 w, ///< Value w at which to evaluate + psF64 x, ///< Value x at which to evaluate + psF64 y, ///< Value y at which to evaluate + psF64 z ///< Value z at which to evaluate ); @@ -396,6 +396,6 @@ psS32 n; ///< The number of spline polynomials psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials - float *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. - float *domains; ///< The boundaries between each spline piece. Size is n+1. + psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. + psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. } psSpline1D; @@ -403,12 +403,12 @@ psSpline1D *psSpline1DAlloc(psS32 n, psS32 order, - float min, - float max); + psF32 min, + psF32 max); psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, psS32 order); -float psSpline1DEval( - float x, +psF32 psSpline1DEval( + psF32 x, const psSpline1D *spline ); @@ -427,8 +427,8 @@ psScalar *x); -float p_psNRSpline1DEval(psSpline1D *spline, +psF32 p_psNRSpline1DEval(psSpline1D *spline, const psVector* restrict x, const psVector* restrict y, - float X); + psF32 X); /* \} */// End of MathGroup Functions Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 2787) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 2788) @@ -9,14 +9,10 @@ * @author GLG, MHPCC * - * @version $Revision: 1.95 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-17 00:18:31 $ + * @version $Revision: 1.96 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * * XXX: must follow coding name standards on local functions. - * - * XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via - * LMM) is not addressed here. We are waiting for subsequent SDRs - * which will redefine the LMM functions. * */ @@ -64,10 +60,10 @@ XXX: Use a static vector. *****************************************************************************/ -void psBuildSums1D(double x, +void psBuildSums1D(psF64 x, psS32 polyOrder, psVector* sums) { psS32 i = 0; - double xSum = 0.0; + psF64 xSum = 0.0; if (sums == NULL) { @@ -100,5 +96,5 @@ XXX: do an F64 version? *****************************************************************************/ -float *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) +psF32 *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) const psVector* restrict y) ///< Coordinates { @@ -108,12 +104,12 @@ psS32 i; psS32 k; - float sig; - float p; + psF32 sig; + psF32 p; psS32 n = y->n; - float *u = (float *) psAlloc(n * sizeof(float)); - float *derivs2 = (float *) psAlloc(n * sizeof(float)); - float *X = (float *) & (x->data.F32[0]); - float *Y = (float *) & (y->data.F32[0]); - float qn; + psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *X = (psF32 *) & (x->data.F32[0]); + psF32 *Y = (psF32 *) & (y->data.F32[0]); + psF32 qn; // XXX: The second derivatives at the endpoints, undefined in the SDR, @@ -168,8 +164,8 @@ *****************************************************************************/ /* -float p_psNRSpline1DEval(psSpline1D *spline, +psF32 p_psNRSpline1DEval(psSpline1D *spline, const psVector* restrict x, const psVector* restrict y, - float X) + psF32 X) { PS_PTR_CHECK_NULL(spline, NAN); @@ -185,10 +181,10 @@ psS32 klo; psS32 khi; - float H; - float A; - float B; - float C; - float D; - float Y; + psF32 H; + psF32 A; + psF32 B; + psF32 C; + psF32 D; + psF32 Y; n = spline->n; @@ -258,8 +254,8 @@ "---- psVectorFitSpline1D() begin ----\n"); psS32 numSplines = (y->n)-1; - float tmp; - float H; + psF32 tmp; + psF32 H; psS32 i; - float slope; + psF32 slope; psVector *x32 = NULL; psVector *y32 = NULL; @@ -421,8 +417,8 @@ psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2Gauss1D() begin ----\n"); - float x; - int i; - float mean = params->data.F32[0]; - float stdev = params->data.F32[1]; + psF32 x; + psS32 i; + psF32 mean = params->data.F32[0]; + psF32 stdev = params->data.F32[1]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -443,5 +439,5 @@ for (i=0;in;i++) { x = ((psVector *) (coords->data[i]))->data.F32[0]; - float tmp = (x - mean) * psGaussian(x, mean, stdev, false); + psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false); deriv->data.F32[i][0] = tmp / (stdev * stdev); tmp = (x - mean) * (x - mean) * @@ -491,10 +487,10 @@ PS_PTR_CHECK_NULL(params, NULL); - double normalization = params->data.F32[0]; - double x0 = params->data.F32[1]; - double y0 = params->data.F32[2]; - double sigmaX = params->data.F32[3]; - double sigmaY = params->data.F32[4]; - double theta = params->data.F32[5]; + psF64 normalization = params->data.F32[0]; + psF64 x0 = params->data.F32[1]; + psF64 y0 = params->data.F32[2]; + psF64 sigmaX = params->data.F32[3]; + psF64 sigmaY = params->data.F32[4]; + psF64 theta = params->data.F32[5]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -502,5 +498,5 @@ deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); } else { - // XXX: Check size of derivative + PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL); } @@ -508,14 +504,14 @@ "---- psMinimizeLMChi2Gauss2D() begin ----\n"); - for (int i=0;in;i++) { - double x = ((psVector *) coords->data[i])->data.F32[0]; - double y = ((psVector *) coords->data[i])->data.F32[0]; - - double u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); - double v = (x-x0)*cos(theta) + (y-y0)*sin(theta); - - double flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + - v*v/(2.0 * sigmaY * sigmaY)))/ - (2.0 * M_PI * sigmaX * sigmaY); + for (psS32 i=0;in;i++) { + psF64 x = ((psVector *) coords->data[i])->data.F32[0]; + psF64 y = ((psVector *) coords->data[i])->data.F32[0]; + + psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); + psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta); + + psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + + v*v/(2.0 * sigmaY * sigmaY)))/ + (2.0 * M_PI * sigmaX * sigmaY); out->data.F32[i] = flux; @@ -584,6 +580,4 @@ PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL); } - - // XXX: Generate code that modifies covar matrix if not-NULL. if (covar != NULL) { PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL); @@ -610,17 +604,10 @@ psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); - - // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *)); - // for (i=0;idata.F64[j][k] = - (double) ((1.0 + lamda) * alpha->data.F32[j][k]); + (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]); } else { - A->data.F64[j][k] = (double) alpha->data.F32[j][k]; + A->data.F64[j][k] = (psF64) alpha->data.F32[j][k]; } } @@ -741,5 +728,5 @@ } else { newParams->data.F32[i] = params->data.F32[i] - - (float) paramDeltasF64->data.F64[i]; + (psF32) paramDeltasF64->data.F64[i]; } } @@ -777,5 +764,5 @@ // We already masked params. for (i=0;idata.F32[i] = (float) newParams->data.F32[i]; + params->data.F32[i] = (psF32) newParams->data.F32[i]; } lamda*= 0.1; @@ -826,6 +813,6 @@ psS32 k; psS32 n = x->n; - double fac; - double sum; + psF64 fac; + psF64 sum; PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64); psScalar *fScalar; @@ -836,8 +823,8 @@ // variable declarations. I retain them here to maintain coherence // with the NR code. - double min = -1.0; - double max = 1.0; - double bma = 0.5 * (max-min); // 1 - double bpa = 0.5 * (max+min); // 0 + psF64 min = -1.0; + psF64 max = 1.0; + psF64 bma = 0.5 * (max-min); // 1 + psF64 bpa = 0.5 * (max+min); // 0 // In this loop, we first calculate the values of X for which the @@ -851,6 +838,6 @@ for (psS32 i=0;in;j++) { @@ -879,5 +866,5 @@ for (k=0;kdata.F64[k] * - cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n)); + cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n)); } @@ -1075,5 +1062,5 @@ *****************************************************************************/ psMinimization *psMinimizationAlloc(psS32 maxIter, - float tol) + psF32 tol) { PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL); @@ -1092,5 +1079,5 @@ // LINE to it. We assume BASEMASK is non-null. #define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if (BASEMASK->data.U8[i] == 0) { \ OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ @@ -1102,5 +1089,5 @@ #define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \ BOOL_VAR = true; \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \ BOOL_VAR = false; \ @@ -1111,5 +1098,5 @@ #define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \ BOOL_VAR = true; \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \ BOOL_VAR = false; \ @@ -1152,17 +1139,17 @@ psMinimizePowellFunc func) { - float a = 0.0; - float b = 0.0; - float c = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; psS32 iter = 100; - float aDir = 0.0; - float cDir = 0.0; - float new_aDir = 0.0; - float new_cDir = 0.0; + psF32 aDir = 0.0; + psF32 cDir = 0.0; + psF32 new_aDir = 0.0; + psF32 new_cDir = 0.0; psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); - float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; + psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; psVector *tmp = NULL; psBool boolLineIsNull = true; @@ -1301,15 +1288,15 @@ psMinimizePowellFunc func) { - float a = 0.0; - float b = 0.0; - float c = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; psS32 iter = 0; PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32); psBool boolLineIsNull = true; - float prevMin = 0.0; - int countMin = 0; + psF32 prevMin = 0.0; + psS32 countMin = 0; psTrace(".psLib.dataManip.p_psDetermineBracket", 4, @@ -1403,5 +1390,5 @@ *****************************************************************************/ #define PS_LINEMIN_MAX_ITERATIONS 30 -float p_psLineMin(psMinimization *min, +psF32 p_psLineMin(psMinimization *min, psVector *params, psVector *line, @@ -1423,13 +1410,13 @@ PS_PTR_CHECK_NULL(func, NAN); psVector *bracket; - float a = 0.0; - float b = 0.0; - float c = 0.0; - float n = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; - float fn = 0.0; - float mul = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 n = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; + psF32 fn = 0.0; + psF32 mul = 0.0; PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32); PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32); @@ -1544,5 +1531,5 @@ This routine must minimize a possibly multi-dimensional function. The function to be minimized "func" is: - float func(psVector *params, psArray *coords) + psF32 func(psVector *params, psArray *coords) The "params" are the parameters of the function which are varied. The data points at which the function is varied are in the argument "coords" which is @@ -1578,10 +1565,10 @@ psVector *myParamMask = NULL; psMinimization dummyMin; - float mul = 0.0; - float baseFuncVal = 0.0; - float currFuncVal = 0.0; + psF32 mul = 0.0; + psF32 baseFuncVal = 0.0; + psF32 currFuncVal = 0.0; psS32 biggestIter = 0; - float biggestDiff = 0.0; - int iterationNumber = 0; + psF32 biggestDiff = 0.0; + psS32 iterationNumber = 0; psTrace(".psLib.dataManip.psMinimizePowell", 4, @@ -1711,9 +1698,9 @@ } } - float fqp = func(pQP, coords); - float term1 = (baseFuncVal - currFuncVal) - biggestDiff; + psF32 fqp = func(pQP, coords); + psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff; term1*= term1; term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp); - float term2 = baseFuncVal - fqp; + psF32 term2 = baseFuncVal - fqp; term2*= term2 * biggestDiff; if (term1 < term2) { @@ -1760,8 +1747,8 @@ PS_PTR_CHECK_NULL(params, NULL); - float x; - int i; - float mean = params->data.F32[0]; - float stdev = params->data.F32[1]; + psF32 x; + psS32 i; + psF32 mean = params->data.F32[0]; + psF32 stdev = params->data.F32[1]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -1783,5 +1770,5 @@ XXX: This is F32 only *****************************************************************************/ -float myPowellChi2Func(const psVector *params, +psF32 myPowellChi2Func(const psVector *params, const psArray *coords) { @@ -1794,6 +1781,6 @@ PS_PTR_CHECK_NULL(coords, NAN); - float chi2 = 0.0; - float d; + psF32 chi2 = 0.0; + psF32 d; psS32 i; psVector *tmp; Index: /trunk/psLib/src/dataManip/psMinimize.h =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.h (revision 2787) +++ /trunk/psLib/src/dataManip/psMinimize.h (revision 2788) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-01 23:57:08 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -40,13 +40,13 @@ { psS32 maxIter; ///< Convergence limit - float tol; ///< Error Tolerance - float value; ///< Value of function at minimum + psF32 tol; ///< Error Tolerance + psF32 value; ///< Value of function at minimum psS32 iter; ///< Number of iterations required - float lastDelta; ///< The last difference for the fit + psF32 lastDelta; ///< The last difference for the fit } psMinimization; psMinimization *psMinimizationAlloc(psS32 maxIter, - float tol); + psF32 tol); /** Derive a polynomial fit. @@ -91,5 +91,5 @@ typedef -float (*psMinimizePowellFunc)(const psVector *params, +psF32 (*psMinimizePowellFunc)(const psVector *params, const psArray *coords); Index: /trunk/psLib/src/dataManip/psStats.c =================================================================== --- /trunk/psLib/src/dataManip/psStats.c (revision 2787) +++ /trunk/psLib/src/dataManip/psStats.c (revision 2788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.107 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-22 00:54:28 $ + * @version $Revision: 1.108 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -70,5 +70,5 @@ /*****************************************************************************/ -psBool p_psGetStatValue(const psStats* stats, double *value) +psBool p_psGetStatValue(const psStats* stats, psF64 *value) { @@ -145,13 +145,13 @@ this routine sets stats->sampleMean to NAN. *****************************************************************************/ -int p_psVectorSampleMean(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorSampleMean(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float mean = 0.0; // The mean + psF32 mean = 0.0; // The mean psS32 count = 0; // # of points in this mean @@ -171,5 +171,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -185,5 +185,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -199,5 +199,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -207,5 +207,5 @@ mean += myVector->data.F32[i]; } - mean /= (float)myVector->n; + mean /= (psF32)myVector->n; } } @@ -283,11 +283,11 @@ this routine sets stats->max to NAN. *****************************************************************************/ -int p_psVectorMax(const psVector* restrict myVector, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorMax(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float max = -PS_MAX_F32; // The calculated maximum + psF32 max = -PS_MAX_F32; // The calculated maximum psS32 empty = true; // Does this vector have valid elements? @@ -348,11 +348,11 @@ this routine sets stats->min to NAN. *****************************************************************************/ -int p_psVectorMin(const psVector* restrict myVector, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorMin(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float min = PS_MAX_F32; // The calculated maximum + psF32 min = PS_MAX_F32; // The calculated maximum psS32 empty = true; // Does this vector have valid elements? @@ -612,5 +612,5 @@ *****************************************************************************/ psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, - float sigma) + psF32 sigma) { PS_PTR_CHECK_NULL(robustHistogram, NULL); @@ -619,6 +619,6 @@ psS32 i = 0; // Loop index variable psS32 j = 0; // Loop index variable - float iMid; - float jMid; + psF32 iMid; + psF32 jMid; psS32 numBins = robustHistogram->nums->n; psS32 numBounds = robustHistogram->bounds->n; @@ -626,6 +626,6 @@ psS32 jMin = 0; psS32 jMax = 0; - float firstBound = robustHistogram->bounds->data.F32[0]; - float lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + psF32 firstBound = robustHistogram->bounds->data.F32[0]; + psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; psScalar x; @@ -795,9 +795,9 @@ psS32 i = 0; // Loop index variable psS32 countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable + psF32 countFloat = 0.0; // # of data points being used + psF32 mean = 0.0; // The mean + psF32 diff = 0.0; // Used in calculating stdev + psF32 sumSquares = 0.0; // temporary variable + psF32 sumDiffs = 0.0; // temporary variable // This procedure requires the mean. If it has not been already @@ -865,5 +865,5 @@ psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", countInt); } else { - countFloat = (float)countInt; + countFloat = (psF32)countInt; stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); } @@ -889,9 +889,9 @@ psS32 i = 0; // Loop index variable psS32 countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable + psF32 countFloat = 0.0; // # of data points being used + psF32 mean = 0.0; // The mean + psF32 diff = 0.0; // Used in calculating stdev + psF32 sumSquares = 0.0; // temporary variable + psF32 sumDiffs = 0.0; // temporary variable // psF32 sum1; // psF32 sum2; @@ -982,5 +982,5 @@ stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor)); } else { - countFloat = (float)countInt; + countFloat = (psF32)countInt; stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); @@ -1003,16 +1003,16 @@ -2: warning *****************************************************************************/ -int p_psVectorClippedStats(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorClippedStats(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable psS32 j = 0; // Loop index variable - float clippedMean = 0.0; // self-explanatory - float clippedStdev = 0.0; // self-explanatory - float oldStanMean = 0.0; // Temporary variable - float oldStanStdev = 0.0; // Temporary variable + psF32 clippedMean = 0.0; // self-explanatory + psF32 clippedStdev = 0.0; // self-explanatory + psF32 oldStanMean = 0.0; // Temporary variable + psF32 oldStanStdev = 0.0; // Temporary variable psVector* tmpMask = NULL; // Temporary vector @@ -1240,8 +1240,8 @@ XXX: Create a 2nd-order polynomial version and solve for X analytically. *****************************************************************************/ -float p_ps1DPolyMedian(psPolynomial1D* myPoly, - float rangeLow, - float rangeHigh, - float getThisValue) +psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly, + psF32 rangeLow, + psF32 rangeHigh, + psF32 getThisValue) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1250,9 +1250,9 @@ // falls within the range of y-values of the polynomial "myPoly" in the // specified x-range (rangeLow:rangeHigh). - float fLo = psPolynomial1DEval( + psF32 fLo = psPolynomial1DEval( myPoly, rangeLow ); - float fHi = psPolynomial1DEval( + psF32 fHi = psPolynomial1DEval( myPoly, rangeHigh @@ -1266,7 +1266,7 @@ psS32 numIterations = 0; - float midpoint = 0.0; - float oldMidpoint = 1.0; - float f = 0.0; + psF32 midpoint = 0.0; + psF32 oldMidpoint = 1.0; + psF32 f = 0.0; while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) { @@ -1309,8 +1309,8 @@ tests to ensure that binNum is within acceptable ranges for both vectors. *****************************************************************************/ -float fitQuadraticSearchForYThenReturnX(psVector *xVec, +psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, psVector *yVec, psS32 binNum, - float yVal) + psF32 yVal) { PS_VECTOR_CHECK_NULL(xVec, NAN); @@ -1331,11 +1331,11 @@ psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - float tmpFloat; + psF32 tmpFloat; if ((binNum > 0) && (binNum < (yVec->n - 2))) { // The general case. We have all three points. - x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); - x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); - x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); + x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); y->data.F64[0] = yVec->data.F32[binNum - 1]; y->data.F64[1] = yVec->data.F32[binNum]; @@ -1346,5 +1346,5 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, - (double)yVal,y->data.F64[2],y->data.F64[0]); + (psF64)yVal,y->data.F64[2],y->data.F64[0]); } yErr->data.F64[0] = 1.0; @@ -1425,13 +1425,13 @@ XXX: Check for errors in psLib routines that we call. *****************************************************************************/ -int p_psVectorRobustStats(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorRobustStats(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psHistogram* robustHistogram = NULL; psVector* robustHistogramVector = NULL; - float binSize = 0.0; // Size of the histogram bins + psF32 binSize = 0.0; // Size of the histogram bins psS32 LQBinNum = -1; // Bin num for lower quartile psS32 UQBinNum = -1; // Bin num for upper quartile @@ -1439,17 +1439,17 @@ psS32 i = 0; // Loop index variable psS32 modeBinNum = 0; - float modeBinCount = 0.0; - float dL = 0.0; + psF32 modeBinCount = 0.0; + psF32 dL = 0.0; psS32 numBins = 0; - float myMean = 0.0; - float myStdev = 0.0; - float countFloat = 0.0; - float diff = 0.0; - float sumSquares = 0.0; - float sumDiffs = 0.0; + psF32 myMean = 0.0; + psF32 myStdev = 0.0; + psF32 countFloat = 0.0; + psF32 diff = 0.0; + psF32 sumSquares = 0.0; + psF32 sumDiffs = 0.0; psVector* cumulativeRobustSums = NULL; - float sumRobust = 0.0; - float sumN50 = 0.0; - float sumNfit = 0.0; + psF32 sumRobust = 0.0; + psF32 sumN50 = 0.0; + psF32 sumNfit = 0.0; psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; @@ -1460,5 +1460,5 @@ // that by 10.0; //XXX: add errors - int rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); + psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); if (rc != 0) { psError(PS_ERR_UNEXPECTED_NULL, @@ -1614,5 +1614,5 @@ for (i=0;in;i++) { myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); - ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i; + ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i; y->data.F32[i] = robustHistogramVector->data.F32[i]; } @@ -1770,5 +1770,5 @@ The histogram structure *****************************************************************************/ -psHistogram* psHistogramAlloc(float lower, float upper, psS32 n) +psHistogram* psHistogramAlloc(psF32 lower, psF32 upper, psS32 n) { PS_INT_CHECK_POSITIVE(n, NULL); @@ -1777,5 +1777,5 @@ psS32 i = 0; // Loop index variable psHistogram* newHist = NULL; // The new histogram structure - float binSize = 0.0; // The histogram bin size + psF32 binSize = 0.0; // The histogram bin size // Allocate memory for the new histogram structure. If there are N @@ -1787,10 +1787,10 @@ // Calculate the bounds for each bin. - binSize = (upper - lower) / (float)n; + binSize = (upper - lower) / (psF32)n; // XXX: Is the following necessary? It prevents the max data point // from being in a non-existant bin. binSize += FLT_EPSILON; for (i = 0; i < n + 1; i++) { - newHist->bounds->data.F32[i] = lower + (binSize * (float)i); + newHist->bounds->data.F32[i] = lower + (binSize * (psF32)i); } @@ -1857,4 +1857,17 @@ } +/***************************************************************************** +UpdateHistogramBins(binNum, out, data, error): This routine is to be used when +updating the histogram in the presence of errors in the input data. We treat +the data point as a boxcar PDF and update a range of points surrounding the +histogram bin which contains the point. The width of that boxcar is defined +as 2.35 * error. Inputs: + binNum: the bin number of the data point in the histogram + out: the histogram structure + data: the data point value + error: the error in that data point + +XXX: Must test this. + *****************************************************************************/ psS32 UpdateHistogramBins(psS32 binNum, psHistogram* out, @@ -1863,20 +1876,68 @@ { PS_PTR_CHECK_NULL(out, -1); - PS_INT_CHECK_RANGE(binNum, 0, out->nums->n-1, -2); - /* - psF32 width = 2.35 * error; - psF32 centerBinWidth = out->bounds->data.F32[binNum+1] - out->bounds->data.F32[binNum] - psF32 boxcarCenter = (out->bounds->data.F32[binNum] + out->bounds->data.F32[binNum+1]) / 2.0; - - if (width <= centerBinWidth) { - out->nums->data.F32[binNum]+= 1.0; - } else { - out->nums->data.F32[binNum]+= centerBinWidth / width; - // XXX: walk to the left, adding fractional values. - // XXX: walk to the right, adding fractional values. - - - } - */ + PS_PTR_CHECK_NULL(out->bounds, -1); + PS_PTR_CHECK_NULL(out->nums, -1); + PS_INT_CHECK_RANGE(binNum, 0, ((out->nums->n)-1), -2); + PS_FLOAT_COMPARE(0.0, error, -3); + PS_FLOAT_CHECK_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4); + + psF32 boxcarWidth = 2.35 * error; + psF32 boxcarCenter = (out->bounds->data.F32[binNum] + + out->bounds->data.F32[binNum+1]) / 2.0; + psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); + psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); + psS32 bin; + psS32 boxcarLeftBinNum; + psS32 boxcarRightBinNum; + + // Determine the left endpoint of the boxcar for the PDF. + for (bin=binNum ; bin >= 0 ; bin--) { + if (out->nums->data.F32[bin] <= boxcarLeft) { + boxcarLeftBinNum = bin; + break; + } + } + + // Determine the right endpoint of the boxcar for the PDF. + for (bin=binNum ; bin < out->nums->n ; bin++) { + if (out->nums->data.F32[bin] >= boxcarRight) { + boxcarRightBinNum = bin; + break; + } + } + + // + // If the boxcar fits entirely inside this bin, then simply add 1.0 to the + // bin and return. + // + if (boxcarLeftBinNum == boxcarRightBinNum) { + out->nums->data.F32[binNum]+= 1.0; + return(0); + } + + // + // If we get here, multiple bins must be updated. We handle the left + // endpoint, and right endpoint differently. + // + out->nums->data.F32[boxcarLeftBinNum]+= + (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth; + + // + // Loop through the center bins, if any. + // + for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) { + out->nums->data.F32[bin]+= + (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth; + } + + // + // Handle the right endpoint differently. + // + out->nums->data.F32[boxcarRightBinNum]+= + (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth; + + // + // Return 0 on success. + // return(0); } @@ -1906,5 +1967,8 @@ { PS_PTR_CHECK_NULL(out, NULL); + PS_VECTOR_CHECK_NULL(out->bounds, NULL); PS_VECTOR_CHECK_TYPE(out->bounds, PS_TYPE_F32, NULL); + PS_INT_CHECK_NON_NEGATIVE(out->bounds->n, NULL); + PS_VECTOR_CHECK_NULL(out->nums, NULL); PS_VECTOR_CHECK_TYPE(out->nums, PS_TYPE_F32, NULL); PS_INT_CHECK_NON_NEGATIVE(out->nums->n, NULL); @@ -1920,16 +1984,24 @@ psS32 i = 0; // Loop index variable - float binSize = 0.0; // Histogram bin size + psF32 binSize = 0.0; // Histogram bin size psS32 binNum = 0; // A temporary bin number psS32 numBins = 0; // The total number of bins - psS32 tmp = 0; psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; - psVector* inF32; - psS32 mustFreeTmp = 1; + psVector* inF32 = NULL; + psVector* errorsF32 = NULL; + psS32 mustFreeVectorIn = 1; + psS32 mustFreeVectorErrors = 1; + + // Convert input and errors vectors to F32 if necessary. inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { inF32 = (psVector *) in; - mustFreeTmp = 0; + mustFreeVectorIn = 0; + } + errorsF32 = p_psConvertToF32((psVector *) errors); + if (errorsF32 == NULL) { + errorsF32 = (psVector *) errors; + mustFreeVectorErrors = 0; } @@ -1950,9 +2022,9 @@ binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); - if (errors != NULL) { + if (errorsF32 != NULL) { // XXX: Check return codes. UpdateHistogramBins(binNum, out, inF32->data.F32[i], - errors->data.F32[i]); + errorsF32->data.F32[i]); } else { // XXX: This if-statement really shouldn't be necessary. @@ -1969,16 +2041,16 @@ // correct bin number requires a bit more work. tmpScalar.data.F32 = inF32->data.F32[i]; - tmp = p_psVectorBinDisect(out->bounds, &tmpScalar); - if (tmp < 0) { + binNum = p_psVectorBinDisect(out->bounds, &tmpScalar); + if (binNum < 0) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: psVectorHistogram(): element outside histogram bounds.\n"); } else { - if (errors != NULL) { + if (errorsF32 != NULL) { // XXX: Check return codes. - UpdateHistogramBins(tmp, out, + UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i]); } else { - (out->nums->data.F32[tmp])+= 1.0; + (out->nums->data.F32[binNum])+= 1.0; } } @@ -1988,6 +2060,9 @@ } - if (mustFreeTmp == 1) { + if (mustFreeVectorIn == 1) { psFree(inF32); + } + if (mustFreeVectorErrors == 1) { + psFree(errorsF32); } return (out); @@ -2015,26 +2090,51 @@ tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.S8[i]; + tmp->data.F32[i] = (psF32)in->data.S8[i]; + } + } else if (in->type.type == PS_TYPE_S16) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32) in->data.S16[i]; + } + } else if (in->type.type == PS_TYPE_S32) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.S32[i]; + } + } else if (in->type.type == PS_TYPE_S64) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.S64[i]; + } + } else if (in->type.type == PS_TYPE_U8) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.U8[i]; } } else if (in->type.type == PS_TYPE_U16) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.U16[i]; - } - } else if (in->type.type == PS_TYPE_U8) { + tmp->data.F32[i] = (psF32)in->data.U16[i]; + } + } else if (in->type.type == PS_TYPE_U32) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.U8[i]; + tmp->data.F32[i] = (psF32)in->data.U32[i]; + } + } else if (in->type.type == PS_TYPE_U64) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.U64[i]; } } else if (in->type.type == PS_TYPE_F64) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.F64[i]; + tmp->data.F32[i] = (psF32)in->data.F64[i]; } } else if (in->type.type == PS_TYPE_F32) { // do nothing } else { - char* strType; - PS_TYPE_NAME(strType,in->type.type); + psS8* strType; + PS_TYPE_NAME(strType, in->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED, @@ -2074,6 +2174,6 @@ } - psVector* inF32; - psVector* errorsF32; + psVector* inF32 = NULL; + psVector* errorsF32 = NULL; psS32 mustFreeVectorIn = 1; psS32 mustFreeVectorErrors = 1; @@ -2178,2 +2278,3 @@ return (stats); } + Index: /trunk/psLib/src/dataManip/psStats.h =================================================================== --- /trunk/psLib/src/dataManip/psStats.h (revision 2787) +++ /trunk/psLib/src/dataManip/psStats.h (revision 2788) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:42:07 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -56,25 +56,25 @@ typedef struct { - 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 + 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; ///< psS32 robustNfit; ///< - double clippedMean; ///< Nsigma clipped mean - double clippedStdev; ///< standard deviation after clipping + psF64 clippedMean; ///< Nsigma clipped mean + psF64 clippedStdev; ///< standard deviation after clipping psS32 clippedNvalues; ///< ??? - double clipSigma; ///< Nsigma used for clipping; user input + psF64 clipSigma; ///< Nsigma used for clipping; user input psS32 clipIter; ///< Number of clipping iterations; user input - double min; ///< minimum data value in array - double max; ///< maximum data value in array - double binsize; ///< + psF64 min; ///< minimum data value in array + psF64 max; ///< maximum data value in array + psF64 binsize; ///< psStatsOptions options; ///< bitmask of calculated values } @@ -132,6 +132,6 @@ */ psHistogram* psHistogramAlloc( - float lower, ///< Lower limit for the bins - float upper, ///< Upper limit for the bins + psF32 lower, ///< Lower limit for the bins + psF32 upper, ///< Upper limit for the bins psS32 n ///< Number of bins ); @@ -172,5 +172,5 @@ ///< the statistic struct to operate on - double *value + psF64 *value ///< if return is true, this is set to the specified statistic value by stats->options ); Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 2787) +++ /trunk/psLib/src/math/psConstants.h (revision 2788) @@ -6,8 +6,14 @@ * @author GLG, MHPCC * - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:42:07 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + * XXX: Add parenthesis around all arguments so that these macros can be + * called with complex expressions. + * + * + * */ Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 2787) +++ /trunk/psLib/src/math/psMinimize.c (revision 2788) @@ -9,14 +9,10 @@ * @author GLG, MHPCC * - * @version $Revision: 1.95 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-17 00:18:31 $ + * @version $Revision: 1.96 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * * XXX: must follow coding name standards on local functions. - * - * XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via - * LMM) is not addressed here. We are waiting for subsequent SDRs - * which will redefine the LMM functions. * */ @@ -64,10 +60,10 @@ XXX: Use a static vector. *****************************************************************************/ -void psBuildSums1D(double x, +void psBuildSums1D(psF64 x, psS32 polyOrder, psVector* sums) { psS32 i = 0; - double xSum = 0.0; + psF64 xSum = 0.0; if (sums == NULL) { @@ -100,5 +96,5 @@ XXX: do an F64 version? *****************************************************************************/ -float *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) +psF32 *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) const psVector* restrict y) ///< Coordinates { @@ -108,12 +104,12 @@ psS32 i; psS32 k; - float sig; - float p; + psF32 sig; + psF32 p; psS32 n = y->n; - float *u = (float *) psAlloc(n * sizeof(float)); - float *derivs2 = (float *) psAlloc(n * sizeof(float)); - float *X = (float *) & (x->data.F32[0]); - float *Y = (float *) & (y->data.F32[0]); - float qn; + psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *X = (psF32 *) & (x->data.F32[0]); + psF32 *Y = (psF32 *) & (y->data.F32[0]); + psF32 qn; // XXX: The second derivatives at the endpoints, undefined in the SDR, @@ -168,8 +164,8 @@ *****************************************************************************/ /* -float p_psNRSpline1DEval(psSpline1D *spline, +psF32 p_psNRSpline1DEval(psSpline1D *spline, const psVector* restrict x, const psVector* restrict y, - float X) + psF32 X) { PS_PTR_CHECK_NULL(spline, NAN); @@ -185,10 +181,10 @@ psS32 klo; psS32 khi; - float H; - float A; - float B; - float C; - float D; - float Y; + psF32 H; + psF32 A; + psF32 B; + psF32 C; + psF32 D; + psF32 Y; n = spline->n; @@ -258,8 +254,8 @@ "---- psVectorFitSpline1D() begin ----\n"); psS32 numSplines = (y->n)-1; - float tmp; - float H; + psF32 tmp; + psF32 H; psS32 i; - float slope; + psF32 slope; psVector *x32 = NULL; psVector *y32 = NULL; @@ -421,8 +417,8 @@ psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2Gauss1D() begin ----\n"); - float x; - int i; - float mean = params->data.F32[0]; - float stdev = params->data.F32[1]; + psF32 x; + psS32 i; + psF32 mean = params->data.F32[0]; + psF32 stdev = params->data.F32[1]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -443,5 +439,5 @@ for (i=0;in;i++) { x = ((psVector *) (coords->data[i]))->data.F32[0]; - float tmp = (x - mean) * psGaussian(x, mean, stdev, false); + psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false); deriv->data.F32[i][0] = tmp / (stdev * stdev); tmp = (x - mean) * (x - mean) * @@ -491,10 +487,10 @@ PS_PTR_CHECK_NULL(params, NULL); - double normalization = params->data.F32[0]; - double x0 = params->data.F32[1]; - double y0 = params->data.F32[2]; - double sigmaX = params->data.F32[3]; - double sigmaY = params->data.F32[4]; - double theta = params->data.F32[5]; + psF64 normalization = params->data.F32[0]; + psF64 x0 = params->data.F32[1]; + psF64 y0 = params->data.F32[2]; + psF64 sigmaX = params->data.F32[3]; + psF64 sigmaY = params->data.F32[4]; + psF64 theta = params->data.F32[5]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -502,5 +498,5 @@ deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); } else { - // XXX: Check size of derivative + PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL); } @@ -508,14 +504,14 @@ "---- psMinimizeLMChi2Gauss2D() begin ----\n"); - for (int i=0;in;i++) { - double x = ((psVector *) coords->data[i])->data.F32[0]; - double y = ((psVector *) coords->data[i])->data.F32[0]; - - double u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); - double v = (x-x0)*cos(theta) + (y-y0)*sin(theta); - - double flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + - v*v/(2.0 * sigmaY * sigmaY)))/ - (2.0 * M_PI * sigmaX * sigmaY); + for (psS32 i=0;in;i++) { + psF64 x = ((psVector *) coords->data[i])->data.F32[0]; + psF64 y = ((psVector *) coords->data[i])->data.F32[0]; + + psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); + psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta); + + psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + + v*v/(2.0 * sigmaY * sigmaY)))/ + (2.0 * M_PI * sigmaX * sigmaY); out->data.F32[i] = flux; @@ -584,6 +580,4 @@ PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL); } - - // XXX: Generate code that modifies covar matrix if not-NULL. if (covar != NULL) { PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL); @@ -610,17 +604,10 @@ psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); - - // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *)); - // for (i=0;idata.F64[j][k] = - (double) ((1.0 + lamda) * alpha->data.F32[j][k]); + (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]); } else { - A->data.F64[j][k] = (double) alpha->data.F32[j][k]; + A->data.F64[j][k] = (psF64) alpha->data.F32[j][k]; } } @@ -741,5 +728,5 @@ } else { newParams->data.F32[i] = params->data.F32[i] - - (float) paramDeltasF64->data.F64[i]; + (psF32) paramDeltasF64->data.F64[i]; } } @@ -777,5 +764,5 @@ // We already masked params. for (i=0;idata.F32[i] = (float) newParams->data.F32[i]; + params->data.F32[i] = (psF32) newParams->data.F32[i]; } lamda*= 0.1; @@ -826,6 +813,6 @@ psS32 k; psS32 n = x->n; - double fac; - double sum; + psF64 fac; + psF64 sum; PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64); psScalar *fScalar; @@ -836,8 +823,8 @@ // variable declarations. I retain them here to maintain coherence // with the NR code. - double min = -1.0; - double max = 1.0; - double bma = 0.5 * (max-min); // 1 - double bpa = 0.5 * (max+min); // 0 + psF64 min = -1.0; + psF64 max = 1.0; + psF64 bma = 0.5 * (max-min); // 1 + psF64 bpa = 0.5 * (max+min); // 0 // In this loop, we first calculate the values of X for which the @@ -851,6 +838,6 @@ for (psS32 i=0;in;j++) { @@ -879,5 +866,5 @@ for (k=0;kdata.F64[k] * - cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n)); + cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n)); } @@ -1075,5 +1062,5 @@ *****************************************************************************/ psMinimization *psMinimizationAlloc(psS32 maxIter, - float tol) + psF32 tol) { PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL); @@ -1092,5 +1079,5 @@ // LINE to it. We assume BASEMASK is non-null. #define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if (BASEMASK->data.U8[i] == 0) { \ OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ @@ -1102,5 +1089,5 @@ #define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \ BOOL_VAR = true; \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \ BOOL_VAR = false; \ @@ -1111,5 +1098,5 @@ #define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \ BOOL_VAR = true; \ -for (int i=0;in;i++) { \ +for (psS32 i=0;in;i++) { \ if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \ BOOL_VAR = false; \ @@ -1152,17 +1139,17 @@ psMinimizePowellFunc func) { - float a = 0.0; - float b = 0.0; - float c = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; psS32 iter = 100; - float aDir = 0.0; - float cDir = 0.0; - float new_aDir = 0.0; - float new_cDir = 0.0; + psF32 aDir = 0.0; + psF32 cDir = 0.0; + psF32 new_aDir = 0.0; + psF32 new_cDir = 0.0; psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); - float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; + psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; psVector *tmp = NULL; psBool boolLineIsNull = true; @@ -1301,15 +1288,15 @@ psMinimizePowellFunc func) { - float a = 0.0; - float b = 0.0; - float c = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; psS32 iter = 0; PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32); psBool boolLineIsNull = true; - float prevMin = 0.0; - int countMin = 0; + psF32 prevMin = 0.0; + psS32 countMin = 0; psTrace(".psLib.dataManip.p_psDetermineBracket", 4, @@ -1403,5 +1390,5 @@ *****************************************************************************/ #define PS_LINEMIN_MAX_ITERATIONS 30 -float p_psLineMin(psMinimization *min, +psF32 p_psLineMin(psMinimization *min, psVector *params, psVector *line, @@ -1423,13 +1410,13 @@ PS_PTR_CHECK_NULL(func, NAN); psVector *bracket; - float a = 0.0; - float b = 0.0; - float c = 0.0; - float n = 0.0; - float fa = 0.0; - float fb = 0.0; - float fc = 0.0; - float fn = 0.0; - float mul = 0.0; + psF32 a = 0.0; + psF32 b = 0.0; + psF32 c = 0.0; + psF32 n = 0.0; + psF32 fa = 0.0; + psF32 fb = 0.0; + psF32 fc = 0.0; + psF32 fn = 0.0; + psF32 mul = 0.0; PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32); PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32); @@ -1544,5 +1531,5 @@ This routine must minimize a possibly multi-dimensional function. The function to be minimized "func" is: - float func(psVector *params, psArray *coords) + psF32 func(psVector *params, psArray *coords) The "params" are the parameters of the function which are varied. The data points at which the function is varied are in the argument "coords" which is @@ -1578,10 +1565,10 @@ psVector *myParamMask = NULL; psMinimization dummyMin; - float mul = 0.0; - float baseFuncVal = 0.0; - float currFuncVal = 0.0; + psF32 mul = 0.0; + psF32 baseFuncVal = 0.0; + psF32 currFuncVal = 0.0; psS32 biggestIter = 0; - float biggestDiff = 0.0; - int iterationNumber = 0; + psF32 biggestDiff = 0.0; + psS32 iterationNumber = 0; psTrace(".psLib.dataManip.psMinimizePowell", 4, @@ -1711,9 +1698,9 @@ } } - float fqp = func(pQP, coords); - float term1 = (baseFuncVal - currFuncVal) - biggestDiff; + psF32 fqp = func(pQP, coords); + psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff; term1*= term1; term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp); - float term2 = baseFuncVal - fqp; + psF32 term2 = baseFuncVal - fqp; term2*= term2 * biggestDiff; if (term1 < term2) { @@ -1760,8 +1747,8 @@ PS_PTR_CHECK_NULL(params, NULL); - float x; - int i; - float mean = params->data.F32[0]; - float stdev = params->data.F32[1]; + psF32 x; + psS32 i; + psF32 mean = params->data.F32[0]; + psF32 stdev = params->data.F32[1]; psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); @@ -1783,5 +1770,5 @@ XXX: This is F32 only *****************************************************************************/ -float myPowellChi2Func(const psVector *params, +psF32 myPowellChi2Func(const psVector *params, const psArray *coords) { @@ -1794,6 +1781,6 @@ PS_PTR_CHECK_NULL(coords, NAN); - float chi2 = 0.0; - float d; + psF32 chi2 = 0.0; + psF32 d; psS32 i; psVector *tmp; Index: /trunk/psLib/src/math/psMinimize.h =================================================================== --- /trunk/psLib/src/math/psMinimize.h (revision 2787) +++ /trunk/psLib/src/math/psMinimize.h (revision 2788) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-01 23:57:08 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -40,13 +40,13 @@ { psS32 maxIter; ///< Convergence limit - float tol; ///< Error Tolerance - float value; ///< Value of function at minimum + psF32 tol; ///< Error Tolerance + psF32 value; ///< Value of function at minimum psS32 iter; ///< Number of iterations required - float lastDelta; ///< The last difference for the fit + psF32 lastDelta; ///< The last difference for the fit } psMinimization; psMinimization *psMinimizationAlloc(psS32 maxIter, - float tol); + psF32 tol); /** Derive a polynomial fit. @@ -91,5 +91,5 @@ typedef -float (*psMinimizePowellFunc)(const psVector *params, +psF32 (*psMinimizePowellFunc)(const psVector *params, const psArray *coords); Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 2787) +++ /trunk/psLib/src/math/psPolynomial.c (revision 2788) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.74 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:42:07 $ + * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -55,5 +55,5 @@ static void dPolynomial4DFree(psDPolynomial4D* myPoly); static void spline1DFree(psSpline1D *tmpSpline); -static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x); +static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); @@ -277,9 +277,9 @@ *****************************************************************************/ -static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly) +static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) { psS32 loop_x = 0; - float polySum = 0.0; - float xSum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, @@ -307,5 +307,5 @@ // XXX: How does the mask vector effect Crenshaw's formula? // XXX: We assume that x is scaled between -1.0 and 1.0; -static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) { // PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0); @@ -314,5 +314,5 @@ psS32 n; psS32 i; - float tmp; + psF32 tmp; n = myPoly->n; @@ -337,5 +337,5 @@ psS32 n; psS32 i; - float tmp; + psF32 tmp; psPolynomial1D **chebPolys = NULL; @@ -355,5 +355,5 @@ } -static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -361,7 +361,7 @@ psS32 loop_x = 0; psS32 loop_y = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -379,5 +379,5 @@ } -static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -386,5 +386,5 @@ psS32 loop_y = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -414,13 +414,13 @@ } -static float ordPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; + psF32 zSum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -442,5 +442,5 @@ } -static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) { psS32 loop_x = 0; @@ -448,5 +448,5 @@ psS32 loop_z = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -483,5 +483,5 @@ } -static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -489,9 +489,9 @@ psS32 loop_y = 0; psS32 loop_z = 0; - float polySum = 0.0; - float wSum = 1.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; + psF32 polySum = 0.0; + psF32 wSum = 1.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; + psF32 zSum = 1.0; for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { @@ -517,5 +517,5 @@ } -static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -524,5 +524,5 @@ psS32 loop_z = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -568,9 +568,9 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) { psS32 loop_x = 0; - double polySum = 0.0; - double xSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; for (loop_x = 0; loop_x < myPoly->n; loop_x++) { @@ -586,10 +586,10 @@ // XXX: You can do this without having to psAlloc() vector d. // XXX: How does the mask vector effect Crenshaw's formula? -static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) { psVector *d; psS32 n; psS32 i; - double tmp; + psF64 tmp; n = myPoly->n; @@ -611,11 +611,11 @@ } -static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -633,10 +633,10 @@ } -static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -667,13 +667,13 @@ } -static double dOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -695,5 +695,5 @@ } -static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) { psS32 loop_x = 0; @@ -701,5 +701,5 @@ psS32 loop_z = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -736,5 +736,5 @@ } -static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -742,9 +742,9 @@ psS32 loop_y = 0; psS32 loop_z = 0; - double polySum = 0.0; - double wSum = 1.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; + psF64 polySum = 0.0; + psF64 wSum = 1.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { @@ -770,5 +770,5 @@ } -static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -777,5 +777,5 @@ psS32 loop_z = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -828,5 +828,5 @@ *****************************************************************************/ #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ -static float fullInterpolate1D##TYPE(ps##TYPE *domain, \ +static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ ps##TYPE *range, \ psS32 n, \ @@ -896,9 +896,9 @@ LaGrange interpolation. *****************************************************************************/ -static float interpolate1DF32(float *domain, - float *range, +static psF32 interpolate1DF32(psF32 *domain, + psF32 *range, psS32 n, psS32 order, - float x) + psF32 x) { psS32 binNum; @@ -941,7 +941,7 @@ evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] *****************************************************************************/ -float psGaussian(float x, float mean, float sigma, psBool normal) -{ - float tmp = 1.0; +psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal) +{ + psF32 tmp = 1.0; psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, @@ -967,5 +967,5 @@ XXX: There is no way to seed the random generator. *****************************************************************************/ -psVector* p_psGaussianDev(float mean, float sigma, psS32 Npts) +psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) { PS_INT_CHECK_NON_NEGATIVE(Npts, NULL); @@ -1007,7 +1007,7 @@ newPoly->type = type; newPoly->n = n; - newPoly->coeff = (float *)psAlloc(n * sizeof(float)); - newPoly->coeffErr = (float *)psAlloc(n * sizeof(float)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); + 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; @@ -1036,11 +1036,11 @@ newPoly->nY = nY; - newPoly->coeff = (float **)psAlloc(nX * sizeof(float *)); - newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + 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] = (float *)psAlloc(nY * sizeof(float)); - newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + 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++) { @@ -1075,15 +1075,15 @@ newPoly->nZ = nZ; - newPoly->coeff = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + 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] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + 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] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + 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)); } } @@ -1124,19 +1124,19 @@ newPoly->nZ = nZ; - newPoly->coeff = (float ****)psAlloc(nW * sizeof(float ***)); - newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***)); - newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); + newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); + newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); } } @@ -1157,5 +1157,5 @@ } -float psPolynomial1DEval(const psPolynomial1D* myPoly, float x) +psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1195,5 +1195,5 @@ } -float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, float y) +psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1250,5 +1250,5 @@ } -float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, float z) +psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1316,5 +1316,5 @@ } -float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, float z) +psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1406,7 +1406,7 @@ newPoly->type = type; newPoly->n = n; - newPoly->coeff = (double *)psAlloc(n * sizeof(double)); - newPoly->coeffErr = (double *)psAlloc(n * sizeof(double)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); + 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; @@ -1435,11 +1435,11 @@ newPoly->nY = nY; - newPoly->coeff = (double **)psAlloc(nX * sizeof(double *)); - newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + 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] = (double *)psAlloc(nY * sizeof(double)); - newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + 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++) { @@ -1474,15 +1474,15 @@ newPoly->nZ = nZ; - newPoly->coeff = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + 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] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + 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] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + 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)); } } @@ -1523,19 +1523,19 @@ newPoly->nZ = nZ; - newPoly->coeff = (double ****)psAlloc(nW * sizeof(double ***)); - newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***)); - newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); + newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); + newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); + newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); + newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); } } @@ -1557,5 +1557,5 @@ -double psDPolynomial1DEval(const psDPolynomial1D* myPoly, double x) +psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1598,7 +1598,7 @@ -double psDPolynomial2DEval(const psDPolynomial2D* myPoly, - double x, - double y) +psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly, + psF64 x, + psF64 y) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1655,8 +1655,8 @@ -double psDPolynomial3DEval(const psDPolynomial3D* myPoly, - double x, - double y, - double z) +psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly, + psF64 x, + psF64 y, + psF64 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1724,9 +1724,9 @@ } -double psDPolynomial4DEval(const psDPolynomial4D* myPoly, - double w, - double x, - double y, - double z) +psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly, + psF64 w, + psF64 x, + psF64 y, + psF64 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1810,6 +1810,6 @@ // psS32 n; // psPolynomial1D **spline; -// float *p_psDeriv2; -// float *domains; +// psF32 *p_psDeriv2; +// psF32 *domains; //} psSpline1D; @@ -1822,6 +1822,6 @@ psSpline1D *psSpline1DAlloc(psS32 numSplines, psS32 order, - float min, - float max) + psF32 min, + psF32 max) { PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL); @@ -1831,6 +1831,6 @@ psSpline1D *tmp = NULL; psS32 i; - float tmpDomain; - float width; + psF32 tmpDomain; + psF32 width; tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); @@ -1845,6 +1845,6 @@ tmp->p_psDeriv2 = NULL; - tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); - width = (max - min) / ((float) numSplines); + tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32)); + width = (max - min) / ((psF32) numSplines); (tmp->domains)[0] = min; @@ -1884,5 +1884,5 @@ } - tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); + tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32)); for (i=0;in;i++) { @@ -1977,5 +1977,5 @@ { PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3); - char* strType; + psS8* strType; switch (x->type.type) { @@ -2074,10 +2074,10 @@ domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); - psScalar *tmpScalar = psScalarAlloc((double) + psScalar *tmpScalar = psScalarAlloc((psF64) interpolate1DF32(domain32->data.F32, range32->data.F32, domain32->n, order, - (float) x->data.F64), PS_TYPE_F64); + (psF32) x->data.F64), PS_TYPE_F64); psFree(range32); psFree(domain32); @@ -2089,5 +2089,5 @@ } else { - char* strType; + psS8* strType; PS_TYPE_NAME(strType,x->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, @@ -2115,6 +2115,6 @@ the spline fit functions require F32 and F64. *****************************************************************************/ -float psSpline1DEval( - float x, +psF32 psSpline1DEval( + psF32 x, const psSpline1D *spline ) @@ -2172,10 +2172,10 @@ for (i=0;in;i++) { tmpVector->data.F32[i] = psSpline1DEval( - (float) x->data.F64[i], + (psF32) x->data.F64[i], spline ); } } else { - char* strType; + psS8* strType; PS_TYPE_NAME(strType,x->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, Index: /trunk/psLib/src/math/psPolynomial.h =================================================================== --- /trunk/psLib/src/math/psPolynomial.h (revision 2787) +++ /trunk/psLib/src/math/psPolynomial.h (revision 2788) @@ -12,6 +12,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-12-02 21:12:51 $ +* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -39,10 +39,10 @@ * \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 stddev, ///< Standard deviation for the Gaussian + * @return psF32 value on the gaussian curve given the input parameters + */ +psF32 psGaussian( + psF32 x, ///< Value at which to evaluate + psF32 mean, ///< Mean for the Gaussian + psF32 stddev, ///< Standard deviation for the Gaussian psBool normal ///< Indicates whether result should be normalized ); @@ -55,6 +55,6 @@ */ psVector* p_psGaussianDev( - float mean, ///< The mean of the Gaussian - float sigma, ///< The sigma of the Gaussian + psF32 mean, ///< The mean of the Gaussian + psF32 sigma, ///< The sigma of the Gaussian psS32 Npts ///< The size of the vector ); @@ -70,7 +70,7 @@ psPolynomialType type; ///< Polynomial type psS32 n; ///< Number of terms - float *coeff; ///< Coefficients - float *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + psF32 *coeff; ///< Coefficients + psF32 *coeffErr; ///< Error in coefficients + psU8 *mask; ///< Coefficient mask } psPolynomial1D; @@ -82,7 +82,7 @@ psS32 nX; ///< Number of terms in x psS32 nY; ///< Number of terms in y - float **coeff; ///< Coefficients - float **coeffErr; ///< Error in coefficients - char **mask; ///< Coefficients mask + psF32 **coeff; ///< Coefficients + psF32 **coeffErr; ///< Error in coefficients + psU8 **mask; ///< Coefficients mask } psPolynomial2D; @@ -95,7 +95,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - float ***coeff; ///< Coefficients - float ***coeffErr; ///< Error in coefficients - char ***mask; ///< Coefficients mask + psF32 ***coeff; ///< Coefficients + psF32 ***coeffErr; ///< Error in coefficients + psU8 ***mask; ///< Coefficients mask } psPolynomial3D; @@ -109,7 +109,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - float ****coeff; ///< Coefficients - float ****coeffErr; ///< Error in coefficients - char ****mask; ///< Coefficients mask + psF32 ****coeff; ///< Coefficients + psF32 ****coeffErr; ///< Error in coefficients + psU8 ****mask; ///< Coefficients mask } psPolynomial4D; @@ -160,42 +160,42 @@ /** Evaluates a 1-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial1DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial1DEval( const psPolynomial1D* myPoly, ///< Coefficients for the polynomial - float x ///< location at which to evaluate + psF32 x ///< location at which to evaluate ); /** Evaluates a 2-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial2DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial2DEval( const psPolynomial2D* myPoly, ///< Coefficients for the polynomial - float x, ///< x location at which to evaluate - float y ///< y location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y ///< y location at which to evaluate ); /** Evaluates a 3-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial3DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial3DEval( const psPolynomial3D* myPoly, ///< Coefficients for the polynomial - float x, ///< x location at which to evaluate - float y, ///< y location at which to evaluate - float z ///< z location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y, ///< y location at which to evaluate + psF32 z ///< z location at which to evaluate ); /** Evaluates a 4-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial4DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial4DEval( const psPolynomial4D* myPoly, ///< Coefficients for the polynomial - float w, ///< w location at which to evaluate - float x, ///< x location at which to evaluate - float y, ///< y location at which to evaluate - float z ///< z location at which to evaluate + psF32 w, ///< w location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y, ///< y location at which to evaluate + psF32 z ///< z location at which to evaluate ); @@ -235,7 +235,7 @@ psPolynomialType type; ///< Polynomial type psS32 n; ///< Number of terms - double *coeff; ///< Coefficients - double *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + psF64 *coeff; ///< Coefficients + psF64 *coeffErr; ///< Error in coefficients + psU8 *mask; ///< Coefficient mask } psDPolynomial1D; @@ -247,7 +247,7 @@ psS32 nX; ///< Number of terms in x psS32 nY; ///< Number of terms in y - double **coeff; ///< Coefficients - double **coeffErr; ///< Error in coefficients - char **mask; ///< Coefficients mask + psF64 **coeff; ///< Coefficients + psF64 **coeffErr; ///< Error in coefficients + psU8 **mask; ///< Coefficients mask } psDPolynomial2D; @@ -260,7 +260,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - double ***coeff; ///< Coefficients - double ***coeffErr; ///< Error in coefficients - char ***mask; ///< Coefficient mask + psF64 ***coeff; ///< Coefficients + psF64 ***coeffErr; ///< Error in coefficients + psU8 ***mask; ///< Coefficient mask } psDPolynomial3D; @@ -274,7 +274,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - double ****coeff; ///< Coefficients - double ****coeffErr; ///< Error in coefficients - char ****mask; ///< Coefficients mask + psF64 ****coeff; ///< Coefficients + psF64 ****coeffErr; ///< Error in coefficients + psU8 ****mask; ///< Coefficients mask } psDPolynomial4D; @@ -324,42 +324,42 @@ /** Evaluates a double-precision 1-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -double psDPolynomial1DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial1DEval( const psDPolynomial1D* myPoly, ///< Coefficients for the polynomial - double x ///< Value at which to evaluate + psF64 x ///< Value at which to evaluate ); /** Evaluates a double-precision 2-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -double psDPolynomial2DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial2DEval( const psDPolynomial2D* myPoly, ///< Coefficients for the polynomial - double x, ///< Value x at which to evaluate - double y ///< Value y at which to evaluate + 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 float result of polynomial at given location - */ -double psDPolynomial3DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial3DEval( const psDPolynomial3D* myPoly, ///< Coefficients for the polynomial - double x, ///< Value x at which to evaluate - double y, ///< Value y at which to evaluate - double z ///< Value z at which to evaluate + 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 float result of polynomial at given location - */ -double psDPolynomial4DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial4DEval( const psDPolynomial4D* myPoly, ///< Coefficients for the polynomial - double w, ///< Value w at which to evaluate - double x, ///< Value x at which to evaluate - double y, ///< Value y at which to evaluate - double z ///< Value z at which to evaluate + psF64 w, ///< Value w at which to evaluate + psF64 x, ///< Value x at which to evaluate + psF64 y, ///< Value y at which to evaluate + psF64 z ///< Value z at which to evaluate ); @@ -396,6 +396,6 @@ psS32 n; ///< The number of spline polynomials psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials - float *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. - float *domains; ///< The boundaries between each spline piece. Size is n+1. + psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. + psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. } psSpline1D; @@ -403,12 +403,12 @@ psSpline1D *psSpline1DAlloc(psS32 n, psS32 order, - float min, - float max); + psF32 min, + psF32 max); psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, psS32 order); -float psSpline1DEval( - float x, +psF32 psSpline1DEval( + psF32 x, const psSpline1D *spline ); @@ -427,8 +427,8 @@ psScalar *x); -float p_psNRSpline1DEval(psSpline1D *spline, +psF32 p_psNRSpline1DEval(psSpline1D *spline, const psVector* restrict x, const psVector* restrict y, - float X); + psF32 X); /* \} */// End of MathGroup Functions Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 2787) +++ /trunk/psLib/src/math/psSpline.c (revision 2788) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.74 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:42:07 $ + * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -55,5 +55,5 @@ static void dPolynomial4DFree(psDPolynomial4D* myPoly); static void spline1DFree(psSpline1D *tmpSpline); -static psS32 vectorBinDisectF32(float *bins,psS32 numBins,float x); +static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); @@ -277,9 +277,9 @@ *****************************************************************************/ -static float ordPolynomial1DEval(float x, const psPolynomial1D* myPoly) +static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) { psS32 loop_x = 0; - float polySum = 0.0; - float xSum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, @@ -307,5 +307,5 @@ // XXX: How does the mask vector effect Crenshaw's formula? // XXX: We assume that x is scaled between -1.0 and 1.0; -static float chebPolynomial1DEval(float x, const psPolynomial1D* myPoly) +static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) { // PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0); @@ -314,5 +314,5 @@ psS32 n; psS32 i; - float tmp; + psF32 tmp; n = myPoly->n; @@ -337,5 +337,5 @@ psS32 n; psS32 i; - float tmp; + psF32 tmp; psPolynomial1D **chebPolys = NULL; @@ -355,5 +355,5 @@ } -static float ordPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -361,7 +361,7 @@ psS32 loop_x = 0; psS32 loop_y = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -379,5 +379,5 @@ } -static float chebPolynomial2DEval(float x, float y, const psPolynomial2D* myPoly) +static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -386,5 +386,5 @@ psS32 loop_y = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -414,13 +414,13 @@ } -static float ordPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; - float polySum = 0.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; + psF32 polySum = 0.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; + psF32 zSum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -442,5 +442,5 @@ } -static float chebPolynomial3DEval(float x, float y, float z, const psPolynomial3D* myPoly) +static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) { psS32 loop_x = 0; @@ -448,5 +448,5 @@ psS32 loop_z = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -483,5 +483,5 @@ } -static float ordPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -489,9 +489,9 @@ psS32 loop_y = 0; psS32 loop_z = 0; - float polySum = 0.0; - float wSum = 1.0; - float xSum = 1.0; - float ySum = 1.0; - float zSum = 1.0; + psF32 polySum = 0.0; + psF32 wSum = 1.0; + psF32 xSum = 1.0; + psF32 ySum = 1.0; + psF32 zSum = 1.0; for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { @@ -517,5 +517,5 @@ } -static float chebPolynomial4DEval(float w, float x, float y, float z, const psPolynomial4D* myPoly) +static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -524,5 +524,5 @@ psS32 loop_z = 0; psS32 i = 0; - float polySum = 0.0; + psF32 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -568,9 +568,9 @@ Polynomial coefficients will be accessed in [w][x][y][z] fashion. *****************************************************************************/ -static double dOrdPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) { psS32 loop_x = 0; - double polySum = 0.0; - double xSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; for (loop_x = 0; loop_x < myPoly->n; loop_x++) { @@ -586,10 +586,10 @@ // XXX: You can do this without having to psAlloc() vector d. // XXX: How does the mask vector effect Crenshaw's formula? -static double dChebPolynomial1DEval(double x, const psDPolynomial1D* myPoly) +static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) { psVector *d; psS32 n; psS32 i; - double tmp; + psF64 tmp; n = myPoly->n; @@ -611,11 +611,11 @@ } -static double dOrdPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -633,10 +633,10 @@ } -static double dChebPolynomial2DEval(double x, double y, const psDPolynomial2D* myPoly) +static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -667,13 +667,13 @@ } -static double dOrdPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) { psS32 loop_x = 0; psS32 loop_y = 0; psS32 loop_z = 0; - double polySum = 0.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; + psF64 polySum = 0.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { @@ -695,5 +695,5 @@ } -static double dChebPolynomial3DEval(double x, double y, double z, const psDPolynomial3D* myPoly) +static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) { psS32 loop_x = 0; @@ -701,5 +701,5 @@ psS32 loop_z = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -736,5 +736,5 @@ } -static double dOrdPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -742,9 +742,9 @@ psS32 loop_y = 0; psS32 loop_z = 0; - double polySum = 0.0; - double wSum = 1.0; - double xSum = 1.0; - double ySum = 1.0; - double zSum = 1.0; + psF64 polySum = 0.0; + psF64 wSum = 1.0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + psF64 zSum = 1.0; for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { @@ -770,5 +770,5 @@ } -static double dChebPolynomial4DEval(double w, double x, double y, double z, const psDPolynomial4D* myPoly) +static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) { psS32 loop_w = 0; @@ -777,5 +777,5 @@ psS32 loop_z = 0; psS32 i = 0; - double polySum = 0.0; + psF64 polySum = 0.0; psPolynomial1D* *chebPolys = NULL; psS32 maxChebyPoly = 0; @@ -828,5 +828,5 @@ *****************************************************************************/ #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ -static float fullInterpolate1D##TYPE(ps##TYPE *domain, \ +static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ ps##TYPE *range, \ psS32 n, \ @@ -896,9 +896,9 @@ LaGrange interpolation. *****************************************************************************/ -static float interpolate1DF32(float *domain, - float *range, +static psF32 interpolate1DF32(psF32 *domain, + psF32 *range, psS32 n, psS32 order, - float x) + psF32 x) { psS32 binNum; @@ -941,7 +941,7 @@ evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] *****************************************************************************/ -float psGaussian(float x, float mean, float sigma, psBool normal) -{ - float tmp = 1.0; +psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal) +{ + psF32 tmp = 1.0; psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, @@ -967,5 +967,5 @@ XXX: There is no way to seed the random generator. *****************************************************************************/ -psVector* p_psGaussianDev(float mean, float sigma, psS32 Npts) +psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) { PS_INT_CHECK_NON_NEGATIVE(Npts, NULL); @@ -1007,7 +1007,7 @@ newPoly->type = type; newPoly->n = n; - newPoly->coeff = (float *)psAlloc(n * sizeof(float)); - newPoly->coeffErr = (float *)psAlloc(n * sizeof(float)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); + 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; @@ -1036,11 +1036,11 @@ newPoly->nY = nY; - newPoly->coeff = (float **)psAlloc(nX * sizeof(float *)); - newPoly->coeffErr = (float **)psAlloc(nX * sizeof(float *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + 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] = (float *)psAlloc(nY * sizeof(float)); - newPoly->coeffErr[x] = (float *)psAlloc(nY * sizeof(float)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + 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++) { @@ -1075,15 +1075,15 @@ newPoly->nZ = nZ; - newPoly->coeff = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->coeffErr = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + 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] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + 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] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + 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)); } } @@ -1124,19 +1124,19 @@ newPoly->nZ = nZ; - newPoly->coeff = (float ****)psAlloc(nW * sizeof(float ***)); - newPoly->coeffErr = (float ****)psAlloc(nW * sizeof(float ***)); - newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); + newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); + newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->coeffErr[w] = (float ***)psAlloc(nX * sizeof(float **)); - newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); + newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->coeffErr[w][x] = (float **)psAlloc(nY * sizeof(float *)); - newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); + newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->coeffErr[w][x][y] = (float *)psAlloc(nZ * sizeof(float)); - newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); + newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); } } @@ -1157,5 +1157,5 @@ } -float psPolynomial1DEval(const psPolynomial1D* myPoly, float x) +psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1195,5 +1195,5 @@ } -float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, float y) +psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1250,5 +1250,5 @@ } -float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, float z) +psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1316,5 +1316,5 @@ } -float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, float z) +psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1406,7 +1406,7 @@ newPoly->type = type; newPoly->n = n; - newPoly->coeff = (double *)psAlloc(n * sizeof(double)); - newPoly->coeffErr = (double *)psAlloc(n * sizeof(double)); - newPoly->mask = (char *)psAlloc(n * sizeof(char)); + 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; @@ -1435,11 +1435,11 @@ newPoly->nY = nY; - newPoly->coeff = (double **)psAlloc(nX * sizeof(double *)); - newPoly->coeffErr = (double **)psAlloc(nX * sizeof(double *)); - newPoly->mask = (char **)psAlloc(nX * sizeof(char *)); + 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] = (double *)psAlloc(nY * sizeof(double)); - newPoly->coeffErr[x] = (double *)psAlloc(nY * sizeof(double)); - newPoly->mask[x] = (char *)psAlloc(nY * sizeof(char)); + 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++) { @@ -1474,15 +1474,15 @@ newPoly->nZ = nZ; - newPoly->coeff = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->coeffErr = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->mask = (char ***)psAlloc(nX * sizeof(char **)); + 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] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->mask[x] = (char **)psAlloc(nY * sizeof(char *)); + 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] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->mask[x][y] = (char *)psAlloc(nZ * sizeof(char)); + 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)); } } @@ -1523,19 +1523,19 @@ newPoly->nZ = nZ; - newPoly->coeff = (double ****)psAlloc(nW * sizeof(double ***)); - newPoly->coeffErr = (double ****)psAlloc(nW * sizeof(double ***)); - newPoly->mask = (char ****)psAlloc(nW * sizeof(char ***)); + newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); + newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); + newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); for (w = 0; w < nW; w++) { - newPoly->coeff[w] = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->coeffErr[w] = (double ***)psAlloc(nX * sizeof(double **)); - newPoly->mask[w] = (char ***)psAlloc(nX * sizeof(char **)); + newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); + newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); for (x = 0; x < nX; x++) { - newPoly->coeff[w][x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->coeffErr[w][x] = (double **)psAlloc(nY * sizeof(double *)); - newPoly->mask[w][x] = (char **)psAlloc(nY * sizeof(char *)); + newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); + newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); for (y = 0; y < nY; y++) { - newPoly->coeff[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->coeffErr[w][x][y] = (double *)psAlloc(nZ * sizeof(double)); - newPoly->mask[w][x][y] = (char *)psAlloc(nZ * sizeof(char)); + newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); + newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); + newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); } } @@ -1557,5 +1557,5 @@ -double psDPolynomial1DEval(const psDPolynomial1D* myPoly, double x) +psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1598,7 +1598,7 @@ -double psDPolynomial2DEval(const psDPolynomial2D* myPoly, - double x, - double y) +psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly, + psF64 x, + psF64 y) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1655,8 +1655,8 @@ -double psDPolynomial3DEval(const psDPolynomial3D* myPoly, - double x, - double y, - double z) +psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly, + psF64 x, + psF64 y, + psF64 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1724,9 +1724,9 @@ } -double psDPolynomial4DEval(const psDPolynomial4D* myPoly, - double w, - double x, - double y, - double z) +psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly, + psF64 w, + psF64 x, + psF64 y, + psF64 z) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1810,6 +1810,6 @@ // psS32 n; // psPolynomial1D **spline; -// float *p_psDeriv2; -// float *domains; +// psF32 *p_psDeriv2; +// psF32 *domains; //} psSpline1D; @@ -1822,6 +1822,6 @@ psSpline1D *psSpline1DAlloc(psS32 numSplines, psS32 order, - float min, - float max) + psF32 min, + psF32 max) { PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL); @@ -1831,6 +1831,6 @@ psSpline1D *tmp = NULL; psS32 i; - float tmpDomain; - float width; + psF32 tmpDomain; + psF32 width; tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); @@ -1845,6 +1845,6 @@ tmp->p_psDeriv2 = NULL; - tmp->domains = (float *) psAlloc((numSplines+1) * sizeof(float)); - width = (max - min) / ((float) numSplines); + tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32)); + width = (max - min) / ((psF32) numSplines); (tmp->domains)[0] = min; @@ -1884,5 +1884,5 @@ } - tmp->domains = (float *) psAlloc((bounds->n) * sizeof(float)); + tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32)); for (i=0;in;i++) { @@ -1977,5 +1977,5 @@ { PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3); - char* strType; + psS8* strType; switch (x->type.type) { @@ -2074,10 +2074,10 @@ domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); - psScalar *tmpScalar = psScalarAlloc((double) + psScalar *tmpScalar = psScalarAlloc((psF64) interpolate1DF32(domain32->data.F32, range32->data.F32, domain32->n, order, - (float) x->data.F64), PS_TYPE_F64); + (psF32) x->data.F64), PS_TYPE_F64); psFree(range32); psFree(domain32); @@ -2089,5 +2089,5 @@ } else { - char* strType; + psS8* strType; PS_TYPE_NAME(strType,x->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, @@ -2115,6 +2115,6 @@ the spline fit functions require F32 and F64. *****************************************************************************/ -float psSpline1DEval( - float x, +psF32 psSpline1DEval( + psF32 x, const psSpline1D *spline ) @@ -2172,10 +2172,10 @@ for (i=0;in;i++) { tmpVector->data.F32[i] = psSpline1DEval( - (float) x->data.F64[i], + (psF32) x->data.F64[i], spline ); } } else { - char* strType; + psS8* strType; PS_TYPE_NAME(strType,x->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, Index: /trunk/psLib/src/math/psSpline.h =================================================================== --- /trunk/psLib/src/math/psSpline.h (revision 2787) +++ /trunk/psLib/src/math/psSpline.h (revision 2788) @@ -12,6 +12,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-12-02 21:12:51 $ +* @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -39,10 +39,10 @@ * \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 stddev, ///< Standard deviation for the Gaussian + * @return psF32 value on the gaussian curve given the input parameters + */ +psF32 psGaussian( + psF32 x, ///< Value at which to evaluate + psF32 mean, ///< Mean for the Gaussian + psF32 stddev, ///< Standard deviation for the Gaussian psBool normal ///< Indicates whether result should be normalized ); @@ -55,6 +55,6 @@ */ psVector* p_psGaussianDev( - float mean, ///< The mean of the Gaussian - float sigma, ///< The sigma of the Gaussian + psF32 mean, ///< The mean of the Gaussian + psF32 sigma, ///< The sigma of the Gaussian psS32 Npts ///< The size of the vector ); @@ -70,7 +70,7 @@ psPolynomialType type; ///< Polynomial type psS32 n; ///< Number of terms - float *coeff; ///< Coefficients - float *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + psF32 *coeff; ///< Coefficients + psF32 *coeffErr; ///< Error in coefficients + psU8 *mask; ///< Coefficient mask } psPolynomial1D; @@ -82,7 +82,7 @@ psS32 nX; ///< Number of terms in x psS32 nY; ///< Number of terms in y - float **coeff; ///< Coefficients - float **coeffErr; ///< Error in coefficients - char **mask; ///< Coefficients mask + psF32 **coeff; ///< Coefficients + psF32 **coeffErr; ///< Error in coefficients + psU8 **mask; ///< Coefficients mask } psPolynomial2D; @@ -95,7 +95,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - float ***coeff; ///< Coefficients - float ***coeffErr; ///< Error in coefficients - char ***mask; ///< Coefficients mask + psF32 ***coeff; ///< Coefficients + psF32 ***coeffErr; ///< Error in coefficients + psU8 ***mask; ///< Coefficients mask } psPolynomial3D; @@ -109,7 +109,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - float ****coeff; ///< Coefficients - float ****coeffErr; ///< Error in coefficients - char ****mask; ///< Coefficients mask + psF32 ****coeff; ///< Coefficients + psF32 ****coeffErr; ///< Error in coefficients + psU8 ****mask; ///< Coefficients mask } psPolynomial4D; @@ -160,42 +160,42 @@ /** Evaluates a 1-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial1DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial1DEval( const psPolynomial1D* myPoly, ///< Coefficients for the polynomial - float x ///< location at which to evaluate + psF32 x ///< location at which to evaluate ); /** Evaluates a 2-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial2DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial2DEval( const psPolynomial2D* myPoly, ///< Coefficients for the polynomial - float x, ///< x location at which to evaluate - float y ///< y location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y ///< y location at which to evaluate ); /** Evaluates a 3-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial3DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial3DEval( const psPolynomial3D* myPoly, ///< Coefficients for the polynomial - float x, ///< x location at which to evaluate - float y, ///< y location at which to evaluate - float z ///< z location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y, ///< y location at which to evaluate + psF32 z ///< z location at which to evaluate ); /** Evaluates a 4-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -float psPolynomial4DEval( + * @return psF32 result of polynomial at given location + */ +psF32 psPolynomial4DEval( const psPolynomial4D* myPoly, ///< Coefficients for the polynomial - float w, ///< w location at which to evaluate - float x, ///< x location at which to evaluate - float y, ///< y location at which to evaluate - float z ///< z location at which to evaluate + psF32 w, ///< w location at which to evaluate + psF32 x, ///< x location at which to evaluate + psF32 y, ///< y location at which to evaluate + psF32 z ///< z location at which to evaluate ); @@ -235,7 +235,7 @@ psPolynomialType type; ///< Polynomial type psS32 n; ///< Number of terms - double *coeff; ///< Coefficients - double *coeffErr; ///< Error in coefficients - char *mask; ///< Coefficient mask + psF64 *coeff; ///< Coefficients + psF64 *coeffErr; ///< Error in coefficients + psU8 *mask; ///< Coefficient mask } psDPolynomial1D; @@ -247,7 +247,7 @@ psS32 nX; ///< Number of terms in x psS32 nY; ///< Number of terms in y - double **coeff; ///< Coefficients - double **coeffErr; ///< Error in coefficients - char **mask; ///< Coefficients mask + psF64 **coeff; ///< Coefficients + psF64 **coeffErr; ///< Error in coefficients + psU8 **mask; ///< Coefficients mask } psDPolynomial2D; @@ -260,7 +260,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - double ***coeff; ///< Coefficients - double ***coeffErr; ///< Error in coefficients - char ***mask; ///< Coefficient mask + psF64 ***coeff; ///< Coefficients + psF64 ***coeffErr; ///< Error in coefficients + psU8 ***mask; ///< Coefficient mask } psDPolynomial3D; @@ -274,7 +274,7 @@ psS32 nY; ///< Number of terms in y psS32 nZ; ///< Number of terms in z - double ****coeff; ///< Coefficients - double ****coeffErr; ///< Error in coefficients - char ****mask; ///< Coefficients mask + psF64 ****coeff; ///< Coefficients + psF64 ****coeffErr; ///< Error in coefficients + psU8 ****mask; ///< Coefficients mask } psDPolynomial4D; @@ -324,42 +324,42 @@ /** Evaluates a double-precision 1-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -double psDPolynomial1DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial1DEval( const psDPolynomial1D* myPoly, ///< Coefficients for the polynomial - double x ///< Value at which to evaluate + psF64 x ///< Value at which to evaluate ); /** Evaluates a double-precision 2-D polynomial at specific coordinates. * - * @return float result of polynomial at given location - */ -double psDPolynomial2DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial2DEval( const psDPolynomial2D* myPoly, ///< Coefficients for the polynomial - double x, ///< Value x at which to evaluate - double y ///< Value y at which to evaluate + 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 float result of polynomial at given location - */ -double psDPolynomial3DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial3DEval( const psDPolynomial3D* myPoly, ///< Coefficients for the polynomial - double x, ///< Value x at which to evaluate - double y, ///< Value y at which to evaluate - double z ///< Value z at which to evaluate + 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 float result of polynomial at given location - */ -double psDPolynomial4DEval( + * @return psF32 result of polynomial at given location + */ +psF64 psDPolynomial4DEval( const psDPolynomial4D* myPoly, ///< Coefficients for the polynomial - double w, ///< Value w at which to evaluate - double x, ///< Value x at which to evaluate - double y, ///< Value y at which to evaluate - double z ///< Value z at which to evaluate + psF64 w, ///< Value w at which to evaluate + psF64 x, ///< Value x at which to evaluate + psF64 y, ///< Value y at which to evaluate + psF64 z ///< Value z at which to evaluate ); @@ -396,6 +396,6 @@ psS32 n; ///< The number of spline polynomials psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials - float *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. - float *domains; ///< The boundaries between each spline piece. Size is n+1. + psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. + psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. } psSpline1D; @@ -403,12 +403,12 @@ psSpline1D *psSpline1DAlloc(psS32 n, psS32 order, - float min, - float max); + psF32 min, + psF32 max); psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, psS32 order); -float psSpline1DEval( - float x, +psF32 psSpline1DEval( + psF32 x, const psSpline1D *spline ); @@ -427,8 +427,8 @@ psScalar *x); -float p_psNRSpline1DEval(psSpline1D *spline, +psF32 p_psNRSpline1DEval(psSpline1D *spline, const psVector* restrict x, const psVector* restrict y, - float X); + psF32 X); /* \} */// End of MathGroup Functions Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 2787) +++ /trunk/psLib/src/math/psStats.c (revision 2788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.107 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-22 00:54:28 $ + * @version $Revision: 1.108 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -70,5 +70,5 @@ /*****************************************************************************/ -psBool p_psGetStatValue(const psStats* stats, double *value) +psBool p_psGetStatValue(const psStats* stats, psF64 *value) { @@ -145,13 +145,13 @@ this routine sets stats->sampleMean to NAN. *****************************************************************************/ -int p_psVectorSampleMean(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorSampleMean(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float mean = 0.0; // The mean + psF32 mean = 0.0; // The mean psS32 count = 0; // # of points in this mean @@ -171,5 +171,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -185,5 +185,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -199,5 +199,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -207,5 +207,5 @@ mean += myVector->data.F32[i]; } - mean /= (float)myVector->n; + mean /= (psF32)myVector->n; } } @@ -283,11 +283,11 @@ this routine sets stats->max to NAN. *****************************************************************************/ -int p_psVectorMax(const psVector* restrict myVector, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorMax(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float max = -PS_MAX_F32; // The calculated maximum + psF32 max = -PS_MAX_F32; // The calculated maximum psS32 empty = true; // Does this vector have valid elements? @@ -348,11 +348,11 @@ this routine sets stats->min to NAN. *****************************************************************************/ -int p_psVectorMin(const psVector* restrict myVector, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorMin(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float min = PS_MAX_F32; // The calculated maximum + psF32 min = PS_MAX_F32; // The calculated maximum psS32 empty = true; // Does this vector have valid elements? @@ -612,5 +612,5 @@ *****************************************************************************/ psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, - float sigma) + psF32 sigma) { PS_PTR_CHECK_NULL(robustHistogram, NULL); @@ -619,6 +619,6 @@ psS32 i = 0; // Loop index variable psS32 j = 0; // Loop index variable - float iMid; - float jMid; + psF32 iMid; + psF32 jMid; psS32 numBins = robustHistogram->nums->n; psS32 numBounds = robustHistogram->bounds->n; @@ -626,6 +626,6 @@ psS32 jMin = 0; psS32 jMax = 0; - float firstBound = robustHistogram->bounds->data.F32[0]; - float lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + psF32 firstBound = robustHistogram->bounds->data.F32[0]; + psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; psScalar x; @@ -795,9 +795,9 @@ psS32 i = 0; // Loop index variable psS32 countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable + psF32 countFloat = 0.0; // # of data points being used + psF32 mean = 0.0; // The mean + psF32 diff = 0.0; // Used in calculating stdev + psF32 sumSquares = 0.0; // temporary variable + psF32 sumDiffs = 0.0; // temporary variable // This procedure requires the mean. If it has not been already @@ -865,5 +865,5 @@ psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", countInt); } else { - countFloat = (float)countInt; + countFloat = (psF32)countInt; stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); } @@ -889,9 +889,9 @@ psS32 i = 0; // Loop index variable psS32 countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable + psF32 countFloat = 0.0; // # of data points being used + psF32 mean = 0.0; // The mean + psF32 diff = 0.0; // Used in calculating stdev + psF32 sumSquares = 0.0; // temporary variable + psF32 sumDiffs = 0.0; // temporary variable // psF32 sum1; // psF32 sum2; @@ -982,5 +982,5 @@ stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor)); } else { - countFloat = (float)countInt; + countFloat = (psF32)countInt; stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); @@ -1003,16 +1003,16 @@ -2: warning *****************************************************************************/ -int p_psVectorClippedStats(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorClippedStats(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable psS32 j = 0; // Loop index variable - float clippedMean = 0.0; // self-explanatory - float clippedStdev = 0.0; // self-explanatory - float oldStanMean = 0.0; // Temporary variable - float oldStanStdev = 0.0; // Temporary variable + psF32 clippedMean = 0.0; // self-explanatory + psF32 clippedStdev = 0.0; // self-explanatory + psF32 oldStanMean = 0.0; // Temporary variable + psF32 oldStanStdev = 0.0; // Temporary variable psVector* tmpMask = NULL; // Temporary vector @@ -1240,8 +1240,8 @@ XXX: Create a 2nd-order polynomial version and solve for X analytically. *****************************************************************************/ -float p_ps1DPolyMedian(psPolynomial1D* myPoly, - float rangeLow, - float rangeHigh, - float getThisValue) +psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly, + psF32 rangeLow, + psF32 rangeHigh, + psF32 getThisValue) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1250,9 +1250,9 @@ // falls within the range of y-values of the polynomial "myPoly" in the // specified x-range (rangeLow:rangeHigh). - float fLo = psPolynomial1DEval( + psF32 fLo = psPolynomial1DEval( myPoly, rangeLow ); - float fHi = psPolynomial1DEval( + psF32 fHi = psPolynomial1DEval( myPoly, rangeHigh @@ -1266,7 +1266,7 @@ psS32 numIterations = 0; - float midpoint = 0.0; - float oldMidpoint = 1.0; - float f = 0.0; + psF32 midpoint = 0.0; + psF32 oldMidpoint = 1.0; + psF32 f = 0.0; while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) { @@ -1309,8 +1309,8 @@ tests to ensure that binNum is within acceptable ranges for both vectors. *****************************************************************************/ -float fitQuadraticSearchForYThenReturnX(psVector *xVec, +psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, psVector *yVec, psS32 binNum, - float yVal) + psF32 yVal) { PS_VECTOR_CHECK_NULL(xVec, NAN); @@ -1331,11 +1331,11 @@ psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - float tmpFloat; + psF32 tmpFloat; if ((binNum > 0) && (binNum < (yVec->n - 2))) { // The general case. We have all three points. - x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); - x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); - x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); + x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); y->data.F64[0] = yVec->data.F32[binNum - 1]; y->data.F64[1] = yVec->data.F32[binNum]; @@ -1346,5 +1346,5 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, - (double)yVal,y->data.F64[2],y->data.F64[0]); + (psF64)yVal,y->data.F64[2],y->data.F64[0]); } yErr->data.F64[0] = 1.0; @@ -1425,13 +1425,13 @@ XXX: Check for errors in psLib routines that we call. *****************************************************************************/ -int p_psVectorRobustStats(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorRobustStats(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psHistogram* robustHistogram = NULL; psVector* robustHistogramVector = NULL; - float binSize = 0.0; // Size of the histogram bins + psF32 binSize = 0.0; // Size of the histogram bins psS32 LQBinNum = -1; // Bin num for lower quartile psS32 UQBinNum = -1; // Bin num for upper quartile @@ -1439,17 +1439,17 @@ psS32 i = 0; // Loop index variable psS32 modeBinNum = 0; - float modeBinCount = 0.0; - float dL = 0.0; + psF32 modeBinCount = 0.0; + psF32 dL = 0.0; psS32 numBins = 0; - float myMean = 0.0; - float myStdev = 0.0; - float countFloat = 0.0; - float diff = 0.0; - float sumSquares = 0.0; - float sumDiffs = 0.0; + psF32 myMean = 0.0; + psF32 myStdev = 0.0; + psF32 countFloat = 0.0; + psF32 diff = 0.0; + psF32 sumSquares = 0.0; + psF32 sumDiffs = 0.0; psVector* cumulativeRobustSums = NULL; - float sumRobust = 0.0; - float sumN50 = 0.0; - float sumNfit = 0.0; + psF32 sumRobust = 0.0; + psF32 sumN50 = 0.0; + psF32 sumNfit = 0.0; psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; @@ -1460,5 +1460,5 @@ // that by 10.0; //XXX: add errors - int rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); + psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); if (rc != 0) { psError(PS_ERR_UNEXPECTED_NULL, @@ -1614,5 +1614,5 @@ for (i=0;in;i++) { myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); - ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i; + ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i; y->data.F32[i] = robustHistogramVector->data.F32[i]; } @@ -1770,5 +1770,5 @@ The histogram structure *****************************************************************************/ -psHistogram* psHistogramAlloc(float lower, float upper, psS32 n) +psHistogram* psHistogramAlloc(psF32 lower, psF32 upper, psS32 n) { PS_INT_CHECK_POSITIVE(n, NULL); @@ -1777,5 +1777,5 @@ psS32 i = 0; // Loop index variable psHistogram* newHist = NULL; // The new histogram structure - float binSize = 0.0; // The histogram bin size + psF32 binSize = 0.0; // The histogram bin size // Allocate memory for the new histogram structure. If there are N @@ -1787,10 +1787,10 @@ // Calculate the bounds for each bin. - binSize = (upper - lower) / (float)n; + binSize = (upper - lower) / (psF32)n; // XXX: Is the following necessary? It prevents the max data point // from being in a non-existant bin. binSize += FLT_EPSILON; for (i = 0; i < n + 1; i++) { - newHist->bounds->data.F32[i] = lower + (binSize * (float)i); + newHist->bounds->data.F32[i] = lower + (binSize * (psF32)i); } @@ -1857,4 +1857,17 @@ } +/***************************************************************************** +UpdateHistogramBins(binNum, out, data, error): This routine is to be used when +updating the histogram in the presence of errors in the input data. We treat +the data point as a boxcar PDF and update a range of points surrounding the +histogram bin which contains the point. The width of that boxcar is defined +as 2.35 * error. Inputs: + binNum: the bin number of the data point in the histogram + out: the histogram structure + data: the data point value + error: the error in that data point + +XXX: Must test this. + *****************************************************************************/ psS32 UpdateHistogramBins(psS32 binNum, psHistogram* out, @@ -1863,20 +1876,68 @@ { PS_PTR_CHECK_NULL(out, -1); - PS_INT_CHECK_RANGE(binNum, 0, out->nums->n-1, -2); - /* - psF32 width = 2.35 * error; - psF32 centerBinWidth = out->bounds->data.F32[binNum+1] - out->bounds->data.F32[binNum] - psF32 boxcarCenter = (out->bounds->data.F32[binNum] + out->bounds->data.F32[binNum+1]) / 2.0; - - if (width <= centerBinWidth) { - out->nums->data.F32[binNum]+= 1.0; - } else { - out->nums->data.F32[binNum]+= centerBinWidth / width; - // XXX: walk to the left, adding fractional values. - // XXX: walk to the right, adding fractional values. - - - } - */ + PS_PTR_CHECK_NULL(out->bounds, -1); + PS_PTR_CHECK_NULL(out->nums, -1); + PS_INT_CHECK_RANGE(binNum, 0, ((out->nums->n)-1), -2); + PS_FLOAT_COMPARE(0.0, error, -3); + PS_FLOAT_CHECK_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4); + + psF32 boxcarWidth = 2.35 * error; + psF32 boxcarCenter = (out->bounds->data.F32[binNum] + + out->bounds->data.F32[binNum+1]) / 2.0; + psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); + psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); + psS32 bin; + psS32 boxcarLeftBinNum; + psS32 boxcarRightBinNum; + + // Determine the left endpoint of the boxcar for the PDF. + for (bin=binNum ; bin >= 0 ; bin--) { + if (out->nums->data.F32[bin] <= boxcarLeft) { + boxcarLeftBinNum = bin; + break; + } + } + + // Determine the right endpoint of the boxcar for the PDF. + for (bin=binNum ; bin < out->nums->n ; bin++) { + if (out->nums->data.F32[bin] >= boxcarRight) { + boxcarRightBinNum = bin; + break; + } + } + + // + // If the boxcar fits entirely inside this bin, then simply add 1.0 to the + // bin and return. + // + if (boxcarLeftBinNum == boxcarRightBinNum) { + out->nums->data.F32[binNum]+= 1.0; + return(0); + } + + // + // If we get here, multiple bins must be updated. We handle the left + // endpoint, and right endpoint differently. + // + out->nums->data.F32[boxcarLeftBinNum]+= + (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth; + + // + // Loop through the center bins, if any. + // + for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) { + out->nums->data.F32[bin]+= + (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth; + } + + // + // Handle the right endpoint differently. + // + out->nums->data.F32[boxcarRightBinNum]+= + (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth; + + // + // Return 0 on success. + // return(0); } @@ -1906,5 +1967,8 @@ { PS_PTR_CHECK_NULL(out, NULL); + PS_VECTOR_CHECK_NULL(out->bounds, NULL); PS_VECTOR_CHECK_TYPE(out->bounds, PS_TYPE_F32, NULL); + PS_INT_CHECK_NON_NEGATIVE(out->bounds->n, NULL); + PS_VECTOR_CHECK_NULL(out->nums, NULL); PS_VECTOR_CHECK_TYPE(out->nums, PS_TYPE_F32, NULL); PS_INT_CHECK_NON_NEGATIVE(out->nums->n, NULL); @@ -1920,16 +1984,24 @@ psS32 i = 0; // Loop index variable - float binSize = 0.0; // Histogram bin size + psF32 binSize = 0.0; // Histogram bin size psS32 binNum = 0; // A temporary bin number psS32 numBins = 0; // The total number of bins - psS32 tmp = 0; psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; - psVector* inF32; - psS32 mustFreeTmp = 1; + psVector* inF32 = NULL; + psVector* errorsF32 = NULL; + psS32 mustFreeVectorIn = 1; + psS32 mustFreeVectorErrors = 1; + + // Convert input and errors vectors to F32 if necessary. inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { inF32 = (psVector *) in; - mustFreeTmp = 0; + mustFreeVectorIn = 0; + } + errorsF32 = p_psConvertToF32((psVector *) errors); + if (errorsF32 == NULL) { + errorsF32 = (psVector *) errors; + mustFreeVectorErrors = 0; } @@ -1950,9 +2022,9 @@ binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); - if (errors != NULL) { + if (errorsF32 != NULL) { // XXX: Check return codes. UpdateHistogramBins(binNum, out, inF32->data.F32[i], - errors->data.F32[i]); + errorsF32->data.F32[i]); } else { // XXX: This if-statement really shouldn't be necessary. @@ -1969,16 +2041,16 @@ // correct bin number requires a bit more work. tmpScalar.data.F32 = inF32->data.F32[i]; - tmp = p_psVectorBinDisect(out->bounds, &tmpScalar); - if (tmp < 0) { + binNum = p_psVectorBinDisect(out->bounds, &tmpScalar); + if (binNum < 0) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: psVectorHistogram(): element outside histogram bounds.\n"); } else { - if (errors != NULL) { + if (errorsF32 != NULL) { // XXX: Check return codes. - UpdateHistogramBins(tmp, out, + UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i]); } else { - (out->nums->data.F32[tmp])+= 1.0; + (out->nums->data.F32[binNum])+= 1.0; } } @@ -1988,6 +2060,9 @@ } - if (mustFreeTmp == 1) { + if (mustFreeVectorIn == 1) { psFree(inF32); + } + if (mustFreeVectorErrors == 1) { + psFree(errorsF32); } return (out); @@ -2015,26 +2090,51 @@ tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.S8[i]; + tmp->data.F32[i] = (psF32)in->data.S8[i]; + } + } else if (in->type.type == PS_TYPE_S16) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32) in->data.S16[i]; + } + } else if (in->type.type == PS_TYPE_S32) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.S32[i]; + } + } else if (in->type.type == PS_TYPE_S64) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.S64[i]; + } + } else if (in->type.type == PS_TYPE_U8) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.U8[i]; } } else if (in->type.type == PS_TYPE_U16) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.U16[i]; - } - } else if (in->type.type == PS_TYPE_U8) { + tmp->data.F32[i] = (psF32)in->data.U16[i]; + } + } else if (in->type.type == PS_TYPE_U32) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.U8[i]; + tmp->data.F32[i] = (psF32)in->data.U32[i]; + } + } else if (in->type.type == PS_TYPE_U64) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.U64[i]; } } else if (in->type.type == PS_TYPE_F64) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.F64[i]; + tmp->data.F32[i] = (psF32)in->data.F64[i]; } } else if (in->type.type == PS_TYPE_F32) { // do nothing } else { - char* strType; - PS_TYPE_NAME(strType,in->type.type); + psS8* strType; + PS_TYPE_NAME(strType, in->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED, @@ -2074,6 +2174,6 @@ } - psVector* inF32; - psVector* errorsF32; + psVector* inF32 = NULL; + psVector* errorsF32 = NULL; psS32 mustFreeVectorIn = 1; psS32 mustFreeVectorErrors = 1; @@ -2178,2 +2278,3 @@ return (stats); } + Index: /trunk/psLib/src/math/psStats.h =================================================================== --- /trunk/psLib/src/math/psStats.h (revision 2787) +++ /trunk/psLib/src/math/psStats.h (revision 2788) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-21 20:42:07 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -56,25 +56,25 @@ typedef struct { - 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 + 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; ///< psS32 robustNfit; ///< - double clippedMean; ///< Nsigma clipped mean - double clippedStdev; ///< standard deviation after clipping + psF64 clippedMean; ///< Nsigma clipped mean + psF64 clippedStdev; ///< standard deviation after clipping psS32 clippedNvalues; ///< ??? - double clipSigma; ///< Nsigma used for clipping; user input + psF64 clipSigma; ///< Nsigma used for clipping; user input psS32 clipIter; ///< Number of clipping iterations; user input - double min; ///< minimum data value in array - double max; ///< maximum data value in array - double binsize; ///< + psF64 min; ///< minimum data value in array + psF64 max; ///< maximum data value in array + psF64 binsize; ///< psStatsOptions options; ///< bitmask of calculated values } @@ -132,6 +132,6 @@ */ psHistogram* psHistogramAlloc( - float lower, ///< Lower limit for the bins - float upper, ///< Upper limit for the bins + psF32 lower, ///< Lower limit for the bins + psF32 upper, ///< Upper limit for the bins psS32 n ///< Number of bins ); @@ -172,5 +172,5 @@ ///< the statistic struct to operate on - double *value + psF64 *value ///< if return is true, this is set to the specified statistic value by stats->options );