IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 26, 2009, 1:59:32 PM (17 years ago)
Author:
beaumont
Message:

merged with trunk

Location:
branches/cnb_branches/cnb_branch_20090301
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branches/cnb_branch_20090301

  • branches/cnb_branches/cnb_branch_20090301/psLib/src/math/psStats.c

    r23594 r24244  
    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; }
    121 
    122 #define TRACE "psLib.math"
    123129
    124130/*****************************************************************************/
    … …  
    156162static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords);
    157163// 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);
     164// static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
     165static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    159166
    160167/******************************************************************************
    … …  
    182189                                 psStats* stats)
    183190{
    184     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    185 
    186191    long count = 0;                     // Number of points contributing to this mean
    187192    psF32 mean = 0.0;                   // The mean
    … …  
    222227    }
    223228    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__);
     229
     230    if (!isnan(mean)) {
     231        stats->results |= PS_STAT_SAMPLE_MEAN;
     232    }
    232233    return true;
    233234}
    … …  
    252253        )
    253254{
    254     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    255255    psF32 max = -PS_MAX_F32;            // The calculated maximum
    256256    psF32 min = PS_MAX_F32;             // The calculated minimum
    … …  
    289289        stats->results |= PS_STAT_MAX;
    290290    }
    291     psTrace(TRACE, 4, "---- %s(%d) end ----\n", __func__, numValid);
    292291    return numValid;
    293292}
    … …  
    303302                               psStats* stats)
    304303{
    305     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    306 
    307304    bool useRange = stats->options & PS_STAT_USE_RANGE;
    308305    psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference the vector
    … …  
    334331
    335332    if (count == 0) {
    336         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n");
     333        psLogMsg(TRACE, PS_LOG_DETAIL, "No valid data in input vector.\n");
    337334        stats->sampleUQ = NAN;
    338335        stats->sampleLQ = NAN;
    … …  
    343340    // Sort the temporary vector.
    344341    if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
     342        // an error in psVectorSort is a serious error
    345343        psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
    346344        stats->sampleUQ = NAN;
    … …  
    364362    stats->results |= PS_STAT_SAMPLE_QUARTILE;
    365363
    366     // Return "true" on success.
    367     psTrace(TRACE, 4, "---- %s(true) end ----\n", __func__);
    368364    return true;
    369365}
    … …  
    390386                              psStats* stats)
    391387{
    392     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    393 
    394388    // This procedure requires the mean.  If it has not been already
    395389    // calculated, then call vectorSampleMean()
    … …  
    400394    // If the mean is NAN, then generate a warning and set the stdev to NAN.
    401395    if (isnan(stats->sampleMean)) {
    402         psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): sample mean is NAN.  "
    403                 "Setting stats->sampleStdev = NAN.\n");
     396        psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.\n");
    404397        stats->sampleStdev = NAN;
    405398        return true;
    … …  
    445438        // Assume that the user knows what he's doing when he masks out everything --> no error.
    446439        stats->sampleStdev = NAN;
    447         psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld).  "
    448                 "Setting stats->sampleStdev = NAN.\n", count);
     440        psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): no valid psVector elements (%ld). Setting stats->sampleStdev = NAN.\n", count);
    449441        return true;
    450442    }
    451443    if (count == 1) {
    452444        stats->sampleStdev = 0.0;
    453         psTrace(TRACE, 5, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld).  "
     445        psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld).  "
    454446                "Setting stats->sampleStdev = 0.0.\n", count);
    455447        return true;
    … …  
    462454    }
    463455    stats->results |= PS_STAT_SAMPLE_STDEV;
    464 
    465     psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    466456
    467457    return true;
    … …  
    473463                                psStats* stats)
    474464{
    475     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    476 
    477465    // This procedure requires the mean and standard deviation
    478466    if (!(stats->results & PS_STAT_SAMPLE_MEAN)) {
    … …  
    480468    }
    481469    if (isnan(stats->sampleMean)) {
    482         psTrace(TRACE, 5, "WARNING: vectorSampleMoments(): sample mean is NAN.\n");
     470        psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleMoments(): sample mean is NAN.\n");
    483471        goto SAMPLE_MOMENTS_BAD;
    484472    }
    … …  
    487475    }
    488476    if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) {
    489         psTrace(TRACE, 5, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n");
     477        psLogMsg(TRACE, PS_LOG_DETAIL, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n");
    490478        goto SAMPLE_MOMENTS_BAD;
    491479    }
    … …  
    534522    stats->results |= PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS;
    535523
    536     psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    537 
    538524    return true;
    539525
    … …  
    542528    stats->sampleSkewness = NAN;
    543529    stats->sampleKurtosis = NAN;
    544     stats->results |= PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS;
    545530    return true;
    546531}
    … …  
    565550    )
    566551{
    567     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    568     psTrace(TRACE, 4, "Trace level is %d\n", psTraceGetLevel("psLib.math"));
    569 
    570552    // Ensure that stats->clipIter is within the proper range.
    571553    PS_ASSERT_INT_WITHIN_RANGE(stats->clipIter,
    … …  
    601583    vectorSampleMedian(myVector, tmpMask, maskVal, stats);
    602584    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;
     585        psLogMsg(TRACE, PS_LOG_DETAIL, "Call to vectorSampleMedian returned NAN\n");
     586        return true;
    606587    }
    607588    psTrace(TRACE, 6, "The initial sample median is %f\n", stats->sampleMedian);
    … …  
    610591    vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
    611592    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;
     593        psLogMsg(TRACE, PS_LOG_DETAIL, "Call to vectorSampleStdev returned NAN\n");
     594        return true;
    615595    }
    616596    psTrace(TRACE, 6, "The initial sample stdev is %f\n", stats->sampleStdev);
    … …  
    669649        if (isnan(stats->sampleMean) || isnan(stats->sampleStdev)) {
    670650            iter = stats->clipIter;
    671             psError(PS_ERR_UNKNOWN, true, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
    672             return false;
     651            psLogMsg(TRACE, PS_LOG_DETAIL, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n");
     652            clippedMean = NAN;
     653            clippedStdev = NAN;
     654            return true;
    673655        } else {
    674656            clippedMean = stats->sampleMean;
    … …  
    693675    psTrace(TRACE, 6, "The final clipped stdev is %f\n", clippedStdev);
    694676
    695     psTrace(TRACE, 4, "---- %s(true) end ----\n", __func__);
    696677    return true;
    697678}
    … …  
    722703                              psStats* stats)
    723704{
    724     psTrace(TRACE, 4, "---- %s() begin ----\n", __func__);
    725705    if (psTraceGetLevel("psLib.math") >= 8) {
    726706        PS_VECTOR_PRINT_F32(myVector);
    … …  
    767747        if (numValid == 0 || isnan(min) || isnan(max)) {
    768748            // Data range calculation failed
    769             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     749            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
    770750            goto escape;
    771751        }
    … …  
    779759            stats->robustLQ = min;
    780760            stats->robustN50 = numValid;
     761            // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits:
    781762            stats->results |= PS_STAT_ROBUST_MEDIAN;
    782763            stats->results |= PS_STAT_ROBUST_STDEV;
    783764            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__);
     765            psLogMsg(TRACE, PS_LOG_DETAIL, "All data points have the same value: %f.\n", min);
    786766            psFree(mask);
    787767            psFree(statsMinMax);
    … …  
    809789        // XXXXX we need to consider this step if errors -> variance
    810790        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
     791            // if psVectorHistogram returns false, we have a programming error
    811792            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n");
    812             goto escape;
     793            psFree(histogram);
     794            psFree(cumulative);
     795            psFree(statsMinMax);
     796            psFree(mask);
     797
     798            return false;
    813799        }
    814800        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    837823        PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);
    838824
    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);
     825        psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]);
     826
     827        // ADD step 3: Interpolate to the exact 50% position in bin units
     828        stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     829        // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     830        // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
     831
     832        // convert bin to bin value: this is the robust histogram median.
    844833        if (isnan(stats->robustMedian)) {
    845             psError(PS_ERR_UNKNOWN, false,
    846                     "Failed to fit a quadratic and calculate the 50-percent position.\n");
     834            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a quadratic and calculate the 50-percent position.\n");
    847835            goto escape;
    848836        }
    … …  
    866854
    867855        if ((binLo < 0) || (binHi < 0)) {
    868             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
     856            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n");
    869857            goto escape;
    870858        }
    … …  
    962950                stats->results |= PS_STAT_ROBUST_STDEV;
    963951                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__);
     952                psLogMsg(TRACE, PS_LOG_DETAIL, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
    966953                psFree(mask);
    967954                psFree(statsMinMax);
    968                 return false;
     955                return true;
    969956            }
    970957        } else {
    … …  
    988975    // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
    989976    // 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);
     977    psF32 binLo25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
     978    psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
    992979    if (isnan(binLo25F32) || isnan(binHi25F32)) {
    993         psError(PS_ERR_UNKNOWN, false,
    994                 "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
     980        psLogMsg(TRACE, PS_LOG_DETAIL, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n");
    995981        goto escape;
    996982    }
    … …  
    10201006    stats->results |= PS_STAT_ROBUST_QUARTILE;
    10211007
    1022     psTrace(TRACE, 4, "---- %s(0) end  ----\n", __func__);
    10231008    return true;
    10241009
    … …  
    10331018    stats->results |= PS_STAT_ROBUST_QUARTILE;
    10341019
    1035     psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1036 
    10371020    psFree(histogram);
    10381021    psFree(cumulative);
    … …  
    10401023    psFree(mask);
    10411024
    1042     return false;
     1025    return true;
    10431026}
    10441027
    … …  
    10571040    // calculated, then call vectorSampleMean()
    10581041    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1059         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1042        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1043            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1044            return false;
     1045        }
    10601046    }
    10611047
    … …  
    10641050        stats->fittedStdev = NAN;
    10651051        stats->fittedStdev = NAN;
    1066         psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1067         return false;
     1052        return true;
    10681053    }
    10691054
    … …  
    10961081        float max = statsMinMax->max;
    10971082        if (numValid == 0 || isnan(min) || isnan(max)) {
    1098             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1083            psTrace(TRACE, 5, "Failed to calculate the min/max of the input vector.\n");
    10991084            psFree(statsMinMax);
    1100             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1101             return false;
     1085            return true;
    11021086        }
    11031087
    … …  
    11111095        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    11121096        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
     1097            // if psVectorHistogram returns false, we have a programming error
    11131098            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n");
    11141099            psFree(histogram);
    … …  
    11751160            PS_VECTOR_PRINT_F32(y);
    11761161        }
     1162       
     1163        // psMinimizeLMChi2 can return false for bad data as well as for serious failures
    11771164        if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
    11781165            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    … …  
    11831170            psFree(histogram);
    11841171            psFree(statsMinMax);
    1185             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1186             return false;
     1172            return true;
    11871173        }
    11881174        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    12351221    // calculated, then call vectorSampleMean()
    12361222    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1237         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1223        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1224            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1225            return false;
     1226        }
    12381227    }
    12391228
    … …  
    12431232        stats->fittedStdev = NAN;
    12441233        psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1245         return false;
     1234        return true;
    12461235    }
    12471236
    … …  
    12741263        float max = statsMinMax->max;
    12751264        if (numValid == 0 || isnan(min) || isnan(max)) {
    1276             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1265            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
    12771266            psFree(statsMinMax);
    12781267            psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1279             return false;
     1268            return true;
    12801269        }
    12811270
    … …  
    13541343        // psVectorInit (fitMask, 0);
    13551344
     1345        // XXX not sure if these should result in errors or not...
    13561346        if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
    13571347            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
    … …  
    13661356
    13671357        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]);
     1358            psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    13701359            psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
    13711360            psFree(x);
    … …  
    14291418    // calculated, then call vectorSampleMean()
    14301419    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1431         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1420        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1421            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1422            return false;
     1423        }
    14321424    }
    14331425
    … …  
    14371429        stats->fittedStdev = NAN;
    14381430        psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1439         return false;
     1431        return true;
    14401432    }
    14411433
    … …  
    14681460        float max = statsMinMax->max;
    14691461        if (numValid == 0 || isnan(min) || isnan(max)) {
    1470             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1462            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
    14711463            psFree(statsMinMax);
    14721464            psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1473             return false;
     1465            return true;
    14741466        }
    14751467
    … …  
    15161508        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    15171509        if (binMin == binMax) {
    1518             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1510            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
    15191511            psFree(statsMinMax);
    1520             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1521             return false;
     1512            return true;
    15221513        }
    15231514
    … …  
    17241715    // calculated, then call vectorSampleMean()
    17251716    if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) {
    1726         vectorRobustStats(myVector, errors, mask, maskVal, stats);
     1717        if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
     1718            psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
     1719            return false;
     1720        }
    17271721    }
    17281722
    … …  
    17311725        stats->fittedStdev = NAN;
    17321726        stats->fittedStdev = NAN;
    1733         psTrace(TRACE, 4, "---- %s() end ----\n", __func__);
    1734         return false;
     1727        return true;
    17351728    }
    17361729
    … …  
    17631756        float max = statsMinMax->max;
    17641757        if (numValid == 0 || isnan(min) || isnan(max)) {
    1765             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1758            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
    17661759            psFree(statsMinMax);
    1767             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1768             return false;
     1760            stats->fittedStdev = NAN;
     1761            stats->fittedStdev = NAN;
     1762            return true;
    17691763        }
    17701764
    … …  
    17801774        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    17811775        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    1782             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics v4.\n");
     1776            psLogMsg(TRACE, PS_LOG_DETAIL, "Unable to generate histogram for fitted statistics v4.\n");
    17831777            psFree(histogram);
    17841778            psFree(statsMinMax);
    1785             return false;
     1779            stats->fittedStdev = NAN;
     1780            stats->fittedStdev = NAN;
     1781            return true;
    17861782        }
    17871783        if (psTraceGetLevel("psLib.math") >= 8) {
    … …  
    18131809        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
    18141810        if (binMin == binMax) {
    1815             psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
     1811            psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to calculate the min/max of the input vector.\n");
    18161812            psFree(statsMinMax);
    1817             psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1818             return false;
     1813            stats->fittedStdev = NAN;
     1814            stats->fittedStdev = NAN;
     1815            return true;
    18191816        }
    18201817
    … …  
    18771874
    18781875            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
     1876            // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
     1877            // we probably should be able to handle the data errors gracefully
    18791878            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    18801879            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
    … …  
    18831882
    18841883            if (!status) {
    1885                 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
     1884                psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a gaussian to the robust histogram.\n");
    18861885                psFree(poly);
    18871886                psFree(histogram);
    18881887                psFree(statsMinMax);
    1889                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1890                 return false;
     1888                stats->fittedStdev = NAN;
     1889                stats->fittedStdev = NAN;
     1890                return true;
    18911891            }
    18921892
    18931893            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]);
     1894                psLogMsg(TRACE, PS_LOG_MINUTIA, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    18961895                psFree(poly);
    18971896                psFree(histogram);
    … …  
    19071906                }
    19081907
    1909                 psError(PS_ERR_UNKNOWN, false, "fit did not converge\n");
    1910                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1911                 return false;
     1908                psLogMsg(TRACE, PS_LOG_DETAIL, "fit did not converge\n");
     1909                stats->fittedStdev = NAN;
     1910                stats->fittedStdev = NAN;
     1911                return true;
    19121912            }
    19131913
    … …  
    19751975
    19761976            if (!status) {
    1977                 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
     1977                psLogMsg(TRACE, PS_LOG_DETAIL, "Failed to fit a gaussian to the robust histogram.\n");
    19781978                psFree(poly);
    19791979                psFree(histogram);
    19801980                psFree(statsMinMax);
    1981                 psTrace(TRACE, 4, "---- %s(false) end  ----\n", __func__);
    1982                 return false;
     1981                stats->fittedStdev = NAN;
     1982                stats->fittedStdev = NAN;
     1983                return true;
    19831984            }
    19841985
    … …  
    21552156psStats* p_psStatsAlloc(const char *file, unsigned int lineno, const char *func, psStatsOptions options)
    21562157{
    2157     psTrace(TRACE, 3,"---- %s() begin  ----\n", __func__);
    2158 
    21592158    psStats *stats = p_psAlloc(file, lineno, func, sizeof(psStats));
    21602159    psMemSetDeallocator(stats, (psFreeFunc)statsFree);
    … …  
    21722171    stats->tmpMask = NULL;
    21732172
    2174     psTrace(TRACE, 3, "---- %s() end  ----\n", __func__);
    21752173    return stats;
    21762174}
    … …  
    22312229                   psVectorMaskType maskVal)
    22322230{
    2233     psTrace(TRACE, 3,"---- %s() begin  ----\n", __func__);
    22342231    PS_ASSERT_PTR_NON_NULL(stats, false);
    22352232    PS_ASSERT_VECTOR_NON_NULL(in, false);
    … …  
    23802377    psFree(errorsF32);
    23812378    psFree(maskVector);
    2382     psTrace(TRACE, 3,"---- %s() end  ----\n", __func__);
    23832379    return status;
    23842380}
    … …  
    27452741    return tmpFloat;
    27462742}
    2747 # endif
    27482743
    27492744/******************************************************************************
    … …  
    28732868    return tmpFloat;
    28742869}
     2870# endif
     2871
     2872/******************************************************************************
     2873fitQuadraticSearchForYThenReturnXusingValues(*xVec, *yVec, binNum, yVal): A general routine
     2874which fits a quadratic to three points and returns the x bin value corresponding to the input
     2875y-value.  This routine takes psVectors of x/y pairs as input, and fits a quadratic to the 3
     2876points surrounding element binNum in the vectors.  This version uses the values of x[i] for the
     2877x coordinates (not the midpoints).  This is appropriate for a cumulative histogram.  It then
     2878determines for what value x does that quadratic f(x) = yVal (the input parameter).
     2879
     2880XXX this function is used a fair amount in an inner loop: the polynomial fitting and evaluation
     2881could easily be done with statically allocated doubles, skipping the psLib versions of
     2882polynomial fitting, etc.
     2883
     2884*****************************************************************************/
     2885static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec,
     2886                                                 psVector *yVec,
     2887                                                 psS32 binNum,
     2888                                                 psF32 yVal
     2889    )
     2890{
     2891    psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
     2892    if (psTraceGetLevel("psLib.math") >= 8) {
     2893        PS_VECTOR_PRINT_F32(xVec);
     2894        PS_VECTOR_PRINT_F32(yVec);
     2895    }
     2896
     2897    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
     2898    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
     2899    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
     2900    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
     2901    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
     2902    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
     2903
     2904    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     2905    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     2906    psF32 tmpFloat = 0.0f;
     2907
     2908    if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
     2909        // The general case.  We have all three points.
     2910        x->data.F64[0] = binNum - 1;
     2911        x->data.F64[1] = binNum;
     2912        x->data.F64[2] = binNum + 1;
     2913        y->data.F64[0] = yVec->data.F32[binNum - 1];
     2914        y->data.F64[1] = yVec->data.F32[binNum];
     2915        y->data.F64[2] = yVec->data.F32[binNum + 1];
     2916        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]);
     2917        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
     2918        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
     2919
     2920        // Ensure that the y value lies within range of the y values.
     2921        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
     2922               ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
     2923            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     2924                    _("Specified yVal, %g, is not within y-range, %g to %g."),
     2925                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
     2926            return NAN;
     2927        }
     2928
     2929        // Ensure that the y values are monotonic.
     2930        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
     2931            ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
     2932            psError(PS_ERR_UNKNOWN, true,
     2933                    "This routine must be called with monotonically increasing or decreasing data points.\n");
     2934            psFree(x);
     2935            psFree(y);
     2936            return NAN;
     2937        }
     2938
     2939        // Determine the coefficients of the polynomial.
     2940        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
     2941        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
     2942            psError(PS_ERR_UNEXPECTED_NULL, false,
     2943                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
     2944                      "Returning NAN."));
     2945            psFree(x);
     2946            psFree(y);
     2947            return NAN;
     2948        }
     2949
     2950        psTrace(TRACE, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
     2951        psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
     2952        psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]);
     2953        psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n",
     2954                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
     2955                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]),
     2956                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2]));
     2957
     2958        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
     2959        float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
     2960        psFree(myPoly);
     2961
     2962        if (isnan(binValue)) {
     2963            psError(PS_ERR_UNEXPECTED_NULL,
     2964                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
     2965            psFree(x);
     2966            psFree(y);
     2967            return(NAN);
     2968        }
     2969
     2970        // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
     2971        assert (binValue >= binNum - 1);
     2972        assert (binValue <= binNum + 1);
     2973
     2974        int fitBin = binValue;
     2975        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
     2976        float dY = binValue - fitBin;
     2977        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
     2978    } else {
     2979        // These are special cases where the bin is at the beginning or end of the vector.
     2980        if (binNum == 0) {
     2981            // We have two points only at the beginning of the vectors x and y.
     2982            // X = (dX/dY)(Y - Yo) + Xo
     2983            float dX = xVec->data.F32[1] - xVec->data.F32[0];
     2984            float dY = yVec->data.F32[1] - yVec->data.F32[0];       
     2985            tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
     2986        } else if (binNum == (xVec->n - 1)) {
     2987            // We have two points only at the end of the vectors x and y.
     2988            // X = (dX/dY)(Y - Yo) + Xo
     2989            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
     2990            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];           
     2991            tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
     2992        }
     2993    }
     2994
     2995    psTrace(TRACE, 6, "FIT: return %f\n", tmpFloat);
     2996    psFree(x);
     2997    psFree(y);
     2998
     2999    return tmpFloat;
     3000}
    28753001
    28763002/******************************************************************************
Note: See TracChangeset for help on using the changeset viewer.