Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 2762) +++ trunk/psLib/src/math/psStats.c (revision 2778) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.104 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-20 21:39:42 $ + * @version $Revision: 1.105 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-21 20:42:07 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -146,4 +146,5 @@ *****************************************************************************/ int p_psVectorSampleMean(const psVector* restrict myVector, + const psVector* restrict errors, const psVector* restrict maskVector, psU32 maskVal, @@ -156,53 +157,112 @@ // If PS_STAT_USE_RANGE is requested, then we enter a slightly different // loop. - if (stats->options & PS_STAT_USE_RANGE) { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - // Check if the data is with the specified range - if (!(maskVal & maskVector->data.U8[i]) && - (stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { + if (errors == NULL) { + if (stats->options & PS_STAT_USE_RANGE) { + if (maskVector != NULL) { + for (i = 0; i < myVector->n; i++) { + // Check if the data is with the specified range + if (!(maskVal & maskVector->data.U8[i]) && + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + mean += myVector->data.F32[i]; + count++; + } + } + if (count != 0) { + mean /= (float)count; + } else { + mean = NAN; + } + + } else { + for (i = 0; i < myVector->n; i++) { + if ((stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + mean += myVector->data.F32[i]; + count++; + } + } + if (count != 0) { + mean /= (float)count; + } else { + mean = NAN; + } + } + } else { + if (maskVector != NULL) { + for (i = 0; i < myVector->n; i++) { + if (!(maskVal & maskVector->data.U8[i])) { + mean += myVector->data.F32[i]; + count++; + } + } + if (count != 0) { + mean /= (float)count; + } else { + mean = NAN; + } + } else { + for (i = 0; i < myVector->n; i++) { mean += myVector->data.F32[i]; - count++; - } - } - if (count != 0) { - mean /= (float)count; + } + mean /= (float)myVector->n; + } + } + } else { + psF32 errorDivisor = 0.0; + if (stats->options & PS_STAT_USE_RANGE) { + if (maskVector != NULL) { + for (i = 0; i < myVector->n; i++) { + // Check if the data is with the specified range + if (!(maskVal & maskVector->data.U8[i]) && + (stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + mean += myVector->data.F32[i]/errors->data.F32[i]; + errorDivisor+= (1.0 / errors->data.F32[i]); + count++; + } + } + if (count != 0) { + mean /= errorDivisor; + } else { + mean = NAN; + } + } else { - mean = NAN; - } - - } else { - for (i = 0; i < myVector->n; i++) { - if ((stats->min <= myVector->data.F32[i]) && - (myVector->data.F32[i] <= stats->max)) { - mean += myVector->data.F32[i]; - count++; - } - } - if (count != 0) { - mean /= (float)count; + for (i = 0; i < myVector->n; i++) { + if ((stats->min <= myVector->data.F32[i]) && + (myVector->data.F32[i] <= stats->max)) { + mean += myVector->data.F32[i]/errors->data.F32[i]; + errorDivisor+= (1.0 / errors->data.F32[i]); + count++; + } + } + if (count != 0) { + mean /= errorDivisor; + } else { + mean = NAN; + } + } + } else { + if (maskVector != NULL) { + for (i = 0; i < myVector->n; i++) { + if (!(maskVal & maskVector->data.U8[i])) { + mean += myVector->data.F32[i]/errors->data.F32[i]; + errorDivisor+= (1.0 / errors->data.F32[i]); + count++; + } + } + if (count != 0) { + mean /= errorDivisor; + } else { + mean = NAN; + } } else { - mean = NAN; - } - } - } else { - if (maskVector != NULL) { - for (i = 0; i < myVector->n; i++) { - if (!(maskVal & maskVector->data.U8[i])) { - mean += myVector->data.F32[i]; - count++; - } - } - if (count != 0) { - mean /= (float)count; - } else { - mean = NAN; - } - } else { - for (i = 0; i < myVector->n; i++) { - mean += myVector->data.F32[i]; - } - mean /= (float)myVector->n; + for (i = 0; i < myVector->n; i++) { + mean += myVector->data.F32[i]/errors->data.F32[i]; + errorDivisor+= (1.0 / errors->data.F32[i]); + } + mean /= errorDivisor; + } } } @@ -349,6 +409,6 @@ /****************************************************************************** p_psVectorNValues(myVector, maskVector, maskVal, stats): This routine returns -"true" if the inputPsVector as 1 or more valid elements (a valid element is an -unmasked element within the specifined min/max range). Otherwise, return +"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". *****************************************************************************/ @@ -447,4 +507,6 @@ median of the input vector. Returns true on success (including if there were no valid input vector elements). + +XXX: Use static vectors for sort arrays. *****************************************************************************/ bool p_psVectorSampleMedian(const psVector* restrict myVector, @@ -522,5 +584,6 @@ } - // Calculate the median exactly. + // Calculate the median exactly. Use the average if the number of samples + // is even. if (0 == (nValues % 2)) { stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues / 2) - 1] + @@ -723,8 +786,9 @@ *****************************************************************************/ -void p_psVectorSampleStdev(const psVector* restrict myVector, - const psVector* restrict maskVector, - psU32 maskVal, - psStats* stats) +void p_psVectorSampleStdevOLD(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable @@ -739,5 +803,5 @@ // calculated, then call p_psVectorSampleMean() if (0 != isnan(stats->sampleMean)) { - p_psVectorSampleMean(myVector, maskVector, maskVal, stats); + p_psVectorSampleMean(myVector, errors, maskVector, maskVal, stats); } // If the mean is NAN, then generate a warning and set the stdev to NAN. @@ -801,19 +865,121 @@ } else { countFloat = (float)countInt; - #ifdef DARWIN - - stats->sampleStdev = (float)sqrt((sumSquares - - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - #else - - stats->sampleStdev = sqrtf((sumSquares - - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - #endif - - } -} - + stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); + } +} /****************************************************************************** -p_psVectorClippedStats(myVector, maskVector, maskVal, stats): calculates the +p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the +stdev of the input vector. +Inputs + myVector + maskVector + maskVal + stats +Returns + NULL + + *****************************************************************************/ +void p_psVectorSampleStdev(const psVector* restrict myVector, + const psVector* restrict errors, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) +{ + psS32 i = 0; // Loop index variable + psS32 countInt = 0; // # of data points being used + float countFloat = 0.0; // # of data points being used + float mean = 0.0; // The mean + float diff = 0.0; // Used in calculating stdev + float sumSquares = 0.0; // temporary variable + float sumDiffs = 0.0; // temporary variable + // psF32 sum1; + // psF32 sum2; + psF32 errorDivisor; + + // 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"); + 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++; + errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); + } + } + } 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++; + errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); + } + } + 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++; + errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); + } + } + } else { + for (i = 0; i < myVector->n; i++) { + diff = myVector->data.F32[i] - mean; + sumSquares += (diff * diff); + sumDiffs += diff; + countInt++; + } + countInt = myVector->n; + errorDivisor+= (1.0 / PS_SQR(errors->data.F32[i])); + } + } + + 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 { + // XXX: The ADD specifies this as the definition of the standard + // deviation if the errors are known. Verify this with IfA: none of + // the data points in the vector are used. Verify that the masks and + // data ranges are used correctly. + if (errors != NULL) { + stats->sampleStdev = (1.0 / PS_SQRT_F32(errorDivisor)); + } else { + countFloat = (float)countInt; + stats->sampleStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); + + } + } +} + +/****************************************************************************** +p_psVectorClippedStats(myVector, errors, maskVector, maskVal, stats): calculates the clipped stats (mean or stdev) of the input vector. @@ -829,4 +995,5 @@ *****************************************************************************/ int p_psVectorClippedStats(const psVector* restrict myVector, + const psVector* restrict errors, const psVector* restrict maskVector, psU32 maskVal, @@ -873,5 +1040,5 @@ // 2. Compute the sample standard deviation. - p_psVectorSampleStdev(myVector, maskVector, maskVal, stats); + p_psVectorSampleStdev(myVector, errors, maskVector, maskVal, stats); if (isnan(stats->sampleStdev)) { psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleStdev returned NAN\n"); @@ -894,14 +1061,22 @@ // 5. Repeat N times: for (i = 0; i < stats->clipIter; i++) { - for (j = 0; j < myVector->n; j++) { - // a) Exclude all values x_i for which |x_i - x| > K * stdev - if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { - tmpMask->data.U8[i] = 0xff; + // a) Exclude all values x_i for which |x_i - x| > K * stdev + if (errors != NULL) { + for (j = 0; j < myVector->n; j++) { + if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * errors->data.F32[j])) { + tmpMask->data.U8[i] = 0xff; + } + } + } else { + for (j = 0; j < myVector->n; j++) { + if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { + tmpMask->data.U8[i] = 0xff; + } } } // b) compute new mean and stdev - p_psVectorSampleMean(myVector, tmpMask, maskVal, stats); - p_psVectorSampleStdev(myVector, tmpMask, maskVal, stats); + p_psVectorSampleMean(myVector, errors, tmpMask, maskVal, stats); + p_psVectorSampleStdev(myVector, errors, tmpMask, maskVal, stats); // If the new mean and stdev are not NAN, then we use them. @@ -1242,4 +1417,5 @@ *****************************************************************************/ int p_psVectorRobustStats(const psVector* restrict myVector, + const psVector* restrict errors, const psVector* restrict maskVector, psU32 maskVal, @@ -1274,5 +1450,6 @@ // by computing the clipped standard deviation of the vector, and dividing // that by 10.0; - int rc = p_psVectorClippedStats(myVector, maskVector, maskVal, tmpStats); + //XXX: add errors + int rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); if (rc != 0) { psError(PS_ERR_UNEXPECTED_NULL, @@ -1333,5 +1510,5 @@ // Populate the histogram array. - psVectorHistogram(robustHistogram, myVector, maskVector, maskVal); + psVectorHistogram(robustHistogram, myVector, errors, maskVector, maskVal); // Smooth the histogram, Gaussian-style. @@ -1416,5 +1593,5 @@ } } - myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); + myStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); // Using the above (myMean, myStdev) as initial estimates, we fit a @@ -1674,8 +1851,8 @@ /***************************************************************************** -psVectorHistogram(out, in, 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. +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: @@ -1686,7 +1863,10 @@ Returns: The histogram structure "out". + +XXX: Waiting for direction on how to use the errors parameter. *****************************************************************************/ psHistogram* psVectorHistogram(psHistogram* out, const psVector* restrict in, + const psVector* restrict errors, const psVector* restrict mask, psU32 maskVal) @@ -1700,4 +1880,8 @@ PS_VECTOR_CHECK_SIZE_EQUAL(in, mask, NULL); PS_VECTOR_CHECK_TYPE(mask, PS_TYPE_U8, NULL); + } + if (errors != NULL) { + PS_VECTOR_CHECK_SIZE_EQUAL(in, errors, NULL); + PS_VECTOR_CHECK_TYPE(errors, in->type.type, NULL); } @@ -1828,4 +2012,5 @@ psStats* psVectorStats(psStats* stats, psVector* in, + psVector* errors, psVector* mask, psU32 maskVal) @@ -1837,11 +2022,23 @@ PS_VECTOR_CHECK_TYPE(mask, PS_TYPE_U8, stats); } + if (errors != NULL) { + PS_VECTOR_CHECK_SIZE_EQUAL(errors, in, stats); + PS_VECTOR_CHECK_TYPE(errors, in->type.type, stats); + } + psVector* inF32; - psS32 mustFreeTmp = 1; + psVector* errorsF32; + psS32 mustFreeVectorIn = 1; + psS32 mustFreeVectorErrors = 1; inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { inF32 = in; - mustFreeTmp = 0; + mustFreeVectorIn = 0; + } + errorsF32 = p_psConvertToF32((psVector *) errors); + if (errorsF32 == NULL) { + errorsF32 = errors; + mustFreeVectorErrors = 0; } @@ -1852,5 +2049,5 @@ // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_MEAN) { - if (0 != p_psVectorSampleMean(inF32, mask, maskVal, stats)) { + if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); @@ -1866,9 +2063,9 @@ // ************************************************************************ if (stats->options & PS_STAT_SAMPLE_STDEV) { - if (0 != p_psVectorSampleMean(inF32, mask, maskVal, stats)) { + if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n"); } - p_psVectorSampleStdev(inF32, mask, maskVal, stats); + p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats); } // ************************************************************************ @@ -1887,5 +2084,5 @@ (stats->options & PS_STAT_ROBUST_STDEV) || (stats->options & PS_STAT_ROBUST_QUARTILE)) { - if (0 != p_psVectorRobustStats(inF32, mask, maskVal, stats)) { + if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED); @@ -1894,11 +2091,12 @@ } + // XXX: Different conditions for return -1 and -2? if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { - if (-1 == p_psVectorClippedStats(inF32, mask, maskVal, stats)) { + if (-1 > p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, - "Failed to calculate clipped statistics for input psVector."); + "Failed to calculate clipped statistics for input psVector.\n"); stats->clippedMean = NAN; stats->clippedStdev = NAN; - } else if (-2 == p_psVectorClippedStats(inF32, mask, maskVal, stats)) { + } else if (-2 == p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats)) { psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector."); stats->clippedMean = NAN; @@ -1925,7 +2123,10 @@ } - if (mustFreeTmp == 1) { + if (mustFreeVectorIn == 1) { psFree(inF32); } + if (mustFreeVectorErrors == 1) { + psFree(errorsF32); + } return (stats); }