Index: trunk/psLib/src/fits/psFitsTable.c =================================================================== --- trunk/psLib/src/fits/psFitsTable.c (revision 30995) +++ trunk/psLib/src/fits/psFitsTable.c (revision 31152) @@ -645,5 +645,12 @@ psMetadata *row = table->data[i]; // The row of interest psMetadataItem *dataItem = psMetadataLookup(row, colSpecItem->name); // Value of interest - memcpy(&columnData->data.U8[i * dataSize], &dataItem->data, dataSize); + if (dataItem) { + memcpy(&columnData->data.U8[i * dataSize], &dataItem->data, dataSize); + } else { + // this element is missing from this row; insert an appropriate-sized place holder + // XXX this should insert a NAN for float / double and an appropriate blank for int types + // XXX for the moment I am putting in 0.0 + memset(&columnData->data.U8[i * dataSize], 0, dataSize); + } } Index: trunk/psLib/src/imageops/psImageBackground.c =================================================================== --- trunk/psLib/src/imageops/psImageBackground.c (revision 30995) +++ trunk/psLib/src/imageops/psImageBackground.c (revision 31152) @@ -13,4 +13,5 @@ #include "psType.h" #include "psAssert.h" +#include "psAbort.h" #include "psRandom.h" #include "psError.h" @@ -23,6 +24,4 @@ } -// XXX allow the user to choose the stats method? -// (SAMPLE_MEAN, CLIPPED_MEAN, ROBUST_MEDIAN, FITTED_MEAN) bool psImageBackground(psStats *stats, psVector **sample, const psImage *image, const psImage *mask, psImageMaskType maskValue, psRandom *rng) { @@ -45,28 +44,59 @@ long nx = image->numCols; long ny = image->numRows; - - psImage *bad = psImageAlloc(nx, ny, PS_TYPE_U8); // Image with bad pixels - psImageInit(bad, 0); - - int Npixels = 0; // Total number of pixels - for (int y = 0; y < ny; y++) { - for (int x = 0; x < nx; x++) { - if (!isfinite(image->data.F32[y][x]) || - (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue)) { - bad->data.U8[y][x] = 0xFF; - } - Npixels++; - } - } - - const int Nsubset = (stats->nSubsample == 0) ? Npixels : PS_MIN(stats->nSubsample, Npixels); // Number of pixels in subset - - psVector *values; // Vector containing subsample + long nPixels = nx * ny; + + long nSubset = (stats->nSubsample == 0) ? (nx * ny) : PS_MIN(stats->nSubsample, (nx * ny)); + + /*** discussion of the sampling and algorithm choices: + + We want stats->nSubsample pixels. How many good pixels do we actually have? + + There are a few domains of interest: + + If nGoodPixels < f1 * nSubset, we should fail [1] + If nSubset >= nGoodPixels, we should just loop over all pixels [2] + If (nSubset < f2 * nGoodPixels) and (nGoodPixels >= f3 * nPixels), we should just use the random sample technique [3] + If (nSubset >= f2 * nGoodPixels) or (nGoodPixel < f3 * nPixels), we should use the Fisher-Yates technique. + + to use Fisher-Yates, we need to have a copy of the pixels so we can re-shuffle. We + should not do that with the whole list unless we need to. + + [1] we just have too few samples to be useful; f1 ~ 0.01? + + [2] we need to select from all pixels to reach the desired sample size + + [3] this avoids a full-image-sized alloc, but only makes sense if nGoodPixels is actually + large. f2 ~ 0.01, f3 ~ 0.1 (4 tries per success on average) + + ***/ + + // count the number of good pixels + long nGoodPixels = 0; + for (long iy = 0; iy < ny; iy++) { + for (long ix = 0; ix < nx; ix++) { + if (!isfinite(image->data.F32[iy][ix])) continue; + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue; + nGoodPixels ++; + } + } + + // XXX This value of 1% is somewhat arbitrary + if (nGoodPixels < 0.01*nSubset) { + if ((nFailures < 3) || (nFailures % 100 == 0)) { + psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", nGoodPixels, nFailures); + } + nFailures ++; + psTrace ("psLib.imageops", 4, "case 1: nGoodPixels < 0.01*nSubset (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels); + return false; + } + + // alloc or recycle the vector containing subsample pixels + psVector *values; if (sample) { - *sample = psVectorRecycle(*sample, Nsubset, PS_TYPE_F32); + *sample = psVectorRecycle(*sample, nSubset, PS_TYPE_F32); values = psMemIncrRefCounter(*sample); values->n = 0; } else { - values = psVectorAllocEmpty(Nsubset, PS_TYPE_F32); + values = psVectorAllocEmpty(nSubset, PS_TYPE_F32); } @@ -75,56 +105,64 @@ float max = -PS_MAX_F32; - // select a subset of the image pixels to measure the stats - long n = 0; // Number of actual pixels in subset - if (Nsubset >= Npixels) { - // if we have an image smaller than Nsubset, just loop over the (good) image pixels - for (int iy = 0; iy < ny; iy++) { - for (int ix = 0; ix < nx; ix++) { - if (bad->data.U8[iy][ix]) { - continue; - } - - float value = image->data.F32[iy][ix]; - min = PS_MIN(value, min); - max = PS_MAX(value, max); - values->data.F32[n] = value; - n++; - } - } - } else { - // Subsample all pixels - // This is not optimal, since there may be a large masked fraction that leaves us with few good pixels. - // In that case, you should have set Nsubset....... - Npixels = nx * ny; - for (long i = 0; i < Nsubset; i++) { + if ((nSubset < 0.01*nGoodPixels) && (nGoodPixels >= 0.1*nPixels)) { + psTrace ("psLib.imageops", 4, "case 3: nSubset < 0.01*nGoodPixels && nGoodPixels > 0.1*nPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels); + // for small subsets, just use simple random sampling (saves a full-image-sized alloc), + // unless we have lots of masked pixels + for (long i = 0; i < nSubset; i++) { double frnd = psRandomUniform(rng); - int pixel = Npixels * frnd; - int ix = pixel % nx; - int iy = pixel / nx; - - if (bad->data.U8[iy][ix]) { - continue; - } + long pixel = nPixels * frnd; + long ix = pixel % nx; + long iy = pixel / nx; + + if (!isfinite(image->data.F32[iy][ix])) continue; + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue; float value = image->data.F32[iy][ix]; min = PS_MIN(value, min); max = PS_MAX(value, max); - values->data.F32[n] = value; - n++; - } - } - - psFree(bad); - - if (n < 0.01*Nsubset) { - if ((nFailures < 3) || (nFailures % 100 == 0)) { - psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", n, nFailures); - } - nFailures ++; - psFree(values); - return false; - } - - values->n = n; + psVectorAppend (values, value); + } + } else if (nSubset >= nGoodPixels) { + psTrace ("psLib.imageops", 4, "case 2: nSubset >= nGoodPixels (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels); + // in this case, we have to select from all masked pixels just to get the desired + // sample size + for (long iy = 0; iy < ny; iy++) { + for (long ix = 0; ix < nx; ix++) { + if (!isfinite(image->data.F32[iy][ix])) continue; + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue; + float value = image->data.F32[iy][ix]; + min = PS_MIN(value, min); + max = PS_MAX(value, max); + psVectorAppend(values, value); + } + } + } else { + psTrace ("psLib.imageops", 4, "case 4: Fisher-Yates (%d x %d : %d : %d : %d)\n", (int) nx, (int) ny, (int) nSubset, (int) nGoodPixels, (int) nPixels); + // use Knuth's version of Fisher-Yates to select a random subset of the pixels, but + // only drawing each value once + + // generate a vector of all pixels which may in theory be selected + psVector *pixelVector = psVectorAllocEmpty(nGoodPixels, PS_TYPE_F32); + for (long iy = 0; iy < ny; iy++) { + for (long ix = 0; ix < nx; ix++) { + if (!isfinite(image->data.F32[iy][ix])) continue; + if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue; + psVectorAppend(pixelVector, image->data.F32[iy][ix]); + } + } + psAssert (nGoodPixels == pixelVector->n, "we must have mis-counted above"); + + // generate the unique random subset + for (int i = 0; i < nSubset; i++) { + double frnd = psRandomUniform(rng); + int pixel = pixelVector->n * frnd; + + psVectorAppend(values, pixelVector->data.F32[pixel]); + pixelVector->data.F32[pixel] = pixelVector->data.F32[pixelVector->n - 1]; + pixelVector->n --; + } + psFree (pixelVector); + } + psAssert (values->n >= 0.01*nSubset, "didn't we test this above?"); if (stats->options & PS_STAT_ROBUST_QUARTILE) { @@ -132,6 +170,6 @@ // XXX this hack is just for testing, drop when I am happy with the psStats version of the values - int imin = stats->min * n; - int imax = stats->max * n; + int imin = stats->min * values->n; + int imax = stats->max * values->n; int npts = imax - imin + 1; @@ -147,7 +185,7 @@ // Subtract the median when we add the numbers, so we don't get numerical problems float median = npts % 2 ? 0.5 * (values->data.F32[npts/2 - 1] + values->data.F32[npts/2]) : - values->data.F32[npts/2]; + values->data.F32[npts/2]; double value = 0; - for (long i = imin; (i <= imax) && (i < n); i++) { + for (long i = imin; (i <= imax) && (i < values->n); i++) { value += values->data.F32[i] - median; } @@ -185,2 +223,32 @@ return true; } + + +/*** old comments + + // Subsample all pixels + // This is not optimal, since there may be a large masked fraction that leaves us with few good pixels. + // In that case, you should have set Nsubset....... + + // 2011/03/22 - MWV: On prompting from Gene, changed the sampling so that it doesn't build and then sort + // a huge array when Nsubset << Npixel. + // O(Npixel) log(Npixel) is unnecessarily painful when Npixel=38 million (one skycell) + // but we only needed Nsubset pixels and Nsubset << Npixel. + + // It also perhaps partially defeats any gain of taking only subsamples if we're sorting the array. + // I should do some performance tests on this anyway. + + // If our chance of collision is low and the gain from not sorting the array is likely to be high, + // then just pick randomnly + // Note that we're guaranteeing Nsubset pixels returned here (up to the limit of Npixels) + + // 2011/03/16 - MWV: Fixed double-sampling of sky pixels. + // The pick-a-random-pixel-at-a-time approach doesn't work well when Npixels is close to Nsubset + // If Nsubset is 0.8*Npixels, then we will get lots of pixels double-counted + // This is correct on average, but isn't the optimal way to estimate the sky level + // Instead we take the following approach: + // 1) Calculate a random ordering of the pixels + // 2) Go through this ordering up to Nsubset to select pixels + // This is O(n log(n)) + + ***/ Index: trunk/psLib/src/math/psMinimizePolyFit.c =================================================================== --- trunk/psLib/src/math/psMinimizePolyFit.c (revision 30995) +++ trunk/psLib/src/math/psMinimizePolyFit.c (revision 31152) @@ -767,6 +767,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); @@ -1211,6 +1211,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); @@ -1621,6 +1621,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); @@ -2055,6 +2055,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 30995) +++ trunk/psLib/src/math/psStats.c (revision 31152) @@ -172,8 +172,4 @@ *****************************************************************************/ -// static prototypes: -static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords); -// static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); -// static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); @@ -1087,13 +1083,16 @@ } -/* - * vectorFittedStats requires guess for fittedMean and fittedStdev - * robustN50 should also be set - */ +/******************** + * perform an asymmetric fit to the population. In development, this was called + * "vectorFittedStats_v4" all versions of fitted stats now resolve to this function (only v4 + * has really been used) vectorFittedStats requires guess for fittedMean and fittedStdev + * robustN50 should also be set gaussian fit is performed using 2D polynomial to ln(y) this + * version follows the upper portion of the distribution until it passes 0.5*peak + ********************/ static bool vectorFittedStats (const psVector* myVector, - const psVector* errors, - psVector* mask, - psVectorMaskType maskVal, - psStats* stats) + const psVector* errors, + psVector* mask, + psVectorMaskType maskVal, + psStats* stats) { @@ -1121,401 +1120,4 @@ stats->results |= PS_STAT_FITTED_MEAN; stats->results |= PS_STAT_FITTED_STDEV; - return true; - } - - float guessStdev = stats->robustStdev; // pass the guess sigma - float guessMean = stats->robustMedian; // pass the guess mean - - psTrace(TRACE, 6, "The guess mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The guess stdev is %f.\n", guessStdev); - - bool done = false; - for (int iteration = 0; !done && (iteration < 2); iteration ++) { - psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - - psF32 binSize = 1; - if (stats->options & PS_STAT_USE_BINSIZE) { - // Set initial bin size to the specified value. - binSize = stats->binsize; - psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { - // construct a histogram with (sigma/10 < binsize < sigma) - // set roughly so that the lowest bins have about 2 cnts - // Nsmallest ~ N50 / (4*dN)) - psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8)); - binSize = guessStdev / dN; - } - - // Determine the min/max of the vector (which prior outliers masked out) - int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values - float min = statsMinMax->min; - float max = statsMinMax->max; - if (numValid == 0 || isnan(min) || isnan(max)) { - psTrace(TRACE, 5, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); - return true; - } - - // Calculate the number of bins. - // XXX can we calculate the binMin, binMax **before** building this histogram? - long numBins = (max - min) / binSize; - psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max); - psTrace(TRACE, 6, "The new bin size is %f.\n", binSize); - psTrace(TRACE, 6, "The numBins is %ld\n", numBins); - - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - // if psVectorHistogram returns false, we have a programming error - psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n"); - psFree(histogram); - psFree(statsMinMax); - return false; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - // Fit a Gaussian to the bins in the requested range about the guess mean - // min,max overrides clipSigma - psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - maxFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->max)) { - maxFitSigma = fabs(stats->max); - } - - psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - minFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minFitSigma = fabs(stats->min); - } - - // select the min and max bins, saturating on the lower and upper end-points - long binMin, binMax; - PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0); - PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0); - - // Generate the variables that will be used in the Gaussian fitting - // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins - psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates - psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates - for (long i = binMin, j = 0; i <= binMax ; i++, j++) { - y->data.F32[j] = histogram->nums->data.F32[i]; - psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value - ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i); - x->data[j] = ordinate; - } - if (psTraceGetLevel("psLib.math") >= 8) { - // XXX: Print the x array somehow. - PS_VECTOR_PRINT_F32(y); - } - psTrace(TRACE, 6, "The clipped numBins is %ld\n", y->n); - psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax); - - // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) - // XXX EAM : why not just fit the normalization as well?? - p_psNormalizeVectorRange(y, 0.0, 1.0); - - // Fit a Gaussian to the data. - psMinimization *minimizer = psMinimizationAlloc(100, 0.01, 1.0); // The minimizer information - psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian - // Initial guess for the mean (index 0) and var (index 1). - params->data.F32[0] = guessMean; - params->data.F32[1] = PS_SQR(guessStdev); - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(params); - PS_VECTOR_PRINT_F32(y); - } - - // psMinimizeLMChi2 can return false for bad data as well as for serious failures - if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(params); - psFree(minimizer); - psFree(x); - psFree(y); - psFree(histogram); - psFree(statsMinMax); - return true; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(params); - } - - guessMean = params->data.F32[0]; - guessStdev = sqrt(params->data.F32[1]); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } - - // Clean up after fitting - psFree (params); - psFree (minimizer); - psFree (x); - psFree (y); - psFree (histogram); - psFree (statsMinMax); - } - - // The fitted mean is the Gaussian mean. - stats->fittedMean = guessMean; - psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean); - - // The fitted standard deviation - stats->fittedStdev = guessStdev; - psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - - stats->results |= PS_STAT_FITTED_MEAN; - stats->results |= PS_STAT_FITTED_STDEV; - - return true; -} - - -/******************** - * vectorFittedStats_v2 requires guess for fittedMean and fittedStdev - * robustN50 should also be set - * gaussian fit is performed using 2D polynomial to ln(y) - ********************/ -static bool vectorFittedStats_v2 (const psVector* myVector, - const psVector* errors, - psVector* mask, - psVectorMaskType maskVal, - psStats* stats) -{ - - // This procedure requires the mean. If it has not been already - // calculated, then call vectorSampleMean() - if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) { - if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } - } - - // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (isnan(stats->robustMedian)) { - stats->fittedMean = NAN; - stats->fittedStdev = NAN; - stats->results |= PS_STAT_FITTED_MEAN_V2; - stats->results |= PS_STAT_FITTED_STDEV_V2; - return true; - } - - if (stats->robustStdev <= FLT_EPSILON) { - stats->fittedMean = stats->robustMedian; - stats->fittedStdev = stats->robustStdev; - stats->results |= PS_STAT_FITTED_MEAN_V2; - stats->results |= PS_STAT_FITTED_STDEV_V2; - return true; - } - - float guessStdev = stats->robustStdev; // pass the guess sigma - float guessMean = stats->robustMedian; // pass the guess mean - - psTrace(TRACE, 6, "The ** starting ** guess mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The ** starting ** guess stdev is %f.\n", guessStdev); - - bool done = false; - for (int iteration = 0; !done && (iteration < 2); iteration ++) { - psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - - psF32 binSize = 1; - if (stats->options & PS_STAT_USE_BINSIZE) { - // Set initial bin size to the specified value. - binSize = stats->binsize; - psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { - // construct a histogram with (sigma/10 < binsize < sigma) - // set roughly so that the lowest bins have about 2 cnts - // Nsmallest ~ N50 / (4*dN)) - psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8)); - binSize = guessStdev / dN; - } - - // Determine the min/max of the vector (which prior outliers masked out) - int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values - float min = statsMinMax->min; - float max = statsMinMax->max; - if (numValid == 0 || isnan(min) || isnan(max)) { - COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return true; - } - - // Calculate the number of bins. - // XXX can we calculate the binMin, binMax **before** building this histogram? - long numBins = (max - min) / binSize; - psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max); - psTrace(TRACE, 6, "The new bin size is %f.\n", binSize); - psTrace(TRACE, 6, "The numBins is %ld\n", numBins); - - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistcs v2.\n"); - psFree(histogram); - psFree(statsMinMax); - return false; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - // Fit a Gaussian to the bins in the requested range about the guess mean - // min,max overrides clipSigma - psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - maxFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->max)) { - maxFitSigma = fabs(stats->max); - } - - psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - minFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minFitSigma = fabs(stats->min); - } - - // select the min and max bins, saturating on the lower and upper end-points - long binMin, binMax; - PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0); - PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0); - - // Generate the variables that will be used in the Gaussian fitting - // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins - psVector *y = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates - psVector *x = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of ordinates - long j = 0; - for (long i = binMin; i <= binMax ; i++) { - if (histogram->nums->data.F32[i] <= 0.0) - continue; - x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); - y->data.F32[j] = log(histogram->nums->data.F32[i]); - // XXX note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) - j++; - } - y->n = x->n = j; - if (psTraceGetLevel("psLib.math") >= 8) { - // XXX: Print the x array somehow. - PS_VECTOR_PRINT_F32(y); - } - psTrace(TRACE, 6, "The clipped numBins is %ld\n", y->n); - psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax); - - // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 - psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - - // XXX how can we protect against some extreme outliers? the robust histogram - // should have already dealt with those, but we are again sensitive to them... - // psStats *fitStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - // fitStats->clipIter = 3.0; - // fitStats->clipSigma = 3.0; - // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_VECTOR_MASK); - // psVectorInit (fitMask, 0); - - // XXX not sure if these should result in errors or not... - if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(x); - psFree(y); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - if (poly->coeff[2] >= 0.0) { - psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psError(PS_ERR_UNKNOWN, false, "fit did not converge\n"); - psFree(x); - psFree(y); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } else { - psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The new guess mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The new guess stdev is %f.\n", guessStdev); - } - - // Clean up after fitting - psFree (x); - psFree (y); - psFree (poly); - // psFree (fitStats); - // psFree (fitMask); - psFree (histogram); - psFree (statsMinMax); - } - - // The fitted mean is the Gaussian mean. - stats->fittedMean = guessMean; - psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean); - - // The fitted standard deviation - stats->fittedStdev = guessStdev; - psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - - stats->results |= PS_STAT_FITTED_MEAN_V2; - stats->results |= PS_STAT_FITTED_STDEV_V2; - - return true; -} - -/******************** - * perform an asymmetric fit to the population - * vectorFittedStats_v3 requires guess for fittedMean and fittedStdev - * robustN50 should also be set - * gaussian fit is performed using 2D polynomial to ln(y) - ********************/ -static bool vectorFittedStats_v3 (const psVector* myVector, - const psVector* errors, - psVector* mask, - psVectorMaskType maskVal, - psStats* stats) -{ - - // This procedure requires the mean. If it has not been already - // calculated, then call vectorSampleMean() - if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) { - if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } - } - - // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (isnan(stats->robustMedian)) { - stats->fittedMean = NAN; - stats->fittedStdev = NAN; - stats->results |= PS_STAT_FITTED_MEAN_V3; - stats->results |= PS_STAT_FITTED_STDEV_V3; - return true; - } - - if (stats->robustStdev <= FLT_EPSILON) { - stats->fittedMean = stats->robustMedian; - stats->fittedStdev = stats->robustStdev; - stats->results |= PS_STAT_FITTED_MEAN_V3; - stats->results |= PS_STAT_FITTED_STDEV_V3; return true; } @@ -1551,310 +1153,4 @@ COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return true; - } - - // Calculate the number of bins. - // XXX can we calculate the binMin, binMax **before** building this histogram? - long numBins = (max - min) / binSize; - psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max); - psTrace(TRACE, 6, "The new bin size is %f.\n", binSize); - psTrace(TRACE, 6, "The numBins is %ld\n", numBins); - - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics v3.\n"); - psFree(histogram); - psFree(statsMinMax); - return false; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - // now fit a Gaussian to the upper and lower halves about the peak independently - - // set the full-range upper and lower limits - psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - maxFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->max)) { - maxFitSigma = fabs(stats->max); - } - - psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - minFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minFitSigma = fabs(stats->min); - } - - // select the min and max bins, saturating on the lower and upper end-points - long binMin, binMax; - PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0); - PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0); - if (binMin == binMax) { - COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); - return true; - } - - // search for mode (peak of histogram within range mean-2sigma - mean+2sigma - long binPeak = binMin; - float valPeak = histogram->nums->data.F32[binPeak]; - for (int i = binMin; i < binMax; i++) { - if (histogram->nums->data.F32[i] > valPeak) { - binPeak = i; - valPeak = histogram->nums->data.F32[binPeak]; - } - psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); - } - psTrace (TRACE, 6, "\n"); - - // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak - - psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin); - psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); - psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binPeak), binPeak); - psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); - - { - // fit the lower half of the distribution - // run down until we drop below 0.25*valPeak - long binS = binMin; - long binE = PS_MIN (binPeak + 3, binMax); - for (int i = binPeak-3; i >= binMin; i--) { - if (histogram->nums->data.F32[i] < 0.25*valPeak) { - binS = i; - break; - } - } - psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binS), binS); - psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binE), binE); - - psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates - psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates - long j = 0; - for (long i = binS; i < binE; i++) { - if (histogram->nums->data.F32[i] <= 0.0) - continue; - x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); - // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) - y->data.F32[j] = log(histogram->nums->data.F32[i]); - j++; - } - y->n = x->n = j; - - // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 - psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); - psFree(x); - psFree(y); - - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - if (poly->coeff[2] >= 0.0) { - psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - - // sometimes, the guessStdev is much too large. in this case, the entire real population - // tends to be found in a single bin. make one attempt to recover by dropping the guessStdev - // down by a jump and trying again - if (iteration == 0) { - guessStdev = 0.25*guessStdev; - psTrace(TRACE, 6, "*** retry, new stdev is %f.\n", guessStdev); - continue; - } - - psError(PS_ERR_UNKNOWN, false, "fit did not converge\n"); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - // calculate lower mean & stdev from parabolic fit -- use this as the result - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } - psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The lower mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev); - - psFree(poly); - } - - // for test, measure the same result for the upper section - { - // fit the upper half of the distribution - // run up until we drop below 0.25*valPeak - long binS = PS_MAX (binPeak - 3, 0); - long binE = binMax; - for (int i = binPeak+3; i < binMax; i++) { - if (histogram->nums->data.F32[i] < 0.25*valPeak) { - binE = i; - break; - } - } - psTrace(TRACE, 6, "Lower bound for upper half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binS), binS); - psTrace(TRACE, 6, "Upper bound for upper half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binE), binE); - - psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates - psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates - long j = 0; - for (long i = binS; i < binE; i++) { - if (histogram->nums->data.F32[i] <= 0.0) - continue; - x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); - // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) - y->data.F32[j] = log(histogram->nums->data.F32[i]); - j++; - } - y->n = x->n = j; - - // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 - psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); - psFree(x); - psFree(y); - - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - // calculate upper mean & stdev from parabolic fit -- ignore this value - float upperStdev = sqrt(-0.5/poly->coeff[2]); - float upperMean = poly->coeff[1]*PS_SQR(upperStdev); -#ifndef PS_NO_TRACE - psTrace(TRACE, 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The upper mean is %f.\n", upperMean); - psTrace(TRACE, 6, "The upper stdev is %f.\n", upperStdev); -#endif - - // if the resulting value is outside of the range binMin - binMax, use the upper value - if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) { - guessMean = upperMean; - guessStdev = upperStdev; - } - - psFree (poly); - } - - // Clean up after fitting - psFree (histogram); - psFree (statsMinMax); - } - - // The fitted mean is the Gaussian mean. - stats->fittedMean = guessMean; - psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean); - - // The fitted standard deviation - stats->fittedStdev = guessStdev; - psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - - stats->results |= PS_STAT_FITTED_MEAN_V3; - stats->results |= PS_STAT_FITTED_STDEV_V3; - - return true; -} - -/******************** - * perform an asymmetric fit to the population - * vectorFittedStats_v4 requires guess for fittedMean and fittedStdev - * robustN50 should also be set - * gaussian fit is performed using 2D polynomial to ln(y) - * this version follows the upper portion of the distribution until it passes 0.5*peak - ********************/ -static bool vectorFittedStats_v4 (const psVector* myVector, - const psVector* errors, - psVector* mask, - psVectorMaskType maskVal, - psStats* stats) -{ - - // This procedure requires the mean. If it has not been already - // calculated, then call vectorSampleMean() - if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) { - if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } - } - - // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (isnan(stats->robustMedian)) { - stats->fittedMean = NAN; - stats->fittedStdev = NAN; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; - return true; - } - - if (stats->robustStdev <= FLT_EPSILON) { - stats->fittedMean = stats->robustMedian; - stats->fittedStdev = stats->robustStdev; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; - return true; - } - - float guessStdev = stats->robustStdev; // pass the guess sigma - float guessMean = stats->robustMedian; // pass the guess mean - - psTrace(TRACE, 6, "The ** starting ** guess mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The ** starting ** guess stdev is %f.\n", guessStdev); - - bool done = false; - for (int iteration = 0; !done && (iteration < 2); iteration ++) { - psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - - psF32 binSize = 1; - if (stats->options & PS_STAT_USE_BINSIZE) { - // Set initial bin size to the specified value. - binSize = stats->binsize; - psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { - // construct a histogram with (sigma/2 < binsize < sigma) - // set roughly so that the lowest bins have about 2 cnts - // Nsmallest ~ N50 / (4*dN)) - psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); - binSize = guessStdev / dN; - } - - // Determine the min/max of the vector (which prior outliers masked out) - int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values - float min = statsMinMax->min; - float max = statsMinMax->max; - if (numValid == 0 || isnan(min) || isnan(max)) { - COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); goto escape; } @@ -1865,6 +1161,6 @@ stats->fittedMean = min; stats->fittedStdev = 0.0; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; } @@ -1933,12 +1229,12 @@ // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak - psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin); psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); + psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak); psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); + float lowfitMean = NAN; + float lowfitStdev = NAN; { // fit the lower half of the distribution @@ -2014,24 +1310,21 @@ // calculate lower mean & stdev from parabolic fit -- use this as the result - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } - psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The lower mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev); + lowfitStdev = sqrt(-0.5/poly->coeff[2]); + lowfitMean = poly->coeff[1]*PS_SQR(lowfitStdev); + + psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + psTrace(TRACE, 6, "The lower mean is %f.\n", lowfitMean); + psTrace(TRACE, 6, "The lower stdev is %f.\n", lowfitStdev); psFree(poly); } - // if we converge on a solution outside the range binMin - binMax, use a more conservative range - float minValue = PS_BIN_MIDPOINT(histogram, binMin); - float maxValue = PS_BIN_MIDPOINT(histogram, binMax - 1); - - if (done && ((guessMean < minValue) || (guessMean > maxValue))) { - psTrace(TRACE, 6, "Inconsistent result, re-trying the fit\n"); - + float fullfitMean = NAN; + float fullfitStdev = NAN; + float minValueSym = NAN; + float maxValueSym = NAN; + + // try the full fit as well: + { // fit a symmetric distribution // run up until we drop below 0.15*valPeak @@ -2085,26 +1378,25 @@ // calculate upper mean & stdev from parabolic fit -- ignore this value - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); + fullfitStdev = sqrt(-0.5/poly->coeff[2]); + fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev); #ifndef PS_NO_TRACE - psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The symmetric mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The symmetric stdev is %f.\n", guessStdev); + psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + psTrace(TRACE, 6, "The symmetric mean is %f.\n", fullfitMean); + psTrace(TRACE, 6, "The symmetric stdev is %f.\n", fullfitStdev); #endif // if we converge on a solution outside the range binMin - binMax, use a more conservative range - float minValueSym = PS_BIN_MIDPOINT(histogram, binS); - float maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1); + minValueSym = PS_BIN_MIDPOINT(histogram, binS); + maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1); // saturate on min or max value - if (guessMean < minValueSym) { - guessMean = minValueSym; + if (fullfitMean < minValueSym) { + fullfitMean = minValueSym; psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); } // saturate on min or max value - if (guessMean > maxValueSym) { - guessMean = maxValueSym; + if (fullfitMean > maxValueSym) { + fullfitMean = maxValueSym; psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); } @@ -2112,4 +1404,24 @@ psFree (poly); } + + // we now have the fullfit and the lowfit mean and stdev values + // accept the fullfit unless minValueSym < lowfitMean < fullfitMean + + if (isfinite(lowfitMean) && isfinite(lowfitStdev) && (lowfitMean < fullfitMean) && (lowfitMean > minValueSym)) { + guessMean = lowfitMean; + guessStdev = lowfitStdev; + } else { + guessMean = fullfitMean; + guessStdev = fullfitStdev; + } + + if (!isfinite(guessMean) || !isfinite(guessStdev)) { + guessMean = stats->robustMedian; + guessStdev = stats->robustStdev; + } + + if (guessStdev > 0.75*stats->robustStdev) { + done = true; + } // Clean up after fitting @@ -2126,6 +1438,6 @@ psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; @@ -2134,6 +1446,6 @@ stats->fittedMean = NAN; stats->fittedStdev = NAN; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; @@ -2442,37 +1754,4 @@ // ************************************************************************ - if (stats->options & (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2)) { - if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V2"); - } - if (!vectorFittedStats_v2(inF32, errorsF32, maskVector, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); - status &= false; - } - } - - // ************************************************************************ - if (stats->options & (PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_STDEV_V3)) { - if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V3"); - } - if (!vectorFittedStats_v3(inF32, errorsF32, maskVector, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); - status &= false; - } - } - - // ************************************************************************ - if (stats->options & (PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4)) { - if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V4"); - } - if (!vectorFittedStats_v4(inF32, errorsF32, maskVector, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); - status &= false; - } - } - - // ************************************************************************ if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) { @@ -2513,16 +1792,16 @@ READ_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); READ_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE); - READ_STAT("FITTED", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); - READ_STAT("FITTED_V2", PS_STAT_FITTED_MEAN_V2); - READ_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN_V2); - READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2); - READ_STAT("FITTED_V3", PS_STAT_FITTED_MEAN_V3); - READ_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN_V3); - READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3); - READ_STAT("FITTED_V4", PS_STAT_FITTED_MEAN_V4); - READ_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN_V4); - READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4); + READ_STAT("FITTED", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); + READ_STAT("FITTED_V2", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV); + READ_STAT("FITTED_V3", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV); + READ_STAT("FITTED_V4", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV); READ_STAT("CLIPPED", PS_STAT_CLIPPED_MEAN); READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); @@ -2554,10 +1833,4 @@ WRITE_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); WRITE_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); - WRITE_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN_V2); - WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2); - WRITE_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN_V3); - WRITE_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3); - WRITE_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN_V4); - WRITE_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4); WRITE_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); @@ -2610,10 +1883,4 @@ case PS_STAT_FITTED_MEAN: case PS_STAT_FITTED_STDEV: - case PS_STAT_FITTED_MEAN_V2: - case PS_STAT_FITTED_STDEV_V2: - case PS_STAT_FITTED_MEAN_V3: - case PS_STAT_FITTED_STDEV_V3: - case PS_STAT_FITTED_MEAN_V4: - case PS_STAT_FITTED_STDEV_V4: case PS_STAT_CLIPPED_MEAN: case PS_STAT_CLIPPED_STDEV: @@ -2631,6 +1898,5 @@ { return options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | - PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2 | - PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_MEAN_V4); + PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN); } @@ -2638,6 +1904,5 @@ { return options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | - PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2 | PS_STAT_FITTED_STDEV_V3 | - PS_STAT_FITTED_STDEV_V4); + PS_STAT_FITTED_STDEV); } @@ -2665,16 +1930,4 @@ case PS_STAT_FITTED_STDEV: return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V2: - return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V2: - return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V3: - return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V3: - return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V4: - return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V4: - return stats->fittedStdev; case PS_STAT_CLIPPED_MEAN: return stats->clippedMean; @@ -3115,37 +2368,2 @@ return tmpFloat; } - -/****************************************************************************** -NOTE: We assume unnormalized gaussians. -*****************************************************************************/ -static psF32 minimizeLMChi2Gauss1D(psVector *deriv, - const psVector *params, - const psVector *coords - ) -{ - psTrace(TRACE, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN); - PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN); - PS_ASSERT_VECTOR_NON_NULL(coords, NAN); - PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN); - PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN); - - psF32 x = coords->data.F32[0]; - psF32 mean = params->data.F32[0]; - psF32 var = params->data.F32[1]; - psF32 dx = (x - mean); - - psF32 gauss = exp (-0.5*PS_SQR(dx)/var); - if (deriv) { - PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN); - PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN); - psF32 tmp = dx * gauss; - deriv->data.F32[0] = tmp / var; - deriv->data.F32[1] = tmp * dx / (var * var); - } - - - psTrace(TRACE, 4, "---- %s() end ----\n", __func__); - return gauss; -} Index: trunk/psLib/src/math/psStats.h =================================================================== --- trunk/psLib/src/math/psStats.h (revision 30995) +++ trunk/psLib/src/math/psStats.h (revision 31152) @@ -43,14 +43,8 @@ PS_STAT_FITTED_MEAN = 0x001000, ///< Fitted Mean PS_STAT_FITTED_STDEV = 0x002000, ///< Fitted Standard Deviation - PS_STAT_FITTED_MEAN_V2 = 0x004000, ///< Fitted Mean - PS_STAT_FITTED_STDEV_V2 = 0x008000, ///< Fitted Standard Deviation - PS_STAT_FITTED_MEAN_V3 = 0x010000, ///< Fitted Mean - PS_STAT_FITTED_STDEV_V3 = 0x020000, ///< Fitted Standard Deviation PS_STAT_CLIPPED_MEAN = 0x040000, ///< Clipped Mean PS_STAT_CLIPPED_STDEV = 0x080000, ///< Clipped Standard Deviation PS_STAT_USE_RANGE = 0x100000, ///< Range PS_STAT_USE_BINSIZE = 0x200000, ///< Binsize - PS_STAT_FITTED_MEAN_V4 = 0x400000, ///< Fitted Mean - PS_STAT_FITTED_STDEV_V4 = 0x800000, ///< Fitted Standard Deviation } psStatsOptions; Index: trunk/psLib/src/sys/psMemory.c =================================================================== --- trunk/psLib/src/sys/psMemory.c (revision 30995) +++ trunk/psLib/src/sys/psMemory.c (revision 31152) @@ -27,4 +27,5 @@ #include #include +#include #if defined(PS_MEM_BACKTRACE) && defined(HAVE_BACKTRACE) @@ -1090,2 +1091,35 @@ return (memBlock1->freeFunc == memBlock2->freeFunc); } + +bool static dumpMemory = false; + +void psMemDumpSetState (bool state) { + dumpMemory = state; +} + +void psMemDump(const char *name) +{ + if (!dumpMemory) return; + + char filename[1024]; // don't make your sub-names too long! + static int num = 0; // Counter, to make files unique and give an idea of sequence + + snprintf (filename, 1024, "memdump_%s_%03d.txt", name, num); + FILE *memFile = fopen(filename, "w"); + + psMemBlock **leaks = NULL; + int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true); + fprintf(memFile, "# MemBlock Size Source\n"); + unsigned long total = 0; // Total memory used + for (int i = 0; i < numLeaks; i++) { + psMemBlock *mb = leaks[i]; + fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, mb->file, mb->lineno); + total += mb->userMemorySize; + } + fclose(memFile); + psFree(leaks); + + // fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, (void *) sbrk(0)); + fprintf(stderr, "Memdump %s %d: Memory use: %ld\n", name, num, total); + num++; +} Index: trunk/psLib/src/sys/psMemory.h =================================================================== --- trunk/psLib/src/sys/psMemory.h (revision 30995) +++ trunk/psLib/src/sys/psMemory.h (revision 31152) @@ -647,4 +647,7 @@ ); +void psMemDumpSetState (bool state); +void psMemDump(const char *name); + // Ensure this is a psLib pointer #define PS_ASSERT_PTR_HEAVY(PTR, RVAL) \ Index: trunk/psLib/src/sys/psString.c =================================================================== --- trunk/psLib/src/sys/psString.c (revision 30995) +++ trunk/psLib/src/sys/psString.c (revision 31152) @@ -456,4 +456,17 @@ +psString psStringFileBasename (const char *fullname) { + + char *file; + + const char *ptr = strrchr (fullname, '/'); + if (ptr) { + file = psStringCopy(ptr + 1); + } else { + file = psStringCopy(fullname); + } + return (file); +} + psString psStringStripCVS(const char *string, const char *tagName) { Index: trunk/psLib/src/sys/psString.h =================================================================== --- trunk/psLib/src/sys/psString.h (revision 30995) +++ trunk/psLib/src/sys/psString.h (revision 31152) @@ -335,4 +335,6 @@ char *psStrcasestr (const char *haystack, const char *needle); +psString psStringFileBasename (const char *fullname); + /// @} #endif // #ifndef PS_STRING_H Index: trunk/psLib/src/types/psArguments.h =================================================================== --- trunk/psLib/src/types/psArguments.h (revision 30995) +++ trunk/psLib/src/types/psArguments.h (revision 31152) @@ -81,5 +81,5 @@ * specific routine called pkgnameVersionLong() is presumed to exist. */ -#define PSARGUMENTS_INSTANTIATE_GENERICS( pkgname, config, argc, argv ) \ +#define PS_ARGUMENTS_GENERIC( pkgname, config, argc, argv ) \ { int N= psArgumentGet (argc, argv, "-version"); \ if (N) { \ @@ -115,5 +115,5 @@ * presumed to exist. */ -#define PSARGUMENTS_INSTANTIATE_THREADSARG( pkgname, config, argc, argv ) \ +#define PS_ARGUMENTS_THREADS( pkgname, config, argc, argv ) \ { int N= psArgumentGet(argc, argv, "-threads"); \ if (N) { \ Index: trunk/psLib/test/math/tap_psStats_Sample_01.c =================================================================== --- trunk/psLib/test/math/tap_psStats_Sample_01.c (revision 30995) +++ trunk/psLib/test/math/tap_psStats_Sample_01.c (revision 31152) @@ -586,5 +586,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); @@ -622,5 +622,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); @@ -657,5 +657,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); @@ -694,5 +694,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); Index: trunk/psLib/test/optime/tap_psStatsTiming.c =================================================================== --- trunk/psLib/test/optime/tap_psStatsTiming.c (revision 30995) +++ trunk/psLib/test/optime/tap_psStatsTiming.c (revision 31152) @@ -678,5 +678,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32); for (int i = 0; i < rnd2->n; i++) @@ -702,5 +702,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32); for (int i = 0; i < rnd2->n; i++) @@ -725,5 +725,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32); for (int i = 0; i < rnd2->n; i++) @@ -790,5 +790,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *sample = psVectorAlloc (1000, PS_TYPE_F32); psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);