Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 6270) +++ trunk/psLib/src/math/psStats.c (revision 6305) @@ -3,5 +3,5 @@ * @ingroup Stat * - * This file will hold the definition of the histogram and stats data + * This file will hold the definitions of the histogram and stats data * structures. It also contains prototypes for procedures which operate * on those data structures. @@ -10,18 +10,12 @@ * * XXX: The following stats members are never used, or set in this code. - * stats->robustN50 * stats->clippedNvalues - * stats->binsize * * XXX: Must do * nSubsample points * use ->min and ->max (PS_STAT_USE_RANGE) - * use ->binsize (PS_STAT_USE_BINSIZE) * - * - * - * - * @version $Revision: 1.164 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-01 02:07:24 $ + * @version $Revision: 1.165 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-02 21:09:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -49,4 +43,5 @@ #include "psPolynomial.h" #include "psConstants.h" +#include "psMathUtils.h" #include "psErrorText.h" @@ -61,5 +56,5 @@ #define PS_CLIPPED_SIGMA_LB 1.0 #define PS_CLIPPED_SIGMA_UB 10.0 -#define PS_POLY_MEDIAN_MAX_ITERATIONS 10 +#define PS_POLY_MEDIAN_MAX_ITERATIONS 30 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \ @@ -847,6 +842,4 @@ // Calculate the quartile points exactly. - // XXX: We should probably do a bin-midpoint if the quartile is not an - // integer. stats->sampleUQ = sortedVector->data.F32[3 * (nValues / 4)]; stats->sampleLQ = sortedVector->data.F32[nValues / 4]; @@ -859,102 +852,4 @@ } -/****************************************************************************** -p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the -stdev of the input vector. -Inputs - myVector - maskVector - maskVal - stats -Returns - NULL - -XXX: remove this - *****************************************************************************/ -void p_psVectorSampleStdevOLD(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psS32 i = 0; // Loop index variable - psS32 countInt = 0; // # of data points being used - psF32 countFloat = 0.0; // # of data points being used - psF32 mean = 0.0; // The mean - psF32 diff = 0.0; // Used in calculating stdev - psF32 sumSquares = 0.0; // temporary variable - psF32 sumDiffs = 0.0; // temporary variable - - // This procedure requires the mean. If it has not been already - // calculated, then call p_psVectorSampleMean() - if (0 != isnan(stats->sampleMean)) { - p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats); - } - // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (0 != isnan(stats->sampleMean)) { - stats->sampleStdev = NAN; - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): p_psVectorSampleMean() reported a NAN mean.\n"); - psTrace(__func__, 4, "---- %s() end ----\n", __func__); - return; - } - - mean = stats->sampleMean; - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); - sumDiffs += diff; - countInt++; - } - } - } else { - for (i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); - sumDiffs += diff; - countInt++; - } - } - countInt = myVector->n; - } - } else { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); - sumDiffs += diff; - countInt++; - } - } - } else { - for (i = 0; i < myVector->n; i++) { - diff = myVector->data.F32[i] - mean; - sumSquares += (diff * diff); - sumDiffs += diff; - countInt++; - } - countInt = myVector->n; - } - } - if (countInt == 0) { - stats->sampleStdev = NAN; - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): no valid psVector elements (%d). Setting stats->sampleStdev = NAN.\n", countInt); - } else if (countInt == 1) { - stats->sampleStdev = 0.0; - psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleStdev(): only one valid psVector elements (%d). Setting stats->sampleStdev = 0.0.\n", countInt); - } else { - countFloat = (psF32)countInt; - stats->sampleStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - } - psTrace(__func__, 4, "---- %s() end ----\n", __func__); -} /****************************************************************************** p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the @@ -1094,9 +989,5 @@ -2: warning -XXX: Do we really need to calculate median and mean? - XXX: Use static vectors for tmpMask. - -XXX: This has not been tested. *****************************************************************************/ psS32 p_psVectorClippedStats(const psVector* myVector, @@ -1223,7 +1114,9 @@ // Otherwise, use the new results and continue. if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) { - // Exit loop. XXX: throw an error/warning here? + // Exit loop. iter = stats->clipIter; - rc = -1; + psError(PS_ERR_UNKNOWN, true, "p_psVectorSampleMean() or p_psVectorSampleStdev() returned a NAN.\n"); + psFree(tmpMask); + return(-1); } else { clippedMean = statsTmp->sampleMean; @@ -1246,116 +1139,6 @@ } -/***************************************************************************** -These macros and functions define the following functions: - - p_psNormalizeVectorRange(myData, low, high) - -That assumes that the low/high arguments are PS_TYPE_F64; the vector myData -can be of any type. Arguments low/high will be converted to the appropriate -type and one of the type-specific functions below will be called: - - p_psNormalizeVectorRangeU8(myData, low, high) - p_psNormalizeVectorRangeU16(myData, low, high) - p_psNormalizeVectorRangeU32(myData, low, high) - p_psNormalizeVectorRangeU64(myData, low, high) - p_psNormalizeVectorRangeS8(myData, low, high) - p_psNormalizeVectorRangeS16(myData, low, high) - p_psNormalizeVectorRangeS32(myData, low, high) - p_psNormalizeVectorRangeS64(myData, low, high) - p_psNormalizeVectorRangeF32(myData, low, high) - p_psNormalizeVectorRangeF64(myData, low, high) - *****************************************************************************/ -#define PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(TYPE) \ -void p_psNormalizeVectorRange##TYPE(psVector* myData, \ - ps##TYPE outLow, \ - ps##TYPE outHigh) \ -{ \ - ps##TYPE min = (ps##TYPE) PS_MAX_##TYPE; \ - ps##TYPE max = (ps##TYPE) -PS_MAX_##TYPE; \ - psS32 i = 0; \ - \ - for (i = 0; i < myData->n; i++) { \ - if (myData->data.TYPE[i] < min) { \ - min = myData->data.TYPE[i]; \ - } \ - if (myData->data.TYPE[i] > max) { \ - max = myData->data.TYPE[i]; \ - } \ - } \ - \ - /* Ensure that max!=min before we divide by (max-min) */ \ - if (max != min) { \ - for (i = 0; i < myData->n; i++) { \ - myData->data.TYPE[i] = (outLow + (myData->data.TYPE[i] - min) * \ - (outHigh - outLow) / (max - min)); \ - } \ - } else { \ - psLogMsg(__func__, PS_LOG_WARN, "WARNING: (max==min). Setting all elements to min.\n"); \ - for (i = 0; i < myData->n; i++) \ - { \ - \ - myData->data.TYPE[i] = outLow; \ - \ - } \ - } \ -} \ - -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U8) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U16) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(U64) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S8) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S16) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(S64) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F32) -PS_FUNC_MACRO_NORMALIZE_VECTOR_RANGE(F64) - -void p_psNormalizeVectorRange(psVector* myData, - psF64 outLow, - psF64 outHigh) -{ - psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - switch (myData->type.type) { - case PS_TYPE_U8: - p_psNormalizeVectorRangeU8(myData, (psU8) outLow, (psU8) outHigh); - break; - case PS_TYPE_U16: - p_psNormalizeVectorRangeU16(myData, (psU16) outLow, (psU16) outHigh); - break; - case PS_TYPE_U32: - p_psNormalizeVectorRangeU32(myData, (psU32) outLow, (psU32) outHigh); - break; - case PS_TYPE_U64: - p_psNormalizeVectorRangeU64(myData, (psU64) outLow, (psU64) outHigh); - break; - case PS_TYPE_S8: - p_psNormalizeVectorRangeS8(myData, (psS8) outLow, (psS8) outHigh); - break; - case PS_TYPE_S16: - p_psNormalizeVectorRangeS16(myData, (psS16) outLow, (psS16) outHigh); - break; - case PS_TYPE_S32: - p_psNormalizeVectorRangeS32(myData, (psS32) outLow, (psS32) outHigh); - break; - case PS_TYPE_S64: - p_psNormalizeVectorRangeS64(myData, (psS64) outLow, (psS64) outHigh); - break; - case PS_TYPE_F32: - p_psNormalizeVectorRangeF32(myData, (psF32) outLow, (psF32) outHigh); - break; - case PS_TYPE_F64: - p_psNormalizeVectorRangeF64(myData, (psF64) outLow, (psF64) outHigh); - break; - case PS_TYPE_C32: - case PS_TYPE_C64: - default: - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - "Unallowable operation: %s has incorrect type.", - myData); - break; - } - psTrace(__func__, 4, "---- %s() end ----\n", __func__); -} + + /****************************************************************************** @@ -1371,8 +1154,9 @@ XXX: Create a 2nd-order polynomial version and solve for X analytically. *****************************************************************************/ -psF32 p_ps1DPolyMedian(psPolynomial1D* myPoly, - psF32 rangeLow, - psF32 rangeHigh, - psF32 getThisValue) +psF32 p_ps1DPolyMedian( + psPolynomial1D* myPoly, + psF32 rangeLow, + psF32 rangeHigh, + psF32 getThisValue) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -1405,4 +1189,5 @@ f = psPolynomial1DEval(myPoly, midpoint); + printf("p_ps1DPolyMedian(): f(%f) is %f\n", midpoint, f); if (fabs(f - getThisValue) <= FLT_EPSILON) { return (midpoint); @@ -1419,4 +1204,65 @@ return (midpoint); } + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +#define PS_SQRT(x) sqrt(x) +static psF32 QuadraticInverse( + psF32 a, + psF32 b, + psF32 c, + psF32 y) +{ + psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a)); + + psF64 x1 = -b/(2.0*a) + tmp; + psF64 x2 = -b/(2.0*a) - tmp; + + psF64 y1 = (a * x1 * x1) + (b * x1) + c; + psF64 y2 = (a * x2 * x2) + (b * x2) + c; + + printf("QuadraticInverse: %fx^2 + %fx + %f\n", a, b, c); + printf("QuadraticInverse: y is %f\n", y); + printf("QuadraticInverse: (x1, x2) is (%f %f)\n", x1, x2); + printf("QuadraticInverse: (y1, y2) is (%f %f)\n", y1, y2); + + return(x1); +} + /****************************************************************************** @@ -1438,8 +1284,9 @@ arbitrary vectors. We should probably test that condition. *****************************************************************************/ -psF32 fitQuadraticSearchForYThenReturnX(psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal) +psF32 fitQuadraticSearchForYThenReturnX( + psVector *xVec, + psVector *yVec, + psS32 binNum, + psF32 yVal) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -1461,5 +1308,4 @@ psVector *y = psVectorAlloc(3, PS_TYPE_F64); psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); - psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); psF32 tmpFloat = 0.0f; @@ -1489,5 +1335,4 @@ if (!(y->data.F64[1] <= y->data.F64[2])) { psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); - psFree(myPoly); psFree(x); psFree(y); @@ -1505,5 +1350,4 @@ if (!(y->data.F64[1] >= y->data.F64[2])) { psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); - psFree(myPoly); psFree(x); psFree(y); @@ -1525,5 +1369,5 @@ // Determine the coefficients of the polynomial. - // myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); + psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x); if (myPoly == NULL) { @@ -1531,5 +1375,4 @@ false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT); - psFree(myPoly); psFree(x); psFree(y); @@ -1541,5 +1384,6 @@ psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]); psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]); - psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n", (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]), + psTrace(__func__, 6, "Fitted y vec is (%f %f %f)\n", + (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]), (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]), (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2])); @@ -1548,5 +1392,11 @@ // polynomial, looking for the value x such that f(x) = yVal psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); + printf("Calling p_ps1DPolyMedian() for the y-value (%.2f) that has an x in the range (%.2f - %.2f)\n", yVal, x->data.F64[0], x->data.F64[2]); tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal); + printf("Cool, tmpFloat is %.2f\n", tmpFloat); + + QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal); + psFree(myPoly); + if (isnan(tmpFloat)) { @@ -1554,5 +1404,4 @@ false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN); - psFree(myPoly); psFree(x); psFree(y); @@ -1581,5 +1430,4 @@ psTrace(__func__, 6, "FIT: return %f\n", tmpFloat); - psFree(myPoly); psFree(x); psFree(y); @@ -1591,9 +1439,10 @@ /****************************************************************************** -XXX: We assume unnormalized gaussians. +NOTE: We assume unnormalized gaussians. *****************************************************************************/ -psF32 psMinimizeLMChi2Gauss1D(psVector *deriv, - const psVector *params, - const psVector *coords) +psF32 psMinimizeLMChi2Gauss1D( + psVector *deriv, + const psVector *params, + const psVector *coords) { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); @@ -1627,4 +1476,5 @@ } +//XXX: Use the psLib function? psF32 LinInterpolate(psF32 x0, psF32 x1, @@ -1643,76 +1493,4 @@ return(x0 + y * (x1 - x0) / (y1 - y0)); } - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - @@ -1754,5 +1532,5 @@ XXX: Review and ensure that all memory is free'ed at premature exits. *****************************************************************************/ -#define INITIAL_NUM_BINS 1000.0 +#define INITIAL_NUM_BINS 500.0 psS32 p_psVectorRobustStats(const psVector* myVector, const psVector* errors, @@ -1770,6 +1548,6 @@ psHistogram *cumulativeRobustHistogram = NULL; psS32 numBins = 0; - psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32); - tmpScalar->type.type = PS_TYPE_F32; + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F32; psS32 totalDataPoints = 0; psS32 rc = 0; @@ -1806,5 +1584,4 @@ psFree(tmpStatsMinMax); psFree(tmpMaskVec); - psFree(tmpScalar); psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); return(1); @@ -1833,5 +1610,4 @@ psFree(tmpStatsMinMax); psFree(tmpMaskVec); - psFree(tmpScalar); psTrace(__func__, 4, "---- %s(0) end ----\n", __func__); @@ -1885,6 +1661,6 @@ binMedian = 0; } else { - tmpScalar->data.F32 = totalDataPoints/2.0; - binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + tmpScalar.data.F32 = totalDataPoints/2.0; + binMedian = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); if (binMedian < 0) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian); @@ -1892,5 +1668,4 @@ psFree(robustHistogram); psFree(cumulativeRobustHistogram); - psFree(tmpScalar); psFree(tmpMaskVec); psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); @@ -1905,4 +1680,5 @@ // XXX: Check return codes. // + printf("ADD: step 3: Interpolate to the exact 50-percent position: this is the robust histogram median.\n"); stats->robustMedian = fitQuadraticSearchForYThenReturnX( *(psVector* *)&cumulativeRobustHistogram->bounds, @@ -1918,20 +1694,20 @@ psS32 binLo = 0; psS32 binHi = 0; - tmpScalar->data.F32 = totalDataPoints * 0.158655f; - if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + tmpScalar.data.F32 = totalDataPoints * 0.158655f; + if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { binLo = 0; } else { // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. - binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); if (binLo > cumulativeRobustHistogram->nums->n-1) { binLo = cumulativeRobustHistogram->nums->n-1; } } - tmpScalar->data.F32 = totalDataPoints * 0.841345f; - if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + tmpScalar.data.F32 = totalDataPoints * 0.841345f; + if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { binHi = 0; } else { // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. - binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + binHi = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); if (binHi > cumulativeRobustHistogram->nums->n-1) { binHi = cumulativeRobustHistogram->nums->n-1; @@ -1947,5 +1723,4 @@ psFree(robustHistogram); psFree(cumulativeRobustHistogram); - psFree(tmpScalar); psFree(tmpMaskVec); psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); @@ -2014,17 +1789,6 @@ if (sigma < (2 * binSize)) { psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize); - psF32 medianLo; - psF32 medianHi; - if ((binMedian - 25) > 0) { - medianLo = robustHistogram->bounds->data.F32[binMedian - 25]; - } else { - medianLo = robustHistogram->bounds->data.F32[0]; - } - if ((binMedian + 25) < robustHistogram->bounds->n) { - medianHi = robustHistogram->bounds->data.F32[binMedian + 25]; - } else { - medianHi = robustHistogram->bounds->data.F32[robustHistogram->bounds->n - 1]; - } - + psF32 medianLo = robustHistogram->bounds->data.F32[PS_MAX(0, (binMedian - 25))]; + psF32 medianHi = robustHistogram->bounds->data.F32[PS_MIN(robustHistogram->bounds->n - 1, (binMedian + 25))]; psTrace(__func__, 6, "Masking data more than 25 bins from the median (%.2f)\n", stats->robustMedian); psTrace(__func__, 6, "The median is at bin number %d. We mask bins outside the bin range (%d:%d)\n", binMedian, binMedian - 25, binMedian + 25); @@ -2052,17 +1816,17 @@ psS32 binHi25 = 0; - tmpScalar->data.F32 = totalDataPoints * 0.25f; - if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + tmpScalar.data.F32 = totalDataPoints * 0.25f; + if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { // XXX: Special case where its in first bin. Must code last bin. binLo25 = 0; } else { - binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); - } - tmpScalar->data.F32 = totalDataPoints * 0.75f; - if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + binLo25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); + } + tmpScalar.data.F32 = totalDataPoints * 0.75f; + if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { // XXX: Special case where its in first bin. Must code last bin. binHi25 = 0; } else { - binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + binHi25 = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, &tmpScalar); } if ((binLo25 < 0) || (binHi25 < 0)) { @@ -2071,5 +1835,4 @@ psFree(robustHistogram); psFree(cumulativeRobustHistogram); - psFree(tmpScalar); psFree(tmpMaskVec); psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); @@ -2084,4 +1847,5 @@ // XXX: Check for errors. // + printf("ADD: step 8: Interpolate for the lower quartile positions.\n"); psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( *(psVector* *)&cumulativeRobustHistogram->bounds, @@ -2089,4 +1853,5 @@ binLo25, totalDataPoints * 0.25f); + printf("ADD: step 8: Interpolate for the upper quartile positions.\n"); psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( *(psVector* *)&cumulativeRobustHistogram->bounds, @@ -2133,5 +1898,4 @@ psFree(robustHistogram); psFree(cumulativeRobustHistogram); - psFree(tmpScalar); psFree(tmpMaskVec); psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); @@ -2170,16 +1934,16 @@ psS32 binMin = 0; psS32 binMax = 0; - tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma); - if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) { + tmpScalar.data.F32 = stats->robustMedian - (2.0 * sigma); + if (tmpScalar.data.F32 <= newHistogram->bounds->data.F32[0]) { binMin = 0; } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); - } - - tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma); - if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); + } + + tmpScalar.data.F32 = stats->robustMedian + (2.0 + sigma); + if (tmpScalar.data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { binMax = newHistogram->bounds->n-1; } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); } if ((binMin < 0) || (binMax < 0)) { @@ -2207,16 +1971,16 @@ // histogram median. // - tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma); - if (tmpScalar->data.F32 < tmpStatsMinMax->min) { + tmpScalar.data.F32 = stats->robustMedian - (20.0 * sigma); + if (tmpScalar.data.F32 < tmpStatsMinMax->min) { binMin = 0; } else { - binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); // XXX: check for errors here. } - tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma); - if (tmpScalar->data.F32 > tmpStatsMinMax->max) { + tmpScalar.data.F32 = stats->robustMedian + (20.0 * sigma); + if (tmpScalar.data.F32 > tmpStatsMinMax->max) { binMax = newHistogramSmoothed->n - 1; } else { - binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, &tmpScalar); // XXX: check for errors here. } @@ -2275,5 +2039,4 @@ psFree(robustHistogram); psFree(cumulativeRobustHistogram); - psFree(tmpScalar); psFree(newHistogram); psFree(x); @@ -2314,5 +2077,4 @@ psFree(tmpStatsMinMax); psFree(cumulativeRobustHistogram); - psFree(tmpScalar); psFree(tmpMaskVec); psFree(robustHistogram); @@ -2321,12 +2083,4 @@ return(0); } - - - - - - - - /*****************************************************************************/ @@ -2362,5 +2116,5 @@ newStruct->robustUQ = NAN; newStruct->robustLQ = NAN; - newStruct->robustN50 = -1; // XXX: This is never used + newStruct->robustN50 = -1; newStruct->fittedMean = NAN; newStruct->fittedStdev = NAN; @@ -2373,5 +2127,5 @@ newStruct->min = NAN; newStruct->max = NAN; - newStruct->binsize = NAN; // XXX: This is never used + newStruct->binsize = NAN; newStruct->nSubsample = 100000; newStruct->options = options;