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

#include "psMemory.h"
#include "psVector.h"
#include "psTrace.h"
#include "psError.h"
#include "psAbort.h"
#include "psStats.h"

#include "float.h"
#define DEFAULT_ROBUST_SIZE_THRESHOLD 30000   // Vectors that are large than this
// 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 MYMAXFLOAT HUGE

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

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->sampleLimit = DEFAULT_ROBUST_SIZE_THRESHOLD;
    newStruct->robustMean = NAN;
    newStruct->robustMedian = NAN;
    newStruct->robustMode = NAN;
    newStruct->robustStdev = NAN;
    newStruct->robustUQ = NAN;
    newStruct->robustLQ = NAN;
    newStruct->robustN50 = NAN;
    newStruct->robustNfit = NAN;
    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);
    }

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

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

    return(newHist);
}

/******************************************************************************
psHistogramAlloc(lower, upper, n): allocate a non-uniform histogram structure
with the specifed bounds.  The number of elements in the bounds vector is n.
Therefore, the number of bins is n-1.
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);
    }

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


/*****************************************************************************
psHistogramVector(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 *psHistogramVector(psHistogram *out,
                               const psVector *restrict in,
                               const psVector *restrict mask,
                               int maskVal)
{
    int i = 0;                                  // Loop index variable
    int j = 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

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

    // 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);
                    (out->nums->data.S32[binNum])++;
                    // If this is a non-uniform histogram, determining the correct
                    // bin number requires a bit more work.
                } else {
                    // NOTE: This is slow.  Put a smarter algorithm here to
                    // find the correct bin number (bin search, probably)
                    for (j=0;j<(out->bounds->n)-1;j++) {
                        if ((out->bounds->data.S32[j] <= in->data.F32[i]) &&
                                (in->data.F32[i] <= out->bounds->data.S32[j+1])) {
                            (out->nums->data.S32[j])++;
                        }
                    }
                }
            }
        }
    }
    return(out);
}

/*****************************************************************************
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
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 = -MYMAXFLOAT;                    // 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 = MYMAXFLOAT;                     // 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 {
        rangeMin = stats->min;
        rangeMax = stats->max;
        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
    psStats *stats2 = NULL;                     // Temporary stats structure


    // Determine if the number of data points exceed a threshold which will
    // cause to generate robust stats, as opposed to exact stats.

    if (myVector->n > stats->sampleLimit) {
        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");

        // 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.
    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.
 *****************************************************************************/
void p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
                                  float sigma)
{
    int i = 0;                                  // Loop index variable
    int j = 0;                                  // Loop index variable
    float tmpf = 0.0;                           // Temporary variable
    float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients

    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
            tmpf = (float) sqrt(2.0f * M_PI * sigma * sigma);
            #else

            tmpf = sqrtf(2.0f * M_PI * sigma * sigma);
            #endif

            gaussianCoefs[i] = (float) exp( (-((float) (i-GAUSS_WIDTH)) *
                                             ((float) (i-GAUSS_WIDTH))) /
                                            (2.0f * sigma * sigma)) / tmpf;
        } 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.
             */
            return;
        }
    }

    for(i=0;i<robustHistogram->nums->n;i++) {
        for (j=-GAUSS_WIDTH;j<=+GAUSS_WIDTH;j++) {
            if (((j+i) >= 0) && ((j+i) < robustHistogram->nums->n)) {
                robustHistogram->nums->data.S32[j+i]+=
                    (gaussianCoefs[j+GAUSS_WIDTH] *
                     (float) robustHistogram->nums->data.S32[j+i]);
            }
        }
    }
}

/******************************************************************************
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 if the number of data points exceed a threshold which will
    // cause to generate robust stats, as opposed to exact stats.
    if (myVector->n > stats->sampleLimit) {
        psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
        psStats *stats2 = NULL;
        // Calculate the robust quartiles.
        stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
        p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);

        // Store the robust quartiles into the sample quartile members.
        stats->sampleUQ = stats2->robustUQ;
        stats->sampleLQ = stats2->robustLQ;

        // 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 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);
    // NOTE: This is the
}


/******************************************************************************
p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
stdev of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 
NOTE: the mean is always calculated exactly.  Robust means are never
calculated in this routine.
 *****************************************************************************/
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_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;
    float binSize = 0.0;                        // Size of the histogram bins
    float sigmaE = 0.0;
    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;
    int maxBinCount = 0;
    float dL = 0.0;
    int numBins = 0;
    psStats *tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV|PS_STAT_CLIPPED_MEAN);

    // NOTE: The SDRS states that the sample quartiles must be used to
    // determine the initial bin sizes.  However, the sample quartiles are
    // calculated based on a full sort of the data set, regardless of the
    // size of the data set.  We should consult with IfA to ensure that this
    // is really required.
    /*
        if (isnan(stats->sampleUQ) ||
            isnan(stats->sampleLQ)) {
            stats->options = stats->options | PS_STAT_SAMPLE_QUARTILE;
            p_psVectorSampleQuartiles(myVector,
                                      maskVector,
                                      maskVal,
                                      stats);
        }
    */
    // Compute the initial bin size of the robust histogram.
    p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats);
    binSize = stats->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;
        }
        psFree(tmpStats);
        psFree(robustHistogram);
        return;
    }

    // Detemine 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 (yes, 2 is necessary, not 1).  Also,
    // if we get here, we know that binSize != 0.0.
    numBins = 2 + (int) ((stats->max - stats->min) / binSize);

    robustHistogram = psHistogramAlloc(stats->min,
                                       stats->max,
                                       numBins);
    // Populate the histogram array.
    psHistogramVector(robustHistogram, myVector, maskVector, maskVal);

    // Smooth the histogram.
    p_psVectorsmoothHistGaussian(robustHistogram, sigmaE/4.0f);

    LQBinNum = -1;
    UQBinNum = -1;
    for (i=0;i<robustHistogram->nums->n;i++) {
        if ((robustHistogram->nums->data.S32[i] <= stats->sampleLQ) &&
                (stats->sampleLQ <= robustHistogram->nums->data.S32[i])) {
            LQBinNum = i;
        }

        if ((robustHistogram->nums->data.S32[i] <= stats->sampleUQ) &&
                (stats->sampleUQ <= robustHistogram->nums->data.S32[i])) {
            UQBinNum = i;
        }
    }

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

    dL = (stats->robustUQ - stats->robustLQ) / 8.0;

    // Fit a Gaussian to the bins in the range MODE-dL to Mode+dL
    // NOTE: This step is dependent on the functions in psMinimize.c being
    // implemented.  Currently, they are not.

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

    if  (stats->options & PS_STAT_ROBUST_MEDIAN) {
        stats->robustMedian = 0.0;
    }
    if  (stats->options & PS_STAT_ROBUST_MODE) {
        stats->robustMode = maxBinNum;
    }
    if  (stats->options & PS_STAT_ROBUST_STDEV) {
        stats->robustStdev = 0.0;
    }
    if  (stats->options & PS_STAT_ROBUST_QUARTILE) {
        stats->robustUQ = 0.0;
        stats->robustLQ = 0.0;
    }
    stats->robustNfit = 0.0;
    stats->robustN50 = 0.0;
    psFree(tmpStats);
    psFree(robustHistogram);
}


/******************************************************************************
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)
{
    // NOTE: Verify that this is the correct action.
    if (in == NULL) {
        return(stats);
    }
    if (stats == NULL) {
        return(NULL);
    }

    // Ensure that the data is of type PS_TYPE_F32.  Eventually, more data
    // types will be implemented.
    if (in->type.type != PS_TYPE_F32) {
        psAbort(__func__,
                "Only data type PS_TYPE_F32 is currently supported (0x%x).",
                in->type.type);
    }

    // Ensure that the mask vector is of the proper size and 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");
        }
    }

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEAN) {
        p_psVectorSampleMean(in, mask, maskVal, stats);
    }

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
        p_psVectorSampleMedian(in, 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(in, mask, maskVal, stats);
        p_psVectorSampleStdev(in, mask, maskVal, stats);
    }

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
        p_psVectorSampleQuartiles(in, 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(in, mask, maskVal, stats);
    }

    if ((stats->options & PS_STAT_CLIPPED_MEAN) ||
            (stats->options & PS_STAT_CLIPPED_STDEV)) {
        p_psVectorClippedStats(in, mask, maskVal, stats);
    }

    // ************************************************************************
    if (stats->options & PS_STAT_MAX) {
        p_psVectorMax(in, mask, maskVal, stats);
    }

    // ************************************************************************
    if (stats->options & PS_STAT_MIN) {
        p_psVectorMin(in, mask, maskVal, stats);
    }

    return(stats);
}
