/** @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
 *
 *  XXX: The following stats members are never used, or set in this code.
 *      stats->clippedNvalues
 *
 * XXX: Must do
 * nSubsample points
 * use ->min and ->max (PS_STAT_USE_RANGE)
 *
 *  @version $Revision: 1.192 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-12-01 22:41:33 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

#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 "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().  Should we perform the sanity
checks in each routine anyway?
 
XXX: For many of these private stats routines, what should be done if 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)?
Currently we set the value to NAN.
 
XXX: Optimization: many routines have an "empty" boolean variable which keeps
track of whether or not the vector has any valid elements.  This code can
possibly be optimized away.
 *****************************************************************************/
/******************************************************************************
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.
 *****************************************************************************/
# if (0)
    static bool vectorSampleMean(const psVector* myVector,
                                 const psVector* errors,
                                 const psVector* maskVector,
                                 psMaskType maskVal,
                                 psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);

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

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

    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different loop.
    if (!errors) {
        if (stats->options & PS_STAT_USE_RANGE) {
            if (maskVector) {
                // No errors, check range, check mask
                psU8 *maskData = maskVector->data.U8;
                for (long i = 0; i < numData; i++) {
                    // Check if the data is with the specified range
                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
                        mean += data[i];
                        count++;
                    }
                }
            } else {
                // No errors, check range, no mask
                for (long i = 0; i < numData; i++) {
                    if (stats->min <= data[i] && data[i] <= stats->max) {
                        mean += data[i];
                        count++;
                    }
                }
            }
        } else {
            if (maskVector) {
                // No errors, no range check, check mask
                psU8 *maskData = maskVector->data.U8;
                for (long i = 0; i < numData; i++) {
                    if (!(maskVal & maskData[i])) {
                        mean += data[i];
                        count++;
                    }
                }

            } else {
                // No errors, no range check, no mask
                for (long i = 0; i < numData; i++) {
                    mean += data[i];
                }
                count = numData;
            }
        }

        if (count > 0) {
            mean /= (psF32)count;
        } else {
            mean = NAN;
        }

    } else {
        psF32 sumWeights = 0.0;         // The sum of the weights
        psF32 *errorsData = errors->data.F32;
        if (stats->options & PS_STAT_USE_RANGE) {
            if (maskVector) {
                psU8 *maskData = maskVector->data.U8;
                for (long i = 0; i < numData; i++) {
                    // Check if the data is with the specified range
                    if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) {
                        float weight = 1.0 / PS_SQR(errorsData[i]);
                        mean += data[i] * weight;
                        sumWeights += weight;
                        count++;
                    }
                }
            } else {
                for (long i = 0; i < myVector->n; i++) {
                    if (stats->min <= data[i] && data[i] <= stats->max) {
                        float weight = 1.0 / PS_SQR(errorsData[i]);
                        mean += data[i] * weight;
                        sumWeights += weight;
                        count++;
                    }
                }
            }
        } else {
            if (maskVector) {
                psU8 *maskData = maskVector->data.U8;
                for (long i = 0; i < myVector->n; i++) {
                    if (!(maskVal & maskData[i])) {
                        float weight = 1.0 / PS_SQR(errorsData[i]);
                        mean += data[i] * weight;
                        sumWeights += weight;
                        count++;
                    }
                }
            } else {
                for (long i = 0; i < myVector->n; i++) {
                    float weight = 1.0 / PS_SQR(errorsData[i]);
                    mean += data[i] * weight;
                    sumWeights += weight;
                }
                count = myVector->n;
            }
        }

        if (count > 0) {
            mean /= sumWeights;
        } else {
            mean = NAN;
        }
    }

    stats->sampleMean = mean;
    if (isnan(mean)) {
        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
        return false;
    }

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

/******************************************************************************
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.
 
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.067 sec   0.073 sec
(mask: 1, range: 0): 0.098 sec  0.102 sec
(mask: 0, range: 0): 0.136 sec  0.162 sec
(mask: 1, range: 0): 0.177 sec  0.181 sec
 *****************************************************************************/
static bool vectorSampleMean(const psVector* myVector,
                             const psVector* errors,
                             const psVector* maskVector,
                             psMaskType maskVal,
                             psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);

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

    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;

    // If errors exist, we enter a slightly different loop.
    // XXX collapse to a single loop?
    if (!errors) {
        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;
            mean += data[i];
            count++;
        }
        mean = (count > 0) ? mean / (psF32) count : NAN;
    } else {
        psF32 sumWeights = 0.0;         // The sum of the weights
        psF32 *errorsData = 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;

            float weight = 1.0 / PS_SQR(errorsData[i]);
            mean += data[i] * weight;
            sumWeights += weight;
            count++;
        }
        mean = (count > 0) ? mean / sumWeights : NAN;
    }

    stats->sampleMean = mean;
    if (isnan(mean)) {
        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
        return false;
    }

    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);
    }

    if (numValid == 0) {
        stats->max = NAN;
        stats->min = NAN;
    } else {
        stats->max = max;
        stats->min = min;
    }
    psTrace("psLib.math", 4, "---- %s(%d) end ----\n", __func__, numValid);
    return numValid;
}

#if 0
/******************************************************************************
vectorCheckNonEmpty(myVector, maskVector, maskVal, stats): This routine returns
"true" if the inputPsVector has 1 or more valid elements (a valid element is an
unmasked element within the specified min/max range).  Otherwise, return
"false".
 *****************************************************************************/
static bool vectorCheckNonEmpty(const psVector* myVector,
                                const psVector* maskVector,
                                psMaskType maskVal,
                                psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(myVector, false);
    PS_ASSERT_PTR_NON_NULL(stats, false);

    if (stats->options & PS_STAT_USE_RANGE) {
        if (maskVector) {
            for (long i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
                    return true;
                }
            }
        } else {
            for (long i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
                    return true;
                }
            }
        }
    } else {
        if (maskVector != NULL) {
            for (long i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
                    return true;
                }
            }
        } else {
            if (myVector->n > 0) {
                psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
                return true;
            }
        }
    }
    psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
    return(false);
}
#endif

/******************************************************************************
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).
 *****************************************************************************/
static bool vectorSampleMedian(const psVector* myVector,
                               const psVector* maskVector,
                               psMaskType maskVal,
                               psStats* stats)
{
    psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__);

    // Allocate temporary vectors for the data.
    psVector* vector = psVectorAllocEmpty(myVector->n, PS_TYPE_F32); // Vector containing the unmasked values
    long count = 0;                     // Number of valid entries

    // Determine if we must only use data points within a min/max range.
    if (stats->options & PS_STAT_USE_RANGE) {
        // Store all non-masked data points within the min/max range
        // into the temporary vectors.
        if (maskVector != NULL) {
            for (long i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    vector->data.F32[count] = myVector->data.F32[i];
                    count++;
                }
            }
        } else {
            for (long i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    vector->data.F32[count] = myVector->data.F32[i];
                    count++;
                }
            }
        }
    } else {
        // Store all non-masked data points into the temporary vectors.
        if (maskVector) {
            for (long i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    vector->data.F32[count] = myVector->data.F32[i];
                    count++;
                }
            }
        } else {
            vector = psVectorCopy(vector, myVector, PS_TYPE_F32);
            count = myVector->n;
        }
    }
    vector->n = count;
    if (count == 0) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n");
        psFree(vector);
        stats->sampleMedian = NAN;
        return false;
    }

    // Sort the temporary vector.
    vector = psVectorSort(vector, vector); // Sort in-place (since it's a copy, it's OK)
    if (!vector) {
        psError(PS_ERR_UNEXPECTED_NULL,
                false,
                _("Failed to sort input data."));
        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
        psFree(vector);
        return false;
    }

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

    // Calculate the quartile points exactly.
    stats->sampleUQ = vector->data.F32[3 * (count / 4)];
    stats->sampleLQ = vector->data.F32[count / 4];

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

    // Return "true" on success.
    psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__);
    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.
 *****************************************************************************/
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);
}

/******************************************************************************
vectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
stdev of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 
 *****************************************************************************/
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 (isnan(stats->sampleMean)) {
        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;
    }

    // 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
    if (stats->options & PS_STAT_USE_RANGE) {
        if (maskVector) {
            for (long i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    double diff = myVector->data.F32[i] - mean;
                    sumSquares += PS_SQR(diff);
                    sumDiffs += diff;
                    count++;
                    if (errors) {
                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
                    }
                }
            }
        } else {
            for (long i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
                    double diff = myVector->data.F32[i] - mean;
                    sumSquares += PS_SQR(diff);
                    sumDiffs += diff;
                    count++;
                    if (errors) {
                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
                    }
                }
            }
            count = myVector->n;
        }
    } else {
        if (maskVector) {
            for (long i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    double diff = myVector->data.F32[i] - mean;
                    sumSquares += PS_SQR(diff);
                    sumDiffs += diff;
                    count++;
                    if (errors) {
                        errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
                    }
                }
            }
        } else {
            for (long i = 0; i < myVector->n; i++) {
                double diff = myVector->data.F32[i] - mean;
                sumSquares += PS_SQR(diff);
                sumDiffs += diff;
                count++;
                if (errors) {
                    errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]);
                }
            }
            count = myVector->n;
        }
    }

    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));
    }
    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);

    // Allocate a psStats structure for calculating the mean, median, and
    // stdev.
    psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);

    // 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.
    vectorSampleMedian(myVector, tmpMask, maskVal, statsTmp);
    if (isnan(statsTmp->sampleMedian)) {
        psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleMedian returned NAN\n");
        stats->clippedMean = NAN;
        stats->clippedStdev = NAN;
        psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__);
        psFree(tmpMask);
        psFree(statsTmp);
        return false;
    }
    psTrace("psLib.math", 6, "The initial sample median is %f\n", statsTmp->sampleMedian);

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

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

    // 4. Use the sample stdev as the first estimator of the mean stdev.
    psF32 clippedStdev = statsTmp->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) {
            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
        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;
        }
    }

    // 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.
    if (stats->options & PS_STAT_CLIPPED_MEAN) {
        stats->clippedMean = clippedMean;
        psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean);
    }
    // 8. The last calcuated value of stdev is the clipped stdev.
    if (stats->options & PS_STAT_CLIPPED_STDEV) {
        stats->clippedStdev = clippedStdev;
        psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev);
    }

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

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 0

    psF64 y1 = (a * x1 * x1) + (b * x1) + c;
    psF64 y2 = (a * x2 * x2) + (b * x2) + c;
    printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c);
    printf("QuadraticInverse: y is %f\n", y);
    printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2);
    printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2);
    #endif

    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);
        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);
        if (myPoly == NULL) {
            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;
}

/*
 * 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)
{

    psTrace("psLib.math", 6, "The guess mean  is %f.\n", stats->fittedMean);
    psTrace("psLib.math", 6, "The guess stdev is %f.\n", stats->fittedStdev);

    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 = stats->fittedStdev / 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.
    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 = stats->fittedMean - minFitSigma*stats->fittedStdev;
    if (tmpScalar.data.F32 < min) {
        binMin = 0;
    } else {
        binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar);
    }
    tmpScalar.data.F32 = stats->fittedMean + maxFitSigma*stats->fittedStdev;
    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);

    // 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] = stats->fittedMean;
    params->data.F32[1] = PS_SQR(stats->fittedStdev);
    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(statsMinMax);
        psFree(histogram);
        psFree(x);
        psFree(y);
        psFree(minimizer);
        psFree(params);
        psTrace("psLib.math", 4, "---- %s(false) end  ----\n", __func__);
        return false;
    }
    if (psTraceGetLevel("psLib.math") >= 8) {
        PS_VECTOR_PRINT_F32(params);
    }

    // The fitted mean is the Gaussian mean.
    stats->fittedMean = params->data.F32[0];
    psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean);

    // The fitted standard deviation
    stats->fittedStdev = sqrt(params->data.F32[1]);
    psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev);

    // Clean up after fitting
    psFree (histogram);
    psFree (x);
    psFree (y);
    psFree (minimizer);
    psFree (params);
    psFree (statsMinMax);
    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;
            }
        }
    }

    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__);
            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);
            if (stats->options & PS_STAT_ROBUST_MEDIAN) {
                stats->robustMedian = min;
            }
            if (stats->options & PS_STAT_ROBUST_QUARTILE) {
                stats->robustUQ = min;
                stats->robustLQ = min;
            }
            stats->robustN50 = numValid;
            psFree(statsMinMax);

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

        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);
        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__);
                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__);
            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__);
            return false;
        }
        // XXXX we can probably just use the bin values to get sigma without interpolating...
        // ADD step 4b: Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
        #if 0
        // XXX: I am overriding the ADD for now.  The following code implements the ADD exactly and is having
        // problems fitting a 2nd-order polynomial to data y-values suchs as (1, 1, 100).  Therefore, I
        // commented out the code and am doing simply linear interpolation.

        // Value for the 15.8655% mark
        float binLoF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binLo,
                         totalDataPoints * 0.158655f);
        // Value for the 84.1345% mark
        float binHiF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binHi,
                         totalDataPoints * 0.841345f);
        psTrace("psLib.math", 6, "The exact 15.8655%% and 84.1345%% data point positions are: (%f, %f)\n",
                binLoF32, binHiF32);
        #else
        // This code basically interpolates to find the positions exactly.
        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 = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo];
        float deltaNums;
        float prevPixels;
        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];
        }
        float percentNums = (totalDataPoints * 0.158655f) - prevPixels;
        float base = cumulative->bounds->data.F32[binLo];
        float 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);

        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);
        psTrace("psLib.math", 6,
                "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
                binLoF32, binHiF32);
        #endif

        // 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__);
        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__);
        return false;
    }

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

    // Fitted statistics
    if (stats->options & PS_STAT_FITTED_MEAN || stats->options & PS_STAT_FITTED_STDEV) {
        stats->fittedStdev = stats->robustStdev;  // pass the guess sigma
        stats->fittedMean = stats->robustMedian;  // pass the guess mean
        if (!vectorFittedStats(myVector, errors, mask, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 1)");
            psFree(statsMinMax);
            psFree(mask);

            return false;
        }
        // if there is a large swing in sigma, try a second time
        if ((stats->fittedStdev/stats->robustStdev) < 0.75) {
            if (!vectorFittedStats (myVector, errors, mask, maskVal, stats)) {
                psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 2)");
                psFree(statsMinMax);
                psFree(mask);

                return false;
            }
        }
    }

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

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

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

/* FUNCTION IMPLEMENTATION - PUBLIC                                          */

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

static void histogramFree(psHistogram* myHist);

// 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 );
}



/******************************************************************************
psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure
with the specifed upper and lower limits, and the specifed number of bins.
This routine will also set the bounds for each of the bins.
 
Input:
    lower
    upper
    n
Returns:
    The histogram structure
 *****************************************************************************/
psHistogram* psHistogramAlloc(float lower, float upper, int n)
{
    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
    psTrace("psLib.math", 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n);
    PS_ASSERT_INT_POSITIVE(n, NULL);
    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL);

    // Allocate memory for the new histogram structure.  If there are N bins, then there are N+1 bounds to
    // those bins.
    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
    psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
    psVector* newBounds = psVectorAlloc(n + 1, PS_TYPE_F32);
    newHist->bounds = newBounds;

    // Calculate the bounds for each bin.
    psF32 binSize = (upper - lower) / (psF32)n; // The histogram bin size
    // XXX: Is the following necessary? It prevents the max data point from being in a non-existant bin.
    binSize += FLT_EPSILON;
    for (long i = 0; i < n + 1; i++) {
        newBounds->data.F32[i] = lower + (binSize * (psF32)i);
    }

    // Allocate the bins, and initialize them to zero.
    newHist->nums = psVectorAlloc(n, PS_TYPE_F32);
    psVectorInit(newHist->nums, 0.0);

    // Initialize the other members.
    newHist->minNum = 0;
    newHist->maxNum = 0;
    newHist->uniform = true;

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

/******************************************************************************
psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure
with the specifed bounds.
 
Input:
    bounds
Returns:
    The histogram structure
 *****************************************************************************/
psHistogram* psHistogramAllocGeneric(const psVector* bounds)
{
    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
    PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL);
    PS_ASSERT_LONG_LARGER_THAN_OR_EQUAL(bounds->n, (long)2, NULL);

    // Allocate memory for the new histogram structure.
    psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure
    psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
    psVector* newBounds = psVectorCopy(NULL, bounds, PS_TYPE_F32);
    newHist->bounds = newBounds;

    // Allocate the bins, and initialize them to zero.  If there are N bounds,
    // then there are N-1 bins.
    newHist->nums = psVectorAlloc((bounds->n) - 1, PS_TYPE_F32);
    psVectorInit(newHist->nums, 0.0);

    // Initialize the other members.
    newHist->minNum = 0;
    newHist->maxNum = 0;
    newHist->uniform = false;

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

static void histogramFree(psHistogram* myHist)
{
    psFree(myHist->bounds);
    psFree(myHist->nums);
}


bool psMemCheckHistogram(psPtr ptr)
{
    PS_ASSERT_PTR(ptr, false);
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)histogramFree );
}



/*****************************************************************************
UpdateHistogramBins(binNum, out, data, error): This routine is to be used when
updating the histogram in the presence of errors in the input data.  We treat
the data point as a boxcar PDF and update a range of points surrounding the
histogram bin which contains the point.  The width of that boxcar is defined
as 2.35 * error.
 
XXX: Must test this.
 *****************************************************************************/
static bool UpdateHistogramBins(long binNum, // Bin number of the data point
                                psHistogram* out, // The histogram to be updated
                                psF32 data, // The data point value
                                psF32 error // The error in the data point
                               )
{
    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(out, false);
    PS_ASSERT_PTR_NON_NULL(out->bounds, false);
    PS_ASSERT_PTR_NON_NULL(out->nums, false);
    PS_ASSERT_LONG_WITHIN_RANGE(binNum, (long)0, (long)((out->nums->n)-1), false);
    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, false);
    PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0],
                                 out->bounds->data.F32[(out->bounds->n)-1], false);

    psF32 boxcarWidth = 2.35 * error;   // Width of the boxcar
    psF32 boxcarCenter = (out->bounds->data.F32[binNum] +
                          out->bounds->data.F32[binNum+1]) / 2.0; // Centre of the boxcar
    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); // Left endpoint of the boxcar for the PDF
    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); // Right endpoint of the boxcar for the PDF
    psS32 boxcarLeftBinNum = 0;         // Bin number for left endpoint
    psS32 boxcarRightBinNum = 0;        // Bin number for right endpoint

    // Determine the left endpoint of the boxcar for the PDF.
    for (long bin = binNum; bin >= 0; bin--) {
        if (out->nums->data.F32[bin] <= boxcarLeft) {
            boxcarLeftBinNum = bin;
            break;
        }
    }

    // Determine the right endpoint of the boxcar for the PDF.
    for (long bin = binNum; bin < out->nums->n; bin++) {
        if (out->nums->data.F32[bin] >= boxcarRight) {
            boxcarRightBinNum = bin;
            break;
        }
    }

    // If the boxcar fits entirely inside this bin, then simply add 1.0 to the
    // bin and return.
    if (boxcarLeftBinNum == boxcarRightBinNum) {
        out->nums->data.F32[binNum] += 1.0;
        psTrace("psLib.math", 3, "---- %s(true) end  ----\n", __func__);
        return true;
    }

    // If we get here, multiple bins must be updated.  We handle the left-most endpoint, and right-most
    // endpoints differently.
    out->nums->data.F32[boxcarLeftBinNum] +=
        (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth;

    // Loop through the center bins, if any.
    for (long bin = boxcarLeftBinNum + 1; bin < (boxcarRightBinNum - 1); bin++) {
        out->nums->data.F32[bin] +=
            (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth;
    }

    // Handle the right endpoint differently.
    out->nums->data.F32[boxcarRightBinNum]+=
        (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth;

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


/*****************************************************************************
psVectorHistogram(out, in, errors, mask, maskVal): this procedure takes as
input a preallocated and initialized histogram structure.  It fills the bins
in that histogram structure in accordance with the input data "in" and the,
possibly NULL, mask vector.
 
Inputs:
    out
    in
    mask
    maskVal
Returns:
    The histogram structure "out".
 *****************************************************************************/
psHistogram* psVectorHistogram(psHistogram* out,
                               const psVector* values,
                               const psVector* errors,
                               const psVector* mask,
                               psMaskType maskVal)
{
    psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(out, NULL);
    PS_ASSERT_VECTOR_NON_NULL(out->bounds, NULL);
    PS_ASSERT_VECTOR_TYPE(out->bounds, PS_TYPE_F32, NULL);
    PS_ASSERT_INT_NONNEGATIVE(out->bounds->n, NULL);
    PS_ASSERT_VECTOR_NON_NULL(out->nums, NULL);
    PS_ASSERT_VECTOR_TYPE(out->nums, PS_TYPE_F32, NULL);
    PS_ASSERT_INT_NONNEGATIVE(out->nums->n, NULL);
    PS_ASSERT_VECTOR_NON_NULL(values, out);
    if (mask) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, NULL);
        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
    }
    if (errors) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NULL);
        PS_ASSERT_VECTOR_TYPE(errors, values->type.type, NULL);
    }

    long binNum = 0;                    // A temporary bin number
    long numBins = out->nums->n;        // The total number of bins
    psScalar tmpScalar;
    tmpScalar.type.type = PS_TYPE_F32;

    // Convert input and errors vectors to F32 if necessary.
    psVector* inF32 = NULL;             // F32 version of input vector
    if (values->type.type == PS_TYPE_F32) {
        inF32 = psMemIncrRefCounter((psPtr)values);
    } else {
        inF32 = psVectorCopy(NULL, values, PS_TYPE_F32);
    }
    psVector* errorsF32 = NULL;         // F32 version of errors vector
    if (errors) {
        if (errors->type.type == PS_TYPE_F32) {
            errorsF32 = psMemIncrRefCounter((psPtr)errors);
        } else {
            errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32);
        }
    }

    for (long i = 0; i < inF32->n; i++) {
        // Check if this pixel is masked, and if so, skip it.
        if (!mask || (mask && (!(mask->data.U8[i] & maskVal)))) {
            if (inF32->data.F32[i] < out->bounds->data.F32[0]) {
                // If this pixel is below minimum value, count it, then skip.
                out->minNum++;
            } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) {
                // If this pixel is above maximum value, count it, then skip.
                out->maxNum++;
            } else {
                // If this is a uniform histogram, determining the correct
                // number is trivial.
                if (out->uniform == true) {
                    float binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; // Histogram bin size
                    binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize;
                    if (errorsF32) {
                        if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errorsF32->data.F32[i])) {
                            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
                                     "bins with the errors vector.\n");
                        }
                    } else {
                        // This if-statement really shouldn't be necessary.
                        // However, due to numerical lack of precision, we
                        // occasionally produce a binNum outside the range.
                        if (binNum >= out->nums->n) {
                            binNum = out->nums->n - 1;
                        }
                        (out->nums->data.F32[binNum])+= 1.0;
                    }

                } else {
                    // If this is a non-uniform histogram, determining the
                    // correct bin number requires a bit more work.
                    tmpScalar.data.F32 = inF32->data.F32[i];
                    binNum = p_psVectorBinDisect( *(psVector* *)&out->bounds, &tmpScalar);
                    if (binNum < 0) {
                        psLogMsg(__func__, PS_LOG_WARN,
                                 "WARNING: psVectorHistogram(): element outside histogram bounds.\n");
                    } else {
                        if (errorsF32 != NULL) {
                            if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) {
                                psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
                                         "bins with the errors vector.\n");
                            }
                        } else {
                            out->nums->data.F32[binNum] += 1.0;
                        }
                    }
                }
            }
        }
    }

    psFree(inF32);
    psFree(errorsF32);

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

/******************************************************************************
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).
 *****************************************************************************/
psStats* 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);
    }

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEAN) {
        if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
            stats->sampleMean = NAN;
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) {
        if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMedian() returned an error.\n");
            stats->sampleMedian = NAN;
            stats->sampleUQ = NAN;
            stats->sampleLQ = NAN;
        }
    }

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_STDEV) {
        if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n");
            stats->sampleMean = NAN;
        } else {
            vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats);
        }
    }

    // ************************************************************************
    if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE |
                          PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) {
        if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, _("Failed to calculate the specified statistic."));
            psFree(inF32);
            psFree(errorsF32);
            psTrace("psLib.math", 3,"---- %s(NULL) end  ----\n", __func__);
            return(NULL);
        }
    }

    // ************************************************************************
    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 for input psVector.\n");
            stats->clippedMean = NAN;
            stats->clippedStdev = NAN;
        }
    }

    // ************************************************************************
    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");
        }
    }

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

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_ROBUST_STDEV);
    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_ROBUST_STDEV);
    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_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_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;
}
