IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 7, 2009, 4:08:25 PM (17 years ago)
Author:
Paul Price
Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:
branches/pap
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/pap

  • branches/pap/psLib/src/math/psStats.c

    r23416 r25027  
    2424// unity, even though the standard deviation is not defined in that case (NAN).
    2525
     26// reworking the return values and failure conditions: it should only be an error if the
     27// inputs imply a programming error: eg, NULL data vectors, non-sensical input
     28// parameters.  If the statistic cannot be calculated (0 length vector, 0 range, no
     29// unmasked data values), the statistic should be reported as NAN, but an error should not
     30// be raised.  (TBD: do we need to have a unique field in psStats or a return parameter
     31// that can be checked for an invalid result?)
     32
    2633#ifdef HAVE_CONFIG_H
    27 # include "config.h"
     34#include "config.h"
    2835#endif
    2936
    … …  
    4249#include "psMemory.h"
    4350#include "psAbort.h"
    44 //#include "psImage.h"
    4551#include "psVector.h"
    4652#include "psTrace.h"
    … …  
    6975#define PS_POLY_MEDIAN_MAX_ITERATIONS 30
    7076
     77#define TRACE "psLib.math"
     78
    7179#define MASK_MARK 0x80   // XXX : can we change this? bit to use internally to mark data as bad
    7280#define PS_ROBUST_MAX_ITERATIONS 20     // Maximum number of iterations for robust statistics
    … …  
    8795            break; \
    8896          case PS_BINARY_DISECT_OUTSIDE_RANGE: \
    89             psTrace ("psLib.math", 6, "selected bin outside range"); \
     97            psTrace(TRACE, 6, "selected bin outside range"); \
    9098            if (USE_END == -1) { RESULT = 0; } \
    9199            if (USE_END == +1) { RESULT = VECTOR->n - 1; } \
    … …  
    116124        } \
    117125        Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \
    118         psTrace("psLib.math", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
     126        psTrace(TRACE, 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \
    119127                Xo, Yo, dX, dY, Xt, VALUE); \
    120128        RESULT = Xt; }
    121129
    122 #define TRACE "psLib.math"
     130# define COUNT_WARNING(LIMIT, INTERVAL, ...) { \
     131        static int nCalls = 1; \
     132        if (nCalls < LIMIT) { \
     133            psWarning(__VA_ARGS__); \
     134        } \
     135        if (!(nCalls % INTERVAL)) { \
     136            psWarning(__VA_ARGS__); \
     137            psWarning("(warning raised %d times)", nCalls); \
     138        } \
     139        nCalls ++; \
     140}
     141
    123142
    124143/*****************************************************************************/
    … …  
    156175static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
    157176// static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    158 static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
     177// static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
     178static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    159179
    160180/******************************************************************************
    … …  
    182202                                 psStats* stats)
    183203{
    184     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    185 
    186204    long count = 0;                     // Number of points contributing to this mean
    187205    psF32 mean = 0.0;                   // The mean
    … …  
    199217    for (long i = 0; i < numData; i++) {
    200218        // Check if the data is with the specified range
    201         if (!isfinite(data[i]))
    202             continue;
     219        if (!isfinite(data[i]))
     220            continue;
    203221        if (useRange && (data[i] < stats->min))
    204222            continue;
    … …  
    222240    }
    223241    stats->sampleMean = mean;
    224     if (isnan(mean)) {
    225         // XXX raise an error here?
    226         psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__);
    227         return true;
    228     }
    229 
    230     stats->results |= PS_STAT_SAMPLE_MEAN;
    231     psTrace(TRACE, 4, "---- %s(true) end ----\n", __func__);
     242
     243    if (!isnan(mean)) {
     244        stats->results |= PS_STAT_SAMPLE_MEAN;
     245    }
    232246    return true;
    233247}
    … …  
    252266        )
    253267{
    254     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    255268    psF32 max = -PS_MAX_F32;            // The calculated maximum
    256269    psF32 min = PS_MAX_F32;             // The calculated minimum
    … …  
    265278    for (long i = 0; i < num; i++) {
    266279        // Check if the data is with the specified range
    267         if (!isfinite(vector[i]))
    268             continue;
     280        if (!isfinite(vector[i]))
     281            continue;
    269282        if (useRange && (vector[i] < stats->min))
    270283            continue;
    … …  
    289302        stats->results |= PS_STAT_MAX;
    290303    }
    291     psTrace(TRACE, 4, "---- %s(%d) end ----\n", __func__, numValid);
    292304    return numValid;
    293305}
    … …  
    303315                               psStats* stats)
    304316{
    305     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    306 
    307317    bool useRange = stats->options & PS_STAT_USE_RANGE;
    308318    psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference the vector
    … …  
    319329    // into the temporary vectors.
    320330    for (long i = 0; i < inVector->n; i++) {
    321         if (!isfinite(input[i]))
    322             continue;
     331        if (!isfinite(input[i]))
     332            continue;
    323333        if (useRange && (input[i] < stats->min))
    324334            continue;
    … …  
    334344
    335345    if (count == 0) {
    336         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n");
     346        COUNT_WARNING(10, 100, "No valid data in input vector.\n");
    337347        stats->sampleUQ = NAN;
    338348        stats->sampleLQ = NAN;
    … …  
    343353    // Sort the temporary vector.
    344354    if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
     355        // an error in psVectorSort is a serious error
    345356        psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
    346357        stats->sampleUQ = NAN;
    … …  
    364375    stats->results |= PS_STAT_SAMPLE_QUARTILE;
    365376
    366     // Return "true" on success.
    367     psTrace(TRACE, 4, "---- %s(true) end ----\n", __func__);
    368377    return true;
    369378}
    … …  
    390399                              psStats* stats)
    391400{
    392     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    393 
    394401    // This procedure requires the mean.  If it has not been already
    395402    // calculated, then call vectorSampleMean()
    … …  
    400407    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    401408    if (isnan(stats->sampleMean)) {
    402         psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): sample mean is NAN.  "
    403                 "Setting stats->sampleStdev = NAN.\n");
     409        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.");
    404410        stats->sampleStdev = NAN;
    405411        return true;
    … …  
    419425    for (long i = 0; i < myVector->n; i++) {
    420426        // Check if the data is with the specified range
    421         if (!isfinite(data[i]))
    422             continue;
     427        if (!isfinite(data[i]))
     428            continue;
    423429        if (useRange && (data[i] < stats->min)) {
    424430            continue;
    … …  
    445451        // Assume that the user knows what he's doing when he masks out everything --> no error.
    446452        stats->sampleStdev = NAN;
    447         psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld).  "
    448                 "Setting stats->sampleStdev = NAN.\n", count);
     453        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). Setting stats->sampleStdev = NAN.\n", count);
    449454        return true;
    450455    }
    451456    if (count == 1) {
    452457        stats->sampleStdev = 0.0;
    453         psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld).  "
    454                 "Setting stats->sampleStdev = 0.0.\n", count);
     458        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). Setting stats->sampleStdev = 0.0.\n", count);
    455459        return true;
    456460    }
    … …  
    462466    }
    463467    stats->results |= PS_STAT_SAMPLE_STDEV;
    464 
    465     psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    466468
    467469    return true;
    … …  
    473475                                psStats* stats)
    474476{
    475     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    476 
    477477    // This procedure requires the mean and standard deviation
    478478    if (!(stats->results & PS_STAT_SAMPLE_MEAN)) {
    … …  
    480480    }
    481481    if (isnan(stats->sampleMean)) {
    482         psTrace(TRACE, 5, "WARNING: vectorSampleMoments(): sample mean is NAN.\n");
     482        COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample mean is NAN.\n");
    483483        goto SAMPLE_MOMENTS_BAD;
    484484    }
    … …  
    487487    }
    488488    if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) {
    489         psTrace(TRACE, 5, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n");
     489        COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n");
    490490        goto SAMPLE_MOMENTS_BAD;
    491491    }
    … …  
    502502    for (long i = 0; i < myVector->n; i++) {
    503503        // Check if the data is with the specified range
    504         if (!isfinite(data[i]))
    505             continue;
     504        if (!isfinite(data[i]))
     505            continue;
    506506        if (useRange && (data[i] < stats->min)) {
    507507            continue;
    … …  
    534534    stats->results |= PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS;
    535535
    536     psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    537 
    538536    return true;
    539537
    … …  
    542540    stats->sampleSkewness = NAN;
    543541    stats->sampleKurtosis = NAN;
    544     stats->results |= PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS;
    545542    return true;
    546543}
    … …  
    565562    )
    566563{
    567     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    568     psTrace(TRACE, 4, "Trace level is %d\n", psTraceGetLevel("psLib.math"));
    569 
    570564    // Ensure that stats->clipIter is within the proper range.
    571565    PS_ASSERT_INT_WITHIN_RANGE(stats->clipIter,
    … …  
    601595    vectorSampleMedian(myVector, tmpMask, maskVal, stats);
    602596    if (isnan(stats->sampleMedian)) {
    603         psTrace(TRACE, 5, "Call to vectorSampleMedian returned NAN\n");
    604         psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__);
    605         return false;
     597        COUNT_WARNING(10, 100, "Call to vectorSampleMedian returned NAN\n");
     598        return true;
    606599    }
    607600    psTrace(TRACE, 6, "The initial sample median is %f\n", stats->sampleMedian);
    … …  
    610603    vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
    611604    if (isnan(stats->sampleStdev)) {
    612         psTrace(TRACE, 5, "Call to vectorSampleStdev returned NAN\n");
    613         psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__);
    614         return false;
     605        COUNT_WARNING(10, 100, "Call to vectorSampleStdev returned NAN\n");
     606        return true;
    615607    }
    616608    psTrace(TRACE, 6, "The initial sample stdev is %f\n", stats->sampleStdev);
    … …  
    669661        if (isnan(stats->sampleMean) || isnan(stats->sampleStdev)) {
    670662            iter = stats->clipIter;
    671             psError(PS_ERR_UNKNOWN, true, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
    672             return false;
     663            COUNT_WARNING(10, 100, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
     664            clippedMean = NAN;
     665            clippedStdev = NAN;
     666            return true;
    673667        } else {
    674668            clippedMean = stats->sampleMean;
    … …  
    693687    psTrace(TRACE, 6, "The final clipped stdev is %f\n", clippedStdev);
    694688
    695     psTrace(TRACE, 4, "---- %s(true) end ----\n", __func__);
    696689    return true;
    697690}
    … …  
    722715                              psStats* stats)
    723716{
    724     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    725717    if (psTraceGetLevel("psLib.math") >= 8) {
    726718        PS_VECTOR_PRINT_F32(myVector);
    … …  
    767759        if (numValid == 0 || isnan(min) || isnan(max)) {
    768760            // Data range calculation failed
    769             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     761            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    770762            goto escape;
    771763        }
    … …  
    779771            stats->robustLQ = min;
    780772            stats->robustN50 = numValid;
     773            // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits:
    781774            stats->results |= PS_STAT_ROBUST_MEDIAN;
    782775            stats->results |= PS_STAT_ROBUST_STDEV;
    783776            stats->results |= PS_STAT_ROBUST_QUARTILE;
    784             psTrace(TRACE, 5, "All data points have the same value: %f.\n", min);
    785             psTrace(TRACE, 4, "---- %s(0) end  ----\n", __func__);
     777            COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min);
    786778            psFree(mask);
    787779            psFree(statsMinMax);
    … …  
    809801        // XXXXX we need to consider this step if errors -> variance
    810802        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
     803            // if psVectorHistogram returns false, we have a programming error
    811804            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n");
    812             goto escape;
     805            psFree(histogram);
     806            psFree(cumulative);
     807            psFree(statsMinMax);
     808            psFree(mask);
     809
     810            return false;
    813811        }
    814812        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    837835        PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);
    838836
    839         psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian,
    840                 cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
    841 
    842         // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median.
    843         stats->robustMedian = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     837        psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
     838
     839        // ADD step 3: Interpolate to the exact 50% position in bin units
     840        stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     841        // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     842        // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
     843
     844        // convert bin to bin value: this is the robust histogram median.
    844845        if (isnan(stats->robustMedian)) {
    845             psError(PS_ERR_UNKNOWN, false,
    846                     "Failed to fit a quadratic and calculate the 50-percent position.\n");
     846            COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n");
    847847            goto escape;
    848848        }
    … …  
    866866
    867867        if ((binLo < 0) || (binHi < 0)) {
    868             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
     868            COUNT_WARNING(10, 100, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
    869869            goto escape;
    870870        }
    … …  
    962962                stats->results |= PS_STAT_ROBUST_STDEV;
    963963                stats->results |= PS_STAT_ROBUST_QUARTILE;
    964                 psTrace(TRACE, 5, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
    965                 psTrace(TRACE, 4, "---- %s(0) end  ----\n", __func__);
     964                COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
    966965                psFree(mask);
    967966                psFree(statsMinMax);
    968                 return false;
     967                return true;
    969968            }
    970969        } else {
    … …  
    988987    // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
    989988    // positions.
    990     psF32 binLo25F32 = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
    991     psF32 binHi25F32 = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
     989    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
     990    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
    992991    if (isnan(binLo25F32) || isnan(binHi25F32)) {
    993         psError(PS_ERR_UNKNOWN, false,
    994                 "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
     992        COUNT_WARNING(10, 100, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n");
    995993        goto escape;
    996994    }
    … …  
    10201018    stats->results |= PS_STAT_ROBUST_QUARTILE;
    10211019
    1022     psTrace(TRACE, 4, "---- %s(0) end  ----\n", __func__);
    10231020    return true;
    10241021
    … …  
    10331030    stats->results |= PS_STAT_ROBUST_QUARTILE;
    10341031
    1035     psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1036 
    10371032    psFree(histogram);
    10381033    psFree(cumulative);
    … …  
    10401035    psFree(mask);
    10411036
    1042     return false;
     1037    return true;
    10431038}
    10441039
    … …  
    10571052    // calculated, then call vectorSampleMean()
    10581053    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1059         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1054        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1055            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1056            return false;
     1057        }
    10601058    }
    10611059
    … …  
    10641062        stats->fittedStdev = NAN;
    10651063        stats->fittedStdev = NAN;
    1066         psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1067         return false;
     1064        return true;
    10681065    }
    10691066
    … …  
    10961093        float max = statsMinMax->max;
    10971094        if (numValid == 0 || isnan(min) || isnan(max)) {
    1098             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1095            psTrace(TRACE, 5, "Failed to calculate the min/max of the input vector.\n");
    10991096            psFree(statsMinMax);
    1100             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1101             return false;
     1097            return true;
    11021098        }
    11031099
    … …  
    11111107        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    11121108        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
     1109            // if psVectorHistogram returns false, we have a programming error
    11131110            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n");
    11141111            psFree(histogram);
    … …  
    11751172            PS_VECTOR_PRINT_F32(y);
    11761173        }
     1174
     1175        // psMinimizeLMChi2 can return false for bad data as well as for serious failures
    11771176        if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
    11781177            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    … …  
    11831182            psFree(histogram);
    11841183            psFree(statsMinMax);
    1185             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1186             return false;
     1184            return true;
    11871185        }
    11881186        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    12351233    // calculated, then call vectorSampleMean()
    12361234    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1237         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1235        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1236            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1237            return false;
     1238        }
    12381239    }
    12391240
    … …  
    12431244        stats->fittedStdev = NAN;
    12441245        psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1245         return false;
     1246        return true;
    12461247    }
    12471248
    … …  
    12741275        float max = statsMinMax->max;
    12751276        if (numValid == 0 || isnan(min) || isnan(max)) {
    1276             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1277            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    12771278            psFree(statsMinMax);
    12781279            psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1279             return false;
     1280            return true;
    12801281        }
    12811282
    … …  
    13541355        // psVectorInit (fitMask, 0);
    13551356
     1357        // XXX not sure if these should result in errors or not...
    13561358        if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
    13571359            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    … …  
    13661368
    13671369        if (poly->coeff[2] >= 0.0) {
    1368             psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n",
    1369                     poly->coeff[0], poly->coeff[1], poly->coeff[2]);
     1370            psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    13701371            psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
    13711372            psFree(x);
    … …  
    14291430    // calculated, then call vectorSampleMean()
    14301431    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1431         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1432        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1433            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1434            return false;
     1435        }
    14321436    }
    14331437
    … …  
    14371441        stats->fittedStdev = NAN;
    14381442        psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1439         return false;
     1443        return true;
    14401444    }
    14411445
    … …  
    14681472        float max = statsMinMax->max;
    14691473        if (numValid == 0 || isnan(min) || isnan(max)) {
    1470             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1474            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    14711475            psFree(statsMinMax);
    14721476            psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1473             return false;
     1477            return true;
    14741478        }
    14751479
    … …  
    15161520        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    15171521        if (binMin == binMax) {
    1518             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1522            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    15191523            psFree(statsMinMax);
    1520             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1521             return false;
     1524            return true;
    15221525        }
    15231526
    … …  
    17241727    // calculated, then call vectorSampleMean()
    17251728    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1726         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1729        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1730            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1731            return false;
     1732        }
    17271733    }
    17281734
    17291735    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    1730     if (isnan(stats->robustMedian)) {
    1731         stats->fittedStdev = NAN;
    1732         stats->fittedStdev = NAN;
    1733         psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1734         return false;
    1735     }
     1736    if (isnan(stats->robustMedian)) goto escape;
    17361737
    17371738    float guessStdev = stats->robustStdev;  // pass the guess sigma
    … …  
    17631764        float max = statsMinMax->max;
    17641765        if (numValid == 0 || isnan(min) || isnan(max)) {
    1765             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1766            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    17661767            psFree(statsMinMax);
    1767             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1768             return false;
     1768            goto escape;
     1769        }
     1770
     1771        // If all data points have the same value, then we set the appropriate members of stats and return.
     1772        if (fabs(max - min) <= FLT_EPSILON) {
     1773            COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min);
     1774            stats->fittedMean = min;
     1775            stats->fittedStdev = 0.0;
     1776            stats->results |= PS_STAT_FITTED_MEAN_V4;
     1777            stats->results |= PS_STAT_FITTED_STDEV_V4;
     1778            return true;
    17691779        }
    17701780
    … …  
    17801790        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    17811791        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    1782             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics v4.\n");
     1792            COUNT_WARNING(10, 100, "Unable to generate histogram for fitted statistics v4.\n");
    17831793            psFree(histogram);
    17841794            psFree(statsMinMax);
    1785             return false;
     1795            goto escape;
    17861796        }
    17871797        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    18131823        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    18141824        if (binMin == binMax) {
    1815             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1825            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    18161826            psFree(statsMinMax);
    1817             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1818             return false;
     1827            goto escape;
    18191828        }
    18201829
    … …  
    18771886
    18781887            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
     1888            // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
     1889            // we probably should be able to handle the data errors gracefully
    18791890            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    18801891            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
    … …  
    18831894
    18841895            if (!status) {
    1885                 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
     1896                COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
    18861897                psFree(poly);
    18871898                psFree(histogram);
    18881899                psFree(statsMinMax);
    1889                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1890                 return false;
     1900                goto escape;
    18911901            }
    18921902
    18931903            if (poly->coeff[2] >= 0.0) {
    1894                 psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n",
    1895                         poly->coeff[0], poly->coeff[1], poly->coeff[2]);
     1904                COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    18961905                psFree(poly);
    18971906                psFree(histogram);
    … …  
    19071916                }
    19081917
    1909                 psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
    1910                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1911                 return false;
     1918                COUNT_WARNING(10, 100, "fit did not converge\n");
     1919                goto escape;
    19121920            }
    19131921
    … …  
    19751983
    19761984            if (!status) {
    1977                 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
     1985                COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
    19781986                psFree(poly);
    19791987                psFree(histogram);
    19801988                psFree(statsMinMax);
    1981                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1982                 return false;
     1989                goto escape;
    19831990            }
    19841991
    … …  
    20252032    psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
    20262033
     2034    stats->results |= PS_STAT_FITTED_MEAN_V4;
     2035    stats->results |= PS_STAT_FITTED_STDEV_V4;
     2036
     2037    return true;
     2038
     2039escape:
     2040    stats->fittedMean = NAN;
     2041    stats->fittedStdev = NAN;
    20272042    stats->results |= PS_STAT_FITTED_MEAN_V4;
    20282043    stats->results |= PS_STAT_FITTED_STDEV_V4;
    … …  
    21552170psStats* p_psStatsAlloc(const char *file, unsigned int lineno, const char *func, psStatsOptions options)
    21562171{
    2157     psTrace(TRACE, 3,"---- %s() begin  ----\n", __func__);
    2158 
    21592172    psStats *stats = p_psAlloc(file, lineno, func, sizeof(psStats));
    21602173    psMemSetDeallocator(stats, (psFreeFunc)statsFree);
    … …  
    21722185    stats->tmpMask = NULL;
    21732186
    2174     psTrace(TRACE, 3, "---- %s() end  ----\n", __func__);
    21752187    return stats;
    21762188}
    … …  
    22312243                   psVectorMaskType maskVal)
    22322244{
    2233     psTrace(TRACE, 3,"---- %s() begin  ----\n", __func__);
    22342245    PS_ASSERT_PTR_NON_NULL(stats, false);
    22352246    PS_ASSERT_VECTOR_NON_NULL(in, false);
    … …  
    23802391    psFree(errorsF32);
    23812392    psFree(maskVector);
    2382     psTrace(TRACE, 3,"---- %s() end  ----\n", __func__);
    23832393    return status;
    23842394}
    … …  
    27452755    return tmpFloat;
    27462756}
    2747 # endif
    27482757
    27492758/******************************************************************************
    … …  
    28732882    return tmpFloat;
    28742883}
     2884# endif
     2885
     2886/******************************************************************************
     2887fitQuadraticSearchForYThenReturnXusingValues(*xVec, *yVec, binNum, yVal): A general routine
     2888which fits a quadratic to three points and returns the x bin value corresponding to the input
     2889y-value.  This routine takes psVectors of x/y pairs as input, and fits a quadratic to the 3
     2890points surrounding element binNum in the vectors.  This version uses the values of x[i] for the
     2891x coordinates (not the midpoints).  This is appropriate for a cumulative histogram.  It then
     2892determines for what value x does that quadratic f(x) = yVal (the input parameter).
     2893
     2894XXX this function is used a fair amount in an inner loop: the polynomial fitting and evaluation
     2895could easily be done with statically allocated doubles, skipping the psLib versions of
     2896polynomial fitting, etc.
     2897
     2898*****************************************************************************/
     2899static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec,
     2900                                                 psVector *yVec,
     2901                                                 psS32 binNum,
     2902                                                 psF32 yVal
     2903    )
     2904{
     2905    psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
     2906    if (psTraceGetLevel("psLib.math") >= 8) {
     2907        PS_VECTOR_PRINT_F32(xVec);
     2908        PS_VECTOR_PRINT_F32(yVec);
     2909    }
     2910
     2911    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
     2912    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
     2913    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
     2914    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
     2915    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
     2916    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
     2917
     2918    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     2919    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     2920    psF32 tmpFloat = 0.0f;
     2921
     2922    if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
     2923        // The general case.  We have all three points.
     2924        x->data.F64[0] = binNum - 1;
     2925        x->data.F64[1] = binNum;
     2926        x->data.F64[2] = binNum + 1;
     2927        y->data.F64[0] = yVec->data.F32[binNum - 1];
     2928        y->data.F64[1] = yVec->data.F32[binNum];
     2929        y->data.F64[2] = yVec->data.F32[binNum + 1];
     2930        psTrace(TRACE, 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]);
     2931        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
     2932        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
     2933
     2934        // Ensure that the y value lies within range of the y values.
     2935        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
     2936               ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
     2937            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     2938                    _("Specified yVal, %g, is not within y-range, %g to %g."),
     2939                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
     2940            return NAN;
     2941        }
     2942
     2943        // Ensure that the y values are monotonic.
     2944        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
     2945            ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
     2946            psError(PS_ERR_UNKNOWN, true,
     2947                    "This routine must be called with monotonically increasing or decreasing data points.\n");
     2948            psFree(x);
     2949            psFree(y);
     2950            return NAN;
     2951        }
     2952
     2953        // Determine the coefficients of the polynomial.
     2954        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
     2955        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
     2956            psError(PS_ERR_UNEXPECTED_NULL, false,
     2957                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
     2958                      "Returning NAN."));
     2959            psFree(x);
     2960            psFree(y);
     2961            return NAN;
     2962        }
     2963
     2964        psTrace(TRACE, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
     2965        psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
     2966        psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
     2967        psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n",
     2968                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
     2969                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
     2970                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
     2971
     2972        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
     2973        float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
     2974        psFree(myPoly);
     2975
     2976        if (isnan(binValue)) {
     2977            psError(PS_ERR_UNEXPECTED_NULL,
     2978                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
     2979            psFree(x);
     2980            psFree(y);
     2981            return(NAN);
     2982        }
     2983
     2984        // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
     2985        assert (binValue >= binNum - 1);
     2986        assert (binValue <= binNum + 1);
     2987
     2988        int fitBin = binValue;
     2989        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
     2990        float dY = binValue - fitBin;
     2991        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
     2992    } else {
     2993        // These are special cases where the bin is at the beginning or end of the vector.
     2994        if (binNum == 0) {
     2995            // We have two points only at the beginning of the vectors x and y.
     2996            // X = (dX/dY)(Y - Yo) + Xo
     2997            float dX = xVec->data.F32[1] - xVec->data.F32[0];
     2998            float dY = yVec->data.F32[1] - yVec->data.F32[0];
     2999            if (dY == 0.0) {
     3000                tmpFloat = xVec->data.F32[0];
     3001            } else {
     3002                tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
     3003            }
     3004        } else if (binNum == (xVec->n - 1)) {
     3005            // We have two points only at the end of the vectors x and y.
     3006            // X = (dX/dY)(Y - Yo) + Xo
     3007            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
     3008            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];
     3009            if (dY == 0.0) {
     3010                tmpFloat = xVec->data.F32[binNum-1];
     3011            } else {
     3012                tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
     3013            }
     3014        }
     3015    }
     3016
     3017    psTrace(TRACE, 6, "FIT: return %f\n", tmpFloat);
     3018    psFree(x);
     3019    psFree(y);
     3020
     3021    return tmpFloat;
     3022}
    28753023
    28763024/******************************************************************************
Note: See TracChangeset for help on using the changeset viewer.