/** @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.133 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-06-13 20:18:18 $
 *
 *  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 "psFunctions.h"
#include "psConstants.h"

#include "psDataManipErrors.h"

/*****************************************************************************/
/* DEFINE STATEMENTS                                                         */
/*****************************************************************************/
#define PS_GAUSS_WIDTH 5       // The width of the Gaussian or boxcar smoothing.
#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)
{
    switch (stats->options & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE | PS_STAT_ROBUST_FOR_SAMPLE)) {
    case PS_STAT_SAMPLE_MEAN:
        *value = stats->sampleMean;
        return true;

    case PS_STAT_SAMPLE_MEDIAN:
        *value = stats->sampleMedian;
        return true;

    case PS_STAT_SAMPLE_STDEV:
        *value = stats->sampleStdev;
        return true;

    case PS_STAT_ROBUST_MEAN:
        *value = stats->robustMean;
        return true;

    case PS_STAT_ROBUST_MEDIAN:
        *value = stats->robustMedian;
        return true;

    case PS_STAT_ROBUST_MODE:
        *value = stats->robustMode;
        return true;

    case PS_STAT_ROBUST_STDEV:
        *value = stats->robustStdev;
        return true;

    case PS_STAT_CLIPPED_MEAN:
        *value = stats->clippedMean;
        return true;

    case PS_STAT_CLIPPED_STDEV:
        *value = stats->clippedStdev;
        return true;

    case PS_STAT_MAX:
        *value = stats->max;
        return true;

    case PS_STAT_MIN:
        *value = stats->min;
        return true;

    default:
        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)
{

    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)) {
        return(-1);
    } else {
        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.
 *****************************************************************************/
psS32 p_psVectorMax(const psVector* myVector,
                    const psVector* maskVector,
                    psU32 maskVal,
                    psStats* stats)
{
    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;
        return(1);
    }
    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)
{
    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;
        return(1);
    }
    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)
{
    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)) {
                    return(true);
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((stats->min <= myVector->data.F32[i]) &&
                        (myVector->data.F32[i] <= stats->max)) {
                    return(true);
                }
            }
        }
    } else {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    return(true);
                }
            }
        } else {
            if (myVector->n > 0) {
                return(true);
            }
        }
    }
    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)
{
    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;
        }
    }

    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)
{
    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);
        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.
    return(true);
}

/******************************************************************************
p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
robustHistogram with a Gaussian of width sigma.
 
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* robustHistogram,
                                       psF32 sigma)
{
    PS_ASSERT_PTR_NON_NULL(robustHistogram, NULL);
    PS_ASSERT_PTR_NON_NULL(robustHistogram->bounds, NULL);

    psS32 i = 0;                  // Loop index variable
    psS32 j = 0;                  // Loop index variable
    psF32 iMid;
    psF32 jMid;
    psS32 numBins = robustHistogram->nums->n;
    psS32 numBounds = robustHistogram->bounds->n;
    psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32);
    psS32 jMin = 0;
    psS32 jMax = 0;
    psF32 firstBound = robustHistogram->bounds->data.F32[0];
    psF32 lastBound = robustHistogram->bounds->data.F32[numBounds-1];
    psScalar x;

    x.type.type = PS_TYPE_F32;
    for (i = 0; i < numBins; i++) {
        // Determine the midpoint of bin i.
        iMid = (robustHistogram->bounds->data.F32[i] +
                robustHistogram->bounds->data.F32[i+1]) / 2.0;


        // We determine the bin numbers corresponding to a range of data
        // values surrounding iMid.  The ranges is of size
        // s*PS_GAUSS_WIDTH*sigma

        // YYY: The p_psVectorBinDisect() routine does much of the work of
        // the following conditionals, however, it also reports a warning
        // message.  I don't want the warning message so I reproduce the
        // conditionals here.  Maybe p_psVectorBinDisect() should not produce
        // warnings?

        x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma);
        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
            jMin = p_psVectorBinDisect(robustHistogram->bounds, &x);
            if (jMin < 0) {
                psError(PS_ERR_UNEXPECTED_NULL,
                        false,
                        PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
                return(NULL);
            }
        } else if (x.data.F32 <= firstBound) {
            jMin = 0;
        } else if (x.data.F32 >= lastBound) {
            jMin = robustHistogram->bounds->n - 1;
        }

        x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma);
        if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) {
            jMax = p_psVectorBinDisect(robustHistogram->bounds, &x);
            if (jMax < 0) {
                psError(PS_ERR_UNEXPECTED_NULL,
                        false,
                        PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM);
                return(NULL);
            }
        } else if (x.data.F32 <= firstBound) {
            jMax = 0;
        } else if (x.data.F32 >= lastBound) {
            jMax = robustHistogram->bounds->n - 1;
        }

        smooth->data.F32[i] = 0.0;
        for (j = jMin ; j <= jMax ; j++) {
            jMid = (robustHistogram->bounds->data.F32[j] +
                    robustHistogram->bounds->data.F32[j+1]) / 2.0;
            smooth->data.F32[i] +=
                robustHistogram->nums->data.F32[j] *
                psGaussian(jMid, iMid, sigma, true);
        }
    }

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

        }
    }
}

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

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

/******************************************************************************
p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
polynomials here) 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) == midpoint.  This functions 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)
{
    PS_ASSERT_POLY_NON_NULL(myPoly, NAN);
    PS_FLOAT_COMPARE(rangeLow, rangeHigh, 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,
                PS_ERRORTEXT_psStats_STATS_POLY_MEDIAN_OUT_OF_RANGE);
        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++;
    }
    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.
*****************************************************************************/
psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec,
                                        psVector *yVec,
                                        psS32 binNum,
                                        psF32 yVal)
{
    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);

    //    PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64);
    //    PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64);
    //    PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64);
    //    PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD);
    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];

        // 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]);
        }
        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);
        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);
            return(NAN);
        }
        // Call p_ps1DPolyMedian(), which does a binary search on the
        // polynomial, looking for the value x such that f(x) = 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);
            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]);
        }
    }

    psFree(myPoly);
    psFree(x);
    psFree(y);
    psFree(yErr);
    return(tmpFloat);
}

/******************************************************************************
p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
calculates a variety of robust stat measures:
    PS_STAT_ROBUST_MEAN
    PS_STAT_ROBUST_MEDIAN
    PS_STAT_ROBUST_MODE
    PS_STAT_ROBUST_STDEV
    PS_STAT_ROBUST_QUARTILE
I have included all that computation in a single function, as opposed to
breaking it across several functions for one primary reason:  they all require
the same basic initial processing steps (calculate the histogram, etc.).
 
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    0 on success.
 
XXX: Check for errors in psLib routines that we call.
*****************************************************************************/
psS32 p_psVectorRobustStats(const psVector* myVector,
                            const psVector* errors,
                            const psVector* maskVector,
                            psU32 maskVal,
                            psStats* stats)
{
    psHistogram* robustHistogram = NULL;
    psVector* robustHistogramVector = NULL;
    psF32 binSize = 0.0;        // Size of the histogram bins
    psS32 LQBinNum = -1;          // Bin num for lower quartile
    psS32 UQBinNum = -1;          // Bin num for upper quartile
    psS32 medianBinNum = -1;
    psS32 i = 0;                  // Loop index variable
    psS32 modeBinNum = 0;
    psF32 modeBinCount = 0.0;
    psF32 dL = 0.0;
    psS32 numBins = 0;
    psF32 myMean = 0.0;
    psF32 myStdev = 0.0;
    psF32 countFloat = 0.0;
    psF32 diff = 0.0;
    psF32 sumSquares = 0.0;
    psF32 sumDiffs = 0.0;
    psVector* cumulativeRobustSums = NULL;
    psF32 sumRobust = 0.0;
    psF32 sumN50 = 0.0;
    psF32 sumNfit = 0.0;
    psScalar tmpScalar;
    tmpScalar.type.type = PS_TYPE_F32;
    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);

    // Compute the initial bin size of the robust histogram.  This is done
    // by computing the clipped standard deviation of the vector, and dividing
    // that by 10.0;
    //XXX: add errors
    psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats);
    if (rc != 0) {
        psError(PS_ERR_UNEXPECTED_NULL,
                false,
                PS_ERRORTEXT_psStats_ROBUST_STATS_CLIPPED_STATS);
        return(1);
    }
    binSize = tmpStats->clippedStdev / 10.0f;

    // If stats->clippedStdev == 0.0, then all data elements have the same
    // value.  Therefore, we can set the appropiate results and return.
    if (fabs(binSize) <= FLT_EPSILON) {
        if (stats->options & PS_STAT_ROBUST_MEAN) {
            stats->robustMean = tmpStats->clippedMean;
        }
        if (stats->options & PS_STAT_ROBUST_MEDIAN) {
            stats->robustMedian = tmpStats->clippedMean;
        }
        if (stats->options & PS_STAT_ROBUST_MODE) {
            stats->robustMode = tmpStats->clippedMean;
        }
        if (stats->options & PS_STAT_ROBUST_STDEV) {
            stats->robustStdev = 0.0;
        }
        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
            stats->robustUQ = tmpStats->clippedMean;
            stats->robustLQ = tmpStats->clippedMean;
        }
        // XXX: Set these to the number of unmasked data points?
        stats->robustNfit = 0.0;
        stats->robustN50 = 0.0;
        psFree(tmpStats);
        return(0);
    }

    // Determine minimum and maximum values in the data vector.
    // XXX: remove this conditional?
    if (isnan(tmpStats->min)) {
        if (0 != p_psVectorMin(myVector, maskVector, maskVal, tmpStats)) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: p_psVectorMin(): p_psVectorMin() reported a NAN mean.\n");
            return(1);
        }
    }
    // XXX: remove this conditional?
    if (isnan(tmpStats->max)) {
        if (0 != p_psVectorMax(myVector, maskVector, maskVal, tmpStats)) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: p_psVectorMin(): p_psVectorMax() reported a NAN mean.\n");
            return(1);
        }
    }

    // Create the histogram structure.  NOTE: we can not specify the bin size
    // precisely since the argument to psHistogramAlloc() is the number of
    // bins, not the binSize.  Also, if we get here, we know that
    // binSize != 0.0.
    numBins = (psS32)((tmpStats->max - tmpStats->min) / binSize);
    robustHistogram = psHistogramAlloc(tmpStats->min, tmpStats->max, numBins);

    // Populate the histogram array.
    psVectorHistogram(robustHistogram, myVector, errors, maskVector, maskVal);

    // Smooth the histogram, Gaussian-style.
    robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram,
                            tmpStats->clippedStdev / 4.0f);

    /**************************************************************************
    Determine the median/lower/upper quartile bin numbers.

    We define a vector called "cumulativeRobustSums" where the value at
    index position i is equal to the sum of bins 0:i.  This will be used in
    determining the median and lower/upper quartiles.
    **************************************************************************/
    cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0];
    for (i = 1; i < robustHistogramVector->n; i++) {
        cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] +
                                            robustHistogramVector->data.F32[i];
    }
    sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];

    tmpScalar.data.F32 = sumRobust / 4.0;
    LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
    tmpScalar.data.F32 = 3.0 * sumRobust / 4.0;
    UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);
    tmpScalar.data.F32 = sumRobust / 2.0;
    medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar);

    if ((LQBinNum < 0) || (UQBinNum < 0)) {
        psError(PS_ERR_UNKNOWN, true,
                PS_ERRORTEXT_psStats_ROBUST_QUARTILE_BINS_FAILED);
        return(1);
    }
    if (medianBinNum < 0) {
        psError(PS_ERR_UNKNOWN, true,
                PS_ERRORTEXT_psStats_ROBUST_QUARTILE_BINS_FAILED);
        return(1);
    }
    /**************************************************************************
    Determine the mode in the range LQ:UQ.
    **************************************************************************/
    // Determine the bin with the peak value in the range LQ to UQ.
    modeBinNum = LQBinNum;
    modeBinCount = robustHistogramVector->data.F32[LQBinNum];
    sumN50 = robustHistogram->nums->data.F32[LQBinNum];
    for (i = LQBinNum + 1; i <= UQBinNum; i++) {
        if (robustHistogramVector->data.F32[i] > modeBinCount) {
            modeBinNum = i;
            modeBinCount = robustHistogramVector->data.F32[i];
        }
        sumN50 += robustHistogram->nums->data.F32[i];
    }

    dL = (UQBinNum - LQBinNum) / 4;
    /**************************************************************************
    Determine the mean/stdev for the bins in the range mode-dL to mode+dL
    **************************************************************************/
    // Calculate the mean of the smoothed robust histogram in the range
    // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
    // that bin (this is a non-exact approximation).
    sumNfit = 0.0;
    myMean = 0.0;
    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
        if ((0 <= i) && (i < robustHistogramVector->n)) {
            myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i);
            countFloat += robustHistogramVector->data.F32[i];
        }

        sumNfit += robustHistogram->nums->data.F32[i];
    }
    // XXX: divide by zero?
    myMean /= countFloat;

    // Calculate the stdev of the smoothed robust histogram in the range
    // mode-dL to mode+dL.  We use the midpoint of each bin as the mean for
    // that bin.
    for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
        if ((0 <= i) && (i < robustHistogramVector->n)) {
            diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean;
            sumSquares += diff * diff * robustHistogramVector->data.F32[i];
            sumDiffs += diff * robustHistogramVector->data.F32[i];
        }
    }
    myStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));

    p_psNormalizeVectorRangeF32(robustHistogramVector, 0.0, 1.0);

    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
            (stats->options & PS_STAT_ROBUST_STDEV)) {
        // We fit a 1-D polynomial to the data.
        // XXX: I implement the changes requested in bug 366.  However, this
        // code no longer produces sensible results.
        // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably
        // remove some of theabove code that calculated the initial estimate
        // for the mean and sigma.

        psVector *y = psVectorAlloc(2 * dL + 1, PS_TYPE_F32);
        psVector *x = psVectorAlloc(2 * dL + 1, PS_TYPE_F32);
        for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
            int index = i - modeBinNum + dL;
            // XXX: Should this be the natural log?
            y->data.F32[index] = robustHistogramVector->data.F32[i];
            //            y->data.F32[index] = logf(robustHistogramVector->data.F32[i]);
            x->data.F32[index] = (psF32) index;
        }

        psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD);
        // XXX: What about the NULL x argument?
        tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, y, NULL);
        if (tmpPoly == NULL) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: failed fit a 1D polynomial.\n");
        }
        psF32 polyFitSigma = sqrtf(-0.5 / tmpPoly->coeff[2]);
        psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma);
        // psF32 polyFitNorm = exp(tmpPoly->coedd[0] + PS_SQR(polyFitMean) / (2.0 * PS_SQR(polyFitSigma)));

        if (stats->options & PS_STAT_ROBUST_MEAN) {
            stats->robustMean = polyFitMean;
        }

        if (stats->options & PS_STAT_ROBUST_STDEV) {
            stats->robustStdev = polyFitSigma;
        }
        // For testing only.  This shows that the polynomial fit is successful.
        if (0) {
            printf("--- Results of the 1-D polynomial fit:\n");
            for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
                int index = i - modeBinNum + dL;
                printf("(x, y, poly(x)) is (%.2f, %.2f, %.2f)\n",
                       x->data.F32[index],
                       y->data.F32[index],
                       psPolynomial1DEval(tmpPoly, x->data.F32[index]));
            }
        }

        psFree(x);
        psFree(y);
        psFree(tmpPoly);
    }


    /**************************************************************************
    Set the appropriate members in the output stats struct.
    **************************************************************************/

    if (stats->options & PS_STAT_ROBUST_MODE) {
        stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum);
    }


    // To determine the median, we fit a quadratic y=f(x) to the three bins
    // surrounding the bin containing the median (x is the midpoint of each
    // bin and y is the value of each bin).  Then we figure out what value
    // of x corresponds to f(x) being the median (half of all points).
    if (stats->options & PS_STAT_ROBUST_MEDIAN) {
        // Take a psVector.  Fit a polynomial y = f(x) to the 3 data elements
        // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0.

        //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths.
        //Determine id that is okay.
        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
                                  robustHistogram->bounds,
                                  cumulativeRobustSums,
                                  medianBinNum,
                                  sumRobust/2.0);
    }

    // To determine the quartiles, we fit a quadratic y=f(x) to the three bins
    // surrounding the bin containing LQ/UQ (x is the midpoint of each
    // bin and y is the value of each bin).  Then we figure out what value
    // of x corresponds to f(x) being the LQ/UQ.
    if (stats->options & PS_STAT_ROBUST_QUARTILE) {
        countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1];

        //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths.
        //Determine id that is okay.
        stats->robustLQ = fitQuadraticSearchForYThenReturnX(
                              robustHistogram->bounds,
                              cumulativeRobustSums,
                              LQBinNum,
                              countFloat/4.0);
        //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths.
        //Determine id that is okay.
        stats->robustUQ = fitQuadraticSearchForYThenReturnX(
                              robustHistogram->bounds,
                              cumulativeRobustSums,
                              UQBinNum,
                              3.0 * countFloat/4.0);
    }
    // XXX: I think sumNfit == sumN50 here.
    stats->robustNfit = sumNfit;
    stats->robustN50 = sumN50;

    psFree(tmpStats);
    psFree(robustHistogram);
    psFree(robustHistogramVector);
    psFree(cumulativeRobustSums);
    return(0);
}

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

/* FUNCTION IMPLEMENTATION - PUBLIC                                          */

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

static void histogramFree(psHistogram* myHist);

/******************************************************************************
    psStatsAlloc(): This routine must create a new psStats data structure.
 *****************************************************************************/
psStats* psStatsAlloc(psStatsOptions options)
{
    psStats* newStruct = NULL;

    newStruct = (psStats* ) psAlloc(sizeof(psStats));
    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;

    return (newStruct);
}

/******************************************************************************
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(psF32 lower, psF32 upper, psS32 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, (psFreeFcn) histogramFree);
    newHist->bounds = psVectorAlloc(n + 1, PS_TYPE_F32);
    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;

    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)
{
    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, (psFreeFcn) histogramFree);
    newHist->bounds = psVectorAlloc(bounds->n, PS_TYPE_F32);
    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;

    return (newHist);
}

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

/*****************************************************************************
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)
{
    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;
        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.
    //
    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* in,
                               const psVector* errors,
                               const psVector* mask,
                               psU32 maskVal)
{
    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(in, out);
    if (mask != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(in, mask, NULL);
        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
    }
    if (errors != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(in, errors, NULL);
        PS_ASSERT_VECTOR_TYPE(errors, in->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 *) in);
    if (inF32 == NULL) {
        inF32 = (psVector *) in;
        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(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);
    }
    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)
{
    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);
    }
    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,
                       psU32 maskVal)
{
    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: 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);
    }
    return (stats);
}

