/** @file  psStats.c
 *  \brief basic statistical operations
 *  @ingroup Stat
 *
 *  This file will hold the definitions of the histogram and stats data
 *  structures.  It also contains prototypes for procedures which operate
 *  on those data structures.
 *
 *  @author GLG, MHPCC
 *
 *  XXX: The following stats members are never used, or set in this code.
 *      stats->clippedNvalues
 *
 * XXX: Must do
 * nSubsample points
 * use ->min and ->max (PS_STAT_USE_RANGE)
 *
 *  @version $Revision: 1.167 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-02-03 22:05:22 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

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

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

#include "psErrorText.h"

/*****************************************************************************/
/* DEFINE STATEMENTS                                                         */
/*****************************************************************************/
#define PS_GAUSS_WIDTH 5                // The width of the Gaussian smoothing.
// This corresponds to N in the ADD.
#define PS_CLIPPED_NUM_ITER_LB 1
#define PS_CLIPPED_NUM_ITER_UB 10
#define PS_CLIPPED_SIGMA_LB 1.0
#define PS_CLIPPED_SIGMA_UB 10.0
#define PS_POLY_MEDIAN_MAX_ITERATIONS 30

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

// None

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

// None

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

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

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

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

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

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

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

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

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

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

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

/******************************************************************************
MISC PRIVATE STATISTICAL FUNCTIONS
 
NOTE: it is assumed that any call to these statistical functions will have
been preceded by a call to the psVectorStats() function.  Various sanity tests
will only be performed in psVectorStats().  Should we perform the sanity
checks in each routine anyway?
 
XXX: For many of these private stats routines, what should be done if there
are no acceptable elements in the input vector (if no elements lie within
range, or there are no unmasked elements, or the input vector is NULL)?
Currently we set the value to NAN.
 
XXX: Optimization: many routines have an "empty" boolean variable which keeps
track of whether or not the vector has any valid elements.  This code can
possibly be optimized away.
 *****************************************************************************/
/******************************************************************************
p_psVectorSampleMean(myVector, maskVector, maskVal, stats): calculates the
mean of the input vector.  If there was a problem with the mean calculation,
this routine sets stats->sampleMean to NAN.
 *****************************************************************************/
psS32 p_psVectorSampleMean(const psVector* myVector,
                           const psVector* errors,
                           const psVector* maskVector,
                           psU32 maskVal,
                           psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);

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

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

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

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

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

}

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

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

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

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

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

/******************************************************************************
p_psVectorNValues(myVector, maskVector, maskVal, stats): This routine returns
"true" if the inputPsVector has 1 or more valid elements (a valid element is an
unmasked element within the specified min/max range).  Otherwise, return
"false".
 
XXX: Can you use psVectorCountPixelMask here?
 *****************************************************************************/
bool p_psVectorCheckNonEmpty(const psVector* myVector,
                             const psVector* maskVector,
                             psU32 maskVal,
                             psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    PS_ASSERT_VECTOR_NON_NULL(myVector, -1);
    PS_ASSERT_PTR_NON_NULL(stats, -1);
    psS32 i = 0;                // Loop index variable

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


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

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

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

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

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

    if (nValues <= 0) {
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n");
        return(true);
        stats->sampleMedian = NAN;
    }

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

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

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

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

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

/******************************************************************************
p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input
robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
smoothed data.
 
XXX: Only PS_TYPE_F32 is supported.
 
XXX: Write a general routine which smoothes a psVector.  This routine should
call that.  Is that possible?
 *****************************************************************************/
psVector *p_psVectorSmoothHistGaussian(
    psHistogram *histogram,
    psF32 sigma)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psTrace(__func__, 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma);
    PS_ASSERT_PTR_NON_NULL(histogram, NULL);
    PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL);
    PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL);
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(histogram->nums);
    }

    psS32 numBins = histogram->nums->n;
    psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32);
    psS32 jMin = 0;
    psS32 jMax = 0;

    if (histogram->uniform == false) {
        //
        // We get here if the histogram is non-uniform.  This code is not tested.
        // However, it is also not used anywhere, yet.
        //
        psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n");

        for (psS32 i = 0; i < numBins; i++) {
            // Determine the midpoint of bin i.
            psS32 iMid = PS_BIN_MIDPOINT(histogram, i);

            //
            // We determine the bin numbers (jMin:jMax) corresponding to a
            // range of data values surrounding iMid.  The range is of size:
            // 2*PS_GAUSS_WIDTH*sigma
            psF32 x = iMid - (PS_GAUSS_WIDTH * sigma);
            jMin = i;
            while (jMin > 0) {
                if ((histogram->bounds->data.F32[jMin] <= x) &&
                        (histogram->bounds->data.F32[jMin+1] >= x)) {
                    break;
                }
                jMin--;
            }

            x = iMid + (PS_GAUSS_WIDTH * sigma);
            jMax = i;
            while (jMax < (histogram->bounds->n - 1)) {
                if ((histogram->bounds->data.F32[jMax] <= x) &&
                        (histogram->bounds->data.F32[jMax+1] >= x)) {
                    break;
                }
                jMax++;
            }


            //
            // Loop from jMin to jMax, computing the gaussian of data i.
            //
            smooth->data.F32[i] = 0.0;
            for (psS32 j = jMin ; j <= jMax ; j++) {
                psS32 jMid = PS_BIN_MIDPOINT(histogram, j);
                smooth->data.F32[i] +=
                    histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true);
            }
        }
    } else {
        //
        // We get here if the histogram is uniform.
        //
        for (psS32 i = 0; i < numBins; i++) {
            psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0];
            psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize);

            //
            // XXX: The following is wrong.  However, in practice, the sigma was too
            // large compared to the binSize.  This meant that the smoothing was done
            // over 500 bins in the robust stats algorithm.  This mean that the smoothing
            // took way too long.
            //
            gaussWidth = 10;
            //
            // We determine the bin numbers (jMin:jMax) corresponding to a
            // range of data values surrounding iMid.  The range is of size:
            // 2*PS_GAUSS_WIDTH*sigma
            //
            psS32 jMin = i - gaussWidth;
            if (jMin < 0 ) {
                jMin = 0;
            }
            psS32 jMax = i + gaussWidth;
            if (jMax > (histogram->bounds->n - 1)) {
                jMax = (histogram->bounds->n - 1);
            }

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

    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(smooth);
    }
    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
    return(smooth);
}
/******************************************************************************
p_psVectorSampleQuartiles(myVector, maskVector, maskVal, statsalculates
the upper and/or lower quartiles of the input vector.
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    NULL
 *****************************************************************************/
bool p_psVectorSampleQuartiles(const psVector* myVector,
                               const psVector* maskVector,
                               psU32 maskVal,
                               psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psVector* unsortedVector = NULL;
    psVector* sortedVector = NULL;
    psS32 i = 0;
    psS32 count = 0;
    psS32 nValues = 0;

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

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

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

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

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

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

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

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

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

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

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

/******************************************************************************
p_psVectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the
clipped stats (mean or stdev) of the input vector.
 
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    0 for success.
    -1: error
    -2: warning
 
XXX: Use static vectors for tmpMask.
 *****************************************************************************/
psS32 p_psVectorClippedStats(const psVector* myVector,
                             const psVector* errors,
                             const psVector* maskVector,
                             psU32 maskVal,
                             psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psTrace(__func__, 4, "Trace level is %d\n", psTraceGetLevel(__func__));
    psF32 clippedMean = 0.0;
    psF32 clippedStdev = 0.0;
    psVector* tmpMask = NULL;
    static psStats *statsTmp = NULL;
    psS32 rc = 0;

    // Ensure that stats->clipIter is within the proper range.
    PS_ASSERT_INT_WITHIN_RANGE(stats->clipIter,
                               PS_CLIPPED_NUM_ITER_LB,
                               PS_CLIPPED_NUM_ITER_UB, -1);

    // Ensure that stats->clipSigma is within the proper range.
    PS_ASSERT_INT_WITHIN_RANGE(stats->clipSigma,
                               PS_CLIPPED_SIGMA_LB,
                               PS_CLIPPED_SIGMA_UB, -1);

    // Allocate a psStats structure for calculating the mean, median, and
    // stdev.
    if (statsTmp == NULL) {
        statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
        p_psMemSetPersistent(statsTmp, true);
    } else {
        // Initialize structure if already allocated
        statsTmp->sampleMean = NAN;
        statsTmp->sampleMedian = NAN;
        statsTmp->sampleStdev = NAN;
        statsTmp->sampleUQ = NAN;
        statsTmp->sampleLQ = NAN;
        statsTmp->robustMedian = NAN;
        statsTmp->robustStdev = NAN;
        statsTmp->robustUQ = NAN;
        statsTmp->robustLQ = NAN;
        statsTmp->robustN50 = -1;
        statsTmp->fittedMean = NAN;
        statsTmp->fittedStdev = NAN;
        statsTmp->fittedNfit = -1;
        statsTmp->clippedMean = NAN;
        statsTmp->clippedStdev = NAN;
        statsTmp->clippedNvalues = -1;
        statsTmp->clipSigma = 3.0;
        statsTmp->clipIter = 3;
        statsTmp->min = NAN;
        statsTmp->max = NAN;
        statsTmp->binsize = NAN;
        statsTmp->nSubsample = 100000;
    }

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

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

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

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

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

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

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

        // b) compute new mean and stdev
        p_psVectorSampleMean(myVector, errors, tmpMask, 0xff, statsTmp);
        p_psVectorSampleStdev(myVector, errors, tmpMask, 0xff, statsTmp);
        psTrace(__func__, 6, "The new sample mean is %f\n", statsTmp->sampleMean);
        psTrace(__func__, 6, "The new sample stdev is %f\n", statsTmp->sampleStdev);

        // If the new mean and stdev are NAN, we must exit the loop.
        // Otherwise, use the new results and continue.
        if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) {
            // Exit loop.
            iter = stats->clipIter;
            psError(PS_ERR_UNKNOWN, true, "p_psVectorSampleMean() or p_psVectorSampleStdev() returned a NAN.\n");
            psFree(tmpMask);
            return(-1);
        } else {
            clippedMean = statsTmp->sampleMean;
            clippedStdev = statsTmp->sampleStdev;
        }
    }

    // 7. The last calcuated value of x is the cliped mean.
    if (stats->options & PS_STAT_CLIPPED_MEAN) {
        stats->clippedMean = clippedMean;
        psTrace(__func__, 6, "The final clipped mean is %f\n", clippedMean);
    }
    // 8. The last calcuated value of stdev is the cliped stdev.
    if (stats->options & PS_STAT_CLIPPED_STDEV) {
        stats->clippedStdev = clippedStdev;
        psTrace(__func__, 6, "The final clipped stdev is %f\n", clippedStdev);
    }

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

static psF32 QuadraticInverse(
    psF32 a,
    psF32 b,
    psF32 c,
    psF32 y,
    psF32 xLo,
    psF32 xHi)
{
    psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a));

    psF64 x1 = -b/(2.0*a) + tmp;
    psF64 x2 = -b/(2.0*a) - tmp;

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

    if ((xLo <= x1) && (x1 <= xHi)) {
        return(x1);
    } else if ((xLo <= x2) && (x2 <= xHi)) {
        return(x2);
    } else {
        return(0.5 * (xLo + xHi));
    }
}


/******************************************************************************
fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general
routine which fits a quadratic to three points and returns the x-value
corresponding to the input y-value.  This routine takes psVectors of x/y pairs
as input, and fits a quadratic to the 3 points surrounding element binNum in
the vectors (the midpoint between element i and i+1 is used for x[i]).  It
then determines for what value x does that quadratic f(x) = yVal (the input
parameter).
 
*****************************************************************************/
psF32 fitQuadraticSearchForYThenReturnX(
    psVector *xVec,
    psVector *yVec,
    psS32 binNum,
    psF32 yVal)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    psTrace(__func__, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(xVec);
        PS_VECTOR_PRINT_F32(yVec);
    }

    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
    //    PS_ASSERT_VECTORS_SIZE_EQUAL(xVec, yVec, NAN);
    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (xVec->n - 1), NAN);
    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (yVec->n - 1), NAN);

    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
    psF32 tmpFloat = 0.0f;

    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
        // The general case.  We have all three points.
        x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
        x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1]));
        x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2]));
        y->data.F64[0] = yVec->data.F32[binNum - 1];
        y->data.F64[1] = yVec->data.F32[binNum];
        y->data.F64[2] = yVec->data.F32[binNum + 1];
        psTrace(__func__, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1],
                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
        psTrace(__func__, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
        psTrace(__func__, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);

        //
        // Ensure that the y values are monotonic.
        //
        // XXX: This routine should probably be rewritten in a more general fashion
        // so that the following checks are not necessary.
        //
        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
                ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
        }

        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
                ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
            psError(PS_ERR_UNKNOWN, true, "This routine must be called with monotonically increasing or decreasing data points.\n");
            psFree(x);
            psFree(y);
            psTrace(__func__, 5, "---- %s() end ----\n", __func__);
            return(NAN);
        }

        // Determine the coefficients of the polynomial.
        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);
        if (myPoly == NULL) {
            psError(PS_ERR_UNEXPECTED_NULL, false,
                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT);
            psFree(x);
            psFree(y);
            psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
            return(NAN);
        }
        psTrace(__func__, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
        psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
        psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
        psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n",
                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));

        psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
        tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
        psFree(myPoly);

        if (isnan(tmpFloat)) {
            psError(PS_ERR_UNEXPECTED_NULL,
                    false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
            psFree(x);
            psFree(y);
            psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
            return(NAN);
        }
    } else {
        // These are special cases where the bin is at the beginning or end of the vector.
        if (binNum == 0) {
            // We have two points only at the beginning of the vectors x and y.
            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                              xVec->data.F32[binNum + 1]);
        } else if (binNum == (xVec->n - 1)) {
            // The special case where we have two points only at the end of
            // the vectors x and y.
            // XXX: Is this right?
            tmpFloat = xVec->data.F32[binNum];
        } else if (binNum == (xVec->n - 2)) {
            // XXX: Is this right?
            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                              xVec->data.F32[binNum + 1]);
        }
    }

    psTrace(__func__, 6, "FIT: return %f\n", tmpFloat);
    psFree(x);
    psFree(y);

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

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

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

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

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

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

/******************************************************************************
p_psVectorRobustStats(myVector, maskVector, maskVal, stats): This is the new
version of the robust stats routine.
 
XXX: MUST DO: If the errors in the input values are known, then the same
approach is used, except that the histograms become probability density
functions (PDFs). In this case, the input values are spread out, so that they
do not simply contribute a single unit to the histogram, but rather contribute
a fraction of a value, equivalent to the weight. In the interests of speed, a
boxcar PDF may be used to represent each input value (as opposed to a
Gaussian), where the boxcar width is equal to 2p2 ln 2 times the error and
each input value contributes constant area. Then the robust median and
standard deviation are estimated in the same manner as above.
 
XXX: Check for errors in psLib routines that we call.
 
XXX: Review and ensure that all memory is free'ed at premature exits.
*****************************************************************************/
#define INITIAL_NUM_BINS 1000.0
psS32 p_psVectorRobustStats(const psVector* myVector,
                            const psVector* errors,
                            const psVector* maskVector,
                            psU32 maskVal,
                            psStats* stats)
{
    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
    if (psTraceGetLevel(__func__) >= 8) {
        PS_VECTOR_PRINT_F32(myVector);
    }
    psHistogram *robustHistogram = NULL;
    psHistogram *cumulativeRobustHistogram = NULL;
    psS32 numBins = 0;
    psScalar tmpScalar;
    tmpScalar.type.type = PS_TYPE_F32;
    psS32 totalDataPoints = 0;
    psS32 rc = 0;
    psS32 rcBool = false;
    psVector *tmpMaskVec = psVectorAlloc(myVector->n, PS_TYPE_U8);
    if (maskVector != NULL) {
        for (psS32 i = 0 ; i < myVector->n ; i++) {
            tmpMaskVec->data.U8[i] = maskVector->data.U8[i];
        }
    } else {
        for (psS32 i = 0 ; i < myVector->n ; i++) {
            tmpMaskVec->data.U8[i] = 0;
        }
    }
    psS32 iterate = 1;
    psF32 sigma;
    psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX);

    while (iterate > 0) {
        psTrace(__func__, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
        psF32 binSize = 0.0;
        if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) {
            // Set initial bin size to the specified value.
            binSize = stats->binsize;
            psTrace(__func__, 6, "Setting initial robust bin size to %.2f\n", binSize);
        } else {
            // Determine the bin size of the robust histogram.  This is done
            // by computing the total range of data values and dividing by 1000.0.

            rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax);
            rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax);
            if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
                psFree(tmpStatsMinMax);
                psFree(tmpMaskVec);
                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
                return(1);
            }
            psTrace(__func__, 6, "Data min/max is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
            binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / INITIAL_NUM_BINS;
        }
        psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize);

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

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

        //
        // ADD: Step 0.
        // Construct the histogram with the specified bin size.
        //
        // NOTE: we can not specify the bin size precisely since the argument
        // to psHistogramAlloc() is the number of bins, not the binSize.
        // If we get here, we know that binSize != 0.0.
        //
        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize);
        psTrace(__func__, 6, "Numbins is %d\n", numBins);
        psTrace(__func__, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
        robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
        cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);

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

        //
        // ADD: Step 2.
        // Find the bin which contains the 50% data point.
        //
        totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1];
        psTrace(__func__, 6, "Total data points is %d\n", totalDataPoints);
        psS32 binMedian;
        if (totalDataPoints/2.0 < cumulativeRobustHistogram->nums->data.F32[0]) {
            // Special case: the median is in the first bin.  Perhaps we should recode this.
            // XXX: Must try a special case where the median in in the last bin.
            binMedian = 0;
        } else {
            tmpScalar.data.F32 = totalDataPoints/2.0;
            binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
            if (binMedian < 0) {
                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian);
                psFree(tmpStatsMinMax);
                psFree(robustHistogram);
                psFree(cumulativeRobustHistogram);
                psFree(tmpMaskVec);
                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
                return(1);
            }
        }
        psTrace(__func__, 6, "The median bin is %d (%.2f to %.2f)\n", binMedian, cumulativeRobustHistogram->bounds->data.F32[binMedian], cumulativeRobustHistogram->bounds->data.F32[binMedian+1]);

        //
        // ADD: Step 3.
        // Interpolate to the exact 50% position: this is the robust histogram median.
        //
        stats->robustMedian = fitQuadraticSearchForYThenReturnX(
                                  *(psVector* *)&cumulativeRobustHistogram->bounds,
                                  *(psVector* *)&cumulativeRobustHistogram->nums,
                                  binMedian,
                                  totalDataPoints/2.0);
        if (isnan(stats->robustMedian)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n");
            psFree(tmpStatsMinMax);
            psFree(robustHistogram);
            psFree(cumulativeRobustHistogram);
            psFree(tmpMaskVec);
            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
            return(1);
        }
        psTrace(__func__, 6, "Current robust median is %f\n", stats->robustMedian);

        //
        // ADD: Step 4.
        // Find the bins which contains the 15.8655% and 84.1345% data points.
        //
        psS32 binLo = 0;
        psS32 binHi = 0;
        tmpScalar.data.F32 = totalDataPoints * 0.158655f;
        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
            binLo = 0;
        } else {
            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
            binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
            if (binLo > cumulativeRobustHistogram->nums->n-1) {
                binLo = cumulativeRobustHistogram->nums->n-1;
            }
        }
        tmpScalar.data.F32 = totalDataPoints * 0.841345f;
        if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
            binHi = 0;
        } else {
            // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined.
            binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar);
            if (binHi > cumulativeRobustHistogram->nums->n-1) {
                binHi = cumulativeRobustHistogram->nums->n-1;
            }
        }
        psTrace(__func__, 6, "The 15.8655-percent and 84.1345-percent data point bins are (%d, %d).\n",
                binLo, binHi);
        psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binLo));
        psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binHi));

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

        //
        // ADD: Step 4b.
        // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions.
        //
        // XXX: I am overriding the ADD for now.  The following code implements
        // the ADD exactly and is having problems fitting a 2nd-order polynomial
        // to data y-values suchs as (1, 1, 100).  Therefore, I commented out the
        // code and am doing simply linear interpolation.
        //
        psF32 binLoF32;
        psF32 binHiF32;
        if (0) {
            binLoF32 = fitQuadraticSearchForYThenReturnX(
                           *(psVector* *)&cumulativeRobustHistogram->bounds,
                           *(psVector* *)&cumulativeRobustHistogram->nums,
                           binLo,
                           totalDataPoints * 0.158655f);
            binHiF32 = fitQuadraticSearchForYThenReturnX(
                           *(psVector* *)&cumulativeRobustHistogram->bounds,
                           *(psVector* *)&cumulativeRobustHistogram->nums,
                           binHi,
                           totalDataPoints * 0.841345f);
            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n",
                    binLoF32, binHiF32);
        }

        //
        // This code basically interpolates to find the positions exactly.
        //
        if (1) {
            psTrace(__func__, 6, "binLo is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
                    binLo, cumulativeRobustHistogram->nums->data.F32[binLo],
                    cumulativeRobustHistogram->nums->data.F32[binLo+1]);
            psTrace(__func__, 6, "binHi is %d.  Nums at that bin and the next are (%.2f, %.2f)\n",
                    binHi, cumulativeRobustHistogram->nums->data.F32[binHi],
                    cumulativeRobustHistogram->nums->data.F32[binHi+1]);

            psF32 deltaNums;
            psF32 deltaBounds;
            psF32 prevPixels;
            psF32 percentNums;
            psF32 base;
            deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binLo+1] - cumulativeRobustHistogram->bounds->data.F32[binLo];
            if (binLo == 0) {
                deltaNums = cumulativeRobustHistogram->nums->data.F32[0];
                prevPixels = 0;
            } else {
                deltaNums = cumulativeRobustHistogram->nums->data.F32[binLo] - cumulativeRobustHistogram->nums->data.F32[binLo-1];
                prevPixels = cumulativeRobustHistogram->nums->data.F32[binLo-1];
            }
            percentNums = (totalDataPoints * 0.158655f) - prevPixels;
            base = cumulativeRobustHistogram->bounds->data.F32[binLo];
            binLoF32 = base + (deltaBounds / deltaNums) * percentNums;
            psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                    base, deltaBounds, deltaNums, prevPixels, percentNums);

            deltaBounds = cumulativeRobustHistogram->bounds->data.F32[binHi+1] - cumulativeRobustHistogram->bounds->data.F32[binHi];
            if (binHi == 0) {
                deltaNums = cumulativeRobustHistogram->nums->data.F32[0];
                prevPixels = 0;
            } else {
                deltaNums = cumulativeRobustHistogram->nums->data.F32[binHi] - cumulativeRobustHistogram->nums->data.F32[binHi-1];
                prevPixels = cumulativeRobustHistogram->nums->data.F32[binHi-1];
            }
            percentNums = (totalDataPoints * 0.841345f) - prevPixels;
            base = cumulativeRobustHistogram->bounds->data.F32[binHi];
            binHiF32 = base + (deltaBounds / deltaNums) * percentNums;
            psTrace(__func__, 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n",
                    base, deltaBounds, deltaNums, prevPixels, percentNums);
            psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32);
        }

        //
        // ADD: Step 5.
        // Determine SIGMA as 1/2 of the distance between these positions.
        //
        sigma = (binHiF32 - binLoF32) / 2.0;
        psTrace(__func__, 6, "The current sigma is %f.\n", sigma);
        stats->robustStdev = sigma;

        //
        // ADD: Step 6.
        // If the measured SIGMA is less than 2 times the bin size, exclude
        // points which are more than 25 bins from the median,
        // recalculate the bin size, and perform the algorithm again.
        //
        if (sigma < (2 * binSize)) {
            psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
            psS32 maskLo = PS_MAX(0, (binMedian - 25));
            psS32 maskHi = PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25));
            psF32 medianLo = robustHistogram->bounds->data.F32[maskLo];
            psF32 medianHi = robustHistogram->bounds->data.F32[maskHi];
            psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian);
            psTrace(__func__, 6, "The median is at bin number %d.  We mask bins outside the bin range (%d:%d)\n",
                    binMedian, maskLo, maskHi);
            psTrace(__func__, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
            for (psS32 i = 0 ; i < myVector->n ; i++) {
                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
                    tmpMaskVec->data.U8[i] = 1;
                    psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]);
                }
            }
            psFree(robustHistogram);
            psFree(cumulativeRobustHistogram);
            iterate++;
        } else {
            psTrace(__func__, 6, "*************: No more iteration.  sigma is %f\n", sigma);
            iterate = 0;
        }
    }

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

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

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

    psS32 N50 = 0;
    for (psS32 i = 0 ; i < myVector->n ; i++) {
        // XXX: use maskVal here?
        if ((0 == tmpMaskVec->data.U8[i]) &&
                (binLo25F32 <= myVector->data.F32[i]) &&
                (binHi25F32 >= myVector->data.F32[i])) {
            N50++;
        }
    }
    stats->robustN50 = N50;
    psTrace(__func__, 6, "The robustN50 is %d.\n", N50);

    // ************************************************************************
    if ((stats->options & PS_STAT_FITTED_MEAN) ||
            (stats->options & PS_STAT_FITTED_STDEV)) {
        //
        // New algorithm for calculated mean and stdev.
        //

        // XXX: Where is this documented?
        psF32 dN = (psF32) (0.17 * N50);
        if (dN < 1.0) {
            dN = 1.0;
        } else if (dN > 4.0) {
            dN = 4.0;
        }
        psF32 newBinSize = sigma / dN;

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

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

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

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

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

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

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

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

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

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

        for (psS32 i = binMin ; i <= binMax ; i++) {
            y->data.F32[j] = newHistogramSmoothed->data.F32[i];
            x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
            ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
            j++;
        }
        if (psTraceGetLevel(__func__) >= 8) {
            // XXX: Print the x array somehow.
            PS_VECTOR_PRINT_F32(y);
        }

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

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

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

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

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

        psFree(newHistogramSmoothed);
        psFree(newHistogram);
        psFree(x);
        psFree(y);
        psFree(min);
        psFree(params);
    }

    psFree(tmpStatsMinMax);
    psFree(cumulativeRobustHistogram);
    psFree(tmpMaskVec);
    psFree(robustHistogram);

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

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

/* FUNCTION IMPLEMENTATION - PUBLIC                                          */

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

static void histogramFree(psHistogram* myHist);

static void statsFree(psStats *newStruct)
{
    // There are non dynamic allocated items
}

/******************************************************************************
    psStatsAlloc(): This routine must create a new psStats data structure.
 *****************************************************************************/
psStats* psStatsAlloc(psStatsOptions options)
{
    psTrace(__func__, 3,"---- %s() begin  ----\n", __func__);
    psStats* newStruct = NULL;

    newStruct = (psStats* ) psAlloc(sizeof(psStats));
    psMemSetDeallocator(newStruct, (psFreeFunc)statsFree);
    newStruct->sampleMean = NAN;
    newStruct->sampleMedian = NAN;
    newStruct->sampleStdev = NAN;
    newStruct->sampleUQ = NAN;
    newStruct->sampleLQ = NAN;
    newStruct->robustMedian = NAN;
    newStruct->robustStdev = NAN;
    newStruct->robustUQ = NAN;
    newStruct->robustLQ = NAN;
    newStruct->robustN50 = -1;
    newStruct->fittedMean = NAN;
    newStruct->fittedStdev = NAN;
    newStruct->fittedNfit = -1;
    newStruct->clippedMean = NAN;
    newStruct->clippedStdev = NAN;
    newStruct->clippedNvalues = -1;     // XXX: This is never used
    newStruct->clipSigma = 3.0;
    newStruct->clipIter = 3;
    newStruct->min = NAN;
    newStruct->max = NAN;
    newStruct->binsize = NAN;
    newStruct->nSubsample = 100000;
    newStruct->options = options;

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


bool psMemCheckStats(psPtr ptr)
{
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)statsFree );
}



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

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

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

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

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

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

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

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

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

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

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

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

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

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


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



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

    PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0], out->bounds->data.F32[(out->bounds->n)-1], -4);

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

/******************************************************************************
psVectorStats(myVector, maskVector, maskVal, stats): this is the public API
function which calls the above private stats functions based on what bits
were set in stats->options.
 
Inputs
    myVector
    maskVector
    maskVal
    stats
Returns
    The stats structure.
 
XXX: Should we free stats if the asserts fail?
 *****************************************************************************/
psStats* psVectorStats(psStats* stats,
                       const psVector* in,
                       const psVector* errors,
                       const psVector* mask,
                       psMaskType maskVal)
{
    psTrace(__func__, 3,"---- %s() begin  ----\n", __func__);
    PS_ASSERT_PTR_NON_NULL(stats, NULL);
    PS_ASSERT_VECTOR_NON_NULL(in, NULL);
    if (mask != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, stats);
        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, stats);
    }
    if (errors != NULL) {
        PS_ASSERT_VECTORS_SIZE_EQUAL(errors, in, stats);
        PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats);
    }
    // XXX: Assert that "in" is F64, F32, U16, or S8

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

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

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

    if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) {
        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, stats);
    }

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

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
        if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian() returned an error.\n");
            stats->sampleMedian = NAN;
        }
    }

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

    // ************************************************************************
    if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
        if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) {
            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleQuartiles() returned an error.\n");
            stats->sampleLQ = NAN;
            stats->sampleUQ = NAN;
        }
    }

    // ************************************************************************
    if ((stats->options & PS_STAT_ROBUST_MEDIAN) ||
            (stats->options & PS_STAT_ROBUST_STDEV) ||
            (stats->options & PS_STAT_ROBUST_QUARTILE) ||
            (stats->options & PS_STAT_FITTED_MEAN) ||
            (stats->options & PS_STAT_FITTED_STDEV)) {
        if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED);
            psFree(stats);
            psTrace(__func__, 3,"---- %s(NULL) end  ----\n", __func__);
            return(NULL);
        }
    }

    // ************************************************************************
    if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
        psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats);
        if (rc < 0) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n");
            stats->clippedMean = NAN;
            stats->clippedStdev = NAN;
        }
    }

    // ************************************************************************
    if (stats->options & PS_STAT_MAX) {
        if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector maximum");
            stats->max = NAN;
        }
    }

    // ************************************************************************
    if (stats->options & PS_STAT_MIN) {
        if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) {
            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum");
            stats->min = NAN;
        }
    }

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

