Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 2268) +++ trunk/psLib/src/math/psStats.c (revision 2269) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.83 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-02 19:13:03 $ + * @version $Revision: 1.84 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-03 03:30:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -403,16 +403,18 @@ *****************************************************************************/ void p_psVectorSampleMedian(const psVector* restrict myVector, - const psVector* restrict maskVector, psU32 maskVal, psStats* stats) + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psVector* unsortedVector = NULL; // Temporary vector psVector* sortedVector = NULL; // Temporary vector - psS32 i = 0; // Loop index variable - psS32 count = 0; // # of points in this mean? - psS32 nValues = 0; // # of points in vector - float rangeMin = 0.0; // Exclude data below this - float rangeMax = 0.0; // Exclude date above this + psS32 i = 0; // Loop index variable + psS32 count = 0; // # of points in this mean? + psS32 nValues = 0; // # of points in vector + float rangeMin = 0.0; // Exclude data below this + float rangeMax = 0.0; // Exclude date above this // Determine how many data points fit inside this min/max range - // and are not masked, IF the maskVector is not NULL> + // and are not masked, if the maskVector is not NULL> nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); @@ -551,10 +553,12 @@ *****************************************************************************/ void p_psVectorSampleQuartiles(const psVector* restrict myVector, - const psVector* restrict maskVector, psU32 maskVal, psStats* stats) + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psVector* unsortedVector = NULL; // Temporary vector psVector* sortedVector = NULL; // Temporary vector psS32 i = 0; // Loop index variable - psS32 count = 0; // # of points in this mean? + psS32 count = 0; // # of points in this mean. psS32 nValues = 0; // # data points float rangeMin = 0.0; // Exclude data below this @@ -631,5 +635,7 @@ *****************************************************************************/ void p_psVectorSampleStdev(const psVector* restrict myVector, - const psVector* restrict maskVector, psU32 maskVal, psStats* stats) + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable @@ -713,8 +719,10 @@ stats Returns - NULL - *****************************************************************************/ -void p_psVectorClippedStats(const psVector* restrict myVector, - const psVector* restrict maskVector, psU32 maskVal, psStats* stats) + 0 for success. + *****************************************************************************/ +int p_psVectorClippedStats(const psVector* restrict myVector, + const psVector* restrict maskVector, + psU32 maskVal, + psStats* stats) { psS32 i = 0; // Loop index variable @@ -726,11 +734,13 @@ psVector* tmpMask = NULL; // Temporary vector - // Endure that stats->clipIter is within the proper range. + // Ensure that stats->clipIter is within the proper range. if (!((PS_CLIPPED_NUM_ITER_LB <= stats->clipIter) && (stats->clipIter <= PS_CLIPPED_NUM_ITER_UB))) { psError(__func__, "Unallowed value for clipIter (%d).\n", stats->clipIter); - } - // Endure that stats->clipSigma is within the proper range. + return(-1); + } + // Ensure that stats->clipSigma is within the proper range. if (!((PS_CLIPPED_SIGMA_LB <= stats->clipSigma) && (stats->clipSigma <= PS_CLIPPED_SIGMA_UB))) { psError(__func__, "Unallowed value for clipSigma (%f).\n", stats->clipSigma); + return(-1); } // We allocate a temporary mask vector since during the iterative @@ -798,7 +808,12 @@ psFree(tmpMask); -} - -void p_psNormalizeVectorF32_0(psVector* myData) + return(0); +} + +/***************************************************************************** + *****************************************************************************/ +void p_psNormalizeVectorRangeF32(psVector* myData, + float outLow, + float outHigh) { float min = (float)HUGE; @@ -815,22 +830,11 @@ } - // float range = max - min; - // for (i = 0; i < myData->n; i++) { - // myData->data.F32[i] = (myData->data.F32[i] - min) / range; - // } for (i = 0; i < myData->n; i++) { - myData->data.F32[i] = (myData->data.F32[i] - min) / (max - min); - } -} - - - -/***************************************************************************** -p_psNormalizeVectorF32(myData): this is a private function which normalizes the -elements of a vector to a range between 0.0 and 1.0. - -XXX: how about an arbitrary p_psNormalizeVectorF32(myData, min, max)? - *****************************************************************************/ -void p_psNormalizeVectorF32(psVector* myData) + myData->data.F32[i] = outLow + (myData->data.F32[i] - min) * (outHigh - outLow) / (max - min); + } +} +void p_psNormalizeVectorRangeF64(psVector* myData, + float outLow, + float outHigh) { float min = (float)HUGE; @@ -840,63 +844,28 @@ for (i = 0; i < myData->n; i++) { if (myData->data.F32[i] < min) { - min = myData->data.F32[i]; + min = myData->data.F64[i]; } if (myData->data.F32[i] > max) { - max = myData->data.F32[i]; - } - } - - // float range = max - min; - // for (i = 0; i < myData->n; i++) { - // myData->data.F32[i] = (myData->data.F32[i] - min) / range; - // } + max = myData->data.F64[i]; + } + } + for (i = 0; i < myData->n; i++) { - myData->data.F32[i] = (myData->data.F32[i] - 0.5 * (min + max)) / - (0.5 * (max - min)); + myData->data.F64[i] = outLow + (myData->data.F64[i] - min) * (outHigh - outLow) / (max - min); } } /***************************************************************************** -p_psNormalizeVectorF64(myData): this is a private function which normalizes the -elements of a vector to a range between -1.0 and 1.0. -XXX: 0-1 or -1:1? - -XXX: how about an arbitrary p_psNormalizeVectorF64(myData, min, max)? - *****************************************************************************/ -void p_psNormalizeVectorF64(psVector* myData) -{ - float min = (double)HUGE; - float max = (double)-HUGE; - double range = 0.0; - psS32 i = 0; - - for (i = 0; i < myData->n; i++) { - if (myData->data.F64[i] < min) { - min = myData->data.F64[i]; - } - if (myData->data.F64[i] > max) { - max = myData->data.F64[i]; - } - } - - range = max - min; - for (i = 0; i < myData->n; i++) { - myData->data.F64[i] = -1.0 + 2.0 * - ((myData->data.F64[i] - min) / range); - } - - // for (i = 0; i < myData->n; i++) { - // myData->data.F64[i] = (myData->data.F64[i] - 0.5 * (min + max)) / - // (0.5 * (max - min)); - // } -} - -// XXX: how about an arbitrary p_psNormalizeVector(myData, min, max)? -void p_psNormalizeVector(psVector* myData) +psNormalizeVectorRange(myData, low, high): this is a private function which +normalizes the elements of a vector to a range between low and high. + *****************************************************************************/ +void psNormalizeVectorRange(psVector* myData, + float low, + float high) { if (myData->type.type == PS_TYPE_F32) { - p_psNormalizeVectorF32(myData); + p_psNormalizeVectorRangeF32(myData, low, high); } else if (myData->type.type == PS_TYPE_F64) { - p_psNormalizeVectorF64(myData); + p_psNormalizeVectorRangeF64(myData, low, high); } else { psError(__func__, "Unalowable data type.\n"); @@ -914,4 +883,6 @@ XXX: Terminate when f(x)-getThisValue is within some error tolerance. + +XXX: Solve for X analytically. *****************************************************************************/ float p_ps1DPolyMedian(psPolynomial1D* myPoly, @@ -924,6 +895,4 @@ float oldMidpoint = 1.0; float f = 0.0; - - // printf("p_ps1DPolyMedian(%f, %f, %f) \n", rangeLow, rangeHigh, getThisValue); while (numIterations < PS_POLY_MEDIAN_MAX_ITERATIONS) { @@ -956,4 +925,6 @@ and i+1 is used for x[i]). It then determines for what value x does that quadratic f(x) = yVal (the input parameter). + +XXX: After you fit the polynomial, solve for X analytically. *****************************************************************************/ float fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1034,5 +1005,5 @@ stats Returns - NULL + 0 on success. *****************************************************************************/ int p_psVectorRobustStats(const psVector* restrict myVector, @@ -1119,8 +1090,4 @@ tmpStats->clippedStdev / 4.0f); - // The following was necessary to fit a gaussian to the data, since - // gaussian functions produce data between 0.0 and 1.0. - // p_psNormalizeVectorF32(robustHistogramVector); - /************************************************************************** Determine the median/lower/upper quartile bin numbers. @@ -1207,5 +1174,5 @@ psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - p_psNormalizeVectorF32_0(robustHistogramVector); + p_psNormalizeVectorRange(robustHistogramVector, 0.0, 1.0); for (i=0;in;i++) { myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32); @@ -1620,6 +1587,4 @@ } // ************************************************************************ - // XXX: The Stdev calculation requires the mean. Should we assume the - // mean has already been calculated? Or should we always calculate it? if (stats->options & PS_STAT_SAMPLE_STDEV) { p_psVectorSampleMean(inF32, mask, maskVal, stats); @@ -1643,5 +1608,7 @@ if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { - p_psVectorClippedStats(inF32, mask, maskVal, stats); + if (0 != p_psVectorClippedStats(inF32, mask, maskVal, stats)) { + psError(__func__, "p_psVectorClippedStats() failed.\n"); + } } // ************************************************************************