/** @file  psStats.c
 *  \brief basic statistical operations
 *  @ingroup Stat
 *
 *  This file will hold the definitions of the histogram and stats data
 *  structures.  It also contains prototypes for procedures which operate
 *  on those data structures.
 *
 *  @author GLG (MHPCC), EAM (IfA)
 *
 * XXX: Must do
 * nSubsample points
 * use ->min and ->max (PS_STAT_USE_RANGE)
 *
 *  @version $Revision: 1.199 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2007-01-19 04:32:27 $
 *
 *  Copyright 2006 IfA, University of Hawaii
 */

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <stdlib.h>
#include <stdio.h>
#include <stdarg.h>
#include <float.h>
#include <math.h>
#include <limits.h>
#include <strings.h>

/*****************************************************************************/
/* INCLUDE FILES                                                             */
/*****************************************************************************/
#include "psMemory.h"
#include "psAbort.h"
//#include "psImage.h"
#include "psVector.h"
#include "psTrace.h"
#include "psLogMsg.h"
#include "psError.h"
#include "psHistogram.h"
#include "psStats.h"
#include "psMinimizeLMM.h"
#include "psMinimizePolyFit.h"
#include "psPolynomial.h"
#include "psAssert.h"
#include "psMathUtils.h"
#include "psList.h"
#include "psString.h"



/*****************************************************************************/
/* DEFINE STATEMENTS                                                         */
/*****************************************************************************/
#define PS_GAUSS_WIDTH 3                // The width of the Gaussian smoothing.
#define PS_CLIPPED_NUM_ITER_LB 1 // This corresponds to N in the ADD.
#define PS_CLIPPED_NUM_ITER_UB 10
#define PS_CLIPPED_SIGMA_LB 1.0
#define PS_CLIPPED_SIGMA_UB 10.0
#define PS_POLY_MEDIAN_MAX_ITERATIONS 30
#define MASK_MARK 0x80   // bit to use internally to mark data as bad

#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))

/*****************************************************************************/
/* TYPE DEFINITIONS                                                          */
/*****************************************************************************/

// None

/*****************************************************************************/
/* GLOBAL VARIABLES                                                          */
/*****************************************************************************/

// None

/*****************************************************************************/
/* FILE STATIC VARIABLES                                                     */
/*****************************************************************************/

// None

/******************************************************************************
MISC PRIVATE STATISTICAL FUNCTIONS
 
NOTE: it is assumed that any call to these statistical functions will have
been preceded by a call to the psVectorStats() function.  Various sanity tests
will only be performed in psVectorStats().
 
For many of these private stats routines, it is possible that there are no acceptable elements
in the input vector (if no elements lie within range, or there are no unmasked elements, or the
input vector is NULL).  In such cases, we set the desired value to NAN and return an error.
The calling function may clear the error.
*****************************************************************************/

// static prototypes:
static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);

/******************************************************************************
vectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
mean of the input vector.  If there was a problem with the mean calculation,
this routine sets stats->sampleMean to NAN.
 
using the method below, with a single loop for various options costs only a small amount and is
much easier to debug. running 10000 tests of 1000 point vectors for the two methods gives:
 
                     separate    single          w/errors
(mask: 0, range: 0): 0.067 sec  0.073 sec (10%)  0.072 sec
(mask: 1, range: 0): 0.098 sec  0.102 sec ( 4%)  0.119 sec
(mask: 0, range: 0): 0.136 sec  0.162 sec (20%)  0.170 sec
(mask: 1, range: 0): 0.177 sec  0.181 sec ( 2%)  0.198 sec
 
(effect of errors not tested)
 
To optmize this, use a macro and ifdef in or out the three states (errors, mask, range)
*****************************************************************************/
static bool vectorSampleMean(const psVector* myVector,
                             const psVector* errors,
                             const psVector* maskVector,
                             psMaskType maskVal,
                             psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);

    long count = 0;                     // Number of points contributing to this mean
    psF32 mean = 0.0;                   // The mean
    psF32 weight;

    psF32 *data = myVector->data.F32;   // Dereference
    int numData = myVector->n;          // Number of data points

    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
    bool useRange = stats->options & PS_STAT_USE_RANGE;

    psF32 sumWeights = 0.0;  // The sum of the weights
    psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32;

    for (long i = 0; i < numData; i++) {
        // Check if the data is with the specified range
        if (useRange && (data[i] < stats->min))
            continue;
        if (useRange && (data[i] > stats->max))
            continue;
        if (maskData && (maskData[i] & maskVal))
            continue;
        if (errors) {
            weight = (errorsData[i] == 0) ? 0.0 : PS_SQR(errorsData[i]);
            mean += data[i] * weight;
            sumWeights += weight;
        } else {
            mean += data[i];
        }
        count++;
    }
    if (errors) {
        mean = (count > 0) ? mean / sumWeights : NAN;
    } else {
        mean = (count > 0) ? mean / count : NAN;
    }
    stats->sampleMean = mean;
    if (isnan(mean)) {
        // XXX raise an error here?
        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
        return false;
    }

    stats->results |= PS_STAT_SAMPLE_MEAN;
    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
    return true;
}

/******************************************************************************
vectorMinMax(myVector, maskVector, maskVal, stats): calculates the minimum and maximum of the input vector.
If there was a problem with the calculation, this routine sets stats->max and stats->min to NAN.  Return the
number of valid values in the vector (those not masked or outside the specified range.
XXX : using the method below, with a single loop for various options
      costs only a small amount and is much easier to debug. running 10000 tests
      of 1000 point vectors for the two methods gives:
                     separate    single
(mask: 0, range: 0):  0.101 sec  0.149 sec
(mask: 1, range: 0):  0.125 sec  0.160 sec
(mask: 0, range: 0):  0.179 sec  0.208 sec
(mask: 1, range: 0):  0.200 sec  0.244 sec
*****************************************************************************/
static long vectorMinMax(const psVector* myVector,
                         const psVector* maskVector,
                         psMaskType maskVal,
                         psStats* stats
                        )
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    psF32 max = -PS_MAX_F32;            // The calculated maximum
    psF32 min = PS_MAX_F32;             // The calculated minimum
    psF32 *vector = myVector->data.F32; // Dereference the vector

    int num = myVector->n;              // Number of values
    int numValid = 0;                   // Number of valid values

    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
    bool useRange = stats->options & PS_STAT_USE_RANGE;

    for (long i = 0; i < num; i++) {
        // Check if the data is with the specified range
        if (useRange && (vector[i] < stats->min))
            continue;
        if (useRange && (vector[i] > stats->max))
            continue;
        if (maskData && (maskData[i] & maskVal))
            continue;
        numValid++;
        max = PS_MAX (vector[i], max);
        min = PS_MIN (vector[i], min);
    }

    // XXX save numValid in psStats?
    if (numValid == 0) {
        stats->max = NAN;
        stats->min = NAN;
    } else {
        stats->max = max;
        stats->min = min;
        stats->results |= PS_STAT_MIN;
        stats->results |= PS_STAT_MAX;
    }
    psTrace("psLib.math", 4, "---- %s(%d) end ----\n", __func__, numValid);
    return numValid;
}

/******************************************************************************
vectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the median
and quartiles of the input vector.  Returns true on success (including if there
were no valid input vector elements). Expects F32 vector for input.
*****************************************************************************/
static bool vectorSampleMedian(const psVector* inVector,
                               const psVector* maskVector,
                               psMaskType maskVal,
                               psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);

    bool useRange = stats->options & PS_STAT_USE_RANGE;
    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; // Dereference the vector
    psF32 *input = inVector->data.F32; // Dereference the vector

    // Allocate temporary vectors for the data.
    psVector *outVector = psVectorAllocEmpty(inVector->n, PS_TYPE_F32); // Vector containing the unmasked values
    psF32 *output = outVector->data.F32; // Dereference the vector

    long count = 0;                     // Number of valid entries

    // Store all non-masked data points within the min/max range
    // into the temporary vectors.
    for (long i = 0; i < inVector->n; i++) {
        if (useRange && (input[i] < stats->min))
            continue;
        if (useRange && (input[i] > stats->max))
            continue;
        if (maskData && (maskData[i] & maskVal))
            continue;

        output[count] = input[i];
        count++;
    }
    outVector->n = count;

    if (count == 0) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n");
        psFree(outVector);
        stats->sampleUQ = NAN;
        stats->sampleLQ = NAN;
        stats->sampleMedian = NAN;
        return false;
    }

    // Sort the temporary vector.
    // XXX this is messed up: should psVectorSort free on error or not?
    if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
        psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
        stats->sampleUQ = NAN;
        stats->sampleLQ = NAN;
        stats->sampleMedian = NAN;
        return false;
    }

    // Calculate the median.  Use the average if the number of samples if even.
    if (count % 2 == 0) {
        stats->sampleMedian = 0.5 * (output[(count/2) - 1] + output[count/2]);
    } else {
        stats->sampleMedian = output[count/2];
    }

    // Calculate the quartile points exactly.
    stats->sampleUQ = output[(int)(0.75*count)];
    stats->sampleLQ = output[(int)(0.25*count)];

    stats->results |= PS_STAT_SAMPLE_MEDIAN;
    stats->results |= PS_STAT_SAMPLE_QUARTILE;

    // Free the temporary data structures.
    psFree(outVector);

    // Return "true" on success.
    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
    return true;
}

/******************************************************************************
vectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
stdev of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 
using the method below, with a single loop for various options costs only a small amount and is
much easier to debug. running 10000 tests of 1000 point vectors for the two methods gives:
 
                     single
(mask: 0, range: 0): 0.193 sec
(mask: 1, range: 0): 0.257 sec
(mask: 0, range: 0): 0.349 sec
(mask: 1, range: 0): 0.401 sec
 
*****************************************************************************/
static bool vectorSampleStdev(const psVector* myVector,
                              const psVector* errors,
                              const psVector* maskVector,
                              psMaskType maskVal,
                              psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);

    // This procedure requires the mean.  If it has not been already
    // calculated, then call vectorSampleMean()
    if (!(stats->results & PS_STAT_SAMPLE_MEAN)) {
        vectorSampleMean(myVector, errors, maskVector, maskVal, stats);
    }

    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    if (isnan(stats->sampleMean)) {
        stats->sampleStdev = NAN;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): vectorSampleMean() reported a NAN mean.\n");
        psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
        return false;
    }

    psF32 *data = myVector->data.F32;   // Dereference
    psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8;
    bool useRange = stats->options & PS_STAT_USE_RANGE;
    psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32;

    // Accumulate the sums
    double mean = stats->sampleMean;    // The mean
    double sumSquares = 0.0;            // Sum of the squares
    double sumDiffs = 0.0;              // Sum of the differences
    double errorDivisor = 0.0;          // Division for errors
    long count = 0;                     // Number of data points being used
    for (long i = 0; i < myVector->n; i++) {
        // Check if the data is with the specified range
        if (useRange && (data[i] < stats->min))
            continue;
        if (useRange && (data[i] > stats->max))
            continue;
        if (maskData && (maskData[i] & maskVal))
            continue;

        double diff = data[i] - mean;
        sumSquares += PS_SQR(diff);
        sumDiffs += diff;
        count++;
        if (errors) {
            errorDivisor += 1.0 / PS_SQR(errorsData[i]);
        }
    }

    if (count == 0) {
        stats->sampleStdev = NAN;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld).  Setting stats->sampleStdev = NAN.\n", count);
        return false;
    }
    if (count == 1) {
        stats->sampleStdev = 0.0;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld).  Setting stats->sampleStdev = 0.0.\n", count);
        return false;
    }

    if (errors) {
        stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
    } else {
        stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / (float)(count - 1));
    }
    stats->results |= PS_STAT_SAMPLE_STDEV;

    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);

    return true;
}

/******************************************************************************
vectorClippedStats(myVector, errors, maskInput, maskValInput, stats): calculates the
clipped stats (mean or stdev) of the input vector.
 
Inputs
    myVector
    maskInput
    maskValInput
    stats
Returns
    true for success; false otherwise
*****************************************************************************/
static bool vectorClippedStats(const psVector* myVector,
                               const psVector* errors,
                               psVector* maskInput,
                               psMaskType maskValInput,
                               psStats* stats
                              )
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    psTrace("psLib.math", 4, "Trace level is %d\n", psTraceGetLevel("psLib.math"));

    // Ensure that stats->clipIter is within the proper range.
    PS_ASSERT_INT_WITHIN_RANGE(stats->clipIter,
                               PS_CLIPPED_NUM_ITER_LB,
                               PS_CLIPPED_NUM_ITER_UB, -1);

    // Ensure that stats->clipSigma is within the proper range.
    PS_ASSERT_FLOAT_WITHIN_RANGE(stats->clipSigma,
                                 PS_CLIPPED_SIGMA_LB,
                                 PS_CLIPPED_SIGMA_UB, -1);

    // unless we succeed, these will have NAN values
    stats->clippedMean = NAN;
    stats->clippedStdev = NAN;
    stats->clippedNvalues = 0;

    // We copy the mask vector, to preserve the original
    psMaskType maskVal = MASK_MARK | maskValInput;
    psVector *tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);
    psVectorInit(tmpMask, 0);
    if (maskInput) {
        for (long i = 0; i < myVector->n; i++) {
            if (maskInput->data.U8[i] & maskValInput) {
                tmpMask->data.U8[i] = maskVal;
            }
        }
    }

    // 1. Compute the sample median, which we save for output
    vectorSampleMedian(myVector, tmpMask, maskVal, stats);
    if (isnan(stats->sampleMedian)) {
        psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleMedian returned NAN\n");
        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
        psFree(tmpMask);
        return false;
    }
    psTrace("psLib.math", 6, "The initial sample median is %f\n", stats->sampleMedian);

    // 2. Compute the sample standard deviation, which we save for output
    vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
    if (isnan(stats->sampleStdev)) {
        psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleStdev returned NAN\n");
        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
        psFree(tmpMask);
        return false;
    }
    psTrace("psLib.math", 6, "The initial sample stdev is %f\n", stats->sampleStdev);

    // 3. Use the sample median as the first estimator of the mean X.
    psF32 clippedMean = stats->sampleMedian;

    // 4. Use the sample stdev as the first estimator of the mean stdev.
    psF32 clippedStdev = stats->sampleStdev;

    // 5. Repeat N (stats->clipIter) times:
    long numClipped = 0;                // Number of values clipped
    bool clipped = true;                // Have we clipped anything in this iteration
    for (int iter = 0; iter < stats->clipIter && clipped; iter++) {
        clipped = false;
        psTrace("psLib.math", 6, "------------ Iteration %d ------------\n", iter);
        // a) Exclude all values x_i for which |x_i - x| > K * stdev
        if (errors) {
            // XXXX if we convert errors to variance, this should square the other terms (A*A faster than sqrt(A))
            for (long j = 0; j < myVector->n; j++) {
                if (!tmpMask->data.U8[j] &&
                        fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) {
                    tmpMask->data.U8[j] = 0xff;
                    psTrace("psLib.math", 10, "Clipped %ld: %f +/- %f\n", j,
                            myVector->data.F32[j], errors->data.F32[j]);
                    numClipped++;
                    clipped = true;
                }
            }
        } else {
            for (long j = 0; j < myVector->n; j++) {
                if (!tmpMask->data.U8[j] &&
                        fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
                    tmpMask->data.U8[j] = 0xff;
                    psTrace("psLib.math", 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]);
                    numClipped++;
                    clipped = true;
                }
            }
        }

        // b) compute new mean and stdev
        // Allocate a psStats structure for calculating the mean and stdev.
        psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
        vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp);
        vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
        psTrace("psLib.math", 6, "The new sample mean is %f\n", statsTmp->sampleMean);
        psTrace("psLib.math", 6, "The new sample stdev is %f\n", statsTmp->sampleStdev);

        // If the new mean and stdev are NAN, we must exit the loop.
        // Otherwise, use the new results and continue.
        if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) {
            // Exit loop.
            iter = stats->clipIter;
            psError(PS_ERR_UNKNOWN, true, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
            psFree(tmpMask);
            return false;
        } else {
            clippedMean = statsTmp->sampleMean;
            clippedStdev = statsTmp->sampleStdev;
        }
        psFree(statsTmp);
    }

    // Number of values used in calculation is the total number of data values, minus those we clipped
    stats->clippedNvalues = myVector->n - numClipped;

    // 7. The last calcuated value of x is the clipped mean.
    // 8. The last calcuated value of stdev is the clipped stdev.
    // we always return both stats even if only one was requested
    stats->clippedMean = clippedMean;
    stats->clippedStdev = clippedStdev;

    stats->results |= PS_STAT_CLIPPED_MEAN;
    stats->results |= PS_STAT_CLIPPED_STDEV;

    psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean);
    psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev);

    psFree(tmpMask);
    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
    return true;
}

/******************************************************************************
vectorRobustStats(myVector, maskInput, maskValInput, stats): This is the new
version of the robust stats routine.
 
XXX: MUST DO: If the errors in the input values are known, then the same
approach is used, except that the histograms become probability density
functions (PDFs). In this case, the input values are spread out, so that they
do not simply contribute a single unit to the histogram, but rather contribute
a fraction of a value, equivalent to the weight. In the interests of speed, a
boxcar PDF may be used to represent each input value (as opposed to a
Gaussian), where the boxcar width is equal to 2p2 ln 2 times the error and
each input value contributes constant area. Then the robust median and
standard deviation are estimated in the same manner as above.
 
XXX: Check for errors in psLib routines that we call.
 
XXX: Review and ensure that all memory is free'ed at premature exits.
*****************************************************************************/
#define INITIAL_NUM_BINS 1000.0
static bool vectorRobustStats(const psVector* myVector,
                              const psVector* errors,
                              psVector* maskInput,
                              psMaskType maskValInput,
                              psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    if (psTraceGetLevel("psLib.math") >= 8) {
        PS_VECTOR_PRINT_F32(myVector);
    }

    psScalar tmpScalar;                 // Temporary scalar variable, for p_psVectorBinDisect
    tmpScalar.type.type = PS_TYPE_F32;

    // we need to generate an internal mask and maskVal which ensures a bit is set
    // and tested even if there is no supplied mask (and/or the maskVal is 0)
    // XXX this would be better if we had globally defined mask values
    psMaskType maskVal = MASK_MARK | maskValInput;
    psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_MASK); // The actual mask we will use
    psVectorInit(mask, 0);
    if (maskInput) {
        for (long i = 0; i < myVector->n; i++) {
            if (maskInput->data.U8[i] & maskValInput) {
                mask->data.U8[i] = maskVal;
            }
        }
    }

    // statsMinMax is only applied to a subset of the data points
    psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max
    psHistogram *histogram = NULL;      // Histogram of the data
    psHistogram *cumulative = NULL;     // Cumulative histogram of the data
    float min = NAN, max = NAN;         // Mimimum and maximum values
    float sigma = NAN;                  // The robust standard deviation
    long totalDataPoints = 0;           // Total number of (unmasked) data points

    // Iterate to get the best bin size
    for (int iterate = 1; iterate > 0; iterate++) {
        psTrace("psLib.math", 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);

        // Get the minimum and maximum values
        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of valid values
        min = statsMinMax->min;
        max = statsMinMax->max;
        if (numValid == 0 || isnan(min) || isnan(max)) {
            // Data range calculation failed
            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__);
            psFree(mask);
            return false;
        }
        psTrace("psLib.math", 6, "Data min/max is (%.2f, %.2f)\n", min, max);

        // If all data points have the same value, then we set the appropriate members of stats and return.
        if (fabs(max - min) <= FLT_EPSILON) {
            psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min);
            stats->robustMedian = min;
            stats->robustStdev = 0.0;
            stats->robustUQ = min;
            stats->robustLQ = min;
            stats->robustN50 = numValid;
            psFree(statsMinMax);

            stats->results |= PS_STAT_ROBUST_MEDIAN;
            stats->results |= PS_STAT_ROBUST_STDEV;
            stats->results |= PS_STAT_ROBUST_QUARTILE;

            psTrace("psLib.math", 4, "---- %s(0) end  ----\n", __func__);
            psFree(mask);
            return true;
        }

        float binSize = 0.0;            // Size of bins for histogram
        if ((iterate == 1) && (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 {
            // Determine the bin size of the robust histogram, using the pre-defined number of bins
            binSize = (max - min) / INITIAL_NUM_BINS;
        }
        psTrace("psLib.math", 6, "Initial robust bin size is %.2f\n", binSize);

        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can not specify the bin
        // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize.  If
        // we get here, we know that binSize != 0.0.
        long numBins = (max - min) / binSize; // Number of bins
        psTrace("psLib.math", 6, "Numbins is %ld\n", numBins);
        psTrace("psLib.math", 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
        // Generate the histogram
        histogram = psHistogramAlloc(min, max, numBins);
        // XXXXX we need to consider this step if errors -> variance
        histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal);
        if (psTraceGetLevel("psLib.math") >= 8) {
            PS_VECTOR_PRINT_F32(histogram->bounds);
            PS_VECTOR_PRINT_F32(histogram->nums);
        }

        // ADD step 1: Convert the specific histogram to a cumulative histogram
        cumulative = psHistogramAlloc(min, max, numBins);
        cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
        for (long i = 1; i < histogram->nums->n; i++) {
            cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
        }
        if (psTraceGetLevel("psLib.math") >= 8) {
            PS_VECTOR_PRINT_F32(cumulative->bounds);
            PS_VECTOR_PRINT_F32(cumulative->nums);
        }

        // ADD step 2: Find the bin which contains the 50% data point.
        totalDataPoints = cumulative->nums->data.F32[numBins - 1];
        psTrace("psLib.math", 6, "Total data points is %ld\n", totalDataPoints);
        long binMedian;
        if (totalDataPoints/2.0 < cumulative->nums->data.F32[0]) {
            // Special case: the median is in the first bin.  Perhaps we should recode this.
            // XXX: Must try a special case where the median in in the last bin.
            binMedian = 0;
        } else {
            tmpScalar.data.F32 = totalDataPoints/2.0;
            binMedian = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
            if (binMedian < 0) {
                psError(PS_ERR_UNKNOWN, false,
                        "Failed to calculate the 50 precent data point (%ld).\n", binMedian);
                psFree(statsMinMax);
                psFree(histogram);
                psFree(cumulative);
                psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
                psFree(mask);
                return false;
            }
        }
        psTrace("psLib.math", 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian,
                cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);

        // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median.
        stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
                              binMedian, totalDataPoints/2.0);
        if (isnan(stats->robustMedian)) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to fit a quadratic and calculate the 50-percent position.\n");
            psFree(statsMinMax);
            psFree(histogram);
            psFree(cumulative);
            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
            psFree(mask);
            return false;
        }
        psTrace("psLib.math", 6, "Current robust median is %f\n", stats->robustMedian);

        // ADD step 4: Find the bins which contains the 15.8655% and 84.1345% data points.
        long binLo = 0;                 // Bin containing the 15.8655% mark
        tmpScalar.data.F32 = totalDataPoints * 0.158655f;
        if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) {
            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
            binLo = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
            if (binLo > cumulative->nums->n-1) {
                binLo = cumulative->nums->n-1;
            }
        }

        long binHi = 0;                 // Bin containing the 84.1345% mark
        tmpScalar.data.F32 = totalDataPoints * 0.841345f;
        if (tmpScalar.data.F32 > cumulative->nums->data.F32[0]) {
            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
            binHi = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
            if (binHi > cumulative->nums->n-1) {
                binHi = cumulative->nums->n-1;
            }
        }
        psTrace("psLib.math", 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n",
                binLo, binHi);
        psTrace("psLib.math", 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binLo));
        psTrace("psLib.math", 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binHi));

        if ((binLo < 0) || (binHi < 0)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
            psFree(statsMinMax);
            psFree(histogram);
            psFree(cumulative);
            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
            psFree(mask);
            return false;
        }

        // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma positions.
        psTrace("psLib.math", 6, "binLo is %ld.  Nums at that bin and the next are (%.2f, %.2f)\n",
                binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]);
        psTrace("psLib.math", 6, "binHi is %ld.  Nums at that bin and the next are (%.2f, %.2f)\n",
                binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]);

        float deltaBounds, deltaNums, prevPixels;
        float percentNums, base, binLoF32;

        // find the -1 sigma points with linear interpolation
        deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
        if (binLo == 0) {
            deltaNums = cumulative->nums->data.F32[0];
            prevPixels = 0;
        } else {
            deltaNums = cumulative->nums->data.F32[binLo] - cumulative->nums->data.F32[binLo-1];
            prevPixels = cumulative->nums->data.F32[binLo-1];
        }
        percentNums = (totalDataPoints * 0.158655f) - prevPixels;
        base = cumulative->bounds->data.F32[binLo];
        binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark
        psTrace("psLib.math", 6,
                "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                base, deltaBounds, deltaNums, prevPixels, percentNums);

        // find the +1 sigma points with linear interpolation
        deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi];
        if (binHi == 0) {
            deltaNums = cumulative->nums->data.F32[0];
            prevPixels = 0;
        } else {
            deltaNums = cumulative->nums->data.F32[binHi] - cumulative->nums->data.F32[binHi-1];
            prevPixels = cumulative->nums->data.F32[binHi-1];
        }
        percentNums = (totalDataPoints * 0.841345f) - prevPixels;
        base = cumulative->bounds->data.F32[binHi];
        float binHiF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 84.1345% mark
        psTrace("psLib.math", 6,
                "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                base, deltaBounds, deltaNums, prevPixels, percentNums);

        // report +/- 1 sigma points
        psTrace("psLib.math", 5,
                "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
                binLoF32, binHiF32);

        // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions.
        sigma = (binHiF32 - binLoF32) / 2.0;
        psTrace("psLib.math", 6, "The current sigma is %f.\n", sigma);
        stats->robustStdev = sigma;

        // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more
        // than 25 bins from the median, recalculate the bin size, and perform the algorithm again.
        if (sigma < (2 * binSize)) {
            psTrace("psLib.math", 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
            long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
            long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking
            psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index
            psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index
            psTrace("psLib.math", 6, "Masking data more than 25 bins from the median\n");
            psTrace("psLib.math", 6,
                    "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n",
                    binMedian, maskLo, maskHi);
            psTrace("psLib.math", 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
            for (long i = 0 ; i < myVector->n ; i++) {
                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
                    mask->data.U8[i] |= MASK_MARK;
                    psTrace("psLib.math", 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
                }
            }
            // Free the histograms; they will be recreated on the next iteration, with new bounds
            psFree(histogram);
            psFree(cumulative);
        } else {
            // We've got the bin size correct now
            psTrace("psLib.math", 6, "*************: No more iteration.  sigma is %f\n", sigma);
            iterate = -1;
        }
    }
    psFree(histogram);

    // ADD step 7: Find the bins which contains the 25% and 75% data points.
    long binLo25 = 0;                   // Bin containing the 25% value
    long binHi25 = 0;                   // Bin containing the 75% value

    tmpScalar.data.F32 = totalDataPoints * 0.25f;
    if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) {
        // Special case where its in first bin.  The last bin should be handled automatically.
        binLo25 = 0;
    } else {
        binLo25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
    }
    tmpScalar.data.F32 = totalDataPoints * 0.75f;
    if (tmpScalar.data.F32 <= cumulative->nums->data.F32[0]) {
        // Special case where its in first bin.  The last bin should be handled automatically.
        binHi25 = 0;
    } else {
        binHi25 = 1 + p_psVectorBinDisect(cumulative->nums, &tmpScalar);
    }
    if ((binLo25 < 0) || (binHi25 < 0)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25 and 75 percent data points\n");
        psFree(statsMinMax);
        psFree(cumulative);
        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
        psFree(mask);
        return false;
    }
    psTrace("psLib.math", 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25);

    // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
    // positions.
    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
                       binLo25, totalDataPoints * 0.25f);
    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums,
                       binHi25, totalDataPoints * 0.75f);
    psFree(cumulative);

    if (isnan(binLo25F32) || isnan(binHi25F32)) {
        psError(PS_ERR_UNKNOWN, false,
                "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
        psFree(statsMinMax);
        psTrace("psLib.math", 4, "---- %s(1) end  ----\n", __func__);
        psFree(mask);
        return false;
    }

    stats->results |= PS_STAT_ROBUST_QUARTILE;
    stats->robustLQ = binLo25F32;
    stats->robustUQ = binHi25F32;
    psTrace("psLib.math", 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n",
            binLo25F32, binHi25F32);
    long N50 = 0;
    for (long i = 0 ; i < myVector->n ; i++) {
        if (!mask->data.U8[i] &&
                (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) {
            N50++;
        }
    }
    stats->robustN50 = N50;
    psTrace("psLib.math", 6, "The robustN50 is %ld.\n", N50);

    // Clean up
    psFree(statsMinMax);
    psFree(mask);

    stats->results |= PS_STAT_ROBUST_MEDIAN;
    stats->results |= PS_STAT_ROBUST_STDEV;

    psTrace("psLib.math", 4, "---- %s(0) end  ----\n", __func__);
    return true;
}

/*
 * vectorFittedStats requires guess for fittedMean and fittedStdev
 * robustN50 should also be set
 */
static bool vectorFittedStats (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 guess mean  is %f.\n", guessMean);
    psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev);

    bool done = false;
    for (int iteration = 0; !done && (iteration < 2); iteration ++) {
        psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max

        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/10 < binsize < sigma)
            // set roughly so that the lowest bins have about 2 cnts
            // Nsmallest ~ N50 / (4*dN))
            psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
            binSize = guessStdev / dN;
        }

        // Determine the min/max of the vector (which prior outliers masked out)
        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
        float min = statsMinMax->min;
        float max = statsMinMax->max;
        if (numValid == 0 || isnan(min) || isnan(max)) {
            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);
        }

        long binMin = 0;                // Low bin to check
        long binMax = 0;                // High bin to check

        // Fit a Gaussian to the bins in the requested range about the guess mean
        // min,max overrides clipSigma
        psF32 maxFitSigma = 2.0;
        if (isfinite(stats->clipSigma)) {
            maxFitSigma = fabs(stats->clipSigma);
        }
        if (isfinite(stats->max)) {
            maxFitSigma = fabs(stats->max);
        }

        psF32 minFitSigma = 2.0;
        if (isfinite(stats->clipSigma)) {
            minFitSigma = fabs(stats->clipSigma);
        }
        if (isfinite(stats->min)) {
            minFitSigma = fabs(stats->min);
        }

        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 - 1;
        } else {
            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
        }
        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;
        }

        // Generate the variables that will be used in the Gaussian fitting
        // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
        psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates
        for (long i = binMin, j = 0; i <= binMax ; i++, j++) {
            y->data.F32[j] = histogram->nums->data.F32[i];
            psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value
            ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i);
            x->data[j] = ordinate;
        }
        if (psTraceGetLevel("psLib.math") >= 8) {
            // XXX: Print the x array somehow.
            PS_VECTOR_PRINT_F32(y);
        }
        psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n);
        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);

        // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
        // XXX EAM : why not just fit the normalization as well??
        p_psNormalizeVectorRange(y, 0.0, 1.0);

        // Fit a Gaussian to the data.
        psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
        psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
        // Initial guess for the mean (index 0) and var (index 1).
        params->data.F32[0] = guessMean;
        params->data.F32[1] = PS_SQR(guessStdev);
        if (psTraceGetLevel("psLib.math") >= 8) {
            PS_VECTOR_PRINT_F32(params);
            PS_VECTOR_PRINT_F32(y);
        }
        if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
            psFree(params);
            psFree(minimizer);
            psFree(x);
            psFree(y);
            psFree(histogram);
            psFree(statsMinMax);
            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
            return false;
        }
        if (psTraceGetLevel("psLib.math") >= 8) {
            PS_VECTOR_PRINT_F32(params);
        }

        guessMean = params->data.F32[0];
        guessStdev = sqrt(params->data.F32[1]);
        if (guessStdev > 0.75*stats->robustStdev) {
            done = true;
        }

        // Clean up after fitting
        psFree (params);
        psFree (minimizer);
        psFree (x);
        psFree (y);
        psFree (histogram);
        psFree (statsMinMax);
    }

    // The fitted mean is the Gaussian mean.
    stats->fittedMean = guessMean;
    psTrace("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;
    stats->results |= PS_STAT_FITTED_STDEV;

    return true;
}


/********************
 * vectorFittedStats_v2 requires guess for fittedMean and fittedStdev
 * robustN50 should also be set
 * gaussian fit is performed using 2D polynomial to ln(y)
 ********************/
static bool vectorFittedStats_v2 (const psVector* myVector,
                                  const psVector* errors,
                                  psVector* mask,
                                  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/10 < binsize < sigma)
            // set roughly so that the lowest bins have about 2 cnts
            // Nsmallest ~ N50 / (4*dN))
            psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8));
            binSize = guessStdev / dN;
        }

        // Determine the min/max of the vector (which prior outliers masked out)
        int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values
        float min = statsMinMax->min;
        float max = statsMinMax->max;
        if (numValid == 0 || isnan(min) || isnan(max)) {
            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);
        }

        long binMin = 0;                // Low bin to check
        long binMax = 0;                // High bin to check

        // Fit a Gaussian to the bins in the requested range about the guess mean
        // min,max overrides clipSigma
        psF32 maxFitSigma = 2.0;
        if (isfinite(stats->clipSigma)) {
            maxFitSigma = fabs(stats->clipSigma);
        }
        if (isfinite(stats->max)) {
            maxFitSigma = fabs(stats->max);
        }

        psF32 minFitSigma = 2.0;
        if (isfinite(stats->clipSigma)) {
            minFitSigma = fabs(stats->clipSigma);
        }
        if (isfinite(stats->min)) {
            minFitSigma = fabs(stats->min);
        }

        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 - 1;
        } else {
            binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
        }
        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;
        }

        // Generate the variables that will be used in the Gaussian fitting
        // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins
        psVector *y = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates
        psVector *x = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of ordinates
        long j = 0;
        for (long i = binMin; i <= binMax ; i++) {
            if (histogram->nums->data.F32[i] <= 0.0)
                continue;
            x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i);
            y->data.F32[j] = log(histogram->nums->data.F32[i]);
            // XXX note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2)
            j++;
        }
        y->n = x->n = j;
        if (psTraceGetLevel("psLib.math") >= 8) {
            // XXX: Print the x array somehow.
            PS_VECTOR_PRINT_F32(y);
        }
        psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n);
        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);

        // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
        psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);

        // XXX how can we protect against some extreme outliers?  the robust histogram
        // should have already dealt with those, but we are again sensitive to them...
        // psStats *fitStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
        // fitStats->clipIter = 3.0;
        // fitStats->clipSigma = 3.0;
        // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_U8);
        // psVectorInit (fitMask, 0);

        if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
            psFree(x);
            psFree(y);
            psFree(poly);
            psFree(histogram);
            psFree(statsMinMax);
            psTrace("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);
            psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
            return false;
        }

        guessStdev = sqrt(-0.5/poly->coeff[2]);
        guessMean = poly->coeff[1]*PS_SQR(guessStdev);
        if (guessStdev > 0.75*stats->robustStdev) {
            done = true;
        } else {
            psTrace("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);
        }

        // Clean up after fitting
        psFree (x);
        psFree (y);
        psFree (poly);
        // psFree (fitStats);
        // psFree (fitMask);
        psFree (histogram);
        psFree (statsMinMax);
    }

    // The fitted mean is the Gaussian mean.
    stats->fittedMean = guessMean;
    psTrace("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_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;
}


/******************************************************************************
vectorSmoothHistGaussian(): This routine smoothes the data in the input
robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
smoothed data.
 
XXX this function is unused
*****************************************************************************/
psVector *vectorSmoothHistGaussian(psHistogram *histogram,
                                   psF32 sigma)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    psTrace("psLib.math", 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma);
    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
    if (psTraceGetLevel("psLib.math") >= 8) {
        PS_VECTOR_PRINT_F32(histogram->nums);
    }

    long numBins = histogram->nums->n;  // Number of histogram bins
    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins
    const psVector *bounds = histogram->bounds; // The bounds for the histogram bins

    if (!histogram->uniform) {
        //
        // We get here if the histogram is non-uniform.  This code is not tested.
        // However, it is also not used anywhere, yet.
        //
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");

        for (long i = 0; i < numBins; i++) {
            // Determine the midpoint of bin i.
            float iMid = PS_BIN_MIDPOINT(histogram, i);

            //
            // We determine the bin numbers (jMin:jMax) corresponding to a
            // range of data values surrounding iMid.  The range is of size:
            // 2*PS_GAUSS_WIDTH*sigma
            long jMin, jMax;
            psF32 x = iMid - (PS_GAUSS_WIDTH * sigma);
            for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--)
                ;
            x = iMid + (PS_GAUSS_WIDTH * sigma);
            for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++)
                ;

            //
            // Loop from jMin to jMax, computing the gaussian of data i.
            //
            smooth->data.F32[i] = 0.0;
            for (long j = jMin ; j <= jMax ; j++) {
                float jMid = PS_BIN_MIDPOINT(histogram, j);
                smooth->data.F32[i] +=
                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
            }
        }
    } else {
        //
        // We get here if the histogram is uniform.
        //
        for (long i = 0; i < numBins; i++) {
            psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0];
            psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize);

            //
            // XXX: The following is wrong.  However, in practice, the sigma was too
            // large compared to the binSize.  This meant that the smoothing was done
            // over 500 bins in the robust stats algorithm.  This mean that the smoothing
            // took way too long.
            //
            #if 0

            gaussWidth = 10;
            #endif
            //
            // We determine the bin numbers (jMin:jMax) corresponding to a
            // range of data values surrounding iMid.  The range is of size:
            // 2*PS_GAUSS_WIDTH*sigma
            //
            psS32 jMin = PS_MAX(i - gaussWidth, 0);
            psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1);

            //
            // Loop from jMin to jMax, computing the gaussian of data i.
            //
            smooth->data.F32[i] = 0.0;
            float iMid = PS_BIN_MIDPOINT(histogram, i);
            for (long j = jMin ; j <= jMax ; j++) {
                float jMid = PS_BIN_MIDPOINT(histogram, j);
                smooth->data.F32[i] +=
                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
            }
        }
    }

    if (psTraceGetLevel("psLib.math") >= 8) {
        PS_VECTOR_PRINT_F32(smooth);
    }
    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
    return(smooth);
}

/*****************************************************************************/

/* FUNCTION IMPLEMENTATION - PUBLIC                                          */

/*****************************************************************************/

// We keep statsFree so that we can identify statistics pointers from the memblock
static void statsFree(psStats *newStruct)
{
    // There are non dynamic allocated items
}

/******************************************************************************
    psStatsAlloc(): This routine must create a new psStats data structure.
*****************************************************************************/
psStats* psStatsAlloc(psStatsOptions options)
{
    psTrace("psLib.math", 3,"---- %s() begin  ----\n", __func__);
    psStats* newStruct = NULL;

    newStruct = (psStats* ) psAlloc(sizeof(psStats));
    psMemSetDeallocator(newStruct, (psFreeFunc)statsFree);
    newStruct->sampleMean = NAN;
    newStruct->sampleMedian = NAN;
    newStruct->sampleStdev = NAN;
    newStruct->sampleUQ = NAN;
    newStruct->sampleLQ = NAN;
    newStruct->robustMedian = NAN;
    newStruct->robustStdev = NAN;
    newStruct->robustUQ = NAN;
    newStruct->robustLQ = NAN;
    newStruct->robustN50 = -1;
    newStruct->fittedMean = NAN;
    newStruct->fittedStdev = NAN;
    newStruct->fittedNfit = -1;
    newStruct->clippedMean = NAN;
    newStruct->clippedStdev = NAN;
    newStruct->clippedNvalues = -1;     // XXX: This is never used
    newStruct->clipSigma = 3.0;
    newStruct->clipIter = 3;
    newStruct->min = NAN;
    newStruct->max = NAN;
    newStruct->binsize = NAN;
    newStruct->nSubsample = 100000;
    newStruct->options = options;

    psTrace("psLib.math", 3, "---- %s() end  ----\n", __func__);
    return (newStruct);
}


bool psMemCheckStats(psPtr ptr)
{
    PS_ASSERT_PTR(ptr, false);
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)statsFree );
}

/******************************************************************************
psVectorStats(in, mask, maskVal, stats): this is the public API
function which calls the above private stats functions based on what bits
were set in stats->options.
 
Inputs
    in
    mask
    maskVal
    stats
Returns
    The stats structure.
 
XXX: Should we free stats if the asserts fail? NO; we don't own it (RHL).
*****************************************************************************/
bool psVectorStats(psStats* stats,
                   const psVector* in,
                   const psVector* errors,
                   const psVector* mask,
                   psMaskType maskVal)
{
    psTrace("psLib.math", 3,"---- %s() begin  ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(stats, NULL);
    PS_ASSERT_VECTOR_NON_NULL(in, NULL);
    if (mask) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, stats);
        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, stats);
    }
    if (errors) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(errors, in, stats);
        PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats);
    }

    // Convert types, as necessary
    psVector *inF32 = NULL;             // Input vector of values, F32 version
    if (in->type.type == PS_TYPE_F32) {
        inF32 = psMemIncrRefCounter((psPtr)in);
    } else {
        inF32 = psVectorCopy(NULL, in, PS_TYPE_F32);
    }
    psVector *errorsF32 = NULL;         // Input vector of errors, F32 version
    if (errors) {
        if (errors->type.type == PS_TYPE_F32) {
            errorsF32 = psMemIncrRefCounter((psPtr)errors);
        } else {
            errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32);
        }
    }
    psVector *maskU8 = NULL;            // Input mask vector, U8 version
    if (mask) {
        if (mask->type.type == PS_TYPE_MASK) {
            maskU8 = psMemIncrRefCounter((psPtr)mask);
        } else {
            maskU8 = psVectorCopy(NULL, mask, PS_TYPE_MASK);
        }
    }

    if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->max, stats->min, stats);
    }

    if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) {
        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, stats);
    }

    // init the value of stats->results: this is used internally to check if
    // prior functions have been called
    stats->results = PS_STAT_NONE;
    bool status = true;

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEAN) {
        if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean");
            status &= false;
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) {
        if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median");
            status &= false;
        }
    }

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_STDEV) {
        if (!vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev");
            status &= false;
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) {
        if (vectorMinMax(inF32, maskU8, maskVal, stats) == 0) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum");
            status &= false;
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) {
        if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics"));
            status &= false;
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
        if (!vectorFittedStats(inF32, errorsF32, maskU8, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
            status &= false;
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2)) {
        if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
            psAbort ("stats", "you may not specify both FITTED_MEAN and FITTED_MEAN_V2");
        }
        if (!vectorFittedStats_v2(inF32, errorsF32, maskU8, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics"));
            status &= false;
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_STDEV_V3)) {
        if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
            psAbort ("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)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n");
            status &= false;
        }
    }

    psFree(inF32);
    psFree(errorsF32);
    psFree(maskU8);
    psTrace("psLib.math", 3,"---- %s() end  ----\n", __func__);
    return (status);
}

psStatsOptions psStatsOptionFromString(const char *string)
{
    #define READ_STAT(NAME, SYMBOL) \
    if (strcasecmp(string, NAME) == 0) { \
        return SYMBOL; \
    }

    READ_STAT("MEAN",       PS_STAT_SAMPLE_MEAN);
    READ_STAT("STDEV",      PS_STAT_SAMPLE_STDEV);
    READ_STAT("MEDIAN",     PS_STAT_SAMPLE_MEDIAN);
    READ_STAT("QUARTILE",   PS_STAT_SAMPLE_QUARTILE);
    READ_STAT("SAMPLE_MEAN",     PS_STAT_SAMPLE_MEAN);
    READ_STAT("SAMPLE_STDEV",    PS_STAT_SAMPLE_STDEV);
    READ_STAT("SAMPLE_MEDIAN",   PS_STAT_SAMPLE_MEDIAN);
    READ_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE);
    READ_STAT("ROBUST",          PS_STAT_ROBUST_MEDIAN);
    READ_STAT("ROBUST_MEDIAN",   PS_STAT_ROBUST_MEDIAN);
    READ_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
    READ_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE);
    READ_STAT("FITTED",         PS_STAT_FITTED_MEAN);
    READ_STAT("FITTED_MEAN",    PS_STAT_FITTED_MEAN);
    READ_STAT("FITTED_STDEV",   PS_STAT_FITTED_STDEV);
    READ_STAT("FITTED_V2",       PS_STAT_FITTED_MEAN_V2);
    READ_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
    READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
    READ_STAT("FITTED_V3",       PS_STAT_FITTED_MEAN_V3);
    READ_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN_V3);
    READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3);
    READ_STAT("CLIPPED",         PS_STAT_CLIPPED_MEAN);
    READ_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
    READ_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);

    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to parse statistic: %s\n", string);
    return 0;
}

psString psStatsOptionToString(psStatsOptions option)
{
    psString string = NULL;             // String to return

    #define WRITE_STAT(NAME, SYMBOL) \
    if (option & SYMBOL) { \
        psStringAppend(&string, "%s ", NAME); \
    }

    // Same list as above (for psStatsFromString), but with repeat symbols removed
    WRITE_STAT("SAMPLE_MEAN",     PS_STAT_SAMPLE_MEAN);
    WRITE_STAT("SAMPLE_STDEV",    PS_STAT_SAMPLE_STDEV);
    WRITE_STAT("SAMPLE_MEDIAN",   PS_STAT_SAMPLE_MEDIAN);
    WRITE_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE);
    WRITE_STAT("ROBUST_MEDIAN",   PS_STAT_ROBUST_MEDIAN);
    WRITE_STAT("ROBUST_STDEV",    PS_STAT_ROBUST_STDEV);
    WRITE_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE);
    WRITE_STAT("FITTED_MEAN",     PS_STAT_FITTED_MEAN);
    WRITE_STAT("FITTED_STDEV",    PS_STAT_FITTED_STDEV);
    WRITE_STAT("FITTED_MEAN_V2",  PS_STAT_FITTED_MEAN_V2);
    WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2);
    WRITE_STAT("FITTED_MEAN_V3",  PS_STAT_FITTED_MEAN_V3);
    WRITE_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3);
    WRITE_STAT("CLIPPED_MEAN",    PS_STAT_CLIPPED_MEAN);
    WRITE_STAT("CLIPPED_STDEV",   PS_STAT_CLIPPED_STDEV);

    return string;
}

psStats *psStatsFromString(const char *string)
{
    psList *subStrings = psStringSplit(string, " ,;", false); // List of sub-strings
    if (!subStrings || psListLength(subStrings) == 0) {
        // Nothing here
        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "No string to parse for statistics: %s\n", string);
        psFree(subStrings);
        return NULL;
    }
    psStats *stats = psStatsAlloc(0);   // Generate empty stats structure
    psListIterator *iterator = psListIteratorAlloc(subStrings, PS_LIST_HEAD, false); // Iterator
    psString statString;                // Statistic string, from iteration
    while ((statString = psListGetAndIncrement(iterator))) {
        psStatsOptions option = psStatsOptionFromString(statString);
        if (option == 0) {
            psLogMsg(__func__, PS_LOG_WARN, "Unable to interpret statistic option: %s --- ignored.\n",
                     statString);
            continue;
        }
        stats->options |= option;
    }
    psFree(iterator);
    psFree(subStrings);
    return stats;
}

psString psStatsToString(const psStats *stats)
{
    return psStatsOptionToString(stats->options);
}

inline psStatsOptions psStatsSingleOption(psStatsOptions option)
{
    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) {
    case PS_STAT_SAMPLE_MEAN:
    case PS_STAT_SAMPLE_MEDIAN:
    case PS_STAT_SAMPLE_STDEV:
    case PS_STAT_SAMPLE_QUARTILE:
    case PS_STAT_ROBUST_MEDIAN:
    case PS_STAT_ROBUST_STDEV:
    case PS_STAT_ROBUST_QUARTILE:
    case PS_STAT_FITTED_MEAN:
    case PS_STAT_FITTED_STDEV:
    case PS_STAT_FITTED_MEAN_V2:
    case PS_STAT_FITTED_STDEV_V2:
    case PS_STAT_FITTED_MEAN_V3:
    case PS_STAT_FITTED_STDEV_V3:
    case PS_STAT_CLIPPED_MEAN:
    case PS_STAT_CLIPPED_STDEV:
    case PS_STAT_MAX:
    case PS_STAT_MIN:
        return option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE);
    default:
        return 0;
    }
    psAbort(__func__, "Should never get here.\n");
    return 0;
}

inline double psStatsGetValue(const psStats *stats, psStatsOptions option)
{
    // We could call psStatsSingle to check, but it would be a waste since we effectively do it anyway
    switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) {
    case PS_STAT_SAMPLE_MEAN:
        return stats->sampleMean;
    case PS_STAT_SAMPLE_MEDIAN:
        return stats->sampleMedian;
    case PS_STAT_SAMPLE_STDEV:
        return stats->sampleStdev;
    case PS_STAT_ROBUST_MEDIAN:
        return stats->robustMedian;
    case PS_STAT_ROBUST_STDEV:
        return stats->robustStdev;
    case PS_STAT_FITTED_MEAN:
        return stats->fittedMean;
    case PS_STAT_FITTED_STDEV:
        return stats->fittedStdev;
    case PS_STAT_FITTED_MEAN_V2:
        return stats->fittedMean;
    case PS_STAT_FITTED_STDEV_V2:
        return stats->fittedStdev;
    case PS_STAT_FITTED_MEAN_V3:
        return stats->fittedMean;
    case PS_STAT_FITTED_STDEV_V3:
        return stats->fittedStdev;
    case PS_STAT_CLIPPED_MEAN:
        return stats->clippedMean;
    case PS_STAT_CLIPPED_STDEV:
        return stats->clippedStdev;
    case PS_STAT_MAX:
        return stats->max;
    case PS_STAT_MIN:
        return stats->min;
    case PS_STAT_SAMPLE_QUARTILE:
    case PS_STAT_ROBUST_QUARTILE:
        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot return a single quartile value; "
                "get them yourself.\n");
        return NAN;
    default:
        return NAN;
    }
    psAbort(__func__, "Should never get here.\n");
    return NAN;
}

// other private functions used above

static psF32 QuadraticInverse(psF32 a,
                              psF32 b,
                              psF32 c,
                              psF32 y,
                              psF32 xLo,
                              psF32 xHi
                             )
{
    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));

    psF64 x1 = -b/(2.0*a) + tmp;
    psF64 x2 = -b/(2.0*a) - tmp;

    if (xLo <= x1 && x1 <= xHi) {
        return x1;
    }
    if (xLo <= x2 && x2 <= xHi) {
        return x2;
    }
    return 0.5 * (xLo + xHi);
}


/******************************************************************************
fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
routine which fits a quadratic to three points and returns the x-value
corresponding to the input y-value.  This routine takes psVectors of x/y pairs
as input, and fits a quadratic to the 3 points surrounding element binNum in
the vectors (the midpoint between element i and i+1 is used for x[i]).  It
then determines for what value x does that quadratic f(x) = yVal (the input
parameter).
 
*****************************************************************************/
static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec,
        psVector *yVec,
        psS32 binNum,
        psF32 yVal
                                              )
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    psTrace("psLib.math", 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
    if (psTraceGetLevel("psLib.math") >= 8) {
        PS_VECTOR_PRINT_F32(xVec);
        PS_VECTOR_PRINT_F32(yVec);
    }

    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
    //    PS_ASSERT_VECTORS_SIZE_EQUAL(xVec, yVec, NAN);
    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);

    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    psF32 tmpFloat = 0.0f;

    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
        // The general case.  We have all three points.
        x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
        x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
        y->data.F64[0] = yVec->data.F32[binNum - 1];
        y->data.F64[1] = yVec->data.F32[binNum];
        y->data.F64[2] = yVec->data.F32[binNum + 1];
        psTrace("psLib.math", 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
        psTrace("psLib.math", 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
        psTrace("psLib.math", 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);

        //
        // Ensure that the y value lies within range of the y values.
        //
        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
                ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    _("Specified yVal, %g, is not within y-range, %g to %g."),
                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
        }

        //
        // Ensure that the y values are monotonic.
        //
        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
                ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
            psError(PS_ERR_UNKNOWN, true,
                    "This routine must be called with monotonically increasing or decreasing data points.\n");
            psFree(x);
            psFree(y);
            psTrace("psLib.math", 5, "---- %s() end ----\n", __func__);
            return NAN;
        }

        // Determine the coefficients of the polynomial.
        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
            psError(PS_ERR_UNEXPECTED_NULL, false,
                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  Returning NAN."));
            psFree(x);
            psFree(y);
            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
            return NAN;
        }
        psTrace("psLib.math", 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
        psTrace("psLib.math", 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
        psTrace("psLib.math", 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
        psTrace("psLib.math", 6, "Fitted y vec is (%f %f %f)\n",
                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));

        psTrace("psLib.math", 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal,
                                    x->data.F64[0], x->data.F64[2]);
        psFree(myPoly);

        if (isnan(tmpFloat)) {
            psError(PS_ERR_UNEXPECTED_NULL,
                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
            psFree(x);
            psFree(y);
            psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__);
            return(NAN);
        }
    } else {
        // These are special cases where the bin is at the beginning or end of the vector.
        if (binNum == 0) {
            // We have two points only at the beginning of the vectors x and y.
            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                              xVec->data.F32[binNum + 1]);
        } else if (binNum == (xVec->n - 1)) {
            // The special case where we have two points only at the end of
            // the vectors x and y.
            // XXX: Is this right?
            tmpFloat = xVec->data.F32[binNum];
        } else if (binNum == (xVec->n - 2)) {
            // XXX: Is this right?
            tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
        }
    }

    psTrace("psLib.math", 6, "FIT: return %f\n", tmpFloat);
    psFree(x);
    psFree(y);

    psTrace("psLib.math", 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
    return tmpFloat;
}

/******************************************************************************
NOTE: We assume unnormalized gaussians.
*****************************************************************************/
static psF32 minimizeLMChi2Gauss1D(psVector *deriv,
                                   const psVector *params,
                                   const psVector *coords
                                  )
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN);
    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN);
    PS_ASSERT_VECTOR_NON_NULL(coords, NAN);
    PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN);
    PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN);

    psF32 x = coords->data.F32[0];
    psF32 mean = params->data.F32[0];
    psF32 var = params->data.F32[1];
    psF32 dx = (x - mean);

    psF32 gauss = exp (-0.5*PS_SQR(dx)/var);
    if (deriv) {
        PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN);
        PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN);
        psF32 tmp = dx * gauss;
        deriv->data.F32[0] = tmp / var;
        deriv->data.F32[1] = tmp * dx / (var * var);
    }


    psTrace("psLib.math", 4, "---- %s() end ----\n", __func__);
    return gauss;
}
