Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 6214) +++ /trunk/psLib/src/math/psStats.c (revision 6215) @@ -14,6 +14,14 @@ * stats->binsize * - * @version $Revision: 1.162 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-26 21:10:22 $ + * XXX: Must do + * nSubsample points + * use ->min and ->max (PS_STAT_USE_RANGE) + * use ->binsize (PS_STAT_USE_BINSIZE) + * + * + * + * + * @version $Revision: 1.163 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-26 23:49:11 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -97,16 +105,6 @@ return true; - case PS_STAT_ROBUST_MEAN: - *value = stats->robustMean; - psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return true; - case PS_STAT_ROBUST_MEDIAN: *value = stats->robustMedian; - psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); - return true; - - case PS_STAT_ROBUST_MODE: - *value = stats->robustMode; psTrace(__func__, 4, "---- %s(true) end ----\n", __func__); return true; @@ -449,4 +447,6 @@ unmasked element within the specified min/max range). Otherwise, return "false". + +XXX: Can you use psVectorCountPixelMask here? *****************************************************************************/ bool p_psVectorCheckNonEmpty(const psVector* myVector, @@ -503,4 +503,6 @@ number of non-masked pixels in the vector that fall within the min/max range, if specified. + +XXX: Can you use psVectorCountPixelMask here? *****************************************************************************/ psS32 p_psVectorNValues(const psVector* myVector, @@ -766,5 +768,5 @@ PS_VECTOR_PRINT_F32(smooth); } - psTrace(__func__, 4, "---- %s(psVector) end ----\n", __func__); + psTrace(__func__, 4, "---- %s() end ----\n", __func__); return(smooth); } @@ -858,5 +860,5 @@ /****************************************************************************** -p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the +p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the stdev of the input vector. Inputs @@ -868,4 +870,5 @@ NULL +XXX: remove this *****************************************************************************/ void p_psVectorSampleStdevOLD(const psVector* myVector, @@ -1126,5 +1129,5 @@ p_psMemSetPersistent(statsTmp, true); } else { - // EAM : initialize structure if already allocated + // Initialize structure if already allocated statsTmp->sampleMean = NAN; statsTmp->sampleMedian = NAN; @@ -1132,20 +1135,21 @@ statsTmp->sampleUQ = NAN; statsTmp->sampleLQ = NAN; - statsTmp->robustMean = NAN; statsTmp->robustMedian = NAN; - statsTmp->robustMode = NAN; statsTmp->robustStdev = NAN; statsTmp->robustUQ = NAN; statsTmp->robustLQ = NAN; - statsTmp->robustN50 = -1; // XXX: This is never used - statsTmp->robustNfit = -1; + statsTmp->robustN50 = -1; + statsTmp->fittedMean = NAN; + statsTmp->fittedStdev = NAN; + statsTmp->fittedNfit = -1; statsTmp->clippedMean = NAN; statsTmp->clippedStdev = NAN; - statsTmp->clippedNvalues = -1; // XXX: This is never used + statsTmp->clippedNvalues = -1; statsTmp->clipSigma = 3.0; statsTmp->clipIter = 3; statsTmp->min = NAN; statsTmp->max = NAN; - statsTmp->binsize = NAN; // XXX: This is never used + statsTmp->binsize = NAN; + statsTmp->nSubsample = 100000; } @@ -1688,27 +1692,34 @@ } } - psBool iterate = true; + psS32 iterate = 1; psF32 sigma; psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); - while (iterate) { + while (iterate > 0) { psTrace(__func__, 6, "Iterating on Bin size.\n"); - // - // Determine the bin size of the robust histogram. This is done - // by computing the total range of data values and dividing by 1000.0. - // - rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax); - rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax); - if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(tmpStatsMinMax); - psFree(tmpMaskVec); - psFree(tmpScalar); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max); - psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; - psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize); + psF32 binSize = 0.0; + if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) { + // Set initial bin size to the specified value. + binSize = stats->binsize; + psTrace(__func__, 6, "Setting initial robust bin size to %.2f\n", binSize); + } else { + // Determine the bin size of the robust histogram. This is done + // by computing the total range of data values and dividing by 1000.0. + + rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax); + if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(tmpStatsMinMax); + psFree(tmpMaskVec); + psFree(tmpScalar); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max); + binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; + psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize); + } + psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize); // @@ -1726,6 +1737,5 @@ } // XXX: Set these to the number of unmasked data points? - stats->robustNfit = 0.0; - stats->robustN50 = 0.0; + stats->robustN50 = 0; psFree(tmpStatsMinMax); psFree(tmpMaskVec); @@ -1900,4 +1910,5 @@ sigma = (binHiF32 - binLoF32) / 2.0; psTrace(__func__, 6, "The current sigma is %f.\n", sigma); + stats->robustStdev = sigma; // @@ -1921,8 +1932,8 @@ psFree(robustHistogram); psFree(cumulativeRobustHistogram); - + iterate++; } else { psTrace(__func__, 6, "*************: No more iteration. sigma is %f\n", sigma); - iterate = false; + iterate = 0; } } @@ -1981,11 +1992,7 @@ psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32); - // - // - // New algorithm for calculated mean and stdev. - // - // psS32 N50 = 0; for (psS32 i = 0 ; i < myVector->n ; i++) { + // XXX: use maskVal here? if ((0 == tmpMaskVec->data.U8[i]) && (binLo25F32 <= myVector->data.F32[i]) && @@ -1997,163 +2004,199 @@ psTrace(__func__, 6, "The robustN50 is %d.\n", N50); - psF32 dN = (psF32) (0.17 * N50); - if (dN < 1.0) { - dN = 1.0; - } else if (dN > 4.0) { - dN = 4.0; - } - psF32 newBinSize = sigma / dN; - - rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); - rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); - if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpScalar); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize); - psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max); - psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize); - psTrace(__func__, 6, "The numBins is %d\n", numBins); - - psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); - newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1); - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(newHistogram->nums); - } - - // - // Smooth the resulting histogram with a Gaussian with sigma_s = 1 - // bin. - // - psF32 sigma_s = newBinSize; - - psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s); - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(newHistogramSmoothed); - } - - // - // Find the bin with the peak value in the range 2 sigma of the - // robust histogram median. - // - - psS32 binMin = 0; - psS32 binMax = 0; - tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma); - if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - } - - tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma); - if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { - binMax = newHistogram->bounds->n-1; - } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - } - if ((binMin < 0) || (binMax < 0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n"); - psFree(tmpMaskVec); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpStatsMinMax); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - - psS32 binNum = binMin; - psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum]; - for (psS32 i = binMin+1 ; i <= binMax ; i++) { - if (newHistogramSmoothed->data.F32[i] > binMaxNums) { - binNum = i; - binMaxNums = newHistogramSmoothed->data.F32[i]; - } - } - psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums); - - // - // Fit a Gaussian to the bins in the range 20 sigma of the robust - // histogram median. - // - tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma); - if (tmpScalar->data.F32 < tmpStatsMinMax->min) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - // XXX: check for errors here. - } - tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma); - if (tmpScalar->data.F32 > tmpStatsMinMax->max) { - binMax = newHistogramSmoothed->n - 1; - } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - // XXX: check for errors here. - } - psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psArray *x = psArrayAlloc((1 + (binMax - binMin))); - stats->robustNfit = 0; - psS32 j = 0; - - for (psS32 i = binMin ; i <= binMax ; i++) { - y->data.F32[j] = newHistogramSmoothed->data.F32[i]; - x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32); - ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i); - xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i); - - stats->robustNfit+= newHistogramSmoothed->data.F32[i]; - j++; - } - if (psTraceGetLevel(__func__) >= 8) { - // XXX: Print the x array somehow. - PS_VECTOR_PRINT_F32(y); - } - - // XXX: Use the min/max routines for this - psF32 minY = FLT_MAX; - psF32 maxY = -FLT_MAX; - for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { - if (y->data.F32[i] > maxY) { - maxY = y->data.F32[i]; - } - if (y->data.F32[i] < minY) { - minY = y->data.F32[i]; - } - } - // - // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) - // XXX: Use the normalize routines for this. - // - for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { - y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY); - } - - // - psMinimization *min = psMinimizationAlloc(100, 0.01); - psVector *params = psVectorAlloc(2, PS_TYPE_F32); - // Initial guess for the mean ([0]) and standard dev. - params->data.F32[0] = stats->robustMedian; - params->data.F32[1] = sigma; - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(params); - PS_VECTOR_PRINT_F32(y); - } - psFree(xTmp); - rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D); - - if (rcBool != true) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpScalar); + // ************************************************************************ + if ((stats->options & PS_STAT_FITTED_MEAN) || + (stats->options & PS_STAT_FITTED_STDEV)) { + // + // New algorithm for calculated mean and stdev. + // + + // XXX: Where is this documented? + psF32 dN = (psF32) (0.17 * N50); + if (dN < 1.0) { + dN = 1.0; + } else if (dN > 4.0) { + dN = 4.0; + } + psF32 newBinSize = sigma / dN; + + rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); + if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(tmpMaskVec); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + + numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize); + psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max); + psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize); + psTrace(__func__, 6, "The numBins is %d\n", numBins); + + psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); + newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(newHistogram->nums); + } + + // + // FITTED STATISTICS HERE + // + // Smooth the resulting histogram with a Gaussian with sigma_s = 1 + // bin. + // + psF32 sigma_s = newBinSize; + + psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(newHistogramSmoothed); + } + + // + // Find the bin with the peak value in the range 2 sigma of the + // robust histogram median. + // + + psS32 binMin = 0; + psS32 binMax = 0; + tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma); + if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + } + + tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma); + if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { + binMax = newHistogram->bounds->n-1; + } else { + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + } + if ((binMin < 0) || (binMax < 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n"); + psFree(tmpMaskVec); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpStatsMinMax); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + + psS32 binNum = binMin; + psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum]; + for (psS32 i = binMin+1 ; i <= binMax ; i++) { + if (newHistogramSmoothed->data.F32[i] > binMaxNums) { + binNum = i; + binMaxNums = newHistogramSmoothed->data.F32[i]; + } + } + psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums); + + // + // Fit a Gaussian to the bins in the range 20 sigma of the robust + // histogram median. + // + tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma); + if (tmpScalar->data.F32 < tmpStatsMinMax->min) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + // XXX: check for errors here. + } + tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma); + if (tmpScalar->data.F32 > tmpStatsMinMax->max) { + binMax = newHistogramSmoothed->n - 1; + } else { + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + // XXX: check for errors here. + } + psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psArray *x = psArrayAlloc((1 + (binMax - binMin))); + psS32 j = 0; + + for (psS32 i = binMin ; i <= binMax ; i++) { + y->data.F32[j] = newHistogramSmoothed->data.F32[i]; + x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32); + ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i); + xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i); + j++; + } + if (psTraceGetLevel(__func__) >= 8) { + // XXX: Print the x array somehow. + PS_VECTOR_PRINT_F32(y); + } + + // XXX: Use the min/max routines for this + psF32 minY = FLT_MAX; + psF32 maxY = -FLT_MAX; + for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + if (y->data.F32[i] > maxY) { + maxY = y->data.F32[i]; + } + if (y->data.F32[i] < minY) { + minY = y->data.F32[i]; + } + } + // + // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) + // XXX: Use the normalize routines for this. + // + for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY); + } + + // + psMinimization *min = psMinimizationAlloc(100, 0.01); + psVector *params = psVectorAlloc(2, PS_TYPE_F32); + // Initial guess for the mean ([0]) and standard dev. + params->data.F32[0] = stats->robustMedian; + params->data.F32[1] = sigma; + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(params); + PS_VECTOR_PRINT_F32(y); + } + psFree(xTmp); + rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D); + + if (rcBool != true) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(newHistogram); + psFree(x); + psFree(y); + psFree(min); + psFree(params); + psFree(tmpMaskVec); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(params); + } + + // + // The fitted mean mean_r is derived directly from the fitted + // Gaussian mean. + // + stats->fittedMean = params->data.F32[0]; + psTrace(__func__, 6, "The fitted mean is %f.\n", params->data.F32[0]); + + // + // The fitted standard deviation, SIGMA_r is determined by + // subtracting the smoothing scale in quadrature: + // SIGMA_r^2 = SIGMA^2 - sigma_s^2 + // + stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s)); + psTrace(__func__, 6, "The fitted stdev is %f.\n", stats->fittedStdev); + + psFree(newHistogramSmoothed); psFree(newHistogram); psFree(x); @@ -2161,36 +2204,9 @@ psFree(min); psFree(params); - psFree(tmpMaskVec); - psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); - return(1); - } - if (psTraceGetLevel(__func__) >= 8) { - PS_VECTOR_PRINT_F32(params); - } - - // - // The robust mean mean_r is derived directly from the fitted - // Gaussian mean. - // - stats->robustMean = params->data.F32[0]; - psTrace(__func__, 6, "The robust mean is %f.\n", params->data.F32[0]); - - // - // The robust standard deviation, SIGMA_r is determined by - // subtracting the smoothing scale in quadrature: - // SIGMA_r^2 = SIGMA^2 - sigma_s^2 - // - stats->robustStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s)); - psTrace(__func__, 6, "The robust stdev is %f.\n", stats->robustStdev); - - psFree(newHistogramSmoothed); + } + psFree(tmpStatsMinMax); psFree(cumulativeRobustHistogram); psFree(tmpScalar); - psFree(newHistogram); - psFree(x); - psFree(y); - psFree(min); - psFree(params); psFree(tmpMaskVec); psFree(robustHistogram); @@ -2236,16 +2252,15 @@ newStruct->sampleUQ = NAN; newStruct->sampleLQ = NAN; - newStruct->robustMean = NAN; newStruct->robustMedian = NAN; - newStruct->robustMode = NAN; newStruct->robustStdev = NAN; newStruct->robustUQ = NAN; newStruct->robustLQ = NAN; newStruct->robustN50 = -1; // XXX: This is never used - newStruct->robustNfit = -1; + newStruct->fittedMean = NAN; + newStruct->fittedStdev = NAN; + newStruct->fittedNfit = -1; newStruct->clippedMean = NAN; newStruct->clippedStdev = NAN; newStruct->clippedNvalues = -1; // XXX: This is never used - // XXX: Where do these values come from? newStruct->clipSigma = 3.0; newStruct->clipIter = 3; @@ -2253,7 +2268,8 @@ newStruct->max = NAN; newStruct->binsize = NAN; // XXX: This is never used + newStruct->nSubsample = 100000; newStruct->options = options; - psTrace(__func__, 3, "---- %s(psStats) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (newStruct); } @@ -2318,5 +2334,5 @@ newHist->uniform = true; - psTrace(__func__, 3, "---- %s(psHistogram) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (newHist); } @@ -2363,5 +2379,5 @@ newHist->uniform = false; - psTrace(__func__, 3, "---- %s(psHistogram) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (newHist); } @@ -2594,5 +2610,5 @@ } - psTrace(__func__, 3, "---- %s(psHistogram) end ----\n", __func__); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); return (out); } @@ -2674,5 +2690,5 @@ } - psTrace(__func__, 4,"---- %s(psVector) end ----\n", __func__); + psTrace(__func__, 4,"---- %s() end ----\n", __func__); return (tmp); } @@ -2708,4 +2724,5 @@ PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats); } + // XXX: Assert that "in" is F64, F32, U16, or S8 psVector* inF32 = NULL; @@ -2729,47 +2746,51 @@ } + if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) { + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, stats); + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEAN) { if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); - } - } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n"); + stats->sampleMean = NAN; + } + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEDIAN) { if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleMedian() returned an error.\n"); - } - } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian() returned an error.\n"); + stats->sampleMedian = NAN; + } + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_STDEV) { if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); - } - p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats); - } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n"); + stats->sampleMean = NAN; + } else { + p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats); + } + } + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_QUARTILE) { if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorStats(): p_psVectorSampleQuartiles() returned an error.\n"); - } - } - // Since the various robust stats quantities share much computation, they - // are grouped together in a single private function: - // p_psVectorRobustStats() - if ((stats->options & PS_STAT_ROBUST_MEAN) || - (stats->options & PS_STAT_ROBUST_MEDIAN) || - (stats->options & PS_STAT_ROBUST_MODE) || + psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleQuartiles() returned an error.\n"); + stats->sampleLQ = NAN; + stats->sampleUQ = NAN; + } + } + + // ************************************************************************ + if ((stats->options & PS_STAT_ROBUST_MEDIAN) || (stats->options & PS_STAT_ROBUST_STDEV) || - (stats->options & PS_STAT_ROBUST_QUARTILE)) { + (stats->options & PS_STAT_ROBUST_QUARTILE) || + (stats->options & PS_STAT_FITTED_MEAN) || + (stats->options & PS_STAT_FITTED_STDEV)) { if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, - PS_ERRORTEXT_psStats_STATS_FAILED); - // XXX: Set to NAN - // XXX: Is this the right thing to do? - // XXX: If so, do it for other funcs? + psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED); psFree(stats); psTrace(__func__, 3,"---- %s(NULL) end ----\n", __func__); @@ -2778,31 +2799,26 @@ } - // XXX: Different conditions for return -1 and -2? + // ************************************************************************ if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats); - if (-1 == rc) { - psError(PS_ERR_UNKNOWN, false, - "Failed to calculate clipped statistics for input psVector.\n"); + if (rc < 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n"); stats->clippedMean = NAN; stats->clippedStdev = NAN; - } else if (-2 == rc) { - psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector."); - stats->clippedMean = NAN; - stats->clippedStdev = NAN; - } - } + } + } + // ************************************************************************ if (stats->options & PS_STAT_MAX) { if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, - "Failed to calculate vector maximum"); + psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector maximum"); stats->max = NAN; } } + // ************************************************************************ if (stats->options & PS_STAT_MIN) { if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, - "Failed to calculate vector minimum"); + psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum"); stats->min = NAN; } @@ -2815,5 +2831,5 @@ psFree(errorsF32); } - psTrace(__func__, 3,"---- %s(psStats) end ----\n", __func__); + psTrace(__func__, 3,"---- %s() end ----\n", __func__); return (stats); } Index: /trunk/psLib/src/math/psStats.h =================================================================== --- /trunk/psLib/src/math/psStats.h (revision 6214) +++ /trunk/psLib/src/math/psStats.h (revision 6215) @@ -14,6 +14,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.48 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-08-30 01:14:13 $ + * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-01-26 23:49:11 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -37,19 +37,19 @@ // XXX: Is PS_STAT_ROBUST_FOR_SAMPLE obsolete? typedef enum { - PS_STAT_SAMPLE_MEAN = 0x000001, ///< Sample Mean - PS_STAT_SAMPLE_MEDIAN = 0x000002, ///< Sample Median - PS_STAT_SAMPLE_STDEV = 0x000004, ///< Sample Standard Deviation - PS_STAT_SAMPLE_QUARTILE = 0x000008, ///< Sample Quartile - PS_STAT_ROBUST_MEAN = 0x000010, ///< Robust Mean - PS_STAT_ROBUST_MEDIAN = 0x000020, ///< Robust Median - PS_STAT_ROBUST_MODE = 0x000040, ///< Robust Mode - PS_STAT_ROBUST_STDEV = 0x000080, ///< Robust Standarad Deviation - PS_STAT_ROBUST_QUARTILE = 0x000100, ///< Robust Quartile - PS_STAT_CLIPPED_MEAN = 0x000200, ///< Clipped Mean - PS_STAT_CLIPPED_STDEV = 0x000400, ///< Clipped Standard Deviation - PS_STAT_MAX = 0x000800, ///< Maximum - PS_STAT_MIN = 0x001000, ///< Minumum - PS_STAT_USE_RANGE = 0x002000, ///< Range - PS_STAT_USE_BINSIZE = 0x004000, ///< Binsize + PS_STAT_SAMPLE_MEAN = 0x000001, ///< Sample Mean + PS_STAT_SAMPLE_MEDIAN = 0x000002, ///< Sample Median + PS_STAT_SAMPLE_STDEV = 0x000004, ///< Sample Standard Deviation + PS_STAT_SAMPLE_QUARTILE = 0x000008, ///< Sample Quartile + PS_STAT_ROBUST_MEDIAN = 0x000010, ///< Robust Median + PS_STAT_ROBUST_STDEV = 0x000020, ///< Robust Standarad Deviation + PS_STAT_ROBUST_QUARTILE = 0x000040, ///< Robust Quartile + PS_STAT_FITTED_MEAN = 0x000080, ///< Fitted Mean + PS_STAT_FITTED_STDEV = 0x000100, ///< Fitted Standard Deviation + PS_STAT_CLIPPED_MEAN = 0x000200, ///< Clipped Mean + PS_STAT_CLIPPED_STDEV = 0x000400, ///< Clipped Standard Deviation + PS_STAT_MAX = 0x000800, ///< Maximum + PS_STAT_MIN = 0x001000, ///< Minumum + PS_STAT_USE_RANGE = 0x002000, ///< Range + PS_STAT_USE_BINSIZE = 0x004000, ///< Binsize } psStatsOptions; @@ -62,14 +62,14 @@ double sampleMedian; ///< formal median of sample double sampleStdev; ///< standard deviation of sample + double sampleLQ; ///< lower quartile of sample double sampleUQ; ///< upper quartile of sample - double sampleLQ; ///< lower quartile of sample - double robustMean; ///< robust mean of array double robustMedian; ///< robust median of array - double robustMode; ///< Robust mode of array double robustStdev; ///< robust standard deviation of array double robustUQ; ///< robust upper quartile double robustLQ; ///< robust lower quartile int robustN50; ///< Number of points in Gaussian fit. XXX: This is currently never set. - int robustNfit; ///< The number of points between the quartiles. + double fittedMean; ///< robust mean of data + double fittedStdev; ///< robust standard deviation of data + int fittedNfit; ///< Number of points in Gaussian fit double clippedMean; ///< Nsigma clipped mean double clippedStdev; ///< standard deviation after clipping @@ -80,4 +80,5 @@ double max; ///< maximum data value in array double binsize; ///< binsize for robust fit (input/ouput) + int nSubsample; ///< maxinum number of measurements (input) psStatsOptions options; ///< bitmask of calculated values } Index: /trunk/psLib/test/math/tst_psStats00.c =================================================================== --- /trunk/psLib/test/math/tst_psStats00.c (revision 6214) +++ /trunk/psLib/test/math/tst_psStats00.c (revision 6215) @@ -6,11 +6,11 @@ * do a much simpler test with data types PS_TYPE_U8, PS_TYPE_U16, PS_TYPE_F64. * -* If the psStats,c code every changes such that vectors of different type +* If the psStats,c code ever changes such that vectors of different type * are handled by different routines, then these tests must be extended. * * @author GLG, MHPCC * -* @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-24 01:10:00 $ +* @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-01-26 23:49:11 $ * * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii @@ -22,5 +22,5 @@ #define ERROR_TOL 0.0001 #define N 15 - +#define VERBOSE 1 static psS32 testStatsSampleMeanF32(void); static psS32 testStatsSampleMeanS8(void); @@ -71,5 +71,5 @@ psTraceSetLevel(".", 0); psTraceSetLevel("p_psGetStatValue", 0); - psTraceSetLevel("p_psVectorSampleMean", 0); + psTraceSetLevel("p_psVectorSampleMean", 10); psTraceSetLevel("p_psVectorMax", 0); psTraceSetLevel("p_psVectorMin", 0); @@ -119,4 +119,5 @@ } + printf("HMMM: myStats->options is (0x%x)\n", myStats->options); /*************************************************************************/ /* Call psVectorStats() with no vector mask. */ @@ -125,5 +126,8 @@ mean = myStats->sampleMean; // Verify return value is as expected - if ( fabs(mean - expectedMeanNoMaskF32) > ERROR_TOL ) { + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskF32); + } + if (isnan(myStats->sampleMean) || (fabs(mean - expectedMeanNoMaskF32) > ERROR_TOL)) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", mean, expectedMeanNoMaskF32); @@ -135,5 +139,8 @@ mean = myStats->sampleMean; // Verify return value is as expected - if ( fabs(mean - expectedWeightMeanNoMaskF32) > ERROR_TOL ) { + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanNoMaskF32); + } + if (isnan(myStats->sampleMean) || (fabs(mean - expectedWeightMeanNoMaskF32) > ERROR_TOL)) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", mean, expectedWeightMeanNoMaskF32); @@ -148,4 +155,7 @@ mean = myStats->sampleMean; // Verify return value is as expected + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanRangeNoMaskF32); + } if ( fabs(mean - expectedMeanRangeNoMaskF32) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f", @@ -158,4 +168,7 @@ mean = myStats->sampleMean; // Verify return value is as expected + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanNoMaskRangeF32); + } if ( fabs(mean - expectedWeightMeanNoMaskRangeF32) > ERROR_TOL) { psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f", @@ -170,4 +183,7 @@ myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); mean = myStats->sampleMean; + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanWithMaskF32); + } if ( fabs(mean - expectedMeanWithMaskF32) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", @@ -179,4 +195,7 @@ myStats = psVectorStats(myStats, myVector, myErrors, maskVector, 1); mean = myStats->sampleMean; + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanWithMaskF32); + } if ( fabs(mean - expectedWeightMeanWithMaskF32) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", @@ -190,4 +209,7 @@ mean = myStats->sampleMean; // Verify return value is as expected + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanRangeWithMaskF32); + } if ( fabs(mean - expectedMeanRangeWithMaskF32) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f", @@ -200,4 +222,7 @@ mean = myStats->sampleMean; // Verify return value is as expected + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanWithMaskRangeF32); + } if ( fabs(mean - expectedWeightMeanWithMaskRangeF32) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f", @@ -219,4 +244,7 @@ myStats = psVectorStats(myStats, myVector, NULL, maskVector, 2); mean = myStats->sampleMean; + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanWithMaskF32); + } if (fabs(mean - expectedMeanWithMaskF32) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", @@ -237,4 +265,7 @@ myStats = psVectorStats(myStats, myVector, NULL, maskVector, 4); mean = myStats->sampleMean; + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskF32); + } if (fabs(mean - expectedMeanNoMaskF32) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f", @@ -307,4 +338,7 @@ mean = myStats->sampleMean; // Verify return value is as expected + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskS8); + } if ( fabs(mean - expectedMeanNoMaskS8) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", @@ -347,4 +381,7 @@ mean = myStats->sampleMean; // Verify return value is as expected + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskU16); + } if ( fabs(mean - expectedMeanNoMaskU16) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", @@ -387,4 +424,7 @@ mean = myStats->sampleMean; // Verify return value is as expected + if (VERBOSE) { + printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskF64); + } if ( fabs(mean - expectedMeanNoMaskF64) > ERROR_TOL ) { psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f", Index: /trunk/psLib/test/math/tst_psStats03.c =================================================================== --- /trunk/psLib/test/math/tst_psStats03.c (revision 6214) +++ /trunk/psLib/test/math/tst_psStats03.c (revision 6215) @@ -18,13 +18,13 @@ // We list pertinent psStats.c functions here for debugging ease. // - psTraceSetLevel(".", 0); - psTraceSetLevel("p_psGetStatValue", 0); - psTraceSetLevel("p_psVectorSampleMean", 0); - psTraceSetLevel("p_psVectorCheckNonEmpty", 0); - psTraceSetLevel("p_psVectorNValues", 0); - psTraceSetLevel("p_psVectorSampleMedian", 0); - psTraceSetLevel("psStatsAlloc", 0); - psTraceSetLevel("p_psConvertToF32", 0); - psTraceSetLevel("psVectorStats", 0); + psTraceSetLevel(".", 10); + psTraceSetLevel("p_psGetStatValue", 10); + psTraceSetLevel("p_psVectorSampleMean", 10); + psTraceSetLevel("p_psVectorCheckNonEmpty", 10); + psTraceSetLevel("p_psVectorNValues", 10); + psTraceSetLevel("p_psVectorSampleMedian", 10); + psTraceSetLevel("psStatsAlloc", 10); + psTraceSetLevel("p_psConvertToF32", 10); + psTraceSetLevel("psVectorStats", 10); psStats *myStats = NULL; Index: /trunk/psLib/test/math/tst_psStats07.c =================================================================== --- /trunk/psLib/test/math/tst_psStats07.c (revision 6214) +++ /trunk/psLib/test/math/tst_psStats07.c (revision 6215) @@ -40,5 +40,4 @@ float realUQNoMask = MEAN + ( 0.6 * STDEV ); psS32 realN50NoMask = N / 4; - psS32 realNfitNoMask = N / 4; /*************************************************************************/ @@ -79,9 +78,9 @@ psFree(mySampleStatsWithMask); - myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN | - PS_STAT_ROBUST_MEDIAN | - PS_STAT_ROBUST_MODE | - PS_STAT_ROBUST_STDEV | - PS_STAT_ROBUST_QUARTILE ); + myStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | + PS_STAT_ROBUST_STDEV | + PS_STAT_ROBUST_QUARTILE | + PS_STAT_FITTED_MEAN | + PS_STAT_FITTED_STDEV); // Create a full outliers: myVector->data.F32[N/4] = -1000.0 * MEAN; @@ -101,7 +100,7 @@ printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n", - realMeanNoMask, myStats->robustMean ); - - if ( fabs( myStats->robustMean - realMeanNoMask ) < ( ERROR_TOLERANCE * realMeanNoMask ) ) { + realMeanNoMask, myStats->fittedMean ); + + if ( fabs( myStats->fittedMean - realMeanNoMask ) < ( ERROR_TOLERANCE * realMeanNoMask ) ) { testStatus = true; } else { @@ -157,6 +156,6 @@ printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n", - realStdevNoMask, myStats->robustStdev ); - if ( fabs( myStats->robustStdev - realStdevNoMask ) < ( ERROR_TOLERANCE * realStdevNoMask ) ) { + realStdevNoMask, myStats->fittedStdev ); + if ( fabs( myStats->fittedStdev - realStdevNoMask ) < ( ERROR_TOLERANCE * realStdevNoMask ) ) { testStatus = true; } else { @@ -208,5 +207,4 @@ "psStats functions", "PS_STAT_ROBUST_STATS: robust N50: no vector mask" ); - // XXX: realN50NoMask = myStats->robustN50; @@ -227,25 +225,4 @@ testStatus ); - printPositiveTestHeader( stdout, - "psStats functions", - "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask" ); - - // XXX: - realNfitNoMask = myStats->robustNfit; - - printf( "The expected Nfit was %d; the calculated Nfit was %d\n", - realNfitNoMask, myStats->robustNfit ); - /* XXX: fix - if ( fabs( myStats->robustNfit - realNfitNoMask ) < ( ERROR_TOLERANCE * realNfitNoMask ) ) { - testStatus = true; - } else { - testStatus = false; - globalTestStatus = false; - } - */ - printFooter( stdout, - "psVector functions", - "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask", - testStatus ); /*************************************************************************/ @@ -296,13 +273,12 @@ float realUQWithMask = MEAN; psS32 realN50WithMask = N / 4; - psS32 realNfitWithMask = N / 4; /*************************************************************************/ /* Allocate and initialize data structures */ /*************************************************************************/ - myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN | - PS_STAT_ROBUST_MEDIAN | - PS_STAT_ROBUST_MODE | - PS_STAT_ROBUST_STDEV | - PS_STAT_ROBUST_QUARTILE ); + myStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | + PS_STAT_ROBUST_STDEV | + PS_STAT_ROBUST_QUARTILE | + PS_STAT_FITTED_MEAN | + PS_STAT_FITTED_STDEV); maskVector = psVectorAlloc( N, PS_TYPE_U8 ); @@ -334,6 +310,6 @@ printf( "Called psVectorStats() on a vector with elements masked.\n" ); printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n", - realMeanWithMask, myStats->robustMean ); - if ( fabs( myStats->robustMean - realMeanWithMask ) < ( ERROR_TOLERANCE * realMeanWithMask ) ) { + realMeanWithMask, myStats->fittedMean ); + if ( fabs( myStats->fittedMean - realMeanWithMask ) < ( ERROR_TOLERANCE * realMeanWithMask ) ) { testStatus = true; } else { @@ -389,6 +365,6 @@ printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n", - realStdevWithMask, myStats->robustStdev ); - if ( fabs( myStats->robustStdev - realStdevWithMask ) < ( ERROR_TOLERANCE * realStdevWithMask ) ) { + realStdevWithMask, myStats->fittedStdev ); + if ( fabs( myStats->fittedStdev - realStdevWithMask ) < ( ERROR_TOLERANCE * realStdevWithMask ) ) { testStatus = true; } else { @@ -459,26 +435,4 @@ "PS_STAT_ROBUST_STATS: robust N50: with vector mask", testStatus ); - - - - printPositiveTestHeader( stdout, - "psStats functions", - "PS_STAT_ROBUST_STATS: robust Nfit: with vector mask" ); - - printf( "The expected Nfit was %d; the calculated Nfit was %d\n", - realNfitWithMask, myStats->robustNfit ); - /* XXX: fix - if ( fabs( myStats->robustNfit - realNfitWithMask ) < ( ERROR_TOLERANCE * realNfitWithMask ) ) { - testStatus = true; - } else { - testStatus = false; - globalTestStatus = false; - } - */ - printFooter( stdout, - "psVector functions", - "PS_STAT_ROBUST_STATS: robust Nfit: with vector mask", - testStatus ); -