/** @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
 *
 *  @version $Revision: 1.67 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-10-13 02:40:13 $
 *
 *  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 "psError.h"
#include "psAbort.h"
#include "psStats.h"
#include "psMinimize.h"
#include "psFunctions.h"
#include "psConstants.h"

/*****************************************************************************/
/* DEFINE STATEMENTS                                                         */
/*****************************************************************************/
// will use robust statistical methods.
#define GAUSS_WIDTH 20                     // The width of the Gaussian or boxcar smoothing.
#define CLIPPED_NUM_ITER_LB 1
#define CLIPPED_NUM_ITER_UB 10
#define CLIPPED_SIGMA_LB 1.0
#define CLIPPED_SIGMA_UB 10.0
#define true 1
#define false 0
#define MY_MAX_FLOAT HUGE
#define MAX_ITERATIONS 10

void p_psVectorRobustStats(const psVector* restrict myVector,
                           const psVector* restrict maskVector, unsigned int maskVal, psStats* stats);


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

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

// None

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

// None

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

bool p_psGetStatValue(const psStats* stats, double *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;
    }
}

/*****************************************************************************
p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal
function that simply prints a vector to STDOUT.  Used primarily for
debugging.
 *****************************************************************************/

void p_psVectorPrint(psVector* myVector, psVector* maskVector, unsigned int maskVal, psStats* stats)
{
    int i = 0;                  // Loop index variable.

    for (i = 0; i < myVector->n; i++) {
        if (maskVector != NULL)
            printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]);
        else
            printf("Element %d is %f\n", i, myVector->data.F32[i]);
    }
}

/******************************************************************************
 ******************************************************************************
 ******************************************************************************
    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().
        Is the mask vector the same length as the data vector?
        Is the mask vector of type PS_TYPE_U8?
        Is the stats data structure NULL?
        Is the in data structure NULL?
        Is the in data structure of type PS_TYPE_F32?
 ******************************************************************************
 ******************************************************************************
 *****************************************************************************/

/******************************************************************************
p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
mean of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
ASSUMPTION: the mean is always calculated exactly.  Robust means are never
calculated in this routine.
 *****************************************************************************/

void p_psVectorSampleMean(const psVector* restrict myVector,
                          const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    int i = 0;                  // Loop index variable
    float mean = 0.0;           // The mean
    int count = 0;              // # of points in this mean?
    float rangeMin = 0.0;       // Exclude data below this
    float rangeMax = 0.0;       // Exclude date above this

    // If PS_STAT_USE_RANGE is requested, then we enter a slightly different
    // loop.
    if (stats->options & PS_STAT_USE_RANGE) {
        rangeMin = stats->min;
        rangeMax = stats->max;
        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]) &&
                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    mean += myVector->data.F32[i];
                    count++;
                }
            }
            mean /= (float)count;
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    mean += myVector->data.F32[i];
                    count++;
                }
            }
            mean /= (float)count;
        }
    } else {
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    mean += myVector->data.F32[i];
                    count++;
                }
            }
            mean /= (float)count;
        } else {
            for (i = 0; i < myVector->n; i++) {
                mean += myVector->data.F32[i];
            }
            mean /= (float)myVector->n;
        }
    }

    stats->sampleMean = mean;
}

/******************************************************************************
p_psVectorSampleMax(myVector, maskVector, maskVal, stats): calculates the
max of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 *****************************************************************************/
void p_psVectorMax(const psVector* restrict myVector,
                   const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    int i = 0;                  // Loop index variable
    float max = -MY_MAX_FLOAT;  // The calculated maximum
    float rangeMin = 0.0;       // Exclude data below this
    float rangeMax = 0.0;       // Exclude date above this

    // If PS_STAT_USE_RANGE is requested, then we enter a different loop.
    if (stats->options & PS_STAT_USE_RANGE) {
        rangeMin = stats->min;
        rangeMax = stats->max;
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i])) {
                    if ((myVector->data.F32[i] > max) &&
                            (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                        max = myVector->data.F32[i];
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((myVector->data.F32[i] > max) &&
                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    max = myVector->data.F32[i];
                }
            }
        }
    } 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];
                    }
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if (myVector->data.F32[i] > max) {
                    max = myVector->data.F32[i];
                }
            }
        }
    }

    stats->max = max;
}

/******************************************************************************
p_psVectorSampleMin(myVector, maskVector, maskVal, stats): calculates the
minimum of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 *****************************************************************************/
void p_psVectorMin(const psVector* restrict myVector,
                   const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    int i = 0;                  // Loop index variable
    float min = MY_MAX_FLOAT;   // The calculated maximum
    float rangeMin = 0.0;       // Exclude data below this
    float rangeMax = 0.0;       // Exclude date above this

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

    stats->min = min;
}

/******************************************************************************
p_psVectorNValues(myVector, maskVector, maskVal, stats): calculates the
number of non-masked pixels in the vector that fall within the min/max
range, if given.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 *****************************************************************************/
int p_psVectorNValues(const psVector* restrict myVector,
                      const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    int i = 0;                  // Loop index variable
    int numData = 0;            // The number of data points
    float rangeMin = 0.0;       // Exclude data below this
    float rangeMax = 0.0;       // Exclude date above this

    if (stats->options & PS_STAT_USE_RANGE) {
        rangeMin = stats->min;
        rangeMax = stats->max;
        if (maskVector != NULL) {
            for (i = 0; i < myVector->n; i++) {
                if (!(maskVal & maskVector->data.U8[i]) &&
                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    numData++;
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    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.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 *****************************************************************************/
void p_psVectorSampleMedian(const psVector* restrict myVector,
                            const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    psVector* unsortedVector = NULL;    // Temporary vector
    psVector* sortedVector = NULL;      // Temporary vector
    int i = 0;                  // Loop index variable
    int count = 0;              // # of points in this mean?
    int nValues = 0;            // # of points in vector
    float rangeMin = 0.0;       // Exclude data below this
    float rangeMax = 0.0;       // Exclude date above this

    // Determine if the number of data points exceed a threshold which will
    // cause to generate robust stats, as opposed to exact stats.
    // XXX: This code is no longer used.  We now calculate the median exactly
    // regardless of the vector size.
    /*
     * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     * 
     * // Store the robust quartiles into the sample quartile members. stats->sampleMedian =
     * stats2->robustMedian;
     * 
     * // Free temporary data buffers. psFree(stats2);
     * 
     * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options |
     * PS_STAT_ROBUST_FOR_SAMPLE;
     * 
     * return; } */

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

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

    sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    sortedVector->n = sortedVector->nalloc;

    // Determine if we must only use data points within a min/max range.
    if (stats->options & PS_STAT_USE_RANGE) {
        rangeMin = stats->min;
        rangeMax = stats->max;

        // 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]) &&
                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    unsortedVector->data.F32[count++] = maskVector->data.F32[i];
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    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.
    psVectorSort(sortedVector, unsortedVector);

    // Calculate the median exactly.
    // XXX: Is this the correct action?
    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);
}

/******************************************************************************
    This routine smoothes the data in the input robustHistogram with a
    Gaussian of width sigma.
 
    XXX: Must consult with IfA on the proper values for the Gaussian
    smoothing.  Currently, the Gaussian coefficients are such that only
    the center coefficient is non-zero.  This is because "e" is being
    raised to a very large negative number for all points except the
    center point.
 
    XXX: use a static variable for gaussianCoefs[] and compute them once.
 *****************************************************************************/
psVector* p_psVectorsmoothHistGaussian(psHistogram* robustHistogram,
                                       float sigma)
{
    int i = 0;                  // Loop index variable
    int j = 0;                  // Loop index variable
    float denom = 0.0;          // Temporary variable
    float expo = 0.0;           // Temporary variable
    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients
    psVector* smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32);

    for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) {
        if (fabs(sigma) >= FLT_EPSILON) {
            // If sigma does not equal zero, then we use Gaussian smoothing.
            #ifdef  DARWIN
            denom = (float)sqrt(2.0 * M_PI * sigma * sigma);
            #else

            denom = sqrtf(2.0 * M_PI * sigma * sigma);
            #endif

            expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH));
            expo /= (2.0 * sigma * sigma);
            gaussianCoefs[i] = exp(expo / denom);

            //printf("p_psVectorsmoothHistGaussian(): gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
            // NOTE: Gaussian smoothing just isn't working with low sigma
            // values.  The problem is that the Gaussian coefficients are
            // all zero, except for the middle coefficient, which is exactly
            // one.  Therefore, I'm using boxcar smoothing.
            gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH));
            // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]);
        } else {
            /* If sigma equals zero (all pixels have the same value) the above code will divide by zero.
             * Therefore, we don't need to smooth the data. */
            for (i = 0; i < robustHistogram->nums->n; i++) {
                smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i];
            }
            return (smooth);
        }
    }

    // Perform the actual smoothing.
    for (i = 0; i < robustHistogram->nums->n; i++) {
        smooth->data.F32[i] = 0.0;
        for (j = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) {
            if (((j + i) >= 0) && ((j + i) < smooth->n)) {
                smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] *
                                        (float)robustHistogram->nums->data.S32[j + i]);
            }
        }
    }
    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
 *****************************************************************************/
void p_psVectorSampleQuartiles(const psVector* restrict myVector,
                               const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    psVector* unsortedVector = NULL;    // Temporary vector
    psVector* sortedVector = NULL;      // Temporary vector
    int i = 0;                  // Loop index variable
    int count = 0;              // # of points in this mean?
    int nValues = 0;            // # data points
    float rangeMin = 0.0;       // Exclude data below this
    float rangeMax = 0.0;       // Exclude date above this

    // 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);
    unsortedVector->n = unsortedVector->nalloc;
    sortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
    sortedVector->n = sortedVector->nalloc;

    // Determine if we must only use data points within a min/max range.
    if (stats->options & PS_STAT_USE_RANGE) {
        rangeMin = stats->min;
        rangeMax = stats->max;
        // 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]) &&
                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    unsortedVector->data.F32[count++] = myVector->data.F32[i];
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    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.
    psVectorSort(sortedVector, unsortedVector);

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

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

/******************************************************************************
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* restrict myVector,
                           const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    int i = 0;                  // Loop index variable
    int countInt = 0;           // # of data points being used
    float countFloat = 0.0;     // # of data points being used
    float mean = 0.0;           // The mean
    float diff = 0.0;           // Used in calculating stdev
    float sumSquares = 0.0;     // temporary variable
    float sumDiffs = 0.0;       // temporary variable
    float rangeMin = 0.0;       // Exclude data below this
    float rangeMax = 0.0;       // Exclude date above this

    // This procedure requires the mean.  If it has not been already
    // calculated, then call p_psVectorSampleMean()
    if (0 != isnan(stats->sampleMean)) {
        p_psVectorSampleMean(myVector, maskVector, maskVal, stats);
    }
    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]) &&
                        (rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    diff = myVector->data.F32[i] - mean;
                    sumSquares += (diff * diff);
                    sumDiffs += diff;
                    countInt++;
                }
            }
        } else {
            for (i = 0; i < myVector->n; i++) {
                if ((rangeMin <= myVector->data.F32[i]) && (myVector->data.F32[i] <= rangeMax)) {
                    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;
        }
    }
    countFloat = (float)countInt;

    #ifdef DARWIN

    stats->sampleStdev = (float)sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    #else

    stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    #endif
}

/******************************************************************************
p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the
clipped stats (mean or stdev) of the input vector.
 
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 *****************************************************************************/
void p_psVectorClippedStats(const psVector* restrict myVector,
                            const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    int i = 0;                  // Loop index variable
    int j = 0;                  // Loop index variable
    float clippedMean = 0.0;    // self-explanatory
    float clippedStdev = 0.0;   // self-explanatory
    float oldStanMean = 0.0;    // Temporary variable
    float oldStanStdev = 0.0;   // Temporary variable
    psVector* tmpMask = NULL;   // Temporary vector

    // Endure that stats->clipIter is within the proper range.
    if (!((CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= CLIPPED_NUM_ITER_UB))) {
        psAbort(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter);
    }
    // Endure that stats->clipSigma is within the proper range.
    if (!((CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= CLIPPED_SIGMA_UB))) {
        psAbort(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma);
    }
    // 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);
    tmpMask->n = myVector->n;

    // If we were called with a mask vector, then initialize the temporary
    // mask vector with those values.
    if (maskVector != NULL) {
        for (i = 0; i < tmpMask->n; i++) {
            tmpMask->data.U8[i] = maskVector->data.U8[i];
        }
    }
    // 1. Compute the sample median.
    // NOTE: This seems odd.  Verify with IfA that we want to calculate the
    // median here, not the mean.
    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);

    // 2. Compute the sample standard deviation.
    p_psVectorSampleStdev(myVector, maskVector, maskVal, stats);

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

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

    // Must save the old sampleMean and sampleStdev since the following code
    // block overwrites them.
    oldStanMean = stats->sampleMean;
    oldStanStdev = stats->sampleStdev;

    // 5. Repeat N times:
    for (i = 0; i < stats->clipIter; i++) {
        for (j = 0; j < myVector->n; j++) {
            // a) Exclude all values x_i for which |x_i - x| > K * stdev
            if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
                tmpMask->data.U8[i] = 0xff;
            }
            // b) compute new mean and stdev
            p_psVectorSampleMedian(myVector, tmpMask, maskVal, stats);
            p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats);

            // c) Use the new mean for x
            clippedMean = stats->sampleMean;

            // d) Use the new stdev for stdev
            clippedStdev = stats->sampleStdev;
        }

    }
    stats->sampleMean = oldStanMean;
    stats->sampleStdev = oldStanStdev;

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

/*****************************************************************************
p_psNormalizeVectorF32(myData): this is a private function which normalizes the
elements of a vector to a range between 0.0 and 1.0.
 *****************************************************************************/
void p_psNormalizeVectorF32(psVector* myData)
{
    float min = (float)HUGE;
    float max = (float)-HUGE;
    int i = 0;

    for (i = 0; i < myData->n; i++) {
        if (myData->data.F32[i] < min) {
            min = myData->data.F32[i];
        }
        if (myData->data.F32[i] > max) {
            max = myData->data.F32[i];
        }
    }

    //  float range = max - min;
    //    for (i = 0; i < myData->n; i++) {
    //        myData->data.F32[i] = (myData->data.F32[i] - min) / range;
    //    }
    for (i = 0; i < myData->n; i++) {
        myData->data.F32[i] = (myData->data.F32[i] - 0.5 * (min + max)) /
                              (0.5 * (max - min));
    }
}

/*****************************************************************************
p_psNormalizeVectorF64(myData): this is a private function which normalizes the
elements of a vector to a range between -1.0 and 1.0.
XXX: 0-1 or -1:1?
 *****************************************************************************/
void p_psNormalizeVectorF64(psVector* myData)
{
    float min = (double)HUGE;
    float max = (double)-HUGE;
    double range = 0.0;
    int i = 0;

    for (i = 0; i < myData->n; i++) {
        if (myData->data.F64[i] < min) {
            min = myData->data.F64[i];
        }
        if (myData->data.F64[i] > max) {
            max = myData->data.F64[i];
        }
    }

    range = max - min;
    for (i = 0; i < myData->n; i++) {
        myData->data.F64[i] = -1.0 + 2.0 *
                              ((myData->data.F64[i] - min) / range);
    }

    //    for (i = 0; i < myData->n; i++) {
    //        myData->data.F64[i] = (myData->data.F64[i] - 0.5 * (min + max)) /
    //                              (0.5 * (max - min));
    //    }
}

void p_psNormalizeVector(psVector* myData)
{
    if (myData->type.type == PS_TYPE_F32) {
        p_psNormalizeVectorF32(myData);
    } else if (myData->type.type == PS_TYPE_F64) {
        p_psNormalizeVectorF64(myData);
    } else {
        psError(__func__, "Unalowable data type.\n");
    }
}

/*****************************************************************************
p_psGaussian(myData, myParams): an internal function, used by robust stats,
intended to compute the Gaussian with a specified sigma/mean, at the
specified data point.
 *****************************************************************************/
float p_psGaussian(const psVector* restrict myData, const psVector* restrict myParams)
{
    float x = myData->data.F32[0];
    float mean = myParams->data.F32[0];
    float stdev = myParams->data.F32[1];
    float tmp = exp(-((x - mean) * (x - mean)) / (2.0 * stdev * stdev));

    tmp /= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));

    // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
    return (tmp);
}

/*****************************************************************************
p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which
calculates the specified partial derivative of the above Gaussian function.
 *****************************************************************************/
float p_psGaussianDeriv(const psVector* restrict myData, const psVector* restrict myParams, int whichParam)
{
    float x = myData->data.F32[0];
    float mean = myParams->data.F32[0];
    float stdev = myParams->data.F32[1];
    float tmp = 0.0;

    if (whichParam == 0) {
        // Return the derivative w.r.t. the mean.
        tmp = (x - mean) * p_psGaussian(myData, myParams);
        tmp /= (stdev * stdev);
    } else if (whichParam == 1) {
        // Return the derivative w.r.t. the stdev.
        tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams);
        tmp /= (stdev * stdev * stdev);
    }
    printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp);

    return (tmp);
}

/*****************************************************************************
p_psQuadratic(myData, myParams): an internal function, used by robust stats,
intended to compute a quadratic, with the specified parameters, at the
specified data point.
 *****************************************************************************/
float p_psQuadratic(const psVector* restrict myParams, const psVector* restrict myCoords)
{
    float x = myCoords->data.F32[0];
    float A = myParams->data.F32[0];
    float B = myParams->data.F32[1];
    float C = myParams->data.F32[2];
    float tmp = 0.0;

    tmp = (A * x * x) + (B * x) + C;
    return (tmp);
}

/*****************************************************************************
p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which
calculates the specified partial derivative of the above quadratic function.
 *****************************************************************************/
float p_psQuadraticDeriv(const psVector* restrict myParams,
                         const psVector* restrict myCoords, int whichParamDeriv)
{
    float x = myCoords->data.F32[0];
    float tmp = 0.0;

    if (whichParamDeriv == 0) {
        tmp = x * x;
    } else if (whichParamDeriv == 1) {
        tmp = x;
    } else if (whichParamDeriv == 2) {
        tmp = 1.0;
    }

    return (tmp);
}

/******************************************************************************
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.
 *****************************************************************************/
float p_ps1DPolyMedian(psPolynomial1D* myPoly,
                       float rangeLow,
                       float rangeHigh,
                       float getThisValue)
{
    int numIterations = 0;
    float midpoint = 0.0;
    float oldMidpoint = 1.0;
    float f = 0.0;

    // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue);

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

        f = psPolynomial1DEval(midpoint, myPoly);
        // printf("p_ps1DPolyMedian() iteration %d.  f(%f) is %f\n", numIterations, midpoint, f);
        if (fabs(f - getThisValue) <= FLT_EPSILON) {
            return (midpoint);
        }

        if (f > getThisValue) {
            rangeHigh = midpoint;
        } else {
            rangeLow = midpoint;
        }
        numIterations++;
    }
    return (midpoint);
}

/******************************************************************************
p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given
the specified histogram, with the specified pre-calculated cumulativeSums,
this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum",
and then finds the point x such that f(x) == fitFloat.
 
XXX: This function is currently not being used.
 *****************************************************************************/
float p_psFitQuadratic(psHistogram* histogram,
                       psVector* cumulativeSums,
                       int binNum,
                       float fitFloat)
{
    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);

    if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) {
        x->data.F64[0] = (double)0.5 *
                         (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]);
        x->data.F64[1] = (double)0.5 *
                         (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]);
        x->data.F64[2] = (double)0.5 *
                         (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]);

        y->data.F64[0] = cumulativeSums->data.F32[binNum - 1];
        y->data.F64[1] = cumulativeSums->data.F32[binNum];
        y->data.F64[2] = cumulativeSums->data.F32[binNum + 1];

        if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) {
            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
        }

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

        myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
        return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat));
    } else {
        return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum]));
    }

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

/******************************************************************************
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.)  however there is no currently defined means, in the SDRS, to specify
this computation as a separate step.  If the above robust stat measures were
calcualted in separate functiosn, then much of the initial processing would
be duplicated.
 
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
*****************************************************************************/
void p_psVectorRobustStats(const psVector* restrict myVector,
                           const psVector* restrict maskVector, unsigned int maskVal, psStats* stats)
{
    psHistogram* robustHistogram = NULL;
    psVector* robustHistogramVector = NULL;
    float binSize = 0.0;        // Size of the histogram bins
    int LQBinNum = -1;          // Bin num for lower quartile
    int UQBinNum = -1;          // Bin num for upper quartile
    int i = 0;                  // Loop index variable
    int maxBinNum = 0;
    float maxBinCount = 0.0;
    float dL = 0.0;
    int numBins = 0;
    psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN);

    // psImage* domain;
    // psVector* errors;
    // psVector* data;
    // psVector* initialGuess;
    // psVector* theParams;
    // float chiSq=0.0;
    // float max = -HUGE;
    // float max_pos;
    float myMean = 0.0;
    float myStdev = 0.0;
    float countFloat = 0.0;
    float diff = 0.0;
    float sumSquares = 0.0;
    float sumDiffs = 0.0;
    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);
    psVector* cumulativeRobustSumsFullRange = NULL;
    psVector* cumulativeRobustSumsDlRange = NULL;
    float sumRobust = 0.0;
    float sumN50 = 0.0;
    float sumNfit = 0.0;
    float cumulativeMedian = 0.0;

    // 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;
    p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
    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 = stats->clippedMean;
        }
        if (stats->options & PS_STAT_ROBUST_MEDIAN) {
            stats->robustMedian = stats->clippedMean;
        }
        if (stats->options & PS_STAT_ROBUST_MODE) {
            stats->robustMode = stats->clippedMean;
        }
        if (stats->options & PS_STAT_ROBUST_STDEV) {
            stats->robustStdev = 0.0;
        }
        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
            stats->robustUQ = stats->clippedMean;
            stats->robustLQ = stats->clippedMean;
        }
        // XXX: Set these to the number of unmasked data points?
        stats->robustNfit = 0.0;
        stats->robustN50 = 0.0;
        psFree(tmpStats);
        return;
    }
    // Determine minimum and maximum values in the data vector.
    if (isnan(stats->min)) {
        p_psVectorMin(myVector, maskVector, maskVal, stats);
    }
    if (isnan(stats->max)) {
        p_psVectorMax(myVector, maskVector, maskVal, stats);
    }
    // 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 = (int)((stats->max - stats->min) / binSize);
    robustHistogram = psHistogramAlloc(stats->min, stats->max, numBins);

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

    // Smooth the histogram.
    // XXX: is that the right stdev?
    robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram,
                            tmpStats->clippedStdev / 4.0f);

    // The following was necessary to fit a gaussian to the data, since
    // gaussian functions produce data between 0.0 and 1.0.
    // p_psNormalizeVectorF32(robustHistogramVector);

    /**************************************************************************
    Determine the lower/upper quartiles.
    **************************************************************************/
    // 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
    // now and later in determining the lower/upper quartiles.
    cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0];
    for (i = 1; i < robustHistogramVector->n; i++) {
        cumulativeRobustSumsFullRange->data.F32[i] =
            cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
    }
    sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];

    // Determine the bin number containing the lower quartile point.
    LQBinNum = -1;
    for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) {
        if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) {
            LQBinNum = i;
            break;
        }
    }

    // Determine the bin number containing the upper quartile point.
    UQBinNum = -1;
    for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) {
        if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) {
            UQBinNum = i;
            break;
        }
    }

    if ((LQBinNum == -1) || (UQBinNum == -1)) {
        psAbort(__func__, "Could not determine the robust lower/upper quartiles.");
    }

    /**************************************************************************
    Determine the mode in the range LQ:UQ.
    **************************************************************************/
    // Determine the bin with the peak value in the range LQ to UQ.
    maxBinNum = LQBinNum;
    maxBinCount = robustHistogramVector->data.F32[LQBinNum];
    sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum];
    for (i = LQBinNum + 1; i <= UQBinNum; i++) {
        if (robustHistogramVector->data.F32[i] > maxBinCount) {
            maxBinNum = i;
            maxBinCount = robustHistogramVector->data.F32[i];
        }
        sumN50 += (float)robustHistogram->nums->data.S32[i];
    }

    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
    dL = (UQBinNum - LQBinNum) / 4;

    printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum);

    /**************************************************************************
    Determine the mean/stdev for the bins in the range mode-dL to mode+dL
    **************************************************************************/
    cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    for (i = 0; i < robustHistogramVector->n; i++) {
        cumulativeRobustSumsDlRange->data.F32[i] = 0.0;
    }
    sumNfit = 0.0;
    cumulativeMedian = 0.0;
    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
        if ((0 <= i) && (i < robustHistogramVector->n)) {
            cumulativeRobustSumsDlRange->data.F32[i] =
                cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i];
            cumulativeMedian += robustHistogramVector->data.F32[i];
            sumNfit += (float)robustHistogram->nums->data.S32[i];
        }
    }

    // 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).
    myMean = 0.0;
    for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) {
        if ((0 <= i) && (i < robustHistogramVector->n)) {
            myMean += (robustHistogramVector->data.F32[i]) * 0.5 *
                      (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]);
            countFloat += robustHistogramVector->data.F32[i];
        }
    }
    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 = maxBinNum - dL; i <= maxBinNum + dL; i++) {
        if ((0 <= i) && (i < robustHistogramVector->n)) {
            diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] +
                           robustHistogram->bounds->data.F32[i])) - myMean;
            sumSquares += diff * diff * robustHistogramVector->data.F32[i];
            sumDiffs += diff * robustHistogramVector->data.F32[i];
        }
    }
    myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));

    /**************************************************************************
    Set the appropriate members in the output stats struct.
    **************************************************************************/
    if (stats->options & PS_STAT_ROBUST_MEAN) {
        stats->robustMean = myMean;
    }

    if (stats->options & PS_STAT_ROBUST_MODE) {
        stats->robustMode = 0.5 *
                            (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]);
    }

    if (stats->options & PS_STAT_ROBUST_STDEV) {
        stats->robustStdev = myStdev;
    }
    // To determine the median (and later, the lower/upper quartile), we fit
    // a quadratic to the three bins surrounding the bin containing the median.
    // The quadratic y=f(x) with x being the midpoint of each bin, and y being
    // the cumulative number of data points in all bins up to, and including,
    // this bin.  We then solve the quadratic for

    // XXX: Create a "p_psQuadraticFitAndSolveForX() function.

    if (stats->options & PS_STAT_ROBUST_MEDIAN) {
        if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) {
            x->data.F64[0] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[maxBinNum - 1] +
                              robustHistogram->bounds->data.F32[maxBinNum]);
            x->data.F64[1] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[maxBinNum] +
                              robustHistogram->bounds->data.F32[maxBinNum + 1]);
            x->data.F64[2] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[maxBinNum + 1] +
                              robustHistogram->bounds->data.F32[maxBinNum + 2]);

            y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1];
            y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum];
            y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1];

            // Ensure that cumulativeMedian/2 is actually within the range of the bins
            // we are using.
            cumulativeMedian *= 0.5;
            if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) {
                printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]);
                psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n");
            }
            // XXX: yErr is not currently used by psVectorFitPolynomial1D().  We
            // may have to set this meaningfully later.
            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);
            // Call p_ps1DPolyMedian(), which does a binary search on the
            // polynomial, looking for the value x such that
            // f(x) = cumulativeMedian.
            stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian);
        } else {
            // If the mode is the first/last histogram bin, then simply use
            // the midpoint of that bin.
            stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] +
                                         robustHistogram->bounds->data.F32[maxBinNum]);
        }
    }
    // The lower/upper quartile calculations are very similar to the median
    // calculations.  We fit a quadratic to the array containing the
    // cumulative data points in each bin, then search for x such that
    // f(x) equals the lower/upper quartile exactly.
    //
    if (stats->options & PS_STAT_ROBUST_QUARTILE) {
        countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1];

        if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) {
            x->data.F64[0] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[LQBinNum - 1] +
                              robustHistogram->bounds->data.F32[LQBinNum]);
            x->data.F64[1] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[LQBinNum] +
                              robustHistogram->bounds->data.F32[LQBinNum + 1]);
            x->data.F64[2] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[LQBinNum + 1] +
                              robustHistogram->bounds->data.F32[LQBinNum + 2]);

            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1];
            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum];
            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1];

            if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) {
                psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n");
            }

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

            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
            stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0);

        } else {
            // If the LQ is the first/last histogram bin, then simply use
            // the midpoint of that bin.
            stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] +
                                     robustHistogram->bounds->data.F32[LQBinNum]);
        }

        if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) {
            x->data.F64[0] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[UQBinNum - 1] +
                              robustHistogram->bounds->data.F32[UQBinNum]);
            x->data.F64[1] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[UQBinNum] +
                              robustHistogram->bounds->data.F32[UQBinNum + 1]);
            x->data.F64[2] = (double)0.5 *
                             (robustHistogram->bounds->data.F32[UQBinNum + 1] +
                              robustHistogram->bounds->data.F32[UQBinNum + 2]);

            y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1];
            y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum];
            y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1];

            if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) &&
                    ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) {
                psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n");
            }

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

            myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr);
            stats->robustUQ = p_ps1DPolyMedian(myPoly,
                                               x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0);
        } else {
            // If the UQ is the first/last histogram bin, then simply use
            // the midpoint of that bin.
            stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] +
                                     robustHistogram->bounds->data.F32[UQBinNum]);
        }
    }
    stats->robustNfit = sumNfit;
    stats->robustN50 = sumN50;

    psFree(x);
    psFree(y);
    psFree(yErr);
    psFree(tmpStats);
    psFree(robustHistogram);
    psFree(myPoly);
    psFree(cumulativeRobustSumsFullRange);
    psFree(cumulativeRobustSumsDlRange);
}

/*
void p_ps_FitTheGaussian()
{
    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    // NOTE: This code uses the psMinimize.c functions to perform
    // the fit.  It doesn't quite work, 100% of the time.  For the time being
    // I am commenting this code out, and replacing it with code which
    // calculates the mean directly on the robustHistogram.
    //
    // XXX: Replace this with the LaGrange interpolation stuff.
 
    domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32);
    errors = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    data = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    initialGuess = psVectorAlloc(2, PS_TYPE_F32);
 
    max = -HUGE;
    max_pos = -1;
    for (i=0;i<robustHistogramVector->n;i++) {
        domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] +
                                        robustHistogram->bounds->data.F32[i]);
        data->data.F32[i] = (float) robustHistogramVector->data.F32[i];
        errors->data.F32[i] = 1.0;
        //printf("DATA (%.2f, %.2f)\n", domain->data.F32[i][0], data->data.F32[i]);
        if (data->data.F32[i] > max) {
            max = data->data.F32[i];
            max_pos = domain->data.F32[i][0];
        }
    }
 
    initialGuess->data.F32[0] = max_pos;
    initialGuess->data.F32[1] = 1.0;
    printf("Initial (mean. stdev) is (%.3f, %.3f)\n", initialGuess->data.F32[0], initialGuess->data.F32[1]);
 
    printf("Calling psMinimizeChi2()\n");
    theParams = psMinimizeChi2(p_psGaussian,
                               p_psGaussianDeriv,
                               domain,
                               data,
                               errors,
                               initialGuess,
                               NULL,
                               &chiSq);
    printf("Called psMinimizeChi2()\n");
    printf("p_psVectorRobustStats() is (%f, %f)\n", theParams->data.F32[0], theParams->data.F32[1]);
}
*/

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

/* 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 = 0;
    newStruct->robustNfit = 0;
    newStruct->clippedMean = NAN;
    newStruct->clippedStdev = NAN;
    newStruct->clipSigma = 3.0;
    newStruct->clipIter = 3.0;
    newStruct->min = NAN;
    newStruct->max = NAN;
    newStruct->binsize = NAN;
    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(float lower, float upper, int n)
{
    int i = 0;                  // Loop index variable
    psHistogram* newHist = NULL;        // The new histogram structure
    float binSize = 0.0;        // The histogram bin size

    // NOTE: Verify that this is the correct action.
    if (n == 0) {
        return (NULL);
    }
    if (n < 0) {
        psAbort(__func__, "psHistogramAlloc() called with bin size %d.\n", n);
    }
    // NOTE: Verify that this is the correct action.
    if (lower > upper) {
        return (NULL);
    }
    // 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));
    p_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) / (float)n;
    // NOTE: 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 * (float)i);
    }

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

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

    return (newHist);
}

/******************************************************************************
psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure
with the specifed bounds.
 
Input:
    bounds
Returns:
    The histogram structure
 *****************************************************************************/
psHistogram* psHistogramAllocGeneric(const psVector* restrict bounds)
{
    psHistogram* newHist = NULL;        // The new histogram structure
    int i;                      // Loop index variable

    // NOTE: Verify that this is the correct action.
    if (bounds == NULL) {
        // psAbort(__func__, "psHistogram requested with NULL bounds");
        return (NULL);
    }
    // NOTE: Verify that this is the correct action.
    if (bounds->n <= 1) {
        // psAbort(__func__, "psHistogram requested with NULL bounds");
        return (NULL);
    }

    if (bounds->type.type != PS_TYPE_F32) {
        // psAbort(__func__, "psHistogram request a bound which is not type F32");
        return (NULL);
    }
    // Allocate memory for the new histogram structure.
    newHist = (psHistogram* ) psAlloc(sizeof(psHistogram));
    p_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_U32);
    newHist->nums->n = newHist->nums->nalloc;
    for (i = 0; i < newHist->nums->n; i++) {
        newHist->nums->data.U32[i] = 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);
}

/*****************************************************************************
psVectorHistogram(out, in, 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* restrict in,
                               const psVector* restrict mask,
                               unsigned int maskVal)
{
    int i = 0;                  // Loop index variable
    float binSize = 0.0;        // Histogram bin size
    int binNum = 0;             // A temporary bin number
    int numBins = 0;            // The total number of bins
    int tmp = 0;
    psScalar tmpScalar;
    tmpScalar.type.type = PS_TYPE_F32;

    // NOTE: Verify that this is the correct action.
    if (out == NULL) {
        return (NULL);
    }
    // Check the specified output histogram for type psF32
    if (out->bounds->type.type != PS_TYPE_F32) {
        psAbort(__func__, "Only data type PS_TYPE_F32 for the output.bounds member.");
    }

    if (out->nums->type.type != PS_TYPE_U32) {
        psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member.");
    }
    // NOTE: Verify that this is the correct action.
    if (in == NULL) {
        return (out);
    }

    if (in->type.type != PS_TYPE_F32) {
        psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type);
    }

    if (mask != NULL) {
        if (in->n != mask->n) {
            psAbort(__func__, "Vector data and vector mask are of different sizes.");
        }
        if (mask->type.type != PS_TYPE_U8) {
            psAbort(__func__, "Vector mask must be type PS_TYPE_U8");
        }
    }

    // NOTE: determine the correct action for a variety of other cases:
    // in vector has 0 elements, and histogram structure has zero bins.

    numBins = out->nums->n;
    for (i = 0; i < in->n; i++) {
        // Check if this pixel is masked, and if so, skip it.
        if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) {
            // Check if this pixel is below the minimum value, and if so
            // count it, then skip it.
            if (in->data.F32[i] < out->bounds->data.F32[0]) {
                out->minNum++;
                // Check if this pixel is above the maximum value, and if so
                // count it, then skip it.
            } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) {
                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 = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize);

                    // NOTE: This next if-statement really shouldn't be necessary.
                    // However, do to numerical lack of precision, we occasionally
                    // produce a binNum outside the range of bins.
                    if (binNum >= out->nums->n) {
                        binNum = out->nums->n - 1;
                    }

                    (out->nums->data.S32[binNum])++;

                    // If this is a non-uniform histogram, determining the correct
                    // bin number requires a bit more work.
                } else {
                    tmpScalar.data.F32 = in->data.F32[i];
                    tmp = p_psVectorBinDisect(out->bounds,
                                              &tmpScalar);
                    if (tmp == -1) {
                        //                      XXX: Generate warning message
                    } else {
                        (out->nums->data.S32[tmp])++;
                    }
                }
            }
        }
    }
    return (out);
}

/******************************************************************************
p_psConvertToF32(myVector, maskVector, maskVal, stats): 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 Athlons...
 *****************************************************************************/
psVector* p_psConvertToF32(psStats* stats, psVector* in, psVector* mask, unsigned int maskVal)
{
    int 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] = (float)in->data.S8[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] = (float)in->data.U16[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] = (float)in->data.F64[i];
        }
    } else if (in->type.type == PS_TYPE_F32) {
        // do nothing
    } else {
        psAbort(__func__, "unsupported vector type 0x%x\n", in->type.type);
    }
    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.
 
NOTE: The current strategy is to implement everything assuming that all
input data is of type PS_TYPE_F32.  Once the basic code is in place, we will
macro-ize everything and add PS_TYPE_U16 and PS_TYPE_F64.
 *****************************************************************************/
psStats* psVectorStats(psStats* stats,
                       psVector* in,
                       psVector* mask,
                       unsigned int maskVal)
{
    psVector* inF32;
    int mustFreeTmp = 1;

    // NOTE: Verify that this is the correct action.
    if (in == NULL) {
        return (stats);
    }
    if (stats == NULL) {
        return (NULL);
    }

    inF32 = p_psConvertToF32(stats, in, mask, maskVal);
    if (inF32 == NULL) {
        inF32 = in;
        mustFreeTmp = 0;
    }
    // XXX: Should we abort if (stats->min == stats->max)?
    if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) {
        psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max);
    }

    if (mask != NULL) {
        PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(mask, in);
        PS_CHECK_VECTOR_TYPE_RETURN_NULL(mask, PS_TYPE_U8);
    }
    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEAN) {
        p_psVectorSampleMean(inF32, mask, maskVal, stats);
    }
    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
        p_psVectorSampleMedian(inF32, mask, maskVal, stats);
    }
    // ************************************************************************
    // NOTE: The Stdev calculation requires the mean.  Should we assume the
    // mean has already been calculated? Or should we always calculate it?
    if (stats->options & PS_STAT_SAMPLE_STDEV) {
        p_psVectorSampleMean(inF32, mask, maskVal, stats);
        p_psVectorSampleStdev(inF32, mask, maskVal, stats);
    }
    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
        p_psVectorSampleQuartiles(inF32, mask, maskVal, stats);
    }
    // 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)) {
        p_psVectorRobustStats(inF32, mask, maskVal, stats);
    }

    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
        p_psVectorClippedStats(inF32, mask, maskVal, stats);
    }
    // ************************************************************************
    if (stats->options & PS_STAT_MAX) {
        p_psVectorMax(inF32, mask, maskVal, stats);
    }
    // ************************************************************************
    if (stats->options & PS_STAT_MIN) {
        p_psVectorMin(inF32, mask, maskVal, stats);
    }

    if (mustFreeTmp == 1) {
        psFree(inF32);
    }
    return (stats);
}
