Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 2782) +++ trunk/psLib/src/math/psStats.c (revision 2788) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.107 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-22 00:54:28 $ + * @version $Revision: 1.108 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-22 05:09:32 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -70,5 +70,5 @@ /*****************************************************************************/ -psBool p_psGetStatValue(const psStats* stats, double *value) +psBool p_psGetStatValue(const psStats* stats, psF64 *value) { @@ -145,13 +145,13 @@ this routine sets stats->sampleMean to NAN. *****************************************************************************/ -int p_psVectorSampleMean(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorSampleMean(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float mean = 0.0; // The mean + psF32 mean = 0.0; // The mean psS32 count = 0; // # of points in this mean @@ -171,5 +171,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -185,5 +185,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -199,5 +199,5 @@ } if (count != 0) { - mean /= (float)count; + mean /= (psF32)count; } else { mean = NAN; @@ -207,5 +207,5 @@ mean += myVector->data.F32[i]; } - mean /= (float)myVector->n; + mean /= (psF32)myVector->n; } } @@ -283,11 +283,11 @@ this routine sets stats->max to NAN. *****************************************************************************/ -int p_psVectorMax(const psVector* restrict myVector, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorMax(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float max = -PS_MAX_F32; // The calculated maximum + psF32 max = -PS_MAX_F32; // The calculated maximum psS32 empty = true; // Does this vector have valid elements? @@ -348,11 +348,11 @@ this routine sets stats->min to NAN. *****************************************************************************/ -int p_psVectorMin(const psVector* restrict myVector, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorMin(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable - float min = PS_MAX_F32; // The calculated maximum + psF32 min = PS_MAX_F32; // The calculated maximum psS32 empty = true; // Does this vector have valid elements? @@ -612,5 +612,5 @@ *****************************************************************************/ psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, - float sigma) + psF32 sigma) { PS_PTR_CHECK_NULL(robustHistogram, NULL); @@ -619,6 +619,6 @@ psS32 i = 0; // Loop index variable psS32 j = 0; // Loop index variable - float iMid; - float jMid; + psF32 iMid; + psF32 jMid; psS32 numBins = robustHistogram->nums->n; psS32 numBounds = robustHistogram->bounds->n; @@ -626,6 +626,6 @@ psS32 jMin = 0; psS32 jMax = 0; - float firstBound = robustHistogram->bounds->data.F32[0]; - float lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + psF32 firstBound = robustHistogram->bounds->data.F32[0]; + psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1]; psScalar x; @@ -795,9 +795,9 @@ psS32 i = 0; // Loop index variable psS32 countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable + psF32 countFloat = 0.0; // # of data points being used + psF32 mean = 0.0; // The mean + psF32 diff = 0.0; // Used in calculating stdev + psF32 sumSquares = 0.0; // temporary variable + psF32 sumDiffs = 0.0; // temporary variable // This procedure requires the mean. If it has not been already @@ -865,5 +865,5 @@ psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", countInt); } else { - countFloat = (float)countInt; + countFloat = (psF32)countInt; stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); } @@ -889,9 +889,9 @@ psS32 i = 0; // Loop index variable psS32 countInt = 0; // # of data points being used - float countFloat = 0.0; // # of data points being used - float mean = 0.0; // The mean - float diff = 0.0; // Used in calculating stdev - float sumSquares = 0.0; // temporary variable - float sumDiffs = 0.0; // temporary variable + psF32 countFloat = 0.0; // # of data points being used + psF32 mean = 0.0; // The mean + psF32 diff = 0.0; // Used in calculating stdev + psF32 sumSquares = 0.0; // temporary variable + psF32 sumDiffs = 0.0; // temporary variable // psF32 sum1; // psF32 sum2; @@ -982,5 +982,5 @@ stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor)); } else { - countFloat = (float)countInt; + countFloat = (psF32)countInt; stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); @@ -1003,16 +1003,16 @@ -2: warning *****************************************************************************/ -int p_psVectorClippedStats(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorClippedStats(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable psS32 j = 0; // Loop index variable - float clippedMean = 0.0; // self-explanatory - float clippedStdev = 0.0; // self-explanatory - float oldStanMean = 0.0; // Temporary variable - float oldStanStdev = 0.0; // Temporary variable + psF32 clippedMean = 0.0; // self-explanatory + psF32 clippedStdev = 0.0; // self-explanatory + psF32 oldStanMean = 0.0; // Temporary variable + psF32 oldStanStdev = 0.0; // Temporary variable psVector* tmpMask = NULL; // Temporary vector @@ -1240,8 +1240,8 @@ XXX: Create a 2nd-order polynomial version and solve for X analytically. *****************************************************************************/ -float p_ps1DPolyMedian(psPolynomial1D* myPoly, - float rangeLow, - float rangeHigh, - float getThisValue) +psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly, + psF32 rangeLow, + psF32 rangeHigh, + psF32 getThisValue) { PS_POLY_CHECK_NULL(myPoly, NAN); @@ -1250,9 +1250,9 @@ // falls within the range of y-values of the polynomial "myPoly" in the // specified x-range (rangeLow:rangeHigh). - float fLo = psPolynomial1DEval( + psF32 fLo = psPolynomial1DEval( myPoly, rangeLow ); - float fHi = psPolynomial1DEval( + psF32 fHi = psPolynomial1DEval( myPoly, rangeHigh @@ -1266,7 +1266,7 @@ psS32 numIterations = 0; - float midpoint = 0.0; - float oldMidpoint = 1.0; - float f = 0.0; + psF32 midpoint = 0.0; + psF32 oldMidpoint = 1.0; + psF32 f = 0.0; while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) { @@ -1309,8 +1309,8 @@ tests to ensure that binNum is within acceptable ranges for both vectors. *****************************************************************************/ -float fitQuadraticSearchForYThenReturnX(psVector *xVec, +psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, psVector *yVec, psS32 binNum, - float yVal) + psF32 yVal) { PS_VECTOR_CHECK_NULL(xVec, NAN); @@ -1331,11 +1331,11 @@ psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - float tmpFloat; + psF32 tmpFloat; if ((binNum > 0) && (binNum < (yVec->n - 2))) { // The general case. We have all three points. - x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); - x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); - x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); + x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); y->data.F64[0] = yVec->data.F32[binNum - 1]; y->data.F64[1] = yVec->data.F32[binNum]; @@ -1346,5 +1346,5 @@ psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, - (double)yVal,y->data.F64[2],y->data.F64[0]); + (psF64)yVal,y->data.F64[2],y->data.F64[0]); } yErr->data.F64[0] = 1.0; @@ -1425,13 +1425,13 @@ XXX: Check for errors in psLib routines that we call. *****************************************************************************/ -int p_psVectorRobustStats(const psVector* restrict myVector, - const psVector* restrict errors, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +psS32 p_psVectorRobustStats(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psHistogram* robustHistogram = NULL; psVector* robustHistogramVector = NULL; - float binSize = 0.0; // Size of the histogram bins + psF32 binSize = 0.0; // Size of the histogram bins psS32 LQBinNum = -1; // Bin num for lower quartile psS32 UQBinNum = -1; // Bin num for upper quartile @@ -1439,17 +1439,17 @@ psS32 i = 0; // Loop index variable psS32 modeBinNum = 0; - float modeBinCount = 0.0; - float dL = 0.0; + psF32 modeBinCount = 0.0; + psF32 dL = 0.0; psS32 numBins = 0; - float myMean = 0.0; - float myStdev = 0.0; - float countFloat = 0.0; - float diff = 0.0; - float sumSquares = 0.0; - float sumDiffs = 0.0; + psF32 myMean = 0.0; + psF32 myStdev = 0.0; + psF32 countFloat = 0.0; + psF32 diff = 0.0; + psF32 sumSquares = 0.0; + psF32 sumDiffs = 0.0; psVector* cumulativeRobustSums = NULL; - float sumRobust = 0.0; - float sumN50 = 0.0; - float sumNfit = 0.0; + psF32 sumRobust = 0.0; + psF32 sumN50 = 0.0; + psF32 sumNfit = 0.0; psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; @@ -1460,5 +1460,5 @@ // that by 10.0; //XXX: add errors - int rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); + psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); if (rc != 0) { psError(PS_ERR_UNEXPECTED_NULL, @@ -1614,5 +1614,5 @@ for (i=0;in;i++) { myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); - ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) i; + ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i; y->data.F32[i] = robustHistogramVector->data.F32[i]; } @@ -1770,5 +1770,5 @@ The histogram structure *****************************************************************************/ -psHistogram* psHistogramAlloc(float lower, float upper, psS32 n) +psHistogram* psHistogramAlloc(psF32 lower, psF32 upper, psS32 n) { PS_INT_CHECK_POSITIVE(n, NULL); @@ -1777,5 +1777,5 @@ psS32 i = 0; // Loop index variable psHistogram* newHist = NULL; // The new histogram structure - float binSize = 0.0; // The histogram bin size + psF32 binSize = 0.0; // The histogram bin size // Allocate memory for the new histogram structure. If there are N @@ -1787,10 +1787,10 @@ // Calculate the bounds for each bin. - binSize = (upper - lower) / (float)n; + binSize = (upper - lower) / (psF32)n; // XXX: Is the following necessary? It prevents the max data point // from being in a non-existant bin. binSize += FLT_EPSILON; for (i = 0; i < n + 1; i++) { - newHist->bounds->data.F32[i] = lower + (binSize * (float)i); + newHist->bounds->data.F32[i] = lower + (binSize * (psF32)i); } @@ -1857,4 +1857,17 @@ } +/***************************************************************************** +UpdateHistogramBins(binNum, out, data, error): This routine is to be used when +updating the histogram in the presence of errors in the input data. We treat +the data point as a boxcar PDF and update a range of points surrounding the +histogram bin which contains the point. The width of that boxcar is defined +as 2.35 * error. Inputs: + binNum: the bin number of the data point in the histogram + out: the histogram structure + data: the data point value + error: the error in that data point + +XXX: Must test this. + *****************************************************************************/ psS32 UpdateHistogramBins(psS32 binNum, psHistogram* out, @@ -1863,20 +1876,68 @@ { PS_PTR_CHECK_NULL(out, -1); - PS_INT_CHECK_RANGE(binNum, 0, out->nums->n-1, -2); - /* - psF32 width = 2.35 * error; - psF32 centerBinWidth = out->bounds->data.F32[binNum+1] - out->bounds->data.F32[binNum] - psF32 boxcarCenter = (out->bounds->data.F32[binNum] + out->bounds->data.F32[binNum+1]) / 2.0; - - if (width <= centerBinWidth) { - out->nums->data.F32[binNum]+= 1.0; - } else { - out->nums->data.F32[binNum]+= centerBinWidth / width; - // XXX: walk to the left, adding fractional values. - // XXX: walk to the right, adding fractional values. - - - } - */ + PS_PTR_CHECK_NULL(out->bounds, -1); + PS_PTR_CHECK_NULL(out->nums, -1); + PS_INT_CHECK_RANGE(binNum, 0, ((out->nums->n)-1), -2); + PS_FLOAT_COMPARE(0.0, error, -3); + PS_FLOAT_CHECK_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4); + + psF32 boxcarWidth = 2.35 * error; + psF32 boxcarCenter = (out->bounds->data.F32[binNum] + + out->bounds->data.F32[binNum+1]) / 2.0; + psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); + psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); + psS32 bin; + psS32 boxcarLeftBinNum; + psS32 boxcarRightBinNum; + + // Determine the left endpoint of the boxcar for the PDF. + for (bin=binNum ; bin >= 0 ; bin--) { + if (out->nums->data.F32[bin] <= boxcarLeft) { + boxcarLeftBinNum = bin; + break; + } + } + + // Determine the right endpoint of the boxcar for the PDF. + for (bin=binNum ; bin < out->nums->n ; bin++) { + if (out->nums->data.F32[bin] >= boxcarRight) { + boxcarRightBinNum = bin; + break; + } + } + + // + // If the boxcar fits entirely inside this bin, then simply add 1.0 to the + // bin and return. + // + if (boxcarLeftBinNum == boxcarRightBinNum) { + out->nums->data.F32[binNum]+= 1.0; + return(0); + } + + // + // If we get here, multiple bins must be updated. We handle the left + // endpoint, and right endpoint differently. + // + out->nums->data.F32[boxcarLeftBinNum]+= + (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth; + + // + // Loop through the center bins, if any. + // + for (bin = boxcarLeftBinNum + 1 ; bin < (boxcarRightBinNum - 1) ; bin++) { + out->nums->data.F32[bin]+= + (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth; + } + + // + // Handle the right endpoint differently. + // + out->nums->data.F32[boxcarRightBinNum]+= + (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth; + + // + // Return 0 on success. + // return(0); } @@ -1906,5 +1967,8 @@ { PS_PTR_CHECK_NULL(out, NULL); + PS_VECTOR_CHECK_NULL(out->bounds, NULL); PS_VECTOR_CHECK_TYPE(out->bounds, PS_TYPE_F32, NULL); + PS_INT_CHECK_NON_NEGATIVE(out->bounds->n, NULL); + PS_VECTOR_CHECK_NULL(out->nums, NULL); PS_VECTOR_CHECK_TYPE(out->nums, PS_TYPE_F32, NULL); PS_INT_CHECK_NON_NEGATIVE(out->nums->n, NULL); @@ -1920,16 +1984,24 @@ psS32 i = 0; // Loop index variable - float binSize = 0.0; // Histogram bin size + psF32 binSize = 0.0; // Histogram bin size psS32 binNum = 0; // A temporary bin number psS32 numBins = 0; // The total number of bins - psS32 tmp = 0; psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; - psVector* inF32; - psS32 mustFreeTmp = 1; + psVector* inF32 = NULL; + psVector* errorsF32 = NULL; + psS32 mustFreeVectorIn = 1; + psS32 mustFreeVectorErrors = 1; + + // Convert input and errors vectors to F32 if necessary. inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { inF32 = (psVector *) in; - mustFreeTmp = 0; + mustFreeVectorIn = 0; + } + errorsF32 = p_psConvertToF32((psVector *) errors); + if (errorsF32 == NULL) { + errorsF32 = (psVector *) errors; + mustFreeVectorErrors = 0; } @@ -1950,9 +2022,9 @@ binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); - if (errors != NULL) { + if (errorsF32 != NULL) { // XXX: Check return codes. UpdateHistogramBins(binNum, out, inF32->data.F32[i], - errors->data.F32[i]); + errorsF32->data.F32[i]); } else { // XXX: This if-statement really shouldn't be necessary. @@ -1969,16 +2041,16 @@ // correct bin number requires a bit more work. tmpScalar.data.F32 = inF32->data.F32[i]; - tmp = p_psVectorBinDisect(out->bounds, &tmpScalar); - if (tmp < 0) { + binNum = p_psVectorBinDisect(out->bounds, &tmpScalar); + if (binNum < 0) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: psVectorHistogram(): element outside histogram bounds.\n"); } else { - if (errors != NULL) { + if (errorsF32 != NULL) { // XXX: Check return codes. - UpdateHistogramBins(tmp, out, + UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i]); } else { - (out->nums->data.F32[tmp])+= 1.0; + (out->nums->data.F32[binNum])+= 1.0; } } @@ -1988,6 +2060,9 @@ } - if (mustFreeTmp == 1) { + if (mustFreeVectorIn == 1) { psFree(inF32); + } + if (mustFreeVectorErrors == 1) { + psFree(errorsF32); } return (out); @@ -2015,26 +2090,51 @@ tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.S8[i]; + tmp->data.F32[i] = (psF32)in->data.S8[i]; + } + } else if (in->type.type == PS_TYPE_S16) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32) in->data.S16[i]; + } + } else if (in->type.type == PS_TYPE_S32) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.S32[i]; + } + } else if (in->type.type == PS_TYPE_S64) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.S64[i]; + } + } else if (in->type.type == PS_TYPE_U8) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.U8[i]; } } else if (in->type.type == PS_TYPE_U16) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.U16[i]; - } - } else if (in->type.type == PS_TYPE_U8) { + tmp->data.F32[i] = (psF32)in->data.U16[i]; + } + } else if (in->type.type == PS_TYPE_U32) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.U8[i]; + tmp->data.F32[i] = (psF32)in->data.U32[i]; + } + } else if (in->type.type == PS_TYPE_U64) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (psF32)in->data.U64[i]; } } else if (in->type.type == PS_TYPE_F64) { tmp = psVectorAlloc(in->n, PS_TYPE_F32); for (i = 0; i < in->n; i++) { - tmp->data.F32[i] = (float)in->data.F64[i]; + tmp->data.F32[i] = (psF32)in->data.F64[i]; } } else if (in->type.type == PS_TYPE_F32) { // do nothing } else { - char* strType; - PS_TYPE_NAME(strType,in->type.type); + psS8* strType; + PS_TYPE_NAME(strType, in->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, true, PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED, @@ -2074,6 +2174,6 @@ } - psVector* inF32; - psVector* errorsF32; + psVector* inF32 = NULL; + psVector* errorsF32 = NULL; psS32 mustFreeVectorIn = 1; psS32 mustFreeVectorErrors = 1; @@ -2178,2 +2278,3 @@ return (stats); } +