Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 1292) +++ trunk/psLib/src/math/psStats.c (revision 1293) @@ -131,4 +131,8 @@ mean of the input vector. Inputs + myVector + maskVector + maskVal + stats Returns NULL @@ -404,28 +408,30 @@ float rangeMin = 0.0; // Exclude data below this float rangeMax = 0.0; // Exclude date above this - psStats *stats2 = NULL; // Temporary stats structure // Determine if the number of data points exceed a threshold which will // cause to generate robust stats, as opposed to exact stats. - - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleMedian = stats2->robustMedian; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } + // XXX: This code is no longer used. We now calculate the median exactly + // regardless of the vector size. + /* + if (myVector->n > stats->sampleLimit) { + psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); + + // Calculate the robust quartiles. + stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); + + // Store the robust quartiles into the sample quartile members. + stats->sampleMedian = stats2->robustMedian; + + // Free temporary data buffers. + psFree(stats2); + + // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. + stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; + + return; + } + */ // Determine how many data points fit inside this min/max range @@ -483,4 +489,5 @@ // Calculate the median exactly. + // XXX: Is this the correct action? if (0 == (nValues % 2)) { stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] + @@ -498,4 +505,12 @@ This routine smoothes the data in the input robustHistogram with a Gaussian of width sigma. + + XXX: Must consult with IfA on the proper values for the Gaussian + smoothing. Currently, the Gaussian coefficients are such that only + the center coefficient is non-zero. This is because "e" is being + raised to a very large negative number for all points except the + center point. + + XXX: use a static variable for gaussianCoefs[] and compute them once. *****************************************************************************/ psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram, @@ -541,4 +556,5 @@ } + // Perform the actual smoothing. for(i=0;inums->n;i++) { smooth->data.F32[i] = 0.0; @@ -576,26 +592,4 @@ float rangeMin = 0.0; // Exclude data below this float rangeMax = 0.0; // Exclude date above this - - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - if (myVector->n > stats->sampleLimit) { - psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented."); - psStats *stats2 = NULL; - // Calculate the robust quartiles. - stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE); - p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - - // Store the robust quartiles into the sample quartile members. - stats->sampleUQ = stats2->robustUQ; - stats->sampleLQ = stats2->robustLQ; - - // Free temporary data buffers. - psFree(stats2); - - // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. - stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE; - - return; - } // Determine how many data points fit inside this min/max range @@ -861,4 +855,6 @@ /***************************************************************************** +p_psNormalizeVector(myData): this is a private function which normalizes the +elements of a vector to a range between 0.0 and 1.0. *****************************************************************************/ void p_psNormalizeVector(psVector *myData) @@ -885,4 +881,9 @@ +/***************************************************************************** +p_psGaussian(myData, myParams): an internal function, used by robust stats, +intended to compute the Gaussian with a specified sigma/mean, at the +specified data point. + *****************************************************************************/ float p_psGaussian(const psVector *restrict myData, const psVector *restrict myParams) @@ -894,8 +895,12 @@ tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); - printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); + // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); return(tmp); } +/***************************************************************************** +p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which +calculates the specified partial derivative of the above Gaussian function. + *****************************************************************************/ float p_psGaussianDeriv(const psVector *restrict myData, const psVector *restrict myParams, @@ -922,4 +927,9 @@ +/***************************************************************************** +p_psQuadratic(myData, myParams): an internal function, used by robust stats, +intended to compute a quadratic, with the specified parameters, at the +specified data point. + *****************************************************************************/ float p_psQuadratic(const psVector *restrict myParams, const psVector *restrict myCoords) @@ -935,4 +945,8 @@ } +/***************************************************************************** +p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which +calculates the specified partial derivative of the above quadratic function. + *****************************************************************************/ float p_psQuadraticDeriv(const psVector *restrict myParams, const psVector *restrict myCoords, @@ -959,6 +973,6 @@ [rangeLow, rangeHigh]. It determines the x-value of that polynomial such that f(x) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polnomial is monotonically increasing within that -range. +range and assumes that the polynomial is monotonically increasing or +decreasing within that range. *****************************************************************************/ float p_ps1DPolyMedian(psPolynomial1D *myPoly, @@ -998,49 +1012,58 @@ /****************************************************************************** -p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes -as input a 1-D polynomial of arbitrary order (though we are using 2nd-order -polynomials here) and a range of x-values for which it is defined: -[rangeLow, rangeHigh]. It determines the x-value of that polynomial such -that f(x) == midpoint. This functions uses a binary-search algorithm on the -range and assumes that the polnomial is monotonically increasing within that -range. - *****************************************************************************/ -float p_psFitQuadratic(psHistogram *robustHistogram, +p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given +the specified histogram, with the specified pre-calculated cumulativeSums, +this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum", +and then finds the point x such that f(x) == fitFloat. + +XXX: This function is currently not being used. + *****************************************************************************/ +float p_psFitQuadratic(psHistogram *histogram, + psVector *cumulativeSums, int binNum, float fitFloat) { - /* - if ((binNum > 0) && - (binNum < (robustHistogram->nums->n+1))) { - x->data.F64[0] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum-1] + - robustHistogram->bounds->data.F32[binNum]); - x->data.F64[1] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum] + - robustHistogram->bounds->data.F32[binNum+1]); - x->data.F64[2] = (double) 0.5 * - (robustHistogram->bounds->data.F32[binNum+1] + - robustHistogram->bounds->data.F32[binNum+2]); - - y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1]; - y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum]; - y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1]; - - if (!((y->data.F64[0] <= fitFloat) && - (fitFloat <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); - return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat); - } else { - return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] + - robustHistogram->bounds->data.F32[binNum])); - } - */ + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + psPolynomial1D *myPoly = psPolynomial1DAlloc(2); + + if ((binNum > 0) && + (binNum < (histogram->nums->n+1))) { + x->data.F64[0] = (double) 0.5 * + (histogram->bounds->data.F32[binNum-1] + + histogram->bounds->data.F32[binNum]); + x->data.F64[1] = (double) 0.5 * + (histogram->bounds->data.F32[binNum] + + histogram->bounds->data.F32[binNum+1]); + x->data.F64[2] = (double) 0.5 * + (histogram->bounds->data.F32[binNum+1] + + histogram->bounds->data.F32[binNum+2]); + + y->data.F64[0] = cumulativeSums->data.F32[binNum-1]; + y->data.F64[1] = cumulativeSums->data.F32[binNum]; + y->data.F64[2] = cumulativeSums->data.F32[binNum+1]; + + if (!((y->data.F64[0] <= fitFloat) && + (fitFloat <= y->data.F64[2]))) { + psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); + } + + yErr->data.F64[0] = 1.0; + yErr->data.F64[1] = 1.0; + yErr->data.F64[2] = 1.0; + + myPoly = psGetArrayPolynomial(myPoly, x, y, yErr); + return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], + fitFloat)); + } else { + return(0.5 * (histogram->bounds->data.F32[binNum+1] + + histogram->bounds->data.F32[binNum])); + } + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myPoly); return(0.0); } @@ -1049,25 +1072,24 @@ p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure calculates a variety of robust stat measures: -PS_STAT_ROBUST_MEAN -PS_STAT_ROBUST_MEDIAN -PS_STAT_ROBUST_MODE -PS_STAT_ROBUST_STDEV -PS_STAT_ROBUST_QUARTILE + PS_STAT_ROBUST_MEAN + PS_STAT_ROBUST_MEDIAN + PS_STAT_ROBUST_MODE + PS_STAT_ROBUST_STDEV + PS_STAT_ROBUST_QUARTILE I have included all that computation in a single function, as opposed to -breaking it across several functions for one primary reason: -they all +breaking it across several functions for one primary reason: they all require the same basic initial processing steps (calculate the histogram, -etc.) -however there is no currently defined means, in the SDRS, to -specify this computation as a separate step. If the above robust stat -measures were calcualted in separate functiosn, then much of the initial -processing would be duplicated. +etc.) however there is no currently defined means, in the SDRS, to specify +this computation as a separate step. If the above robust stat measures were +calcualted in separate functiosn, then much of the initial processing would +be duplicated. + Inputs -myVector -maskVector -maskVal -stats + myVector + maskVector + maskVal + stats Returns -NULL + NULL *****************************************************************************/ void p_psVectorRobustStats(const psVector *restrict myVector, @@ -1224,4 +1246,5 @@ sumN50+= (float) robustHistogram->nums->data.S32[i]; } + // XXX: is dL defined as the value at the LQ/UQ, or the bin number? dL = (UQBinNum - LQBinNum) / 4; @@ -1244,5 +1267,5 @@ robustHistogramVector->data.F32[i]; cumulativeMedian+= robustHistogramVector->data.F32[i]; - sumNfit+= (float) robustHistogramVector->data.S32[i]; + sumNfit+= (float) robustHistogram->nums->data.S32[i]; } } @@ -1780,4 +1803,12 @@ } + // XXX: Should we abort if (stats->min == stats->max)? + if (stats->min >= stats->max) { + psAbort(__func__, "psVectorStats() called with range: %f to %f\n", + stats->min, stats->max); + } + + + PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32); if (mask != NULL) {