IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 31152 for trunk/psLib


Ignore:
Timestamp:
Apr 4, 2011, 12:57:08 PM (15 years ago)
Author:
eugene
Message:

consolidate multiple FITTED stats; updates to psImageBackground based on results from MWV; added memory dump API

Location:
trunk/psLib
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib

    • Property svn:ignore

      •  

        old new  
        3434*.bbg
        3535*.da
         36a.out.dSYM

  • trunk/psLib/src/fits/psFitsTable.c

    r28580 r31152  
    645645                    psMetadata *row = table->data[i]; // The row of interest
    646646                    psMetadataItem *dataItem = psMetadataLookup(row, colSpecItem->name); // Value of interest
    647                     memcpy(&columnData->data.U8[i * dataSize], &dataItem->data, dataSize);
     647                    if (dataItem) {
     648                        memcpy(&columnData->data.U8[i * dataSize], &dataItem->data, dataSize);
     649                    } else {
     650                        // this element is missing from this row; insert an appropriate-sized place holder
     651                        // XXX this should insert a NAN for float / double and an appropriate blank for int types
     652                        // XXX for the moment I am putting in 0.0
     653                        memset(&columnData->data.U8[i * dataSize], 0, dataSize);
     654                    }
    648655                }
    649656

  • trunk/psLib/src/imageops/psImageBackground.c

    r27069 r31152  
    1313#include "psType.h"
    1414#include "psAssert.h"
     15#include "psAbort.h"
    1516#include "psRandom.h"
    1617#include "psError.h"
    … …  
    2324}
    2425
    25 // XXX allow the user to choose the stats method?
    26 // (SAMPLE_MEAN, CLIPPED_MEAN, ROBUST_MEDIAN, FITTED_MEAN)
    2726bool psImageBackground(psStats *stats, psVector **sample, const psImage *image, const psImage *mask, psImageMaskType maskValue, psRandom *rng)
    2827{
    … …  
    4544    long nx = image->numCols;
    4645    long ny = image->numRows;
    47 
    48     psImage *bad = psImageAlloc(nx, ny, PS_TYPE_U8); // Image with bad pixels
    49     psImageInit(bad, 0);
    50 
    51     int Npixels = 0;                    // Total number of pixels
    52     for (int y = 0; y < ny; y++) {
    53         for (int x = 0; x < nx; x++) {
    54             if (!isfinite(image->data.F32[y][x]) ||
    55                 (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue)) {
    56                 bad->data.U8[y][x] = 0xFF;
    57             }
    58             Npixels++;
    59         }
    60     }
    61 
    62     const int Nsubset = (stats->nSubsample == 0) ? Npixels : PS_MIN(stats->nSubsample, Npixels); // Number of pixels in subset
    63 
    64     psVector *values;                   // Vector containing subsample
     46    long nPixels = nx * ny;
     47
     48    long nSubset = (stats->nSubsample == 0) ? (nx * ny) : PS_MIN(stats->nSubsample, (nx * ny));
     49
     50    /*** discussion of the sampling and algorithm choices:
     51
     52         We want stats->nSubsample pixels.  How many good pixels do we actually have? 
     53
     54         There are a few domains of interest:
     55
     56         If nGoodPixels < f1 * nSubset, we should fail [1]
     57         If nSubset >= nGoodPixels, we should just loop over all pixels [2]
     58         If (nSubset <  f2 * nGoodPixels) and (nGoodPixels >= f3 * nPixels), we should just use the random sample technique [3]
     59         If (nSubset >= f2 * nGoodPixels) or  (nGoodPixel  <  f3 * nPixels), we should use the Fisher-Yates technique.
     60
     61         to use Fisher-Yates, we need to have a copy of the pixels so we can re-shuffle.  We
     62         should not do that with the whole list unless we need to.
     63
     64         [1] we just have too few samples to be useful; f1 ~ 0.01?
     65
     66         [2] we need to select from all pixels to reach the desired sample size
     67
     68         [3] this avoids a full-image-sized alloc, but only makes sense if nGoodPixels is actually
     69         large.  f2 ~ 0.01, f3 ~ 0.1 (4 tries per success on average)
     70
     71    ***/
     72
     73    // count the number of good pixels
     74    long nGoodPixels = 0;
     75    for (long iy = 0; iy < ny; iy++) {
     76        for (long ix = 0; ix < nx; ix++) {
     77            if (!isfinite(image->data.F32[iy][ix])) continue;
     78            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
     79            nGoodPixels ++;
     80        }
     81    }
     82
     83    // XXX This value of 1% is somewhat arbitrary
     84    if (nGoodPixels < 0.01*nSubset) {
     85        if ((nFailures < 3) || (nFailures % 100 == 0)) {
     86            psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", nGoodPixels, nFailures);
     87        }
     88        nFailures ++;
     89        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);
     90        return false;
     91    }
     92
     93    // alloc or recycle the vector containing subsample pixels
     94    psVector *values;
    6595    if (sample) {
    66         *sample = psVectorRecycle(*sample, Nsubset, PS_TYPE_F32);
     96        *sample = psVectorRecycle(*sample, nSubset, PS_TYPE_F32);
    6797        values = psMemIncrRefCounter(*sample);
    6898        values->n = 0;
    6999    } else {
    70         values = psVectorAllocEmpty(Nsubset, PS_TYPE_F32);
     100        values = psVectorAllocEmpty(nSubset, PS_TYPE_F32);
    71101    }
    72102
    … …  
    75105    float max = -PS_MAX_F32;
    76106
    77     // select a subset of the image pixels to measure the stats
    78     long n = 0;                         // Number of actual pixels in subset
    79     if (Nsubset >= Npixels) {
    80         // if we have an image smaller than Nsubset, just loop over the (good) image pixels
    81         for (int iy = 0; iy < ny; iy++) {
    82             for (int ix = 0; ix < nx; ix++) {
    83                 if (bad->data.U8[iy][ix]) {
    84                     continue;
    85                 }
    86 
    87                 float value = image->data.F32[iy][ix];
    88                 min = PS_MIN(value, min);
    89                 max = PS_MAX(value, max);
    90                 values->data.F32[n] = value;
    91                 n++;
    92             }
    93         }
    94     } else {
    95         // Subsample all pixels
    96         // This is not optimal, since there may be a large masked fraction that leaves us with few good pixels.
    97         // In that case, you should have set Nsubset.......
    98         Npixels = nx * ny;
    99         for (long i = 0; i < Nsubset; i++) {
     107    if ((nSubset < 0.01*nGoodPixels) && (nGoodPixels >= 0.1*nPixels)) {
     108        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);
     109        // for small subsets, just use simple random sampling (saves a full-image-sized alloc),
     110        // unless we have lots of masked pixels
     111        for (long i = 0; i < nSubset; i++) {
    100112            double frnd = psRandomUniform(rng);
    101             int pixel = Npixels * frnd;
    102             int ix = pixel % nx;
    103             int iy = pixel / nx;
    104 
    105             if (bad->data.U8[iy][ix]) {
    106                 continue;
    107             }
     113            long pixel = nPixels * frnd;
     114            long ix = pixel % nx;
     115            long iy = pixel / nx;
     116
     117            if (!isfinite(image->data.F32[iy][ix])) continue;
     118            if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
    108119
    109120            float value = image->data.F32[iy][ix];
    110121            min = PS_MIN(value, min);
    111122            max = PS_MAX(value, max);
    112             values->data.F32[n] = value;
    113             n++;
    114         }
    115     }
    116 
    117     psFree(bad);
    118 
    119     if (n < 0.01*Nsubset) {
    120         if ((nFailures < 3) || (nFailures % 100 == 0)) {
    121             psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", n, nFailures);
    122         }
    123         nFailures ++;
    124         psFree(values);
    125         return false;
    126     }
    127 
    128     values->n = n;
     123            psVectorAppend (values, value);
     124        }
     125    } else if (nSubset >= nGoodPixels) {
     126            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);
     127        // in this case, we have to select from all masked pixels just to get the desired
     128        // sample size
     129        for (long iy = 0; iy < ny; iy++) {
     130            for (long ix = 0; ix < nx; ix++) {
     131                if (!isfinite(image->data.F32[iy][ix])) continue;
     132                if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
     133                float value = image->data.F32[iy][ix];
     134                min = PS_MIN(value, min);
     135                max = PS_MAX(value, max);
     136                psVectorAppend(values, value);
     137            }
     138        }
     139    } else {
     140        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);
     141        // use Knuth's version of Fisher-Yates to select a random subset of the pixels, but
     142        // only drawing each value once
     143           
     144        // generate a vector of all pixels which may in theory be selected
     145        psVector *pixelVector = psVectorAllocEmpty(nGoodPixels, PS_TYPE_F32);
     146        for (long iy = 0; iy < ny; iy++) {
     147            for (long ix = 0; ix < nx; ix++) {
     148                if (!isfinite(image->data.F32[iy][ix])) continue;
     149                if (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue) continue;
     150                psVectorAppend(pixelVector, image->data.F32[iy][ix]);
     151            }
     152        }
     153        psAssert (nGoodPixels == pixelVector->n, "we must have mis-counted above");
     154       
     155        // generate the unique random subset
     156        for (int i = 0; i < nSubset; i++) {
     157            double frnd = psRandomUniform(rng);
     158            int pixel = pixelVector->n * frnd;
     159               
     160            psVectorAppend(values, pixelVector->data.F32[pixel]);
     161            pixelVector->data.F32[pixel] = pixelVector->data.F32[pixelVector->n - 1];
     162            pixelVector->n --;
     163        }
     164        psFree (pixelVector);
     165    }
     166    psAssert (values->n >= 0.01*nSubset, "didn't we test this above?");
    129167
    130168    if (stats->options & PS_STAT_ROBUST_QUARTILE) {
    … …  
    132170        // XXX this hack is just for testing, drop when I am happy with the psStats version of the values
    133171
    134         int imin = stats->min * n;
    135         int imax = stats->max * n;
     172        int imin = stats->min * values->n;
     173        int imax = stats->max * values->n;
    136174        int npts = imax - imin + 1;
    137175
    … …  
    147185        // Subtract the median when we add the numbers, so we don't get numerical problems
    148186        float median = npts % 2 ? 0.5 * (values->data.F32[npts/2 - 1] + values->data.F32[npts/2]) :
    149                        values->data.F32[npts/2];
     187            values->data.F32[npts/2];
    150188        double value = 0;
    151         for (long i = imin; (i <= imax) && (i < n); i++) {
     189        for (long i = imin; (i <= imax) && (i < values->n); i++) {
    152190            value += values->data.F32[i] - median;
    153191        }
    … …  
    185223    return true;
    186224}
     225
     226
     227/***  old comments
     228
     229 // Subsample all pixels
     230 // This is not optimal, since there may be a large masked fraction that leaves us with few good pixels.
     231 // In that case, you should have set Nsubset.......
     232
     233 // 2011/03/22 - MWV:  On prompting from Gene, changed the sampling so that it doesn't build and then sort
     234 //     a huge array when Nsubset << Npixel. 
     235 //   O(Npixel) log(Npixel) is unnecessarily painful when Npixel=38 million (one skycell)
     236 //     but we only needed Nsubset pixels and Nsubset << Npixel.
     237
     238 //   It also perhaps partially defeats any gain of taking only subsamples if we're sorting the array.
     239 //   I should do some performance tests on this anyway.
     240
     241 // If our chance of collision is low and the gain from not sorting the array is likely to be high,
     242 // then just pick randomnly
     243 // Note that we're guaranteeing Nsubset pixels returned here (up to the limit of Npixels)
     244
     245 // 2011/03/16 - MWV:  Fixed double-sampling of sky pixels.
     246 //   The pick-a-random-pixel-at-a-time approach  doesn't work well when Npixels is close to Nsubset
     247 //   If Nsubset is 0.8*Npixels, then we will get lots of pixels double-counted
     248 //   This is correct on average, but isn't the optimal way to estimate the sky level
     249 // Instead we take the following approach:
     250 //   1) Calculate a random ordering of the pixels
     251 //   2) Go through this ordering up to Nsubset to select pixels
     252 // This is O(n log(n))
     253
     254 ***/

  • trunk/psLib/src/math/psMinimizePolyFit.c

    r24086 r31152  
    767767    // XXX enforce consistency?
    768768    // XXX psStatsGetValue() probably has inverted precedence
    769     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);
    770     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);
     769    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);
     770    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV);
    771771    if (!meanOption) {
    772772        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
    … …  
    12111211    // XXX enforce consistency?
    12121212    // XXX psStatsGetValue() probably has inverted precedence
    1213     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);
    1214     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);
     1213    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);
     1214    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV);
    12151215    if (!meanOption) {
    12161216        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
    … …  
    16211621    // XXX enforce consistency?
    16221622    // XXX psStatsGetValue() probably has inverted precedence
    1623     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);
    1624     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);
     1623    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);
     1624    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV);
    16251625    if (!meanOption) {
    16261626        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
    … …  
    20552055    // XXX enforce consistency?
    20562056    // XXX psStatsGetValue() probably has inverted precedence
    2057     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);
    2058     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);
     2057    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);
     2058    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV);
    20592059    if (!meanOption) {
    20602060        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");

  • trunk/psLib/src/math/psStats.c

    r30075 r31152  
    172172*****************************************************************************/
    173173
    174 // static prototypes:
    175 static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
    176 // static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    177 // static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    178174static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    179175
    … …  
    10871083}
    10881084
    1089 /*
    1090  * vectorFittedStats requires guess for fittedMean and fittedStdev
    1091  * robustN50 should also be set
    1092  */
     1085/********************
     1086 * perform an asymmetric fit to the population.  In development, this was called
     1087 * "vectorFittedStats_v4" all versions of fitted stats now resolve to this function (only v4
     1088 * has really been used) vectorFittedStats requires guess for fittedMean and fittedStdev
     1089 * robustN50 should also be set gaussian fit is performed using 2D polynomial to ln(y) this
     1090 * version follows the upper portion of the distribution until it passes 0.5*peak
     1091 ********************/
    10931092static bool vectorFittedStats (const psVector* myVector,
    1094                                const psVector* errors,
    1095                                psVector* mask,
    1096                                psVectorMaskType maskVal,
    1097                                psStats* stats)
     1093                                  const psVector* errors,
     1094                                  psVector* mask,
     1095                                  psVectorMaskType maskVal,
     1096                                  psStats* stats)
    10981097{
    10991098
    … …  
    11211120        stats->results |= PS_STAT_FITTED_MEAN;
    11221121        stats->results |= PS_STAT_FITTED_STDEV;
    1123         return true;
    1124     }
    1125 
    1126     float guessStdev = stats->robustStdev;  // pass the guess sigma
    1127     float guessMean = stats->robustMedian;  // pass the guess mean
    1128 
    1129     psTrace(TRACE, 6, "The guess mean  is %f.\n", guessMean);
    1130     psTrace(TRACE, 6, "The guess stdev is %f.\n", guessStdev);
    1131 
    1132     bool done = false;
    1133     for (int iteration = 0; !done && (iteration < 2); iteration ++) {
    1134         psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
    1135 
    1136         psF32 binSize = 1;
    1137         if (stats->options & PS_STAT_USE_BINSIZE) {
    1138             // Set initial bin size to the specified value.
    1139             binSize = stats->binsize;
    1140             psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize);
    1141         } else {
    1142             // construct a histogram with (sigma/10 < binsize < sigma)
    1143             // set roughly so that the lowest bins have about 2 cnts
    1144             // Nsmallest ~ N50 / (4*dN))
    1145             psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
    1146             binSize = guessStdev / dN;
    1147         }
    1148 
    1149         // Determine the min/max of the vector (which prior outliers masked out)
    1150         int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
    1151         float min = statsMinMax->min;
    1152         float max = statsMinMax->max;
    1153         if (numValid == 0 || isnan(min) || isnan(max)) {
    1154             psTrace(TRACE, 5, "Failed to calculate the min/max of the input vector.\n");
    1155             psFree(statsMinMax);
    1156             return true;
    1157         }
    1158 
    1159         // Calculate the number of bins.
    1160         // XXX can we calculate the binMin, binMax **before** building this histogram?
    1161         long numBins = (max - min) / binSize;
    1162         psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
    1163         psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
    1164         psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
    1165 
    1166         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    1167         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    1168             // if psVectorHistogram returns false, we have a programming error
    1169             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n");
    1170             psFree(histogram);
    1171             psFree(statsMinMax);
    1172             return false;
    1173         }
    1174         if (psTraceGetLevel("psLib.math") >= 8) {
    1175             PS_VECTOR_PRINT_F32(histogram->nums);
    1176         }
    1177 
    1178         // Fit a Gaussian to the bins in the requested range about the guess mean
    1179         // min,max overrides clipSigma
    1180         psF32 maxFitSigma = 2.0;
    1181         if (isfinite(stats->clipSigma)) {
    1182             maxFitSigma = fabs(stats->clipSigma);
    1183         }
    1184         if (isfinite(stats->max)) {
    1185             maxFitSigma = fabs(stats->max);
    1186         }
    1187 
    1188         psF32 minFitSigma = 2.0;
    1189         if (isfinite(stats->clipSigma)) {
    1190             minFitSigma = fabs(stats->clipSigma);
    1191         }
    1192         if (isfinite(stats->min)) {
    1193             minFitSigma = fabs(stats->min);
    1194         }
    1195 
    1196         // select the min and max bins, saturating on the lower and upper end-points
    1197         long binMin, binMax;
    1198         PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0);
    1199         PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    1200 
    1201         // Generate the variables that will be used in the Gaussian fitting
    1202         // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
    1203         psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
    1204         psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates
    1205         for (long i = binMin, j = 0; i <= binMax ; i++, j++) {
    1206             y->data.F32[j] = histogram->nums->data.F32[i];
    1207             psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value
    1208             ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i);
    1209             x->data[j] = ordinate;
    1210         }
    1211         if (psTraceGetLevel("psLib.math") >= 8) {
    1212             // XXX: Print the x array somehow.
    1213             PS_VECTOR_PRINT_F32(y);
    1214         }
    1215         psTrace(TRACE, 6, "The clipped numBins is %ld\n", y->n);
    1216         psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
    1217         psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
    1218 
    1219         // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
    1220         // XXX EAM : why not just fit the normalization as well??
    1221         p_psNormalizeVectorRange(y, 0.0, 1.0);
    1222 
    1223         // Fit a Gaussian to the data.
    1224         psMinimization *minimizer = psMinimizationAlloc(100, 0.01, 1.0); // The minimizer information
    1225         psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
    1226         // Initial guess for the mean (index 0) and var (index 1).
    1227         params->data.F32[0] = guessMean;
    1228         params->data.F32[1] = PS_SQR(guessStdev);
    1229         if (psTraceGetLevel("psLib.math") >= 8) {
    1230             PS_VECTOR_PRINT_F32(params);
    1231             PS_VECTOR_PRINT_F32(y);
    1232         }
    1233 
    1234         // psMinimizeLMChi2 can return false for bad data as well as for serious failures
    1235         if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
    1236             psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    1237             psFree(params);
    1238             psFree(minimizer);
    1239             psFree(x);
    1240             psFree(y);
    1241             psFree(histogram);
    1242             psFree(statsMinMax);
    1243             return true;
    1244         }
    1245         if (psTraceGetLevel("psLib.math") >= 8) {
    1246             PS_VECTOR_PRINT_F32(params);
    1247         }
    1248 
    1249         guessMean = params->data.F32[0];
    1250         guessStdev = sqrt(params->data.F32[1]);
    1251         if (guessStdev > 0.75*stats->robustStdev) {
    1252             done = true;
    1253         }
    1254 
    1255         // Clean up after fitting
    1256         psFree (params);
    1257         psFree (minimizer);
    1258         psFree (x);
    1259         psFree (y);
    1260         psFree (histogram);
    1261         psFree (statsMinMax);
    1262     }
    1263 
    1264     // The fitted mean is the Gaussian mean.
    1265     stats->fittedMean = guessMean;
    1266     psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean);
    1267 
    1268     // The fitted standard deviation
    1269     stats->fittedStdev = guessStdev;
    1270     psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
    1271 
    1272     stats->results |= PS_STAT_FITTED_MEAN;
    1273     stats->results |= PS_STAT_FITTED_STDEV;
    1274 
    1275     return true;
    1276 }
    1277 
    1278 
    1279 /********************
    1280  * vectorFittedStats_v2 requires guess for fittedMean and fittedStdev
    1281  * robustN50 should also be set
    1282  * gaussian fit is performed using 2D polynomial to ln(y)
    1283  ********************/
    1284 static bool vectorFittedStats_v2 (const psVector* myVector,
    1285                                   const psVector* errors,
    1286                                   psVector* mask,
    1287                                   psVectorMaskType maskVal,
    1288                                   psStats* stats)
    1289 {
    1290 
    1291     // This procedure requires the mean.  If it has not been already
    1292     // calculated, then call vectorSampleMean()
    1293     if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1294         if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
    1295             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
    1296             return false;
    1297         }
    1298     }
    1299 
    1300     // If the mean is NAN, then generate a warning and set the stdev to NAN.
    1301     if (isnan(stats->robustMedian)) {
    1302         stats->fittedMean = NAN;
    1303         stats->fittedStdev = NAN;
    1304         stats->results |= PS_STAT_FITTED_MEAN_V2;
    1305         stats->results |= PS_STAT_FITTED_STDEV_V2;
    1306         return true;
    1307     }
    1308 
    1309     if (stats->robustStdev <= FLT_EPSILON) {
    1310         stats->fittedMean = stats->robustMedian;
    1311         stats->fittedStdev = stats->robustStdev;
    1312         stats->results |= PS_STAT_FITTED_MEAN_V2;
    1313         stats->results |= PS_STAT_FITTED_STDEV_V2;
    1314         return true;
    1315     }
    1316 
    1317     float guessStdev = stats->robustStdev;  // pass the guess sigma
    1318     float guessMean = stats->robustMedian;  // pass the guess mean
    1319 
    1320     psTrace(TRACE, 6, "The ** starting ** guess mean  is %f.\n", guessMean);
    1321     psTrace(TRACE, 6, "The ** starting ** guess stdev is %f.\n", guessStdev);
    1322 
    1323     bool done = false;
    1324     for (int iteration = 0; !done && (iteration < 2); iteration ++) {
    1325         psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
    1326 
    1327         psF32 binSize = 1;
    1328         if (stats->options & PS_STAT_USE_BINSIZE) {
    1329             // Set initial bin size to the specified value.
    1330             binSize = stats->binsize;
    1331             psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize);
    1332         } else {
    1333             // construct a histogram with (sigma/10 < binsize < sigma)
    1334             // set roughly so that the lowest bins have about 2 cnts
    1335             // Nsmallest ~ N50 / (4*dN))
    1336             psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
    1337             binSize = guessStdev / dN;
    1338         }
    1339 
    1340         // Determine the min/max of the vector (which prior outliers masked out)
    1341         int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
    1342         float min = statsMinMax->min;
    1343         float max = statsMinMax->max;
    1344         if (numValid == 0 || isnan(min) || isnan(max)) {
    1345             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    1346             psFree(statsMinMax);
    1347             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1348             return true;
    1349         }
    1350 
    1351         // Calculate the number of bins.
    1352         // XXX can we calculate the binMin, binMax **before** building this histogram?
    1353         long numBins = (max - min) / binSize;
    1354         psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
    1355         psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
    1356         psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
    1357 
    1358         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    1359         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    1360             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistcs v2.\n");
    1361             psFree(histogram);
    1362             psFree(statsMinMax);
    1363             return false;
    1364         }
    1365         if (psTraceGetLevel("psLib.math") >= 8) {
    1366             PS_VECTOR_PRINT_F32(histogram->nums);
    1367         }
    1368 
    1369         // Fit a Gaussian to the bins in the requested range about the guess mean
    1370         // min,max overrides clipSigma
    1371         psF32 maxFitSigma = 2.0;
    1372         if (isfinite(stats->clipSigma)) {
    1373             maxFitSigma = fabs(stats->clipSigma);
    1374         }
    1375         if (isfinite(stats->max)) {
    1376             maxFitSigma = fabs(stats->max);
    1377         }
    1378 
    1379         psF32 minFitSigma = 2.0;
    1380         if (isfinite(stats->clipSigma)) {
    1381             minFitSigma = fabs(stats->clipSigma);
    1382         }
    1383         if (isfinite(stats->min)) {
    1384             minFitSigma = fabs(stats->min);
    1385         }
    1386 
    1387         // select the min and max bins, saturating on the lower and upper end-points
    1388         long binMin, binMax;
    1389         PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0);
    1390         PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    1391 
    1392         // Generate the variables that will be used in the Gaussian fitting
    1393         // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
    1394         psVector *y = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
    1395         psVector *x = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of ordinates
    1396         long j = 0;
    1397         for (long i = binMin; i <= binMax ; i++) {
    1398             if (histogram->nums->data.F32[i] <= 0.0)
    1399                 continue;
    1400             x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
    1401             y->data.F32[j] = log(histogram->nums->data.F32[i]);
    1402             // XXX note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
    1403             j++;
    1404         }
    1405         y->n = x->n = j;
    1406         if (psTraceGetLevel("psLib.math") >= 8) {
    1407             // XXX: Print the x array somehow.
    1408             PS_VECTOR_PRINT_F32(y);
    1409         }
    1410         psTrace(TRACE, 6, "The clipped numBins is %ld\n", y->n);
    1411         psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
    1412         psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
    1413 
    1414         // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
    1415         psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    1416 
    1417         // XXX how can we protect against some extreme outliers?  the robust histogram
    1418         // should have already dealt with those, but we are again sensitive to them...
    1419         // psStats *fitStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    1420         // fitStats->clipIter = 3.0;
    1421         // fitStats->clipSigma = 3.0;
    1422         // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_VECTOR_MASK);
    1423         // psVectorInit (fitMask, 0);
    1424 
    1425         // XXX not sure if these should result in errors or not...
    1426         if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
    1427             psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    1428             psFree(x);
    1429             psFree(y);
    1430             psFree(poly);
    1431             psFree(histogram);
    1432             psFree(statsMinMax);
    1433             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1434             return false;
    1435         }
    1436 
    1437         if (poly->coeff[2] >= 0.0) {
    1438             psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    1439             psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
    1440             psFree(x);
    1441             psFree(y);
    1442             psFree(poly);
    1443             psFree(histogram);
    1444             psFree(statsMinMax);
    1445             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1446             return false;
    1447         }
    1448 
    1449         guessStdev = sqrt(-0.5/poly->coeff[2]);
    1450         guessMean = poly->coeff[1]*PS_SQR(guessStdev);
    1451         if (guessStdev > 0.75*stats->robustStdev) {
    1452             done = true;
    1453         } else {
    1454             psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n",
    1455                     poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    1456             psTrace(TRACE, 6, "The new guess mean  is %f.\n", guessMean);
    1457             psTrace(TRACE, 6, "The new guess stdev is %f.\n", guessStdev);
    1458         }
    1459 
    1460         // Clean up after fitting
    1461         psFree (x);
    1462         psFree (y);
    1463         psFree (poly);
    1464         // psFree (fitStats);
    1465         // psFree (fitMask);
    1466         psFree (histogram);
    1467         psFree (statsMinMax);
    1468     }
    1469 
    1470     // The fitted mean is the Gaussian mean.
    1471     stats->fittedMean = guessMean;
    1472     psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean);
    1473 
    1474     // The fitted standard deviation
    1475     stats->fittedStdev = guessStdev;
    1476     psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
    1477 
    1478     stats->results |= PS_STAT_FITTED_MEAN_V2;
    1479     stats->results |= PS_STAT_FITTED_STDEV_V2;
    1480 
    1481     return true;
    1482 }
    1483 
    1484 /********************
    1485  * perform an asymmetric fit to the population
    1486  * vectorFittedStats_v3 requires guess for fittedMean and fittedStdev
    1487  * robustN50 should also be set
    1488  * gaussian fit is performed using 2D polynomial to ln(y)
    1489  ********************/
    1490 static bool vectorFittedStats_v3 (const psVector* myVector,
    1491                                   const psVector* errors,
    1492                                   psVector* mask,
    1493                                   psVectorMaskType maskVal,
    1494                                   psStats* stats)
    1495 {
    1496 
    1497     // This procedure requires the mean.  If it has not been already
    1498     // calculated, then call vectorSampleMean()
    1499     if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1500         if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
    1501             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
    1502             return false;
    1503         }
    1504     }
    1505 
    1506     // If the mean is NAN, then generate a warning and set the stdev to NAN.
    1507     if (isnan(stats->robustMedian)) {
    1508         stats->fittedMean = NAN;
    1509         stats->fittedStdev = NAN;
    1510         stats->results |= PS_STAT_FITTED_MEAN_V3;
    1511         stats->results |= PS_STAT_FITTED_STDEV_V3;
    1512         return true;
    1513     }
    1514 
    1515     if (stats->robustStdev <= FLT_EPSILON) {
    1516         stats->fittedMean = stats->robustMedian;
    1517         stats->fittedStdev = stats->robustStdev;
    1518         stats->results |= PS_STAT_FITTED_MEAN_V3;
    1519         stats->results |= PS_STAT_FITTED_STDEV_V3;
    15201122        return true;
    15211123    }
    … …  
    15511153            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    15521154            psFree(statsMinMax);
    1553             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1554             return true;
    1555         }
    1556 
    1557         // Calculate the number of bins.
    1558         // XXX can we calculate the binMin, binMax **before** building this histogram?
    1559         long numBins = (max - min) / binSize;
    1560         psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
    1561         psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
    1562         psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
    1563 
    1564         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    1565         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    1566             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics v3.\n");
    1567             psFree(histogram);
    1568             psFree(statsMinMax);
    1569             return false;
    1570         }
    1571         if (psTraceGetLevel("psLib.math") >= 8) {
    1572             PS_VECTOR_PRINT_F32(histogram->nums);
    1573         }
    1574 
    1575         // now fit a Gaussian to the upper and lower halves about the peak independently
    1576 
    1577         // set the full-range upper and lower limits
    1578         psF32 maxFitSigma = 2.0;
    1579         if (isfinite(stats->clipSigma)) {
    1580             maxFitSigma = fabs(stats->clipSigma);
    1581         }
    1582         if (isfinite(stats->max)) {
    1583             maxFitSigma = fabs(stats->max);
    1584         }
    1585 
    1586         psF32 minFitSigma = 2.0;
    1587         if (isfinite(stats->clipSigma)) {
    1588             minFitSigma = fabs(stats->clipSigma);
    1589         }
    1590         if (isfinite(stats->min)) {
    1591             minFitSigma = fabs(stats->min);
    1592         }
    1593 
    1594         // select the min and max bins, saturating on the lower and upper end-points
    1595         long binMin, binMax;
    1596         PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0);
    1597         PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    1598         if (binMin == binMax) {
    1599             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    1600             psFree(statsMinMax);
    1601             return true;
    1602         }
    1603 
    1604         // search for mode (peak of histogram within range mean-2sigma - mean+2sigma
    1605         long  binPeak = binMin;
    1606         float valPeak = histogram->nums->data.F32[binPeak];
    1607         for (int i = binMin; i < binMax; i++) {
    1608             if (histogram->nums->data.F32[i] > valPeak) {
    1609                 binPeak = i;
    1610                 valPeak = histogram->nums->data.F32[binPeak];
    1611             }
    1612             psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]);
    1613         }
    1614         psTrace (TRACE, 6, "\n");
    1615 
    1616         // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak
    1617 
    1618         psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin);
    1619         psTrace(TRACE, 6, "The clipped min is %f (%ld)\n",
    1620                 PS_BIN_MIDPOINT(histogram, binMin), binMin);
    1621         psTrace(TRACE, 6, "The clipped max is %f (%ld)\n",
    1622                 PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1);
    1623         psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n",
    1624                 PS_BIN_MIDPOINT(histogram, binPeak), binPeak);
    1625         psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]);
    1626 
    1627         {
    1628             // fit the lower half of the distribution
    1629             // run down until we drop below 0.25*valPeak
    1630             long binS = binMin;
    1631             long binE = PS_MIN (binPeak + 3, binMax);
    1632             for (int i = binPeak-3; i >= binMin; i--) {
    1633                 if (histogram->nums->data.F32[i] < 0.25*valPeak) {
    1634                     binS = i;
    1635                     break;
    1636                 }
    1637             }
    1638             psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n",
    1639                     PS_BIN_MIDPOINT(histogram, binS), binS);
    1640             psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n",
    1641                     PS_BIN_MIDPOINT(histogram, binE), binE);
    1642 
    1643             psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates
    1644             psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates
    1645             long j = 0;
    1646             for (long i = binS; i < binE; i++) {
    1647                 if (histogram->nums->data.F32[i] <= 0.0)
    1648                     continue;
    1649                 x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
    1650                 // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
    1651                 y->data.F32[j] = log(histogram->nums->data.F32[i]);
    1652                 j++;
    1653             }
    1654             y->n = x->n = j;
    1655 
    1656             // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
    1657             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    1658             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
    1659             psFree(x);
    1660             psFree(y);
    1661 
    1662             if (!status) {
    1663                 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    1664                 psFree(poly);
    1665                 psFree(histogram);
    1666                 psFree(statsMinMax);
    1667                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1668                 return false;
    1669             }
    1670 
    1671             if (poly->coeff[2] >= 0.0) {
    1672                 psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n",
    1673                         poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    1674                 psFree(poly);
    1675                 psFree(histogram);
    1676                 psFree(statsMinMax);
    1677 
    1678                 // sometimes, the guessStdev is much too large.  in this case, the entire real population
    1679                 // tends to be found in a single bin.  make one attempt to recover by dropping the guessStdev
    1680                 // down by a jump and trying again
    1681                 if (iteration == 0) {
    1682                     guessStdev = 0.25*guessStdev;
    1683                     psTrace(TRACE, 6, "*** retry, new stdev is %f.\n", guessStdev);
    1684                     continue;
    1685                 }
    1686 
    1687                 psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
    1688                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1689                 return false;
    1690             }
    1691 
    1692             // calculate lower mean & stdev from parabolic fit -- use this as the result
    1693             guessStdev = sqrt(-0.5/poly->coeff[2]);
    1694             guessMean = poly->coeff[1]*PS_SQR(guessStdev);
    1695             if (guessStdev > 0.75*stats->robustStdev) {
    1696                 done = true;
    1697             }
    1698             psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n",
    1699                     poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    1700             psTrace(TRACE, 6, "The lower mean  is %f.\n", guessMean);
    1701             psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev);
    1702 
    1703             psFree(poly);
    1704         }
    1705 
    1706         // for test, measure the same result for the upper section
    1707         {
    1708             // fit the upper half of the distribution
    1709             // run up until we drop below 0.25*valPeak
    1710             long binS = PS_MAX (binPeak - 3, 0);
    1711             long binE = binMax;
    1712             for (int i = binPeak+3; i < binMax; i++) {
    1713                 if (histogram->nums->data.F32[i] < 0.25*valPeak) {
    1714                     binE = i;
    1715                     break;
    1716                 }
    1717             }
    1718             psTrace(TRACE, 6, "Lower bound for upper half: %f (%ld)\n",
    1719                     PS_BIN_MIDPOINT(histogram, binS), binS);
    1720             psTrace(TRACE, 6, "Upper bound for upper half: %f (%ld)\n",
    1721                     PS_BIN_MIDPOINT(histogram, binE), binE);
    1722 
    1723             psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates
    1724             psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates
    1725             long j = 0;
    1726             for (long i = binS; i < binE; i++) {
    1727                 if (histogram->nums->data.F32[i] <= 0.0)
    1728                     continue;
    1729                 x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
    1730                 // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
    1731                 y->data.F32[j] = log(histogram->nums->data.F32[i]);
    1732                 j++;
    1733             }
    1734             y->n = x->n = j;
    1735 
    1736             // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
    1737             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    1738             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
    1739             psFree(x);
    1740             psFree(y);
    1741 
    1742             if (!status) {
    1743                 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    1744                 psFree(poly);
    1745                 psFree(histogram);
    1746                 psFree(statsMinMax);
    1747                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1748                 return false;
    1749             }
    1750 
    1751             // calculate upper mean & stdev from parabolic fit -- ignore this value
    1752             float upperStdev = sqrt(-0.5/poly->coeff[2]);
    1753             float upperMean = poly->coeff[1]*PS_SQR(upperStdev);
    1754 #ifndef PS_NO_TRACE
    1755             psTrace(TRACE, 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n",
    1756                     poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    1757             psTrace(TRACE, 6, "The upper mean  is %f.\n", upperMean);
    1758             psTrace(TRACE, 6, "The upper stdev is %f.\n", upperStdev);
    1759 #endif
    1760 
    1761             // if the resulting value is outside of the range binMin - binMax, use the upper value
    1762             if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) {
    1763                 guessMean = upperMean;
    1764                 guessStdev = upperStdev;
    1765             }
    1766 
    1767             psFree (poly);
    1768         }
    1769 
    1770         // Clean up after fitting
    1771         psFree (histogram);
    1772         psFree (statsMinMax);
    1773     }
    1774 
    1775     // The fitted mean is the Gaussian mean.
    1776     stats->fittedMean = guessMean;
    1777     psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean);
    1778 
    1779     // The fitted standard deviation
    1780     stats->fittedStdev = guessStdev;
    1781     psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
    1782 
    1783     stats->results |= PS_STAT_FITTED_MEAN_V3;
    1784     stats->results |= PS_STAT_FITTED_STDEV_V3;
    1785 
    1786     return true;
    1787 }
    1788 
    1789 /********************
    1790  * perform an asymmetric fit to the population
    1791  * vectorFittedStats_v4 requires guess for fittedMean and fittedStdev
    1792  * robustN50 should also be set
    1793  * gaussian fit is performed using 2D polynomial to ln(y)
    1794  * this version follows the upper portion of the distribution until it passes 0.5*peak
    1795  ********************/
    1796 static bool vectorFittedStats_v4 (const psVector* myVector,
    1797                                   const psVector* errors,
    1798                                   psVector* mask,
    1799                                   psVectorMaskType maskVal,
    1800                                   psStats* stats)
    1801 {
    1802 
    1803     // This procedure requires the mean.  If it has not been already
    1804     // calculated, then call vectorSampleMean()
    1805     if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1806         if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
    1807             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
    1808             return false;
    1809         }
    1810     }
    1811 
    1812     // If the mean is NAN, then generate a warning and set the stdev to NAN.
    1813     if (isnan(stats->robustMedian)) {
    1814         stats->fittedMean = NAN;
    1815         stats->fittedStdev = NAN;
    1816         stats->results |= PS_STAT_FITTED_MEAN_V4;
    1817         stats->results |= PS_STAT_FITTED_STDEV_V4;
    1818         return true;
    1819     }
    1820 
    1821     if (stats->robustStdev <= FLT_EPSILON) {
    1822         stats->fittedMean = stats->robustMedian;
    1823         stats->fittedStdev = stats->robustStdev;
    1824         stats->results |= PS_STAT_FITTED_MEAN_V4;
    1825         stats->results |= PS_STAT_FITTED_STDEV_V4;
    1826         return true;
    1827     }
    1828 
    1829     float guessStdev = stats->robustStdev;  // pass the guess sigma
    1830     float guessMean = stats->robustMedian;  // pass the guess mean
    1831 
    1832     psTrace(TRACE, 6, "The ** starting ** guess mean  is %f.\n", guessMean);
    1833     psTrace(TRACE, 6, "The ** starting ** guess stdev is %f.\n", guessStdev);
    1834 
    1835     bool done = false;
    1836     for (int iteration = 0; !done && (iteration < 2); iteration ++) {
    1837         psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
    1838 
    1839         psF32 binSize = 1;
    1840         if (stats->options & PS_STAT_USE_BINSIZE) {
    1841             // Set initial bin size to the specified value.
    1842             binSize = stats->binsize;
    1843             psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize);
    1844         } else {
    1845             // construct a histogram with (sigma/2 < binsize < sigma)
    1846             // set roughly so that the lowest bins have about 2 cnts
    1847             // Nsmallest ~ N50 / (4*dN))
    1848             psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
    1849             binSize = guessStdev / dN;
    1850         }
    1851 
    1852         // Determine the min/max of the vector (which prior outliers masked out)
    1853         int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
    1854         float min = statsMinMax->min;
    1855         float max = statsMinMax->max;
    1856         if (numValid == 0 || isnan(min) || isnan(max)) {
    1857             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    1858             psFree(statsMinMax);
    18591155            goto escape;
    18601156        }
    … …  
    18651161            stats->fittedMean = min;
    18661162            stats->fittedStdev = 0.0;
    1867             stats->results |= PS_STAT_FITTED_MEAN_V4;
    1868             stats->results |= PS_STAT_FITTED_STDEV_V4;
     1163            stats->results |= PS_STAT_FITTED_MEAN;
     1164            stats->results |= PS_STAT_FITTED_STDEV;
    18691165            return true;
    18701166        }
    … …  
    19331229
    19341230        // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak
    1935 
    19361231        psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin);
    19371232        psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
    1938         psTrace(TRACE, 6, "The clipped max is %f (%ld)\n",
    1939                 PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1);
     1233        psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1);
    19401234        psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak);
    19411235        psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]);
    19421236
     1237        float lowfitMean = NAN;
     1238        float lowfitStdev = NAN;
    19431239        {
    19441240            // fit the lower half of the distribution
    … …  
    20141310
    20151311            // calculate lower mean & stdev from parabolic fit -- use this as the result
    2016             guessStdev = sqrt(-0.5/poly->coeff[2]);
    2017             guessMean = poly->coeff[1]*PS_SQR(guessStdev);
    2018             if (guessStdev > 0.75*stats->robustStdev) {
    2019                 done = true;
    2020             }
    2021             psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n",
    2022                     poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    2023             psTrace(TRACE, 6, "The lower mean  is %f.\n", guessMean);
    2024             psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev);
     1312            lowfitStdev = sqrt(-0.5/poly->coeff[2]);
     1313            lowfitMean  = poly->coeff[1]*PS_SQR(lowfitStdev);
     1314
     1315            psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
     1316            psTrace(TRACE, 6, "The lower mean  is %f.\n", lowfitMean);
     1317            psTrace(TRACE, 6, "The lower stdev is %f.\n", lowfitStdev);
    20251318
    20261319            psFree(poly);
    20271320        }
    20281321
    2029         // if we converge on a solution outside the range binMin - binMax, use a more conservative range
    2030         float minValue = PS_BIN_MIDPOINT(histogram, binMin);
    2031         float maxValue = PS_BIN_MIDPOINT(histogram, binMax - 1);
    2032 
    2033         if (done && ((guessMean < minValue) || (guessMean > maxValue))) {
    2034             psTrace(TRACE, 6, "Inconsistent result, re-trying the fit\n");
    2035 
     1322        float fullfitMean  = NAN;
     1323        float fullfitStdev = NAN;
     1324        float minValueSym  = NAN;
     1325        float maxValueSym  = NAN;
     1326
     1327        // try the full fit as well:
     1328        {
    20361329            // fit a symmetric distribution
    20371330            // run up until we drop below 0.15*valPeak
    … …  
    20851378
    20861379            // calculate upper mean & stdev from parabolic fit -- ignore this value
    2087             guessStdev = sqrt(-0.5/poly->coeff[2]);
    2088             guessMean = poly->coeff[1]*PS_SQR(guessStdev);
     1380            fullfitStdev = sqrt(-0.5/poly->coeff[2]);
     1381            fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
    20891382#ifndef PS_NO_TRACE
    2090             psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n",
    2091                     poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    2092             psTrace(TRACE, 6, "The symmetric mean  is %f.\n", guessMean);
    2093             psTrace(TRACE, 6, "The symmetric stdev is %f.\n", guessStdev);
     1383            psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
     1384            psTrace(TRACE, 6, "The symmetric mean  is %f.\n", fullfitMean);
     1385            psTrace(TRACE, 6, "The symmetric stdev is %f.\n", fullfitStdev);
    20941386#endif
    20951387
    20961388            // if we converge on a solution outside the range binMin - binMax, use a more conservative range
    2097             float minValueSym = PS_BIN_MIDPOINT(histogram, binS);
    2098             float maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1);
     1389            minValueSym = PS_BIN_MIDPOINT(histogram, binS);
     1390            maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1);
    20991391
    21001392            // saturate on min or max value
    2101             if (guessMean < minValueSym) {
    2102                 guessMean = minValueSym;
     1393            if (fullfitMean < minValueSym) {
     1394                fullfitMean = minValueSym;
    21031395                psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
    21041396            }
    21051397
    21061398            // saturate on min or max value
    2107             if (guessMean > maxValueSym) {
    2108                 guessMean = maxValueSym;
     1399            if (fullfitMean > maxValueSym) {
     1400                fullfitMean = maxValueSym;
    21091401                psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
    21101402            }
    … …  
    21121404            psFree (poly);
    21131405        }
     1406
     1407        // we now have the fullfit and the lowfit mean and stdev values
     1408        // accept the fullfit unless minValueSym < lowfitMean < fullfitMean
     1409
     1410        if (isfinite(lowfitMean) && isfinite(lowfitStdev) && (lowfitMean < fullfitMean) && (lowfitMean > minValueSym)) {
     1411            guessMean  = lowfitMean;
     1412            guessStdev = lowfitStdev;
     1413        } else {
     1414            guessMean  = fullfitMean;
     1415            guessStdev = fullfitStdev;
     1416        }
     1417
     1418        if (!isfinite(guessMean) || !isfinite(guessStdev)) {
     1419            guessMean  = stats->robustMedian;
     1420            guessStdev = stats->robustStdev;
     1421        }
     1422
     1423        if (guessStdev > 0.75*stats->robustStdev) {
     1424            done = true;
     1425        }
    21141426
    21151427        // Clean up after fitting
    … …  
    21261438    psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
    21271439
    2128     stats->results |= PS_STAT_FITTED_MEAN_V4;
    2129     stats->results |= PS_STAT_FITTED_STDEV_V4;
     1440    stats->results |= PS_STAT_FITTED_MEAN;
     1441    stats->results |= PS_STAT_FITTED_STDEV;
    21301442
    21311443    return true;
    … …  
    21341446    stats->fittedMean = NAN;
    21351447    stats->fittedStdev = NAN;
    2136     stats->results |= PS_STAT_FITTED_MEAN_V4;
    2137     stats->results |= PS_STAT_FITTED_STDEV_V4;
     1448    stats->results |= PS_STAT_FITTED_MEAN;
     1449    stats->results |= PS_STAT_FITTED_STDEV;
    21381450
    21391451    return true;
    … …  
    24421754
    24431755    // ************************************************************************
    2444     if (stats->options & (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2)) {
    2445         if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
    2446             psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V2");
    2447         }
    2448         if (!vectorFittedStats_v2(inF32, errorsF32, maskVector, maskVal, stats)) {
    2449             psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
    2450             status &= false;
    2451         }
    2452     }
    2453 
    2454     // ************************************************************************
    2455     if (stats->options & (PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_STDEV_V3)) {
    2456         if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
    2457             psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V3");
    2458         }
    2459         if (!vectorFittedStats_v3(inF32, errorsF32, maskVector, maskVal, stats)) {
    2460             psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
    2461             status &= false;
    2462         }
    2463     }
    2464 
    2465     // ************************************************************************
    2466     if (stats->options & (PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4)) {
    2467         if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
    2468             psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V4");
    2469         }
    2470         if (!vectorFittedStats_v4(inF32, errorsF32, maskVector, maskVal, stats)) {
    2471             psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
    2472             status &= false;
    2473         }
    2474     }
    2475 
    2476     // ************************************************************************
    24771756    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
    24781757        if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) {
    … …  
    25131792    READ_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
    25141793    READ_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE);
    2515     READ_STAT("FITTED",         PS_STAT_FITTED_MEAN);
    2516     READ_STAT("FITTED_MEAN",    PS_STAT_FITTED_MEAN);
    2517     READ_STAT("FITTED_STDEV",   PS_STAT_FITTED_STDEV);
    2518     READ_STAT("FITTED_V2",       PS_STAT_FITTED_MEAN_V2);
    2519     READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
    2520     READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
    2521     READ_STAT("FITTED_V3",       PS_STAT_FITTED_MEAN_V3);
    2522     READ_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN_V3);
    2523     READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3);
    2524     READ_STAT("FITTED_V4",       PS_STAT_FITTED_MEAN_V4);
    2525     READ_STAT("FITTED_MEAN_V4",  PS_STAT_FITTED_MEAN_V4);
    2526     READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4);
     1794    READ_STAT("FITTED",          PS_STAT_FITTED_MEAN);
     1795    READ_STAT("FITTED_MEAN",     PS_STAT_FITTED_MEAN);
     1796    READ_STAT("FITTED_STDEV",    PS_STAT_FITTED_STDEV);
     1797    READ_STAT("FITTED_V2",       PS_STAT_FITTED_MEAN);
     1798    READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN);
     1799    READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV);
     1800    READ_STAT("FITTED_V3",       PS_STAT_FITTED_MEAN);
     1801    READ_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN);
     1802    READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV);
     1803    READ_STAT("FITTED_V4",       PS_STAT_FITTED_MEAN);
     1804    READ_STAT("FITTED_MEAN_V4",  PS_STAT_FITTED_MEAN);
     1805    READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV);
    25271806    READ_STAT("CLIPPED",         PS_STAT_CLIPPED_MEAN);
    25281807    READ_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
    … …  
    25541833    WRITE_STAT("FITTED_MEAN",     PS_STAT_FITTED_MEAN);
    25551834    WRITE_STAT("FITTED_STDEV",    PS_STAT_FITTED_STDEV);
    2556     WRITE_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
    2557     WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
    2558     WRITE_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN_V3);
    2559     WRITE_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3);
    2560     WRITE_STAT("FITTED_MEAN_V4",  PS_STAT_FITTED_MEAN_V4);
    2561     WRITE_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4);
    25621835    WRITE_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
    25631836    WRITE_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);
    … …  
    26101883      case PS_STAT_FITTED_MEAN:
    26111884      case PS_STAT_FITTED_STDEV:
    2612       case PS_STAT_FITTED_MEAN_V2:
    2613       case PS_STAT_FITTED_STDEV_V2:
    2614       case PS_STAT_FITTED_MEAN_V3:
    2615       case PS_STAT_FITTED_STDEV_V3:
    2616       case PS_STAT_FITTED_MEAN_V4:
    2617       case PS_STAT_FITTED_STDEV_V4:
    26181885      case PS_STAT_CLIPPED_MEAN:
    26191886      case PS_STAT_CLIPPED_STDEV:
    … …  
    26311898{
    26321899    return options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN |
    2633                       PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2 |
    2634                       PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_MEAN_V4);
     1900                      PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN);
    26351901}
    26361902
    … …  
    26381904{
    26391905    return options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV |
    2640                       PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2 | PS_STAT_FITTED_STDEV_V3 |
    2641                       PS_STAT_FITTED_STDEV_V4);
     1906                      PS_STAT_FITTED_STDEV);
    26421907}
    26431908
    … …  
    26651930      case PS_STAT_FITTED_STDEV:
    26661931        return stats->fittedStdev;
    2667       case PS_STAT_FITTED_MEAN_V2:
    2668         return stats->fittedMean;
    2669       case PS_STAT_FITTED_STDEV_V2:
    2670         return stats->fittedStdev;
    2671       case PS_STAT_FITTED_MEAN_V3:
    2672         return stats->fittedMean;
    2673       case PS_STAT_FITTED_STDEV_V3:
    2674         return stats->fittedStdev;
    2675       case PS_STAT_FITTED_MEAN_V4:
    2676         return stats->fittedMean;
    2677       case PS_STAT_FITTED_STDEV_V4:
    2678         return stats->fittedStdev;
    26791932      case PS_STAT_CLIPPED_MEAN:
    26801933        return stats->clippedMean;
    … …  
    31152368    return tmpFloat;
    31162369}
    3117 
    3118 /******************************************************************************
    3119 NOTE: We assume unnormalized gaussians.
    3120 *****************************************************************************/
    3121 static psF32 minimizeLMChi2Gauss1D(psVector *deriv,
    3122                                    const psVector *params,
    3123                                    const psVector *coords
    3124     )
    3125 {
    3126     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    3127     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    3128     PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN);
    3129     PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN);
    3130     PS_ASSERT_VECTOR_NON_NULL(coords, NAN);
    3131     PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN);
    3132     PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN);
    3133 
    3134     psF32 x = coords->data.F32[0];
    3135     psF32 mean = params->data.F32[0];
    3136     psF32 var = params->data.F32[1];
    3137     psF32 dx = (x - mean);
    3138 
    3139     psF32 gauss = exp (-0.5*PS_SQR(dx)/var);
    3140     if (deriv) {
    3141         PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN);
    3142         PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN);
    3143         psF32 tmp = dx * gauss;
    3144         deriv->data.F32[0] = tmp / var;
    3145         deriv->data.F32[1] = tmp * dx / (var * var);
    3146     }
    3147 
    3148 
    3149     psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    3150     return gauss;
    3151 }

  • trunk/psLib/src/math/psStats.h

    r21183 r31152  
    4343    PS_STAT_FITTED_MEAN     = 0x001000, ///< Fitted Mean
    4444    PS_STAT_FITTED_STDEV    = 0x002000, ///< Fitted Standard Deviation
    45     PS_STAT_FITTED_MEAN_V2  = 0x004000, ///< Fitted Mean
    46     PS_STAT_FITTED_STDEV_V2 = 0x008000, ///< Fitted Standard Deviation
    47     PS_STAT_FITTED_MEAN_V3  = 0x010000, ///< Fitted Mean
    48     PS_STAT_FITTED_STDEV_V3 = 0x020000, ///< Fitted Standard Deviation
    4945    PS_STAT_CLIPPED_MEAN    = 0x040000, ///< Clipped Mean
    5046    PS_STAT_CLIPPED_STDEV   = 0x080000, ///< Clipped Standard Deviation
    5147    PS_STAT_USE_RANGE       = 0x100000, ///< Range
    5248    PS_STAT_USE_BINSIZE     = 0x200000, ///< Binsize
    53     PS_STAT_FITTED_MEAN_V4  = 0x400000, ///< Fitted Mean
    54     PS_STAT_FITTED_STDEV_V4 = 0x800000, ///< Fitted Standard Deviation
    5549} psStatsOptions;
    5650

  • trunk/psLib/src/sys/psMemory.c

    r30595 r31152  
    2727#include <string.h>
    2828#include <assert.h>
     29#include <unistd.h>
    2930
    3031#if defined(PS_MEM_BACKTRACE) && defined(HAVE_BACKTRACE)
    … …  
    10901091    return (memBlock1->freeFunc == memBlock2->freeFunc);
    10911092}
     1093
     1094bool static dumpMemory = false;
     1095
     1096void psMemDumpSetState (bool state) {
     1097    dumpMemory = state;
     1098}
     1099
     1100void psMemDump(const char *name)
     1101{
     1102    if (!dumpMemory) return;
     1103
     1104    char filename[1024];          // don't make your sub-names too long!
     1105    static int num = 0;           // Counter, to make files unique and give an idea of sequence
     1106
     1107    snprintf (filename, 1024, "memdump_%s_%03d.txt", name, num);
     1108    FILE *memFile = fopen(filename, "w");
     1109
     1110    psMemBlock **leaks = NULL;
     1111    int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true);
     1112    fprintf(memFile, "# MemBlock Size Source\n");
     1113    unsigned long total = 0;            // Total memory used
     1114    for (int i = 0; i < numLeaks; i++) {
     1115        psMemBlock *mb = leaks[i];
     1116        fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, mb->file, mb->lineno);
     1117        total += mb->userMemorySize;
     1118    }
     1119    fclose(memFile);
     1120    psFree(leaks);
     1121
     1122    // fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, (void *) sbrk(0));
     1123    fprintf(stderr, "Memdump %s %d: Memory use: %ld\n", name, num, total);
     1124    num++;
     1125}

  • trunk/psLib/src/sys/psMemory.h

    r26892 r31152  
    647647  );
    648648
     649void psMemDumpSetState (bool state);
     650void psMemDump(const char *name);
     651
    649652// Ensure this is a psLib pointer
    650653#define PS_ASSERT_PTR_HEAVY(PTR, RVAL) \

  • trunk/psLib/src/sys/psString.c

    r29932 r31152  
    456456
    457457
     458psString psStringFileBasename (const char *fullname) {
     459 
     460    char *file;
     461
     462    const char *ptr = strrchr (fullname, '/');
     463    if (ptr) {
     464        file = psStringCopy(ptr + 1);
     465    } else {
     466        file = psStringCopy(fullname);
     467    }
     468  return (file);
     469}
     470
    458471psString psStringStripCVS(const char *string, const char *tagName)
    459472{

  • trunk/psLib/src/sys/psString.h

    r29932 r31152  
    335335char *psStrcasestr (const char *haystack, const char *needle);
    336336
     337psString psStringFileBasename (const char *fullname);
     338
    337339/// @}
    338340#endif // #ifndef PS_STRING_H

  • trunk/psLib/src/types/psArguments.h

    r24143 r31152  
    8181 *  specific routine called pkgnameVersionLong() is presumed to exist.
    8282 */
    83 #define PSARGUMENTS_INSTANTIATE_GENERICS( pkgname, config, argc, argv )   \
     83#define PS_ARGUMENTS_GENERIC( pkgname, config, argc, argv )   \
    8484  { int N= psArgumentGet (argc, argv, "-version");                        \
    8585    if (N) {                                                              \
    … …  
    115115 *  presumed to exist.
    116116 */
    117 #define PSARGUMENTS_INSTANTIATE_THREADSARG( pkgname, config, argc, argv )   \
     117#define PS_ARGUMENTS_THREADS( pkgname, config, argc, argv )   \
    118118  { int N= psArgumentGet(argc, argv, "-threads");                           \
    119119    if (N) {                                                                \

  • trunk/psLib/test/math/tap_psStats_Sample_01.c

    r30077 r31152  
    586586        psFree (stats);
    587587
    588         stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE);
     588        stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE);
    589589        stats->binsize = 1.0;
    590590        psVectorStats (stats, y, NULL, NULL, 1);
    … …  
    622622        psFree (stats);
    623623
    624         stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE);
     624        stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE);
    625625        stats->binsize = 1.0;
    626626        psVectorStats (stats, y, NULL, NULL, 1);
    … …  
    657657        psFree (stats);
    658658
    659         stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE);
     659        stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE);
    660660        stats->binsize = 1.0;
    661661        psVectorStats (stats, y, NULL, NULL, 1);
    … …  
    694694        psFree (stats);
    695695
    696         stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE);
     696        stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE);
    697697        stats->binsize = 1.0;
    698698        psVectorStats (stats, y, NULL, NULL, 1);

  • trunk/psLib/test/optime/tap_psStatsTiming.c

    r24572 r31152  
    678678        psMemId id = psMemGetId();
    679679
    680         psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
     680        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
    681681        psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32);
    682682        for (int i = 0; i < rnd2->n; i++)
    … …  
    702702        psMemId id = psMemGetId();
    703703
    704         psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
     704        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
    705705        psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32);
    706706        for (int i = 0; i < rnd2->n; i++)
    … …  
    725725        psMemId id = psMemGetId();
    726726
    727         psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2);
     727        psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV);
    728728        psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32);
    729729        for (int i = 0; i < rnd2->n; i++)
    … …  
    790790        psMemId id = psMemGetId();
    791791
    792         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);
     792        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);
    793793        psVector *sample = psVectorAlloc (1000, PS_TYPE_F32);
    794794        psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);
Note: See TracChangeset for help on using the changeset viewer.