/** @file  psStats.c
 *  \brief basic statistical operations
 *  @ingroup Stats
 *
 *  This file will hold the definition 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->robustN50
*      stats->clippedNvalues
*      stats->binsize
*
*
*
*
*  @version $Revision: 1.149.4.1 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-11-16 22:51:04 $
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <stdarg.h>
#include <float.h>
#include <math.h>

/*****************************************************************************/
/* INCLUDE FILES                                                             */
/*****************************************************************************/
#include "psMemory.h"
#include "psImage.h"
#include "psVector.h"
#include "psTrace.h"
#include "psLogMsg.h"
#include "psError.h"
#include "psStats.h"
#include "psMinimize.h"
#include "psPolynomial.h"
#include "psConstants.h"

#include "psErrorText.h"

/*****************************************************************************/
/* DEFINE STATEMENTS                                                         */
/*****************************************************************************/
#define PS_GAUSS_WIDTH 5                // The width of the Gaussian smoothing.
// This corresponds to N in the ADD.
#define PS_CLIPPED_NUM_ITER_LB 1
#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 10

#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1]))
/*****************************************************************************/
/* TYPE DEFINITIONS                                                          */
/*****************************************************************************/
psVector* p_psConvertToF32(psVector* in);
/*****************************************************************************/
/* GLOBAL VARIABLES                                                          */
/*****************************************************************************/

// None

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

// None

/*****************************************************************************/
/* FUNCTION IMPLEMENTATION - LOCAL                                           */
/*****************************************************************************/

psBool p_psGetStatValue(const psStats* stats, psF64 *value)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    //    switch (stats->options & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE)) {
    switch (stats->options & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) {
    case PS_STAT_SAMPLE_MEAN:
        *value = stats->sampleMean;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_SAMPLE_MEDIAN:
        *value = stats->sampleMedian;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_SAMPLE_STDEV:
        *value = stats->sampleStdev;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_ROBUST_MEAN:
        *value = stats->robustMean;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_ROBUST_MEDIAN:
        *value = stats->robustMedian;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_ROBUST_MODE:
        *value = stats->robustMode;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_ROBUST_STDEV:
        *value = stats->robustStdev;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_CLIPPED_MEAN:
        *value = stats->clippedMean;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_CLIPPED_STDEV:
        *value = stats->clippedStdev;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_MAX:
        *value = stats->max;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    case PS_STAT_MIN:
        *value = stats->min;
        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
        return true;

    default:
        psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
        return false;
    }
}

/******************************************************************************
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.
 *****************************************************************************/
/******************************************************************************
p_psVectorSampleMean(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.
 *****************************************************************************/
psS32 p_psVectorSampleMean(const psVector* myVector,
                           const psVector* errors,
                           const psVector* maskVector,
                           psU32 maskVal,
                           psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);

    psS32 i = 0;                // Loop index variable
    psF32 mean = 0.0;           // The mean
    psS32 count = 0;            // # of points in this mean

    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
    // loop.
    if (errors == NULL) {
        if (stats->options & PS_STAT_USE_RANGE) {
            if (maskVector != NULL) {
                for (i = 0; i < myVector->n; i++) {
                    // Check if the data is with the specified range
                    if (!(maskVal & maskVector->data.U8[i]) &&
                            (stats->min <= myVector->data.F32[i]) &&
                            (myVector->data.F32[i] <= stats->max)) {
                        mean += myVector->data.F32[i];
                        count++;
                    }
                }
                if (count != 0) {
                    mean /= (psF32)count;
                } else {
                    mean = NAN;
                }

            } else {
                for (i = 0; i < myVector->n; i++) {
                    if ((stats->min <= myVector->data.F32[i]) &&
                            (myVector->data.F32[i] <= stats->max)) {
                        mean += myVector->data.F32[i];
                        count++;
                    }
                }
                if (count != 0) {
                    mean /= (psF32)count;
                } else {
                    mean = NAN;
                }
            }
        } else {
            if (maskVector != NULL) {
                for (i = 0; i < myVector->n; i++) {
                    if (!(maskVal & maskVector->data.U8[i])) {
                        mean += myVector->data.F32[i];
                        count++;
                    }
                }
                if (count != 0) {
                    mean /= (psF32)count;
                } else {
                    mean = NAN;
                }
            } else {
                for (i = 0; i < myVector->n; i++) {
                    mean += myVector->data.F32[i];
                }
                mean /= (psF32)myVector->n;
            }
        }
    } else {
        psF32 errorDivisor = 0.0;
        psF32 errorSqr = 0.0;
        if (stats->options & PS_STAT_USE_RANGE) {
            if (maskVector != NULL) {
                for (i = 0; i < myVector->n; i++) {
                    // Check if the data is with the specified range
                    if (!(maskVal & maskVector->data.U8[i]) &&
                            (stats->min <= myVector->data.F32[i]) &&
                            (myVector->data.F32[i] <= stats->max)) {
                        errorSqr = errors->data.F32[i]*errors->data.F32[i];
                        mean += myVector->data.F32[i]/errorSqr;
                        errorDivisor+= (1.0 / errorSqr);
                        count++;
                    }
                }
                if (count != 0) {
                    mean /= errorDivisor;
                } else {
                    mean = NAN;
                }

            } else {
                for (i = 0; i < myVector->n; i++) {
                    if ((stats->min <= myVector->data.F32[i]) &&
                            (myVector->data.F32[i] <= stats->max)) {
                        errorSqr = errors->data.F32[i]*errors->data.F32[i];
                        mean += myVector->data.F32[i]/errorSqr;
                        errorDivisor+= (1.0 / errorSqr);
                        count++;
                    }
                }
                if (count != 0) {
                    mean /= errorDivisor;
                } else {
                    mean = NAN;
                }
            }
        } else {
            if (maskVector != NULL) {
                for (i = 0; i < myVector->n; i++) {
                    if (!(maskVal & maskVector->data.U8[i])) {
                        errorSqr = errors->data.F32[i]*errors->data.F32[i];
                        mean += myVector->data.F32[i]/errorSqr;
                        errorDivisor+= (1.0 / errorSqr);
                        count++;
                    }
                }
                if (count != 0) {
                    mean /= errorDivisor;
                } else {
                    mean = NAN;
                }
            } else {
                for (i = 0; i < myVector->n; i++) {
                    errorSqr = errors->data.F32[i]*errors->data.F32[i];
                    mean += myVector->data.F32[i]/errorSqr;
                    errorDivisor+= (1.0 / errorSqr);
                }
                mean /= errorDivisor;
            }
        }
    }

    stats->sampleMean = mean;
    if (isnan(mean)) {
        psTrace(__func__, 4, "---- %s(-1) end ----\n", __func__);
        return(-1);
    } else {
        psTrace(__func__, 4, "---- %s(0) end ----\n", __func__);
        return(0);
    }

}

/******************************************************************************
p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
max of the input vector.  If there was a problem with the max calculation,
this routine sets stats->max to NAN.
 
XXX: Do we need to factor errors into it?
 *****************************************************************************/
psS32 p_psVectorMax(const psVector* myVector,
                    const psVector* maskVector,
                    psU32 maskVal,
                    psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psS32 i = 0;                // Loop index variable
    psF32 max = -PS_MAX_F32;    // The calculated maximum
    psS32 empty = true;         // Does this vector have valid elements?

    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
    if (stats->options & PS_STAT_USE_RANGE) {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    if ((myVector->data.F32[i] > max) &&
                            (stats->min <= myVector->data.F32[i]) &&
                            (myVector->data.F32[i] <= stats->max)) {
                        max = myVector->data.F32[i];
                        empty = false;
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((myVector->data.F32[i] > max) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    max = myVector->data.F32[i];
                    empty = false;
                }
            }
        }
    } else {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    if (myVector->data.F32[i] > max) {
                        max = myVector->data.F32[i];
                        empty = false;
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if (myVector->data.F32[i] > max) {
                    max = myVector->data.F32[i];
                    empty = false;
                }
            }
        }
    }
    if (empty == false) {
        stats->max = max;
    } else {
        stats->max = NAN;
        psTrace(__func__, 4, "---- %s(1) end ----\n", __func__);
        return(1);
    }
    psTrace(__func__, 4, "---- %s(0) end ----\n", __func__);
    return(0);
}

/******************************************************************************
p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
minimum of the input vector.  If there was a problem with the min calculation,
this routine sets stats->min to NAN.
 *****************************************************************************/
psS32 p_psVectorMin(const psVector* myVector,
                    const psVector* maskVector,
                    psU32 maskVal,
                    psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psS32 i = 0;                // Loop index variable
    psF32 min = PS_MAX_F32;   // The calculated maximum
    psS32 empty = true;         // Does this vector have valid elements?

    if (stats->options & PS_STAT_USE_RANGE) {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    if ((myVector->data.F32[i] < min) &&
                            (stats->min <= myVector->data.F32[i]) &&
                            (myVector->data.F32[i] <= stats->max)) {
                        min = myVector->data.F32[i];
                        empty = false;
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((myVector->data.F32[i] < min) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    min = myVector->data.F32[i];
                    empty = false;
                }
            }
        }
    } else {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    if (myVector->data.F32[i] < min) {
                        min = myVector->data.F32[i];
                        empty = false;
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if (myVector->data.F32[i] < min) {
                    min = myVector->data.F32[i];
                    empty = false;
                }
            }
        }
    }

    if (empty == false) {
        stats->min = min;
    } else {
        stats->min = NAN;
        psTrace(__func__, 4, "---- %s(1) end ----\n", __func__);
        return(1);
    }
    psTrace(__func__, 4, "---- %s(0) end ----\n", __func__);
    return(0);
}

/******************************************************************************
p_psVectorNValues(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".
 *****************************************************************************/
bool p_psVectorCheckNonEmpty(const psVector* myVector,
                             const psVector* maskVector,
                             psU32 maskVal,
                             psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(myVector, -1);
    PS_ASSERT_PTR_NON_NULL(stats, -1);
    psS32 i = 0;                // Loop index variable

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


/******************************************************************************
p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
number of non-masked pixels in the vector that fall within the min/max
range, if specified.
 *****************************************************************************/
psS32 p_psVectorNValues(const psVector* myVector,
                        const psVector* maskVector,
                        psU32 maskVal,
                        psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(myVector, -1);
    PS_ASSERT_PTR_NON_NULL(stats, -1);
    psS32 i = 0;                // Loop index variable
    psS32 numData = 0;          // The number of data points

    if (stats->options & PS_STAT_USE_RANGE) {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    numData++;
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    numData++;
                }
            }
        }
    } else {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    numData++;
                }
            }
        } else {
            numData = myVector->n;
        }
    }

    psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, numData);
    return (numData);
}

/******************************************************************************
p_psVectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the
median of the input vector.  Returns true on success (including if there were
no valid input vector elements).
 
XXX: Use static vectors for sort arrays.
 *****************************************************************************/
bool p_psVectorSampleMedian(const psVector* myVector,
                            const psVector* maskVector,
                            psU32 maskVal,
                            psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psVector* unsortedVector = NULL;    // Temporary vector
    psVector* sortedVector = NULL;      // Temporary vector
    psS32 i = 0;                        // Loop index variable
    psS32 count = 0;
    psS32 nValues = 0;

    // Determine how many data points fit inside this min/max range
    // and are not masked, if the maskVector is not NULL.
    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);

    // XXX: Is a warning message appropriate?  What value should we set the
    // sampleMedian to?  Should we generate an error?
    if (nValues <= 0) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n");
        return(true);
        stats->sampleMedian = NAN;
    }

    // Allocate temporary vectors for the data.
    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);

    // 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.
        count = 0;
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    unsortedVector->data.F32[count] = myVector->data.F32[i];
                    count++;
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    unsortedVector->data.F32[count] = myVector->data.F32[i];
                    count++;
                }
            }
        }
    } else {
        // Store all non-masked data points into the temporary vectors.
        count = 0;
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    unsortedVector->data.F32[count] = myVector->data.F32[i];
                    count++;
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                unsortedVector->data.F32[i] = myVector->data.F32[i];
            }
        }
    }
    // Sort the temporary vectors.
    sortedVector = psVectorSort(sortedVector, unsortedVector);
    if (sortedVector == NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                false,
                PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM);
        psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
        return(false);
    }

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

    // Free the temporary data structures.
    psFree(unsortedVector);
    psFree(sortedVector);

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

/******************************************************************************
p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
smoothed data.
 
XXX: Only PS_TYPE_F32 is supported.
 
XXX: Write a general routine which smoothes a psVector.  This routine should
call that.  Is that possible?
 *****************************************************************************/
psVector *p_psVectorSmoothHistGaussian(psHistogram *histogram,
                                       psF32 sigma)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psTrace(__func__, 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(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(histogram->nums);
    }

    psS32 numBins = histogram->nums->n;
    psS32 numBounds = histogram->bounds->n;
    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32);
    psF32 firstBound = histogram->bounds->data.F32[0];
    psF32 lastBound = histogram->bounds->data.F32[numBounds-1];
    psScalar x;
    x.type.type = PS_TYPE_F32;
    psS32 jMin = 0;
    psS32 jMax = 0;

    if (histogram->uniform == false) {
        //
        // We get here if the histogram is non-uniform.
        //

        for (psS32 i = 0; i < numBins; i++) {
            // Determine the midpoint of bin i.
            psS32 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
            //
            x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
            if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
                jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x);
                if (jMin < 0) {
                    psError(PS_ERR_UNEXPECTED_NULL,
                            false,
                            PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
                    psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
                    return(NULL);
                }
            } else if (x.data.F32 <= firstBound) {
                jMin = 0;
            } else if (x.data.F32 >= lastBound) {
                jMin = histogram->bounds->n - 1;
            }

            x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
            if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
                jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x);
                if (jMax < 0) {
                    psError(PS_ERR_UNEXPECTED_NULL,
                            false,
                            PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
                    psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
                    return(NULL);
                }
            } else if (x.data.F32 <= firstBound) {
                jMax = 0;
            } else if (x.data.F32 >= lastBound) {
                jMax = histogram->bounds->n - 1;
            }

            //
            // Loop from jMin to jMax, computing the gaussian of data i.
            //
            smooth->data.F32[i] = 0.0;
            for (psS32 j = jMin ; j <= jMax ; j++) {
                psS32 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 (psS32 i = 0; i < numBins; i++) {
            psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0];
            psS32 gaussWidth = (psS32) ((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.
            //
            gaussWidth = 10;
            //
            // 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 = i - gaussWidth;
            if (jMin < 0 ) {
                jMin = 0;
            }
            psS32 jMax = i + gaussWidth;
            if (jMax > (histogram->bounds->n - 1)) {
                jMax = (histogram->bounds->n - 1);
            }

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

    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(smooth);
    }
    psTrace(__func__, 4, "---- %s(psVector) end ----\n", __func__);
    return(smooth);
}
/******************************************************************************
p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates
the upper and/or lower quartiles of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 
XXX: Use static vectors.
 *****************************************************************************/
bool p_psVectorSampleQuartiles(const psVector* myVector,
                               const psVector* maskVector,
                               psU32 maskVal,
                               psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psVector* unsortedVector = NULL;    // Temporary vector
    psVector* sortedVector = NULL;      // Temporary vector
    psS32 i = 0;                  // Loop index variable
    psS32 count = 0;              // # of points in this mean.
    psS32 nValues = 0;            // # data points

    // Determine how many data points fit inside this min/max range
    // and are not maxed, IF the maskVector is not NULL.
    nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);

    // Allocate temporary vectors for the data.
    unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);

    // 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.
        count = 0;
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) {
                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                }
            }
        }
    } else {
        // Store all non-masked data points into the temporary vectors.
        count = 0;
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                unsortedVector->data.F32[i] = myVector->data.F32[i];
            }
        }
    }

    // Sort the temporary vectors.
    sortedVector = psVectorSort(sortedVector, unsortedVector);
    if (sortedVector == NULL) {
        psError(PS_ERR_UNEXPECTED_NULL,
                false,
                PS_ERRORTEXT_psStats_STATS_SAMPLE_MEDIAN_SORT_PROBLEM);
        psTrace(__func__, 4, "---- %s(false) end ----\n", __func__);
        return(false);
    }

    // Calculate the quartile points exactly.
    // XXX: We should probably do a bin-midpoint if the quartile is not an
    // integer.
    stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)];
    stats->sampleLQ = sortedVector->data.F32[nValues / 4];

    // Free the temporary data structures.
    psFree(unsortedVector);
    psFree(sortedVector);
    psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
    return(true);
}

/******************************************************************************
p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
stdev of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 
 *****************************************************************************/
void p_psVectorSampleStdevOLD(const psVector* myVector,
                              const psVector* errors,
                              const psVector* maskVector,
                              psU32 maskVal,
                              psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psS32 i = 0;                  // Loop index variable
    psS32 countInt = 0;           // # of data points being used
    psF32 countFloat = 0.0;     // # of data points being used
    psF32 mean = 0.0;           // The mean
    psF32 diff = 0.0;           // Used in calculating stdev
    psF32 sumSquares = 0.0;     // temporary variable
    psF32 sumDiffs = 0.0;       // temporary variable

    // This procedure requires the mean.  If it has not been already
    // calculated, then call p_psVectorSampleMean()
    if (0 != isnan(stats->sampleMean)) {
        p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats);
    }
    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    if (0 != isnan(stats->sampleMean)) {
        stats->sampleStdev = NAN;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n");
        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
        return;
    }

    mean = stats->sampleMean;
    if (stats->options & PS_STAT_USE_RANGE) {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    diff = myVector->data.F32[i] - mean;
                    sumSquares += (diff * diff);
                    sumDiffs += diff;
                    countInt++;
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    diff = myVector->data.F32[i] - mean;
                    sumSquares += (diff * diff);
                    sumDiffs += diff;
                    countInt++;
                }
            }
            countInt = myVector->n;
        }
    } else {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    diff = myVector->data.F32[i] - mean;
                    sumSquares += (diff * diff);
                    sumDiffs += diff;
                    countInt++;
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                diff = myVector->data.F32[i] - mean;
                sumSquares += (diff * diff);
                sumDiffs += diff;
                countInt++;
            }
            countInt = myVector->n;
        }
    }
    if (countInt == 0) {
        stats->sampleStdev = NAN;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d).  Setting stats->sampleStdev = NAN.\n", countInt);
    } else if (countInt == 1) {
        stats->sampleStdev = 0.0;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
    } else {
        countFloat = (psF32)countInt;
        stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    }
    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
}
/******************************************************************************
p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
stdev of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 
 *****************************************************************************/
void p_psVectorSampleStdev(const psVector* myVector,
                           const psVector* errors,
                           const psVector* maskVector,
                           psU32 maskVal,
                           psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psS32 i = 0;                  // Loop index variable
    psS32 countInt = 0;           // # of data points being used
    psF32 countFloat = 0.0;     // # of data points being used
    psF32 mean = 0.0;           // The mean
    psF32 diff = 0.0;           // Used in calculating stdev
    psF32 sumSquares = 0.0;     // temporary variable
    psF32 sumDiffs = 0.0;       // temporary variable
    //    psF32 sum1;
    //    psF32 sum2;
    psF32 errorDivisor = 0.0f;

    // This procedure requires the mean.  If it has not been already
    // calculated, then call p_psVectorSampleMean()
    if (0 != isnan(stats->sampleMean)) {
        p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats);
    }
    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    if (0 != isnan(stats->sampleMean)) {
        stats->sampleStdev = NAN;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n");
        psTrace(__func__, 4, "---- %s() end ----\n", __func__);
        return;
    }

    mean = stats->sampleMean;
    if (stats->options & PS_STAT_USE_RANGE) {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    diff = myVector->data.F32[i] - mean;
                    sumSquares += (diff * diff);
                    sumDiffs += diff;
                    countInt++;
                    if (errors != NULL) {
                        errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    diff = myVector->data.F32[i] - mean;
                    sumSquares += (diff * diff);
                    sumDiffs += diff;
                    countInt++;
                    if (errors != NULL) {
                        errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
                    }
                }
            }
            countInt = myVector->n;
        }
    } else {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    diff = myVector->data.F32[i] - mean;
                    sumSquares += (diff * diff);
                    sumDiffs += diff;
                    countInt++;
                    if (errors != NULL) {
                        errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                diff = myVector->data.F32[i] - mean;
                sumSquares += (diff * diff);
                sumDiffs += diff;
                countInt++;
            }
            countInt = myVector->n;
            if (errors != NULL) {
                errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i]));
            }
        }
    }

    if (countInt == 0) {
        stats->sampleStdev = NAN;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d).  Setting stats->sampleStdev = NAN.\n", countInt);
    } else if (countInt == 1) {
        stats->sampleStdev = 0.0;
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d).  Setting stats->sampleStdev = 0.0.\n", countInt);
    } else {
        // XXX: The ADD specifies this as the definition of the standard
        // deviation if the errors are known.  Verify this with IfA: none of
        // the data points in the vector are used.  Verify that the masks and
        // data ranges are used correctly.
        if (errors != NULL) {
            stats->sampleStdev = (1.0 / sqrtf(errorDivisor));
        } else {
            countFloat = (psF32)countInt;
            stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));

        }
    }
    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
}

/******************************************************************************
p_psVectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the
clipped stats (mean or stdev) of the input vector.
 
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    0 for success.
    -1: error
    -2: warning
 
XXX: Do we really need to calculate median and mean?
 
XXX: Use static vectors for tmpMask.
 
XXX: This has not been tested.
 *****************************************************************************/
psS32 p_psVectorClippedStats(const psVector* myVector,
                             const psVector* errors,
                             const psVector* maskVector,
                             psU32 maskVal,
                             psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psF32 clippedMean = 0.0;    // self-explanatory
    psF32 clippedStdev = 0.0;   // self-explanatory
    psVector* tmpMask = NULL;   // Temporary vector for masks during iterations.
    static psStats *statsTmp = NULL;   // Temporary psStats struct.
    psS32 rc = 0;               // Return code.

    // 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_INT_WITHIN_RANGE(stats->clipSigma,
                               PS_CLIPPED_SIGMA_LB,
                               PS_CLIPPED_SIGMA_UB, -1);

    // Allocate a psStats structure for calculating the mean, median, and
    // stdev.
    if (statsTmp == NULL) {
        statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
        p_psMemSetPersistent(statsTmp, true);
    } else {
        // XXX EAM : initialize structure if already allocated
        statsTmp->sampleMean = NAN;
        statsTmp->sampleMedian = NAN;
        statsTmp->sampleStdev = NAN;
        statsTmp->sampleUQ = NAN;
        statsTmp->sampleLQ = NAN;
        statsTmp->robustMean = NAN;
        statsTmp->robustMedian = NAN;
        statsTmp->robustMode = NAN;
        statsTmp->robustStdev = NAN;
        statsTmp->robustUQ = NAN;
        statsTmp->robustLQ = NAN;
        statsTmp->robustN50 = -1;            // XXX: This is never used
        statsTmp->robustNfit = -1;
        statsTmp->clippedMean = NAN;
        statsTmp->clippedStdev = NAN;
        statsTmp->clippedNvalues = -1;     // XXX: This is never used
        statsTmp->clipSigma = 3.0;
        statsTmp->clipIter = 3;
        statsTmp->min = NAN;
        statsTmp->max = NAN;
        statsTmp->binsize = NAN;          // XXX: This is never used
    }

    // We allocate a temporary mask vector since during the iterative
    // steps that follow, we will be masking off additional data points.
    // However, we do no want to modify the original mask vector.
    tmpMask = psVectorAlloc(myVector->n, PS_TYPE_U8);

    // If we were called with a mask vector, then initialize the temporary
    // mask vector with those values.  Otherwise, initialize to zero.
    if (maskVector != NULL) {
        for (psS32 i = 0; i < tmpMask->n; i++) {
            tmpMask->data.U8[i] = maskVector->data.U8[i];
        }
    } else {
        for (psS32 i = 0; i < tmpMask->n; i++) {
            tmpMask->data.U8[i] = 0;
        }
    }

    // 1. Compute the sample median.
    p_psVectorSampleMedian(myVector, maskVector, maskVal, statsTmp);
    if (isnan(statsTmp->sampleMedian)) {
        psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleMedian returned NAN\n");
        stats->clippedMean = NAN;
        stats->clippedStdev = NAN;
        psTrace(__func__, 4, "---- %s(-2) end ----\n", __func__);
        return(-2);
    }

    // 2. Compute the sample standard deviation.
    p_psVectorSampleStdev(myVector, errors, maskVector, maskVal, statsTmp);
    if (isnan(statsTmp->sampleStdev)) {
        psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleStdev returned NAN\n");
        stats->clippedMean = NAN;
        stats->clippedStdev = NAN;
        psTrace(__func__, 4, "---- %s(-2) end ----\n", __func__);
        return(-2);
    }

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

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

    // 5. Repeat N (stats->clipIter) times:
    for (psS32 iter = 0; iter < stats->clipIter; iter++) {
        // a) Exclude all values x_i for which |x_i - x| > K * stdev
        if (errors != NULL) {
            for (psS32 j = 0; j < myVector->n; j++) {
                if (fabs(myVector->data.F32[j] - clippedMean) >
                        (stats->clipSigma * errors->data.F32[j])) {
                    tmpMask->data.U8[j] = 0xff;
                }
            }
        } else {
            for (psS32 j = 0; j < myVector->n; j++) {
                if (fabs(myVector->data.F32[j] - clippedMean) >
                        (stats->clipSigma * clippedStdev)) {
                    tmpMask->data.U8[j] = 0xff;
                }
            }
        }

        // b) compute new mean and stdev
        p_psVectorSampleMean(myVector, errors, tmpMask, 0xff, statsTmp);
        p_psVectorSampleStdev(myVector, errors, tmpMask, 0xff, statsTmp);

        // 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.  XXX: Should we throw an error/warning here?
            iter = stats->clipIter;
            rc = -1;
        } else {
            clippedMean = statsTmp->sampleMean;
            clippedStdev = statsTmp->sampleStdev;
        }
    }

    // 7. The last calcuated value of x is the cliped mean.
    if (stats->options & PS_STAT_CLIPPED_MEAN) {
        stats->clippedMean = clippedMean;
    }
    // 8. The last calcuated value of stdev is the cliped stdev.
    if (stats->options & PS_STAT_CLIPPED_STDEV) {
        stats->clippedStdev = clippedStdev;
    }

    psFree(tmpMask);
    psTrace(__func__, 4, "---- %s(%d) end ----\n", __func__, rc);
    return(rc);
}

/*****************************************************************************
These macros and functions define the following functions:
 
    p_psNormalizeVectorRange(myData, low, high)
 
That assumes that the low/high arguments are PS_TYPE_F64; the vector myData
can be of any type.  Arguments low/high will be converted to the appropriate
type and one of the type-specific functions below will be called:
 
    p_psNormalizeVectorRangeU8(myData, low, high)
    p_psNormalizeVectorRangeU16(myData, low, high)
    p_psNormalizeVectorRangeU32(myData, low, high)
    p_psNormalizeVectorRangeU64(myData, low, high)
    p_psNormalizeVectorRangeS8(myData, low, high)
    p_psNormalizeVectorRangeS16(myData, low, high)
    p_psNormalizeVectorRangeS32(myData, low, high)
    p_psNormalizeVectorRangeS64(myData, low, high)
    p_psNormalizeVectorRangeF32(myData, low, high)
    p_psNormalizeVectorRangeF64(myData, low, high)
 *****************************************************************************/
#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \
void p_psNormalizeVectorRange##TYPE(psVector* myData, \
                                    ps##TYPE outLow, \
                                    ps##TYPE outHigh) \
{ \
    ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \
    ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \
    psS32 i = 0; \
    \
    for (i = 0; i < myData->n; i++) { \
        if (myData->data.TYPE[i] < min) { \
            min = myData->data.TYPE[i]; \
        } \
        if (myData->data.TYPE[i] > max) { \
            max = myData->data.TYPE[i]; \
        } \
    } \
    \
    /* Ensure that max!=min before we divide by (max-min) */ \
    if (max != min) { \
        for (i = 0; i < myData->n; i++) { \
            myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \
                                    (outHigh - outLow) / (max - min)); \
        } \
    } else { \
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min).  Setting all elements to min.\n"); \
        for (i = 0; i < myData->n; i++) \
        { \
            \
            myData->data.TYPE[i] = outLow; \
            \
        } \
    } \
} \

PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32)
PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64)

void p_psNormalizeVectorRange(psVector* myData,
                              psF64 outLow,
                              psF64 outHigh)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    switch (myData->type.type) {
    case PS_TYPE_U8:
        p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh);
        break;
    case PS_TYPE_U16:
        p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh);
        break;
    case PS_TYPE_U32:
        p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh);
        break;
    case PS_TYPE_U64:
        p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh);
        break;
    case PS_TYPE_S8:
        p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh);
        break;
    case PS_TYPE_S16:
        p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh);
        break;
    case PS_TYPE_S32:
        p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh);
        break;
    case PS_TYPE_S64:
        p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh);
        break;
    case PS_TYPE_F32:
        p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh);
        break;
    case PS_TYPE_F64:
        p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh);
        break;
    case PS_TYPE_C32:
    case PS_TYPE_C64:
    default:
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                "Unallowable operation: %s has incorrect type.",
                myData);
        break;
    }
    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
}

/******************************************************************************
p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, getThisValue): This routine
takes as input a 1-D polynomial of arbitrary order and a range of x-values for
which it is defined:  [rangeLow, rangeHigh].  It determines the x-value of
that polynomial such that f(x) == getThisValue.  This function uses a
binary-search algorithm on the range and assumes that the polynomial is
monotonically increasing or decreasing within that range.
 
XXX: Terminate when f(x)-getThisValue is within some error tolerance.
 
XXX: Create a 2nd-order polynomial version and solve for X analytically.
 *****************************************************************************/
psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly,
                       psF32 rangeLow,
                       psF32 rangeHigh,
                       psF32 getThisValue)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    PS_ASSERT_POLY_NON_NULL(myPoly, NAN);
    PS_ASSERT_FLOAT_LARGER_THAN(rangeHigh, rangeLow, NAN);
    // We ensure that the requested f(y) value, which is getThisValue, is
    // falls within the range of y-values of the polynomial "myPoly" in the
    // specified x-range (rangeLow:rangeHigh).
    psF32 fLo = psPolynomial1DEval(myPoly, rangeLow);
    psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh);
    if (!((fLo <= getThisValue) && (fHi >= getThisValue))) {
        psError(PS_ERR_UNKNOWN, true, "The requested y value (%f) does not fall within the specified range (%f to %f)\n", getThisValue, fLo, fHi);
        psTrace(__func__, 4, "---- %s(NAN) end ----\n", __func__);
        return(NAN);
    }

    psS32 numIterations = 0;
    psF32 midpoint = 0.0;
    psF32 oldMidpoint = 1.0;
    psF32 f = 0.0;

    while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) {
        midpoint = (rangeHigh + rangeLow) / 2.0;
        if (fabs(midpoint - oldMidpoint) <= FLT_EPSILON) {
            return (midpoint);
        }
        oldMidpoint = midpoint;

        f = psPolynomial1DEval(myPoly, midpoint);
        if (fabs(f - getThisValue) <= FLT_EPSILON) {
            return (midpoint);
        }

        if (f > getThisValue) {
            rangeHigh = midpoint;
        } else {
            rangeLow = midpoint;
        }
        numIterations++;
    }
    psTrace(__func__, 4, "---- %s(%f) end ----\n", __func__, midpoint);
    return (midpoint);
}

/******************************************************************************
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).
 
XXX: After you fit the polynomial, solve for X analytically.
 
XXX: the vectors do not have to be the same length.  Must insert the proper
tests to ensure that binNum is within acceptable ranges for both vectors.
 
XXX: This currently assumes that the three points are monotonically increasing
or decreasing: so, it works for the cumulative histogram vectors, but not for
arbitrary vectors.  We should probably test that condition.
*****************************************************************************/
psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
                                        psVector *yVec,
                                        psS32 binNum,
                                        psF32 yVal)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psTrace(__func__, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
    if (psTraceGetLevel(__func__) >= 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, (xVec->n - 1), NAN);
    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (yVec->n - 1), NAN);

    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
    psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);

    psF32 tmpFloat = 0.0f;

    if ((binNum > 0) && (binNum < (yVec->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(__func__, 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(__func__, 6, "x vec is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
        psTrace(__func__, 6, "y vec is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);

        //
        // Ensure that the y values are monotonic.
        //
        // XXX: This routine should probably be rewritten in a more general fashion
        // so that the folloiwng checks are not necessary.
        //
        if (y->data.F64[0] < y->data.F64[1]) {
            if (!(y->data.F64[1] <= y->data.F64[2])) {
                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
                psFree(myPoly);
                psFree(x);
                psFree(y);
                psFree(yErr);
                psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
                return(NAN);
            }
            // Ensure that yVal is within the range of the bins we are using.
            if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) {
                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                        PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
                        (psF64)yVal, y->data.F64[2], y->data.F64[0]);
            }
        } else if (y->data.F64[0] > y->data.F64[1]) {
            if (!(y->data.F64[1] >= y->data.F64[2])) {
                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
                psFree(myPoly);
                psFree(x);
                psFree(y);
                psFree(yErr);
                psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
                return(NAN);
            }
            // Ensure that yVal is within the range of the bins we are using.
            if (!((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) {
                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                        PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
                        (psF64)yVal, y->data.F64[2], y->data.F64[0]);
            }
        }

        yErr->data.F64[0] = 1.0;
        yErr->data.F64[1] = 1.0;
        yErr->data.F64[2] = 1.0;

        // Determine the coefficients of the polynomial.
        //        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);
        if (myPoly == NULL) {
            psError(PS_ERR_UNEXPECTED_NULL,
                    false,
                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT);
            psFree(myPoly);
            psFree(x);
            psFree(y);
            psFree(yErr);
            psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
            return(NAN);
        }
        psTrace(__func__, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
        psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
        psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
        psTrace(__func__, 6, "Fitted y vec is (%.2f %.2f %.2f)\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]));

        // Call p_ps1DPolyMedian(), which does a binary search on the
        // polynomial, looking for the value x such that f(x) = yVal
        psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
        tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal);

        if (isnan(tmpFloat)) {
            psError(PS_ERR_UNEXPECTED_NULL,
                    false,
                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
            psFree(myPoly);
            psFree(x);
            psFree(y);
            psFree(yErr);
            psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
            return(NAN);
        }

    } else {
        // The special case where we have two points only at the beginning of
        // the vectors x and y.
        if (binNum == 0) {
            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(__func__, 6, "FIT: return %f\n", tmpFloat);
    psFree(myPoly);
    psFree(x);
    psFree(y);
    psFree(yErr);

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

/*****************************************************************************
XXX: Is there a psLib function for this?
 *****************************************************************************/
psVector *PsVectorDup(psVector *in)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(in, NULL);
    psVector *out = psVectorAlloc(in->n, in->type.type);

    if (in->type.type == PS_TYPE_F32) {
        for (psS32 i = 0 ; i < in->n ; i++) {
            out->data.F32[i] = in->data.F32[i];
        }
    } else if (in->type.type == PS_TYPE_F64) {
        for (psS32 i = 0 ; i < in->n ; i++) {
            out->data.F64[i] = in->data.F64[i];
        }
    } else {
        psError(PS_ERR_UNKNOWN, true, "Unallowable vector type.\n");
        psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
        return(NULL);
    }
    psTrace(__func__, 4, "---- %s(psVector) end ----\n", __func__);
    return(out);
}

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

    psF32 x = coords->data.F32[0];
    psF32 mean = params->data.F32[0];
    psF32 stdev = params->data.F32[1];

    if (deriv == NULL) {
        deriv = psVectorAlloc(2, PS_TYPE_F32);
    } else {
        PS_ASSERT_VECTOR_SIZE(deriv, 2, NAN);
        PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN);
    }

    psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false);
    deriv->data.F32[0] = tmp / (stdev * stdev);
    tmp = (x - mean) * (x - mean) * psGaussian(x, mean, stdev, 0);
    deriv->data.F32[1] = tmp / (stdev * stdev * stdev);

    // printf("f(x, mean, stdev) = f(%.2f %.2f %.2f) is %.2f\n", x, mean, stdev, psGaussian(x, mean, stdev, false));

    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    return(psGaussian(x, mean, stdev, false));
}

psF32 LinInterpolate(psF32 x0,
                     psF32 x1,
                     psF32 y0,
                     psF32 y1,
                     psF32 y)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psTrace(__func__, 5, "LinInterpolate(%.2f %.2f %.2f %.2f %.2f)\n", x0, x1, y0, y1, y);
    if (y0 == y1) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: y0 == y1.  Cannot interpolate.  Returning NaN.\n");
        return(NAN);
    }

    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    return(x0 + y * (x1 - x0) / (y1 - y0));
}

// XXX: Get rid of this, or put it in psConstants.h
void PRINT_LEAKS()
{
    psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
    psMemCheckCorruption( 1 );
    psS32 currentId = psMemGetId();
    //    psMemBlock** blks;
    //    psS32 nLeaks = psMemCheckLeaks( currentId, &blks, stderr, false );
    psS32 nLeaks = psMemCheckLeaks(currentId,NULL,stderr,false);
    printf("COOL: %d leaks.\n", nLeaks);
    tmpScalar->data.F32 = 0.0;
}

/*
#define PRINT_MEMLEAKS(NUM) \
    printf("A: ---------------------------------------- (ID: %d)\n", NUM); \
    memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false); \
    printf("B: ---------------------------------------- (%d)\n", memLeaks); \
*/
#define PRINT_MEMLEAKS(NUM) \
memLeaks = currentId;


/******************************************************************************
p_psVectorRobustStats(myVector, maskVector, maskVal, 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.
*****************************************************************************/
psS32 p_psVectorRobustStats(const psVector* myVector,
                            const psVector* errors,
                            const psVector* maskVector,
                            psU32 maskVal,
                            psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(myVector);
    }
    //    psS32 currentId = psMemGetId();
    //    psS32 memLeaks = 0;
    psHistogram *robustHistogram = NULL;
    psHistogram *cumulativeRobustHistogram = NULL;
    psS32 numBins = 0;
    psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32);
    tmpScalar->type.type = PS_TYPE_F32;
    psS32 totalDataPoints = 0;
    psS32 rc = 0;
    psS32 rcBool = false;
    psVector *tmpMaskVec = psVectorAlloc(myVector->n, PS_TYPE_U8);
    if (maskVector != NULL) {
        for (psS32 i = 0 ; i < myVector->n ; i++) {
            tmpMaskVec->data.U8[i] = maskVector->data.U8[i];
        }
    } else {
        for (psS32 i = 0 ; i < myVector->n ; i++) {
            tmpMaskVec->data.U8[i] = 0;
        }
    }
    psBool iterate = true;
    psF32 sigma;
    psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX);

    while (iterate) {
        psTrace(__func__, 6, "Iterating on Bin size.\n");
        //
        // Determine the bin size of the robust histogram.  This is done
        // by computing the total range of data values and dividing by 1000.0.
        //
        rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax);
        rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax);
        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
            psFree(tmpStatsMinMax);
            psFree(tmpMaskVec);
            psFree(tmpScalar);
            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
            return(1);
        }
        psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
        psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f;
        psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize);

        //
        // If all data points have the same value, then we set the appropiate
        // members of stats and return.
        //
        if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) {
            psTrace(__func__, 5, "All data points have the same value.\n", binSize);
            if (stats->options & PS_STAT_ROBUST_MEDIAN) {
                stats->robustMedian = tmpStatsMinMax->min;
            }
            if (stats->options & PS_STAT_ROBUST_QUARTILE) {
                stats->robustUQ = tmpStatsMinMax->min;
                stats->robustLQ = tmpStatsMinMax->min;
            }
            // XXX: Set these to the number of unmasked data points?
            stats->robustNfit = 0.0;
            stats->robustN50 = 0.0;
            psFree(tmpStatsMinMax);
            psFree(tmpMaskVec);
            psFree(tmpScalar);

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

        //
        // 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.
        //
        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize);
        psTrace(__func__, 6, "Numbins is %d\n", numBins);
        robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
        cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);

        // Populate the histogram array.
        robustHistogram = psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal);
        if (psTraceGetLevel(__func__) >= 8) {
            PS_VECTOR_PRINT_F32(robustHistogram->bounds);
            PS_VECTOR_PRINT_F32(robustHistogram->nums);
        }
        //
        // ADD: Step 1.
        // Construct the cumulative histogram from the specific histogram
        //
        cumulativeRobustHistogram->nums->data.F32[0] = robustHistogram->nums->data.F32[0];
        for (psS32 i = 1 ; i < robustHistogram->nums->n ; i++) {
            cumulativeRobustHistogram->nums->data.F32[i] = cumulativeRobustHistogram->nums->data.F32[i-1] +
                    robustHistogram->nums->data.F32[i];
        }
        if (psTraceGetLevel(__func__) >= 8) {
            PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->bounds);
            PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->nums);
        }

        //
        // ADD: Step 2.
        // Find the bin which contains the 50% data point.
        //
        totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1];
        psTrace(__func__, 6, "Total data points is %d\n", totalDataPoints);
        psS32 binMedian;
        if (totalDataPoints/2.0 < cumulativeRobustHistogram->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 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
            if (binMedian < 0) {
                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian);
                psFree(tmpStatsMinMax);
                psFree(robustHistogram);
                psFree(cumulativeRobustHistogram);
                psFree(tmpScalar);
                psFree(tmpMaskVec);
                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
                return(1);
            }
        }
        psTrace(__func__, 6, "The median bin is %d\n", binMedian);

        //
        // ADD: Step 3.
        // Interpolate to the exact 50% position: this is the robust histogram median.
        // XXX: Check for errors here!
        //
        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
                                  *(psVector* *)&cumulativeRobustHistogram->bounds,
                                  *(psVector* *)&cumulativeRobustHistogram->nums,
                                  binMedian,
                                  totalDataPoints/2.0);
        psTrace(__func__, 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.
        //
        psS32 binLo = 0;
        psS32 binHi = 0;
        tmpScalar->data.F32 = totalDataPoints * 0.158655f;
        if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
            binLo = 0;
        } else {
            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
            if (binLo > cumulativeRobustHistogram->nums->n-1) {
                binLo = cumulativeRobustHistogram->nums->n-1;
            }
        }
        tmpScalar->data.F32 = totalDataPoints * 0.841345f;
        if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
            binHi = 0;
        } else {
            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
            binHi = 1+ p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
            if (binHi > cumulativeRobustHistogram->nums->n-1) {
                binHi = cumulativeRobustHistogram->nums->n-1;
            }
        }
        psTrace(__func__, 6, "The 15.8655% and 84.1345% data point bins are (%d, %d).\n", binLo, binHi);
        psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binLo));
        psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binHi));

        if ((binLo < 0) || (binHi < 0)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655 and 84.1345 percent data point\n");
            psFree(tmpStatsMinMax);
            psFree(robustHistogram);
            psFree(cumulativeRobustHistogram);
            psFree(tmpScalar);
            psFree(tmpMaskVec);
            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
            return(1);
        }

        //
        // ADD: Step 4b.
        // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
        //
        // 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.
        //
        psF32 binLoF32;
        psF32 binHiF32;
        if (0) {
            binLoF32 = fitQuadraticSearchForYThenReturnX(
                           *(psVector* *)&cumulativeRobustHistogram->bounds,
                           *(psVector* *)&cumulativeRobustHistogram->nums,
                           binLo,
                           totalDataPoints * 0.158655f);
            binHiF32 = fitQuadraticSearchForYThenReturnX(
                           *(psVector* *)&cumulativeRobustHistogram->bounds,
                           *(psVector* *)&cumulativeRobustHistogram->nums,
                           binHi,
                           totalDataPoints * 0.841345f);
            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32);
        }

        if (1) {
            psTrace(__func__, 6, "binLo is %d.  Nums at that bin and the next are (%.2f, %.2f)\n", binLo, cumulativeRobustHistogram->nums->data.F32[binLo], cumulativeRobustHistogram->nums->data.F32[binLo+1]);
            psTrace(__func__, 6, "binHi is %d.  Nums at that bin and the next are (%.2f, %.2f)\n", binHi, cumulativeRobustHistogram->nums->data.F32[binHi], cumulativeRobustHistogram->nums->data.F32[binHi+1]);
            // XXX: Add checks to ensure that these bins are not the first (seg fault).
            binLoF32 = LinInterpolate(
                           cumulativeRobustHistogram->bounds->data.F32[binLo],
                           cumulativeRobustHistogram->bounds->data.F32[binLo+1],
                           cumulativeRobustHistogram->nums->data.F32[binLo-1],
                           cumulativeRobustHistogram->nums->data.F32[binLo],
                           totalDataPoints * 0.158655f);
            binHiF32 = LinInterpolate(
                           cumulativeRobustHistogram->bounds->data.F32[binHi],
                           cumulativeRobustHistogram->bounds->data.F32[binHi+1],
                           cumulativeRobustHistogram->nums->data.F32[binHi-1],
                           cumulativeRobustHistogram->nums->data.F32[binHi],
                           totalDataPoints * 0.841345f);
            // XXX: Check for NANs
            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32);
        }

        //
        // ADD: Step 5.
        // Determine SIGMA as 1/2 of the distance between these positions.
        //
        sigma = (binHiF32 - binLoF32) / 2.0;
        psTrace(__func__, 6, "The current sigma is %f.\n", 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(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
            psF32 medianLo = robustHistogram->bounds->data.F32[binMedian - 25];
            psF32 medianHi = robustHistogram->bounds->data.F32[binMedian + 25];
            psTrace(__func__, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
            for (psS32 i = 0 ; i < myVector->n ; i++) {
                if ((myVector->data.F32[i] < medianLo) ||
                        (myVector->data.F32[i] > medianHi)) {
                    tmpMaskVec->data.U8[i] = 1;
                    psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]);
                }
            }
            psFree(robustHistogram);
            psFree(cumulativeRobustHistogram);

        } else {
            psTrace(__func__, 6, "*************: No more iteration.  sigma is %f\n", sigma);
            iterate = false;
        }
    }

    //
    // ADD: Step 7.
    // Find the bins which contains the 25% and 75% data points.
    //
    psS32 binLo25 = 0;
    psS32 binHi25 = 0;

    tmpScalar->data.F32 = totalDataPoints * 0.25f;
    if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
        // XXX: Special case where its in first bin.  Must code last bin.
        binLo25 = 0;
    } else {
        binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
    }
    tmpScalar->data.F32 = totalDataPoints * 0.75f;
    if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
        // XXX: Special case where its in first bin.  Must code last bin.
        binHi25 = 0;
    } else {
        binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar);
    }
    if ((binLo25 < 0) || (binHi25 < 0)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25 and 75 percent data points\n");
        psFree(tmpStatsMinMax);
        psFree(robustHistogram);
        psFree(cumulativeRobustHistogram);
        psFree(tmpScalar);
        psFree(tmpMaskVec);
        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
        return(1);
    }
    psTrace(__func__, 6, "The 25% and 75% data point bins are (%d, %d).\n", binLo25, binHi25);

    //
    // ADD: Step 8.
    // Interpolate to find these two positions exactly: these are the upper
    // and lower quartile positions.
    // XXX: Check for errors.
    //
    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
                           *(psVector* *)&cumulativeRobustHistogram->bounds,
                           *(psVector* *)&cumulativeRobustHistogram->nums,
                           binLo25,
                           totalDataPoints * 0.25f);
    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
                           *(psVector* *)&cumulativeRobustHistogram->bounds,
                           *(psVector* *)&cumulativeRobustHistogram->nums,
                           binHi25,
                           totalDataPoints * 0.75f);
    stats->robustLQ = binLo25F32;
    stats->robustUQ = binHi25F32;
    psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32);

    //
    //
    // New algorithm for calculated mean and stdev.
    //
    //
    psS32 N50 = 0;
    for (psS32 i = 0 ; i < myVector->n ; i++) {
        if ((0 == tmpMaskVec->data.U8[i]) &&
                (binLo25F32 <= myVector->data.F32[i]) &&
                (binHi25F32 >= myVector->data.F32[i])) {
            N50++;
        }
    }
    stats->robustN50 = N50;
    psTrace(__func__, 6, "The robustN50 is %d.\n", N50);

    psF32 dN = (psF32) (0.17 * N50);
    if (dN < 1.0) {
        dN = 1.0;
    } else if (dN > 4.0) {
        dN = 4.0;
    }
    psF32 newBinSize = sigma / dN;

    rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
    rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
    if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
        psFree(tmpStatsMinMax);
        psFree(robustHistogram);
        psFree(cumulativeRobustHistogram);
        psFree(tmpScalar);
        psFree(tmpMaskVec);
        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
        return(1);
    }

    numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize);
    psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
    psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize);
    psTrace(__func__, 6, "The numBins is %d\n", numBins);

    psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
    newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1);
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(newHistogram->nums);
    }

    //
    // Smooth the resulting histogram with a Gaussian with sigma_s = 1
    // bin.
    //
    psF32 sigma_s = newBinSize;

    psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s);
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(newHistogramSmoothed);
    }

    //
    // Find the bin with the peak value in the range 2 sigma of the
    // robust histogram median.
    //

    psS32 binMin = 0;
    psS32 binMax = 0;
    tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
    if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
        binMin = 0;
    } else {
        binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
    }

    tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
    if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
        binMax = newHistogram->bounds->n-1;
    } else {
        binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
    }
    if ((binMin < 0) || (binMax < 0)) {
        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n");
        psFree(tmpMaskVec);
        psFree(robustHistogram);
        psFree(cumulativeRobustHistogram);
        psFree(tmpStatsMinMax);
        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
        return(1);
    }

    psS32 binNum = binMin;
    psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum];
    for (psS32 i = binMin+1 ; i <= binMax ; i++) {
        if (newHistogramSmoothed->data.F32[i] > binMaxNums) {
            binNum = i;
            binMaxNums = newHistogramSmoothed->data.F32[i];
        }
    }
    psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums);

    //
    // Fit a Gaussian to the bins in the range 20 sigma of the robust
    // histogram median.
    //
    tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
    if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
        binMin = 0;
    } else {
        binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
        // XXX: check for errors here.
    }
    tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
    if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
        binMax = newHistogramSmoothed->n - 1;
    } else {
        binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
        // XXX: check for errors here.
    }
    psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
    psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
    psArray *x = psArrayAlloc((1 + (binMax - binMin)));
    stats->robustNfit = 0;
    psS32 j = 0;

    for (psS32 i = binMin ; i <= binMax ; i++) {
        y->data.F32[j] = newHistogramSmoothed->data.F32[i];
        x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
        ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
        xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i);

        stats->robustNfit+= newHistogramSmoothed->data.F32[i];
        j++;
    }
    if (psTraceGetLevel(__func__) >= 8) {
        // XXX: Print the x array somehow.
        PS_VECTOR_PRINT_F32(y);
    }

    // XXX: Use the min/max routines for this
    psF32 minY = FLT_MAX;
    psF32 maxY = -FLT_MAX;
    for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
        if (y->data.F32[i] > maxY) {
            maxY = y->data.F32[i];
        }
        if (y->data.F32[i] < minY) {
            minY = y->data.F32[i];
        }
    }
    //
    // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
    // XXX: Use the normalize routines for this.
    //
    for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
        y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
    }

    //
    psMinimization *min = psMinimizationAlloc(100, 0.01);
    psVector *params = psVectorAlloc(2, PS_TYPE_F32);
    // Initial guess for the mean ([0]) and standard dev.
    params->data.F32[0] = stats->robustMedian;
    params->data.F32[1] = sigma;
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(params);
        PS_VECTOR_PRINT_F32(y);
    }
    psFree(xTmp);
    rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D);

    if (rcBool != true) {
        psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
        psFree(tmpStatsMinMax);
        psFree(robustHistogram);
        psFree(cumulativeRobustHistogram);
        psFree(tmpScalar);
        psFree(newHistogram);
        psFree(x);
        psFree(y);
        psFree(min);
        psFree(params);
        psFree(tmpMaskVec);
        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
        return(1);
    }
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(params);
    }

    //
    // The robust mean mean_r is derived directly from the fitted
    // Gaussian mean.
    //
    stats->robustMean = params->data.F32[0];
    psTrace(__func__, 6, "The robust mean is %f.\n", params->data.F32[0]);

    //
    // The robust standard deviation, SIGMA_r is determined by
    // subtracting the smoothing scale in quadrature:
    //     SIGMA_r^2 = SIGMA^2 - sigma_s^2
    //
    stats->robustStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s));
    psTrace(__func__, 6, "The robust stdev is %f.\n", stats->robustStdev);

    psFree(newHistogramSmoothed);
    psFree(tmpStatsMinMax);
    psFree(cumulativeRobustHistogram);
    psFree(tmpScalar);
    psFree(newHistogram);
    psFree(x);
    psFree(y);
    psFree(min);
    psFree(params);
    psFree(tmpMaskVec);
    psFree(robustHistogram);

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









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

/* FUNCTION IMPLEMENTATION - PUBLIC                                          */

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

static void histogramFree(psHistogram* myHist);

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(__func__, 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->robustMean = NAN;
    newStruct->robustMedian = NAN;
    newStruct->robustMode = NAN;
    newStruct->robustStdev = NAN;
    newStruct->robustUQ = NAN;
    newStruct->robustLQ = NAN;
    newStruct->robustN50 = -1;            // XXX: This is never used
    newStruct->robustNfit = -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;          // XXX: This is never used
    newStruct->options = options;

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


bool psMemCheckStats(psPtr ptr)
{
    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(__func__, 3, "---- %s() begin  ----\n", __func__);
    psTrace(__func__, 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n);
    PS_ASSERT_INT_POSITIVE(n, NULL);
    PS_FLOAT_COMPARE(lower, upper, NULL);

    psS32 i = 0;                  // Loop index variable
    psHistogram* newHist = NULL;        // The new histogram structure
    psF32 binSize = 0.0;        // The histogram bin size

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

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

    // Allocate the bins, and initialize them to zero.
    newHist->nums = psVectorAlloc(n, PS_TYPE_F32);
    for (i = 0; i < newHist->nums->n; i++) {
        newHist->nums->data.F32[i] = 0.0;
    }

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

    psTrace(__func__, 3, "---- %s(psHistogram) 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(__func__, 3, "---- %s() begin  ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(bounds, NULL);
    PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL);
    PS_INT_COMPARE(2, bounds->n, NULL);

    psHistogram* newHist = NULL;        // The new histogram structure
    psS32 i;                      // Loop index variable

    // Allocate memory for the new histogram structure.
    newHist = (psHistogram* ) psAlloc(sizeof(psHistogram));
    psMemSetDeallocator(newHist, (psFreeFunc) histogramFree);
    newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
    *(int *)&newHist->bounds->n = newHist->bounds->nalloc;
    for (i = 0; i < bounds->n; i++) {
        newHist->bounds->data.F32[i] = bounds->data.F32[i];
    }

    // 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);
    for (i = 0; i < newHist->nums->n; i++) {
        newHist->nums->data.F32[i] = 0.0;
    }

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

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

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


bool psMemCheckHistogram(psPtr ptr)
{
    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.  Inputs:
    binNum: the bin number of the data point in the histogram
    out: the histogram structure
    data: the data point value
    error: the error in that data point
 
XXX: Must test this.
 *****************************************************************************/
psS32 UpdateHistogramBins(psS32 binNum,
                          psHistogram* out,
                          psF32 data,
                          psF32 error)
{
    psTrace(__func__, 3, "---- %s() begin  ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(out, -1);
    PS_ASSERT_PTR_NON_NULL(out->bounds, -1);
    PS_ASSERT_PTR_NON_NULL(out->nums, -1);
    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, ((out->nums->n)-1), -2);
    PS_FLOAT_COMPARE(0.0, error, -3);
    PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4);

    psF32 boxcarWidth = 2.35 * error;
    psF32 boxcarCenter = (out->bounds->data.F32[binNum] +
                          out->bounds->data.F32[binNum+1]) / 2.0;
    psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0);
    psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0);
    psS32 bin;
    psS32 boxcarLeftBinNum = 0;
    psS32 boxcarRightBinNum = 0;

    // Determine the left endpoint of the boxcar for the PDF.
    for (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 (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(__func__, 3, "---- %s(0) end  ----\n", __func__);
        return(0);
    }

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

    //
    // Loop through the center bins, if any.
    //
    for (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;

    //
    // Return 0 on success.
    //
    psTrace(__func__, 3, "---- %s(0) end  ----\n", __func__);
    return(0);
}


/*****************************************************************************
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(__func__, 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 != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, NULL);
        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
    }
    if (errors != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NULL);
        PS_ASSERT_VECTOR_TYPE(errors, values->type.type, NULL);
    }

    psS32 i = 0;                  // Loop index variable
    psF32 binSize = 0.0;          // Histogram bin size
    psS32 binNum = 0;             // A temporary bin number
    psS32 numBins = 0;            // The total number of bins
    psScalar tmpScalar;
    tmpScalar.type.type = PS_TYPE_F32;
    psVector* inF32 = NULL;
    psVector* errorsF32 = NULL;
    psS32 mustFreeVectorIn = 1;
    psS32 mustFreeVectorErrors = 1;

    // Convert input and errors vectors to F32 if necessary.
    inF32 = p_psConvertToF32((psVector *) values);
    if (inF32 == NULL) {
        inF32 = (psVector *) values;
        mustFreeVectorIn = 0;
    }
    errorsF32 = p_psConvertToF32((psVector *) errors);
    if (errorsF32 == NULL) {
        errorsF32 = (psVector *) errors;
        mustFreeVectorErrors = 0;
    }

    numBins = out->nums->n;
    for (i = 0; i < inF32->n; i++) {
        // Check if this pixel is masked, and if so, skip it.
        if ((mask == NULL) || ((mask != NULL) && (!(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) {
                    binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0];
                    binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
                    if (errorsF32 != NULL) {
                        // XXX: Check return codes.
                        UpdateHistogramBins(binNum, out,
                                            inF32->data.F32[i],
                                            errorsF32->data.F32[i]);
                    } else {
                        // XXX: 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) {
                            // XXX: Check return codes.
                            UpdateHistogramBins(binNum, out,
                                                inF32->data.F32[i],
                                                errors->data.F32[i]);
                        } else {
                            (out->nums->data.F32[binNum])+= 1.0;
                        }
                    }
                }
            }
        }
    }

    if (mustFreeVectorIn == 1) {
        psFree(inF32);
    }
    if (mustFreeVectorErrors == 1) {
        psFree(errorsF32);
    }

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

/******************************************************************************
p_psConvertToF32(in): this is the cheap way to support a variety of vector
data types: we simply convert the input vector to F32 at the beginning, and
write all of our functions in F32.  If the vast majority of all vector stat
operations are F32 (or any other single type), then this is probably the
best way to go.  Otherwise, when the algorithms stablize, we will then macro
everything and put type support in the various stat functions.
 
XXX: Should the default data type be F64?  Since we are buying Opterons...
 *****************************************************************************/
psVector* p_psConvertToF32(psVector* in)
{
    psTrace(__func__, 4,"---- %s() begin  ----\n", __func__);
    if (in == NULL) {
        return(NULL);
    }
    psS32 i = 0;
    psVector* tmp = NULL;

    if (in->type.type == PS_TYPE_S8) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.S8[i];
        }
    } else if (in->type.type == PS_TYPE_S16) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32) in->data.S16[i];
        }
    } else if (in->type.type == PS_TYPE_S32) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.S32[i];
        }
    } else if (in->type.type == PS_TYPE_S64) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.S64[i];
        }
    } else if (in->type.type == PS_TYPE_U8) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.U8[i];
        }
    } else if (in->type.type == PS_TYPE_U16) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.U16[i];
        }
    } else if (in->type.type == PS_TYPE_U32) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.U32[i];
        }
    } else if (in->type.type == PS_TYPE_U64) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.U64[i];
        }
    } else if (in->type.type == PS_TYPE_F64) {
        tmp = psVectorAlloc(in->n, PS_TYPE_F32);
        for (i = 0; i < in->n; i++) {
            tmp->data.F32[i] = (psF32)in->data.F64[i];
        }
    } else if (in->type.type == PS_TYPE_F32) {
        // do nothing
    } else {
        char* strType;
        PS_TYPE_NAME(strType, in->type.type);
        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                PS_ERRORTEXT_psStats_VECTOR_TYPE_UNSUPPORTED,
                strType);
    }

    psTrace(__func__, 4,"---- %s(psVector) end  ----\n", __func__);
    return (tmp);
}

/******************************************************************************
psVectorStats(myVector, maskVector, 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
    myVector
    maskVector
    maskVal
    stats
Returns
    The stats structure.
 *****************************************************************************/
psStats* psVectorStats(psStats* stats,
                       const psVector* in,
                       const psVector* errors,
                       const psVector* mask,
                       psMaskType maskVal)
{
    psTrace(__func__, 3,"---- %s() begin  ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(stats, NULL);
    PS_ASSERT_VECTOR_NON_NULL(in, stats);
    if (mask != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, stats);
        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, stats);
    }
    if (errors != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(errors, in, stats);
        PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats);
    }

    psVector* inF32 = NULL;
    psVector* errorsF32 = NULL;
    psS32 mustFreeVectorIn = 1;
    psS32 mustFreeVectorErrors = 1;

    inF32 = p_psConvertToF32((psVector *) in);
    if (inF32 == NULL) {
        inF32 = (psVector *) in;
        mustFreeVectorIn = 0;
    }
    errorsF32 = p_psConvertToF32((psVector *) errors);
    if (errorsF32 == NULL) {
        errorsF32 = (psVector *) errors;
        mustFreeVectorErrors = 0;
    }

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

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEAN) {
        if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n");
        }
    }
    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
        if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: psVectorStats(): p_psVectorSampleMedian() returned an error.\n");
        }
    }
    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_STDEV) {
        if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n");
        }
        p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats);
    }
    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
        if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: psVectorStats(): p_psVectorSampleQuartiles() returned an error.\n");
        }
    }
    // Since the various robust stats quantities share much computation, they
    // are grouped together in a single private function:
    // p_psVectorRobustStats()
    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
            (stats->options & PS_STAT_ROBUST_MODE) ||
            (stats->options & PS_STAT_ROBUST_STDEV) ||
            (stats->options & PS_STAT_ROBUST_QUARTILE)) {
        if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false,
                    PS_ERRORTEXT_psStats_STATS_FAILED);
            // XXX: Set to NAN
            // XXX: Is this the right thing to do?
            // XXX: If so, do it for other funcs?
            psFree(stats);
            psTrace(__func__, 3,"---- %s(NULL) end  ----\n", __func__);
            return(NULL);
        }
    }

    // XXX: Different conditions for return -1 and -2?
    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
        psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats);
        if (-1 == rc) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to calculate clipped statistics for input psVector.\n");
            stats->clippedMean = NAN;
            stats->clippedStdev = NAN;
        } else if (-2 == rc) {
            psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector.");
            stats->clippedMean = NAN;
            stats->clippedStdev = NAN;
        }
    }
    // ************************************************************************
    if (stats->options & PS_STAT_MAX) {
        if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to calculate vector maximum");
            stats->max = NAN;
        }
    }
    // ************************************************************************
    if (stats->options & PS_STAT_MIN) {
        if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to calculate vector minimum");
            stats->min = NAN;
        }
    }

    if (mustFreeVectorIn == 1) {
        psFree(inF32);
    }
    if (mustFreeVectorErrors == 1) {
        psFree(errorsF32);
    }
    psTrace(__func__, 3,"---- %s(psStats) end  ----\n", __func__);
    return (stats);
}

