Changeset 24987 for trunk/psLib/src/math/psStats.c
- Timestamp:
- Aug 3, 2009, 2:13:09 PM (17 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/math/psStats.c (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/math/psStats.c
r24884 r24987 129 129 130 130 # define COUNT_WARNING(LIMIT, INTERVAL, ...) { \ 131 132 133 134 135 136 131 static int nCalls = 1; \ 132 if (nCalls < LIMIT) { \ 133 psWarning(__VA_ARGS__); \ 134 } \ 135 if (!(nCalls % INTERVAL)) { \ 136 psWarning(__VA_ARGS__); \ 137 137 psWarning("(warning raised %d times)", nCalls); \ 138 139 138 } \ 139 nCalls ++; \ 140 140 } 141 141 142 142 143 143 /*****************************************************************************/ … … 217 217 for (long i = 0; i < numData; i++) { 218 218 // Check if the data is with the specified range 219 if (!isfinite(data[i])) 220 219 if (!isfinite(data[i])) 220 continue; 221 221 if (useRange && (data[i] < stats->min)) 222 222 continue; … … 242 242 243 243 if (!isnan(mean)) { 244 244 stats->results |= PS_STAT_SAMPLE_MEAN; 245 245 } 246 246 return true; … … 278 278 for (long i = 0; i < num; i++) { 279 279 // Check if the data is with the specified range 280 if (!isfinite(vector[i])) 281 280 if (!isfinite(vector[i])) 281 continue; 282 282 if (useRange && (vector[i] < stats->min)) 283 283 continue; … … 329 329 // into the temporary vectors. 330 330 for (long i = 0; i < inVector->n; i++) { 331 if (!isfinite(input[i])) 332 331 if (!isfinite(input[i])) 332 continue; 333 333 if (useRange && (input[i] < stats->min)) 334 334 continue; … … 353 353 // Sort the temporary vector. 354 354 if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK) 355 355 // an error in psVectorSort is a serious error 356 356 psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data.")); 357 357 stats->sampleUQ = NAN; … … 407 407 // If the mean is NAN, then generate a warning and set the stdev to NAN. 408 408 if (isnan(stats->sampleMean)) { 409 409 COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN."); 410 410 stats->sampleStdev = NAN; 411 411 return true; … … 425 425 for (long i = 0; i < myVector->n; i++) { 426 426 // Check if the data is with the specified range 427 if (!isfinite(data[i])) 428 427 if (!isfinite(data[i])) 428 continue; 429 429 if (useRange && (data[i] < stats->min)) { 430 430 continue; … … 502 502 for (long i = 0; i < myVector->n; i++) { 503 503 // Check if the data is with the specified range 504 if (!isfinite(data[i])) 505 504 if (!isfinite(data[i])) 505 continue; 506 506 if (useRange && (data[i] < stats->min)) { 507 507 continue; … … 771 771 stats->robustLQ = min; 772 772 stats->robustN50 = numValid; 773 773 // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits: 774 774 stats->results |= PS_STAT_ROBUST_MEDIAN; 775 775 stats->results |= PS_STAT_ROBUST_STDEV; … … 801 801 // XXXXX we need to consider this step if errors -> variance 802 802 if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { 803 803 // if psVectorHistogram returns false, we have a programming error 804 804 psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n"); 805 806 807 808 809 810 805 psFree(histogram); 806 psFree(cumulative); 807 psFree(statsMinMax); 808 psFree(mask); 809 810 return false; 811 811 } 812 812 if (psTraceGetLevel("psLib.math") >= 8) { … … 838 838 839 839 // ADD step 3: Interpolate to the exact 50% position in bin units 840 841 842 843 844 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. 845 845 if (isnan(stats->robustMedian)) { 846 846 COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n"); … … 1054 1054 if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { 1055 1055 psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); 1056 1057 1056 return false; 1057 } 1058 1058 } 1059 1059 … … 1107 1107 psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) 1108 1108 if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { 1109 1109 // if psVectorHistogram returns false, we have a programming error 1110 1110 psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n"); 1111 1111 psFree(histogram); … … 1172 1172 PS_VECTOR_PRINT_F32(y); 1173 1173 } 1174 1175 1174 1175 // psMinimizeLMChi2 can return false for bad data as well as for serious failures 1176 1176 if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { 1177 1177 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); … … 1235 1235 if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { 1236 1236 psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); 1237 1238 1237 return false; 1238 } 1239 1239 } 1240 1240 … … 1355 1355 // psVectorInit (fitMask, 0); 1356 1356 1357 1357 // XXX not sure if these should result in errors or not... 1358 1358 if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) { 1359 1359 psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); … … 1432 1432 if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { 1433 1433 psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); 1434 1435 1434 return false; 1435 } 1436 1436 } 1437 1437 … … 1729 1729 if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { 1730 1730 psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); 1731 1732 1731 return false; 1732 } 1733 1733 } 1734 1734 … … 1766 1766 COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); 1767 1767 psFree(statsMinMax); 1768 1768 goto escape; 1769 1769 } 1770 1770 … … 1774 1774 stats->fittedMean = min; 1775 1775 stats->fittedStdev = 0.0; 1776 1777 1776 stats->results |= PS_STAT_FITTED_MEAN_V4; 1777 stats->results |= PS_STAT_FITTED_STDEV_V4; 1778 1778 return true; 1779 1779 } … … 1793 1793 psFree(histogram); 1794 1794 psFree(statsMinMax); 1795 1795 goto escape; 1796 1796 } 1797 1797 if (psTraceGetLevel("psLib.math") >= 8) { … … 1825 1825 COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); 1826 1826 psFree(statsMinMax); 1827 1827 goto escape; 1828 1828 } 1829 1829 … … 1886 1886 1887 1887 // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 1888 1889 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 1890 1890 psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); 1891 1891 bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); … … 1898 1898 psFree(histogram); 1899 1899 psFree(statsMinMax); 1900 1900 goto escape; 1901 1901 } 1902 1902 … … 1917 1917 1918 1918 COUNT_WARNING(10, 100, "fit did not converge\n"); 1919 1919 goto escape; 1920 1920 } 1921 1921 … … 1987 1987 psFree(histogram); 1988 1988 psFree(statsMinMax); 1989 1989 goto escape; 1990 1990 } 1991 1991 … … 2898 2898 *****************************************************************************/ 2899 2899 static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, 2900 2901 2902 2900 psVector *yVec, 2901 psS32 binNum, 2902 psF32 yVal 2903 2903 ) 2904 2904 { … … 2982 2982 } 2983 2983 2984 2985 2986 2987 2988 2989 2990 2991 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 2992 } else { 2993 2993 // These are special cases where the bin is at the beginning or end of the vector. 2994 2994 if (binNum == 0) { 2995 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 tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; 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 } 3000 3004 } else if (binNum == (xVec->n - 1)) { 3001 3005 // We have two points only at the end of the vectors x and y. 3002 // X = (dX/dY)(Y - Yo) + Xo 3003 float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1]; 3004 float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1]; 3005 tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1]; 3006 } 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 } 3007 3015 } 3008 3016
Note:
See TracChangeset
for help on using the changeset viewer.
