IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30863


Ignore:
Timestamp:
Mar 10, 2011, 5:18:40 PM (15 years ago)
Author:
eugene
Message:

once more address issues with the fitted mean measurement: if the lower-half fit is poor, it can give a result higher than the full sample fit; only accept the lower-half fit if it is lower than the full fit and higher than the lower range of the lower-half data sample; deprecate all of the old named versions of the fitted mean (eg, FITTED_MEAN_V2,V3, etc). All of these names now select the single best-tested version of the algorithm

Location:
branches/eam_branches/ipp-20110213/psLib
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110213/psLib/src/math/psMinimizePolyFit.c

    r24086 r30863  
    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");

  • branches/eam_branches/ipp-20110213/psLib/src/math/psStats.c

    r30075 r30863  
    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        }
    … …  
    19341230        // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak
    19351231
     1232        float clippedMean = PS_BIN_MIDPOINT(histogram, binPeak);
     1233
    19361234        psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin);
    19371235        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);
    1940         psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak);
     1236        psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1);
     1237        psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", clippedMean, binPeak);
    19411238        psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]);
    19421239
     1240        float lowfitMean = NAN;
     1241        float lowfitStdev = NAN;
    19431242        {
    19441243            // fit the lower half of the distribution
    … …  
    20141313
    20151314            // 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);
     1315            lowfitStdev = sqrt(-0.5/poly->coeff[2]);
     1316            lowfitMean  = poly->coeff[1]*PS_SQR(lowfitStdev);
     1317
     1318            psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
     1319            psTrace(TRACE, 6, "The lower mean  is %f.\n", lowfitMean);
     1320            psTrace(TRACE, 6, "The lower stdev is %f.\n", lowfitStdev);
    20251321
    20261322            psFree(poly);
    20271323        }
    20281324
    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 
     1325        float fullfitMean  = NAN;
     1326        float fullfitStdev = NAN;
     1327        float minValueSym  = NAN;
     1328        float maxValueSym  = NAN;
     1329
     1330        // try the full fit as well:
     1331        {
    20361332            // fit a symmetric distribution
    20371333            // run up until we drop below 0.15*valPeak
    … …  
    20851381
    20861382            // 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);
     1383            fullfitStdev = sqrt(-0.5/poly->coeff[2]);
     1384            fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
    20891385#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);
     1386            psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
     1387            psTrace(TRACE, 6, "The symmetric mean  is %f.\n", fullfitMean);
     1388            psTrace(TRACE, 6, "The symmetric stdev is %f.\n", fullfitStdev);
    20941389#endif
    20951390
    20961391            // 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);
     1392            minValueSym = PS_BIN_MIDPOINT(histogram, binS);
     1393            maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1);
    20991394
    21001395            // saturate on min or max value
    2101             if (guessMean < minValueSym) {
    2102                 guessMean = minValueSym;
     1396            if (fullfitMean < minValueSym) {
     1397                fullfitMean = minValueSym;
    21031398                psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
    21041399            }
    21051400
    21061401            // saturate on min or max value
    2107             if (guessMean > maxValueSym) {
    2108                 guessMean = maxValueSym;
     1402            if (fullfitMean > maxValueSym) {
     1403                fullfitMean = maxValueSym;
    21091404                psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean);
    21101405            }
    … …  
    21121407            psFree (poly);
    21131408        }
     1409
     1410        // we now have the fullfit and the lowfit mean and stdev values
     1411        // accept the fullfit unless minValueSym < lowfitMean < fullfitMean
     1412
     1413        if (isfinite(lowfitMean) && isfinite(lowfitStdev) && (lowfitMean < fullfitMean) && (lowfitMean > minValueSym)) {
     1414            guessMean  = lowfitMean;
     1415            guessStdev = lowfitStdev;
     1416        } else {
     1417            guessMean  = fullfitMean;
     1418            guessStdev = fullfitStdev;
     1419        }
     1420
     1421        if (!isfinite(guessMean) || !isfinite(guessStdev)) {
     1422            guessMean  = stats->robustMedian;
     1423            guessStdev = stats->robustStdev;
     1424        }
     1425
     1426        if (guessStdev > 0.75*stats->robustStdev) {
     1427            done = true;
     1428        }
    21141429
    21151430        // Clean up after fitting
    … …  
    21261441    psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
    21271442
    2128     stats->results |= PS_STAT_FITTED_MEAN_V4;
    2129     stats->results |= PS_STAT_FITTED_STDEV_V4;
     1443    stats->results |= PS_STAT_FITTED_MEAN;
     1444    stats->results |= PS_STAT_FITTED_STDEV;
    21301445
    21311446    return true;
    … …  
    21341449    stats->fittedMean = NAN;
    21351450    stats->fittedStdev = NAN;
    2136     stats->results |= PS_STAT_FITTED_MEAN_V4;
    2137     stats->results |= PS_STAT_FITTED_STDEV_V4;
     1451    stats->results |= PS_STAT_FITTED_MEAN;
     1452    stats->results |= PS_STAT_FITTED_STDEV;
    21381453
    21391454    return true;
    … …  
    24421757
    24431758    // ************************************************************************
    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     // ************************************************************************
    24771759    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
    24781760        if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) {
    … …  
    25131795    READ_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
    25141796    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);
     1797    READ_STAT("FITTED",          PS_STAT_FITTED_MEAN);
     1798    READ_STAT("FITTED_MEAN",     PS_STAT_FITTED_MEAN);
     1799    READ_STAT("FITTED_STDEV",    PS_STAT_FITTED_STDEV);
     1800    READ_STAT("FITTED_V2",       PS_STAT_FITTED_MEAN);
     1801    READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN);
     1802    READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV);
     1803    READ_STAT("FITTED_V3",       PS_STAT_FITTED_MEAN);
     1804    READ_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN);
     1805    READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV);
     1806    READ_STAT("FITTED_V4",       PS_STAT_FITTED_MEAN);
     1807    READ_STAT("FITTED_MEAN_V4",  PS_STAT_FITTED_MEAN);
     1808    READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV);
    25271809    READ_STAT("CLIPPED",         PS_STAT_CLIPPED_MEAN);
    25281810    READ_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
    … …  
    25541836    WRITE_STAT("FITTED_MEAN",     PS_STAT_FITTED_MEAN);
    25551837    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);
    25621838    WRITE_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
    25631839    WRITE_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);
    … …  
    26101886      case PS_STAT_FITTED_MEAN:
    26111887      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:
    26181888      case PS_STAT_CLIPPED_MEAN:
    26191889      case PS_STAT_CLIPPED_STDEV:
    … …  
    26311901{
    26321902    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);
     1903                      PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN);
    26351904}
    26361905
    … …  
    26381907{
    26391908    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);
     1909                      PS_STAT_FITTED_STDEV);
    26421910}
    26431911
    … …  
    26651933      case PS_STAT_FITTED_STDEV:
    26661934        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;
    26791935      case PS_STAT_CLIPPED_MEAN:
    26801936        return stats->clippedMean;
    … …  
    31152371    return tmpFloat;
    31162372}
    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 }

  • branches/eam_branches/ipp-20110213/psLib/src/math/psStats.h

    r21183 r30863  
    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

  • branches/eam_branches/ipp-20110213/psLib/test/math/tap_psStats_Sample_01.c

    r30077 r30863  
    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);

  • branches/eam_branches/ipp-20110213/psLib/test/optime/tap_psStatsTiming.c

    r24572 r30863  
    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.