Index: /trunk/psLib/src/math/psStats.c
===================================================================
--- /trunk/psLib/src/math/psStats.c	(revision 11154)
+++ /trunk/psLib/src/math/psStats.c	(revision 11155)
@@ -13,6 +13,6 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.198 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-01-09 22:38:53 $
+ *  @version $Revision: 1.199 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-01-19 04:32:27 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -1132,6 +1132,6 @@
     float guessMean = stats->robustMedian;  // pass the guess mean
 
-    psTrace("psLib.math", 6, "The guess mean  is %f.\n", guessMean);
-    psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev);
+    psTrace("psLib.math", 6, "The ** starting ** guess mean  is %f.\n", guessMean);
+    psTrace("psLib.math", 6, "The ** starting ** guess stdev is %f.\n", guessStdev);
 
     bool done = false;
@@ -1258,6 +1258,16 @@
             psFree(y);
             psFree(poly);
-            // psFree(fitStats);
-            // psFree(fitMask);
+            psFree(histogram);
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        if (poly->coeff[2] >= 0.0) {
+            psTrace("psLib.math", 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+            psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
+            psFree(x);
+            psFree(y);
+            psFree(poly);
             psFree(histogram);
             psFree(statsMinMax);
@@ -1271,4 +1281,5 @@
             done = true;
         } else {
+            psTrace("psLib.math", 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
             psTrace("psLib.math", 6, "The new guess mean  is %f.\n", guessMean);
             psTrace("psLib.math", 6, "The new guess stdev is %f.\n", guessStdev);
@@ -1295,4 +1306,283 @@
     stats->results |= PS_STAT_FITTED_MEAN_V2;
     stats->results |= PS_STAT_FITTED_STDEV_V2;
+
+    return true;
+}
+
+/********************
+ * perform an asymmetric fit to the population
+ * vectorFittedStats_v3 requires guess for fittedMean and fittedStdev
+ * robustN50 should also be set
+ * gaussian fit is performed using 2D polynomial to ln(y)
+ ********************/
+static bool vectorFittedStats_v3 (const psVector* myVector,
+                                  const psVector* errors,
+                                  psVector* mask,
+                                  psMaskType maskVal,
+                                  psStats* stats)
+{
+
+    // This procedure requires the mean.  If it has not been already
+    // calculated, then call vectorSampleMean()
+    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
+        vectorRobustStats(myVector, errors, mask, maskVal, stats);
+    }
+
+    // If the mean is NAN, then generate a warning and set the stdev to NAN.
+    if (isnan(stats->robustMedian)) {
+        stats->fittedStdev = NAN;
+        stats->fittedStdev = NAN;
+        psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
+        return false;
+    }
+
+    float guessStdev = stats->robustStdev;  // pass the guess sigma
+    float guessMean = stats->robustMedian;  // pass the guess mean
+
+    psTrace("psLib.math", 6, "The ** starting ** guess mean  is %f.\n", guessMean);
+    psTrace("psLib.math", 6, "The ** starting ** guess stdev is %f.\n", guessStdev);
+
+    bool done = false;
+    for (int iteration = 0; !done && (iteration < 2); iteration ++) {
+        psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
+
+        psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
+        tmpScalar.type.type = PS_TYPE_F32;
+
+        psF32 binSize = 1;
+        if (stats->options & PS_STAT_USE_BINSIZE) {
+            // Set initial bin size to the specified value.
+            binSize = stats->binsize;
+            psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize);
+        } else {
+            // construct a histogram with (sigma/2 < binsize < sigma)
+            // set roughly so that the lowest bins have about 2 cnts
+            // Nsmallest ~ N50 / (4*dN))
+            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
+            binSize = guessStdev / dN;
+        }
+
+        // Determine the min/max of the vector (which prior outliers masked out)
+        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
+        float min = statsMinMax->min;
+        float max = statsMinMax->max;
+        if (numValid == 0 || isnan(min) || isnan(max)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(statsMinMax);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // Calculate the number of bins.
+        // XXX can we calculate the binMin, binMax **before** building this histogram?
+        long numBins = (max - min) / binSize;
+        psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max);
+        psTrace("psLib.math", 6, "The new bin size is %f.\n", binSize);
+        psTrace("psLib.math", 6, "The numBins is %ld\n", numBins);
+
+        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
+        histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
+        if (psTraceGetLevel("psLib.math") >= 8) {
+            PS_VECTOR_PRINT_F32(histogram->nums);
+        }
+
+        // now fit a Gaussian to the upper and lower halves about the peak independently
+
+        long binMin, binMax;
+        float upperMean, upperStdev;
+
+        // set the full-range upper and lower limits
+        psF32 maxFitSigma = 2.0;
+        if (isfinite(stats->clipSigma)) {
+            maxFitSigma = fabs(stats->clipSigma);
+        }
+        if (isfinite(stats->max)) {
+            maxFitSigma = fabs(stats->max);
+        }
+
+        psF32 minFitSigma = 2.0;
+        if (isfinite(stats->clipSigma)) {
+            minFitSigma = fabs(stats->clipSigma);
+        }
+        if (isfinite(stats->min)) {
+            minFitSigma = fabs(stats->min);
+        }
+
+        tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev;
+        if (tmpScalar.data.F32 < min) {
+            binMin = 0;
+        } else {
+            binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
+        }
+        tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev;
+        if (tmpScalar.data.F32 > max) {
+            binMax = histogram->bounds->n;
+        } else {
+            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar) + 1;
+        }
+        if (binMin < 0 || binMax < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma);
+            psFree(statsMinMax);
+            psFree(histogram);
+            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+            return false;
+        }
+
+        // search for mode (peak of histogram within range mean-2sigma - mean+2sigma
+        long  binPeak = binMin;
+        float valPeak = histogram->nums->data.F32[binPeak];
+        for (int i = binMin; i < binMax; i++) {
+            if (histogram->nums->data.F32[i] > valPeak) {
+                binPeak = i;
+                valPeak = histogram->nums->data.F32[binPeak];
+            }
+        }
+
+        // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak
+
+        psTrace("psLib.math", 6, "The clipped numBins is %ld\n", binMax - binMin);
+        psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin);
+        psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax);
+        psTrace("psLib.math", 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak);
+
+        {
+            // fit the lower half of the distribution
+            // run down until we drop below 0.25*valPeak
+            long binS = binMin;
+            long binE = PS_MIN (binPeak + 3, binMax);
+            for (int i = binPeak-3; i >= binMin; i--) {
+                if (histogram->nums->data.F32[i] < 0.25*valPeak) {
+                    binS = i;
+                    break;
+                }
+            }
+            psTrace("psLib.math", 6, "Lower bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS);
+            psTrace("psLib.math", 6, "Upper bound for lower half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE);
+
+            psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates
+            psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates
+            long j = 0;
+            for (long i = binS; i < binE; i++) {
+                if (histogram->nums->data.F32[i] <= 0.0)
+                    continue;
+                x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
+                y->data.F32[j] = log(histogram->nums->data.F32[i]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
+                j++;
+            }
+            y->n = x->n = j;
+
+            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
+            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+            psFree(x);
+            psFree(y);
+
+            if (!status) {
+                psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
+                psFree(poly);
+                psFree(histogram);
+                psFree(statsMinMax);
+                psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+                return false;
+            }
+
+            if (poly->coeff[2] >= 0.0) {
+                psTrace("psLib.math", 6, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+                psFree(poly);
+                psFree(histogram);
+                psFree(statsMinMax);
+
+                // sometimes, the guessStdev is much too large.  in this case, the entire real population tends to be found
+                // in a single bin.  make one attempt to recover by dropping the guessStdev down by a jump and trying again
+                if (iteration == 0) {
+                    guessStdev = 0.25*guessStdev;
+                    psTrace("psLib.math", 6, "*** retry stdev is %f.\n", guessStdev);
+                    continue;
+                }
+
+                psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
+                psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+                return false;
+            }
+
+            // calculate lower mean & stdev from parabolic fit -- use this as the result
+            guessStdev = sqrt(-0.5/poly->coeff[2]);
+            guessMean = poly->coeff[1]*PS_SQR(guessStdev);
+            if (guessStdev > 0.75*stats->robustStdev) {
+                done = true;
+            }
+            psTrace("psLib.math", 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+            psTrace("psLib.math", 6, "The lower mean  is %f.\n", guessMean);
+            psTrace("psLib.math", 6, "The lower stdev is %f.\n", guessStdev);
+
+            psFree(poly);
+        }
+
+        {
+            // fit the upper half of the distribution
+            // run up until we drop below 0.25*valPeak
+            long binS = PS_MAX (binPeak - 3, 0);
+            long binE = binMax;
+            for (int i = binPeak+3; i < binMax; i++) {
+                if (histogram->nums->data.F32[i] < 0.25*valPeak) {
+                    binE = i;
+                    break;
+                }
+            }
+            psTrace("psLib.math", 6, "Lower bound for upper half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binS), binS);
+            psTrace("psLib.math", 6, "Upper bound for upper half: %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binE), binE);
+
+            psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates
+            psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates
+            long j = 0;
+            for (long i = binS; i < binE; i++) {
+                if (histogram->nums->data.F32[i] <= 0.0)
+                    continue;
+                x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
+                y->data.F32[j] = log(histogram->nums->data.F32[i]); // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
+                j++;
+            }
+            y->n = x->n = j;
+
+            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
+            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
+            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+            psFree(x);
+            psFree(y);
+
+            if (!status) {
+                psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
+                psFree(poly);
+                psFree(histogram);
+                psFree(statsMinMax);
+                psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
+                return false;
+            }
+
+            // calculate upper mean & stdev from parabolic fit -- use this as the result
+            upperStdev = sqrt(-0.5/poly->coeff[2]);
+            upperMean = poly->coeff[1]*PS_SQR(upperStdev);
+            psTrace("psLib.math", 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+            psTrace("psLib.math", 6, "The upper mean  is %f.\n", upperMean);
+            psTrace("psLib.math", 6, "The upper stdev is %f.\n", upperStdev);
+
+            psFree (poly);
+        }
+
+        // Clean up after fitting
+        psFree (histogram);
+        psFree (statsMinMax);
+    }
+
+    // The fitted mean is the Gaussian mean.
+    stats->fittedMean = guessMean;
+    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);
+
+    // The fitted standard deviation
+    stats->fittedStdev = guessStdev;
+    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);
+
+    stats->results |= PS_STAT_FITTED_MEAN_V3;
+    stats->results |= PS_STAT_FITTED_STDEV_V3;
 
     return true;
@@ -1589,4 +1879,15 @@
 
     // ************************************************************************
+    if (stats->options & (PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_STDEV_V3)) {
+        if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
+            psAbort ("stats", "you may not specify both FITTED_MEAN and FITTED_MEAN_V3");
+        }
+        if (!vectorFittedStats_v3(inF32, errorsF32, maskU8, maskVal, stats)) {
+            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
+            status &= false;
+        }
+    }
+
+    // ************************************************************************
     if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
         if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
@@ -1628,4 +1929,7 @@
     READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
     READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
+    READ_STAT("FITTED_V3",       PS_STAT_FITTED_MEAN_V3);
+    READ_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN_V3);
+    READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3);
     READ_STAT("CLIPPED",         PS_STAT_CLIPPED_MEAN);
     READ_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
@@ -1657,4 +1961,6 @@
     WRITE_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
     WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
+    WRITE_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN_V3);
+    WRITE_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3);
     WRITE_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
     WRITE_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);
@@ -1708,4 +2014,6 @@
     case PS_STAT_FITTED_MEAN_V2:
     case PS_STAT_FITTED_STDEV_V2:
+    case PS_STAT_FITTED_MEAN_V3:
+    case PS_STAT_FITTED_STDEV_V3:
     case PS_STAT_CLIPPED_MEAN:
     case PS_STAT_CLIPPED_STDEV:
@@ -1742,4 +2050,8 @@
     case PS_STAT_FITTED_STDEV_V2:
         return stats->fittedStdev;
+    case PS_STAT_FITTED_MEAN_V3:
+        return stats->fittedMean;
+    case PS_STAT_FITTED_STDEV_V3:
+        return stats->fittedStdev;
     case PS_STAT_CLIPPED_MEAN:
         return stats->clippedMean;
Index: /trunk/psLib/src/math/psStats.h
===================================================================
--- /trunk/psLib/src/math/psStats.h	(revision 11154)
+++ /trunk/psLib/src/math/psStats.h	(revision 11155)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.57 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-12-08 11:38:54 $
+ *  @version $Revision: 1.58 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-01-19 04:32:27 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -50,4 +50,6 @@
     PS_STAT_FITTED_MEAN_V2  = 0x008000, ///< Fitted Mean
     PS_STAT_FITTED_STDEV_V2 = 0x010000, ///< Fitted Standard Deviation
+    PS_STAT_FITTED_MEAN_V3  = 0x020000, ///< Fitted Mean
+    PS_STAT_FITTED_STDEV_V3 = 0x040000, ///< Fitted Standard Deviation
 } psStatsOptions;