Changeset 2788
- Timestamp:
- Dec 21, 2004, 7:09:32 PM (22 years ago)
- Location:
- trunk/psLib/src
- Files:
-
- 16 edited
-
dataManip/psConstants.h (modified) (1 diff)
-
dataManip/psFunctions.c (modified) (58 diffs)
-
dataManip/psFunctions.h (modified) (16 diffs)
-
dataManip/psMinimize.c (modified) (36 diffs)
-
dataManip/psMinimize.h (modified) (3 diffs)
-
dataManip/psStats.c (modified) (40 diffs)
-
dataManip/psStats.h (modified) (4 diffs)
-
math/psConstants.h (modified) (1 diff)
-
math/psMinimize.c (modified) (36 diffs)
-
math/psMinimize.h (modified) (3 diffs)
-
math/psPolynomial.c (modified) (58 diffs)
-
math/psPolynomial.h (modified) (16 diffs)
-
math/psSpline.c (modified) (58 diffs)
-
math/psSpline.h (modified) (16 diffs)
-
math/psStats.c (modified) (40 diffs)
-
math/psStats.h (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/dataManip/psConstants.h
r2778 r2788 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1.4 8$ $Name: not supported by cvs2svn $9 * @date $Date: 2004-12-2 1 20:42:07$8 * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2004-12-22 05:09:32 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 12 * 13 * XXX: Add parenthesis around all arguments so that these macros can be 14 * called with complex expressions. 15 * 16 * 17 * 12 18 */ 13 19 -
trunk/psLib/src/dataManip/psFunctions.c
r2778 r2788 7 7 * polynomials. It also contains a Gaussian functions. 8 8 * 9 * @version $Revision: 1.7 4$ $Name: not supported by cvs2svn $10 * @date $Date: 2004-12-2 1 20:42:07$9 * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-12-22 05:09:32 $ 11 11 * 12 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 55 55 static void dPolynomial4DFree(psDPolynomial4D* myPoly); 56 56 static void spline1DFree(psSpline1D *tmpSpline); 57 static psS32 vectorBinDisectF32( float *bins,psS32 numBins,floatx);57 static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); 58 58 static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); 59 59 … … 277 277 278 278 *****************************************************************************/ 279 static float ordPolynomial1DEval(floatx, const psPolynomial1D* myPoly)279 static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) 280 280 { 281 281 psS32 loop_x = 0; 282 floatpolySum = 0.0;283 floatxSum = 1.0;282 psF32 polySum = 0.0; 283 psF32 xSum = 1.0; 284 284 285 285 psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, … … 307 307 // XXX: How does the mask vector effect Crenshaw's formula? 308 308 // XXX: We assume that x is scaled between -1.0 and 1.0; 309 static float chebPolynomial1DEval(floatx, const psPolynomial1D* myPoly)309 static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) 310 310 { 311 311 // PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0); … … 314 314 psS32 n; 315 315 psS32 i; 316 floattmp;316 psF32 tmp; 317 317 318 318 n = myPoly->n; … … 337 337 psS32 n; 338 338 psS32 i; 339 floattmp;339 psF32 tmp; 340 340 psPolynomial1D **chebPolys = NULL; 341 341 … … 355 355 } 356 356 357 static float ordPolynomial2DEval(float x, floaty, const psPolynomial2D* myPoly)357 static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) 358 358 { 359 359 PS_POLY_CHECK_NULL(myPoly, NAN); … … 361 361 psS32 loop_x = 0; 362 362 psS32 loop_y = 0; 363 floatpolySum = 0.0;364 floatxSum = 1.0;365 floatySum = 1.0;363 psF32 polySum = 0.0; 364 psF32 xSum = 1.0; 365 psF32 ySum = 1.0; 366 366 367 367 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 379 379 } 380 380 381 static float chebPolynomial2DEval(float x, floaty, const psPolynomial2D* myPoly)381 static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) 382 382 { 383 383 PS_POLY_CHECK_NULL(myPoly, NAN); … … 386 386 psS32 loop_y = 0; 387 387 psS32 i = 0; 388 floatpolySum = 0.0;388 psF32 polySum = 0.0; 389 389 psPolynomial1D* *chebPolys = NULL; 390 390 psS32 maxChebyPoly = 0; … … 414 414 } 415 415 416 static float ordPolynomial3DEval(float x, float y, floatz, const psPolynomial3D* myPoly)416 static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) 417 417 { 418 418 psS32 loop_x = 0; 419 419 psS32 loop_y = 0; 420 420 psS32 loop_z = 0; 421 floatpolySum = 0.0;422 floatxSum = 1.0;423 floatySum = 1.0;424 floatzSum = 1.0;421 psF32 polySum = 0.0; 422 psF32 xSum = 1.0; 423 psF32 ySum = 1.0; 424 psF32 zSum = 1.0; 425 425 426 426 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 442 442 } 443 443 444 static float chebPolynomial3DEval(float x, float y, floatz, const psPolynomial3D* myPoly)444 static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) 445 445 { 446 446 psS32 loop_x = 0; … … 448 448 psS32 loop_z = 0; 449 449 psS32 i = 0; 450 floatpolySum = 0.0;450 psF32 polySum = 0.0; 451 451 psPolynomial1D* *chebPolys = NULL; 452 452 psS32 maxChebyPoly = 0; … … 483 483 } 484 484 485 static float ordPolynomial4DEval(float w, float x, float y, floatz, const psPolynomial4D* myPoly)485 static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) 486 486 { 487 487 psS32 loop_w = 0; … … 489 489 psS32 loop_y = 0; 490 490 psS32 loop_z = 0; 491 floatpolySum = 0.0;492 floatwSum = 1.0;493 floatxSum = 1.0;494 floatySum = 1.0;495 floatzSum = 1.0;491 psF32 polySum = 0.0; 492 psF32 wSum = 1.0; 493 psF32 xSum = 1.0; 494 psF32 ySum = 1.0; 495 psF32 zSum = 1.0; 496 496 497 497 for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { … … 517 517 } 518 518 519 static float chebPolynomial4DEval(float w, float x, float y, floatz, const psPolynomial4D* myPoly)519 static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) 520 520 { 521 521 psS32 loop_w = 0; … … 524 524 psS32 loop_z = 0; 525 525 psS32 i = 0; 526 floatpolySum = 0.0;526 psF32 polySum = 0.0; 527 527 psPolynomial1D* *chebPolys = NULL; 528 528 psS32 maxChebyPoly = 0; … … 568 568 Polynomial coefficients will be accessed in [w][x][y][z] fashion. 569 569 *****************************************************************************/ 570 static double dOrdPolynomial1DEval(doublex, const psDPolynomial1D* myPoly)570 static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) 571 571 { 572 572 psS32 loop_x = 0; 573 doublepolySum = 0.0;574 doublexSum = 1.0;573 psF64 polySum = 0.0; 574 psF64 xSum = 1.0; 575 575 576 576 for (loop_x = 0; loop_x < myPoly->n; loop_x++) { … … 586 586 // XXX: You can do this without having to psAlloc() vector d. 587 587 // XXX: How does the mask vector effect Crenshaw's formula? 588 static double dChebPolynomial1DEval(doublex, const psDPolynomial1D* myPoly)588 static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) 589 589 { 590 590 psVector *d; 591 591 psS32 n; 592 592 psS32 i; 593 doubletmp;593 psF64 tmp; 594 594 595 595 n = myPoly->n; … … 611 611 } 612 612 613 static double dOrdPolynomial2DEval(double x, doubley, const psDPolynomial2D* myPoly)613 static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) 614 614 { 615 615 psS32 loop_x = 0; 616 616 psS32 loop_y = 0; 617 doublepolySum = 0.0;618 doublexSum = 1.0;619 doubleySum = 1.0;617 psF64 polySum = 0.0; 618 psF64 xSum = 1.0; 619 psF64 ySum = 1.0; 620 620 621 621 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 633 633 } 634 634 635 static double dChebPolynomial2DEval(double x, doubley, const psDPolynomial2D* myPoly)635 static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) 636 636 { 637 637 psS32 loop_x = 0; 638 638 psS32 loop_y = 0; 639 639 psS32 i = 0; 640 doublepolySum = 0.0;640 psF64 polySum = 0.0; 641 641 psPolynomial1D* *chebPolys = NULL; 642 642 psS32 maxChebyPoly = 0; … … 667 667 } 668 668 669 static double dOrdPolynomial3DEval(double x, double y, doublez, const psDPolynomial3D* myPoly)669 static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) 670 670 { 671 671 psS32 loop_x = 0; 672 672 psS32 loop_y = 0; 673 673 psS32 loop_z = 0; 674 doublepolySum = 0.0;675 doublexSum = 1.0;676 doubleySum = 1.0;677 doublezSum = 1.0;674 psF64 polySum = 0.0; 675 psF64 xSum = 1.0; 676 psF64 ySum = 1.0; 677 psF64 zSum = 1.0; 678 678 679 679 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 695 695 } 696 696 697 static double dChebPolynomial3DEval(double x, double y, doublez, const psDPolynomial3D* myPoly)697 static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) 698 698 { 699 699 psS32 loop_x = 0; … … 701 701 psS32 loop_z = 0; 702 702 psS32 i = 0; 703 doublepolySum = 0.0;703 psF64 polySum = 0.0; 704 704 psPolynomial1D* *chebPolys = NULL; 705 705 psS32 maxChebyPoly = 0; … … 736 736 } 737 737 738 static double dOrdPolynomial4DEval(double w, double x, double y, doublez, const psDPolynomial4D* myPoly)738 static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) 739 739 { 740 740 psS32 loop_w = 0; … … 742 742 psS32 loop_y = 0; 743 743 psS32 loop_z = 0; 744 doublepolySum = 0.0;745 doublewSum = 1.0;746 doublexSum = 1.0;747 doubleySum = 1.0;748 doublezSum = 1.0;744 psF64 polySum = 0.0; 745 psF64 wSum = 1.0; 746 psF64 xSum = 1.0; 747 psF64 ySum = 1.0; 748 psF64 zSum = 1.0; 749 749 750 750 for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { … … 770 770 } 771 771 772 static double dChebPolynomial4DEval(double w, double x, double y, doublez, const psDPolynomial4D* myPoly)772 static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) 773 773 { 774 774 psS32 loop_w = 0; … … 777 777 psS32 loop_z = 0; 778 778 psS32 i = 0; 779 doublepolySum = 0.0;779 psF64 polySum = 0.0; 780 780 psPolynomial1D* *chebPolys = NULL; 781 781 psS32 maxChebyPoly = 0; … … 828 828 *****************************************************************************/ 829 829 #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ 830 static floatfullInterpolate1D##TYPE(ps##TYPE *domain, \830 static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ 831 831 ps##TYPE *range, \ 832 832 psS32 n, \ … … 896 896 LaGrange interpolation. 897 897 *****************************************************************************/ 898 static float interpolate1DF32(float*domain,899 float*range,898 static psF32 interpolate1DF32(psF32 *domain, 899 psF32 *range, 900 900 psS32 n, 901 901 psS32 order, 902 floatx)902 psF32 x) 903 903 { 904 904 psS32 binNum; … … 941 941 evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] 942 942 *****************************************************************************/ 943 float psGaussian(float x, float mean, floatsigma, psBool normal)944 { 945 floattmp = 1.0;943 psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal) 944 { 945 psF32 tmp = 1.0; 946 946 947 947 psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, … … 967 967 XXX: There is no way to seed the random generator. 968 968 *****************************************************************************/ 969 psVector* p_psGaussianDev( float mean, floatsigma, psS32 Npts)969 psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) 970 970 { 971 971 PS_INT_CHECK_NON_NEGATIVE(Npts, NULL); … … 1007 1007 newPoly->type = type; 1008 1008 newPoly->n = n; 1009 newPoly->coeff = ( float *)psAlloc(n * sizeof(float));1010 newPoly->coeffErr = ( float *)psAlloc(n * sizeof(float));1011 newPoly->mask = ( char *)psAlloc(n * sizeof(char));1009 newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32)); 1010 newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32)); 1011 newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); 1012 1012 for (i = 0; i < n; i++) { 1013 1013 newPoly->coeff[i] = 0.0; … … 1036 1036 newPoly->nY = nY; 1037 1037 1038 newPoly->coeff = ( float **)psAlloc(nX * sizeof(float*));1039 newPoly->coeffErr = ( float **)psAlloc(nX * sizeof(float*));1040 newPoly->mask = ( char **)psAlloc(nX * sizeof(char*));1038 newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *)); 1039 newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *)); 1040 newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); 1041 1041 for (x = 0; x < nX; x++) { 1042 newPoly->coeff[x] = ( float *)psAlloc(nY * sizeof(float));1043 newPoly->coeffErr[x] = ( float *)psAlloc(nY * sizeof(float));1044 newPoly->mask[x] = ( char *)psAlloc(nY * sizeof(char));1042 newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); 1043 newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); 1044 newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); 1045 1045 } 1046 1046 for (x = 0; x < nX; x++) { … … 1075 1075 newPoly->nZ = nZ; 1076 1076 1077 newPoly->coeff = ( float ***)psAlloc(nX * sizeof(float**));1078 newPoly->coeffErr = ( float ***)psAlloc(nX * sizeof(float**));1079 newPoly->mask = ( char ***)psAlloc(nX * sizeof(char**));1077 newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1078 newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1079 newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1080 1080 for (x = 0; x < nX; x++) { 1081 newPoly->coeff[x] = ( float **)psAlloc(nY * sizeof(float*));1082 newPoly->coeffErr[x] = ( float **)psAlloc(nY * sizeof(float*));1083 newPoly->mask[x] = ( char **)psAlloc(nY * sizeof(char*));1081 newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1082 newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1083 newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1084 1084 for (y = 0; y < nY; y++) { 1085 newPoly->coeff[x][y] = ( float *)psAlloc(nZ * sizeof(float));1086 newPoly->coeffErr[x][y] = ( float *)psAlloc(nZ * sizeof(float));1087 newPoly->mask[x][y] = ( char *)psAlloc(nZ * sizeof(char));1085 newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1086 newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1087 newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1088 1088 } 1089 1089 } … … 1124 1124 newPoly->nZ = nZ; 1125 1125 1126 newPoly->coeff = ( float ****)psAlloc(nW * sizeof(float***));1127 newPoly->coeffErr = ( float ****)psAlloc(nW * sizeof(float***));1128 newPoly->mask = ( char ****)psAlloc(nW * sizeof(char***));1126 newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); 1127 newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); 1128 newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); 1129 1129 for (w = 0; w < nW; w++) { 1130 newPoly->coeff[w] = ( float ***)psAlloc(nX * sizeof(float**));1131 newPoly->coeffErr[w] = ( float ***)psAlloc(nX * sizeof(float**));1132 newPoly->mask[w] = ( char ***)psAlloc(nX * sizeof(char**));1130 newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1131 newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1132 newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1133 1133 for (x = 0; x < nX; x++) { 1134 newPoly->coeff[w][x] = ( float **)psAlloc(nY * sizeof(float*));1135 newPoly->coeffErr[w][x] = ( float **)psAlloc(nY * sizeof(float*));1136 newPoly->mask[w][x] = ( char **)psAlloc(nY * sizeof(char*));1134 newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1135 newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1136 newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1137 1137 for (y = 0; y < nY; y++) { 1138 newPoly->coeff[w][x][y] = ( float *)psAlloc(nZ * sizeof(float));1139 newPoly->coeffErr[w][x][y] = ( float *)psAlloc(nZ * sizeof(float));1140 newPoly->mask[w][x][y] = ( char *)psAlloc(nZ * sizeof(char));1138 newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1139 newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1140 newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1141 1141 } 1142 1142 } … … 1157 1157 } 1158 1158 1159 float psPolynomial1DEval(const psPolynomial1D* myPoly, floatx)1159 psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x) 1160 1160 { 1161 1161 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1195 1195 } 1196 1196 1197 float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, floaty)1197 psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y) 1198 1198 { 1199 1199 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1250 1250 } 1251 1251 1252 float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, floatz)1252 psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z) 1253 1253 { 1254 1254 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1316 1316 } 1317 1317 1318 float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, floatz)1318 psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z) 1319 1319 { 1320 1320 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1406 1406 newPoly->type = type; 1407 1407 newPoly->n = n; 1408 newPoly->coeff = ( double *)psAlloc(n * sizeof(double));1409 newPoly->coeffErr = ( double *)psAlloc(n * sizeof(double));1410 newPoly->mask = ( char *)psAlloc(n * sizeof(char));1408 newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64)); 1409 newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64)); 1410 newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); 1411 1411 for (i = 0; i < n; i++) { 1412 1412 newPoly->coeff[i] = 0.0; … … 1435 1435 newPoly->nY = nY; 1436 1436 1437 newPoly->coeff = ( double **)psAlloc(nX * sizeof(double*));1438 newPoly->coeffErr = ( double **)psAlloc(nX * sizeof(double*));1439 newPoly->mask = ( char **)psAlloc(nX * sizeof(char*));1437 newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *)); 1438 newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *)); 1439 newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); 1440 1440 for (x = 0; x < nX; x++) { 1441 newPoly->coeff[x] = ( double *)psAlloc(nY * sizeof(double));1442 newPoly->coeffErr[x] = ( double *)psAlloc(nY * sizeof(double));1443 newPoly->mask[x] = ( char *)psAlloc(nY * sizeof(char));1441 newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); 1442 newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); 1443 newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); 1444 1444 } 1445 1445 for (x = 0; x < nX; x++) { … … 1474 1474 newPoly->nZ = nZ; 1475 1475 1476 newPoly->coeff = ( double ***)psAlloc(nX * sizeof(double**));1477 newPoly->coeffErr = ( double ***)psAlloc(nX * sizeof(double**));1478 newPoly->mask = ( char ***)psAlloc(nX * sizeof(char**));1476 newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1477 newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1478 newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1479 1479 for (x = 0; x < nX; x++) { 1480 newPoly->coeff[x] = ( double **)psAlloc(nY * sizeof(double*));1481 newPoly->coeffErr[x] = ( double **)psAlloc(nY * sizeof(double*));1482 newPoly->mask[x] = ( char **)psAlloc(nY * sizeof(char*));1480 newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1481 newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1482 newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1483 1483 for (y = 0; y < nY; y++) { 1484 newPoly->coeff[x][y] = ( double *)psAlloc(nZ * sizeof(double));1485 newPoly->coeffErr[x][y] = ( double *)psAlloc(nZ * sizeof(double));1486 newPoly->mask[x][y] = ( char *)psAlloc(nZ * sizeof(char));1484 newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1485 newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1486 newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1487 1487 } 1488 1488 } … … 1523 1523 newPoly->nZ = nZ; 1524 1524 1525 newPoly->coeff = ( double ****)psAlloc(nW * sizeof(double***));1526 newPoly->coeffErr = ( double ****)psAlloc(nW * sizeof(double***));1527 newPoly->mask = ( char ****)psAlloc(nW * sizeof(char***));1525 newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); 1526 newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); 1527 newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); 1528 1528 for (w = 0; w < nW; w++) { 1529 newPoly->coeff[w] = ( double ***)psAlloc(nX * sizeof(double**));1530 newPoly->coeffErr[w] = ( double ***)psAlloc(nX * sizeof(double**));1531 newPoly->mask[w] = ( char ***)psAlloc(nX * sizeof(char**));1529 newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1530 newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1531 newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1532 1532 for (x = 0; x < nX; x++) { 1533 newPoly->coeff[w][x] = ( double **)psAlloc(nY * sizeof(double*));1534 newPoly->coeffErr[w][x] = ( double **)psAlloc(nY * sizeof(double*));1535 newPoly->mask[w][x] = ( char **)psAlloc(nY * sizeof(char*));1533 newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1534 newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1535 newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1536 1536 for (y = 0; y < nY; y++) { 1537 newPoly->coeff[w][x][y] = ( double *)psAlloc(nZ * sizeof(double));1538 newPoly->coeffErr[w][x][y] = ( double *)psAlloc(nZ * sizeof(double));1539 newPoly->mask[w][x][y] = ( char *)psAlloc(nZ * sizeof(char));1537 newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1538 newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1539 newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1540 1540 } 1541 1541 } … … 1557 1557 1558 1558 1559 double psDPolynomial1DEval(const psDPolynomial1D* myPoly, doublex)1559 psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x) 1560 1560 { 1561 1561 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1598 1598 1599 1599 1600 doublepsDPolynomial2DEval(const psDPolynomial2D* myPoly,1601 doublex,1602 doubley)1600 psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly, 1601 psF64 x, 1602 psF64 y) 1603 1603 { 1604 1604 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1655 1655 1656 1656 1657 doublepsDPolynomial3DEval(const psDPolynomial3D* myPoly,1658 doublex,1659 doubley,1660 doublez)1657 psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly, 1658 psF64 x, 1659 psF64 y, 1660 psF64 z) 1661 1661 { 1662 1662 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1724 1724 } 1725 1725 1726 doublepsDPolynomial4DEval(const psDPolynomial4D* myPoly,1727 doublew,1728 doublex,1729 doubley,1730 doublez)1726 psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly, 1727 psF64 w, 1728 psF64 x, 1729 psF64 y, 1730 psF64 z) 1731 1731 { 1732 1732 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1810 1810 // psS32 n; 1811 1811 // psPolynomial1D **spline; 1812 // float*p_psDeriv2;1813 // float*domains;1812 // psF32 *p_psDeriv2; 1813 // psF32 *domains; 1814 1814 //} psSpline1D; 1815 1815 … … 1822 1822 psSpline1D *psSpline1DAlloc(psS32 numSplines, 1823 1823 psS32 order, 1824 floatmin,1825 floatmax)1824 psF32 min, 1825 psF32 max) 1826 1826 { 1827 1827 PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL); … … 1831 1831 psSpline1D *tmp = NULL; 1832 1832 psS32 i; 1833 floattmpDomain;1834 floatwidth;1833 psF32 tmpDomain; 1834 psF32 width; 1835 1835 1836 1836 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); … … 1845 1845 tmp->p_psDeriv2 = NULL; 1846 1846 1847 tmp->domains = ( float *) psAlloc((numSplines+1) * sizeof(float));1848 width = (max - min) / (( float) numSplines);1847 tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32)); 1848 width = (max - min) / ((psF32) numSplines); 1849 1849 1850 1850 (tmp->domains)[0] = min; … … 1884 1884 } 1885 1885 1886 tmp->domains = ( float *) psAlloc((bounds->n) * sizeof(float));1886 tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32)); 1887 1887 1888 1888 for (i=0;i<bounds->n;i++) { … … 1977 1977 { 1978 1978 PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3); 1979 char* strType;1979 psS8* strType; 1980 1980 1981 1981 switch (x->type.type) { … … 2074 2074 domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); 2075 2075 2076 psScalar *tmpScalar = psScalarAlloc(( double)2076 psScalar *tmpScalar = psScalarAlloc((psF64) 2077 2077 interpolate1DF32(domain32->data.F32, 2078 2078 range32->data.F32, 2079 2079 domain32->n, 2080 2080 order, 2081 ( float) x->data.F64), PS_TYPE_F64);2081 (psF32) x->data.F64), PS_TYPE_F64); 2082 2082 psFree(range32); 2083 2083 psFree(domain32); … … 2089 2089 2090 2090 } else { 2091 char* strType;2091 psS8* strType; 2092 2092 PS_TYPE_NAME(strType,x->type.type); 2093 2093 psError(PS_ERR_BAD_PARAMETER_TYPE, … … 2115 2115 the spline fit functions require F32 and F64. 2116 2116 *****************************************************************************/ 2117 floatpsSpline1DEval(2118 floatx,2117 psF32 psSpline1DEval( 2118 psF32 x, 2119 2119 const psSpline1D *spline 2120 2120 ) … … 2172 2172 for (i=0;i<x->n;i++) { 2173 2173 tmpVector->data.F32[i] = psSpline1DEval( 2174 ( float) x->data.F64[i],2174 (psF32) x->data.F64[i], 2175 2175 spline 2176 2176 ); 2177 2177 } 2178 2178 } else { 2179 char* strType;2179 psS8* strType; 2180 2180 PS_TYPE_NAME(strType,x->type.type); 2181 2181 psError(PS_ERR_BAD_PARAMETER_TYPE, -
trunk/psLib/src/dataManip/psFunctions.h
r2600 r2788 12 12 * @author GLG, MHPCC 13 13 * 14 * @version $Revision: 1.3 7$ $Name: not supported by cvs2svn $15 * @date $Date: 2004-12- 02 21:12:51$14 * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ 15 * @date $Date: 2004-12-22 05:09:32 $ 16 16 * 17 17 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 39 39 * \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] 40 40 * 41 * @return floatvalue on the gaussian curve given the input parameters42 */ 43 floatpsGaussian(44 floatx, ///< Value at which to evaluate45 floatmean, ///< Mean for the Gaussian46 floatstddev, ///< Standard deviation for the Gaussian41 * @return psF32 value on the gaussian curve given the input parameters 42 */ 43 psF32 psGaussian( 44 psF32 x, ///< Value at which to evaluate 45 psF32 mean, ///< Mean for the Gaussian 46 psF32 stddev, ///< Standard deviation for the Gaussian 47 47 psBool normal ///< Indicates whether result should be normalized 48 48 ); … … 55 55 */ 56 56 psVector* p_psGaussianDev( 57 floatmean, ///< The mean of the Gaussian58 floatsigma, ///< The sigma of the Gaussian57 psF32 mean, ///< The mean of the Gaussian 58 psF32 sigma, ///< The sigma of the Gaussian 59 59 psS32 Npts ///< The size of the vector 60 60 ); … … 70 70 psPolynomialType type; ///< Polynomial type 71 71 psS32 n; ///< Number of terms 72 float*coeff; ///< Coefficients73 float*coeffErr; ///< Error in coefficients74 char*mask; ///< Coefficient mask72 psF32 *coeff; ///< Coefficients 73 psF32 *coeffErr; ///< Error in coefficients 74 psU8 *mask; ///< Coefficient mask 75 75 } 76 76 psPolynomial1D; … … 82 82 psS32 nX; ///< Number of terms in x 83 83 psS32 nY; ///< Number of terms in y 84 float**coeff; ///< Coefficients85 float**coeffErr; ///< Error in coefficients86 char**mask; ///< Coefficients mask84 psF32 **coeff; ///< Coefficients 85 psF32 **coeffErr; ///< Error in coefficients 86 psU8 **mask; ///< Coefficients mask 87 87 } 88 88 psPolynomial2D; … … 95 95 psS32 nY; ///< Number of terms in y 96 96 psS32 nZ; ///< Number of terms in z 97 float***coeff; ///< Coefficients98 float***coeffErr; ///< Error in coefficients99 char***mask; ///< Coefficients mask97 psF32 ***coeff; ///< Coefficients 98 psF32 ***coeffErr; ///< Error in coefficients 99 psU8 ***mask; ///< Coefficients mask 100 100 } 101 101 psPolynomial3D; … … 109 109 psS32 nY; ///< Number of terms in y 110 110 psS32 nZ; ///< Number of terms in z 111 float****coeff; ///< Coefficients112 float****coeffErr; ///< Error in coefficients113 char****mask; ///< Coefficients mask111 psF32 ****coeff; ///< Coefficients 112 psF32 ****coeffErr; ///< Error in coefficients 113 psU8 ****mask; ///< Coefficients mask 114 114 } 115 115 psPolynomial4D; … … 160 160 /** Evaluates a 1-D polynomial at specific coordinates. 161 161 * 162 * @return floatresult of polynomial at given location163 */ 164 floatpsPolynomial1DEval(162 * @return psF32 result of polynomial at given location 163 */ 164 psF32 psPolynomial1DEval( 165 165 const psPolynomial1D* myPoly, ///< Coefficients for the polynomial 166 floatx ///< location at which to evaluate166 psF32 x ///< location at which to evaluate 167 167 ); 168 168 169 169 /** Evaluates a 2-D polynomial at specific coordinates. 170 170 * 171 * @return floatresult of polynomial at given location172 */ 173 floatpsPolynomial2DEval(171 * @return psF32 result of polynomial at given location 172 */ 173 psF32 psPolynomial2DEval( 174 174 const psPolynomial2D* myPoly, ///< Coefficients for the polynomial 175 floatx, ///< x location at which to evaluate176 floaty ///< y location at which to evaluate175 psF32 x, ///< x location at which to evaluate 176 psF32 y ///< y location at which to evaluate 177 177 ); 178 178 179 179 /** Evaluates a 3-D polynomial at specific coordinates. 180 180 * 181 * @return floatresult of polynomial at given location182 */ 183 floatpsPolynomial3DEval(181 * @return psF32 result of polynomial at given location 182 */ 183 psF32 psPolynomial3DEval( 184 184 const psPolynomial3D* myPoly, ///< Coefficients for the polynomial 185 floatx, ///< x location at which to evaluate186 floaty, ///< y location at which to evaluate187 floatz ///< z location at which to evaluate185 psF32 x, ///< x location at which to evaluate 186 psF32 y, ///< y location at which to evaluate 187 psF32 z ///< z location at which to evaluate 188 188 ); 189 189 190 190 /** Evaluates a 4-D polynomial at specific coordinates. 191 191 * 192 * @return floatresult of polynomial at given location193 */ 194 floatpsPolynomial4DEval(192 * @return psF32 result of polynomial at given location 193 */ 194 psF32 psPolynomial4DEval( 195 195 const psPolynomial4D* myPoly, ///< Coefficients for the polynomial 196 floatw, ///< w location at which to evaluate197 floatx, ///< x location at which to evaluate198 floaty, ///< y location at which to evaluate199 floatz ///< z location at which to evaluate196 psF32 w, ///< w location at which to evaluate 197 psF32 x, ///< x location at which to evaluate 198 psF32 y, ///< y location at which to evaluate 199 psF32 z ///< z location at which to evaluate 200 200 ); 201 201 … … 235 235 psPolynomialType type; ///< Polynomial type 236 236 psS32 n; ///< Number of terms 237 double*coeff; ///< Coefficients238 double*coeffErr; ///< Error in coefficients239 char*mask; ///< Coefficient mask237 psF64 *coeff; ///< Coefficients 238 psF64 *coeffErr; ///< Error in coefficients 239 psU8 *mask; ///< Coefficient mask 240 240 } 241 241 psDPolynomial1D; … … 247 247 psS32 nX; ///< Number of terms in x 248 248 psS32 nY; ///< Number of terms in y 249 double**coeff; ///< Coefficients250 double**coeffErr; ///< Error in coefficients251 char**mask; ///< Coefficients mask249 psF64 **coeff; ///< Coefficients 250 psF64 **coeffErr; ///< Error in coefficients 251 psU8 **mask; ///< Coefficients mask 252 252 } 253 253 psDPolynomial2D; … … 260 260 psS32 nY; ///< Number of terms in y 261 261 psS32 nZ; ///< Number of terms in z 262 double***coeff; ///< Coefficients263 double***coeffErr; ///< Error in coefficients264 char***mask; ///< Coefficient mask262 psF64 ***coeff; ///< Coefficients 263 psF64 ***coeffErr; ///< Error in coefficients 264 psU8 ***mask; ///< Coefficient mask 265 265 } 266 266 psDPolynomial3D; … … 274 274 psS32 nY; ///< Number of terms in y 275 275 psS32 nZ; ///< Number of terms in z 276 double****coeff; ///< Coefficients277 double****coeffErr; ///< Error in coefficients278 char****mask; ///< Coefficients mask276 psF64 ****coeff; ///< Coefficients 277 psF64 ****coeffErr; ///< Error in coefficients 278 psU8 ****mask; ///< Coefficients mask 279 279 } 280 280 psDPolynomial4D; … … 324 324 /** Evaluates a double-precision 1-D polynomial at specific coordinates. 325 325 * 326 * @return floatresult of polynomial at given location327 */ 328 doublepsDPolynomial1DEval(326 * @return psF32 result of polynomial at given location 327 */ 328 psF64 psDPolynomial1DEval( 329 329 const psDPolynomial1D* myPoly, ///< Coefficients for the polynomial 330 doublex ///< Value at which to evaluate330 psF64 x ///< Value at which to evaluate 331 331 ); 332 332 333 333 /** Evaluates a double-precision 2-D polynomial at specific coordinates. 334 334 * 335 * @return floatresult of polynomial at given location336 */ 337 doublepsDPolynomial2DEval(335 * @return psF32 result of polynomial at given location 336 */ 337 psF64 psDPolynomial2DEval( 338 338 const psDPolynomial2D* myPoly, ///< Coefficients for the polynomial 339 doublex, ///< Value x at which to evaluate340 doubley ///< Value y at which to evaluate339 psF64 x, ///< Value x at which to evaluate 340 psF64 y ///< Value y at which to evaluate 341 341 ); 342 342 343 343 /** Evaluates a double-precision 3-D polynomial at specific coordinates. 344 344 * 345 * @return floatresult of polynomial at given location346 */ 347 doublepsDPolynomial3DEval(345 * @return psF32 result of polynomial at given location 346 */ 347 psF64 psDPolynomial3DEval( 348 348 const psDPolynomial3D* myPoly, ///< Coefficients for the polynomial 349 doublex, ///< Value x at which to evaluate350 doubley, ///< Value y at which to evaluate351 doublez ///< Value z at which to evaluate349 psF64 x, ///< Value x at which to evaluate 350 psF64 y, ///< Value y at which to evaluate 351 psF64 z ///< Value z at which to evaluate 352 352 ); 353 353 354 354 /** Evaluates a double-precision 4-D polynomial at specific coordinates. 355 355 * 356 * @return floatresult of polynomial at given location357 */ 358 doublepsDPolynomial4DEval(356 * @return psF32 result of polynomial at given location 357 */ 358 psF64 psDPolynomial4DEval( 359 359 const psDPolynomial4D* myPoly, ///< Coefficients for the polynomial 360 doublew, ///< Value w at which to evaluate361 doublex, ///< Value x at which to evaluate362 doubley, ///< Value y at which to evaluate363 doublez ///< Value z at which to evaluate360 psF64 w, ///< Value w at which to evaluate 361 psF64 x, ///< Value x at which to evaluate 362 psF64 y, ///< Value y at which to evaluate 363 psF64 z ///< Value z at which to evaluate 364 364 ); 365 365 … … 396 396 psS32 n; ///< The number of spline polynomials 397 397 psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials 398 float*p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1.399 float*domains; ///< The boundaries between each spline piece. Size is n+1.398 psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. 399 psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. 400 400 } 401 401 psSpline1D; … … 403 403 psSpline1D *psSpline1DAlloc(psS32 n, 404 404 psS32 order, 405 floatmin,406 floatmax);405 psF32 min, 406 psF32 max); 407 407 408 408 psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, 409 409 psS32 order); 410 410 411 floatpsSpline1DEval(412 floatx,411 psF32 psSpline1DEval( 412 psF32 x, 413 413 const psSpline1D *spline 414 414 ); … … 427 427 psScalar *x); 428 428 429 floatp_psNRSpline1DEval(psSpline1D *spline,429 psF32 p_psNRSpline1DEval(psSpline1D *spline, 430 430 const psVector* restrict x, 431 431 const psVector* restrict y, 432 floatX);432 psF32 X); 433 433 434 434 /* \} */// End of MathGroup Functions -
trunk/psLib/src/dataManip/psMinimize.c
r2741 r2788 9 9 * @author GLG, MHPCC 10 10 * 11 * @version $Revision: 1.9 5$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-12- 17 00:18:31$11 * @version $Revision: 1.96 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-12-22 05:09:32 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 15 15 * 16 16 * XXX: must follow coding name standards on local functions. 17 *18 * XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via19 * LMM) is not addressed here. We are waiting for subsequent SDRs20 * which will redefine the LMM functions.21 17 * 22 18 */ … … 64 60 XXX: Use a static vector. 65 61 *****************************************************************************/ 66 void psBuildSums1D( doublex,62 void psBuildSums1D(psF64 x, 67 63 psS32 polyOrder, 68 64 psVector* sums) 69 65 { 70 66 psS32 i = 0; 71 doublexSum = 0.0;67 psF64 xSum = 0.0; 72 68 73 69 if (sums == NULL) { … … 100 96 XXX: do an F64 version? 101 97 *****************************************************************************/ 102 float*CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices)98 psF32 *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) 103 99 const psVector* restrict y) ///< Coordinates 104 100 { … … 108 104 psS32 i; 109 105 psS32 k; 110 floatsig;111 floatp;106 psF32 sig; 107 psF32 p; 112 108 psS32 n = y->n; 113 float *u = (float *) psAlloc(n * sizeof(float));114 float *derivs2 = (float *) psAlloc(n * sizeof(float));115 float *X = (float*) & (x->data.F32[0]);116 float *Y = (float*) & (y->data.F32[0]);117 floatqn;109 psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); 110 psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); 111 psF32 *X = (psF32 *) & (x->data.F32[0]); 112 psF32 *Y = (psF32 *) & (y->data.F32[0]); 113 psF32 qn; 118 114 119 115 // XXX: The second derivatives at the endpoints, undefined in the SDR, … … 168 164 *****************************************************************************/ 169 165 /* 170 floatp_psNRSpline1DEval(psSpline1D *spline,166 psF32 p_psNRSpline1DEval(psSpline1D *spline, 171 167 const psVector* restrict x, 172 168 const psVector* restrict y, 173 floatX)169 psF32 X) 174 170 { 175 171 PS_PTR_CHECK_NULL(spline, NAN); … … 185 181 psS32 klo; 186 182 psS32 khi; 187 floatH;188 floatA;189 floatB;190 floatC;191 floatD;192 floatY;183 psF32 H; 184 psF32 A; 185 psF32 B; 186 psF32 C; 187 psF32 D; 188 psF32 Y; 193 189 194 190 n = spline->n; … … 258 254 "---- psVectorFitSpline1D() begin ----\n"); 259 255 psS32 numSplines = (y->n)-1; 260 floattmp;261 floatH;256 psF32 tmp; 257 psF32 H; 262 258 psS32 i; 263 floatslope;259 psF32 slope; 264 260 psVector *x32 = NULL; 265 261 psVector *y32 = NULL; … … 421 417 psTrace(".psLib.dataManip.psMinimize", 4, 422 418 "---- psMinimizeLMChi2Gauss1D() begin ----\n"); 423 floatx;424 inti;425 floatmean = params->data.F32[0];426 floatstdev = params->data.F32[1];419 psF32 x; 420 psS32 i; 421 psF32 mean = params->data.F32[0]; 422 psF32 stdev = params->data.F32[1]; 427 423 psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); 428 424 … … 443 439 for (i=0;i<coords->n;i++) { 444 440 x = ((psVector *) (coords->data[i]))->data.F32[0]; 445 floattmp = (x - mean) * psGaussian(x, mean, stdev, false);441 psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false); 446 442 deriv->data.F32[i][0] = tmp / (stdev * stdev); 447 443 tmp = (x - mean) * (x - mean) * … … 491 487 PS_PTR_CHECK_NULL(params, NULL); 492 488 493 doublenormalization = params->data.F32[0];494 doublex0 = params->data.F32[1];495 doubley0 = params->data.F32[2];496 doublesigmaX = params->data.F32[3];497 doublesigmaY = params->data.F32[4];498 doubletheta = params->data.F32[5];489 psF64 normalization = params->data.F32[0]; 490 psF64 x0 = params->data.F32[1]; 491 psF64 y0 = params->data.F32[2]; 492 psF64 sigmaX = params->data.F32[3]; 493 psF64 sigmaY = params->data.F32[4]; 494 psF64 theta = params->data.F32[5]; 499 495 psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); 500 496 … … 502 498 deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); 503 499 } else { 504 // XXX: Check size of derivative500 PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL); 505 501 } 506 502 … … 508 504 "---- psMinimizeLMChi2Gauss2D() begin ----\n"); 509 505 510 for ( inti=0;i<coords->n;i++) {511 doublex = ((psVector *) coords->data[i])->data.F32[0];512 doubley = ((psVector *) coords->data[i])->data.F32[0];513 514 doubleu = - (x-x0)*cos(theta) + (y-y0)*sin(theta);515 doublev = (x-x0)*cos(theta) + (y-y0)*sin(theta);516 517 doubleflux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +518 v*v/(2.0 * sigmaY * sigmaY)))/519 (2.0 * M_PI * sigmaX * sigmaY);506 for (psS32 i=0;i<coords->n;i++) { 507 psF64 x = ((psVector *) coords->data[i])->data.F32[0]; 508 psF64 y = ((psVector *) coords->data[i])->data.F32[0]; 509 510 psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); 511 psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta); 512 513 psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + 514 v*v/(2.0 * sigmaY * sigmaY)))/ 515 (2.0 * M_PI * sigmaX * sigmaY); 520 516 out->data.F32[i] = flux; 521 517 … … 584 580 PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL); 585 581 } 586 587 // XXX: Generate code that modifies covar matrix if not-NULL.588 582 if (covar != NULL) { 589 583 PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL); … … 610 604 psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); 611 605 psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); 612 613 // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));614 // for (i=0;i<numData;i++) {615 // deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);616 // }617 606 psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32); 618 619 psVector *currValueVec; 620 psVector *newValueVec; 621 622 float currChi2 = 0.0; 623 float newChi2 = 0.0; 624 float lamda = 0.00005; 607 psVector *currValueVec = NULL; 608 psVector *newValueVec = NULL; 609 psF32 currChi2 = 0.0; 610 psF32 newChi2 = 0.0; 611 psF32 lamda = 0.00005; 625 612 lamda = 0.05; 626 613 … … 711 698 if (j == k) { 712 699 A->data.F64[j][k] = 713 ( double) ((1.0 + lamda) * alpha->data.F32[j][k]);700 (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]); 714 701 } else { 715 A->data.F64[j][k] = ( double) alpha->data.F32[j][k];702 A->data.F64[j][k] = (psF64) alpha->data.F32[j][k]; 716 703 } 717 704 } … … 741 728 } else { 742 729 newParams->data.F32[i] = params->data.F32[i] - 743 ( float) paramDeltasF64->data.F64[i];730 (psF32) paramDeltasF64->data.F64[i]; 744 731 } 745 732 } … … 777 764 // We already masked params. 778 765 for (i=0;i<numParams;i++) { 779 params->data.F32[i] = ( float) newParams->data.F32[i];766 params->data.F32[i] = (psF32) newParams->data.F32[i]; 780 767 } 781 768 lamda*= 0.1; … … 826 813 psS32 k; 827 814 psS32 n = x->n; 828 doublefac;829 doublesum;815 psF64 fac; 816 psF64 sum; 830 817 PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64); 831 818 psScalar *fScalar; … … 836 823 // variable declarations. I retain them here to maintain coherence 837 824 // with the NR code. 838 doublemin = -1.0;839 doublemax = 1.0;840 doublebma = 0.5 * (max-min); // 1841 doublebpa = 0.5 * (max+min); // 0825 psF64 min = -1.0; 826 psF64 max = 1.0; 827 psF64 bma = 0.5 * (max-min); // 1 828 psF64 bpa = 0.5 * (max+min); // 0 842 829 843 830 // In this loop, we first calculate the values of X for which the … … 851 838 for (psS32 i=0;i<n;i++) { 852 839 // NR 5.8.4 853 double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));854 doubleX = (Y + bma + bpa) - 1.0;840 psF64 Y = cos(M_PI * (0.5 + ((psF32) i)) / ((psF32) n)); 841 psF64 X = (Y + bma + bpa) - 1.0; 855 842 tmpScalar.data.F64 = X; 856 843 … … 873 860 // coefficients of the Chebyshev polynomial: NR 5.8.7. 874 861 875 fac = 2.0/(( float) n);862 fac = 2.0/((psF32) n); 876 863 // XXX: is this loop bound correct? 877 864 for (j=0;j<myPoly->n;j++) { … … 879 866 for (k=0;k<n;k++) { 880 867 sum+= f->data.F64[k] * 881 cos(M_PI * (( float) j) * (0.5 + ((float) k)) / ((float) n));868 cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n)); 882 869 } 883 870 … … 1075 1062 *****************************************************************************/ 1076 1063 psMinimization *psMinimizationAlloc(psS32 maxIter, 1077 floattol)1064 psF32 tol) 1078 1065 { 1079 1066 PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL); … … 1092 1079 // LINE to it. We assume BASEMASK is non-null. 1093 1080 #define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ 1094 for ( inti=0;i<BASE->n;i++) { \1081 for (psS32 i=0;i<BASE->n;i++) { \ 1095 1082 if (BASEMASK->data.U8[i] == 0) { \ 1096 1083 OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ … … 1102 1089 #define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \ 1103 1090 BOOL_VAR = true; \ 1104 for ( inti=0;i<IN->n;i++) { \1091 for (psS32 i=0;i<IN->n;i++) { \ 1105 1092 if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \ 1106 1093 BOOL_VAR = false; \ … … 1111 1098 #define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \ 1112 1099 BOOL_VAR = true; \ 1113 for ( inti=0;i<IN->n;i++) { \1100 for (psS32 i=0;i<IN->n;i++) { \ 1114 1101 if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \ 1115 1102 BOOL_VAR = false; \ … … 1152 1139 psMinimizePowellFunc func) 1153 1140 { 1154 floata = 0.0;1155 floatb = 0.0;1156 floatc = 0.0;1157 floatfa = 0.0;1158 floatfb = 0.0;1159 floatfc = 0.0;1141 psF32 a = 0.0; 1142 psF32 b = 0.0; 1143 psF32 c = 0.0; 1144 psF32 fa = 0.0; 1145 psF32 fb = 0.0; 1146 psF32 fc = 0.0; 1160 1147 psS32 iter = 100; 1161 floataDir = 0.0;1162 floatcDir = 0.0;1163 floatnew_aDir = 0.0;1164 floatnew_cDir = 0.0;1148 psF32 aDir = 0.0; 1149 psF32 cDir = 0.0; 1150 psF32 new_aDir = 0.0; 1151 psF32 new_cDir = 0.0; 1165 1152 psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); 1166 floatstepSize = PS_DETERMINE_BRACKET_STEP_SIZE;1153 psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; 1167 1154 psVector *tmp = NULL; 1168 1155 psBool boolLineIsNull = true; … … 1301 1288 psMinimizePowellFunc func) 1302 1289 { 1303 floata = 0.0;1304 floatb = 0.0;1305 floatc = 0.0;1306 floatfa = 0.0;1307 floatfb = 0.0;1308 floatfc = 0.0;1290 psF32 a = 0.0; 1291 psF32 b = 0.0; 1292 psF32 c = 0.0; 1293 psF32 fa = 0.0; 1294 psF32 fb = 0.0; 1295 psF32 fc = 0.0; 1309 1296 psS32 iter = 0; 1310 1297 PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32); 1311 1298 psBool boolLineIsNull = true; 1312 floatprevMin = 0.0;1313 intcountMin = 0;1299 psF32 prevMin = 0.0; 1300 psS32 countMin = 0; 1314 1301 1315 1302 psTrace(".psLib.dataManip.p_psDetermineBracket", 4, … … 1403 1390 *****************************************************************************/ 1404 1391 #define PS_LINEMIN_MAX_ITERATIONS 30 1405 floatp_psLineMin(psMinimization *min,1392 psF32 p_psLineMin(psMinimization *min, 1406 1393 psVector *params, 1407 1394 psVector *line, … … 1423 1410 PS_PTR_CHECK_NULL(func, NAN); 1424 1411 psVector *bracket; 1425 floata = 0.0;1426 floatb = 0.0;1427 floatc = 0.0;1428 floatn = 0.0;1429 floatfa = 0.0;1430 floatfb = 0.0;1431 floatfc = 0.0;1432 floatfn = 0.0;1433 floatmul = 0.0;1412 psF32 a = 0.0; 1413 psF32 b = 0.0; 1414 psF32 c = 0.0; 1415 psF32 n = 0.0; 1416 psF32 fa = 0.0; 1417 psF32 fb = 0.0; 1418 psF32 fc = 0.0; 1419 psF32 fn = 0.0; 1420 psF32 mul = 0.0; 1434 1421 PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32); 1435 1422 PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32); … … 1544 1531 This routine must minimize a possibly multi-dimensional function. The 1545 1532 function to be minimized "func" is: 1546 floatfunc(psVector *params, psArray *coords)1533 psF32 func(psVector *params, psArray *coords) 1547 1534 The "params" are the parameters of the function which are varied. The data 1548 1535 points at which the function is varied are in the argument "coords" which is … … 1578 1565 psVector *myParamMask = NULL; 1579 1566 psMinimization dummyMin; 1580 floatmul = 0.0;1581 floatbaseFuncVal = 0.0;1582 floatcurrFuncVal = 0.0;1567 psF32 mul = 0.0; 1568 psF32 baseFuncVal = 0.0; 1569 psF32 currFuncVal = 0.0; 1583 1570 psS32 biggestIter = 0; 1584 floatbiggestDiff = 0.0;1585 intiterationNumber = 0;1571 psF32 biggestDiff = 0.0; 1572 psS32 iterationNumber = 0; 1586 1573 1587 1574 psTrace(".psLib.dataManip.psMinimizePowell", 4, … … 1711 1698 } 1712 1699 } 1713 floatfqp = func(pQP, coords);1714 floatterm1 = (baseFuncVal - currFuncVal) - biggestDiff;1700 psF32 fqp = func(pQP, coords); 1701 psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff; 1715 1702 term1*= term1; 1716 1703 term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp); 1717 floatterm2 = baseFuncVal - fqp;1704 psF32 term2 = baseFuncVal - fqp; 1718 1705 term2*= term2 * biggestDiff; 1719 1706 if (term1 < term2) { … … 1760 1747 PS_PTR_CHECK_NULL(params, NULL); 1761 1748 1762 floatx;1763 inti;1764 floatmean = params->data.F32[0];1765 floatstdev = params->data.F32[1];1749 psF32 x; 1750 psS32 i; 1751 psF32 mean = params->data.F32[0]; 1752 psF32 stdev = params->data.F32[1]; 1766 1753 psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); 1767 1754 … … 1783 1770 XXX: This is F32 only 1784 1771 *****************************************************************************/ 1785 floatmyPowellChi2Func(const psVector *params,1772 psF32 myPowellChi2Func(const psVector *params, 1786 1773 const psArray *coords) 1787 1774 { … … 1794 1781 PS_PTR_CHECK_NULL(coords, NAN); 1795 1782 1796 floatchi2 = 0.0;1797 floatd;1783 psF32 chi2 = 0.0; 1784 psF32 d; 1798 1785 psS32 i; 1799 1786 psVector *tmp; -
trunk/psLib/src/dataManip/psMinimize.h
r2250 r2788 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.3 6$ $Name: not supported by cvs2svn $11 * @date $Date: 2004-1 1-01 23:57:08$10 * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2004-12-22 05:09:32 $ 12 12 * 13 13 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 40 40 { 41 41 psS32 maxIter; ///< Convergence limit 42 floattol; ///< Error Tolerance43 floatvalue; ///< Value of function at minimum42 psF32 tol; ///< Error Tolerance 43 psF32 value; ///< Value of function at minimum 44 44 psS32 iter; ///< Number of iterations required 45 floatlastDelta; ///< The last difference for the fit45 psF32 lastDelta; ///< The last difference for the fit 46 46 } 47 47 psMinimization; 48 48 49 49 psMinimization *psMinimizationAlloc(psS32 maxIter, 50 floattol);50 psF32 tol); 51 51 52 52 /** Derive a polynomial fit. … … 91 91 92 92 typedef 93 float(*psMinimizePowellFunc)(const psVector *params,93 psF32 (*psMinimizePowellFunc)(const psVector *params, 94 94 const psArray *coords); 95 95 -
trunk/psLib/src/dataManip/psStats.c
r2782 r2788 9 9 * @author GLG, MHPCC 10 10 * 11 * @version $Revision: 1.10 7$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-12-22 0 0:54:28$11 * @version $Revision: 1.108 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-12-22 05:09:32 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 70 70 /*****************************************************************************/ 71 71 72 psBool p_psGetStatValue(const psStats* stats, double*value)72 psBool p_psGetStatValue(const psStats* stats, psF64 *value) 73 73 { 74 74 … … 145 145 this routine sets stats->sampleMean to NAN. 146 146 *****************************************************************************/ 147 intp_psVectorSampleMean(const psVector* restrict myVector,148 const psVector* restrict errors,149 const psVector* restrict maskVector,150 psU32 maskVal,151 psStats* stats)147 psS32 p_psVectorSampleMean(const psVector* restrict myVector, 148 const psVector* restrict errors, 149 const psVector* restrict maskVector, 150 psU32 maskVal, 151 psStats* stats) 152 152 { 153 153 154 154 psS32 i = 0; // Loop index variable 155 floatmean = 0.0; // The mean155 psF32 mean = 0.0; // The mean 156 156 psS32 count = 0; // # of points in this mean 157 157 … … 171 171 } 172 172 if (count != 0) { 173 mean /= ( float)count;173 mean /= (psF32)count; 174 174 } else { 175 175 mean = NAN; … … 185 185 } 186 186 if (count != 0) { 187 mean /= ( float)count;187 mean /= (psF32)count; 188 188 } else { 189 189 mean = NAN; … … 199 199 } 200 200 if (count != 0) { 201 mean /= ( float)count;201 mean /= (psF32)count; 202 202 } else { 203 203 mean = NAN; … … 207 207 mean += myVector->data.F32[i]; 208 208 } 209 mean /= ( float)myVector->n;209 mean /= (psF32)myVector->n; 210 210 } 211 211 } … … 283 283 this routine sets stats->max to NAN. 284 284 *****************************************************************************/ 285 intp_psVectorMax(const psVector* restrict myVector,286 const psVector* restrict maskVector,287 psU32 maskVal,288 psStats* stats)285 psS32 p_psVectorMax(const psVector* restrict myVector, 286 const psVector* restrict maskVector, 287 psU32 maskVal, 288 psStats* stats) 289 289 { 290 290 psS32 i = 0; // Loop index variable 291 floatmax = -PS_MAX_F32; // The calculated maximum291 psF32 max = -PS_MAX_F32; // The calculated maximum 292 292 psS32 empty = true; // Does this vector have valid elements? 293 293 … … 348 348 this routine sets stats->min to NAN. 349 349 *****************************************************************************/ 350 intp_psVectorMin(const psVector* restrict myVector,351 const psVector* restrict maskVector,352 psU32 maskVal,353 psStats* stats)350 psS32 p_psVectorMin(const psVector* restrict myVector, 351 const psVector* restrict maskVector, 352 psU32 maskVal, 353 psStats* stats) 354 354 { 355 355 psS32 i = 0; // Loop index variable 356 floatmin = PS_MAX_F32; // The calculated maximum356 psF32 min = PS_MAX_F32; // The calculated maximum 357 357 psS32 empty = true; // Does this vector have valid elements? 358 358 … … 612 612 *****************************************************************************/ 613 613 psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, 614 floatsigma)614 psF32 sigma) 615 615 { 616 616 PS_PTR_CHECK_NULL(robustHistogram, NULL); … … 619 619 psS32 i = 0; // Loop index variable 620 620 psS32 j = 0; // Loop index variable 621 floatiMid;622 floatjMid;621 psF32 iMid; 622 psF32 jMid; 623 623 psS32 numBins = robustHistogram->nums->n; 624 624 psS32 numBounds = robustHistogram->bounds->n; … … 626 626 psS32 jMin = 0; 627 627 psS32 jMax = 0; 628 floatfirstBound = robustHistogram->bounds->data.F32[0];629 floatlastBound = robustHistogram->bounds->data.F32[numBounds-1];628 psF32 firstBound = robustHistogram->bounds->data.F32[0]; 629 psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; 630 630 psScalar x; 631 631 … … 795 795 psS32 i = 0; // Loop index variable 796 796 psS32 countInt = 0; // # of data points being used 797 floatcountFloat = 0.0; // # of data points being used798 floatmean = 0.0; // The mean799 floatdiff = 0.0; // Used in calculating stdev800 floatsumSquares = 0.0; // temporary variable801 floatsumDiffs = 0.0; // temporary variable797 psF32 countFloat = 0.0; // # of data points being used 798 psF32 mean = 0.0; // The mean 799 psF32 diff = 0.0; // Used in calculating stdev 800 psF32 sumSquares = 0.0; // temporary variable 801 psF32 sumDiffs = 0.0; // temporary variable 802 802 803 803 // This procedure requires the mean. If it has not been already … … 865 865 psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", countInt); 866 866 } else { 867 countFloat = ( float)countInt;867 countFloat = (psF32)countInt; 868 868 stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); 869 869 } … … 889 889 psS32 i = 0; // Loop index variable 890 890 psS32 countInt = 0; // # of data points being used 891 floatcountFloat = 0.0; // # of data points being used892 floatmean = 0.0; // The mean893 floatdiff = 0.0; // Used in calculating stdev894 floatsumSquares = 0.0; // temporary variable895 floatsumDiffs = 0.0; // temporary variable891 psF32 countFloat = 0.0; // # of data points being used 892 psF32 mean = 0.0; // The mean 893 psF32 diff = 0.0; // Used in calculating stdev 894 psF32 sumSquares = 0.0; // temporary variable 895 psF32 sumDiffs = 0.0; // temporary variable 896 896 // psF32 sum1; 897 897 // psF32 sum2; … … 982 982 stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor)); 983 983 } else { 984 countFloat = ( float)countInt;984 countFloat = (psF32)countInt; 985 985 stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); 986 986 … … 1003 1003 -2: warning 1004 1004 *****************************************************************************/ 1005 intp_psVectorClippedStats(const psVector* restrict myVector,1006 const psVector* restrict errors,1007 const psVector* restrict maskVector,1008 psU32 maskVal,1009 psStats* stats)1005 psS32 p_psVectorClippedStats(const psVector* restrict myVector, 1006 const psVector* restrict errors, 1007 const psVector* restrict maskVector, 1008 psU32 maskVal, 1009 psStats* stats) 1010 1010 { 1011 1011 psS32 i = 0; // Loop index variable 1012 1012 psS32 j = 0; // Loop index variable 1013 floatclippedMean = 0.0; // self-explanatory1014 floatclippedStdev = 0.0; // self-explanatory1015 floatoldStanMean = 0.0; // Temporary variable1016 floatoldStanStdev = 0.0; // Temporary variable1013 psF32 clippedMean = 0.0; // self-explanatory 1014 psF32 clippedStdev = 0.0; // self-explanatory 1015 psF32 oldStanMean = 0.0; // Temporary variable 1016 psF32 oldStanStdev = 0.0; // Temporary variable 1017 1017 psVector* tmpMask = NULL; // Temporary vector 1018 1018 … … 1240 1240 XXX: Create a 2nd-order polynomial version and solve for X analytically. 1241 1241 *****************************************************************************/ 1242 floatp_ps1DPolyMedian(psPolynomial1D* myPoly,1243 floatrangeLow,1244 floatrangeHigh,1245 floatgetThisValue)1242 psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly, 1243 psF32 rangeLow, 1244 psF32 rangeHigh, 1245 psF32 getThisValue) 1246 1246 { 1247 1247 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1250 1250 // falls within the range of y-values of the polynomial "myPoly" in the 1251 1251 // specified x-range (rangeLow:rangeHigh). 1252 floatfLo = psPolynomial1DEval(1252 psF32 fLo = psPolynomial1DEval( 1253 1253 myPoly, 1254 1254 rangeLow 1255 1255 ); 1256 floatfHi = psPolynomial1DEval(1256 psF32 fHi = psPolynomial1DEval( 1257 1257 myPoly, 1258 1258 rangeHigh … … 1266 1266 1267 1267 psS32 numIterations = 0; 1268 floatmidpoint = 0.0;1269 floatoldMidpoint = 1.0;1270 floatf = 0.0;1268 psF32 midpoint = 0.0; 1269 psF32 oldMidpoint = 1.0; 1270 psF32 f = 0.0; 1271 1271 1272 1272 while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) { … … 1309 1309 tests to ensure that binNum is within acceptable ranges for both vectors. 1310 1310 *****************************************************************************/ 1311 floatfitQuadraticSearchForYThenReturnX(psVector *xVec,1311 psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, 1312 1312 psVector *yVec, 1313 1313 psS32 binNum, 1314 floatyVal)1314 psF32 yVal) 1315 1315 { 1316 1316 PS_VECTOR_CHECK_NULL(xVec, NAN); … … 1331 1331 psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); 1332 1332 1333 floattmpFloat;1333 psF32 tmpFloat; 1334 1334 1335 1335 if ((binNum > 0) && (binNum < (yVec->n - 2))) { 1336 1336 // The general case. We have all three points. 1337 x->data.F64[0] = ( double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));1338 x->data.F64[1] = ( double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));1339 x->data.F64[2] = ( double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));1337 x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); 1338 x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); 1339 x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); 1340 1340 y->data.F64[0] = yVec->data.F32[binNum - 1]; 1341 1341 y->data.F64[1] = yVec->data.F32[binNum]; … … 1346 1346 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 1347 1347 PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, 1348 ( double)yVal,y->data.F64[2],y->data.F64[0]);1348 (psF64)yVal,y->data.F64[2],y->data.F64[0]); 1349 1349 } 1350 1350 yErr->data.F64[0] = 1.0; … … 1425 1425 XXX: Check for errors in psLib routines that we call. 1426 1426 *****************************************************************************/ 1427 intp_psVectorRobustStats(const psVector* restrict myVector,1428 const psVector* restrict errors,1429 const psVector* restrict maskVector,1430 psU32 maskVal,1431 psStats* stats)1427 psS32 p_psVectorRobustStats(const psVector* restrict myVector, 1428 const psVector* restrict errors, 1429 const psVector* restrict maskVector, 1430 psU32 maskVal, 1431 psStats* stats) 1432 1432 { 1433 1433 psHistogram* robustHistogram = NULL; 1434 1434 psVector* robustHistogramVector = NULL; 1435 floatbinSize = 0.0; // Size of the histogram bins1435 psF32 binSize = 0.0; // Size of the histogram bins 1436 1436 psS32 LQBinNum = -1; // Bin num for lower quartile 1437 1437 psS32 UQBinNum = -1; // Bin num for upper quartile … … 1439 1439 psS32 i = 0; // Loop index variable 1440 1440 psS32 modeBinNum = 0; 1441 floatmodeBinCount = 0.0;1442 floatdL = 0.0;1441 psF32 modeBinCount = 0.0; 1442 psF32 dL = 0.0; 1443 1443 psS32 numBins = 0; 1444 floatmyMean = 0.0;1445 floatmyStdev = 0.0;1446 floatcountFloat = 0.0;1447 floatdiff = 0.0;1448 floatsumSquares = 0.0;1449 floatsumDiffs = 0.0;1444 psF32 myMean = 0.0; 1445 psF32 myStdev = 0.0; 1446 psF32 countFloat = 0.0; 1447 psF32 diff = 0.0; 1448 psF32 sumSquares = 0.0; 1449 psF32 sumDiffs = 0.0; 1450 1450 psVector* cumulativeRobustSums = NULL; 1451 floatsumRobust = 0.0;1452 floatsumN50 = 0.0;1453 floatsumNfit = 0.0;1451 psF32 sumRobust = 0.0; 1452 psF32 sumN50 = 0.0; 1453 psF32 sumNfit = 0.0; 1454 1454 psScalar tmpScalar; 1455 1455 tmpScalar.type.type = PS_TYPE_F32; … … 1460 1460 // that by 10.0; 1461 1461 //XXX: add errors 1462 intrc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats);1462 psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); 1463 1463 if (rc != 0) { 1464 1464 psError(PS_ERR_UNEXPECTED_NULL, … … 1614 1614 for (i=0;i<robustHistogramVector->n;i++) { 1615 1615 myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); 1616 ((psVector *) (myCoords->data[i]))->data.F32[0] = ( float) i;1616 ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i; 1617 1617 y->data.F32[i] = robustHistogramVector->data.F32[i]; 1618 1618 } … … 1770 1770 The histogram structure 1771 1771 *****************************************************************************/ 1772 psHistogram* psHistogramAlloc( float lower, floatupper, psS32 n)1772 psHistogram* psHistogramAlloc(psF32 lower, psF32 upper, psS32 n) 1773 1773 { 1774 1774 PS_INT_CHECK_POSITIVE(n, NULL); … … 1777 1777 psS32 i = 0; // Loop index variable 1778 1778 psHistogram* newHist = NULL; // The new histogram structure 1779 floatbinSize = 0.0; // The histogram bin size1779 psF32 binSize = 0.0; // The histogram bin size 1780 1780 1781 1781 // Allocate memory for the new histogram structure. If there are N … … 1787 1787 1788 1788 // Calculate the bounds for each bin. 1789 binSize = (upper - lower) / ( float)n;1789 binSize = (upper - lower) / (psF32)n; 1790 1790 // XXX: Is the following necessary? It prevents the max data point 1791 1791 // from being in a non-existant bin. 1792 1792 binSize += FLT_EPSILON; 1793 1793 for (i = 0; i < n + 1; i++) { 1794 newHist->bounds->data.F32[i] = lower + (binSize * ( float)i);1794 newHist->bounds->data.F32[i] = lower + (binSize * (psF32)i); 1795 1795 } 1796 1796 … … 1857 1857 } 1858 1858 1859 /***************************************************************************** 1860 UpdateHistogramBins(binNum, out, data, error): This routine is to be used when 1861 updating the histogram in the presence of errors in the input data. We treat 1862 the data point as a boxcar PDF and update a range of points surrounding the 1863 histogram bin which contains the point. The width of that boxcar is defined 1864 as 2.35 * error. Inputs: 1865 binNum: the bin number of the data point in the histogram 1866 out: the histogram structure 1867 data: the data point value 1868 error: the error in that data point 1869 1870 XXX: Must test this. 1871 *****************************************************************************/ 1859 1872 psS32 UpdateHistogramBins(psS32 binNum, 1860 1873 psHistogram* out, … … 1863 1876 { 1864 1877 PS_PTR_CHECK_NULL(out, -1); 1865 PS_INT_CHECK_RANGE(binNum, 0, out->nums->n-1, -2); 1866 /* 1867 psF32 width = 2.35 * error; 1868 psF32 centerBinWidth = out->bounds->data.F32[binNum+1] - out->bounds->data.F32[binNum] 1869 psF32 boxcarCenter = (out->bounds->data.F32[binNum] + out->bounds->data.F32[binNum+1]) / 2.0; 1870 1871 if (width <= centerBinWidth) { 1872 out->nums->data.F32[binNum]+= 1.0; 1873 } else { 1874 out->nums->data.F32[binNum]+= centerBinWidth / width; 1875 // XXX: walk to the left, adding fractional values. 1876 // XXX: walk to the right, adding fractional values. 1877 1878 1879 } 1880 */ 1878 PS_PTR_CHECK_NULL(out->bounds, -1); 1879 PS_PTR_CHECK_NULL(out->nums, -1); 1880 PS_INT_CHECK_RANGE(binNum, 0, ((out->nums->n)-1), -2); 1881 PS_FLOAT_COMPARE(0.0, error, -3); 1882 PS_FLOAT_CHECK_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4); 1883 1884 psF32 boxcarWidth = 2.35 * error; 1885 psF32 boxcarCenter = (out->bounds->data.F32[binNum] + 1886 out->bounds->data.F32[binNum+1]) / 2.0; 1887 psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); 1888 psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); 1889 psS32 bin; 1890 psS32 boxcarLeftBinNum; 1891 psS32 boxcarRightBinNum; 1892 1893 // Determine the left endpoint of the boxcar for the PDF. 1894 for (bin=binNum ; bin >= 0 ; bin--) { 1895 if (out->nums->data.F32[bin] <= boxcarLeft) { 1896 boxcarLeftBinNum = bin; 1897 break; 1898 } 1899 } 1900 1901 // Determine the right endpoint of the boxcar for the PDF. 1902 for (bin=binNum ; bin < out->nums->n ; bin++) { 1903 if (out->nums->data.F32[bin] >= boxcarRight) { 1904 boxcarRightBinNum = bin; 1905 break; 1906 } 1907 } 1908 1909 // 1910 // If the boxcar fits entirely inside this bin, then simply add 1.0 to the 1911 // bin and return. 1912 // 1913 if (boxcarLeftBinNum == boxcarRightBinNum) { 1914 out->nums->data.F32[binNum]+= 1.0; 1915 return(0); 1916 } 1917 1918 // 1919 // If we get here, multiple bins must be updated. We handle the left 1920 // endpoint, and right endpoint differently. 1921 // 1922 out->nums->data.F32[boxcarLeftBinNum]+= 1923 (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth; 1924 1925 // 1926 // Loop through the center bins, if any. 1927 // 1928 for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) { 1929 out->nums->data.F32[bin]+= 1930 (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth; 1931 } 1932 1933 // 1934 // Handle the right endpoint differently. 1935 // 1936 out->nums->data.F32[boxcarRightBinNum]+= 1937 (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth; 1938 1939 // 1940 // Return 0 on success. 1941 // 1881 1942 return(0); 1882 1943 } … … 1906 1967 { 1907 1968 PS_PTR_CHECK_NULL(out, NULL); 1969 PS_VECTOR_CHECK_NULL(out->bounds, NULL); 1908 1970 PS_VECTOR_CHECK_TYPE(out->bounds, PS_TYPE_F32, NULL); 1971 PS_INT_CHECK_NON_NEGATIVE(out->bounds->n, NULL); 1972 PS_VECTOR_CHECK_NULL(out->nums, NULL); 1909 1973 PS_VECTOR_CHECK_TYPE(out->nums, PS_TYPE_F32, NULL); 1910 1974 PS_INT_CHECK_NON_NEGATIVE(out->nums->n, NULL); … … 1920 1984 1921 1985 psS32 i = 0; // Loop index variable 1922 floatbinSize = 0.0; // Histogram bin size1986 psF32 binSize = 0.0; // Histogram bin size 1923 1987 psS32 binNum = 0; // A temporary bin number 1924 1988 psS32 numBins = 0; // The total number of bins 1925 psS32 tmp = 0;1926 1989 psScalar tmpScalar; 1927 1990 tmpScalar.type.type = PS_TYPE_F32; 1928 psVector* inF32; 1929 psS32 mustFreeTmp = 1; 1991 psVector* inF32 = NULL; 1992 psVector* errorsF32 = NULL; 1993 psS32 mustFreeVectorIn = 1; 1994 psS32 mustFreeVectorErrors = 1; 1995 1996 // Convert input and errors vectors to F32 if necessary. 1930 1997 inF32 = p_psConvertToF32((psVector *) in); 1931 1998 if (inF32 == NULL) { 1932 1999 inF32 = (psVector *) in; 1933 mustFreeTmp = 0; 2000 mustFreeVectorIn = 0; 2001 } 2002 errorsF32 = p_psConvertToF32((psVector *) errors); 2003 if (errorsF32 == NULL) { 2004 errorsF32 = (psVector *) errors; 2005 mustFreeVectorErrors = 0; 1934 2006 } 1935 2007 … … 1950 2022 binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; 1951 2023 binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); 1952 if (errors != NULL) {2024 if (errorsF32 != NULL) { 1953 2025 // XXX: Check return codes. 1954 2026 UpdateHistogramBins(binNum, out, 1955 2027 inF32->data.F32[i], 1956 errors ->data.F32[i]);2028 errorsF32->data.F32[i]); 1957 2029 } else { 1958 2030 // XXX: This if-statement really shouldn't be necessary. … … 1969 2041 // correct bin number requires a bit more work. 1970 2042 tmpScalar.data.F32 = inF32->data.F32[i]; 1971 tmp= p_psVectorBinDisect(out->bounds, &tmpScalar);1972 if ( tmp< 0) {2043 binNum = p_psVectorBinDisect(out->bounds, &tmpScalar); 2044 if (binNum < 0) { 1973 2045 psLogMsg(__func__, PS_LOG_WARN, 1974 2046 "WARNING: psVectorHistogram(): element outside histogram bounds.\n"); 1975 2047 } else { 1976 if (errors != NULL) {2048 if (errorsF32 != NULL) { 1977 2049 // XXX: Check return codes. 1978 UpdateHistogramBins( tmp, out,2050 UpdateHistogramBins(binNum, out, 1979 2051 inF32->data.F32[i], 1980 2052 errors->data.F32[i]); 1981 2053 } else { 1982 (out->nums->data.F32[ tmp])+= 1.0;2054 (out->nums->data.F32[binNum])+= 1.0; 1983 2055 } 1984 2056 } … … 1988 2060 } 1989 2061 1990 if (mustFree Tmp== 1) {2062 if (mustFreeVectorIn == 1) { 1991 2063 psFree(inF32); 2064 } 2065 if (mustFreeVectorErrors == 1) { 2066 psFree(errorsF32); 1992 2067 } 1993 2068 return (out); … … 2015 2090 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2016 2091 for (i = 0; i < in->n; i++) { 2017 tmp->data.F32[i] = (float)in->data.S8[i]; 2092 tmp->data.F32[i] = (psF32)in->data.S8[i]; 2093 } 2094 } else if (in->type.type == PS_TYPE_S16) { 2095 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2096 for (i = 0; i < in->n; i++) { 2097 tmp->data.F32[i] = (psF32) in->data.S16[i]; 2098 } 2099 } else if (in->type.type == PS_TYPE_S32) { 2100 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2101 for (i = 0; i < in->n; i++) { 2102 tmp->data.F32[i] = (psF32)in->data.S32[i]; 2103 } 2104 } else if (in->type.type == PS_TYPE_S64) { 2105 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2106 for (i = 0; i < in->n; i++) { 2107 tmp->data.F32[i] = (psF32)in->data.S64[i]; 2108 } 2109 } else if (in->type.type == PS_TYPE_U8) { 2110 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2111 for (i = 0; i < in->n; i++) { 2112 tmp->data.F32[i] = (psF32)in->data.U8[i]; 2018 2113 } 2019 2114 } else if (in->type.type == PS_TYPE_U16) { 2020 2115 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2021 2116 for (i = 0; i < in->n; i++) { 2022 tmp->data.F32[i] = ( float)in->data.U16[i];2023 } 2024 } else if (in->type.type == PS_TYPE_U 8) {2117 tmp->data.F32[i] = (psF32)in->data.U16[i]; 2118 } 2119 } else if (in->type.type == PS_TYPE_U32) { 2025 2120 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2026 2121 for (i = 0; i < in->n; i++) { 2027 tmp->data.F32[i] = (float)in->data.U8[i]; 2122 tmp->data.F32[i] = (psF32)in->data.U32[i]; 2123 } 2124 } else if (in->type.type == PS_TYPE_U64) { 2125 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2126 for (i = 0; i < in->n; i++) { 2127 tmp->data.F32[i] = (psF32)in->data.U64[i]; 2028 2128 } 2029 2129 } else if (in->type.type == PS_TYPE_F64) { 2030 2130 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2031 2131 for (i = 0; i < in->n; i++) { 2032 tmp->data.F32[i] = ( float)in->data.F64[i];2132 tmp->data.F32[i] = (psF32)in->data.F64[i]; 2033 2133 } 2034 2134 } else if (in->type.type == PS_TYPE_F32) { 2035 2135 // do nothing 2036 2136 } else { 2037 char* strType;2038 PS_TYPE_NAME(strType, in->type.type);2137 psS8* strType; 2138 PS_TYPE_NAME(strType, in->type.type); 2039 2139 psError(PS_ERR_BAD_PARAMETER_TYPE, true, 2040 2140 PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED, … … 2074 2174 } 2075 2175 2076 psVector* inF32 ;2077 psVector* errorsF32 ;2176 psVector* inF32 = NULL; 2177 psVector* errorsF32 = NULL; 2078 2178 psS32 mustFreeVectorIn = 1; 2079 2179 psS32 mustFreeVectorErrors = 1; … … 2178 2278 return (stats); 2179 2279 } 2280 -
trunk/psLib/src/dataManip/psStats.h
r2778 r2788 10 10 * @author George Gusciora, MHPCC 11 11 * 12 * @version $Revision: 1.3 6$ $Name: not supported by cvs2svn $13 * @date $Date: 2004-12-2 1 20:42:07$12 * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2004-12-22 05:09:32 $ 14 14 * 15 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 56 56 typedef struct 57 57 { 58 doublesampleMean; ///< formal mean of sample59 doublesampleMedian; ///< formal median of sample60 doublesampleStdev; ///< standard deviation of sample61 doublesampleUQ; ///< upper quartile of sample62 doublesampleLQ; ///< lower quartile of sample63 doublerobustMean; ///< robust mean of array64 doublerobustMedian; ///< robust median of array65 doublerobustMode; ///< Robust mode of array66 doublerobustStdev; ///< robust standard deviation of array67 doublerobustUQ; ///< robust upper quartile68 doublerobustLQ; ///< robust lower quartile58 psF64 sampleMean; ///< formal mean of sample 59 psF64 sampleMedian; ///< formal median of sample 60 psF64 sampleStdev; ///< standard deviation of sample 61 psF64 sampleUQ; ///< upper quartile of sample 62 psF64 sampleLQ; ///< lower quartile of sample 63 psF64 robustMean; ///< robust mean of array 64 psF64 robustMedian; ///< robust median of array 65 psF64 robustMode; ///< Robust mode of array 66 psF64 robustStdev; ///< robust standard deviation of array 67 psF64 robustUQ; ///< robust upper quartile 68 psF64 robustLQ; ///< robust lower quartile 69 69 psS32 robustN50; ///< 70 70 psS32 robustNfit; ///< 71 doubleclippedMean; ///< Nsigma clipped mean72 doubleclippedStdev; ///< standard deviation after clipping71 psF64 clippedMean; ///< Nsigma clipped mean 72 psF64 clippedStdev; ///< standard deviation after clipping 73 73 psS32 clippedNvalues; ///< ??? 74 doubleclipSigma; ///< Nsigma used for clipping; user input74 psF64 clipSigma; ///< Nsigma used for clipping; user input 75 75 psS32 clipIter; ///< Number of clipping iterations; user input 76 doublemin; ///< minimum data value in array77 doublemax; ///< maximum data value in array78 doublebinsize; ///<76 psF64 min; ///< minimum data value in array 77 psF64 max; ///< maximum data value in array 78 psF64 binsize; ///< 79 79 psStatsOptions options; ///< bitmask of calculated values 80 80 } … … 132 132 */ 133 133 psHistogram* psHistogramAlloc( 134 floatlower, ///< Lower limit for the bins135 floatupper, ///< Upper limit for the bins134 psF32 lower, ///< Lower limit for the bins 135 psF32 upper, ///< Upper limit for the bins 136 136 psS32 n ///< Number of bins 137 137 ); … … 172 172 ///< the statistic struct to operate on 173 173 174 double*value174 psF64 *value 175 175 ///< if return is true, this is set to the specified statistic value by stats->options 176 176 ); -
trunk/psLib/src/math/psConstants.h
r2778 r2788 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1.4 8$ $Name: not supported by cvs2svn $9 * @date $Date: 2004-12-2 1 20:42:07$8 * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2004-12-22 05:09:32 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 12 * 13 * XXX: Add parenthesis around all arguments so that these macros can be 14 * called with complex expressions. 15 * 16 * 17 * 12 18 */ 13 19 -
trunk/psLib/src/math/psMinimize.c
r2741 r2788 9 9 * @author GLG, MHPCC 10 10 * 11 * @version $Revision: 1.9 5$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-12- 17 00:18:31$11 * @version $Revision: 1.96 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-12-22 05:09:32 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii 15 15 * 16 16 * XXX: must follow coding name standards on local functions. 17 *18 * XXX: Section 4.5.1.1 (predefined functions for Gauss minimization via19 * LMM) is not addressed here. We are waiting for subsequent SDRs20 * which will redefine the LMM functions.21 17 * 22 18 */ … … 64 60 XXX: Use a static vector. 65 61 *****************************************************************************/ 66 void psBuildSums1D( doublex,62 void psBuildSums1D(psF64 x, 67 63 psS32 polyOrder, 68 64 psVector* sums) 69 65 { 70 66 psS32 i = 0; 71 doublexSum = 0.0;67 psF64 xSum = 0.0; 72 68 73 69 if (sums == NULL) { … … 100 96 XXX: do an F64 version? 101 97 *****************************************************************************/ 102 float*CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices)98 psF32 *CalculateSecondDerivs(const psVector* restrict x, ///< Ordinates (or NULL to just use the indices) 103 99 const psVector* restrict y) ///< Coordinates 104 100 { … … 108 104 psS32 i; 109 105 psS32 k; 110 floatsig;111 floatp;106 psF32 sig; 107 psF32 p; 112 108 psS32 n = y->n; 113 float *u = (float *) psAlloc(n * sizeof(float));114 float *derivs2 = (float *) psAlloc(n * sizeof(float));115 float *X = (float*) & (x->data.F32[0]);116 float *Y = (float*) & (y->data.F32[0]);117 floatqn;109 psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); 110 psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); 111 psF32 *X = (psF32 *) & (x->data.F32[0]); 112 psF32 *Y = (psF32 *) & (y->data.F32[0]); 113 psF32 qn; 118 114 119 115 // XXX: The second derivatives at the endpoints, undefined in the SDR, … … 168 164 *****************************************************************************/ 169 165 /* 170 floatp_psNRSpline1DEval(psSpline1D *spline,166 psF32 p_psNRSpline1DEval(psSpline1D *spline, 171 167 const psVector* restrict x, 172 168 const psVector* restrict y, 173 floatX)169 psF32 X) 174 170 { 175 171 PS_PTR_CHECK_NULL(spline, NAN); … … 185 181 psS32 klo; 186 182 psS32 khi; 187 floatH;188 floatA;189 floatB;190 floatC;191 floatD;192 floatY;183 psF32 H; 184 psF32 A; 185 psF32 B; 186 psF32 C; 187 psF32 D; 188 psF32 Y; 193 189 194 190 n = spline->n; … … 258 254 "---- psVectorFitSpline1D() begin ----\n"); 259 255 psS32 numSplines = (y->n)-1; 260 floattmp;261 floatH;256 psF32 tmp; 257 psF32 H; 262 258 psS32 i; 263 floatslope;259 psF32 slope; 264 260 psVector *x32 = NULL; 265 261 psVector *y32 = NULL; … … 421 417 psTrace(".psLib.dataManip.psMinimize", 4, 422 418 "---- psMinimizeLMChi2Gauss1D() begin ----\n"); 423 floatx;424 inti;425 floatmean = params->data.F32[0];426 floatstdev = params->data.F32[1];419 psF32 x; 420 psS32 i; 421 psF32 mean = params->data.F32[0]; 422 psF32 stdev = params->data.F32[1]; 427 423 psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); 428 424 … … 443 439 for (i=0;i<coords->n;i++) { 444 440 x = ((psVector *) (coords->data[i]))->data.F32[0]; 445 floattmp = (x - mean) * psGaussian(x, mean, stdev, false);441 psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false); 446 442 deriv->data.F32[i][0] = tmp / (stdev * stdev); 447 443 tmp = (x - mean) * (x - mean) * … … 491 487 PS_PTR_CHECK_NULL(params, NULL); 492 488 493 doublenormalization = params->data.F32[0];494 doublex0 = params->data.F32[1];495 doubley0 = params->data.F32[2];496 doublesigmaX = params->data.F32[3];497 doublesigmaY = params->data.F32[4];498 doubletheta = params->data.F32[5];489 psF64 normalization = params->data.F32[0]; 490 psF64 x0 = params->data.F32[1]; 491 psF64 y0 = params->data.F32[2]; 492 psF64 sigmaX = params->data.F32[3]; 493 psF64 sigmaY = params->data.F32[4]; 494 psF64 theta = params->data.F32[5]; 499 495 psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); 500 496 … … 502 498 deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); 503 499 } else { 504 // XXX: Check size of derivative500 PS_IMAGE_CHECK_SIZE(deriv, 6, coords->n, NULL); 505 501 } 506 502 … … 508 504 "---- psMinimizeLMChi2Gauss2D() begin ----\n"); 509 505 510 for ( inti=0;i<coords->n;i++) {511 doublex = ((psVector *) coords->data[i])->data.F32[0];512 doubley = ((psVector *) coords->data[i])->data.F32[0];513 514 doubleu = - (x-x0)*cos(theta) + (y-y0)*sin(theta);515 doublev = (x-x0)*cos(theta) + (y-y0)*sin(theta);516 517 doubleflux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) +518 v*v/(2.0 * sigmaY * sigmaY)))/519 (2.0 * M_PI * sigmaX * sigmaY);506 for (psS32 i=0;i<coords->n;i++) { 507 psF64 x = ((psVector *) coords->data[i])->data.F32[0]; 508 psF64 y = ((psVector *) coords->data[i])->data.F32[0]; 509 510 psF64 u = - (x-x0)*cos(theta) + (y-y0)*sin(theta); 511 psF64 v = (x-x0)*cos(theta) + (y-y0)*sin(theta); 512 513 psF64 flux = normalization * exp(-( u*u/(2.0 * sigmaX * sigmaX) + 514 v*v/(2.0 * sigmaY * sigmaY)))/ 515 (2.0 * M_PI * sigmaX * sigmaY); 520 516 out->data.F32[i] = flux; 521 517 … … 584 580 PS_VECTOR_CHECK_SIZE_EQUAL(y, yErr, NULL); 585 581 } 586 587 // XXX: Generate code that modifies covar matrix if not-NULL.588 582 if (covar != NULL) { 589 583 PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL); … … 610 604 psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); 611 605 psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); 612 613 // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *));614 // for (i=0;i<numData;i++) {615 // deriv[i] = psVectorAlloc(numParams, PS_TYPE_F32);616 // }617 606 psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32); 618 619 psVector *currValueVec; 620 psVector *newValueVec; 621 622 float currChi2 = 0.0; 623 float newChi2 = 0.0; 624 float lamda = 0.00005; 607 psVector *currValueVec = NULL; 608 psVector *newValueVec = NULL; 609 psF32 currChi2 = 0.0; 610 psF32 newChi2 = 0.0; 611 psF32 lamda = 0.00005; 625 612 lamda = 0.05; 626 613 … … 711 698 if (j == k) { 712 699 A->data.F64[j][k] = 713 ( double) ((1.0 + lamda) * alpha->data.F32[j][k]);700 (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]); 714 701 } else { 715 A->data.F64[j][k] = ( double) alpha->data.F32[j][k];702 A->data.F64[j][k] = (psF64) alpha->data.F32[j][k]; 716 703 } 717 704 } … … 741 728 } else { 742 729 newParams->data.F32[i] = params->data.F32[i] - 743 ( float) paramDeltasF64->data.F64[i];730 (psF32) paramDeltasF64->data.F64[i]; 744 731 } 745 732 } … … 777 764 // We already masked params. 778 765 for (i=0;i<numParams;i++) { 779 params->data.F32[i] = ( float) newParams->data.F32[i];766 params->data.F32[i] = (psF32) newParams->data.F32[i]; 780 767 } 781 768 lamda*= 0.1; … … 826 813 psS32 k; 827 814 psS32 n = x->n; 828 doublefac;829 doublesum;815 psF64 fac; 816 psF64 sum; 830 817 PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64); 831 818 psScalar *fScalar; … … 836 823 // variable declarations. I retain them here to maintain coherence 837 824 // with the NR code. 838 doublemin = -1.0;839 doublemax = 1.0;840 doublebma = 0.5 * (max-min); // 1841 doublebpa = 0.5 * (max+min); // 0825 psF64 min = -1.0; 826 psF64 max = 1.0; 827 psF64 bma = 0.5 * (max-min); // 1 828 psF64 bpa = 0.5 * (max+min); // 0 842 829 843 830 // In this loop, we first calculate the values of X for which the … … 851 838 for (psS32 i=0;i<n;i++) { 852 839 // NR 5.8.4 853 double Y = cos(M_PI * (0.5 + ((float) i)) / ((float) n));854 doubleX = (Y + bma + bpa) - 1.0;840 psF64 Y = cos(M_PI * (0.5 + ((psF32) i)) / ((psF32) n)); 841 psF64 X = (Y + bma + bpa) - 1.0; 855 842 tmpScalar.data.F64 = X; 856 843 … … 873 860 // coefficients of the Chebyshev polynomial: NR 5.8.7. 874 861 875 fac = 2.0/(( float) n);862 fac = 2.0/((psF32) n); 876 863 // XXX: is this loop bound correct? 877 864 for (j=0;j<myPoly->n;j++) { … … 879 866 for (k=0;k<n;k++) { 880 867 sum+= f->data.F64[k] * 881 cos(M_PI * (( float) j) * (0.5 + ((float) k)) / ((float) n));868 cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n)); 882 869 } 883 870 … … 1075 1062 *****************************************************************************/ 1076 1063 psMinimization *psMinimizationAlloc(psS32 maxIter, 1077 floattol)1064 psF32 tol) 1078 1065 { 1079 1066 PS_INT_CHECK_NON_NEGATIVE(maxIter, NULL); … … 1092 1079 // LINE to it. We assume BASEMASK is non-null. 1093 1080 #define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ 1094 for ( inti=0;i<BASE->n;i++) { \1081 for (psS32 i=0;i<BASE->n;i++) { \ 1095 1082 if (BASEMASK->data.U8[i] == 0) { \ 1096 1083 OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ … … 1102 1089 #define PS_VECTOR_F32_CHECK_ZERO_VECTOR(IN, BOOL_VAR) \ 1103 1090 BOOL_VAR = true; \ 1104 for ( inti=0;i<IN->n;i++) { \1091 for (psS32 i=0;i<IN->n;i++) { \ 1105 1092 if (fabs(IN->data.F32[i]) >= FLT_EPSILON) { \ 1106 1093 BOOL_VAR = false; \ … … 1111 1098 #define PS_VECTOR_WITH_MASK_F32_CHECK_ZERO_VECTOR(IN, INMASK, BOOL_VAR) \ 1112 1099 BOOL_VAR = true; \ 1113 for ( inti=0;i<IN->n;i++) { \1100 for (psS32 i=0;i<IN->n;i++) { \ 1114 1101 if ((INMASK->data.U8[i] == 0) && (fabs(IN->data.F32[i]) >= FLT_EPSILON)) { \ 1115 1102 BOOL_VAR = false; \ … … 1152 1139 psMinimizePowellFunc func) 1153 1140 { 1154 floata = 0.0;1155 floatb = 0.0;1156 floatc = 0.0;1157 floatfa = 0.0;1158 floatfb = 0.0;1159 floatfc = 0.0;1141 psF32 a = 0.0; 1142 psF32 b = 0.0; 1143 psF32 c = 0.0; 1144 psF32 fa = 0.0; 1145 psF32 fb = 0.0; 1146 psF32 fc = 0.0; 1160 1147 psS32 iter = 100; 1161 floataDir = 0.0;1162 floatcDir = 0.0;1163 floatnew_aDir = 0.0;1164 floatnew_cDir = 0.0;1148 psF32 aDir = 0.0; 1149 psF32 cDir = 0.0; 1150 psF32 new_aDir = 0.0; 1151 psF32 new_cDir = 0.0; 1165 1152 psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); 1166 floatstepSize = PS_DETERMINE_BRACKET_STEP_SIZE;1153 psF32 stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; 1167 1154 psVector *tmp = NULL; 1168 1155 psBool boolLineIsNull = true; … … 1301 1288 psMinimizePowellFunc func) 1302 1289 { 1303 floata = 0.0;1304 floatb = 0.0;1305 floatc = 0.0;1306 floatfa = 0.0;1307 floatfb = 0.0;1308 floatfc = 0.0;1290 psF32 a = 0.0; 1291 psF32 b = 0.0; 1292 psF32 c = 0.0; 1293 psF32 fa = 0.0; 1294 psF32 fb = 0.0; 1295 psF32 fc = 0.0; 1309 1296 psS32 iter = 0; 1310 1297 PS_VECTOR_GEN_STATIC_RECYCLED(tmp, params->n, PS_TYPE_F32); 1311 1298 psBool boolLineIsNull = true; 1312 floatprevMin = 0.0;1313 intcountMin = 0;1299 psF32 prevMin = 0.0; 1300 psS32 countMin = 0; 1314 1301 1315 1302 psTrace(".psLib.dataManip.p_psDetermineBracket", 4, … … 1403 1390 *****************************************************************************/ 1404 1391 #define PS_LINEMIN_MAX_ITERATIONS 30 1405 floatp_psLineMin(psMinimization *min,1392 psF32 p_psLineMin(psMinimization *min, 1406 1393 psVector *params, 1407 1394 psVector *line, … … 1423 1410 PS_PTR_CHECK_NULL(func, NAN); 1424 1411 psVector *bracket; 1425 floata = 0.0;1426 floatb = 0.0;1427 floatc = 0.0;1428 floatn = 0.0;1429 floatfa = 0.0;1430 floatfb = 0.0;1431 floatfc = 0.0;1432 floatfn = 0.0;1433 floatmul = 0.0;1412 psF32 a = 0.0; 1413 psF32 b = 0.0; 1414 psF32 c = 0.0; 1415 psF32 n = 0.0; 1416 psF32 fa = 0.0; 1417 psF32 fb = 0.0; 1418 psF32 fc = 0.0; 1419 psF32 fn = 0.0; 1420 psF32 mul = 0.0; 1434 1421 PS_VECTOR_GEN_STATIC_RECYCLED(tmpa, params->n, PS_TYPE_F32); 1435 1422 PS_VECTOR_GEN_STATIC_RECYCLED(tmpb, params->n, PS_TYPE_F32); … … 1544 1531 This routine must minimize a possibly multi-dimensional function. The 1545 1532 function to be minimized "func" is: 1546 floatfunc(psVector *params, psArray *coords)1533 psF32 func(psVector *params, psArray *coords) 1547 1534 The "params" are the parameters of the function which are varied. The data 1548 1535 points at which the function is varied are in the argument "coords" which is … … 1578 1565 psVector *myParamMask = NULL; 1579 1566 psMinimization dummyMin; 1580 floatmul = 0.0;1581 floatbaseFuncVal = 0.0;1582 floatcurrFuncVal = 0.0;1567 psF32 mul = 0.0; 1568 psF32 baseFuncVal = 0.0; 1569 psF32 currFuncVal = 0.0; 1583 1570 psS32 biggestIter = 0; 1584 floatbiggestDiff = 0.0;1585 intiterationNumber = 0;1571 psF32 biggestDiff = 0.0; 1572 psS32 iterationNumber = 0; 1586 1573 1587 1574 psTrace(".psLib.dataManip.psMinimizePowell", 4, … … 1711 1698 } 1712 1699 } 1713 floatfqp = func(pQP, coords);1714 floatterm1 = (baseFuncVal - currFuncVal) - biggestDiff;1700 psF32 fqp = func(pQP, coords); 1701 psF32 term1 = (baseFuncVal - currFuncVal) - biggestDiff; 1715 1702 term1*= term1; 1716 1703 term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp); 1717 floatterm2 = baseFuncVal - fqp;1704 psF32 term2 = baseFuncVal - fqp; 1718 1705 term2*= term2 * biggestDiff; 1719 1706 if (term1 < term2) { … … 1760 1747 PS_PTR_CHECK_NULL(params, NULL); 1761 1748 1762 floatx;1763 inti;1764 floatmean = params->data.F32[0];1765 floatstdev = params->data.F32[1];1749 psF32 x; 1750 psS32 i; 1751 psF32 mean = params->data.F32[0]; 1752 psF32 stdev = params->data.F32[1]; 1766 1753 psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); 1767 1754 … … 1783 1770 XXX: This is F32 only 1784 1771 *****************************************************************************/ 1785 floatmyPowellChi2Func(const psVector *params,1772 psF32 myPowellChi2Func(const psVector *params, 1786 1773 const psArray *coords) 1787 1774 { … … 1794 1781 PS_PTR_CHECK_NULL(coords, NAN); 1795 1782 1796 floatchi2 = 0.0;1797 floatd;1783 psF32 chi2 = 0.0; 1784 psF32 d; 1798 1785 psS32 i; 1799 1786 psVector *tmp; -
trunk/psLib/src/math/psMinimize.h
r2250 r2788 8 8 * @author GLG, MHPCC 9 9 * 10 * @version $Revision: 1.3 6$ $Name: not supported by cvs2svn $11 * @date $Date: 2004-1 1-01 23:57:08$10 * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ 11 * @date $Date: 2004-12-22 05:09:32 $ 12 12 * 13 13 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 40 40 { 41 41 psS32 maxIter; ///< Convergence limit 42 floattol; ///< Error Tolerance43 floatvalue; ///< Value of function at minimum42 psF32 tol; ///< Error Tolerance 43 psF32 value; ///< Value of function at minimum 44 44 psS32 iter; ///< Number of iterations required 45 floatlastDelta; ///< The last difference for the fit45 psF32 lastDelta; ///< The last difference for the fit 46 46 } 47 47 psMinimization; 48 48 49 49 psMinimization *psMinimizationAlloc(psS32 maxIter, 50 floattol);50 psF32 tol); 51 51 52 52 /** Derive a polynomial fit. … … 91 91 92 92 typedef 93 float(*psMinimizePowellFunc)(const psVector *params,93 psF32 (*psMinimizePowellFunc)(const psVector *params, 94 94 const psArray *coords); 95 95 -
trunk/psLib/src/math/psPolynomial.c
r2778 r2788 7 7 * polynomials. It also contains a Gaussian functions. 8 8 * 9 * @version $Revision: 1.7 4$ $Name: not supported by cvs2svn $10 * @date $Date: 2004-12-2 1 20:42:07$9 * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-12-22 05:09:32 $ 11 11 * 12 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 55 55 static void dPolynomial4DFree(psDPolynomial4D* myPoly); 56 56 static void spline1DFree(psSpline1D *tmpSpline); 57 static psS32 vectorBinDisectF32( float *bins,psS32 numBins,floatx);57 static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); 58 58 static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); 59 59 … … 277 277 278 278 *****************************************************************************/ 279 static float ordPolynomial1DEval(floatx, const psPolynomial1D* myPoly)279 static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) 280 280 { 281 281 psS32 loop_x = 0; 282 floatpolySum = 0.0;283 floatxSum = 1.0;282 psF32 polySum = 0.0; 283 psF32 xSum = 1.0; 284 284 285 285 psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, … … 307 307 // XXX: How does the mask vector effect Crenshaw's formula? 308 308 // XXX: We assume that x is scaled between -1.0 and 1.0; 309 static float chebPolynomial1DEval(floatx, const psPolynomial1D* myPoly)309 static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) 310 310 { 311 311 // PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0); … … 314 314 psS32 n; 315 315 psS32 i; 316 floattmp;316 psF32 tmp; 317 317 318 318 n = myPoly->n; … … 337 337 psS32 n; 338 338 psS32 i; 339 floattmp;339 psF32 tmp; 340 340 psPolynomial1D **chebPolys = NULL; 341 341 … … 355 355 } 356 356 357 static float ordPolynomial2DEval(float x, floaty, const psPolynomial2D* myPoly)357 static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) 358 358 { 359 359 PS_POLY_CHECK_NULL(myPoly, NAN); … … 361 361 psS32 loop_x = 0; 362 362 psS32 loop_y = 0; 363 floatpolySum = 0.0;364 floatxSum = 1.0;365 floatySum = 1.0;363 psF32 polySum = 0.0; 364 psF32 xSum = 1.0; 365 psF32 ySum = 1.0; 366 366 367 367 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 379 379 } 380 380 381 static float chebPolynomial2DEval(float x, floaty, const psPolynomial2D* myPoly)381 static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) 382 382 { 383 383 PS_POLY_CHECK_NULL(myPoly, NAN); … … 386 386 psS32 loop_y = 0; 387 387 psS32 i = 0; 388 floatpolySum = 0.0;388 psF32 polySum = 0.0; 389 389 psPolynomial1D* *chebPolys = NULL; 390 390 psS32 maxChebyPoly = 0; … … 414 414 } 415 415 416 static float ordPolynomial3DEval(float x, float y, floatz, const psPolynomial3D* myPoly)416 static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) 417 417 { 418 418 psS32 loop_x = 0; 419 419 psS32 loop_y = 0; 420 420 psS32 loop_z = 0; 421 floatpolySum = 0.0;422 floatxSum = 1.0;423 floatySum = 1.0;424 floatzSum = 1.0;421 psF32 polySum = 0.0; 422 psF32 xSum = 1.0; 423 psF32 ySum = 1.0; 424 psF32 zSum = 1.0; 425 425 426 426 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 442 442 } 443 443 444 static float chebPolynomial3DEval(float x, float y, floatz, const psPolynomial3D* myPoly)444 static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) 445 445 { 446 446 psS32 loop_x = 0; … … 448 448 psS32 loop_z = 0; 449 449 psS32 i = 0; 450 floatpolySum = 0.0;450 psF32 polySum = 0.0; 451 451 psPolynomial1D* *chebPolys = NULL; 452 452 psS32 maxChebyPoly = 0; … … 483 483 } 484 484 485 static float ordPolynomial4DEval(float w, float x, float y, floatz, const psPolynomial4D* myPoly)485 static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) 486 486 { 487 487 psS32 loop_w = 0; … … 489 489 psS32 loop_y = 0; 490 490 psS32 loop_z = 0; 491 floatpolySum = 0.0;492 floatwSum = 1.0;493 floatxSum = 1.0;494 floatySum = 1.0;495 floatzSum = 1.0;491 psF32 polySum = 0.0; 492 psF32 wSum = 1.0; 493 psF32 xSum = 1.0; 494 psF32 ySum = 1.0; 495 psF32 zSum = 1.0; 496 496 497 497 for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { … … 517 517 } 518 518 519 static float chebPolynomial4DEval(float w, float x, float y, floatz, const psPolynomial4D* myPoly)519 static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) 520 520 { 521 521 psS32 loop_w = 0; … … 524 524 psS32 loop_z = 0; 525 525 psS32 i = 0; 526 floatpolySum = 0.0;526 psF32 polySum = 0.0; 527 527 psPolynomial1D* *chebPolys = NULL; 528 528 psS32 maxChebyPoly = 0; … … 568 568 Polynomial coefficients will be accessed in [w][x][y][z] fashion. 569 569 *****************************************************************************/ 570 static double dOrdPolynomial1DEval(doublex, const psDPolynomial1D* myPoly)570 static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) 571 571 { 572 572 psS32 loop_x = 0; 573 doublepolySum = 0.0;574 doublexSum = 1.0;573 psF64 polySum = 0.0; 574 psF64 xSum = 1.0; 575 575 576 576 for (loop_x = 0; loop_x < myPoly->n; loop_x++) { … … 586 586 // XXX: You can do this without having to psAlloc() vector d. 587 587 // XXX: How does the mask vector effect Crenshaw's formula? 588 static double dChebPolynomial1DEval(doublex, const psDPolynomial1D* myPoly)588 static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) 589 589 { 590 590 psVector *d; 591 591 psS32 n; 592 592 psS32 i; 593 doubletmp;593 psF64 tmp; 594 594 595 595 n = myPoly->n; … … 611 611 } 612 612 613 static double dOrdPolynomial2DEval(double x, doubley, const psDPolynomial2D* myPoly)613 static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) 614 614 { 615 615 psS32 loop_x = 0; 616 616 psS32 loop_y = 0; 617 doublepolySum = 0.0;618 doublexSum = 1.0;619 doubleySum = 1.0;617 psF64 polySum = 0.0; 618 psF64 xSum = 1.0; 619 psF64 ySum = 1.0; 620 620 621 621 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 633 633 } 634 634 635 static double dChebPolynomial2DEval(double x, doubley, const psDPolynomial2D* myPoly)635 static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) 636 636 { 637 637 psS32 loop_x = 0; 638 638 psS32 loop_y = 0; 639 639 psS32 i = 0; 640 doublepolySum = 0.0;640 psF64 polySum = 0.0; 641 641 psPolynomial1D* *chebPolys = NULL; 642 642 psS32 maxChebyPoly = 0; … … 667 667 } 668 668 669 static double dOrdPolynomial3DEval(double x, double y, doublez, const psDPolynomial3D* myPoly)669 static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) 670 670 { 671 671 psS32 loop_x = 0; 672 672 psS32 loop_y = 0; 673 673 psS32 loop_z = 0; 674 doublepolySum = 0.0;675 doublexSum = 1.0;676 doubleySum = 1.0;677 doublezSum = 1.0;674 psF64 polySum = 0.0; 675 psF64 xSum = 1.0; 676 psF64 ySum = 1.0; 677 psF64 zSum = 1.0; 678 678 679 679 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 695 695 } 696 696 697 static double dChebPolynomial3DEval(double x, double y, doublez, const psDPolynomial3D* myPoly)697 static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) 698 698 { 699 699 psS32 loop_x = 0; … … 701 701 psS32 loop_z = 0; 702 702 psS32 i = 0; 703 doublepolySum = 0.0;703 psF64 polySum = 0.0; 704 704 psPolynomial1D* *chebPolys = NULL; 705 705 psS32 maxChebyPoly = 0; … … 736 736 } 737 737 738 static double dOrdPolynomial4DEval(double w, double x, double y, doublez, const psDPolynomial4D* myPoly)738 static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) 739 739 { 740 740 psS32 loop_w = 0; … … 742 742 psS32 loop_y = 0; 743 743 psS32 loop_z = 0; 744 doublepolySum = 0.0;745 doublewSum = 1.0;746 doublexSum = 1.0;747 doubleySum = 1.0;748 doublezSum = 1.0;744 psF64 polySum = 0.0; 745 psF64 wSum = 1.0; 746 psF64 xSum = 1.0; 747 psF64 ySum = 1.0; 748 psF64 zSum = 1.0; 749 749 750 750 for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { … … 770 770 } 771 771 772 static double dChebPolynomial4DEval(double w, double x, double y, doublez, const psDPolynomial4D* myPoly)772 static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) 773 773 { 774 774 psS32 loop_w = 0; … … 777 777 psS32 loop_z = 0; 778 778 psS32 i = 0; 779 doublepolySum = 0.0;779 psF64 polySum = 0.0; 780 780 psPolynomial1D* *chebPolys = NULL; 781 781 psS32 maxChebyPoly = 0; … … 828 828 *****************************************************************************/ 829 829 #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ 830 static floatfullInterpolate1D##TYPE(ps##TYPE *domain, \830 static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ 831 831 ps##TYPE *range, \ 832 832 psS32 n, \ … … 896 896 LaGrange interpolation. 897 897 *****************************************************************************/ 898 static float interpolate1DF32(float*domain,899 float*range,898 static psF32 interpolate1DF32(psF32 *domain, 899 psF32 *range, 900 900 psS32 n, 901 901 psS32 order, 902 floatx)902 psF32 x) 903 903 { 904 904 psS32 binNum; … … 941 941 evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] 942 942 *****************************************************************************/ 943 float psGaussian(float x, float mean, floatsigma, psBool normal)944 { 945 floattmp = 1.0;943 psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal) 944 { 945 psF32 tmp = 1.0; 946 946 947 947 psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, … … 967 967 XXX: There is no way to seed the random generator. 968 968 *****************************************************************************/ 969 psVector* p_psGaussianDev( float mean, floatsigma, psS32 Npts)969 psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) 970 970 { 971 971 PS_INT_CHECK_NON_NEGATIVE(Npts, NULL); … … 1007 1007 newPoly->type = type; 1008 1008 newPoly->n = n; 1009 newPoly->coeff = ( float *)psAlloc(n * sizeof(float));1010 newPoly->coeffErr = ( float *)psAlloc(n * sizeof(float));1011 newPoly->mask = ( char *)psAlloc(n * sizeof(char));1009 newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32)); 1010 newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32)); 1011 newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); 1012 1012 for (i = 0; i < n; i++) { 1013 1013 newPoly->coeff[i] = 0.0; … … 1036 1036 newPoly->nY = nY; 1037 1037 1038 newPoly->coeff = ( float **)psAlloc(nX * sizeof(float*));1039 newPoly->coeffErr = ( float **)psAlloc(nX * sizeof(float*));1040 newPoly->mask = ( char **)psAlloc(nX * sizeof(char*));1038 newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *)); 1039 newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *)); 1040 newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); 1041 1041 for (x = 0; x < nX; x++) { 1042 newPoly->coeff[x] = ( float *)psAlloc(nY * sizeof(float));1043 newPoly->coeffErr[x] = ( float *)psAlloc(nY * sizeof(float));1044 newPoly->mask[x] = ( char *)psAlloc(nY * sizeof(char));1042 newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); 1043 newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); 1044 newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); 1045 1045 } 1046 1046 for (x = 0; x < nX; x++) { … … 1075 1075 newPoly->nZ = nZ; 1076 1076 1077 newPoly->coeff = ( float ***)psAlloc(nX * sizeof(float**));1078 newPoly->coeffErr = ( float ***)psAlloc(nX * sizeof(float**));1079 newPoly->mask = ( char ***)psAlloc(nX * sizeof(char**));1077 newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1078 newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1079 newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1080 1080 for (x = 0; x < nX; x++) { 1081 newPoly->coeff[x] = ( float **)psAlloc(nY * sizeof(float*));1082 newPoly->coeffErr[x] = ( float **)psAlloc(nY * sizeof(float*));1083 newPoly->mask[x] = ( char **)psAlloc(nY * sizeof(char*));1081 newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1082 newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1083 newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1084 1084 for (y = 0; y < nY; y++) { 1085 newPoly->coeff[x][y] = ( float *)psAlloc(nZ * sizeof(float));1086 newPoly->coeffErr[x][y] = ( float *)psAlloc(nZ * sizeof(float));1087 newPoly->mask[x][y] = ( char *)psAlloc(nZ * sizeof(char));1085 newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1086 newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1087 newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1088 1088 } 1089 1089 } … … 1124 1124 newPoly->nZ = nZ; 1125 1125 1126 newPoly->coeff = ( float ****)psAlloc(nW * sizeof(float***));1127 newPoly->coeffErr = ( float ****)psAlloc(nW * sizeof(float***));1128 newPoly->mask = ( char ****)psAlloc(nW * sizeof(char***));1126 newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); 1127 newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); 1128 newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); 1129 1129 for (w = 0; w < nW; w++) { 1130 newPoly->coeff[w] = ( float ***)psAlloc(nX * sizeof(float**));1131 newPoly->coeffErr[w] = ( float ***)psAlloc(nX * sizeof(float**));1132 newPoly->mask[w] = ( char ***)psAlloc(nX * sizeof(char**));1130 newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1131 newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1132 newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1133 1133 for (x = 0; x < nX; x++) { 1134 newPoly->coeff[w][x] = ( float **)psAlloc(nY * sizeof(float*));1135 newPoly->coeffErr[w][x] = ( float **)psAlloc(nY * sizeof(float*));1136 newPoly->mask[w][x] = ( char **)psAlloc(nY * sizeof(char*));1134 newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1135 newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1136 newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1137 1137 for (y = 0; y < nY; y++) { 1138 newPoly->coeff[w][x][y] = ( float *)psAlloc(nZ * sizeof(float));1139 newPoly->coeffErr[w][x][y] = ( float *)psAlloc(nZ * sizeof(float));1140 newPoly->mask[w][x][y] = ( char *)psAlloc(nZ * sizeof(char));1138 newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1139 newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1140 newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1141 1141 } 1142 1142 } … … 1157 1157 } 1158 1158 1159 float psPolynomial1DEval(const psPolynomial1D* myPoly, floatx)1159 psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x) 1160 1160 { 1161 1161 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1195 1195 } 1196 1196 1197 float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, floaty)1197 psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y) 1198 1198 { 1199 1199 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1250 1250 } 1251 1251 1252 float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, floatz)1252 psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z) 1253 1253 { 1254 1254 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1316 1316 } 1317 1317 1318 float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, floatz)1318 psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z) 1319 1319 { 1320 1320 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1406 1406 newPoly->type = type; 1407 1407 newPoly->n = n; 1408 newPoly->coeff = ( double *)psAlloc(n * sizeof(double));1409 newPoly->coeffErr = ( double *)psAlloc(n * sizeof(double));1410 newPoly->mask = ( char *)psAlloc(n * sizeof(char));1408 newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64)); 1409 newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64)); 1410 newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); 1411 1411 for (i = 0; i < n; i++) { 1412 1412 newPoly->coeff[i] = 0.0; … … 1435 1435 newPoly->nY = nY; 1436 1436 1437 newPoly->coeff = ( double **)psAlloc(nX * sizeof(double*));1438 newPoly->coeffErr = ( double **)psAlloc(nX * sizeof(double*));1439 newPoly->mask = ( char **)psAlloc(nX * sizeof(char*));1437 newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *)); 1438 newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *)); 1439 newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); 1440 1440 for (x = 0; x < nX; x++) { 1441 newPoly->coeff[x] = ( double *)psAlloc(nY * sizeof(double));1442 newPoly->coeffErr[x] = ( double *)psAlloc(nY * sizeof(double));1443 newPoly->mask[x] = ( char *)psAlloc(nY * sizeof(char));1441 newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); 1442 newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); 1443 newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); 1444 1444 } 1445 1445 for (x = 0; x < nX; x++) { … … 1474 1474 newPoly->nZ = nZ; 1475 1475 1476 newPoly->coeff = ( double ***)psAlloc(nX * sizeof(double**));1477 newPoly->coeffErr = ( double ***)psAlloc(nX * sizeof(double**));1478 newPoly->mask = ( char ***)psAlloc(nX * sizeof(char**));1476 newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1477 newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1478 newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1479 1479 for (x = 0; x < nX; x++) { 1480 newPoly->coeff[x] = ( double **)psAlloc(nY * sizeof(double*));1481 newPoly->coeffErr[x] = ( double **)psAlloc(nY * sizeof(double*));1482 newPoly->mask[x] = ( char **)psAlloc(nY * sizeof(char*));1480 newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1481 newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1482 newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1483 1483 for (y = 0; y < nY; y++) { 1484 newPoly->coeff[x][y] = ( double *)psAlloc(nZ * sizeof(double));1485 newPoly->coeffErr[x][y] = ( double *)psAlloc(nZ * sizeof(double));1486 newPoly->mask[x][y] = ( char *)psAlloc(nZ * sizeof(char));1484 newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1485 newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1486 newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1487 1487 } 1488 1488 } … … 1523 1523 newPoly->nZ = nZ; 1524 1524 1525 newPoly->coeff = ( double ****)psAlloc(nW * sizeof(double***));1526 newPoly->coeffErr = ( double ****)psAlloc(nW * sizeof(double***));1527 newPoly->mask = ( char ****)psAlloc(nW * sizeof(char***));1525 newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); 1526 newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); 1527 newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); 1528 1528 for (w = 0; w < nW; w++) { 1529 newPoly->coeff[w] = ( double ***)psAlloc(nX * sizeof(double**));1530 newPoly->coeffErr[w] = ( double ***)psAlloc(nX * sizeof(double**));1531 newPoly->mask[w] = ( char ***)psAlloc(nX * sizeof(char**));1529 newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1530 newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1531 newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1532 1532 for (x = 0; x < nX; x++) { 1533 newPoly->coeff[w][x] = ( double **)psAlloc(nY * sizeof(double*));1534 newPoly->coeffErr[w][x] = ( double **)psAlloc(nY * sizeof(double*));1535 newPoly->mask[w][x] = ( char **)psAlloc(nY * sizeof(char*));1533 newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1534 newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1535 newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1536 1536 for (y = 0; y < nY; y++) { 1537 newPoly->coeff[w][x][y] = ( double *)psAlloc(nZ * sizeof(double));1538 newPoly->coeffErr[w][x][y] = ( double *)psAlloc(nZ * sizeof(double));1539 newPoly->mask[w][x][y] = ( char *)psAlloc(nZ * sizeof(char));1537 newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1538 newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1539 newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1540 1540 } 1541 1541 } … … 1557 1557 1558 1558 1559 double psDPolynomial1DEval(const psDPolynomial1D* myPoly, doublex)1559 psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x) 1560 1560 { 1561 1561 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1598 1598 1599 1599 1600 doublepsDPolynomial2DEval(const psDPolynomial2D* myPoly,1601 doublex,1602 doubley)1600 psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly, 1601 psF64 x, 1602 psF64 y) 1603 1603 { 1604 1604 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1655 1655 1656 1656 1657 doublepsDPolynomial3DEval(const psDPolynomial3D* myPoly,1658 doublex,1659 doubley,1660 doublez)1657 psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly, 1658 psF64 x, 1659 psF64 y, 1660 psF64 z) 1661 1661 { 1662 1662 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1724 1724 } 1725 1725 1726 doublepsDPolynomial4DEval(const psDPolynomial4D* myPoly,1727 doublew,1728 doublex,1729 doubley,1730 doublez)1726 psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly, 1727 psF64 w, 1728 psF64 x, 1729 psF64 y, 1730 psF64 z) 1731 1731 { 1732 1732 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1810 1810 // psS32 n; 1811 1811 // psPolynomial1D **spline; 1812 // float*p_psDeriv2;1813 // float*domains;1812 // psF32 *p_psDeriv2; 1813 // psF32 *domains; 1814 1814 //} psSpline1D; 1815 1815 … … 1822 1822 psSpline1D *psSpline1DAlloc(psS32 numSplines, 1823 1823 psS32 order, 1824 floatmin,1825 floatmax)1824 psF32 min, 1825 psF32 max) 1826 1826 { 1827 1827 PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL); … … 1831 1831 psSpline1D *tmp = NULL; 1832 1832 psS32 i; 1833 floattmpDomain;1834 floatwidth;1833 psF32 tmpDomain; 1834 psF32 width; 1835 1835 1836 1836 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); … … 1845 1845 tmp->p_psDeriv2 = NULL; 1846 1846 1847 tmp->domains = ( float *) psAlloc((numSplines+1) * sizeof(float));1848 width = (max - min) / (( float) numSplines);1847 tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32)); 1848 width = (max - min) / ((psF32) numSplines); 1849 1849 1850 1850 (tmp->domains)[0] = min; … … 1884 1884 } 1885 1885 1886 tmp->domains = ( float *) psAlloc((bounds->n) * sizeof(float));1886 tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32)); 1887 1887 1888 1888 for (i=0;i<bounds->n;i++) { … … 1977 1977 { 1978 1978 PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3); 1979 char* strType;1979 psS8* strType; 1980 1980 1981 1981 switch (x->type.type) { … … 2074 2074 domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); 2075 2075 2076 psScalar *tmpScalar = psScalarAlloc(( double)2076 psScalar *tmpScalar = psScalarAlloc((psF64) 2077 2077 interpolate1DF32(domain32->data.F32, 2078 2078 range32->data.F32, 2079 2079 domain32->n, 2080 2080 order, 2081 ( float) x->data.F64), PS_TYPE_F64);2081 (psF32) x->data.F64), PS_TYPE_F64); 2082 2082 psFree(range32); 2083 2083 psFree(domain32); … … 2089 2089 2090 2090 } else { 2091 char* strType;2091 psS8* strType; 2092 2092 PS_TYPE_NAME(strType,x->type.type); 2093 2093 psError(PS_ERR_BAD_PARAMETER_TYPE, … … 2115 2115 the spline fit functions require F32 and F64. 2116 2116 *****************************************************************************/ 2117 floatpsSpline1DEval(2118 floatx,2117 psF32 psSpline1DEval( 2118 psF32 x, 2119 2119 const psSpline1D *spline 2120 2120 ) … … 2172 2172 for (i=0;i<x->n;i++) { 2173 2173 tmpVector->data.F32[i] = psSpline1DEval( 2174 ( float) x->data.F64[i],2174 (psF32) x->data.F64[i], 2175 2175 spline 2176 2176 ); 2177 2177 } 2178 2178 } else { 2179 char* strType;2179 psS8* strType; 2180 2180 PS_TYPE_NAME(strType,x->type.type); 2181 2181 psError(PS_ERR_BAD_PARAMETER_TYPE, -
trunk/psLib/src/math/psPolynomial.h
r2600 r2788 12 12 * @author GLG, MHPCC 13 13 * 14 * @version $Revision: 1.3 7$ $Name: not supported by cvs2svn $15 * @date $Date: 2004-12- 02 21:12:51$14 * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ 15 * @date $Date: 2004-12-22 05:09:32 $ 16 16 * 17 17 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 39 39 * \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] 40 40 * 41 * @return floatvalue on the gaussian curve given the input parameters42 */ 43 floatpsGaussian(44 floatx, ///< Value at which to evaluate45 floatmean, ///< Mean for the Gaussian46 floatstddev, ///< Standard deviation for the Gaussian41 * @return psF32 value on the gaussian curve given the input parameters 42 */ 43 psF32 psGaussian( 44 psF32 x, ///< Value at which to evaluate 45 psF32 mean, ///< Mean for the Gaussian 46 psF32 stddev, ///< Standard deviation for the Gaussian 47 47 psBool normal ///< Indicates whether result should be normalized 48 48 ); … … 55 55 */ 56 56 psVector* p_psGaussianDev( 57 floatmean, ///< The mean of the Gaussian58 floatsigma, ///< The sigma of the Gaussian57 psF32 mean, ///< The mean of the Gaussian 58 psF32 sigma, ///< The sigma of the Gaussian 59 59 psS32 Npts ///< The size of the vector 60 60 ); … … 70 70 psPolynomialType type; ///< Polynomial type 71 71 psS32 n; ///< Number of terms 72 float*coeff; ///< Coefficients73 float*coeffErr; ///< Error in coefficients74 char*mask; ///< Coefficient mask72 psF32 *coeff; ///< Coefficients 73 psF32 *coeffErr; ///< Error in coefficients 74 psU8 *mask; ///< Coefficient mask 75 75 } 76 76 psPolynomial1D; … … 82 82 psS32 nX; ///< Number of terms in x 83 83 psS32 nY; ///< Number of terms in y 84 float**coeff; ///< Coefficients85 float**coeffErr; ///< Error in coefficients86 char**mask; ///< Coefficients mask84 psF32 **coeff; ///< Coefficients 85 psF32 **coeffErr; ///< Error in coefficients 86 psU8 **mask; ///< Coefficients mask 87 87 } 88 88 psPolynomial2D; … … 95 95 psS32 nY; ///< Number of terms in y 96 96 psS32 nZ; ///< Number of terms in z 97 float***coeff; ///< Coefficients98 float***coeffErr; ///< Error in coefficients99 char***mask; ///< Coefficients mask97 psF32 ***coeff; ///< Coefficients 98 psF32 ***coeffErr; ///< Error in coefficients 99 psU8 ***mask; ///< Coefficients mask 100 100 } 101 101 psPolynomial3D; … … 109 109 psS32 nY; ///< Number of terms in y 110 110 psS32 nZ; ///< Number of terms in z 111 float****coeff; ///< Coefficients112 float****coeffErr; ///< Error in coefficients113 char****mask; ///< Coefficients mask111 psF32 ****coeff; ///< Coefficients 112 psF32 ****coeffErr; ///< Error in coefficients 113 psU8 ****mask; ///< Coefficients mask 114 114 } 115 115 psPolynomial4D; … … 160 160 /** Evaluates a 1-D polynomial at specific coordinates. 161 161 * 162 * @return floatresult of polynomial at given location163 */ 164 floatpsPolynomial1DEval(162 * @return psF32 result of polynomial at given location 163 */ 164 psF32 psPolynomial1DEval( 165 165 const psPolynomial1D* myPoly, ///< Coefficients for the polynomial 166 floatx ///< location at which to evaluate166 psF32 x ///< location at which to evaluate 167 167 ); 168 168 169 169 /** Evaluates a 2-D polynomial at specific coordinates. 170 170 * 171 * @return floatresult of polynomial at given location172 */ 173 floatpsPolynomial2DEval(171 * @return psF32 result of polynomial at given location 172 */ 173 psF32 psPolynomial2DEval( 174 174 const psPolynomial2D* myPoly, ///< Coefficients for the polynomial 175 floatx, ///< x location at which to evaluate176 floaty ///< y location at which to evaluate175 psF32 x, ///< x location at which to evaluate 176 psF32 y ///< y location at which to evaluate 177 177 ); 178 178 179 179 /** Evaluates a 3-D polynomial at specific coordinates. 180 180 * 181 * @return floatresult of polynomial at given location182 */ 183 floatpsPolynomial3DEval(181 * @return psF32 result of polynomial at given location 182 */ 183 psF32 psPolynomial3DEval( 184 184 const psPolynomial3D* myPoly, ///< Coefficients for the polynomial 185 floatx, ///< x location at which to evaluate186 floaty, ///< y location at which to evaluate187 floatz ///< z location at which to evaluate185 psF32 x, ///< x location at which to evaluate 186 psF32 y, ///< y location at which to evaluate 187 psF32 z ///< z location at which to evaluate 188 188 ); 189 189 190 190 /** Evaluates a 4-D polynomial at specific coordinates. 191 191 * 192 * @return floatresult of polynomial at given location193 */ 194 floatpsPolynomial4DEval(192 * @return psF32 result of polynomial at given location 193 */ 194 psF32 psPolynomial4DEval( 195 195 const psPolynomial4D* myPoly, ///< Coefficients for the polynomial 196 floatw, ///< w location at which to evaluate197 floatx, ///< x location at which to evaluate198 floaty, ///< y location at which to evaluate199 floatz ///< z location at which to evaluate196 psF32 w, ///< w location at which to evaluate 197 psF32 x, ///< x location at which to evaluate 198 psF32 y, ///< y location at which to evaluate 199 psF32 z ///< z location at which to evaluate 200 200 ); 201 201 … … 235 235 psPolynomialType type; ///< Polynomial type 236 236 psS32 n; ///< Number of terms 237 double*coeff; ///< Coefficients238 double*coeffErr; ///< Error in coefficients239 char*mask; ///< Coefficient mask237 psF64 *coeff; ///< Coefficients 238 psF64 *coeffErr; ///< Error in coefficients 239 psU8 *mask; ///< Coefficient mask 240 240 } 241 241 psDPolynomial1D; … … 247 247 psS32 nX; ///< Number of terms in x 248 248 psS32 nY; ///< Number of terms in y 249 double**coeff; ///< Coefficients250 double**coeffErr; ///< Error in coefficients251 char**mask; ///< Coefficients mask249 psF64 **coeff; ///< Coefficients 250 psF64 **coeffErr; ///< Error in coefficients 251 psU8 **mask; ///< Coefficients mask 252 252 } 253 253 psDPolynomial2D; … … 260 260 psS32 nY; ///< Number of terms in y 261 261 psS32 nZ; ///< Number of terms in z 262 double***coeff; ///< Coefficients263 double***coeffErr; ///< Error in coefficients264 char***mask; ///< Coefficient mask262 psF64 ***coeff; ///< Coefficients 263 psF64 ***coeffErr; ///< Error in coefficients 264 psU8 ***mask; ///< Coefficient mask 265 265 } 266 266 psDPolynomial3D; … … 274 274 psS32 nY; ///< Number of terms in y 275 275 psS32 nZ; ///< Number of terms in z 276 double****coeff; ///< Coefficients277 double****coeffErr; ///< Error in coefficients278 char****mask; ///< Coefficients mask276 psF64 ****coeff; ///< Coefficients 277 psF64 ****coeffErr; ///< Error in coefficients 278 psU8 ****mask; ///< Coefficients mask 279 279 } 280 280 psDPolynomial4D; … … 324 324 /** Evaluates a double-precision 1-D polynomial at specific coordinates. 325 325 * 326 * @return floatresult of polynomial at given location327 */ 328 doublepsDPolynomial1DEval(326 * @return psF32 result of polynomial at given location 327 */ 328 psF64 psDPolynomial1DEval( 329 329 const psDPolynomial1D* myPoly, ///< Coefficients for the polynomial 330 doublex ///< Value at which to evaluate330 psF64 x ///< Value at which to evaluate 331 331 ); 332 332 333 333 /** Evaluates a double-precision 2-D polynomial at specific coordinates. 334 334 * 335 * @return floatresult of polynomial at given location336 */ 337 doublepsDPolynomial2DEval(335 * @return psF32 result of polynomial at given location 336 */ 337 psF64 psDPolynomial2DEval( 338 338 const psDPolynomial2D* myPoly, ///< Coefficients for the polynomial 339 doublex, ///< Value x at which to evaluate340 doubley ///< Value y at which to evaluate339 psF64 x, ///< Value x at which to evaluate 340 psF64 y ///< Value y at which to evaluate 341 341 ); 342 342 343 343 /** Evaluates a double-precision 3-D polynomial at specific coordinates. 344 344 * 345 * @return floatresult of polynomial at given location346 */ 347 doublepsDPolynomial3DEval(345 * @return psF32 result of polynomial at given location 346 */ 347 psF64 psDPolynomial3DEval( 348 348 const psDPolynomial3D* myPoly, ///< Coefficients for the polynomial 349 doublex, ///< Value x at which to evaluate350 doubley, ///< Value y at which to evaluate351 doublez ///< Value z at which to evaluate349 psF64 x, ///< Value x at which to evaluate 350 psF64 y, ///< Value y at which to evaluate 351 psF64 z ///< Value z at which to evaluate 352 352 ); 353 353 354 354 /** Evaluates a double-precision 4-D polynomial at specific coordinates. 355 355 * 356 * @return floatresult of polynomial at given location357 */ 358 doublepsDPolynomial4DEval(356 * @return psF32 result of polynomial at given location 357 */ 358 psF64 psDPolynomial4DEval( 359 359 const psDPolynomial4D* myPoly, ///< Coefficients for the polynomial 360 doublew, ///< Value w at which to evaluate361 doublex, ///< Value x at which to evaluate362 doubley, ///< Value y at which to evaluate363 doublez ///< Value z at which to evaluate360 psF64 w, ///< Value w at which to evaluate 361 psF64 x, ///< Value x at which to evaluate 362 psF64 y, ///< Value y at which to evaluate 363 psF64 z ///< Value z at which to evaluate 364 364 ); 365 365 … … 396 396 psS32 n; ///< The number of spline polynomials 397 397 psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials 398 float*p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1.399 float*domains; ///< The boundaries between each spline piece. Size is n+1.398 psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. 399 psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. 400 400 } 401 401 psSpline1D; … … 403 403 psSpline1D *psSpline1DAlloc(psS32 n, 404 404 psS32 order, 405 floatmin,406 floatmax);405 psF32 min, 406 psF32 max); 407 407 408 408 psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, 409 409 psS32 order); 410 410 411 floatpsSpline1DEval(412 floatx,411 psF32 psSpline1DEval( 412 psF32 x, 413 413 const psSpline1D *spline 414 414 ); … … 427 427 psScalar *x); 428 428 429 floatp_psNRSpline1DEval(psSpline1D *spline,429 psF32 p_psNRSpline1DEval(psSpline1D *spline, 430 430 const psVector* restrict x, 431 431 const psVector* restrict y, 432 floatX);432 psF32 X); 433 433 434 434 /* \} */// End of MathGroup Functions -
trunk/psLib/src/math/psSpline.c
r2778 r2788 7 7 * polynomials. It also contains a Gaussian functions. 8 8 * 9 * @version $Revision: 1.7 4$ $Name: not supported by cvs2svn $10 * @date $Date: 2004-12-2 1 20:42:07$9 * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-12-22 05:09:32 $ 11 11 * 12 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 55 55 static void dPolynomial4DFree(psDPolynomial4D* myPoly); 56 56 static void spline1DFree(psSpline1D *tmpSpline); 57 static psS32 vectorBinDisectF32( float *bins,psS32 numBins,floatx);57 static psS32 vectorBinDisectF32(psF32 *bins,psS32 numBins,psF32 x); 58 58 static psS32 vectorBinDisectS32(psS32 *bins,psS32 numBins,psS32 x); 59 59 … … 277 277 278 278 *****************************************************************************/ 279 static float ordPolynomial1DEval(floatx, const psPolynomial1D* myPoly)279 static psF32 ordPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) 280 280 { 281 281 psS32 loop_x = 0; 282 floatpolySum = 0.0;283 floatxSum = 1.0;282 psF32 polySum = 0.0; 283 psF32 xSum = 1.0; 284 284 285 285 psTrace(".psLib.dataManip.psFunctions.ordPolynomial1DEval", 4, … … 307 307 // XXX: How does the mask vector effect Crenshaw's formula? 308 308 // XXX: We assume that x is scaled between -1.0 and 1.0; 309 static float chebPolynomial1DEval(floatx, const psPolynomial1D* myPoly)309 static psF32 chebPolynomial1DEval(psF32 x, const psPolynomial1D* myPoly) 310 310 { 311 311 // PS_FLOAT_CHECK_RANGE(x, -1.0, 1.0, 0.0); … … 314 314 psS32 n; 315 315 psS32 i; 316 floattmp;316 psF32 tmp; 317 317 318 318 n = myPoly->n; … … 337 337 psS32 n; 338 338 psS32 i; 339 floattmp;339 psF32 tmp; 340 340 psPolynomial1D **chebPolys = NULL; 341 341 … … 355 355 } 356 356 357 static float ordPolynomial2DEval(float x, floaty, const psPolynomial2D* myPoly)357 static psF32 ordPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) 358 358 { 359 359 PS_POLY_CHECK_NULL(myPoly, NAN); … … 361 361 psS32 loop_x = 0; 362 362 psS32 loop_y = 0; 363 floatpolySum = 0.0;364 floatxSum = 1.0;365 floatySum = 1.0;363 psF32 polySum = 0.0; 364 psF32 xSum = 1.0; 365 psF32 ySum = 1.0; 366 366 367 367 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 379 379 } 380 380 381 static float chebPolynomial2DEval(float x, floaty, const psPolynomial2D* myPoly)381 static psF32 chebPolynomial2DEval(psF32 x, psF32 y, const psPolynomial2D* myPoly) 382 382 { 383 383 PS_POLY_CHECK_NULL(myPoly, NAN); … … 386 386 psS32 loop_y = 0; 387 387 psS32 i = 0; 388 floatpolySum = 0.0;388 psF32 polySum = 0.0; 389 389 psPolynomial1D* *chebPolys = NULL; 390 390 psS32 maxChebyPoly = 0; … … 414 414 } 415 415 416 static float ordPolynomial3DEval(float x, float y, floatz, const psPolynomial3D* myPoly)416 static psF32 ordPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) 417 417 { 418 418 psS32 loop_x = 0; 419 419 psS32 loop_y = 0; 420 420 psS32 loop_z = 0; 421 floatpolySum = 0.0;422 floatxSum = 1.0;423 floatySum = 1.0;424 floatzSum = 1.0;421 psF32 polySum = 0.0; 422 psF32 xSum = 1.0; 423 psF32 ySum = 1.0; 424 psF32 zSum = 1.0; 425 425 426 426 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 442 442 } 443 443 444 static float chebPolynomial3DEval(float x, float y, floatz, const psPolynomial3D* myPoly)444 static psF32 chebPolynomial3DEval(psF32 x, psF32 y, psF32 z, const psPolynomial3D* myPoly) 445 445 { 446 446 psS32 loop_x = 0; … … 448 448 psS32 loop_z = 0; 449 449 psS32 i = 0; 450 floatpolySum = 0.0;450 psF32 polySum = 0.0; 451 451 psPolynomial1D* *chebPolys = NULL; 452 452 psS32 maxChebyPoly = 0; … … 483 483 } 484 484 485 static float ordPolynomial4DEval(float w, float x, float y, floatz, const psPolynomial4D* myPoly)485 static psF32 ordPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) 486 486 { 487 487 psS32 loop_w = 0; … … 489 489 psS32 loop_y = 0; 490 490 psS32 loop_z = 0; 491 floatpolySum = 0.0;492 floatwSum = 1.0;493 floatxSum = 1.0;494 floatySum = 1.0;495 floatzSum = 1.0;491 psF32 polySum = 0.0; 492 psF32 wSum = 1.0; 493 psF32 xSum = 1.0; 494 psF32 ySum = 1.0; 495 psF32 zSum = 1.0; 496 496 497 497 for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { … … 517 517 } 518 518 519 static float chebPolynomial4DEval(float w, float x, float y, floatz, const psPolynomial4D* myPoly)519 static psF32 chebPolynomial4DEval(psF32 w, psF32 x, psF32 y, psF32 z, const psPolynomial4D* myPoly) 520 520 { 521 521 psS32 loop_w = 0; … … 524 524 psS32 loop_z = 0; 525 525 psS32 i = 0; 526 floatpolySum = 0.0;526 psF32 polySum = 0.0; 527 527 psPolynomial1D* *chebPolys = NULL; 528 528 psS32 maxChebyPoly = 0; … … 568 568 Polynomial coefficients will be accessed in [w][x][y][z] fashion. 569 569 *****************************************************************************/ 570 static double dOrdPolynomial1DEval(doublex, const psDPolynomial1D* myPoly)570 static psF64 dOrdPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) 571 571 { 572 572 psS32 loop_x = 0; 573 doublepolySum = 0.0;574 doublexSum = 1.0;573 psF64 polySum = 0.0; 574 psF64 xSum = 1.0; 575 575 576 576 for (loop_x = 0; loop_x < myPoly->n; loop_x++) { … … 586 586 // XXX: You can do this without having to psAlloc() vector d. 587 587 // XXX: How does the mask vector effect Crenshaw's formula? 588 static double dChebPolynomial1DEval(doublex, const psDPolynomial1D* myPoly)588 static psF64 dChebPolynomial1DEval(psF64 x, const psDPolynomial1D* myPoly) 589 589 { 590 590 psVector *d; 591 591 psS32 n; 592 592 psS32 i; 593 doubletmp;593 psF64 tmp; 594 594 595 595 n = myPoly->n; … … 611 611 } 612 612 613 static double dOrdPolynomial2DEval(double x, doubley, const psDPolynomial2D* myPoly)613 static psF64 dOrdPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) 614 614 { 615 615 psS32 loop_x = 0; 616 616 psS32 loop_y = 0; 617 doublepolySum = 0.0;618 doublexSum = 1.0;619 doubleySum = 1.0;617 psF64 polySum = 0.0; 618 psF64 xSum = 1.0; 619 psF64 ySum = 1.0; 620 620 621 621 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 633 633 } 634 634 635 static double dChebPolynomial2DEval(double x, doubley, const psDPolynomial2D* myPoly)635 static psF64 dChebPolynomial2DEval(psF64 x, psF64 y, const psDPolynomial2D* myPoly) 636 636 { 637 637 psS32 loop_x = 0; 638 638 psS32 loop_y = 0; 639 639 psS32 i = 0; 640 doublepolySum = 0.0;640 psF64 polySum = 0.0; 641 641 psPolynomial1D* *chebPolys = NULL; 642 642 psS32 maxChebyPoly = 0; … … 667 667 } 668 668 669 static double dOrdPolynomial3DEval(double x, double y, doublez, const psDPolynomial3D* myPoly)669 static psF64 dOrdPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) 670 670 { 671 671 psS32 loop_x = 0; 672 672 psS32 loop_y = 0; 673 673 psS32 loop_z = 0; 674 doublepolySum = 0.0;675 doublexSum = 1.0;676 doubleySum = 1.0;677 doublezSum = 1.0;674 psF64 polySum = 0.0; 675 psF64 xSum = 1.0; 676 psF64 ySum = 1.0; 677 psF64 zSum = 1.0; 678 678 679 679 for (loop_x = 0; loop_x < myPoly->nX; loop_x++) { … … 695 695 } 696 696 697 static double dChebPolynomial3DEval(double x, double y, doublez, const psDPolynomial3D* myPoly)697 static psF64 dChebPolynomial3DEval(psF64 x, psF64 y, psF64 z, const psDPolynomial3D* myPoly) 698 698 { 699 699 psS32 loop_x = 0; … … 701 701 psS32 loop_z = 0; 702 702 psS32 i = 0; 703 doublepolySum = 0.0;703 psF64 polySum = 0.0; 704 704 psPolynomial1D* *chebPolys = NULL; 705 705 psS32 maxChebyPoly = 0; … … 736 736 } 737 737 738 static double dOrdPolynomial4DEval(double w, double x, double y, doublez, const psDPolynomial4D* myPoly)738 static psF64 dOrdPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) 739 739 { 740 740 psS32 loop_w = 0; … … 742 742 psS32 loop_y = 0; 743 743 psS32 loop_z = 0; 744 doublepolySum = 0.0;745 doublewSum = 1.0;746 doublexSum = 1.0;747 doubleySum = 1.0;748 doublezSum = 1.0;744 psF64 polySum = 0.0; 745 psF64 wSum = 1.0; 746 psF64 xSum = 1.0; 747 psF64 ySum = 1.0; 748 psF64 zSum = 1.0; 749 749 750 750 for (loop_w = 0; loop_w < myPoly->nW; loop_w++) { … … 770 770 } 771 771 772 static double dChebPolynomial4DEval(double w, double x, double y, doublez, const psDPolynomial4D* myPoly)772 static psF64 dChebPolynomial4DEval(psF64 w, psF64 x, psF64 y, psF64 z, const psDPolynomial4D* myPoly) 773 773 { 774 774 psS32 loop_w = 0; … … 777 777 psS32 loop_z = 0; 778 778 psS32 i = 0; 779 doublepolySum = 0.0;779 psF64 polySum = 0.0; 780 780 psPolynomial1D* *chebPolys = NULL; 781 781 psS32 maxChebyPoly = 0; … … 828 828 *****************************************************************************/ 829 829 #define FUNC_MACRO_FULL_INTERPOLATE_1D(TYPE) \ 830 static floatfullInterpolate1D##TYPE(ps##TYPE *domain, \830 static psF32 fullInterpolate1D##TYPE(ps##TYPE *domain, \ 831 831 ps##TYPE *range, \ 832 832 psS32 n, \ … … 896 896 LaGrange interpolation. 897 897 *****************************************************************************/ 898 static float interpolate1DF32(float*domain,899 float*range,898 static psF32 interpolate1DF32(psF32 *domain, 899 psF32 *range, 900 900 psS32 n, 901 901 psS32 order, 902 floatx)902 psF32 x) 903 903 { 904 904 psS32 binNum; … … 941 941 evaluated Gaussian is: \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] 942 942 *****************************************************************************/ 943 float psGaussian(float x, float mean, floatsigma, psBool normal)944 { 945 floattmp = 1.0;943 psF32 psGaussian(psF32 x, psF32 mean, psF32 sigma, psBool normal) 944 { 945 psF32 tmp = 1.0; 946 946 947 947 psTrace(".psLib.dataManip.psFunctions.psGaussian", 4, … … 967 967 XXX: There is no way to seed the random generator. 968 968 *****************************************************************************/ 969 psVector* p_psGaussianDev( float mean, floatsigma, psS32 Npts)969 psVector* p_psGaussianDev(psF32 mean, psF32 sigma, psS32 Npts) 970 970 { 971 971 PS_INT_CHECK_NON_NEGATIVE(Npts, NULL); … … 1007 1007 newPoly->type = type; 1008 1008 newPoly->n = n; 1009 newPoly->coeff = ( float *)psAlloc(n * sizeof(float));1010 newPoly->coeffErr = ( float *)psAlloc(n * sizeof(float));1011 newPoly->mask = ( char *)psAlloc(n * sizeof(char));1009 newPoly->coeff = (psF32 *)psAlloc(n * sizeof(psF32)); 1010 newPoly->coeffErr = (psF32 *)psAlloc(n * sizeof(psF32)); 1011 newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); 1012 1012 for (i = 0; i < n; i++) { 1013 1013 newPoly->coeff[i] = 0.0; … … 1036 1036 newPoly->nY = nY; 1037 1037 1038 newPoly->coeff = ( float **)psAlloc(nX * sizeof(float*));1039 newPoly->coeffErr = ( float **)psAlloc(nX * sizeof(float*));1040 newPoly->mask = ( char **)psAlloc(nX * sizeof(char*));1038 newPoly->coeff = (psF32 **)psAlloc(nX * sizeof(psF32 *)); 1039 newPoly->coeffErr = (psF32 **)psAlloc(nX * sizeof(psF32 *)); 1040 newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); 1041 1041 for (x = 0; x < nX; x++) { 1042 newPoly->coeff[x] = ( float *)psAlloc(nY * sizeof(float));1043 newPoly->coeffErr[x] = ( float *)psAlloc(nY * sizeof(float));1044 newPoly->mask[x] = ( char *)psAlloc(nY * sizeof(char));1042 newPoly->coeff[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); 1043 newPoly->coeffErr[x] = (psF32 *)psAlloc(nY * sizeof(psF32)); 1044 newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); 1045 1045 } 1046 1046 for (x = 0; x < nX; x++) { … … 1075 1075 newPoly->nZ = nZ; 1076 1076 1077 newPoly->coeff = ( float ***)psAlloc(nX * sizeof(float**));1078 newPoly->coeffErr = ( float ***)psAlloc(nX * sizeof(float**));1079 newPoly->mask = ( char ***)psAlloc(nX * sizeof(char**));1077 newPoly->coeff = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1078 newPoly->coeffErr = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1079 newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1080 1080 for (x = 0; x < nX; x++) { 1081 newPoly->coeff[x] = ( float **)psAlloc(nY * sizeof(float*));1082 newPoly->coeffErr[x] = ( float **)psAlloc(nY * sizeof(float*));1083 newPoly->mask[x] = ( char **)psAlloc(nY * sizeof(char*));1081 newPoly->coeff[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1082 newPoly->coeffErr[x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1083 newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1084 1084 for (y = 0; y < nY; y++) { 1085 newPoly->coeff[x][y] = ( float *)psAlloc(nZ * sizeof(float));1086 newPoly->coeffErr[x][y] = ( float *)psAlloc(nZ * sizeof(float));1087 newPoly->mask[x][y] = ( char *)psAlloc(nZ * sizeof(char));1085 newPoly->coeff[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1086 newPoly->coeffErr[x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1087 newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1088 1088 } 1089 1089 } … … 1124 1124 newPoly->nZ = nZ; 1125 1125 1126 newPoly->coeff = ( float ****)psAlloc(nW * sizeof(float***));1127 newPoly->coeffErr = ( float ****)psAlloc(nW * sizeof(float***));1128 newPoly->mask = ( char ****)psAlloc(nW * sizeof(char***));1126 newPoly->coeff = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); 1127 newPoly->coeffErr = (psF32 ****)psAlloc(nW * sizeof(psF32 ***)); 1128 newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); 1129 1129 for (w = 0; w < nW; w++) { 1130 newPoly->coeff[w] = ( float ***)psAlloc(nX * sizeof(float**));1131 newPoly->coeffErr[w] = ( float ***)psAlloc(nX * sizeof(float**));1132 newPoly->mask[w] = ( char ***)psAlloc(nX * sizeof(char**));1130 newPoly->coeff[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1131 newPoly->coeffErr[w] = (psF32 ***)psAlloc(nX * sizeof(psF32 **)); 1132 newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1133 1133 for (x = 0; x < nX; x++) { 1134 newPoly->coeff[w][x] = ( float **)psAlloc(nY * sizeof(float*));1135 newPoly->coeffErr[w][x] = ( float **)psAlloc(nY * sizeof(float*));1136 newPoly->mask[w][x] = ( char **)psAlloc(nY * sizeof(char*));1134 newPoly->coeff[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1135 newPoly->coeffErr[w][x] = (psF32 **)psAlloc(nY * sizeof(psF32 *)); 1136 newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1137 1137 for (y = 0; y < nY; y++) { 1138 newPoly->coeff[w][x][y] = ( float *)psAlloc(nZ * sizeof(float));1139 newPoly->coeffErr[w][x][y] = ( float *)psAlloc(nZ * sizeof(float));1140 newPoly->mask[w][x][y] = ( char *)psAlloc(nZ * sizeof(char));1138 newPoly->coeff[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1139 newPoly->coeffErr[w][x][y] = (psF32 *)psAlloc(nZ * sizeof(psF32)); 1140 newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1141 1141 } 1142 1142 } … … 1157 1157 } 1158 1158 1159 float psPolynomial1DEval(const psPolynomial1D* myPoly, floatx)1159 psF32 psPolynomial1DEval(const psPolynomial1D* myPoly, psF32 x) 1160 1160 { 1161 1161 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1195 1195 } 1196 1196 1197 float psPolynomial2DEval(const psPolynomial2D* myPoly, float x, floaty)1197 psF32 psPolynomial2DEval(const psPolynomial2D* myPoly, psF32 x, psF32 y) 1198 1198 { 1199 1199 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1250 1250 } 1251 1251 1252 float psPolynomial3DEval(const psPolynomial3D* myPoly, float x, float y, floatz)1252 psF32 psPolynomial3DEval(const psPolynomial3D* myPoly, psF32 x, psF32 y, psF32 z) 1253 1253 { 1254 1254 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1316 1316 } 1317 1317 1318 float psPolynomial4DEval(const psPolynomial4D* myPoly, float w, float x, float y, floatz)1318 psF32 psPolynomial4DEval(const psPolynomial4D* myPoly, psF32 w, psF32 x, psF32 y, psF32 z) 1319 1319 { 1320 1320 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1406 1406 newPoly->type = type; 1407 1407 newPoly->n = n; 1408 newPoly->coeff = ( double *)psAlloc(n * sizeof(double));1409 newPoly->coeffErr = ( double *)psAlloc(n * sizeof(double));1410 newPoly->mask = ( char *)psAlloc(n * sizeof(char));1408 newPoly->coeff = (psF64 *)psAlloc(n * sizeof(psF64)); 1409 newPoly->coeffErr = (psF64 *)psAlloc(n * sizeof(psF64)); 1410 newPoly->mask = (psU8 *)psAlloc(n * sizeof(psU8)); 1411 1411 for (i = 0; i < n; i++) { 1412 1412 newPoly->coeff[i] = 0.0; … … 1435 1435 newPoly->nY = nY; 1436 1436 1437 newPoly->coeff = ( double **)psAlloc(nX * sizeof(double*));1438 newPoly->coeffErr = ( double **)psAlloc(nX * sizeof(double*));1439 newPoly->mask = ( char **)psAlloc(nX * sizeof(char*));1437 newPoly->coeff = (psF64 **)psAlloc(nX * sizeof(psF64 *)); 1438 newPoly->coeffErr = (psF64 **)psAlloc(nX * sizeof(psF64 *)); 1439 newPoly->mask = (psU8 **)psAlloc(nX * sizeof(psU8 *)); 1440 1440 for (x = 0; x < nX; x++) { 1441 newPoly->coeff[x] = ( double *)psAlloc(nY * sizeof(double));1442 newPoly->coeffErr[x] = ( double *)psAlloc(nY * sizeof(double));1443 newPoly->mask[x] = ( char *)psAlloc(nY * sizeof(char));1441 newPoly->coeff[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); 1442 newPoly->coeffErr[x] = (psF64 *)psAlloc(nY * sizeof(psF64)); 1443 newPoly->mask[x] = (psU8 *)psAlloc(nY * sizeof(psU8)); 1444 1444 } 1445 1445 for (x = 0; x < nX; x++) { … … 1474 1474 newPoly->nZ = nZ; 1475 1475 1476 newPoly->coeff = ( double ***)psAlloc(nX * sizeof(double**));1477 newPoly->coeffErr = ( double ***)psAlloc(nX * sizeof(double**));1478 newPoly->mask = ( char ***)psAlloc(nX * sizeof(char**));1476 newPoly->coeff = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1477 newPoly->coeffErr = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1478 newPoly->mask = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1479 1479 for (x = 0; x < nX; x++) { 1480 newPoly->coeff[x] = ( double **)psAlloc(nY * sizeof(double*));1481 newPoly->coeffErr[x] = ( double **)psAlloc(nY * sizeof(double*));1482 newPoly->mask[x] = ( char **)psAlloc(nY * sizeof(char*));1480 newPoly->coeff[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1481 newPoly->coeffErr[x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1482 newPoly->mask[x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1483 1483 for (y = 0; y < nY; y++) { 1484 newPoly->coeff[x][y] = ( double *)psAlloc(nZ * sizeof(double));1485 newPoly->coeffErr[x][y] = ( double *)psAlloc(nZ * sizeof(double));1486 newPoly->mask[x][y] = ( char *)psAlloc(nZ * sizeof(char));1484 newPoly->coeff[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1485 newPoly->coeffErr[x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1486 newPoly->mask[x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1487 1487 } 1488 1488 } … … 1523 1523 newPoly->nZ = nZ; 1524 1524 1525 newPoly->coeff = ( double ****)psAlloc(nW * sizeof(double***));1526 newPoly->coeffErr = ( double ****)psAlloc(nW * sizeof(double***));1527 newPoly->mask = ( char ****)psAlloc(nW * sizeof(char***));1525 newPoly->coeff = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); 1526 newPoly->coeffErr = (psF64 ****)psAlloc(nW * sizeof(psF64 ***)); 1527 newPoly->mask = (psU8 ****)psAlloc(nW * sizeof(psU8 ***)); 1528 1528 for (w = 0; w < nW; w++) { 1529 newPoly->coeff[w] = ( double ***)psAlloc(nX * sizeof(double**));1530 newPoly->coeffErr[w] = ( double ***)psAlloc(nX * sizeof(double**));1531 newPoly->mask[w] = ( char ***)psAlloc(nX * sizeof(char**));1529 newPoly->coeff[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1530 newPoly->coeffErr[w] = (psF64 ***)psAlloc(nX * sizeof(psF64 **)); 1531 newPoly->mask[w] = (psU8 ***)psAlloc(nX * sizeof(psU8 **)); 1532 1532 for (x = 0; x < nX; x++) { 1533 newPoly->coeff[w][x] = ( double **)psAlloc(nY * sizeof(double*));1534 newPoly->coeffErr[w][x] = ( double **)psAlloc(nY * sizeof(double*));1535 newPoly->mask[w][x] = ( char **)psAlloc(nY * sizeof(char*));1533 newPoly->coeff[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1534 newPoly->coeffErr[w][x] = (psF64 **)psAlloc(nY * sizeof(psF64 *)); 1535 newPoly->mask[w][x] = (psU8 **)psAlloc(nY * sizeof(psU8 *)); 1536 1536 for (y = 0; y < nY; y++) { 1537 newPoly->coeff[w][x][y] = ( double *)psAlloc(nZ * sizeof(double));1538 newPoly->coeffErr[w][x][y] = ( double *)psAlloc(nZ * sizeof(double));1539 newPoly->mask[w][x][y] = ( char *)psAlloc(nZ * sizeof(char));1537 newPoly->coeff[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1538 newPoly->coeffErr[w][x][y] = (psF64 *)psAlloc(nZ * sizeof(psF64)); 1539 newPoly->mask[w][x][y] = (psU8 *)psAlloc(nZ * sizeof(psU8)); 1540 1540 } 1541 1541 } … … 1557 1557 1558 1558 1559 double psDPolynomial1DEval(const psDPolynomial1D* myPoly, doublex)1559 psF64 psDPolynomial1DEval(const psDPolynomial1D* myPoly, psF64 x) 1560 1560 { 1561 1561 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1598 1598 1599 1599 1600 doublepsDPolynomial2DEval(const psDPolynomial2D* myPoly,1601 doublex,1602 doubley)1600 psF64 psDPolynomial2DEval(const psDPolynomial2D* myPoly, 1601 psF64 x, 1602 psF64 y) 1603 1603 { 1604 1604 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1655 1655 1656 1656 1657 doublepsDPolynomial3DEval(const psDPolynomial3D* myPoly,1658 doublex,1659 doubley,1660 doublez)1657 psF64 psDPolynomial3DEval(const psDPolynomial3D* myPoly, 1658 psF64 x, 1659 psF64 y, 1660 psF64 z) 1661 1661 { 1662 1662 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1724 1724 } 1725 1725 1726 doublepsDPolynomial4DEval(const psDPolynomial4D* myPoly,1727 doublew,1728 doublex,1729 doubley,1730 doublez)1726 psF64 psDPolynomial4DEval(const psDPolynomial4D* myPoly, 1727 psF64 w, 1728 psF64 x, 1729 psF64 y, 1730 psF64 z) 1731 1731 { 1732 1732 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1810 1810 // psS32 n; 1811 1811 // psPolynomial1D **spline; 1812 // float*p_psDeriv2;1813 // float*domains;1812 // psF32 *p_psDeriv2; 1813 // psF32 *domains; 1814 1814 //} psSpline1D; 1815 1815 … … 1822 1822 psSpline1D *psSpline1DAlloc(psS32 numSplines, 1823 1823 psS32 order, 1824 floatmin,1825 floatmax)1824 psF32 min, 1825 psF32 max) 1826 1826 { 1827 1827 PS_INT_CHECK_NON_NEGATIVE(numSplines, NULL); … … 1831 1831 psSpline1D *tmp = NULL; 1832 1832 psS32 i; 1833 floattmpDomain;1834 floatwidth;1833 psF32 tmpDomain; 1834 psF32 width; 1835 1835 1836 1836 tmp = (psSpline1D *) psAlloc(sizeof(psSpline1D)); … … 1845 1845 tmp->p_psDeriv2 = NULL; 1846 1846 1847 tmp->domains = ( float *) psAlloc((numSplines+1) * sizeof(float));1848 width = (max - min) / (( float) numSplines);1847 tmp->domains = (psF32 *) psAlloc((numSplines+1) * sizeof(psF32)); 1848 width = (max - min) / ((psF32) numSplines); 1849 1849 1850 1850 (tmp->domains)[0] = min; … … 1884 1884 } 1885 1885 1886 tmp->domains = ( float *) psAlloc((bounds->n) * sizeof(float));1886 tmp->domains = (psF32 *) psAlloc((bounds->n) * sizeof(psF32)); 1887 1887 1888 1888 for (i=0;i<bounds->n;i++) { … … 1977 1977 { 1978 1978 PS_PTR_CHECK_TYPE_EQUAL(x, bins, -3); 1979 char* strType;1979 psS8* strType; 1980 1980 1981 1981 switch (x->type.type) { … … 2074 2074 domain32 = psVectorCopy(domain32, domain, PS_TYPE_F32); 2075 2075 2076 psScalar *tmpScalar = psScalarAlloc(( double)2076 psScalar *tmpScalar = psScalarAlloc((psF64) 2077 2077 interpolate1DF32(domain32->data.F32, 2078 2078 range32->data.F32, 2079 2079 domain32->n, 2080 2080 order, 2081 ( float) x->data.F64), PS_TYPE_F64);2081 (psF32) x->data.F64), PS_TYPE_F64); 2082 2082 psFree(range32); 2083 2083 psFree(domain32); … … 2089 2089 2090 2090 } else { 2091 char* strType;2091 psS8* strType; 2092 2092 PS_TYPE_NAME(strType,x->type.type); 2093 2093 psError(PS_ERR_BAD_PARAMETER_TYPE, … … 2115 2115 the spline fit functions require F32 and F64. 2116 2116 *****************************************************************************/ 2117 floatpsSpline1DEval(2118 floatx,2117 psF32 psSpline1DEval( 2118 psF32 x, 2119 2119 const psSpline1D *spline 2120 2120 ) … … 2172 2172 for (i=0;i<x->n;i++) { 2173 2173 tmpVector->data.F32[i] = psSpline1DEval( 2174 ( float) x->data.F64[i],2174 (psF32) x->data.F64[i], 2175 2175 spline 2176 2176 ); 2177 2177 } 2178 2178 } else { 2179 char* strType;2179 psS8* strType; 2180 2180 PS_TYPE_NAME(strType,x->type.type); 2181 2181 psError(PS_ERR_BAD_PARAMETER_TYPE, -
trunk/psLib/src/math/psSpline.h
r2600 r2788 12 12 * @author GLG, MHPCC 13 13 * 14 * @version $Revision: 1.3 7$ $Name: not supported by cvs2svn $15 * @date $Date: 2004-12- 02 21:12:51$14 * @version $Revision: 1.38 $ $Name: not supported by cvs2svn $ 15 * @date $Date: 2004-12-22 05:09:32 $ 16 16 * 17 17 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 39 39 * \f[ exp(-\frac{(x-mean)^2}{2\sigma^2}) \f] 40 40 * 41 * @return floatvalue on the gaussian curve given the input parameters42 */ 43 floatpsGaussian(44 floatx, ///< Value at which to evaluate45 floatmean, ///< Mean for the Gaussian46 floatstddev, ///< Standard deviation for the Gaussian41 * @return psF32 value on the gaussian curve given the input parameters 42 */ 43 psF32 psGaussian( 44 psF32 x, ///< Value at which to evaluate 45 psF32 mean, ///< Mean for the Gaussian 46 psF32 stddev, ///< Standard deviation for the Gaussian 47 47 psBool normal ///< Indicates whether result should be normalized 48 48 ); … … 55 55 */ 56 56 psVector* p_psGaussianDev( 57 floatmean, ///< The mean of the Gaussian58 floatsigma, ///< The sigma of the Gaussian57 psF32 mean, ///< The mean of the Gaussian 58 psF32 sigma, ///< The sigma of the Gaussian 59 59 psS32 Npts ///< The size of the vector 60 60 ); … … 70 70 psPolynomialType type; ///< Polynomial type 71 71 psS32 n; ///< Number of terms 72 float*coeff; ///< Coefficients73 float*coeffErr; ///< Error in coefficients74 char*mask; ///< Coefficient mask72 psF32 *coeff; ///< Coefficients 73 psF32 *coeffErr; ///< Error in coefficients 74 psU8 *mask; ///< Coefficient mask 75 75 } 76 76 psPolynomial1D; … … 82 82 psS32 nX; ///< Number of terms in x 83 83 psS32 nY; ///< Number of terms in y 84 float**coeff; ///< Coefficients85 float**coeffErr; ///< Error in coefficients86 char**mask; ///< Coefficients mask84 psF32 **coeff; ///< Coefficients 85 psF32 **coeffErr; ///< Error in coefficients 86 psU8 **mask; ///< Coefficients mask 87 87 } 88 88 psPolynomial2D; … … 95 95 psS32 nY; ///< Number of terms in y 96 96 psS32 nZ; ///< Number of terms in z 97 float***coeff; ///< Coefficients98 float***coeffErr; ///< Error in coefficients99 char***mask; ///< Coefficients mask97 psF32 ***coeff; ///< Coefficients 98 psF32 ***coeffErr; ///< Error in coefficients 99 psU8 ***mask; ///< Coefficients mask 100 100 } 101 101 psPolynomial3D; … … 109 109 psS32 nY; ///< Number of terms in y 110 110 psS32 nZ; ///< Number of terms in z 111 float****coeff; ///< Coefficients112 float****coeffErr; ///< Error in coefficients113 char****mask; ///< Coefficients mask111 psF32 ****coeff; ///< Coefficients 112 psF32 ****coeffErr; ///< Error in coefficients 113 psU8 ****mask; ///< Coefficients mask 114 114 } 115 115 psPolynomial4D; … … 160 160 /** Evaluates a 1-D polynomial at specific coordinates. 161 161 * 162 * @return floatresult of polynomial at given location163 */ 164 floatpsPolynomial1DEval(162 * @return psF32 result of polynomial at given location 163 */ 164 psF32 psPolynomial1DEval( 165 165 const psPolynomial1D* myPoly, ///< Coefficients for the polynomial 166 floatx ///< location at which to evaluate166 psF32 x ///< location at which to evaluate 167 167 ); 168 168 169 169 /** Evaluates a 2-D polynomial at specific coordinates. 170 170 * 171 * @return floatresult of polynomial at given location172 */ 173 floatpsPolynomial2DEval(171 * @return psF32 result of polynomial at given location 172 */ 173 psF32 psPolynomial2DEval( 174 174 const psPolynomial2D* myPoly, ///< Coefficients for the polynomial 175 floatx, ///< x location at which to evaluate176 floaty ///< y location at which to evaluate175 psF32 x, ///< x location at which to evaluate 176 psF32 y ///< y location at which to evaluate 177 177 ); 178 178 179 179 /** Evaluates a 3-D polynomial at specific coordinates. 180 180 * 181 * @return floatresult of polynomial at given location182 */ 183 floatpsPolynomial3DEval(181 * @return psF32 result of polynomial at given location 182 */ 183 psF32 psPolynomial3DEval( 184 184 const psPolynomial3D* myPoly, ///< Coefficients for the polynomial 185 floatx, ///< x location at which to evaluate186 floaty, ///< y location at which to evaluate187 floatz ///< z location at which to evaluate185 psF32 x, ///< x location at which to evaluate 186 psF32 y, ///< y location at which to evaluate 187 psF32 z ///< z location at which to evaluate 188 188 ); 189 189 190 190 /** Evaluates a 4-D polynomial at specific coordinates. 191 191 * 192 * @return floatresult of polynomial at given location193 */ 194 floatpsPolynomial4DEval(192 * @return psF32 result of polynomial at given location 193 */ 194 psF32 psPolynomial4DEval( 195 195 const psPolynomial4D* myPoly, ///< Coefficients for the polynomial 196 floatw, ///< w location at which to evaluate197 floatx, ///< x location at which to evaluate198 floaty, ///< y location at which to evaluate199 floatz ///< z location at which to evaluate196 psF32 w, ///< w location at which to evaluate 197 psF32 x, ///< x location at which to evaluate 198 psF32 y, ///< y location at which to evaluate 199 psF32 z ///< z location at which to evaluate 200 200 ); 201 201 … … 235 235 psPolynomialType type; ///< Polynomial type 236 236 psS32 n; ///< Number of terms 237 double*coeff; ///< Coefficients238 double*coeffErr; ///< Error in coefficients239 char*mask; ///< Coefficient mask237 psF64 *coeff; ///< Coefficients 238 psF64 *coeffErr; ///< Error in coefficients 239 psU8 *mask; ///< Coefficient mask 240 240 } 241 241 psDPolynomial1D; … … 247 247 psS32 nX; ///< Number of terms in x 248 248 psS32 nY; ///< Number of terms in y 249 double**coeff; ///< Coefficients250 double**coeffErr; ///< Error in coefficients251 char**mask; ///< Coefficients mask249 psF64 **coeff; ///< Coefficients 250 psF64 **coeffErr; ///< Error in coefficients 251 psU8 **mask; ///< Coefficients mask 252 252 } 253 253 psDPolynomial2D; … … 260 260 psS32 nY; ///< Number of terms in y 261 261 psS32 nZ; ///< Number of terms in z 262 double***coeff; ///< Coefficients263 double***coeffErr; ///< Error in coefficients264 char***mask; ///< Coefficient mask262 psF64 ***coeff; ///< Coefficients 263 psF64 ***coeffErr; ///< Error in coefficients 264 psU8 ***mask; ///< Coefficient mask 265 265 } 266 266 psDPolynomial3D; … … 274 274 psS32 nY; ///< Number of terms in y 275 275 psS32 nZ; ///< Number of terms in z 276 double****coeff; ///< Coefficients277 double****coeffErr; ///< Error in coefficients278 char****mask; ///< Coefficients mask276 psF64 ****coeff; ///< Coefficients 277 psF64 ****coeffErr; ///< Error in coefficients 278 psU8 ****mask; ///< Coefficients mask 279 279 } 280 280 psDPolynomial4D; … … 324 324 /** Evaluates a double-precision 1-D polynomial at specific coordinates. 325 325 * 326 * @return floatresult of polynomial at given location327 */ 328 doublepsDPolynomial1DEval(326 * @return psF32 result of polynomial at given location 327 */ 328 psF64 psDPolynomial1DEval( 329 329 const psDPolynomial1D* myPoly, ///< Coefficients for the polynomial 330 doublex ///< Value at which to evaluate330 psF64 x ///< Value at which to evaluate 331 331 ); 332 332 333 333 /** Evaluates a double-precision 2-D polynomial at specific coordinates. 334 334 * 335 * @return floatresult of polynomial at given location336 */ 337 doublepsDPolynomial2DEval(335 * @return psF32 result of polynomial at given location 336 */ 337 psF64 psDPolynomial2DEval( 338 338 const psDPolynomial2D* myPoly, ///< Coefficients for the polynomial 339 doublex, ///< Value x at which to evaluate340 doubley ///< Value y at which to evaluate339 psF64 x, ///< Value x at which to evaluate 340 psF64 y ///< Value y at which to evaluate 341 341 ); 342 342 343 343 /** Evaluates a double-precision 3-D polynomial at specific coordinates. 344 344 * 345 * @return floatresult of polynomial at given location346 */ 347 doublepsDPolynomial3DEval(345 * @return psF32 result of polynomial at given location 346 */ 347 psF64 psDPolynomial3DEval( 348 348 const psDPolynomial3D* myPoly, ///< Coefficients for the polynomial 349 doublex, ///< Value x at which to evaluate350 doubley, ///< Value y at which to evaluate351 doublez ///< Value z at which to evaluate349 psF64 x, ///< Value x at which to evaluate 350 psF64 y, ///< Value y at which to evaluate 351 psF64 z ///< Value z at which to evaluate 352 352 ); 353 353 354 354 /** Evaluates a double-precision 4-D polynomial at specific coordinates. 355 355 * 356 * @return floatresult of polynomial at given location357 */ 358 doublepsDPolynomial4DEval(356 * @return psF32 result of polynomial at given location 357 */ 358 psF64 psDPolynomial4DEval( 359 359 const psDPolynomial4D* myPoly, ///< Coefficients for the polynomial 360 doublew, ///< Value w at which to evaluate361 doublex, ///< Value x at which to evaluate362 doubley, ///< Value y at which to evaluate363 doublez ///< Value z at which to evaluate360 psF64 w, ///< Value w at which to evaluate 361 psF64 x, ///< Value x at which to evaluate 362 psF64 y, ///< Value y at which to evaluate 363 psF64 z ///< Value z at which to evaluate 364 364 ); 365 365 … … 396 396 psS32 n; ///< The number of spline polynomials 397 397 psPolynomial1D **spline; ///< An array of n pointers to the spline polynomials 398 float*p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1.399 float*domains; ///< The boundaries between each spline piece. Size is n+1.398 psF32 *p_psDeriv2; ///< For cubic splines, the second derivative at each domain point. Size is n+1. 399 psF32 *domains; ///< The boundaries between each spline piece. Size is n+1. 400 400 } 401 401 psSpline1D; … … 403 403 psSpline1D *psSpline1DAlloc(psS32 n, 404 404 psS32 order, 405 floatmin,406 floatmax);405 psF32 min, 406 psF32 max); 407 407 408 408 psSpline1D *psSpline1DAllocGeneric(const psVector *bounds, 409 409 psS32 order); 410 410 411 floatpsSpline1DEval(412 floatx,411 psF32 psSpline1DEval( 412 psF32 x, 413 413 const psSpline1D *spline 414 414 ); … … 427 427 psScalar *x); 428 428 429 floatp_psNRSpline1DEval(psSpline1D *spline,429 psF32 p_psNRSpline1DEval(psSpline1D *spline, 430 430 const psVector* restrict x, 431 431 const psVector* restrict y, 432 floatX);432 psF32 X); 433 433 434 434 /* \} */// End of MathGroup Functions -
trunk/psLib/src/math/psStats.c
r2782 r2788 9 9 * @author GLG, MHPCC 10 10 * 11 * @version $Revision: 1.10 7$ $Name: not supported by cvs2svn $12 * @date $Date: 2004-12-22 0 0:54:28$11 * @version $Revision: 1.108 $ $Name: not supported by cvs2svn $ 12 * @date $Date: 2004-12-22 05:09:32 $ 13 13 * 14 14 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 70 70 /*****************************************************************************/ 71 71 72 psBool p_psGetStatValue(const psStats* stats, double*value)72 psBool p_psGetStatValue(const psStats* stats, psF64 *value) 73 73 { 74 74 … … 145 145 this routine sets stats->sampleMean to NAN. 146 146 *****************************************************************************/ 147 intp_psVectorSampleMean(const psVector* restrict myVector,148 const psVector* restrict errors,149 const psVector* restrict maskVector,150 psU32 maskVal,151 psStats* stats)147 psS32 p_psVectorSampleMean(const psVector* restrict myVector, 148 const psVector* restrict errors, 149 const psVector* restrict maskVector, 150 psU32 maskVal, 151 psStats* stats) 152 152 { 153 153 154 154 psS32 i = 0; // Loop index variable 155 floatmean = 0.0; // The mean155 psF32 mean = 0.0; // The mean 156 156 psS32 count = 0; // # of points in this mean 157 157 … … 171 171 } 172 172 if (count != 0) { 173 mean /= ( float)count;173 mean /= (psF32)count; 174 174 } else { 175 175 mean = NAN; … … 185 185 } 186 186 if (count != 0) { 187 mean /= ( float)count;187 mean /= (psF32)count; 188 188 } else { 189 189 mean = NAN; … … 199 199 } 200 200 if (count != 0) { 201 mean /= ( float)count;201 mean /= (psF32)count; 202 202 } else { 203 203 mean = NAN; … … 207 207 mean += myVector->data.F32[i]; 208 208 } 209 mean /= ( float)myVector->n;209 mean /= (psF32)myVector->n; 210 210 } 211 211 } … … 283 283 this routine sets stats->max to NAN. 284 284 *****************************************************************************/ 285 intp_psVectorMax(const psVector* restrict myVector,286 const psVector* restrict maskVector,287 psU32 maskVal,288 psStats* stats)285 psS32 p_psVectorMax(const psVector* restrict myVector, 286 const psVector* restrict maskVector, 287 psU32 maskVal, 288 psStats* stats) 289 289 { 290 290 psS32 i = 0; // Loop index variable 291 floatmax = -PS_MAX_F32; // The calculated maximum291 psF32 max = -PS_MAX_F32; // The calculated maximum 292 292 psS32 empty = true; // Does this vector have valid elements? 293 293 … … 348 348 this routine sets stats->min to NAN. 349 349 *****************************************************************************/ 350 intp_psVectorMin(const psVector* restrict myVector,351 const psVector* restrict maskVector,352 psU32 maskVal,353 psStats* stats)350 psS32 p_psVectorMin(const psVector* restrict myVector, 351 const psVector* restrict maskVector, 352 psU32 maskVal, 353 psStats* stats) 354 354 { 355 355 psS32 i = 0; // Loop index variable 356 floatmin = PS_MAX_F32; // The calculated maximum356 psF32 min = PS_MAX_F32; // The calculated maximum 357 357 psS32 empty = true; // Does this vector have valid elements? 358 358 … … 612 612 *****************************************************************************/ 613 613 psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, 614 floatsigma)614 psF32 sigma) 615 615 { 616 616 PS_PTR_CHECK_NULL(robustHistogram, NULL); … … 619 619 psS32 i = 0; // Loop index variable 620 620 psS32 j = 0; // Loop index variable 621 floatiMid;622 floatjMid;621 psF32 iMid; 622 psF32 jMid; 623 623 psS32 numBins = robustHistogram->nums->n; 624 624 psS32 numBounds = robustHistogram->bounds->n; … … 626 626 psS32 jMin = 0; 627 627 psS32 jMax = 0; 628 floatfirstBound = robustHistogram->bounds->data.F32[0];629 floatlastBound = robustHistogram->bounds->data.F32[numBounds-1];628 psF32 firstBound = robustHistogram->bounds->data.F32[0]; 629 psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; 630 630 psScalar x; 631 631 … … 795 795 psS32 i = 0; // Loop index variable 796 796 psS32 countInt = 0; // # of data points being used 797 floatcountFloat = 0.0; // # of data points being used798 floatmean = 0.0; // The mean799 floatdiff = 0.0; // Used in calculating stdev800 floatsumSquares = 0.0; // temporary variable801 floatsumDiffs = 0.0; // temporary variable797 psF32 countFloat = 0.0; // # of data points being used 798 psF32 mean = 0.0; // The mean 799 psF32 diff = 0.0; // Used in calculating stdev 800 psF32 sumSquares = 0.0; // temporary variable 801 psF32 sumDiffs = 0.0; // temporary variable 802 802 803 803 // This procedure requires the mean. If it has not been already … … 865 865 psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", countInt); 866 866 } else { 867 countFloat = ( float)countInt;867 countFloat = (psF32)countInt; 868 868 stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); 869 869 } … … 889 889 psS32 i = 0; // Loop index variable 890 890 psS32 countInt = 0; // # of data points being used 891 floatcountFloat = 0.0; // # of data points being used892 floatmean = 0.0; // The mean893 floatdiff = 0.0; // Used in calculating stdev894 floatsumSquares = 0.0; // temporary variable895 floatsumDiffs = 0.0; // temporary variable891 psF32 countFloat = 0.0; // # of data points being used 892 psF32 mean = 0.0; // The mean 893 psF32 diff = 0.0; // Used in calculating stdev 894 psF32 sumSquares = 0.0; // temporary variable 895 psF32 sumDiffs = 0.0; // temporary variable 896 896 // psF32 sum1; 897 897 // psF32 sum2; … … 982 982 stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor)); 983 983 } else { 984 countFloat = ( float)countInt;984 countFloat = (psF32)countInt; 985 985 stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); 986 986 … … 1003 1003 -2: warning 1004 1004 *****************************************************************************/ 1005 intp_psVectorClippedStats(const psVector* restrict myVector,1006 const psVector* restrict errors,1007 const psVector* restrict maskVector,1008 psU32 maskVal,1009 psStats* stats)1005 psS32 p_psVectorClippedStats(const psVector* restrict myVector, 1006 const psVector* restrict errors, 1007 const psVector* restrict maskVector, 1008 psU32 maskVal, 1009 psStats* stats) 1010 1010 { 1011 1011 psS32 i = 0; // Loop index variable 1012 1012 psS32 j = 0; // Loop index variable 1013 floatclippedMean = 0.0; // self-explanatory1014 floatclippedStdev = 0.0; // self-explanatory1015 floatoldStanMean = 0.0; // Temporary variable1016 floatoldStanStdev = 0.0; // Temporary variable1013 psF32 clippedMean = 0.0; // self-explanatory 1014 psF32 clippedStdev = 0.0; // self-explanatory 1015 psF32 oldStanMean = 0.0; // Temporary variable 1016 psF32 oldStanStdev = 0.0; // Temporary variable 1017 1017 psVector* tmpMask = NULL; // Temporary vector 1018 1018 … … 1240 1240 XXX: Create a 2nd-order polynomial version and solve for X analytically. 1241 1241 *****************************************************************************/ 1242 floatp_ps1DPolyMedian(psPolynomial1D* myPoly,1243 floatrangeLow,1244 floatrangeHigh,1245 floatgetThisValue)1242 psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly, 1243 psF32 rangeLow, 1244 psF32 rangeHigh, 1245 psF32 getThisValue) 1246 1246 { 1247 1247 PS_POLY_CHECK_NULL(myPoly, NAN); … … 1250 1250 // falls within the range of y-values of the polynomial "myPoly" in the 1251 1251 // specified x-range (rangeLow:rangeHigh). 1252 floatfLo = psPolynomial1DEval(1252 psF32 fLo = psPolynomial1DEval( 1253 1253 myPoly, 1254 1254 rangeLow 1255 1255 ); 1256 floatfHi = psPolynomial1DEval(1256 psF32 fHi = psPolynomial1DEval( 1257 1257 myPoly, 1258 1258 rangeHigh … … 1266 1266 1267 1267 psS32 numIterations = 0; 1268 floatmidpoint = 0.0;1269 floatoldMidpoint = 1.0;1270 floatf = 0.0;1268 psF32 midpoint = 0.0; 1269 psF32 oldMidpoint = 1.0; 1270 psF32 f = 0.0; 1271 1271 1272 1272 while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) { … … 1309 1309 tests to ensure that binNum is within acceptable ranges for both vectors. 1310 1310 *****************************************************************************/ 1311 floatfitQuadraticSearchForYThenReturnX(psVector *xVec,1311 psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, 1312 1312 psVector *yVec, 1313 1313 psS32 binNum, 1314 floatyVal)1314 psF32 yVal) 1315 1315 { 1316 1316 PS_VECTOR_CHECK_NULL(xVec, NAN); … … 1331 1331 psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); 1332 1332 1333 floattmpFloat;1333 psF32 tmpFloat; 1334 1334 1335 1335 if ((binNum > 0) && (binNum < (yVec->n - 2))) { 1336 1336 // The general case. We have all three points. 1337 x->data.F64[0] = ( double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));1338 x->data.F64[1] = ( double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));1339 x->data.F64[2] = ( double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));1337 x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); 1338 x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); 1339 x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); 1340 1340 y->data.F64[0] = yVec->data.F32[binNum - 1]; 1341 1341 y->data.F64[1] = yVec->data.F32[binNum]; … … 1346 1346 psError(PS_ERR_BAD_PARAMETER_VALUE, true, 1347 1347 PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, 1348 ( double)yVal,y->data.F64[2],y->data.F64[0]);1348 (psF64)yVal,y->data.F64[2],y->data.F64[0]); 1349 1349 } 1350 1350 yErr->data.F64[0] = 1.0; … … 1425 1425 XXX: Check for errors in psLib routines that we call. 1426 1426 *****************************************************************************/ 1427 intp_psVectorRobustStats(const psVector* restrict myVector,1428 const psVector* restrict errors,1429 const psVector* restrict maskVector,1430 psU32 maskVal,1431 psStats* stats)1427 psS32 p_psVectorRobustStats(const psVector* restrict myVector, 1428 const psVector* restrict errors, 1429 const psVector* restrict maskVector, 1430 psU32 maskVal, 1431 psStats* stats) 1432 1432 { 1433 1433 psHistogram* robustHistogram = NULL; 1434 1434 psVector* robustHistogramVector = NULL; 1435 floatbinSize = 0.0; // Size of the histogram bins1435 psF32 binSize = 0.0; // Size of the histogram bins 1436 1436 psS32 LQBinNum = -1; // Bin num for lower quartile 1437 1437 psS32 UQBinNum = -1; // Bin num for upper quartile … … 1439 1439 psS32 i = 0; // Loop index variable 1440 1440 psS32 modeBinNum = 0; 1441 floatmodeBinCount = 0.0;1442 floatdL = 0.0;1441 psF32 modeBinCount = 0.0; 1442 psF32 dL = 0.0; 1443 1443 psS32 numBins = 0; 1444 floatmyMean = 0.0;1445 floatmyStdev = 0.0;1446 floatcountFloat = 0.0;1447 floatdiff = 0.0;1448 floatsumSquares = 0.0;1449 floatsumDiffs = 0.0;1444 psF32 myMean = 0.0; 1445 psF32 myStdev = 0.0; 1446 psF32 countFloat = 0.0; 1447 psF32 diff = 0.0; 1448 psF32 sumSquares = 0.0; 1449 psF32 sumDiffs = 0.0; 1450 1450 psVector* cumulativeRobustSums = NULL; 1451 floatsumRobust = 0.0;1452 floatsumN50 = 0.0;1453 floatsumNfit = 0.0;1451 psF32 sumRobust = 0.0; 1452 psF32 sumN50 = 0.0; 1453 psF32 sumNfit = 0.0; 1454 1454 psScalar tmpScalar; 1455 1455 tmpScalar.type.type = PS_TYPE_F32; … … 1460 1460 // that by 10.0; 1461 1461 //XXX: add errors 1462 intrc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats);1462 psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); 1463 1463 if (rc != 0) { 1464 1464 psError(PS_ERR_UNEXPECTED_NULL, … … 1614 1614 for (i=0;i<robustHistogramVector->n;i++) { 1615 1615 myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); 1616 ((psVector *) (myCoords->data[i]))->data.F32[0] = ( float) i;1616 ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i; 1617 1617 y->data.F32[i] = robustHistogramVector->data.F32[i]; 1618 1618 } … … 1770 1770 The histogram structure 1771 1771 *****************************************************************************/ 1772 psHistogram* psHistogramAlloc( float lower, floatupper, psS32 n)1772 psHistogram* psHistogramAlloc(psF32 lower, psF32 upper, psS32 n) 1773 1773 { 1774 1774 PS_INT_CHECK_POSITIVE(n, NULL); … … 1777 1777 psS32 i = 0; // Loop index variable 1778 1778 psHistogram* newHist = NULL; // The new histogram structure 1779 floatbinSize = 0.0; // The histogram bin size1779 psF32 binSize = 0.0; // The histogram bin size 1780 1780 1781 1781 // Allocate memory for the new histogram structure. If there are N … … 1787 1787 1788 1788 // Calculate the bounds for each bin. 1789 binSize = (upper - lower) / ( float)n;1789 binSize = (upper - lower) / (psF32)n; 1790 1790 // XXX: Is the following necessary? It prevents the max data point 1791 1791 // from being in a non-existant bin. 1792 1792 binSize += FLT_EPSILON; 1793 1793 for (i = 0; i < n + 1; i++) { 1794 newHist->bounds->data.F32[i] = lower + (binSize * ( float)i);1794 newHist->bounds->data.F32[i] = lower + (binSize * (psF32)i); 1795 1795 } 1796 1796 … … 1857 1857 } 1858 1858 1859 /***************************************************************************** 1860 UpdateHistogramBins(binNum, out, data, error): This routine is to be used when 1861 updating the histogram in the presence of errors in the input data. We treat 1862 the data point as a boxcar PDF and update a range of points surrounding the 1863 histogram bin which contains the point. The width of that boxcar is defined 1864 as 2.35 * error. Inputs: 1865 binNum: the bin number of the data point in the histogram 1866 out: the histogram structure 1867 data: the data point value 1868 error: the error in that data point 1869 1870 XXX: Must test this. 1871 *****************************************************************************/ 1859 1872 psS32 UpdateHistogramBins(psS32 binNum, 1860 1873 psHistogram* out, … … 1863 1876 { 1864 1877 PS_PTR_CHECK_NULL(out, -1); 1865 PS_INT_CHECK_RANGE(binNum, 0, out->nums->n-1, -2); 1866 /* 1867 psF32 width = 2.35 * error; 1868 psF32 centerBinWidth = out->bounds->data.F32[binNum+1] - out->bounds->data.F32[binNum] 1869 psF32 boxcarCenter = (out->bounds->data.F32[binNum] + out->bounds->data.F32[binNum+1]) / 2.0; 1870 1871 if (width <= centerBinWidth) { 1872 out->nums->data.F32[binNum]+= 1.0; 1873 } else { 1874 out->nums->data.F32[binNum]+= centerBinWidth / width; 1875 // XXX: walk to the left, adding fractional values. 1876 // XXX: walk to the right, adding fractional values. 1877 1878 1879 } 1880 */ 1878 PS_PTR_CHECK_NULL(out->bounds, -1); 1879 PS_PTR_CHECK_NULL(out->nums, -1); 1880 PS_INT_CHECK_RANGE(binNum, 0, ((out->nums->n)-1), -2); 1881 PS_FLOAT_COMPARE(0.0, error, -3); 1882 PS_FLOAT_CHECK_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4); 1883 1884 psF32 boxcarWidth = 2.35 * error; 1885 psF32 boxcarCenter = (out->bounds->data.F32[binNum] + 1886 out->bounds->data.F32[binNum+1]) / 2.0; 1887 psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); 1888 psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); 1889 psS32 bin; 1890 psS32 boxcarLeftBinNum; 1891 psS32 boxcarRightBinNum; 1892 1893 // Determine the left endpoint of the boxcar for the PDF. 1894 for (bin=binNum ; bin >= 0 ; bin--) { 1895 if (out->nums->data.F32[bin] <= boxcarLeft) { 1896 boxcarLeftBinNum = bin; 1897 break; 1898 } 1899 } 1900 1901 // Determine the right endpoint of the boxcar for the PDF. 1902 for (bin=binNum ; bin < out->nums->n ; bin++) { 1903 if (out->nums->data.F32[bin] >= boxcarRight) { 1904 boxcarRightBinNum = bin; 1905 break; 1906 } 1907 } 1908 1909 // 1910 // If the boxcar fits entirely inside this bin, then simply add 1.0 to the 1911 // bin and return. 1912 // 1913 if (boxcarLeftBinNum == boxcarRightBinNum) { 1914 out->nums->data.F32[binNum]+= 1.0; 1915 return(0); 1916 } 1917 1918 // 1919 // If we get here, multiple bins must be updated. We handle the left 1920 // endpoint, and right endpoint differently. 1921 // 1922 out->nums->data.F32[boxcarLeftBinNum]+= 1923 (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth; 1924 1925 // 1926 // Loop through the center bins, if any. 1927 // 1928 for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) { 1929 out->nums->data.F32[bin]+= 1930 (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth; 1931 } 1932 1933 // 1934 // Handle the right endpoint differently. 1935 // 1936 out->nums->data.F32[boxcarRightBinNum]+= 1937 (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth; 1938 1939 // 1940 // Return 0 on success. 1941 // 1881 1942 return(0); 1882 1943 } … … 1906 1967 { 1907 1968 PS_PTR_CHECK_NULL(out, NULL); 1969 PS_VECTOR_CHECK_NULL(out->bounds, NULL); 1908 1970 PS_VECTOR_CHECK_TYPE(out->bounds, PS_TYPE_F32, NULL); 1971 PS_INT_CHECK_NON_NEGATIVE(out->bounds->n, NULL); 1972 PS_VECTOR_CHECK_NULL(out->nums, NULL); 1909 1973 PS_VECTOR_CHECK_TYPE(out->nums, PS_TYPE_F32, NULL); 1910 1974 PS_INT_CHECK_NON_NEGATIVE(out->nums->n, NULL); … … 1920 1984 1921 1985 psS32 i = 0; // Loop index variable 1922 floatbinSize = 0.0; // Histogram bin size1986 psF32 binSize = 0.0; // Histogram bin size 1923 1987 psS32 binNum = 0; // A temporary bin number 1924 1988 psS32 numBins = 0; // The total number of bins 1925 psS32 tmp = 0;1926 1989 psScalar tmpScalar; 1927 1990 tmpScalar.type.type = PS_TYPE_F32; 1928 psVector* inF32; 1929 psS32 mustFreeTmp = 1; 1991 psVector* inF32 = NULL; 1992 psVector* errorsF32 = NULL; 1993 psS32 mustFreeVectorIn = 1; 1994 psS32 mustFreeVectorErrors = 1; 1995 1996 // Convert input and errors vectors to F32 if necessary. 1930 1997 inF32 = p_psConvertToF32((psVector *) in); 1931 1998 if (inF32 == NULL) { 1932 1999 inF32 = (psVector *) in; 1933 mustFreeTmp = 0; 2000 mustFreeVectorIn = 0; 2001 } 2002 errorsF32 = p_psConvertToF32((psVector *) errors); 2003 if (errorsF32 == NULL) { 2004 errorsF32 = (psVector *) errors; 2005 mustFreeVectorErrors = 0; 1934 2006 } 1935 2007 … … 1950 2022 binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; 1951 2023 binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); 1952 if (errors != NULL) {2024 if (errorsF32 != NULL) { 1953 2025 // XXX: Check return codes. 1954 2026 UpdateHistogramBins(binNum, out, 1955 2027 inF32->data.F32[i], 1956 errors ->data.F32[i]);2028 errorsF32->data.F32[i]); 1957 2029 } else { 1958 2030 // XXX: This if-statement really shouldn't be necessary. … … 1969 2041 // correct bin number requires a bit more work. 1970 2042 tmpScalar.data.F32 = inF32->data.F32[i]; 1971 tmp= p_psVectorBinDisect(out->bounds, &tmpScalar);1972 if ( tmp< 0) {2043 binNum = p_psVectorBinDisect(out->bounds, &tmpScalar); 2044 if (binNum < 0) { 1973 2045 psLogMsg(__func__, PS_LOG_WARN, 1974 2046 "WARNING: psVectorHistogram(): element outside histogram bounds.\n"); 1975 2047 } else { 1976 if (errors != NULL) {2048 if (errorsF32 != NULL) { 1977 2049 // XXX: Check return codes. 1978 UpdateHistogramBins( tmp, out,2050 UpdateHistogramBins(binNum, out, 1979 2051 inF32->data.F32[i], 1980 2052 errors->data.F32[i]); 1981 2053 } else { 1982 (out->nums->data.F32[ tmp])+= 1.0;2054 (out->nums->data.F32[binNum])+= 1.0; 1983 2055 } 1984 2056 } … … 1988 2060 } 1989 2061 1990 if (mustFree Tmp== 1) {2062 if (mustFreeVectorIn == 1) { 1991 2063 psFree(inF32); 2064 } 2065 if (mustFreeVectorErrors == 1) { 2066 psFree(errorsF32); 1992 2067 } 1993 2068 return (out); … … 2015 2090 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2016 2091 for (i = 0; i < in->n; i++) { 2017 tmp->data.F32[i] = (float)in->data.S8[i]; 2092 tmp->data.F32[i] = (psF32)in->data.S8[i]; 2093 } 2094 } else if (in->type.type == PS_TYPE_S16) { 2095 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2096 for (i = 0; i < in->n; i++) { 2097 tmp->data.F32[i] = (psF32) in->data.S16[i]; 2098 } 2099 } else if (in->type.type == PS_TYPE_S32) { 2100 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2101 for (i = 0; i < in->n; i++) { 2102 tmp->data.F32[i] = (psF32)in->data.S32[i]; 2103 } 2104 } else if (in->type.type == PS_TYPE_S64) { 2105 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2106 for (i = 0; i < in->n; i++) { 2107 tmp->data.F32[i] = (psF32)in->data.S64[i]; 2108 } 2109 } else if (in->type.type == PS_TYPE_U8) { 2110 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2111 for (i = 0; i < in->n; i++) { 2112 tmp->data.F32[i] = (psF32)in->data.U8[i]; 2018 2113 } 2019 2114 } else if (in->type.type == PS_TYPE_U16) { 2020 2115 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2021 2116 for (i = 0; i < in->n; i++) { 2022 tmp->data.F32[i] = ( float)in->data.U16[i];2023 } 2024 } else if (in->type.type == PS_TYPE_U 8) {2117 tmp->data.F32[i] = (psF32)in->data.U16[i]; 2118 } 2119 } else if (in->type.type == PS_TYPE_U32) { 2025 2120 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2026 2121 for (i = 0; i < in->n; i++) { 2027 tmp->data.F32[i] = (float)in->data.U8[i]; 2122 tmp->data.F32[i] = (psF32)in->data.U32[i]; 2123 } 2124 } else if (in->type.type == PS_TYPE_U64) { 2125 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2126 for (i = 0; i < in->n; i++) { 2127 tmp->data.F32[i] = (psF32)in->data.U64[i]; 2028 2128 } 2029 2129 } else if (in->type.type == PS_TYPE_F64) { 2030 2130 tmp = psVectorAlloc(in->n, PS_TYPE_F32); 2031 2131 for (i = 0; i < in->n; i++) { 2032 tmp->data.F32[i] = ( float)in->data.F64[i];2132 tmp->data.F32[i] = (psF32)in->data.F64[i]; 2033 2133 } 2034 2134 } else if (in->type.type == PS_TYPE_F32) { 2035 2135 // do nothing 2036 2136 } else { 2037 char* strType;2038 PS_TYPE_NAME(strType, in->type.type);2137 psS8* strType; 2138 PS_TYPE_NAME(strType, in->type.type); 2039 2139 psError(PS_ERR_BAD_PARAMETER_TYPE, true, 2040 2140 PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED, … … 2074 2174 } 2075 2175 2076 psVector* inF32 ;2077 psVector* errorsF32 ;2176 psVector* inF32 = NULL; 2177 psVector* errorsF32 = NULL; 2078 2178 psS32 mustFreeVectorIn = 1; 2079 2179 psS32 mustFreeVectorErrors = 1; … … 2178 2278 return (stats); 2179 2279 } 2280 -
trunk/psLib/src/math/psStats.h
r2778 r2788 10 10 * @author George Gusciora, MHPCC 11 11 * 12 * @version $Revision: 1.3 6$ $Name: not supported by cvs2svn $13 * @date $Date: 2004-12-2 1 20:42:07$12 * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2004-12-22 05:09:32 $ 14 14 * 15 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 56 56 typedef struct 57 57 { 58 doublesampleMean; ///< formal mean of sample59 doublesampleMedian; ///< formal median of sample60 doublesampleStdev; ///< standard deviation of sample61 doublesampleUQ; ///< upper quartile of sample62 doublesampleLQ; ///< lower quartile of sample63 doublerobustMean; ///< robust mean of array64 doublerobustMedian; ///< robust median of array65 doublerobustMode; ///< Robust mode of array66 doublerobustStdev; ///< robust standard deviation of array67 doublerobustUQ; ///< robust upper quartile68 doublerobustLQ; ///< robust lower quartile58 psF64 sampleMean; ///< formal mean of sample 59 psF64 sampleMedian; ///< formal median of sample 60 psF64 sampleStdev; ///< standard deviation of sample 61 psF64 sampleUQ; ///< upper quartile of sample 62 psF64 sampleLQ; ///< lower quartile of sample 63 psF64 robustMean; ///< robust mean of array 64 psF64 robustMedian; ///< robust median of array 65 psF64 robustMode; ///< Robust mode of array 66 psF64 robustStdev; ///< robust standard deviation of array 67 psF64 robustUQ; ///< robust upper quartile 68 psF64 robustLQ; ///< robust lower quartile 69 69 psS32 robustN50; ///< 70 70 psS32 robustNfit; ///< 71 doubleclippedMean; ///< Nsigma clipped mean72 doubleclippedStdev; ///< standard deviation after clipping71 psF64 clippedMean; ///< Nsigma clipped mean 72 psF64 clippedStdev; ///< standard deviation after clipping 73 73 psS32 clippedNvalues; ///< ??? 74 doubleclipSigma; ///< Nsigma used for clipping; user input74 psF64 clipSigma; ///< Nsigma used for clipping; user input 75 75 psS32 clipIter; ///< Number of clipping iterations; user input 76 doublemin; ///< minimum data value in array77 doublemax; ///< maximum data value in array78 doublebinsize; ///<76 psF64 min; ///< minimum data value in array 77 psF64 max; ///< maximum data value in array 78 psF64 binsize; ///< 79 79 psStatsOptions options; ///< bitmask of calculated values 80 80 } … … 132 132 */ 133 133 psHistogram* psHistogramAlloc( 134 floatlower, ///< Lower limit for the bins135 floatupper, ///< Upper limit for the bins134 psF32 lower, ///< Lower limit for the bins 135 psF32 upper, ///< Upper limit for the bins 136 136 psS32 n ///< Number of bins 137 137 ); … … 172 172 ///< the statistic struct to operate on 173 173 174 double*value174 psF64 *value 175 175 ///< if return is true, this is set to the specified statistic value by stats->options 176 176 );
Note:
See TracChangeset
for help on using the changeset viewer.
