IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 28, 2013, 4:42:12 PM (13 years ago)
Author:
eugene
Message:

merge changes from trunk (czw stats fixes)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c

    r34703 r36256  
    140140}
    141141
     142// Debug information
     143#define CZW 0
    142144
    143145/*****************************************************************************/
    … …  
    172174*****************************************************************************/
    173175
    174 static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
     176// static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
     177static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal);
    175178
    176179/******************************************************************************
    … …  
    229232        }
    230233        count++;
     234
    231235    }
    232236    if (errors) {
    … …  
    793797        } else {
    794798            // Determine the bin size of the robust histogram, using the pre-defined number of bins
    795             binSize = (max - min) / INITIAL_NUM_BINS;
     799            binSize = (max - min) / INITIAL_NUM_BINS;
    796800        }
    797801        psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize);
    … …  
    876880        cumulative = psHistogramAlloc(min, max, numBins);
    877881        cumulative->nums->data.F32[0] = histogram->nums->data.F32[0];
    878         for (long i = 1; i < histogram->nums->n; i++) {
    879             cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i];
    880             cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i];
    881         }
     882        cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[1];
     883
     884        // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins).
     885        long Nc = 1;  // track the current bin of cumulative
     886        // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the
     887        // upper end of the current histogram bin
     888        for (long i = 1; i < histogram->nums->n - 1; i++) {
     889            if (histogram->nums->data.F32[i] == 0.0) continue;
     890            cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i];
     891            cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i+1];
     892            Nc ++;
     893        }
     894        long Nlast = Nc - 1;  // last valid cumulative bin
     895        for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled.
     896            cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast];
     897            cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0;
     898        }
     899       
    882900        if (psTraceGetLevel("psLib.math") >= 8) {
    883901            PS_VECTOR_PRINT_F32(cumulative->bounds);
    … …  
    895913
    896914        // ADD step 3: Interpolate to the exact 50% position in bin units
    897         stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    898         // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     915        // stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
     916        // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    899917        // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
     918        // There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there.
     919        // Instead, do a 5-point linear fit.
     920        stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
    900921
    901922        // convert bin to bin value: this is the robust histogram median.
    … …  
    912933        PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, 0);
    913934        PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, 0);
    914         PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022481f, 0);
    915         PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977519f, 0);
    916 
     935        PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0);
     936        PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0);
     937       
     938       
    917939        psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n",
    918940                binLo, binHi);
    … …  
    926948            goto escape;
    927949        }
    928 
     950   
    929951        // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma
    930952        // positions.
    … …  
    938960        // (extrapolation should not be needed and will result in errors)
    939961        float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32;
     962#if (0)
    940963        PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo,
    941964                            totalDataPoints * 0.158655f);
    … …  
    947970                            totalDataPoints * 0.691462f);
    948971        PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4,
    949                             totalDataPoints * 0.022481f);
     972                            totalDataPoints * 0.022750f);
    950973        PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4,
    951                             totalDataPoints * 0.977519f);
    952 
     974                            totalDataPoints * 0.977250f);
     975#else
     976        binLoF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo, totalDataPoints * 0.158655);
     977        binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345);         
     978        binL2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL2, totalDataPoints * 0.308538);
     979        binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462);         
     980        binL4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL4, totalDataPoints * 0.022750);
     981        binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250);
     982#endif 
    953983        // report +/- 1 sigma points
    954984        psTrace(TRACE, 5,
    … …  
    959989                binL2F32, binH2F32);
    960990        psTrace(TRACE, 5,
    961                 "The exact 02.22481 and 97.7519 percent data point positions are: (%f, %f)\n",
     991                "The exact 02.2275 and 97.7250 percent data point positions are: (%f, %f)\n",
    962992                binL4F32, binH4F32);
    963993
     994        // If some of the fits failed, attempt to fix this
     995        if (!isfinite(binLoF32) && isfinite(binHiF32)) { binLoF32 = -1.0 * binHiF32; }
     996        if (!isfinite(binHiF32) && isfinite(binLoF32)) { binHiF32 = -1.0 * binLoF32; }
     997        if (!isfinite(binL2F32) && isfinite(binH2F32)) { binL2F32 = -1.0 * binH2F32; }
     998        if (!isfinite(binH2F32) && isfinite(binL2F32)) { binH2F32 = -1.0 * binL2F32; }
     999        if (!isfinite(binL4F32) && isfinite(binH4F32)) { binL4F32 = -1.0 * binH4F32; }
     1000        if (!isfinite(binH4F32) && isfinite(binL4F32)) { binH4F32 = -1.0 * binL4F32; }
     1001       
    9641002        // ADD step 5: Determine SIGMA as the distance between binL2 and binH2 (+/- 0.5 sigma)
     1003
     1004
    9651005        float sigma1 = (binH2F32 - binL2F32);
    9661006        float sigma2 = (binHiF32 - binLoF32) / 2.0;
    9671007        float sigma4 = (binH4F32 - binL4F32) / 4.0;
    9681008
     1009        // Fix again?
     1010        if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4)) { sigma1 = (sigma2 + sigma4) / 2.0; }
     1011        if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4)) { sigma2 = (sigma1 + sigma4) / 2.0; }
     1012        if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1)) { sigma4 = (sigma2 + sigma1) / 2.0; }
     1013       
    9691014        // take the smallest of the three: if we have a clump with wide outliers, sigma2 and
    9701015        // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2
    9711016        // will be biased high.
    972         sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
     1017        //        sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4));
     1018        // CZW: Instead, take the median.  Taking the MIN forces a bias on unbiased data.
     1019        //      It seems like occasionally getting the wrong answer on a complex distribution
     1020        //      is more acceptable than always getting the wrong answer for simple ones.
     1021
     1022       
     1023        sigma = PS_MAX( PS_MIN(sigma1,sigma2),
     1024                        PS_MIN( PS_MAX(sigma1,sigma2),
     1025                                sigma4));
    9731026
    9741027        psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1);
    … …  
    9771030
    9781031        psTrace(TRACE, 6, "The current sigma is %f.\n", sigma);
    979         stats->robustStdev = sigma;
     1032        //        stats->robustStdev = sigma;
     1033        stats->robustStdev = sigma;
     1034
     1035#if (CZW && 0)
     1036        // Skewness check: Find least biased sample for each pair.
     1037        sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian,
     1038                              stats->robustMedian - binL2F32);
     1039        sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian,
     1040                              stats->robustMedian - binLoF32);
     1041        sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian,
     1042                              stats->robustMedian - binL4F32);
     1043        // Kurtosis check: Take median sample as the solution.
     1044        stats->robustStdev = PS_MAX( PS_MIN(sigma1,sigma2),
     1045                                     PS_MIN( PS_MAX(sigma1,sigma2),
     1046                                             sigma4));
     1047#endif
     1048
     1049       
     1050#if (CZW)
     1051        printf("CZW: bad sigma?: %f %f  %f %f  %f %f  %f %f %f  %f\n",
     1052               binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32,
     1053               sigma1,sigma2,sigma4,sigma);
     1054       
     1055        printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
     1056               iterate,
     1057           stats->robustMedian,stats->robustStdev,
     1058           fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]),
     1059
     1060           cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2],
     1061           cumulative->bounds->data.F32[binMedian-1],
     1062           cumulative->bounds->data.F32[binMedian],
     1063           cumulative->bounds->data.F32[binMedian+1],
     1064           cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3],
     1065
     1066           cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2],
     1067           cumulative->nums->data.F32[binMedian-1],
     1068           cumulative->nums->data.F32[binMedian],
     1069           cumulative->nums->data.F32[binMedian+1],
     1070           cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
     1071        //      PS_VECTOR_PRINT_F32(histogram->bounds);
     1072        //      PS_VECTOR_PRINT_F32(histogram->nums);
     1073        //      PS_VECTOR_PRINT_F32(cumulative->bounds);
     1074        //      PS_VECTOR_PRINT_F32(cumulative->nums);
     1075#endif
    9801076
    9811077        // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more
    … …  
    9831079        if (sigma < (3.0 * binSize)) {
    9841080            psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize);
    985             long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
    986             long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking
    987             psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index
    988             psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index
     1081
     1082            // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins
     1083            psF32 medianLo = stats->robustMedian - 25*binSize;
     1084            psF32 medianHi = stats->robustMedian + 25*binSize;
     1085
     1086            // long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region
     1087            // long maskHi = PS_MIN(cumulative->bounds->n - 1, (binMedian + 25)); // High index for masking
     1088            // psF32 medianLo = cumulative->bounds->data.F32[maskLo]; // Value at low index
     1089            // psF32 medianHi = cumulative->bounds->data.F32[maskHi]; // Value at high index
    9891090            psTrace(TRACE, 6, "Masking data more than 25 bins from the median\n");
    990             psTrace(TRACE, 6,
    991                     "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n",
    992                     binMedian, maskLo, maskHi);
     1091            // psTrace(TRACE, 6, "The median is at bin number %ld.  We mask bins outside the bin range (%ld:%ld)\n", binMedian, maskLo, maskHi);
    9931092            psTrace(TRACE, 6, "Masking data outside (%f %f)\n", medianLo, medianHi);
     1093            int Nmasked = 0;
    9941094            for (long i = 0 ; i < myVector->n ; i++) {
    9951095                if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) {
    996                     // XXXX is this correct?  is MASK_MARK safe?
     1096                    if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue;
    9971097                    mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK;
    9981098                    psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]);
     1099                    Nmasked ++;
    9991100                }
    10001101            }
     1102
     1103            if (Nmasked == 0) {
     1104                // no significant change to the sigma & binsize -- we are done here
     1105                iterate = -1;
     1106                continue;
     1107            }
    10011108
    10021109            // Free the histograms; they will be recreated on the next iteration, with new bounds
    … …  
    10301137        }
    10311138    }
    1032 
     1139   
    10331140    // XXX test lines while studying algorithm errors
    10341141    // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld  : %f %f %f\n",
    … …  
    10401147    PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0);
    10411148    PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0);
    1042     psTrace(TRACE, 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25);
     1149    psTrace(TRACE, 6, "The 25-percent and 75-percent data point bins are (%ld, %ld).\n", binLo25, binHi25);
    10431150
    10441151    // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile
    10451152    // positions.
    1046     psF32 binLo25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
    1047     psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
     1153    psF32 binLo25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f);
     1154    psF32 binHi25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f);
    10481155    if (isnan(binLo25F32) || isnan(binHi25F32)) {
    1049         COUNT_WARNING(10, 100, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n");
     1156        COUNT_WARNING(10, 100, "could not determine the robustUQ or LQ: fitLinearSearchForYThenReturnBin() returned a NAN.\n");
    10501157        goto escape;
    10511158    }
    … …  
    11001207 * "vectorFittedStats_v4" all versions of fitted stats now resolve to this function (only v4
    11011208 * has really been used) vectorFittedStats requires guess for fittedMean and fittedStdev
    1102  * robustN50 should also be set gaussian fit is performed using 2D polynomial to ln(y) this
     1209 * robustN50 should also be set gaussian fit is performed using 1D polynomial to ln(y) this
    11031210 * version follows the upper portion of the distribution until it passes 0.5*peak
    11041211 ********************/
    … …  
    11351242        return true;
    11361243    }
    1137 
     1244    if (myVector->n < 1) { printf("There are no elements in this vector.\n"); abort(); }
    11381245    float guessStdev = stats->robustStdev;  // pass the guess sigma
    11391246    float guessMean = stats->robustMedian;  // pass the guess mean
    … …  
    11551262            // set roughly so that the lowest bins have about 2 cnts
    11561263            // Nsmallest ~ N50 / (4*dN))
    1157             psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8));
    1158             binSize = guessStdev / dN;
     1264            psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); 
     1265           binSize = guessStdev / dN;
    11591266        }
    11601267
    … …  
    11821289        // XXX can we calculate the binMin, binMax **before** building this histogram?
    11831290        // if the range is too absurd, adjust numBins & binSize
     1291        // We no longer want to reset the binSize here, as it can cause odd things.  Better to select
     1292        // a number of bins, and then set the min/max values to put those bins sanely around the mean.
    11841293        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
    1185         binSize = (max - min) / (float) numBins;
     1294        //        binSize = (max - min) / (float) numBins;
    11861295        psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
    11871296        psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
    11881297        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
    11891298
     1299
     1300#define FITTED_CLIPPING_NUM 5.0
     1301        if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev) {
     1302          min = guessMean - FITTED_CLIPPING_NUM * guessStdev;
     1303        }
     1304        if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev) {
     1305          max = guessMean + FITTED_CLIPPING_NUM * guessStdev;
     1306        }
     1307       
    11901308        psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
    11911309        if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
    … …  
    12221340        PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, 0);
    12231341        PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0);
     1342
    12241343        if (binMin == binMax) {
    12251344            COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
    … …  
    12481367        psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]);
    12491368
     1369       
    12501370        float lowfitMean = NAN;
    12511371        float lowfitStdev = NAN;
    … …  
    12851405            }
    12861406            y->n = x->n = j;
    1287 
     1407           
    12881408            // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2
    12891409            // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
    … …  
    12971417                psErrorClear();
    12981418                COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
     1419
    12991420                psFree(poly);
    13001421                psFree(histogram);
    … …  
    13781499            psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
    13791500            bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
     1501#if (CZW && 0)
     1502            for (long i = 0; i < x->n; i++) {
     1503              printf("CZW: Dcheck: %ld %f %f %f\n",
     1504                     i,x->data.F32[i],y->data.F32[i],
     1505                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
     1506                     poly->coeff[2] * pow(x->data.F32[i],2));
     1507            }
     1508#endif
    13801509            psFree(x);
    13811510            psFree(y);
    … …  
    13931522            fullfitStdev = sqrt(-0.5/poly->coeff[2]);
    13941523            fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev);
     1524
    13951525#ifndef PS_NO_TRACE
    13961526            psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
    … …  
    14151545            }
    14161546
     1547           
    14171548            psFree (poly);
    14181549        }
    … …  
    14371568            done = true;
    14381569        }
     1570
     1571       
     1572#if (CZW && 1)
     1573        printf("CZW IN FITTED: iter   %d %f \n"
     1574               "               low    %f %f \n"
     1575               "               full   %f %f \n"
     1576               "               robust %f %f \n"
     1577               "               final  %f %f\n",
     1578               iteration,minValueSym,
     1579               lowfitMean,lowfitStdev,
     1580               fullfitMean,fullfitStdev,
     1581               stats->robustMedian,stats->robustStdev,
     1582               guessMean,guessStdev);
     1583#endif
    14391584
    14401585        // Clean up after fitting
    … …  
    19652110// other private functions used above
    19662111
     2112# if (0)
    19672113static psF32 QuadraticInverse(psF32 a,
    19682114                              psF32 b,
    … …  
    19862132    return 0.5 * (xLo + xHi);
    19872133}
     2134
     2135static psF32 LinearInverse(psF32 a,
     2136                           psF32 b,
     2137                           psF32 y,
     2138                           psF32 xLo,
     2139                           psF32 xHi
     2140    )
     2141{
     2142    psF64 x = (y - b) / a;
     2143
     2144    if (xLo <= x && x <= xHi) {
     2145        return x;
     2146    }
     2147    return 0.5 * (xLo + xHi);
     2148}
     2149# endif
    19882150
    19892151# if (0)
    … …  
    22422404    return tmpFloat;
    22432405}
    2244 # endif
    22452406
    22462407/******************************************************************************
    … …  
    22762437    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
    22772438
    2278     psVector *x = psVectorAlloc(3, PS_TYPE_F64);
    2279     psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     2439    //    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     2440    //    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     2441    psVector *x = psVectorAlloc(5, PS_TYPE_F64);
     2442    psVector *y = psVectorAlloc(5, PS_TYPE_F64);
    22802443    psF32 tmpFloat = 0.0f;
    22812444
    2282     if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
     2445    //    if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) {
     2446    if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) {
    22832447        // The general case.  We have all three points.
    2284         x->data.F64[0] = binNum - 1;
    2285         x->data.F64[1] = binNum;
    2286         x->data.F64[2] = binNum + 1;
    2287         y->data.F64[0] = yVec->data.F32[binNum - 1];
    2288         y->data.F64[1] = yVec->data.F32[binNum];
    2289         y->data.F64[2] = yVec->data.F32[binNum + 1];
     2448      //        x->data.F64[0] = binNum - 1;
     2449      //        x->data.F64[1] = binNum;
     2450      //        x->data.F64[2] = binNum + 1;
     2451      x->data.F64[0] = xVec->data.F32[binNum - 2];
     2452      x->data.F64[1] = xVec->data.F32[binNum - 1];
     2453      x->data.F64[2] = xVec->data.F32[binNum + 0];
     2454      x->data.F64[3] = xVec->data.F32[binNum + 1];
     2455      x->data.F64[4] = xVec->data.F32[binNum + 2];
     2456        y->data.F64[0] = yVec->data.F32[binNum - 2];
     2457        y->data.F64[1] = yVec->data.F32[binNum - 1];
     2458        y->data.F64[2] = yVec->data.F32[binNum + 0];
     2459        y->data.F64[3] = yVec->data.F32[binNum + 1];
     2460        y->data.F64[4] = yVec->data.F32[binNum + 2];
    22902461        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]);
    22912462        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
    22922463        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
    22932464
     2465
    22942466        // Ensure that the y value lies within range of the y values.
    2295         if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
    2296                ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
     2467        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
     2468               ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) {
    22972469            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    22982470                    _("Specified yVal, %g, is not within y-range, %g to %g."),
    … …  
    23312503
    23322504        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
    2333         float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]);
     2505        float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[4]);
    23342506        psFree(myPoly);
    23352507
    … …  
    23412513            return(NAN);
    23422514        }
    2343 
     2515       
    23442516        // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
    2345         assert (binValue >= binNum - 1);
    2346         assert (binValue <= binNum + 1);
    2347 
    2348         int fitBin = binValue;
    2349         float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
    2350         float dY = binValue - fitBin;
    2351         tmpFloat = xVec->data.F32[fitBin] + dY * dX;
     2517        //      assert (binValue >= binNum - 1);
     2518        //      assert (binValue <= binNum + 1);
     2519
     2520        //      int fitBin = binValue;
     2521        //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
     2522        //        float dY = binValue - fitBin;
     2523        //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
     2524        tmpFloat = binValue;
     2525       
    23522526    } else {
    23532527        // These are special cases where the bin is at the beginning or end of the vector.
    … …  
    23812555    return tmpFloat;
    23822556}
     2557# endif
     2558
     2559
     2560/******************************************************************************
     2561fitQuadraticSearchForYThenReturnXusingValues(*xVec, *yVec, binNum, yVal): A general routine
     2562which fits a quadratic to three points and returns the x bin value corresponding to the input
     2563y-value.  This routine takes psVectors of x/y pairs as input, and fits a quadratic to the 3
     2564points surrounding element binNum in the vectors.  This version uses the values of x[i] for the
     2565x coordinates (not the midpoints).  This is appropriate for a cumulative histogram.  It then
     2566determines for what value x does that quadratic f(x) = yVal (the input parameter).
     2567
     2568XXX this function is used a fair amount in an inner loop: the polynomial fitting and evaluation
     2569could easily be done with statically allocated doubles, skipping the psLib versions of
     2570polynomial fitting, etc.
     2571
     2572*****************************************************************************/
     2573static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec,
     2574                                              psVector *yVec,
     2575                                              psS32 binNum,
     2576                                              psF32 yVal
     2577    )
     2578{
     2579
     2580# if (1)
     2581# define HALF_SIZE 2
     2582  double Sx = 0.0;
     2583
     2584  double Sy = 0.0;
     2585  double Sxx = 0.0;
     2586  double Sxy = 0.0;
     2587  double deltaY = 0.0;
     2588  int N = 0;
     2589
     2590  for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) {
     2591    if ((u >= 0)&&(u < yVec->n)) {
     2592      if (u+1 < xVec->n) {
     2593        Sx += yVec->data.F32[u];
     2594        Sxx += PS_SQR(yVec->data.F32[u]);
     2595
     2596        deltaY = xVec->data.F32[u];
     2597        //deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]);
     2598        Sy += deltaY;
     2599        Sxy += yVec->data.F32[u] * deltaY;
     2600        N += 1;
     2601      }
     2602    }
     2603  }
     2604  double Det = N * Sxx - Sx * Sx;
     2605  if (Det == 0.0) return NAN;
     2606  if (N == 0) return NAN;
     2607
     2608  double C0 = (Sy*Sxx - Sx*Sxy) / Det;
     2609  double C1 = (Sxy*N - Sx*Sy) / Det;
     2610 
     2611  double value = C0 + yVal*C1;
     2612  return value;
     2613 
     2614 
     2615# else
     2616    psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal);
     2617    if (psTraceGetLevel("psLib.math") >= 8) {
     2618        PS_VECTOR_PRINT_F32(xVec);
     2619        PS_VECTOR_PRINT_F32(yVec);
     2620    }
     2621
     2622    PS_ASSERT_VECTOR_NON_NULL(xVec, NAN);
     2623    PS_ASSERT_VECTOR_NON_NULL(yVec, NAN);
     2624    PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN);
     2625    PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN);
     2626    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN);
     2627    PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN);
     2628
     2629    //    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     2630    //    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     2631    psVector *x = psVectorAlloc(5, PS_TYPE_F64);
     2632    psVector *y = psVectorAlloc(5, PS_TYPE_F64);
     2633    psF32 tmpFloat = 0.0f;
     2634
     2635    if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) {
     2636        x->data.F64[0] = xVec->data.F32[binNum - 2];
     2637        x->data.F64[1] = xVec->data.F32[binNum - 1];
     2638        x->data.F64[2] = xVec->data.F32[binNum + 0];
     2639        x->data.F64[3] = xVec->data.F32[binNum + 1];
     2640        x->data.F64[4] = xVec->data.F32[binNum + 2];
     2641
     2642        y->data.F64[0] = yVec->data.F32[binNum - 2];
     2643        y->data.F64[1] = yVec->data.F32[binNum - 1];
     2644        y->data.F64[2] = yVec->data.F32[binNum + 0];
     2645        y->data.F64[3] = yVec->data.F32[binNum + 1];
     2646        y->data.F64[4] = yVec->data.F32[binNum + 2];
     2647        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]);
     2648        psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
     2649        psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
     2650
     2651        // Ensure that the y value lies within range of the y values.
     2652        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) ||
     2653               ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) {
     2654            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     2655                    _("Specified yVal, %g, is not within y-range, %g to %g."),
     2656                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
     2657            return NAN;
     2658        }
     2659
     2660        // Ensure that the y values are monotonic.
     2661        if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) ||
     2662            ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) {
     2663            psError(PS_ERR_UNKNOWN, true,
     2664                    "This routine must be called with monotonically increasing or decreasing data points.\n");
     2665            psFree(x);
     2666            psFree(y);
     2667            return NAN;
     2668        }
     2669
     2670        // Determine the coefficients of the polynomial.
     2671        psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1);
     2672        if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) {
     2673            psError(PS_ERR_UNEXPECTED_NULL, false,
     2674                    _("Failed to fit a 1-dimensional polynomial to the three specified data points.  "
     2675                      "Returning NAN."));
     2676            psFree(x);
     2677            psFree(y);
     2678            return NAN;
     2679        }
     2680
     2681        psTrace(TRACE, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]);
     2682        psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]);
     2683        psTrace(TRACE, 6, "Fitted y vec is (%f %f)\n",
     2684                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]),
     2685                (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]));
     2686
     2687        psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal);
     2688        float binValue = LinearInverse(myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[4]);
     2689        psFree(myPoly);
     2690
     2691        if (isnan(binValue)) {
     2692            psError(PS_ERR_UNEXPECTED_NULL,
     2693                    false, _("Failed to determine the median of the fitted polynomial.  Returning NAN."));
     2694            psFree(x);
     2695            psFree(y);
     2696            return(NAN);
     2697        }
     2698       
     2699        // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1
     2700        //      assert (binValue >= binNum - 1);
     2701        //      assert (binValue <= binNum + 1);
     2702
     2703        //      int fitBin = binValue;
     2704        //        float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin];
     2705        //        float dY = binValue - fitBin;
     2706        //        tmpFloat = xVec->data.F32[fitBin] + dY * dX;
     2707        tmpFloat = binValue;
     2708               
     2709       
     2710    } else {
     2711        // These are special cases where the bin is at the beginning or end of the vector.
     2712        if (binNum == 0) {
     2713            // We have two points only at the beginning of the vectors x and y.
     2714            // X = (dX/dY)(Y - Yo) + Xo
     2715            float dX = xVec->data.F32[1] - xVec->data.F32[0];
     2716            float dY = yVec->data.F32[1] - yVec->data.F32[0];
     2717            if (dY == 0.0) {
     2718                tmpFloat = xVec->data.F32[0];
     2719            } else {
     2720                tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
     2721            }
     2722        } else if (binNum == (xVec->n - 1)) {
     2723            // We have two points only at the end of the vectors x and y.
     2724            // X = (dX/dY)(Y - Yo) + Xo
     2725            float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1];
     2726            float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1];
     2727            if (dY == 0.0) {
     2728                tmpFloat = xVec->data.F32[binNum-1];
     2729            } else {
     2730                tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
     2731            }
     2732        }
     2733    }
     2734
     2735    psTrace(TRACE, 6, "FIT: return %f\n", tmpFloat);
     2736    psFree(x);
     2737    psFree(y);
     2738
     2739    return tmpFloat;
     2740# endif
     2741}
Note: See TracChangeset for help on using the changeset viewer.