Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 10475) +++ trunk/psLib/src/math/psStats.c (revision 10550) @@ -7,8 +7,5 @@ * on those data structures. * - * @author GLG, MHPCC - * - * XXX: The following stats members are never used, or set in this code. - * stats->clippedNvalues + * @author GLG (MHPCC), EAM (IfA) * * XXX: Must do @@ -16,8 +13,8 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.193 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-05 20:05:57 $ + * @version $Revision: 1.194 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-08 11:38:54 $ * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * Copyright 2006 IfA, University of Hawaii */ @@ -35,9 +32,10 @@ #include "psMemory.h" #include "psAbort.h" -#include "psImage.h" +//#include "psImage.h" #include "psVector.h" #include "psTrace.h" #include "psLogMsg.h" #include "psError.h" +#include "psHistogram.h" #include "psStats.h" #include "psMinimizeLMM.h" @@ -64,4 +62,5 @@ #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \ (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) + /*****************************************************************************/ /* TYPE DEFINITIONS */ @@ -87,147 +86,15 @@ NOTE: it is assumed that any call to these statistical functions will have been preceded by a call to the psVectorStats() function. Various sanity tests -will only be performed in psVectorStats(). Should we perform the sanity -checks in each routine anyway? +will only be performed in psVectorStats(). -XXX: For many of these private stats routines, what should be done if there -are no acceptable elements in the input vector (if no elements lie within -range, or there are no unmasked elements, or the input vector is NULL)? -Currently we set the value to NAN. - -XXX: Optimization: many routines have an "empty" boolean variable which keeps -track of whether or not the vector has any valid elements. This code can -possibly be optimized away. - *****************************************************************************/ -/****************************************************************************** -vectorSampleMean(myVector, maskVector, maskVal, stats): calculates the -mean of the input vector. If there was a problem with the mean calculation, -this routine sets stats->sampleMean to NAN. - *****************************************************************************/ -# if (0) - static bool vectorSampleMean(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats) -{ - psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - - psF32 mean = 0.0; // The mean - long count = 0; // Number of points contributing to this mean - - psF32 *data = myVector->data.F32; // Dereference - int numData = myVector->n; // Number of data points - - // If PS_STAT_USE_RANGE is requested, then we enter a slightly different loop. - if (!errors) { - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector) { - // No errors, check range, check mask - psU8 *maskData = maskVector->data.U8; - for (long i = 0; i < numData; i++) { - // Check if the data is with the specified range - if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) { - mean += data[i]; - count++; - } - } - } else { - // No errors, check range, no mask - for (long i = 0; i < numData; i++) { - if (stats->min <= data[i] && data[i] <= stats->max) { - mean += data[i]; - count++; - } - } - } - } else { - if (maskVector) { - // No errors, no range check, check mask - psU8 *maskData = maskVector->data.U8; - for (long i = 0; i < numData; i++) { - if (!(maskVal & maskData[i])) { - mean += data[i]; - count++; - } - } - - } else { - // No errors, no range check, no mask - for (long i = 0; i < numData; i++) { - mean += data[i]; - } - count = numData; - } - } - - if (count > 0) { - mean /= (psF32)count; - } else { - mean = NAN; - } - - } else { - psF32 sumWeights = 0.0; // The sum of the weights - psF32 *errorsData = errors->data.F32; - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector) { - psU8 *maskData = maskVector->data.U8; - for (long i = 0; i < numData; i++) { - // Check if the data is with the specified range - if (!(maskVal & maskData[i]) && stats->min <= data[i] && data[i] <= stats->max) { - float weight = 1.0 / PS_SQR(errorsData[i]); - mean += data[i] * weight; - sumWeights += weight; - count++; - } - } - } else { - for (long i = 0; i < myVector->n; i++) { - if (stats->min <= data[i] && data[i] <= stats->max) { - float weight = 1.0 / PS_SQR(errorsData[i]); - mean += data[i] * weight; - sumWeights += weight; - count++; - } - } - } - } else { - if (maskVector) { - psU8 *maskData = maskVector->data.U8; - for (long i = 0; i < myVector->n; i++) { - if (!(maskVal & maskData[i])) { - float weight = 1.0 / PS_SQR(errorsData[i]); - mean += data[i] * weight; - sumWeights += weight; - count++; - } - } - } else { - for (long i = 0; i < myVector->n; i++) { - float weight = 1.0 / PS_SQR(errorsData[i]); - mean += data[i] * weight; - sumWeights += weight; - } - count = myVector->n; - } - } - - if (count > 0) { - mean /= sumWeights; - } else { - mean = NAN; - } - } - - stats->sampleMean = mean; - if (isnan(mean)) { - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return false; - } - - psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); - return true; -} -#endif +For many of these private stats routines, it is possible that there are no acceptable elements +in the input vector (if no elements lie within range, or there are no unmasked elements, or the +input vector is NULL). In such cases, we set the desired value to NAN and return an error. +The calling function may clear the error. +*****************************************************************************/ + +// static prototypes: +static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords); +static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); /****************************************************************************** @@ -236,13 +103,17 @@ this routine sets stats->sampleMean to NAN. -XXX : using the method below, with a single loop for various options - costs only a small amount and is much easier to debug. running 10000 tests - of 1000 point vectors for the two methods gives: - separate single -(mask: 0, range: 0): 0.067 sec 0.073 sec -(mask: 1, range: 0): 0.098 sec 0.102 sec -(mask: 0, range: 0): 0.136 sec 0.162 sec -(mask: 1, range: 0): 0.177 sec 0.181 sec - *****************************************************************************/ +using the method below, with a single loop for various options costs only a small amount and is +much easier to debug. running 10000 tests of 1000 point vectors for the two methods gives: + + separate single w/errors +(mask: 0, range: 0): 0.067 sec 0.073 sec (10%) 0.072 sec +(mask: 1, range: 0): 0.098 sec 0.102 sec ( 4%) 0.119 sec +(mask: 0, range: 0): 0.136 sec 0.162 sec (20%) 0.170 sec +(mask: 1, range: 0): 0.177 sec 0.181 sec ( 2%) 0.198 sec + +(effect of errors not tested) + +To optmize this, use a macro and ifdef in or out the three states (errors, mask, range) +*****************************************************************************/ static bool vectorSampleMean(const psVector* myVector, const psVector* errors, @@ -253,6 +124,7 @@ psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); + long count = 0; // Number of points contributing to this mean psF32 mean = 0.0; // The mean - long count = 0; // Number of points contributing to this mean + psF32 weight; psF32 *data = myVector->data.F32; // Dereference @@ -262,45 +134,37 @@ bool useRange = stats->options & PS_STAT_USE_RANGE; - // If errors exist, we enter a slightly different loop. - // XXX collapse to a single loop? - if (!errors) { - for (long i = 0; i < numData; i++) { - // Check if the data is with the specified range - if (useRange && (data[i] < stats->min)) - continue; - if (useRange && (data[i] > stats->max)) - continue; - if (maskData && (maskData[i] & maskVal)) - continue; - mean += data[i]; - count++; - } - mean = (count > 0) ? mean / (psF32) count : NAN; - } else { - psF32 sumWeights = 0.0; // The sum of the weights - psF32 *errorsData = errors->data.F32; - for (long i = 0; i < numData; i++) { - // Check if the data is with the specified range - if (useRange && (data[i] < stats->min)) - continue; - if (useRange && (data[i] > stats->max)) - continue; - if (maskData && (maskData[i] & maskVal)) - continue; - - float weight = 1.0 / PS_SQR(errorsData[i]); + psF32 sumWeights = 0.0; // The sum of the weights + psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32; + + for (long i = 0; i < numData; i++) { + // Check if the data is with the specified range + if (useRange && (data[i] < stats->min)) + continue; + if (useRange && (data[i] > stats->max)) + continue; + if (maskData && (maskData[i] & maskVal)) + continue; + if (errors) { + weight = (errorsData[i] == 0) ? 0.0 : PS_SQR(errorsData[i]); mean += data[i] * weight; sumWeights += weight; - count++; - } + } else { + mean += data[i]; + } + count++; + } + if (errors) { mean = (count > 0) ? mean / sumWeights : NAN; - } - + } else { + mean = (count > 0) ? mean / count : NAN; + } stats->sampleMean = mean; if (isnan(mean)) { + // XXX raise an error here? psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); return false; } + stats->results |= PS_STAT_SAMPLE_MEAN; psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); return true; @@ -319,5 +183,5 @@ (mask: 0, range: 0): 0.179 sec 0.208 sec (mask: 1, range: 0): 0.200 sec 0.244 sec - *****************************************************************************/ +*****************************************************************************/ static long vectorMinMax(const psVector* myVector, const psVector* maskVector, @@ -350,4 +214,5 @@ } + // XXX save numValid in psStats? if (numValid == 0) { stats->max = NAN; @@ -356,70 +221,17 @@ stats->max = max; stats->min = min; + stats->results |= PS_STAT_MIN; + stats->results |= PS_STAT_MAX; } psTrace("psLib.math", 4, "---- %s(%d) end ----\n", __func__, numValid); return numValid; } - -#if 0 -/****************************************************************************** -vectorCheckNonEmpty(myVector, maskVector, maskVal, stats): This routine returns -"true" if the inputPsVector has 1 or more valid elements (a valid element is an -unmasked element within the specified min/max range). Otherwise, return -"false". - *****************************************************************************/ -static bool vectorCheckNonEmpty(const psVector* myVector, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats) -{ - psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(myVector, false); - PS_ASSERT_PTR_NON_NULL(stats, false); - - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector) { - for (long i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); - return true; - } - } - } else { - for (long i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); - return true; - } - } - } - } else { - if (maskVector != NULL) { - for (long i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); - return true; - } - } - } else { - if (myVector->n > 0) { - psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); - return true; - } - } - } - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return(false); -} -#endif /****************************************************************************** vectorSampleMedian(myVector, maskVector, maskVal, stats): calculates the median and quartiles of the input vector. Returns true on success (including if there -were no valid input vector elements). - *****************************************************************************/ -static bool vectorSampleMedian(const psVector* myVector, +were no valid input vector elements). Expects F32 vector for input. +*****************************************************************************/ +static bool vectorSampleMedian(const psVector* inVector, const psVector* maskVector, psMaskType maskVal, @@ -428,48 +240,34 @@ psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); + bool useRange = stats->options & PS_STAT_USE_RANGE; + psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; // Dereference the vector + psF32 *input = inVector->data.F32; // Dereference the vector + // Allocate temporary vectors for the data. - psVector* vector = psVectorAllocEmpty(myVector->n, PS_TYPE_F32); // Vector containing the unmasked values + psVector *outVector = psVectorAllocEmpty(inVector->n, PS_TYPE_F32); // Vector containing the unmasked values + psF32 *output = outVector->data.F32; // Dereference the vector + long count = 0; // Number of valid entries - // Determine if we must only use data points within a min/max range. - if (stats->options & PS_STAT_USE_RANGE) { - // Store all non-masked data points within the min/max range - // into the temporary vectors. - if (maskVector != NULL) { - for (long i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - vector->data.F32[count] = myVector->data.F32[i]; - count++; - } - } - } else { - for (long i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - vector->data.F32[count] = myVector->data.F32[i]; - count++; - } - } - } - } else { - // Store all non-masked data points into the temporary vectors. - if (maskVector) { - for (long i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - vector->data.F32[count] = myVector->data.F32[i]; - count++; - } - } - } else { - vector = psVectorCopy(vector, myVector, PS_TYPE_F32); - count = myVector->n; - } - } - vector->n = count; + // Store all non-masked data points within the min/max range + // into the temporary vectors. + for (long i = 0; i < inVector->n; i++) { + if (useRange && (input[i] < stats->min)) + continue; + if (useRange && (input[i] > stats->max)) + continue; + if (maskData && (maskData[i] & maskVal)) + continue; + + output[count] = input[i]; + count++; + } + outVector->n = count; + if (count == 0) { psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No valid data in input vector.\n"); - psFree(vector); + psFree(outVector); + stats->sampleUQ = NAN; + stats->sampleLQ = NAN; stats->sampleMedian = NAN; return false; @@ -477,133 +275,33 @@ // Sort the temporary vector. - vector = psVectorSort(vector, vector); // Sort in-place (since it's a copy, it's OK) - if (!vector) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - _("Failed to sort input data.")); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - psFree(vector); + // XXX this is messed up: should psVectorSort free on error or not? + if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK) + psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data.")); + stats->sampleUQ = NAN; + stats->sampleLQ = NAN; + stats->sampleMedian = NAN; return false; } - // Calculate the median exactly. Use the average if the number of samples - // is even. + // Calculate the median. Use the average if the number of samples if even. if (count % 2 == 0) { - stats->sampleMedian = 0.5 * (vector->data.F32[(count / 2) - 1] + - vector->data.F32[count / 2]); + stats->sampleMedian = 0.5 * (output[(count/2) - 1] + output[count/2]); } else { - stats->sampleMedian = vector->data.F32[count / 2]; + stats->sampleMedian = output[count/2]; } // Calculate the quartile points exactly. - stats->sampleUQ = vector->data.F32[3 * (count / 4)]; - stats->sampleLQ = vector->data.F32[count / 4]; + stats->sampleUQ = output[(int)(0.75*count)]; + stats->sampleLQ = output[(int)(0.25*count)]; + + stats->results |= PS_STAT_SAMPLE_MEDIAN; + stats->results |= PS_STAT_SAMPLE_QUARTILE; // Free the temporary data structures. - psFree(vector); + psFree(outVector); // Return "true" on success. psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); return true; -} - -/****************************************************************************** -vectorSmoothHistGaussian(): This routine smoothes the data in the input -robustHistogram with a Gaussian of width sigma. It returns a psVector of the -smoothed data. - *****************************************************************************/ -psVector *vectorSmoothHistGaussian(psHistogram *histogram, - psF32 sigma) -{ - psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - psTrace("psLib.math", 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma); - PS_ASSERT_PTR_NON_NULL(histogram, NULL); - PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); - PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - long numBins = histogram->nums->n; // Number of histogram bins - psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins - const psVector *bounds = histogram->bounds; // The bounds for the histogram bins - - if (!histogram->uniform) { - // - // We get here if the histogram is non-uniform. This code is not tested. - // However, it is also not used anywhere, yet. - // - psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n"); - - for (long i = 0; i < numBins; i++) { - // Determine the midpoint of bin i. - float iMid = PS_BIN_MIDPOINT(histogram, i); - - // - // We determine the bin numbers (jMin:jMax) corresponding to a - // range of data values surrounding iMid. The range is of size: - // 2*PS_GAUSS_WIDTH*sigma - long jMin, jMax; - psF32 x = iMid - (PS_GAUSS_WIDTH * sigma); - for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--) - ; - x = iMid + (PS_GAUSS_WIDTH * sigma); - for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++) - ; - - // - // Loop from jMin to jMax, computing the gaussian of data i. - // - smooth->data.F32[i] = 0.0; - for (long j = jMin ; j <= jMax ; j++) { - float jMid = PS_BIN_MIDPOINT(histogram, j); - smooth->data.F32[i] += - histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); - } - } - } else { - // - // We get here if the histogram is uniform. - // - for (long i = 0; i < numBins; i++) { - psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0]; - psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize); - - // - // XXX: The following is wrong. However, in practice, the sigma was too - // large compared to the binSize. This meant that the smoothing was done - // over 500 bins in the robust stats algorithm. This mean that the smoothing - // took way too long. - // - #if 0 - - gaussWidth = 10; - #endif - // - // We determine the bin numbers (jMin:jMax) corresponding to a - // range of data values surrounding iMid. The range is of size: - // 2*PS_GAUSS_WIDTH*sigma - // - psS32 jMin = PS_MAX(i - gaussWidth, 0); - psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1); - - // - // Loop from jMin to jMax, computing the gaussian of data i. - // - smooth->data.F32[i] = 0.0; - float iMid = PS_BIN_MIDPOINT(histogram, i); - for (long j = jMin ; j <= jMax ; j++) { - float jMid = PS_BIN_MIDPOINT(histogram, j); - smooth->data.F32[i] += - histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); - } - } - } - - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(smooth); - } - psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return(smooth); } @@ -619,5 +317,14 @@ NULL - *****************************************************************************/ +using the method below, with a single loop for various options costs only a small amount and is +much easier to debug. running 10000 tests of 1000 point vectors for the two methods gives: + + single +(mask: 0, range: 0): 0.193 sec +(mask: 1, range: 0): 0.257 sec +(mask: 0, range: 0): 0.349 sec +(mask: 1, range: 0): 0.401 sec + +*****************************************************************************/ static bool vectorSampleStdev(const psVector* myVector, const psVector* errors, @@ -630,7 +337,8 @@ // This procedure requires the mean. If it has not been already // calculated, then call vectorSampleMean() - if (isnan(stats->sampleMean)) { + if (!(stats->results & PS_STAT_SAMPLE_MEAN)) { vectorSampleMean(myVector, errors, maskVector, maskVal, stats); } + // If the mean is NAN, then generate a warning and set the stdev to NAN. if (isnan(stats->sampleMean)) { @@ -641,4 +349,9 @@ } + psF32 *data = myVector->data.F32; // Dereference + psU8 *maskData = (maskVector == NULL) ? NULL : maskVector->data.U8; + bool useRange = stats->options & PS_STAT_USE_RANGE; + psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32; + // Accumulate the sums double mean = stats->sampleMean; // The mean @@ -647,57 +360,19 @@ double errorDivisor = 0.0; // Division for errors long count = 0; // Number of data points being used - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector) { - for (long i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i]) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - double diff = myVector->data.F32[i] - mean; - sumSquares += PS_SQR(diff); - sumDiffs += diff; - count++; - if (errors) { - errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); - } - } - } - } else { - for (long i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && (myVector->data.F32[i] <= stats->max)) { - double diff = myVector->data.F32[i] - mean; - sumSquares += PS_SQR(diff); - sumDiffs += diff; - count++; - if (errors) { - errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); - } - } - } - count = myVector->n; - } - } else { - if (maskVector) { - for (long i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - double diff = myVector->data.F32[i] - mean; - sumSquares += PS_SQR(diff); - sumDiffs += diff; - count++; - if (errors) { - errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); - } - } - } - } else { - for (long i = 0; i < myVector->n; i++) { - double diff = myVector->data.F32[i] - mean; - sumSquares += PS_SQR(diff); - sumDiffs += diff; - count++; - if (errors) { - errorDivisor += 1.0 / PS_SQR(errors->data.F32[i]); - } - } - count = myVector->n; + for (long i = 0; i < myVector->n; i++) { + // Check if the data is with the specified range + if (useRange && (data[i] < stats->min)) + continue; + if (useRange && (data[i] > stats->max)) + continue; + if (maskData && (maskData[i] & maskVal)) + continue; + + double diff = data[i] - mean; + sumSquares += PS_SQR(diff); + sumDiffs += diff; + count++; + if (errors) { + errorDivisor += 1.0 / PS_SQR(errorsData[i]); } } @@ -717,7 +392,8 @@ stats->sampleStdev = (1.0 / sqrtf(errorDivisor)); } else { - stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / - (float)(count - 1)); - } + stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / (float)(count - 1)); + } + stats->results |= PS_STAT_SAMPLE_STDEV; + psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); @@ -736,5 +412,5 @@ Returns true for success; false otherwise - *****************************************************************************/ +*****************************************************************************/ static bool vectorClippedStats(const psVector* myVector, const psVector* errors, @@ -757,7 +433,8 @@ PS_CLIPPED_SIGMA_UB, -1); - // Allocate a psStats structure for calculating the mean, median, and - // stdev. - psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + // unless we succeed, these will have NAN values + stats->clippedMean = NAN; + stats->clippedStdev = NAN; + stats->clippedNvalues = 0; // We copy the mask vector, to preserve the original @@ -773,35 +450,29 @@ } - // 1. Compute the sample median. - vectorSampleMedian(myVector, tmpMask, maskVal, statsTmp); - if (isnan(statsTmp->sampleMedian)) { + // 1. Compute the sample median, which we save for output + vectorSampleMedian(myVector, tmpMask, maskVal, stats); + if (isnan(stats->sampleMedian)) { psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleMedian returned NAN\n"); - stats->clippedMean = NAN; - stats->clippedStdev = NAN; psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); psFree(tmpMask); - psFree(statsTmp); return false; } - psTrace("psLib.math", 6, "The initial sample median is %f\n", statsTmp->sampleMedian); - - // 2. Compute the sample standard deviation. - vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp); - if (isnan(statsTmp->sampleStdev)) { + psTrace("psLib.math", 6, "The initial sample median is %f\n", stats->sampleMedian); + + // 2. Compute the sample standard deviation, which we save for output + vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats); + if (isnan(stats->sampleStdev)) { psLogMsg(__func__, PS_LOG_WARN, "Call to vectorSampleStdev returned NAN\n"); - stats->clippedMean = NAN; - stats->clippedStdev = NAN; psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); psFree(tmpMask); - psFree(statsTmp); return false; } - psTrace("psLib.math", 6, "The initial sample stdev is %f\n", statsTmp->sampleStdev); + psTrace("psLib.math", 6, "The initial sample stdev is %f\n", stats->sampleStdev); // 3. Use the sample median as the first estimator of the mean X. - psF32 clippedMean = statsTmp->sampleMedian; + psF32 clippedMean = stats->sampleMedian; // 4. Use the sample stdev as the first estimator of the mean stdev. - psF32 clippedStdev = statsTmp->sampleStdev; + psF32 clippedStdev = stats->sampleStdev; // 5. Repeat N (stats->clipIter) times: @@ -836,4 +507,6 @@ // b) compute new mean and stdev + // Allocate a psStats structure for calculating the mean and stdev. + psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp); vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp); @@ -853,4 +526,5 @@ clippedStdev = statsTmp->sampleStdev; } + psFree(statsTmp); } @@ -859,375 +533,19 @@ // 7. The last calcuated value of x is the clipped mean. - if (stats->options & PS_STAT_CLIPPED_MEAN) { - stats->clippedMean = clippedMean; - psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean); - } // 8. The last calcuated value of stdev is the clipped stdev. - if (stats->options & PS_STAT_CLIPPED_STDEV) { - stats->clippedStdev = clippedStdev; - psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev); - } + // we always return both stats even if only one was requested + stats->clippedMean = clippedMean; + stats->clippedStdev = clippedStdev; + + stats->results |= PS_STAT_CLIPPED_MEAN; + stats->results |= PS_STAT_CLIPPED_STDEV; + + psTrace("psLib.math", 6, "The final clipped mean is %f\n", clippedMean); + psTrace("psLib.math", 6, "The final clipped stdev is %f\n", clippedStdev); psFree(tmpMask); - psFree(statsTmp); psTrace("psLib.math", 4, "---- %s(true) end ----\n", __func__); return true; } - -static psF32 QuadraticInverse(psF32 a, - psF32 b, - psF32 c, - psF32 y, - psF32 xLo, - psF32 xHi - ) -{ - 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; - - #if 0 - - 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); - #endif - - if (xLo <= x1 && x1 <= xHi) { - return x1; - } - if (xLo <= x2 && x2 <= xHi) { - return x2; - } - return 0.5 * (xLo + xHi); -} - - -/****************************************************************************** -fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general -routine which fits a quadratic to three points and returns the x-value -corresponding to the input y-value. This routine takes psVectors of x/y pairs -as input, and fits a quadratic to the 3 points surrounding element binNum in -the vectors (the midpoint between element i and i+1 is used for x[i]). It -then determines for what value x does that quadratic f(x) = yVal (the input -parameter). - -*****************************************************************************/ -static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal - ) -{ - psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - psTrace("psLib.math", 5, "binNum, yVal is (%d, %f)\n", binNum, yVal); - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(xVec); - PS_VECTOR_PRINT_F32(yVec); - } - - PS_ASSERT_VECTOR_NON_NULL(xVec, NAN); - PS_ASSERT_VECTOR_NON_NULL(yVec, NAN); - PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN); - PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN); - // PS_ASSERT_VECTORS_SIZE_EQUAL(xVec, yVec, NAN); - PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN); - PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN); - - psVector *x = psVectorAlloc(3, PS_TYPE_F64); - psVector *y = psVectorAlloc(3, PS_TYPE_F64); - psF32 tmpFloat = 0.0f; - - if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) { - // The general case. We have all three points. - x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); - x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); - x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); - y->data.F64[0] = yVec->data.F32[binNum - 1]; - y->data.F64[1] = yVec->data.F32[binNum]; - y->data.F64[2] = yVec->data.F32[binNum + 1]; - psTrace("psLib.math", 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]); - psTrace("psLib.math", 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); - psTrace("psLib.math", 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); - - // - // Ensure that the y value lies within range of the y values. - // - if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || - ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified yVal, %g, is not within y-range, %g to %g."), - (psF64)yVal, y->data.F64[0], y->data.F64[2]); - } - - // - // Ensure that the y values are monotonic. - // - if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || - ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { - psError(PS_ERR_UNKNOWN, true, - "This routine must be called with monotonically increasing or decreasing data points.\n"); - psFree(x); - psFree(y); - psTrace("psLib.math", 5, "---- %s() end ----\n", __func__); - return NAN; - } - - // Determine the coefficients of the polynomial. - psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x); - if (myPoly == NULL) { - psError(PS_ERR_UNEXPECTED_NULL, false, - _("Failed to fit a 1-dimensional polynomial to the three specified data points. Returning NAN.")); - psFree(x); - psFree(y); - psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__); - return NAN; - } - psTrace("psLib.math", 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]); - psTrace("psLib.math", 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]); - psTrace("psLib.math", 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]); - psTrace("psLib.math", 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])); - - psTrace("psLib.math", 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); - tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, - x->data.F64[0], x->data.F64[2]); - psFree(myPoly); - - if (isnan(tmpFloat)) { - psError(PS_ERR_UNEXPECTED_NULL, - false, _("Failed to determine the median of the fitted polynomial. Returning NAN.")); - psFree(x); - psFree(y); - psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - } else { - // These are special cases where the bin is at the beginning or end of the vector. - if (binNum == 0) { - // We have two points only at the beginning of the vectors x and y. - tmpFloat = 0.5 * (xVec->data.F32[binNum] + - xVec->data.F32[binNum + 1]); - } else if (binNum == (xVec->n - 1)) { - // The special case where we have two points only at the end of - // the vectors x and y. - // XXX: Is this right? - tmpFloat = xVec->data.F32[binNum]; - } else if (binNum == (xVec->n - 2)) { - // XXX: Is this right? - tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]); - } - } - - psTrace("psLib.math", 6, "FIT: return %f\n", tmpFloat); - psFree(x); - psFree(y); - - psTrace("psLib.math", 5, "---- %s(%f) end ----\n", __func__, tmpFloat); - return tmpFloat; -} - -/****************************************************************************** -NOTE: We assume unnormalized gaussians. - *****************************************************************************/ -static psF32 minimizeLMChi2Gauss1D(psVector *deriv, - const psVector *params, - const psVector *coords - ) -{ - psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN); - PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN); - PS_ASSERT_VECTOR_NON_NULL(coords, NAN); - PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN); - PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN); - - psF32 x = coords->data.F32[0]; - psF32 mean = params->data.F32[0]; - psF32 var = params->data.F32[1]; - psF32 dx = (x - mean); - - psF32 gauss = exp (-0.5*PS_SQR(dx)/var); - if (deriv) { - PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN); - PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN); - psF32 tmp = dx * gauss; - deriv->data.F32[0] = tmp / var; - deriv->data.F32[1] = tmp * dx / (var * var); - } - - - psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return gauss; -} - -/* - * vectorFittedStats requires guess for fittedMean and fittedStdev - * robustN50 should also be set - */ -static bool vectorFittedStats (const psVector* myVector, - const psVector* errors, - psVector* mask, - psMaskType maskVal, - psStats* stats) -{ - - psTrace("psLib.math", 6, "The guess mean is %f.\n", stats->fittedMean); - psTrace("psLib.math", 6, "The guess stdev is %f.\n", stats->fittedStdev); - - psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - - psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect - tmpScalar.type.type = PS_TYPE_F32; - - psF32 binSize = 1; - if (stats->options & PS_STAT_USE_BINSIZE) { - // Set initial bin size to the specified value. - binSize = stats->binsize; - psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { - // construct a histogram with (sigma/10 < binsize < sigma) - // set roughly so that the lowest bins have about 2 cnts - // Nsmallest ~ N50 / (4*dN)) - psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8)); - binSize = stats->fittedStdev / dN; - } - - // Determine the min/max of the vector (which prior outliers masked out) - int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values - float min = statsMinMax->min; - float max = statsMinMax->max; - if (numValid == 0 || isnan(min) || isnan(max)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return false; - } - - // Calculate the number of bins. - long numBins = (max - min) / binSize; - psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max); - psTrace("psLib.math", 6, "The new bin size is %f.\n", binSize); - psTrace("psLib.math", 6, "The numBins is %ld\n", numBins); - - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal); - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - long binMin = 0; // Low bin to check - long binMax = 0; // High bin to check - - // Fit a Gaussian to the bins in the requested range about the guess mean - // min,max overrides clipSigma - psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - maxFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->max)) { - maxFitSigma = fabs(stats->max); - } - - psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - minFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minFitSigma = fabs(stats->min); - } - - tmpScalar.data.F32 = stats->fittedMean - minFitSigma*stats->fittedStdev; - if (tmpScalar.data.F32 < min) { - binMin = 0; - } else { - binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); - } - tmpScalar.data.F32 = stats->fittedMean + maxFitSigma*stats->fittedStdev; - if (tmpScalar.data.F32 > max) { - binMax = histogram->bounds->n - 1; - } else { - binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar); - } - if (binMin < 0 || binMax < 0) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma); - psFree(statsMinMax); - psFree(histogram); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return false; - } - - // Generate the variables that will be used in the Gaussian fitting - // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins - psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates - psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates - for (long i = binMin, j = 0; i <= binMax ; i++, j++) { - y->data.F32[j] = histogram->nums->data.F32[i]; - psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value - ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i); - x->data[j] = ordinate; - } - if (psTraceGetLevel("psLib.math") >= 8) { - // XXX: Print the x array somehow. - PS_VECTOR_PRINT_F32(y); - } - psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n); - - // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) - // XXX EAM : why not just fit the normalization as well?? - p_psNormalizeVectorRange(y, 0.0, 1.0); - - // Fit a Gaussian to the data. - psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information - psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian - // Initial guess for the mean (index 0) and var (index 1). - params->data.F32[0] = stats->fittedMean; - params->data.F32[1] = PS_SQR(stats->fittedStdev); - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(params); - PS_VECTOR_PRINT_F32(y); - } - if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(statsMinMax); - psFree(histogram); - psFree(x); - psFree(y); - psFree(minimizer); - psFree(params); - psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); - return false; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(params); - } - - // The fitted mean is the Gaussian mean. - stats->fittedMean = params->data.F32[0]; - psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean); - - // The fitted standard deviation - stats->fittedStdev = sqrt(params->data.F32[1]); - psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev); - - // Clean up after fitting - psFree (histogram); - psFree (x); - psFree (y); - psFree (minimizer); - psFree (params); - psFree (statsMinMax); - return true; -} - /****************************************************************************** @@ -1278,4 +596,5 @@ } + // statsMinMax is only applied to a subset of the data points psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max psHistogram *histogram = NULL; // Histogram of the data @@ -1287,7 +606,5 @@ // Iterate to get the best bin size for (int iterate = 1; iterate > 0; iterate++) { - psTrace("psLib.math", 6, - "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", - iterate); + psTrace("psLib.math", 6, "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", iterate); // Get the minimum and maximum values @@ -1308,11 +625,7 @@ if (fabs(max - min) <= FLT_EPSILON) { psTrace("psLib.math", 5, "All data points have the same value: %f.\n", min); - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = min; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = min; - stats->robustLQ = min; - } + stats->robustMedian = min; + stats->robustUQ = min; + stats->robustLQ = min; stats->robustN50 = numValid; psFree(statsMinMax); @@ -1433,21 +746,6 @@ return false; } - // XXXX we can probably just use the bin values to get sigma without interpolating... - // ADD step 4b: Interpolate Sigma to find these two positions exactly: these are the 1sigma positions. - #if 0 - // XXX: I am overriding the ADD for now. The following code implements the ADD exactly and is having - // problems fitting a 2nd-order polynomial to data y-values suchs as (1, 1, 100). Therefore, I - // commented out the code and am doing simply linear interpolation. - - // Value for the 15.8655% mark - float binLoF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binLo, - totalDataPoints * 0.158655f); - // Value for the 84.1345% mark - float binHiF32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, binHi, - totalDataPoints * 0.841345f); - psTrace("psLib.math", 6, "The exact 15.8655%% and 84.1345%% data point positions are: (%f, %f)\n", - binLoF32, binHiF32); - #else - // This code basically interpolates to find the positions exactly. + + // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma positions. psTrace("psLib.math", 6, "binLo is %ld. Nums at that bin and the next are (%.2f, %.2f)\n", binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]); @@ -1455,7 +753,9 @@ binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]); - float deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo]; - float deltaNums; - float prevPixels; + float deltaBounds, deltaNums, prevPixels; + float percentNums, base, binLoF32; + + // find the -1 sigma points with linear interpolation + deltaBounds = cumulative->bounds->data.F32[binLo+1] - cumulative->bounds->data.F32[binLo]; if (binLo == 0) { deltaNums = cumulative->nums->data.F32[0]; @@ -1465,11 +765,12 @@ prevPixels = cumulative->nums->data.F32[binLo-1]; } - float percentNums = (totalDataPoints * 0.158655f) - prevPixels; - float base = cumulative->bounds->data.F32[binLo]; - float binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark + percentNums = (totalDataPoints * 0.158655f) - prevPixels; + base = cumulative->bounds->data.F32[binLo]; + binLoF32 = base + (deltaBounds / deltaNums) * percentNums; // Value for the 15.8655% mark psTrace("psLib.math", 6, "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", base, deltaBounds, deltaNums, prevPixels, percentNums); + // find the +1 sigma points with linear interpolation deltaBounds = cumulative->bounds->data.F32[binHi+1] - cumulative->bounds->data.F32[binHi]; if (binHi == 0) { @@ -1486,8 +787,9 @@ "(base, deltaBounds, deltaNums, prevPixels, percentNums) is (%.2f %.2f %.2f %.2f %.2f)\n", base, deltaBounds, deltaNums, prevPixels, percentNums); - psTrace("psLib.math", 6, + + // report +/- 1 sigma points + psTrace("psLib.math", 5, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32); - #endif // ADD step 5: Determine SIGMA as 1/2 of the distance between these positions. @@ -1585,35 +887,507 @@ psTrace("psLib.math", 6, "The robustN50 is %ld.\n", N50); - // Fitted statistics - if (stats->options & PS_STAT_FITTED_MEAN || stats->options & PS_STAT_FITTED_STDEV) { - stats->fittedStdev = stats->robustStdev; // pass the guess sigma - stats->fittedMean = stats->robustMedian; // pass the guess mean - if (!vectorFittedStats(myVector, errors, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 1)"); - psFree(statsMinMax); - psFree(mask); - - return false; - } - // if there is a large swing in sigma, try a second time - if ((stats->fittedStdev/stats->robustStdev) < 0.75) { - if (!vectorFittedStats (myVector, errors, mask, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "Unable to estimate statistics (pass 2)"); - psFree(statsMinMax); - psFree(mask); - - return false; - } - } - } - // Clean up psFree(statsMinMax); psFree(mask); + stats->results |= PS_STAT_ROBUST_MEDIAN; + stats->results |= PS_STAT_ROBUST_STDEV; + psTrace("psLib.math", 4, "---- %s(0) end ----\n", __func__); return true; } +/* + * vectorFittedStats requires guess for fittedMean and fittedStdev + * robustN50 should also be set + */ +static bool vectorFittedStats (const psVector* myVector, + const psVector* errors, + psVector* mask, + psMaskType maskVal, + psStats* stats) +{ + + // This procedure requires the mean. If it has not been already + // calculated, then call vectorSampleMean() + if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) { + vectorRobustStats(myVector, errors, mask, maskVal, stats); + } + + // If the mean is NAN, then generate a warning and set the stdev to NAN. + if (isnan(stats->robustMedian)) { + stats->fittedStdev = NAN; + stats->fittedStdev = NAN; + psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); + return false; + } + + float guessStdev = stats->robustStdev; // pass the guess sigma + float guessMean = stats->robustMedian; // pass the guess mean + + psTrace("psLib.math", 6, "The guess mean is %f.\n", guessMean); + psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev); + + bool done = false; + for (int iteration = 0; !done && (iteration < 2); iteration ++) { + psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max + + psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect + tmpScalar.type.type = PS_TYPE_F32; + + psF32 binSize = 1; + if (stats->options & PS_STAT_USE_BINSIZE) { + // Set initial bin size to the specified value. + binSize = stats->binsize; + psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize); + } else { + // construct a histogram with (sigma/10 < binsize < sigma) + // set roughly so that the lowest bins have about 2 cnts + // Nsmallest ~ N50 / (4*dN)) + psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8)); + binSize = guessStdev / dN; + } + + // Determine the min/max of the vector (which prior outliers masked out) + int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values + float min = statsMinMax->min; + float max = statsMinMax->max; + if (numValid == 0 || isnan(min) || isnan(max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(statsMinMax); + psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // Calculate the number of bins. + // XXX can we calculate the binMin, binMax **before** building this histogram? + long numBins = (max - min) / binSize; + psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max); + psTrace("psLib.math", 6, "The new bin size is %f.\n", binSize); + psTrace("psLib.math", 6, "The numBins is %ld\n", numBins); + + psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) + histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal); + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(histogram->nums); + } + + long binMin = 0; // Low bin to check + long binMax = 0; // High bin to check + + // Fit a Gaussian to the bins in the requested range about the guess mean + // min,max overrides clipSigma + psF32 maxFitSigma = 2.0; + if (isfinite(stats->clipSigma)) { + maxFitSigma = fabs(stats->clipSigma); + } + if (isfinite(stats->max)) { + maxFitSigma = fabs(stats->max); + } + + psF32 minFitSigma = 2.0; + if (isfinite(stats->clipSigma)) { + minFitSigma = fabs(stats->clipSigma); + } + if (isfinite(stats->min)) { + minFitSigma = fabs(stats->min); + } + + tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev; + if (tmpScalar.data.F32 < min) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); + } + tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev; + if (tmpScalar.data.F32 > max) { + binMax = histogram->bounds->n - 1; + } else { + binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar); + } + if (binMin < 0 || binMax < 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma); + psFree(statsMinMax); + psFree(histogram); + psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // Generate the variables that will be used in the Gaussian fitting + // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins + psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates + psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates + for (long i = binMin, j = 0; i <= binMax ; i++, j++) { + y->data.F32[j] = histogram->nums->data.F32[i]; + psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value + ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i); + x->data[j] = ordinate; + } + if (psTraceGetLevel("psLib.math") >= 8) { + // XXX: Print the x array somehow. + PS_VECTOR_PRINT_F32(y); + } + psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n); + psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); + psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax); + + // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) + // XXX EAM : why not just fit the normalization as well?? + p_psNormalizeVectorRange(y, 0.0, 1.0); + + // Fit a Gaussian to the data. + psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information + psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian + // Initial guess for the mean (index 0) and var (index 1). + params->data.F32[0] = guessMean; + params->data.F32[1] = PS_SQR(guessStdev); + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(params); + PS_VECTOR_PRINT_F32(y); + } + if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(params); + psFree(minimizer); + psFree(x); + psFree(y); + psFree(histogram); + psFree(statsMinMax); + psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); + return false; + } + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(params); + } + + guessMean = params->data.F32[0]; + guessStdev = sqrt(params->data.F32[1]); + if (guessStdev > 0.75*stats->robustStdev) { + done = true; + } + + // Clean up after fitting + psFree (params); + psFree (minimizer); + psFree (x); + psFree (y); + psFree (histogram); + psFree (statsMinMax); + } + + // The fitted mean is the Gaussian mean. + stats->fittedMean = guessMean; + psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean); + + // The fitted standard deviation + stats->fittedStdev = guessStdev; + psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev); + + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; + + return true; +} + + +/******************** + * vectorFittedStats_v2 requires guess for fittedMean and fittedStdev + * robustN50 should also be set + * gaussian fit is performed using 2D polynomial to ln(y) + ********************/ +static bool vectorFittedStats_v2 (const psVector* myVector, + const psVector* errors, + psVector* mask, + psMaskType maskVal, + psStats* stats) +{ + + // This procedure requires the mean. If it has not been already + // calculated, then call vectorSampleMean() + if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) { + vectorRobustStats(myVector, errors, mask, maskVal, stats); + } + + // If the mean is NAN, then generate a warning and set the stdev to NAN. + if (isnan(stats->robustMedian)) { + stats->fittedStdev = NAN; + stats->fittedStdev = NAN; + psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); + return false; + } + + float guessStdev = stats->robustStdev; // pass the guess sigma + float guessMean = stats->robustMedian; // pass the guess mean + + psTrace("psLib.math", 6, "The guess mean is %f.\n", guessMean); + psTrace("psLib.math", 6, "The guess stdev is %f.\n", guessStdev); + + bool done = false; + for (int iteration = 0; !done && (iteration < 2); iteration ++) { + psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max + + psScalar tmpScalar; // Temporary scalar variable, for p_psVectorBinDisect + tmpScalar.type.type = PS_TYPE_F32; + + psF32 binSize = 1; + if (stats->options & PS_STAT_USE_BINSIZE) { + // Set initial bin size to the specified value. + binSize = stats->binsize; + psTrace("psLib.math", 6, "Setting initial robust bin size to %.2f\n", binSize); + } else { + // construct a histogram with (sigma/10 < binsize < sigma) + // set roughly so that the lowest bins have about 2 cnts + // Nsmallest ~ N50 / (4*dN)) + psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8)); + binSize = guessStdev / dN; + } + + // Determine the min/max of the vector (which prior outliers masked out) + int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values + float min = statsMinMax->min; + float max = statsMinMax->max; + if (numValid == 0 || isnan(min) || isnan(max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(statsMinMax); + psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // Calculate the number of bins. + // XXX can we calculate the binMin, binMax **before** building this histogram? + long numBins = (max - min) / binSize; + psTrace("psLib.math", 6, "The new min/max values are (%f, %f).\n", min, max); + psTrace("psLib.math", 6, "The new bin size is %f.\n", binSize); + psTrace("psLib.math", 6, "The numBins is %ld\n", numBins); + + psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) + histogram = psVectorHistogram(histogram, myVector, errors, mask, maskVal); + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(histogram->nums); + } + + long binMin = 0; // Low bin to check + long binMax = 0; // High bin to check + + // Fit a Gaussian to the bins in the requested range about the guess mean + // min,max overrides clipSigma + psF32 maxFitSigma = 2.0; + if (isfinite(stats->clipSigma)) { + maxFitSigma = fabs(stats->clipSigma); + } + if (isfinite(stats->max)) { + maxFitSigma = fabs(stats->max); + } + + psF32 minFitSigma = 2.0; + if (isfinite(stats->clipSigma)) { + minFitSigma = fabs(stats->clipSigma); + } + if (isfinite(stats->min)) { + minFitSigma = fabs(stats->min); + } + + tmpScalar.data.F32 = guessMean - minFitSigma*guessStdev; + if (tmpScalar.data.F32 < min) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect(histogram->bounds, &tmpScalar); + } + tmpScalar.data.F32 = guessMean + maxFitSigma*guessStdev; + if (tmpScalar.data.F32 > max) { + binMax = histogram->bounds->n - 1; + } else { + binMax = p_psVectorBinDisect(histogram->bounds, &tmpScalar); + } + if (binMin < 0 || binMax < 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the -%f sigma : +%f sigma bins\n", minFitSigma, maxFitSigma); + psFree(statsMinMax); + psFree(histogram); + psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); + return false; + } + + // Generate the variables that will be used in the Gaussian fitting + // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins + psVector *y = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates + psVector *x = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of ordinates + long j = 0; + for (long i = binMin; i <= binMax ; i++) { + if (histogram->nums->data.F32[i] <= 0.0) + continue; + x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); + y->data.F32[j] = log(histogram->nums->data.F32[i]); + // XXX note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) + j++; + } + y->n = x->n = j; + if (psTraceGetLevel("psLib.math") >= 8) { + // XXX: Print the x array somehow. + PS_VECTOR_PRINT_F32(y); + } + psTrace("psLib.math", 6, "The clipped numBins is %ld\n", y->n); + psTrace("psLib.math", 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); + psTrace("psLib.math", 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax); + + // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + + // XXX how can we protect against some extreme outliers? the robust histogram + // should have already dealt with those, but we are again sensitive to them... + // psStats *fitStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + // fitStats->clipIter = 3.0; + // fitStats->clipSigma = 3.0; + // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_U8); + // psVectorInit (fitMask, 0); + + if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(x); + psFree(y); + psFree(poly); + // psFree(fitStats); + // psFree(fitMask); + psFree(histogram); + psFree(statsMinMax); + psTrace("psLib.math", 4, "---- %s(false) end ----\n", __func__); + return false; + } + + guessStdev = sqrt(-0.5/poly->coeff[2]); + guessMean = poly->coeff[1]*PS_SQR(guessStdev); + if (guessStdev > 0.75*stats->robustStdev) { + done = true; + } else { + psTrace("psLib.math", 6, "The new guess mean is %f.\n", guessMean); + psTrace("psLib.math", 6, "The new guess stdev is %f.\n", guessStdev); + } + + // Clean up after fitting + psFree (x); + psFree (y); + psFree (poly); + // psFree (fitStats); + // psFree (fitMask); + psFree (histogram); + psFree (statsMinMax); + } + + // The fitted mean is the Gaussian mean. + stats->fittedMean = guessMean; + psTrace("psLib.math", 6, "The fitted mean is %f.\n", stats->fittedMean); + + // The fitted standard deviation + stats->fittedStdev = guessStdev; + psTrace("psLib.math", 6, "The fitted stdev is %f.\n", stats->fittedStdev); + + stats->results |= PS_STAT_FITTED_MEAN_V2; + stats->results |= PS_STAT_FITTED_STDEV_V2; + + return true; +} + + +/****************************************************************************** +vectorSmoothHistGaussian(): This routine smoothes the data in the input +robustHistogram with a Gaussian of width sigma. It returns a psVector of the +smoothed data. + +XXX this function is unused +*****************************************************************************/ +psVector *vectorSmoothHistGaussian(psHistogram *histogram, + psF32 sigma) +{ + psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); + psTrace("psLib.math", 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma); + PS_ASSERT_PTR_NON_NULL(histogram, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); + PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(histogram->nums); + } + + long numBins = histogram->nums->n; // Number of histogram bins + psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins + const psVector *bounds = histogram->bounds; // The bounds for the histogram bins + + if (!histogram->uniform) { + // + // We get here if the histogram is non-uniform. This code is not tested. + // However, it is also not used anywhere, yet. + // + psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSmoothHistGaussian() on non-uniform histograms has not been tested or used.\n"); + + for (long i = 0; i < numBins; i++) { + // Determine the midpoint of bin i. + float iMid = PS_BIN_MIDPOINT(histogram, i); + + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + long jMin, jMax; + psF32 x = iMid - (PS_GAUSS_WIDTH * sigma); + for (jMin = i; jMin > 0 && bounds->data.F32[jMin] > x; jMin--) + ; + x = iMid + (PS_GAUSS_WIDTH * sigma); + for (jMax = i; jMax < bounds->n - 1 && bounds->data.F32[jMax + 1] > x; jMax++) + ; + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + for (long j = jMin ; j <= jMax ; j++) { + float jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } else { + // + // We get here if the histogram is uniform. + // + for (long i = 0; i < numBins; i++) { + psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0]; + psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize); + + // + // XXX: The following is wrong. However, in practice, the sigma was too + // large compared to the binSize. This meant that the smoothing was done + // over 500 bins in the robust stats algorithm. This mean that the smoothing + // took way too long. + // + #if 0 + + gaussWidth = 10; + #endif + // + // We determine the bin numbers (jMin:jMax) corresponding to a + // range of data values surrounding iMid. The range is of size: + // 2*PS_GAUSS_WIDTH*sigma + // + psS32 jMin = PS_MAX(i - gaussWidth, 0); + psS32 jMax = PS_MIN(i + gaussWidth, bounds->n - 1); + + // + // Loop from jMin to jMax, computing the gaussian of data i. + // + smooth->data.F32[i] = 0.0; + float iMid = PS_BIN_MIDPOINT(histogram, i); + for (long j = jMin ; j <= jMax ; j++) { + float jMid = PS_BIN_MIDPOINT(histogram, j); + smooth->data.F32[i] += + histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); + } + } + } + + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(smooth); + } + psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); + return(smooth); +} + /*****************************************************************************/ @@ -1621,6 +1395,4 @@ /*****************************************************************************/ - -static void histogramFree(psHistogram* myHist); // We keep statsFree so that we can identify statistics pointers from the memblock @@ -1632,5 +1404,5 @@ /****************************************************************************** psStatsAlloc(): This routine must create a new psStats data structure. - *****************************************************************************/ +*****************************************************************************/ psStats* psStatsAlloc(psStatsOptions options) { @@ -1675,299 +1447,4 @@ } - - -/****************************************************************************** -psHistogramAlloc(lower, upper, n): allocate a uniform histogram structure -with the specifed upper and lower limits, and the specifed number of bins. -This routine will also set the bounds for each of the bins. - -Input: - lower - upper - n -Returns: - The histogram structure - *****************************************************************************/ -psHistogram* psHistogramAlloc(float lower, float upper, int n) -{ - psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - psTrace("psLib.math", 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n); - PS_ASSERT_INT_POSITIVE(n, NULL); - PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL); - - // Allocate memory for the new histogram structure. If there are N bins, then there are N+1 bounds to - // those bins. - psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure - psMemSetDeallocator(newHist, (psFreeFunc) histogramFree); - psVector* newBounds = psVectorAlloc(n + 1, PS_TYPE_F32); - newHist->bounds = newBounds; - - // Calculate the bounds for each bin. - psF32 binSize = (upper - lower) / (psF32)n; // The histogram bin size - // XXX: Is the following necessary? It prevents the max data point from being in a non-existant bin. - binSize += FLT_EPSILON; - for (long i = 0; i < n + 1; i++) { - newBounds->data.F32[i] = lower + (binSize * (psF32)i); - } - - // Allocate the bins, and initialize them to zero. - newHist->nums = psVectorAlloc(n, PS_TYPE_F32); - psVectorInit(newHist->nums, 0.0); - - // Initialize the other members. - newHist->minNum = 0; - newHist->maxNum = 0; - newHist->uniform = true; - - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return newHist; -} - -/****************************************************************************** -psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure -with the specifed bounds. - -Input: - bounds -Returns: - The histogram structure - *****************************************************************************/ -psHistogram* psHistogramAllocGeneric(const psVector* bounds) -{ - psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(bounds, NULL); - PS_ASSERT_VECTOR_TYPE(bounds, PS_TYPE_F32, NULL); - PS_ASSERT_LONG_LARGER_THAN_OR_EQUAL(bounds->n, (long)2, NULL); - - // Allocate memory for the new histogram structure. - psHistogram *newHist = (psHistogram* ) psAlloc(sizeof(psHistogram)); // The new histogram structure - psMemSetDeallocator(newHist, (psFreeFunc) histogramFree); - psVector* newBounds = psVectorCopy(NULL, bounds, PS_TYPE_F32); - newHist->bounds = newBounds; - - // Allocate the bins, and initialize them to zero. If there are N bounds, - // then there are N-1 bins. - newHist->nums = psVectorAlloc((bounds->n) - 1, PS_TYPE_F32); - psVectorInit(newHist->nums, 0.0); - - // Initialize the other members. - newHist->minNum = 0; - newHist->maxNum = 0; - newHist->uniform = false; - - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return (newHist); -} - -static void histogramFree(psHistogram* myHist) -{ - psFree(myHist->bounds); - psFree(myHist->nums); -} - - -bool psMemCheckHistogram(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc)histogramFree ); -} - - - -/***************************************************************************** -UpdateHistogramBins(binNum, out, data, error): This routine is to be used when -updating the histogram in the presence of errors in the input data. We treat -the data point as a boxcar PDF and update a range of points surrounding the -histogram bin which contains the point. The width of that boxcar is defined -as 2.35 * error. - -XXX: Must test this. - *****************************************************************************/ -static bool UpdateHistogramBins(long binNum, // Bin number of the data point - psHistogram* out, // The histogram to be updated - psF32 data, // The data point value - psF32 error // The error in the data point - ) -{ - psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(out, false); - PS_ASSERT_PTR_NON_NULL(out->bounds, false); - PS_ASSERT_PTR_NON_NULL(out->nums, false); - PS_ASSERT_LONG_WITHIN_RANGE(binNum, (long)0, (long)((out->nums->n)-1), false); - PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(error, 0.0, false); - PS_ASSERT_FLOAT_WITHIN_RANGE(data, out->bounds->data.F32[0], - out->bounds->data.F32[(out->bounds->n)-1], false); - - psF32 boxcarWidth = 2.35 * error; // Width of the boxcar - psF32 boxcarCenter = (out->bounds->data.F32[binNum] + - out->bounds->data.F32[binNum+1]) / 2.0; // Centre of the boxcar - psF32 boxcarLeft = boxcarCenter - (boxcarWidth / 2.0); // Left endpoint of the boxcar for the PDF - psF32 boxcarRight = boxcarCenter + (boxcarWidth / 2.0); // Right endpoint of the boxcar for the PDF - psS32 boxcarLeftBinNum = 0; // Bin number for left endpoint - psS32 boxcarRightBinNum = 0; // Bin number for right endpoint - - // Determine the left endpoint of the boxcar for the PDF. - for (long bin = binNum; bin >= 0; bin--) { - if (out->nums->data.F32[bin] <= boxcarLeft) { - boxcarLeftBinNum = bin; - break; - } - } - - // Determine the right endpoint of the boxcar for the PDF. - for (long bin = binNum; bin < out->nums->n; bin++) { - if (out->nums->data.F32[bin] >= boxcarRight) { - boxcarRightBinNum = bin; - break; - } - } - - // If the boxcar fits entirely inside this bin, then simply add 1.0 to the - // bin and return. - if (boxcarLeftBinNum == boxcarRightBinNum) { - out->nums->data.F32[binNum] += 1.0; - psTrace("psLib.math", 3, "---- %s(true) end ----\n", __func__); - return true; - } - - // If we get here, multiple bins must be updated. We handle the left-most endpoint, and right-most - // endpoints differently. - out->nums->data.F32[boxcarLeftBinNum] += - (out->bounds->data.F32[boxcarLeftBinNum+1] - boxcarLeft) / boxcarWidth; - - // Loop through the center bins, if any. - for (long bin = boxcarLeftBinNum + 1; bin < (boxcarRightBinNum - 1); bin++) { - out->nums->data.F32[bin] += - (out->bounds->data.F32[bin+1] - out->bounds->data.F32[bin]) / boxcarWidth; - } - - // Handle the right endpoint differently. - out->nums->data.F32[boxcarRightBinNum]+= - (boxcarRight - out->bounds->data.F32[boxcarRightBinNum]) / boxcarWidth; - - psTrace("psLib.math", 3, "---- %s(true) end ----\n", __func__); - return true; -} - - -/***************************************************************************** -psVectorHistogram(out, in, errors, mask, maskVal): this procedure takes as -input a preallocated and initialized histogram structure. It fills the bins -in that histogram structure in accordance with the input data "in" and the, -possibly NULL, mask vector. - -Inputs: - out - in - mask - maskVal -Returns: - The histogram structure "out". - *****************************************************************************/ -psHistogram* psVectorHistogram(psHistogram* out, - const psVector* values, - const psVector* errors, - const psVector* mask, - psMaskType maskVal) -{ - psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_PTR_NON_NULL(out, NULL); - PS_ASSERT_VECTOR_NON_NULL(out->bounds, NULL); - PS_ASSERT_VECTOR_TYPE(out->bounds, PS_TYPE_F32, NULL); - PS_ASSERT_INT_NONNEGATIVE(out->bounds->n, NULL); - PS_ASSERT_VECTOR_NON_NULL(out->nums, NULL); - PS_ASSERT_VECTOR_TYPE(out->nums, PS_TYPE_F32, NULL); - PS_ASSERT_INT_NONNEGATIVE(out->nums->n, NULL); - PS_ASSERT_VECTOR_NON_NULL(values, out); - if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(values, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - if (errors) { - PS_ASSERT_VECTORS_SIZE_EQUAL(values, errors, NULL); - PS_ASSERT_VECTOR_TYPE(errors, values->type.type, NULL); - } - - long binNum = 0; // A temporary bin number - long numBins = out->nums->n; // The total number of bins - psScalar tmpScalar; - tmpScalar.type.type = PS_TYPE_F32; - - // Convert input and errors vectors to F32 if necessary. - psVector* inF32 = NULL; // F32 version of input vector - if (values->type.type == PS_TYPE_F32) { - inF32 = psMemIncrRefCounter((psPtr)values); - } else { - inF32 = psVectorCopy(NULL, values, PS_TYPE_F32); - } - psVector* errorsF32 = NULL; // F32 version of errors vector - if (errors) { - if (errors->type.type == PS_TYPE_F32) { - errorsF32 = psMemIncrRefCounter((psPtr)errors); - } else { - errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32); - } - } - - for (long i = 0; i < inF32->n; i++) { - // Check if this pixel is masked, and if so, skip it. - if (!mask || (mask && (!(mask->data.U8[i] & maskVal)))) { - if (inF32->data.F32[i] < out->bounds->data.F32[0]) { - // If this pixel is below minimum value, count it, then skip. - out->minNum++; - } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) { - // If this pixel is above maximum value, count it, then skip. - out->maxNum++; - } else { - // If this is a uniform histogram, determining the correct - // number is trivial. - if (out->uniform == true) { - float binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; // Histogram bin size - binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize; - if (errorsF32) { - if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errorsF32->data.F32[i])) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram " - "bins with the errors vector.\n"); - } - } else { - // This if-statement really shouldn't be necessary. - // However, due to numerical lack of precision, we - // occasionally produce a binNum outside the range. - if (binNum >= out->nums->n) { - binNum = out->nums->n - 1; - } - (out->nums->data.F32[binNum])+= 1.0; - } - - } else { - // If this is a non-uniform histogram, determining the - // correct bin number requires a bit more work. - tmpScalar.data.F32 = inF32->data.F32[i]; - binNum = p_psVectorBinDisect( *(psVector* *)&out->bounds, &tmpScalar); - if (binNum < 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psVectorHistogram(): element outside histogram bounds.\n"); - } else { - if (errorsF32 != NULL) { - if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram " - "bins with the errors vector.\n"); - } - } else { - out->nums->data.F32[binNum] += 1.0; - } - } - } - } - } - } - - psFree(inF32); - psFree(errorsF32); - - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return (out); -} - /****************************************************************************** psVectorStats(in, mask, maskVal, stats): this is the public API @@ -1984,11 +1461,10 @@ XXX: Should we free stats if the asserts fail? NO; we don't own it (RHL). - *****************************************************************************/ -psStats* psVectorStats( - psStats* stats, - const psVector* in, - const psVector* errors, - const psVector* mask, - psMaskType maskVal) +*****************************************************************************/ +bool psVectorStats(psStats* stats, + const psVector* in, + const psVector* errors, + const psVector* mask, + psMaskType maskVal) { psTrace("psLib.math", 3,"---- %s() begin ----\n", __func__); @@ -2036,9 +1512,14 @@ } + // init the value of stats->results: this is used internally to check if + // prior functions have been called + stats->results = PS_STAT_NONE; + bool status = true; + // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEAN) { if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n"); - stats->sampleMean = NAN; + psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean"); + status &= false; } } @@ -2047,8 +1528,6 @@ if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) { if (!vectorSampleMedian(inF32, maskU8, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMedian() returned an error.\n"); - stats->sampleMedian = NAN; - stats->sampleUQ = NAN; - stats->sampleLQ = NAN; + psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median"); + status &= false; } } @@ -2056,30 +1535,7 @@ // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_STDEV) { - if (!vectorSampleMean(inF32, errorsF32, maskU8, maskVal, stats)) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: vectorSampleMean() returned an error.\n"); - stats->sampleMean = NAN; - } else { - vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats); - } - } - - // ************************************************************************ - if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE | - PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, _("Failed to calculate the specified statistic.")); - psFree(inF32); - psFree(errorsF32); - psTrace("psLib.math", 3,"---- %s(NULL) end ----\n", __func__); - return(NULL); - } - } - - // ************************************************************************ - if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { - if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n"); - stats->clippedMean = NAN; - stats->clippedStdev = NAN; + if (!vectorSampleStdev(inF32, errorsF32, maskU8, maskVal, stats)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev"); + status &= false; } } @@ -2089,4 +1545,40 @@ if (vectorMinMax(inF32, maskU8, maskVal, stats) == 0) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum and maximum"); + status &= false; + } + } + + // ************************************************************************ + if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) { + if (!vectorRobustStats(inF32, errorsF32, maskU8, maskVal, stats)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics")); + status &= false; + } + } + + // ************************************************************************ + if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { + if (!vectorFittedStats(inF32, errorsF32, maskU8, maskVal, stats)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); + status &= false; + } + } + + // ************************************************************************ + if (stats->options & (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2)) { + if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { + psAbort ("stats", "you may not specify both FITTED_MEAN and FITTED_MEAN_V2"); + } + if (!vectorFittedStats_v2(inF32, errorsF32, maskU8, maskVal, stats)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); + status &= false; + } + } + + // ************************************************************************ + if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { + if (!vectorClippedStats(inF32, errorsF32, maskU8, maskVal, stats)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n"); + status &= false; } } @@ -2096,5 +1588,5 @@ psFree(maskU8); psTrace("psLib.math", 3,"---- %s() end ----\n", __func__); - return (stats); + return (status); } @@ -2106,8 +1598,8 @@ } - READ_STAT("MEAN", PS_STAT_SAMPLE_MEAN); - READ_STAT("STDEV", PS_STAT_SAMPLE_STDEV); - READ_STAT("MEDIAN", PS_STAT_SAMPLE_MEDIAN); - READ_STAT("QUARTILE", PS_STAT_SAMPLE_QUARTILE); + READ_STAT("MEAN", PS_STAT_SAMPLE_MEAN); + READ_STAT("STDEV", PS_STAT_SAMPLE_STDEV); + READ_STAT("MEDIAN", PS_STAT_SAMPLE_MEDIAN); + READ_STAT("QUARTILE", PS_STAT_SAMPLE_QUARTILE); READ_STAT("SAMPLE_MEAN", PS_STAT_SAMPLE_MEAN); READ_STAT("SAMPLE_STDEV", PS_STAT_SAMPLE_STDEV); @@ -2118,10 +1610,13 @@ READ_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); READ_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE); - READ_STAT("FITTED", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_STDEV", PS_STAT_ROBUST_STDEV); - READ_STAT("CLIPPED", PS_STAT_CLIPPED_MEAN); - READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); - READ_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); + READ_STAT("FITTED", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); + READ_STAT("FITTED_V2", PS_STAT_FITTED_MEAN_V2); + READ_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN_V2); + READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2); + READ_STAT("CLIPPED", PS_STAT_CLIPPED_MEAN); + READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); + READ_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to parse statistic: %s\n", string); @@ -2146,8 +1641,10 @@ WRITE_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); WRITE_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE); - WRITE_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); - WRITE_STAT("FITTED_STDEV", PS_STAT_ROBUST_STDEV); - WRITE_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); - WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); + WRITE_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); + WRITE_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); + WRITE_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN_V2); + WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2); + WRITE_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); + WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); return string; @@ -2197,4 +1694,6 @@ case PS_STAT_FITTED_MEAN: case PS_STAT_FITTED_STDEV: + case PS_STAT_FITTED_MEAN_V2: + case PS_STAT_FITTED_STDEV_V2: case PS_STAT_CLIPPED_MEAN: case PS_STAT_CLIPPED_STDEV: @@ -2227,4 +1726,8 @@ case PS_STAT_FITTED_STDEV: return stats->fittedStdev; + case PS_STAT_FITTED_MEAN_V2: + return stats->fittedMean; + case PS_STAT_FITTED_STDEV_V2: + return stats->fittedStdev; case PS_STAT_CLIPPED_MEAN: return stats->clippedMean; @@ -2246,2 +1749,190 @@ return NAN; } + +// other private functions used above + +static psF32 QuadraticInverse(psF32 a, + psF32 b, + psF32 c, + psF32 y, + psF32 xLo, + psF32 xHi + ) +{ + 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; + + if (xLo <= x1 && x1 <= xHi) { + return x1; + } + if (xLo <= x2 && x2 <= xHi) { + return x2; + } + return 0.5 * (xLo + xHi); +} + + +/****************************************************************************** +fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general +routine which fits a quadratic to three points and returns the x-value +corresponding to the input y-value. This routine takes psVectors of x/y pairs +as input, and fits a quadratic to the 3 points surrounding element binNum in +the vectors (the midpoint between element i and i+1 is used for x[i]). It +then determines for what value x does that quadratic f(x) = yVal (the input +parameter). + +*****************************************************************************/ +static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, + psVector *yVec, + psS32 binNum, + psF32 yVal + ) +{ + psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); + psTrace("psLib.math", 5, "binNum, yVal is (%d, %f)\n", binNum, yVal); + if (psTraceGetLevel("psLib.math") >= 8) { + PS_VECTOR_PRINT_F32(xVec); + PS_VECTOR_PRINT_F32(yVec); + } + + PS_ASSERT_VECTOR_NON_NULL(xVec, NAN); + PS_ASSERT_VECTOR_NON_NULL(yVec, NAN); + PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN); + PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN); + // PS_ASSERT_VECTORS_SIZE_EQUAL(xVec, yVec, NAN); + PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN); + PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN); + + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psF32 tmpFloat = 0.0f; + + if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) { + // The general case. We have all three points. + x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); + x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + y->data.F64[0] = yVec->data.F32[binNum - 1]; + y->data.F64[1] = yVec->data.F32[binNum]; + y->data.F64[2] = yVec->data.F32[binNum + 1]; + psTrace("psLib.math", 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]); + psTrace("psLib.math", 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); + psTrace("psLib.math", 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); + + // + // Ensure that the y value lies within range of the y values. + // + if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || + ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified yVal, %g, is not within y-range, %g to %g."), + (psF64)yVal, y->data.F64[0], y->data.F64[2]); + } + + // + // Ensure that the y values are monotonic. + // + if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || + ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { + psError(PS_ERR_UNKNOWN, true, + "This routine must be called with monotonically increasing or decreasing data points.\n"); + psFree(x); + psFree(y); + psTrace("psLib.math", 5, "---- %s() end ----\n", __func__); + return NAN; + } + + // Determine the coefficients of the polynomial. + psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); + myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x); + if (myPoly == NULL) { + psError(PS_ERR_UNEXPECTED_NULL, false, + _("Failed to fit a 1-dimensional polynomial to the three specified data points. Returning NAN.")); + psFree(x); + psFree(y); + psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__); + return NAN; + } + psTrace("psLib.math", 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]); + psTrace("psLib.math", 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]); + psTrace("psLib.math", 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]); + psTrace("psLib.math", 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])); + + psTrace("psLib.math", 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); + tmpFloat = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, + x->data.F64[0], x->data.F64[2]); + psFree(myPoly); + + if (isnan(tmpFloat)) { + psError(PS_ERR_UNEXPECTED_NULL, + false, _("Failed to determine the median of the fitted polynomial. Returning NAN.")); + psFree(x); + psFree(y); + psTrace("psLib.math", 5, "---- %s(NAN) end ----\n", __func__); + return(NAN); + } + } else { + // These are special cases where the bin is at the beginning or end of the vector. + if (binNum == 0) { + // We have two points only at the beginning of the vectors x and y. + tmpFloat = 0.5 * (xVec->data.F32[binNum] + + xVec->data.F32[binNum + 1]); + } else if (binNum == (xVec->n - 1)) { + // The special case where we have two points only at the end of + // the vectors x and y. + // XXX: Is this right? + tmpFloat = xVec->data.F32[binNum]; + } else if (binNum == (xVec->n - 2)) { + // XXX: Is this right? + tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]); + } + } + + psTrace("psLib.math", 6, "FIT: return %f\n", tmpFloat); + psFree(x); + psFree(y); + + psTrace("psLib.math", 5, "---- %s(%f) end ----\n", __func__, tmpFloat); + return tmpFloat; +} + +/****************************************************************************** +NOTE: We assume unnormalized gaussians. +*****************************************************************************/ +static psF32 minimizeLMChi2Gauss1D(psVector *deriv, + const psVector *params, + const psVector *coords + ) +{ + psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN); + PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN); + PS_ASSERT_VECTOR_NON_NULL(coords, NAN); + PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN); + PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN); + + psF32 x = coords->data.F32[0]; + psF32 mean = params->data.F32[0]; + psF32 var = params->data.F32[1]; + psF32 dx = (x - mean); + + psF32 gauss = exp (-0.5*PS_SQR(dx)/var); + if (deriv) { + PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN); + PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN); + psF32 tmp = dx * gauss; + deriv->data.F32[0] = tmp / var; + deriv->data.F32[1] = tmp * dx / (var * var); + } + + + psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); + return gauss; +}