Changeset 6215

Timestamp:

Jan 26, 2006, 1:49:11 PM (20 years ago)

Author:

gusciora

Message:

Made the requested changes to the psStats structure and added the fitted
statistics.

Location:

trunk/psLib

Files:

• src/math/psStats.c (modified) (26 diffs)
• src/math/psStats.h (modified) (4 diffs)
• test/math/tst_psStats00.c (modified) (17 diffs)
• test/math/tst_psStats03.c (modified) (1 diff)
• test/math/tst_psStats07.c (modified) (10 diffs)

trunk/psLib/src/math/psStats.c

@@ -14 +14 @@
  *      stats->binsize
  *
- *  @version $Revision: 1.162 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-26 21:10:22 $
+ * XXX: Must do
+ * nSubsample points
+ * use ->min and ->max (PS_STAT_USE_RANGE)
+ * use ->binsize (PS_STAT_USE_BINSIZE)
+ *
+ *
+ *
+ *
+ *  @version $Revision: 1.163 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-26 23:49:11 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -97 +105 @@
         return true;
 
-    case PS_STAT_ROBUST_MEAN:
-        *value = stats->robustMean;
-        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-        return true;
-
     case PS_STAT_ROBUST_MEDIAN:
         *value = stats->robustMedian;
-        psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
-        return true;
-
-    case PS_STAT_ROBUST_MODE:
-        *value = stats->robustMode;
         psTrace(__func__, 4, "---- %s(true) end ----\n", __func__);
         return true;
@@ -449 +447 @@
 unmasked element within the specified min/max range).  Otherwise, return
 "false".
+ 
+XXX: Can you use psVectorCountPixelMask here?
  *****************************************************************************/
 bool p_psVectorCheckNonEmpty(const psVector* myVector,
@@ -503 +503 @@
 number of non-masked pixels in the vector that fall within the min/max
 range, if specified.
+ 
+XXX: Can you use psVectorCountPixelMask here?
  *****************************************************************************/
 psS32 p_psVectorNValues(const psVector* myVector,
@@ -766 +768 @@
         PS_VECTOR_PRINT_F32(smooth);
     }
-    psTrace(__func__, 4, "---- %s(psVector) end ----\n", __func__);
+    psTrace(__func__, 4, "---- %s() end ----\n", __func__);
     return(smooth);
 }
@@ -858 +860 @@
 
 /******************************************************************************
-p_psVectorSampleStdev(myVector, maskVector, maskVal, stats): calculates the
+p_psVectorSampleStdevOLD(myVector, maskVector, maskVal, stats): calculates the
 stdev of the input vector.
 Inputs
@@ -868 +870 @@
     NULL
  
+XXX: remove this
  *****************************************************************************/
 void p_psVectorSampleStdevOLD(const psVector* myVector,
@@ -1126 +1129 @@
         p_psMemSetPersistent(statsTmp, true);
     } else {
-        // EAM : initialize structure if already allocated
+        // Initialize structure if already allocated
         statsTmp->sampleMean = NAN;
         statsTmp->sampleMedian = NAN;
@@ -1132 +1135 @@
         statsTmp->sampleUQ = NAN;
         statsTmp->sampleLQ = NAN;
-        statsTmp->robustMean = NAN;
         statsTmp->robustMedian = NAN;
-        statsTmp->robustMode = NAN;
         statsTmp->robustStdev = NAN;
         statsTmp->robustUQ = NAN;
         statsTmp->robustLQ = NAN;
-        statsTmp->robustN50 = -1;            // XXX: This is never used
-        statsTmp->robustNfit = -1;
+        statsTmp->robustN50 = -1;
+        statsTmp->fittedMean = NAN;
+        statsTmp->fittedStdev = NAN;
+        statsTmp->fittedNfit = -1;
         statsTmp->clippedMean = NAN;
         statsTmp->clippedStdev = NAN;
-        statsTmp->clippedNvalues = -1;     // XXX: This is never used
+        statsTmp->clippedNvalues = -1;
         statsTmp->clipSigma = 3.0;
         statsTmp->clipIter = 3;
         statsTmp->min = NAN;
         statsTmp->max = NAN;
-        statsTmp->binsize = NAN;          // XXX: This is never used
+        statsTmp->binsize = NAN;
+        statsTmp->nSubsample = 100000;
     }
 
@@ -1688 +1692 @@
         }
     }
-    psBool iterate = true;
+    psS32 iterate = 1;
     psF32 sigma;
     psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX);
 
-    while (iterate) {
+    while (iterate > 0) {
         psTrace(__func__, 6, "Iterating on Bin size.\n");
-        //
-        // Determine the bin size of the robust histogram.  This is done
-        // by computing the total range of data values and dividing by 1000.0.
-        //
-        rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax);
-        rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax);
-        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
-            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-            psFree(tmpStatsMinMax);
-            psFree(tmpMaskVec);
-            psFree(tmpScalar);
-            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-            return(1);
-        }
-        psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
-        psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f;
-        psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize);
+        psF32 binSize = 0.0;
+        if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) {
+            // Set initial bin size to the specified value.
+            binSize = stats->binsize;
+            psTrace(__func__, 6, "Setting initial robust bin size to %.2f\n", binSize);
+        } else {
+            // Determine the bin size of the robust histogram.  This is done
+            // by computing the total range of data values and dividing by 1000.0.
+
+            rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax);
+            rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax);
+            if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
+                psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+                psFree(tmpStatsMinMax);
+                psFree(tmpMaskVec);
+                psFree(tmpScalar);
+                psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+                return(1);
+            }
+            psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
+            binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f;
+            psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize);
+        }
+        psTrace(__func__, 6, "Initial robust bin size is %.2f\n", binSize);
 
         //
@@ -1726 +1737 @@
             }
             // XXX: Set these to the number of unmasked data points?
-            stats->robustNfit = 0.0;
-            stats->robustN50 = 0.0;
+            stats->robustN50 = 0;
             psFree(tmpStatsMinMax);
             psFree(tmpMaskVec);
@@ -1900 +1910 @@
         sigma = (binHiF32 - binLoF32) / 2.0;
         psTrace(__func__, 6, "The current sigma is %f.\n", sigma);
+        stats->robustStdev = sigma;
 
         //
@@ -1921 +1932 @@
             psFree(robustHistogram);
             psFree(cumulativeRobustHistogram);
-
+            iterate++;
         } else {
             psTrace(__func__, 6, "*************: No more iteration.  sigma is %f\n", sigma);
-            iterate = false;
+            iterate = 0;
         }
     }
@@ -1981 +1992 @@
     psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32);
 
-    //
-    //
-    // New algorithm for calculated mean and stdev.
-    //
-    //
     psS32 N50 = 0;
     for (psS32 i = 0 ; i < myVector->n ; i++) {
+        // XXX: use maskVal here?
         if ((0 == tmpMaskVec->data.U8[i]) &&
                 (binLo25F32 <= myVector->data.F32[i]) &&
@@ -1997 +2004 @@
     psTrace(__func__, 6, "The robustN50 is %d.\n", N50);
 
-    psF32 dN = (psF32) (0.17 * N50);
-    if (dN < 1.0) {
-        dN = 1.0;
-    } else if (dN > 4.0) {
-        dN = 4.0;
-    }
-    psF32 newBinSize = sigma / dN;
-
-    rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
-    rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
-    if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
-        psFree(tmpStatsMinMax);
-        psFree(robustHistogram);
-        psFree(cumulativeRobustHistogram);
-        psFree(tmpScalar);
-        psFree(tmpMaskVec);
-        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-        return(1);
-    }
-
-    numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize);
-    psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
-    psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize);
-    psTrace(__func__, 6, "The numBins is %d\n", numBins);
-
-    psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
-    newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1);
-    if (psTraceGetLevel(__func__) >= 8) {
-        PS_VECTOR_PRINT_F32(newHistogram->nums);
-    }
-
-    //
-    // Smooth the resulting histogram with a Gaussian with sigma_s = 1
-    // bin.
-    //
-    psF32 sigma_s = newBinSize;
-
-    psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s);
-    if (psTraceGetLevel(__func__) >= 8) {
-        PS_VECTOR_PRINT_F32(newHistogramSmoothed);
-    }
-
-    //
-    // Find the bin with the peak value in the range 2 sigma of the
-    // robust histogram median.
-    //
-
-    psS32 binMin = 0;
-    psS32 binMax = 0;
-    tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
-    if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
-        binMin = 0;
-    } else {
-        binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
-    }
-
-    tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
-    if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
-        binMax = newHistogram->bounds->n-1;
-    } else {
-        binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
-    }
-    if ((binMin < 0) || (binMax < 0)) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n");
-        psFree(tmpMaskVec);
-        psFree(robustHistogram);
-        psFree(cumulativeRobustHistogram);
-        psFree(tmpStatsMinMax);
-        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-        return(1);
-    }
-
-    psS32 binNum = binMin;
-    psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum];
-    for (psS32 i = binMin+1 ; i <= binMax ; i++) {
-        if (newHistogramSmoothed->data.F32[i] > binMaxNums) {
-            binNum = i;
-            binMaxNums = newHistogramSmoothed->data.F32[i];
-        }
-    }
-    psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums);
-
-    //
-    // Fit a Gaussian to the bins in the range 20 sigma of the robust
-    // histogram median.
-    //
-    tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
-    if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
-        binMin = 0;
-    } else {
-        binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
-        // XXX: check for errors here.
-    }
-    tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
-    if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
-        binMax = newHistogramSmoothed->n - 1;
-    } else {
-        binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
-        // XXX: check for errors here.
-    }
-    psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
-    psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
-    psArray *x = psArrayAlloc((1 + (binMax - binMin)));
-    stats->robustNfit = 0;
-    psS32 j = 0;
-
-    for (psS32 i = binMin ; i <= binMax ; i++) {
-        y->data.F32[j] = newHistogramSmoothed->data.F32[i];
-        x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
-        ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
-        xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i);
-
-        stats->robustNfit+= newHistogramSmoothed->data.F32[i];
-        j++;
-    }
-    if (psTraceGetLevel(__func__) >= 8) {
-        // XXX: Print the x array somehow.
-        PS_VECTOR_PRINT_F32(y);
-    }
-
-    // XXX: Use the min/max routines for this
-    psF32 minY = FLT_MAX;
-    psF32 maxY = -FLT_MAX;
-    for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
-        if (y->data.F32[i] > maxY) {
-            maxY = y->data.F32[i];
-        }
-        if (y->data.F32[i] < minY) {
-            minY = y->data.F32[i];
-        }
-    }
-    //
-    // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
-    // XXX: Use the normalize routines for this.
-    //
-    for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
-        y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
-    }
-
-    //
-    psMinimization *min = psMinimizationAlloc(100, 0.01);
-    psVector *params = psVectorAlloc(2, PS_TYPE_F32);
-    // Initial guess for the mean ([0]) and standard dev.
-    params->data.F32[0] = stats->robustMedian;
-    params->data.F32[1] = sigma;
-    if (psTraceGetLevel(__func__) >= 8) {
-        PS_VECTOR_PRINT_F32(params);
-        PS_VECTOR_PRINT_F32(y);
-    }
-    psFree(xTmp);
-    rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D);
-
-    if (rcBool != true) {
-        psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
-        psFree(tmpStatsMinMax);
-        psFree(robustHistogram);
-        psFree(cumulativeRobustHistogram);
-        psFree(tmpScalar);
+    // ************************************************************************
+    if ((stats->options & PS_STAT_FITTED_MEAN) ||
+            (stats->options & PS_STAT_FITTED_STDEV)) {
+        //
+        // New algorithm for calculated mean and stdev.
+        //
+
+        // XXX: Where is this documented?
+        psF32 dN = (psF32) (0.17 * N50);
+        if (dN < 1.0) {
+            dN = 1.0;
+        } else if (dN > 4.0) {
+            dN = 4.0;
+        }
+        psF32 newBinSize = sigma / dN;
+
+        rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
+        rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax);
+        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(tmpStatsMinMax);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpScalar);
+            psFree(tmpMaskVec);
+            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+            return(1);
+        }
+
+        numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize);
+        psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max);
+        psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize);
+        psTrace(__func__, 6, "The numBins is %d\n", numBins);
+
+        psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins);
+        newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1);
+        if (psTraceGetLevel(__func__) >= 8) {
+            PS_VECTOR_PRINT_F32(newHistogram->nums);
+        }
+
+        //
+        // FITTED STATISTICS HERE
+        //
+        // Smooth the resulting histogram with a Gaussian with sigma_s = 1
+        // bin.
+        //
+        psF32 sigma_s = newBinSize;
+
+        psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s);
+        if (psTraceGetLevel(__func__) >= 8) {
+            PS_VECTOR_PRINT_F32(newHistogramSmoothed);
+        }
+
+        //
+        // Find the bin with the peak value in the range 2 sigma of the
+        // robust histogram median.
+        //
+
+        psS32 binMin = 0;
+        psS32 binMax = 0;
+        tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma);
+        if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) {
+            binMin = 0;
+        } else {
+            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+        }
+
+        tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma);
+        if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) {
+            binMax = newHistogram->bounds->n-1;
+        } else {
+            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+        }
+        if ((binMin < 0) || (binMax < 0)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n");
+            psFree(tmpMaskVec);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpStatsMinMax);
+            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+            return(1);
+        }
+
+        psS32 binNum = binMin;
+        psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum];
+        for (psS32 i = binMin+1 ; i <= binMax ; i++) {
+            if (newHistogramSmoothed->data.F32[i] > binMaxNums) {
+                binNum = i;
+                binMaxNums = newHistogramSmoothed->data.F32[i];
+            }
+        }
+        psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums);
+
+        //
+        // Fit a Gaussian to the bins in the range 20 sigma of the robust
+        // histogram median.
+        //
+        tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma);
+        if (tmpScalar->data.F32 < tmpStatsMinMax->min) {
+            binMin = 0;
+        } else {
+            binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            // XXX: check for errors here.
+        }
+        tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma);
+        if (tmpScalar->data.F32 > tmpStatsMinMax->max) {
+            binMax = newHistogramSmoothed->n - 1;
+        } else {
+            binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar);
+            // XXX: check for errors here.
+        }
+        psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
+        psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32);
+        psArray *x = psArrayAlloc((1 + (binMax - binMin)));
+        psS32 j = 0;
+
+        for (psS32 i = binMin ; i <= binMax ; i++) {
+            y->data.F32[j] = newHistogramSmoothed->data.F32[i];
+            x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32);
+            ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i);
+            xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i);
+            j++;
+        }
+        if (psTraceGetLevel(__func__) >= 8) {
+            // XXX: Print the x array somehow.
+            PS_VECTOR_PRINT_F32(y);
+        }
+
+        // XXX: Use the min/max routines for this
+        psF32 minY = FLT_MAX;
+        psF32 maxY = -FLT_MAX;
+        for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
+            if (y->data.F32[i] > maxY) {
+                maxY = y->data.F32[i];
+            }
+            if (y->data.F32[i] < minY) {
+                minY = y->data.F32[i];
+            }
+        }
+        //
+        // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
+        // XXX: Use the normalize routines for this.
+        //
+        for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) {
+            y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
+        }
+
+        //
+        psMinimization *min = psMinimizationAlloc(100, 0.01);
+        psVector *params = psVectorAlloc(2, PS_TYPE_F32);
+        // Initial guess for the mean ([0]) and standard dev.
+        params->data.F32[0] = stats->robustMedian;
+        params->data.F32[1] = sigma;
+        if (psTraceGetLevel(__func__) >= 8) {
+            PS_VECTOR_PRINT_F32(params);
+            PS_VECTOR_PRINT_F32(y);
+        }
+        psFree(xTmp);
+        rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D);
+
+        if (rcBool != true) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
+            psFree(tmpStatsMinMax);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpScalar);
+            psFree(newHistogram);
+            psFree(x);
+            psFree(y);
+            psFree(min);
+            psFree(params);
+            psFree(tmpMaskVec);
+            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+            return(1);
+        }
+        if (psTraceGetLevel(__func__) >= 8) {
+            PS_VECTOR_PRINT_F32(params);
+        }
+
+        //
+        // The fitted mean mean_r is derived directly from the fitted
+        // Gaussian mean.
+        //
+        stats->fittedMean = params->data.F32[0];
+        psTrace(__func__, 6, "The fitted mean is %f.\n", params->data.F32[0]);
+
+        //
+        // The fitted standard deviation, SIGMA_r is determined by
+        // subtracting the smoothing scale in quadrature:
+        //     SIGMA_r^2 = SIGMA^2 - sigma_s^2
+        //
+        stats->fittedStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s));
+        psTrace(__func__, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
+
+        psFree(newHistogramSmoothed);
         psFree(newHistogram);
         psFree(x);
@@ -2161 +2204 @@
         psFree(min);
         psFree(params);
-        psFree(tmpMaskVec);
-        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
-        return(1);
-    }
-    if (psTraceGetLevel(__func__) >= 8) {
-        PS_VECTOR_PRINT_F32(params);
-    }
-
-    //
-    // The robust mean mean_r is derived directly from the fitted
-    // Gaussian mean.
-    //
-    stats->robustMean = params->data.F32[0];
-    psTrace(__func__, 6, "The robust mean is %f.\n", params->data.F32[0]);
-
-    //
-    // The robust standard deviation, SIGMA_r is determined by
-    // subtracting the smoothing scale in quadrature:
-    //     SIGMA_r^2 = SIGMA^2 - sigma_s^2
-    //
-    stats->robustStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s));
-    psTrace(__func__, 6, "The robust stdev is %f.\n", stats->robustStdev);
-
-    psFree(newHistogramSmoothed);
+    }
+
     psFree(tmpStatsMinMax);
     psFree(cumulativeRobustHistogram);
     psFree(tmpScalar);
-    psFree(newHistogram);
-    psFree(x);
-    psFree(y);
-    psFree(min);
-    psFree(params);
     psFree(tmpMaskVec);
     psFree(robustHistogram);
@@ -2236 +2252 @@
     newStruct->sampleUQ = NAN;
     newStruct->sampleLQ = NAN;
-    newStruct->robustMean = NAN;
     newStruct->robustMedian = NAN;
-    newStruct->robustMode = NAN;
     newStruct->robustStdev = NAN;
     newStruct->robustUQ = NAN;
     newStruct->robustLQ = NAN;
     newStruct->robustN50 = -1;            // XXX: This is never used
-    newStruct->robustNfit = -1;
+    newStruct->fittedMean = NAN;
+    newStruct->fittedStdev = NAN;
+    newStruct->fittedNfit = -1;
     newStruct->clippedMean = NAN;
     newStruct->clippedStdev = NAN;
     newStruct->clippedNvalues = -1;     // XXX: This is never used
-    // XXX: Where do these values come from?
     newStruct->clipSigma = 3.0;
     newStruct->clipIter = 3;
@@ -2253 +2268 @@
     newStruct->max = NAN;
     newStruct->binsize = NAN;          // XXX: This is never used
+    newStruct->nSubsample = 100000;
     newStruct->options = options;
 
-    psTrace(__func__, 3, "---- %s(psStats) end  ----\n", __func__);
+    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
     return (newStruct);
 }
@@ -2318 +2334 @@
     newHist->uniform = true;
 
-    psTrace(__func__, 3, "---- %s(psHistogram) end  ----\n", __func__);
+    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
     return (newHist);
 }
@@ -2363 +2379 @@
     newHist->uniform = false;
 
-    psTrace(__func__, 3, "---- %s(psHistogram) end  ----\n", __func__);
+    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
     return (newHist);
 }
@@ -2594 +2610 @@
     }
 
-    psTrace(__func__, 3, "---- %s(psHistogram) end  ----\n", __func__);
+    psTrace(__func__, 3, "---- %s() end  ----\n", __func__);
     return (out);
 }
@@ -2674 +2690 @@
     }
 
-    psTrace(__func__, 4,"---- %s(psVector) end  ----\n", __func__);
+    psTrace(__func__, 4,"---- %s() end  ----\n", __func__);
     return (tmp);
 }
@@ -2708 +2724 @@
         PS_ASSERT_VECTOR_TYPE(errors, in->type.type, stats);
     }
+    // XXX: Assert that "in" is F64, F32, U16, or S8
 
     psVector* inF32 = NULL;
@@ -2729 +2746 @@
     }
 
+    if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) {
+        PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, stats);
+    }
+
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_MEAN) {
         if (0 != p_psVectorSampleMean(inF32, errorsF32, mask, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: psVectorStats(): p_psVectorSampleMean() returned an error.\n");
-        }
-    }
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n");
+            stats->sampleMean = NAN;
+        }
+    }
+
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_MEDIAN) {
         if (false == p_psVectorSampleMedian(inF32, mask, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: psVectorStats(): p_psVectorSampleMedian() returned an error.\n");
-        }
-    }
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian() returned an error.\n");
+            stats->sampleMedian = NAN;
+        }
+    }
+
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_STDEV) {
         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, errorsF32, mask, maskVal, stats);
-    }
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMean() returned an error.\n");
+            stats->sampleMean = NAN;
+        } else {
+            p_psVectorSampleStdev(inF32, errorsF32, mask, maskVal, stats);
+        }
+    }
+
     // ************************************************************************
     if (stats->options & PS_STAT_SAMPLE_QUARTILE) {
         if (false == p_psVectorSampleQuartiles(inF32, mask, maskVal, stats)) {
-            psLogMsg(__func__, PS_LOG_WARN,
-                     "WARNING: psVectorStats(): p_psVectorSampleQuartiles() returned an error.\n");
-        }
-    }
-    // Since the various robust stats quantities share much computation, they
-    // are grouped together in a single private function:
-    // p_psVectorRobustStats()
-    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
-            (stats->options & PS_STAT_ROBUST_MEDIAN) ||
-            (stats->options & PS_STAT_ROBUST_MODE) ||
+            psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleQuartiles() returned an error.\n");
+            stats->sampleLQ = NAN;
+            stats->sampleUQ = NAN;
+        }
+    }
+
+    // ************************************************************************
+    if ((stats->options & PS_STAT_ROBUST_MEDIAN) ||
             (stats->options & PS_STAT_ROBUST_STDEV) ||
-            (stats->options & PS_STAT_ROBUST_QUARTILE)) {
+            (stats->options & PS_STAT_ROBUST_QUARTILE) ||
+            (stats->options & PS_STAT_FITTED_MEAN) ||
+            (stats->options & PS_STAT_FITTED_STDEV)) {
         if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false,
-                    PS_ERRORTEXT_psStats_STATS_FAILED);
-            // XXX: Set to NAN
-            // XXX: Is this the right thing to do?
-            // XXX: If so, do it for other funcs?
+            psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED);
             psFree(stats);
             psTrace(__func__, 3,"---- %s(NULL) end  ----\n", __func__);
@@ -2778 +2799 @@
     }
 
-    // XXX: Different conditions for return -1 and -2?
+    // ************************************************************************
     if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) {
         psS32 rc = p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats);
-        if (-1 == rc) {
-            psError(PS_ERR_UNKNOWN, false,
-                    "Failed to calculate clipped statistics for input psVector.\n");
+        if (rc < 0) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics for input psVector.\n");
             stats->clippedMean = NAN;
             stats->clippedStdev = NAN;
-        } else if (-2 == rc) {
-            psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector.");
-            stats->clippedMean = NAN;
-            stats->clippedStdev = NAN;
-        }
-    }
+        }
+    }
+
     // ************************************************************************
     if (stats->options & PS_STAT_MAX) {
         if (0 != p_psVectorMax(inF32, mask, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false,
-                    "Failed to calculate vector maximum");
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector maximum");
             stats->max = NAN;
         }
     }
+
     // ************************************************************************
     if (stats->options & PS_STAT_MIN) {
         if (0 != p_psVectorMin(inF32, mask, maskVal, stats)) {
-            psError(PS_ERR_UNKNOWN, false,
-                    "Failed to calculate vector minimum");
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector minimum");
             stats->min = NAN;
         }
@@ -2815 +2831 @@
         psFree(errorsF32);
     }
-    psTrace(__func__, 3,"---- %s(psStats) end  ----\n", __func__);
+    psTrace(__func__, 3,"---- %s() end  ----\n", __func__);
     return (stats);
 }

trunk/psLib/src/math/psStats.h

@@ -14 +14 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-08-30 01:14:13 $
+ *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-26 23:49:11 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -37 +37 @@
 // XXX: Is PS_STAT_ROBUST_FOR_SAMPLE obsolete?
 typedef enum {
-    PS_STAT_SAMPLE_MEAN = 0x000001,         ///< Sample Mean
-    PS_STAT_SAMPLE_MEDIAN = 0x000002,       ///< Sample Median
-    PS_STAT_SAMPLE_STDEV = 0x000004,        ///< Sample Standard Deviation
-    PS_STAT_SAMPLE_QUARTILE = 0x000008,     ///< Sample Quartile
-    PS_STAT_ROBUST_MEAN = 0x000010,         ///< Robust Mean
-    PS_STAT_ROBUST_MEDIAN = 0x000020,       ///< Robust Median
-    PS_STAT_ROBUST_MODE = 0x000040,         ///< Robust Mode
-    PS_STAT_ROBUST_STDEV = 0x000080,        ///< Robust Standarad Deviation
-    PS_STAT_ROBUST_QUARTILE = 0x000100,     ///< Robust Quartile
-    PS_STAT_CLIPPED_MEAN = 0x000200,        ///< Clipped Mean
-    PS_STAT_CLIPPED_STDEV = 0x000400,       ///< Clipped Standard Deviation
-    PS_STAT_MAX =  0x000800,                ///< Maximum
-    PS_STAT_MIN =  0x001000,                ///< Minumum
-    PS_STAT_USE_RANGE =  0x002000,          ///< Range
-    PS_STAT_USE_BINSIZE = 0x004000,         ///< Binsize
+    PS_STAT_SAMPLE_MEAN     = 0x000001, ///< Sample Mean
+    PS_STAT_SAMPLE_MEDIAN   = 0x000002, ///< Sample Median
+    PS_STAT_SAMPLE_STDEV    = 0x000004, ///< Sample Standard Deviation
+    PS_STAT_SAMPLE_QUARTILE = 0x000008, ///< Sample Quartile
+    PS_STAT_ROBUST_MEDIAN   = 0x000010, ///< Robust Median
+    PS_STAT_ROBUST_STDEV    = 0x000020, ///< Robust Standarad Deviation
+    PS_STAT_ROBUST_QUARTILE = 0x000040, ///< Robust Quartile
+    PS_STAT_FITTED_MEAN     = 0x000080, ///< Fitted Mean
+    PS_STAT_FITTED_STDEV    = 0x000100, ///< Fitted Standard Deviation
+    PS_STAT_CLIPPED_MEAN    = 0x000200, ///< Clipped Mean
+    PS_STAT_CLIPPED_STDEV   = 0x000400, ///< Clipped Standard Deviation
+    PS_STAT_MAX             = 0x000800, ///< Maximum
+    PS_STAT_MIN             = 0x001000, ///< Minumum
+    PS_STAT_USE_RANGE       = 0x002000, ///< Range
+    PS_STAT_USE_BINSIZE     = 0x004000, ///< Binsize
 } psStatsOptions;
 
@@ -62 +62 @@
     double sampleMedian;               ///< formal median of sample
     double sampleStdev;                ///< standard deviation of sample
+    double sampleLQ;                   ///< lower quartile of sample
     double sampleUQ;                   ///< upper quartile of sample
-    double sampleLQ;                   ///< lower quartile of sample
-    double robustMean;                 ///< robust mean of array
     double robustMedian;               ///< robust median of array
-    double robustMode;                 ///< Robust mode of array
     double robustStdev;                ///< robust standard deviation of array
     double robustUQ;                   ///< robust upper quartile
     double robustLQ;                   ///< robust lower quartile
     int robustN50;                     ///< Number of points in Gaussian fit.  XXX: This is currently never set.
-    int robustNfit;                    ///< The number of points between the quartiles.
+    double fittedMean;                 ///< robust mean of data
+    double fittedStdev;                ///< robust standard deviation of data
+    int fittedNfit;                    ///< Number of points in Gaussian fit
     double clippedMean;                ///< Nsigma clipped mean
     double clippedStdev;               ///< standard deviation after clipping
@@ -80 +80 @@
     double max;                        ///< maximum data value in array
     double binsize;                    ///< binsize for robust fit (input/ouput)
+    int nSubsample;                    ///< maxinum number of measurements (input)
     psStatsOptions options;            ///< bitmask of calculated values
 }

trunk/psLib/test/math/tst_psStats00.c

@@ -6 +6 @@
 *  do a much simpler test with data types PS_TYPE_U8, PS_TYPE_U16, PS_TYPE_F64.
 *
-*  If the psStats,c code every changes such that vectors of different type
+*  If the psStats,c code ever changes such that vectors of different type
 *  are handled by different routines, then these tests must be extended.
 *
 *  @author GLG, MHPCC
 * 
-*  @version $Revision: 1.3 $  $Name: not supported by cvs2svn $
-*  @date $Date: 2005-09-24 01:10:00 $
+*  @version $Revision: 1.4 $  $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-26 23:49:11 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii
@@ -22 +22 @@
 #define ERROR_TOL  0.0001
 #define N 15
-
+#define VERBOSE 1
 static psS32 testStatsSampleMeanF32(void);
 static psS32 testStatsSampleMeanS8(void);
@@ -71 +71 @@
     psTraceSetLevel(".", 0);
     psTraceSetLevel("p_psGetStatValue", 0);
-    psTraceSetLevel("p_psVectorSampleMean", 0);
+    psTraceSetLevel("p_psVectorSampleMean", 10);
     psTraceSetLevel("p_psVectorMax", 0);
     psTraceSetLevel("p_psVectorMin", 0);
@@ -119 +119 @@
     }
 
+    printf("HMMM: myStats->options is (0x%x)\n", myStats->options);
     /*************************************************************************/
     /*  Call psVectorStats() with no vector mask.                    */
@@ -125 +126 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
-    if ( fabs(mean - expectedMeanNoMaskF32) > ERROR_TOL ) {
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskF32);
+    }
+    if (isnan(myStats->sampleMean) || (fabs(mean - expectedMeanNoMaskF32) > ERROR_TOL)) {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",
                 mean, expectedMeanNoMaskF32);
@@ -135 +139 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
-    if ( fabs(mean - expectedWeightMeanNoMaskF32) > ERROR_TOL ) {
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanNoMaskF32);
+    }
+    if (isnan(myStats->sampleMean) || (fabs(mean - expectedWeightMeanNoMaskF32) > ERROR_TOL)) {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",
                 mean, expectedWeightMeanNoMaskF32);
@@ -148 +155 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanRangeNoMaskF32);
+    }
     if ( fabs(mean - expectedMeanRangeNoMaskF32) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f",
@@ -158 +168 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanNoMaskRangeF32);
+    }
     if ( fabs(mean - expectedWeightMeanNoMaskRangeF32) > ERROR_TOL) {
         psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f",
@@ -170 +183 @@
     myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1);
     mean = myStats->sampleMean;
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanWithMaskF32);
+    }
     if ( fabs(mean - expectedMeanWithMaskF32) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",
@@ -179 +195 @@
     myStats = psVectorStats(myStats, myVector, myErrors, maskVector, 1);
     mean = myStats->sampleMean;
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanWithMaskF32);
+    }
     if ( fabs(mean - expectedWeightMeanWithMaskF32) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",
@@ -190 +209 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanRangeWithMaskF32);
+    }
     if ( fabs(mean - expectedMeanRangeWithMaskF32) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f",
@@ -200 +222 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedWeightMeanWithMaskRangeF32);
+    }
     if ( fabs(mean - expectedWeightMeanWithMaskRangeF32) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f",
@@ -219 +244 @@
     myStats = psVectorStats(myStats, myVector, NULL, maskVector, 2);
     mean = myStats->sampleMean;
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanWithMaskF32);
+    }
     if (fabs(mean - expectedMeanWithMaskF32) > ERROR_TOL )  {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",
@@ -237 +265 @@
     myStats = psVectorStats(myStats, myVector, NULL, maskVector, 4);
     mean = myStats->sampleMean;
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskF32);
+    }
     if (fabs(mean - expectedMeanNoMaskF32) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Return value %f not as expected %f",
@@ -307 +338 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskS8);
+    }
     if ( fabs(mean - expectedMeanNoMaskS8) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",
@@ -347 +381 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskU16);
+    }
     if ( fabs(mean - expectedMeanNoMaskU16) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",
@@ -387 +424 @@
     mean = myStats->sampleMean;
     // Verify return value is as expected
+    if (VERBOSE) {
+        printf("psVectorStats() returned %.2f: expected was %.2f\n", mean, expectedMeanNoMaskF64);
+    }
     if ( fabs(mean - expectedMeanNoMaskF64) > ERROR_TOL ) {
         psError(PS_ERR_UNKNOWN,true,"Returned value %f not as expected %f",

trunk/psLib/test/math/tst_psStats03.c

@@ -18 +18 @@
     // We list pertinent psStats.c functions here for debugging ease.
     //
-    psTraceSetLevel(".", 0);
-    psTraceSetLevel("p_psGetStatValue", 0);
-    psTraceSetLevel("p_psVectorSampleMean", 0);
-    psTraceSetLevel("p_psVectorCheckNonEmpty", 0);
-    psTraceSetLevel("p_psVectorNValues", 0);
-    psTraceSetLevel("p_psVectorSampleMedian", 0);
-    psTraceSetLevel("psStatsAlloc", 0);
-    psTraceSetLevel("p_psConvertToF32", 0);
-    psTraceSetLevel("psVectorStats", 0);
+    psTraceSetLevel(".", 10);
+    psTraceSetLevel("p_psGetStatValue", 10);
+    psTraceSetLevel("p_psVectorSampleMean", 10);
+    psTraceSetLevel("p_psVectorCheckNonEmpty", 10);
+    psTraceSetLevel("p_psVectorNValues", 10);
+    psTraceSetLevel("p_psVectorSampleMedian", 10);
+    psTraceSetLevel("psStatsAlloc", 10);
+    psTraceSetLevel("p_psConvertToF32", 10);
+    psTraceSetLevel("psVectorStats", 10);
 
     psStats *myStats    = NULL;

trunk/psLib/test/math/tst_psStats07.c

@@ -40 +40 @@
     float realUQNoMask = MEAN + ( 0.6 * STDEV );
     psS32 realN50NoMask = N / 4;
-    psS32 realNfitNoMask = N / 4;
 
     /*************************************************************************/
@@ -79 +78 @@
     psFree(mySampleStatsWithMask);
 
-    myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN |
-                            PS_STAT_ROBUST_MEDIAN |
-                            PS_STAT_ROBUST_MODE |
-                            PS_STAT_ROBUST_STDEV |
-                            PS_STAT_ROBUST_QUARTILE );
+    myStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN |
+                           PS_STAT_ROBUST_STDEV |
+                           PS_STAT_ROBUST_QUARTILE |
+                           PS_STAT_FITTED_MEAN |
+                           PS_STAT_FITTED_STDEV);
     // Create a full outliers:
     myVector->data.F32[N/4] = -1000.0 * MEAN;
@@ -101 +100 @@
 
     printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n",
-            realMeanNoMask, myStats->robustMean );
-
-    if ( fabs( myStats->robustMean - realMeanNoMask ) < ( ERROR_TOLERANCE * realMeanNoMask ) ) {
+            realMeanNoMask, myStats->fittedMean );
+
+    if ( fabs( myStats->fittedMean - realMeanNoMask ) < ( ERROR_TOLERANCE * realMeanNoMask ) ) {
         testStatus = true;
     } else {
@@ -157 +156 @@
 
     printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n",
-            realStdevNoMask, myStats->robustStdev );
-    if ( fabs( myStats->robustStdev - realStdevNoMask ) < ( ERROR_TOLERANCE * realStdevNoMask ) ) {
+            realStdevNoMask, myStats->fittedStdev );
+    if ( fabs( myStats->fittedStdev - realStdevNoMask ) < ( ERROR_TOLERANCE * realStdevNoMask ) ) {
         testStatus = true;
     } else {
@@ -208 +207 @@
                              "psStats functions",
                              "PS_STAT_ROBUST_STATS: robust N50: no vector mask" );
-
     // XXX:
     realN50NoMask = myStats->robustN50;
@@ -227 +225 @@
                  testStatus );
 
-    printPositiveTestHeader( stdout,
-                             "psStats functions",
-                             "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask" );
-
-    // XXX:
-    realNfitNoMask = myStats->robustNfit;
-
-    printf( "The expected Nfit was %d; the calculated Nfit was %d\n",
-            realNfitNoMask, myStats->robustNfit );
-    /* XXX: fix
-        if ( fabs( myStats->robustNfit - realNfitNoMask ) < ( ERROR_TOLERANCE * realNfitNoMask ) ) {
-            testStatus = true;
-        } else {
-            testStatus = false;
-            globalTestStatus = false;
-        }
-    */
-    printFooter( stdout,
-                 "psVector functions",
-                 "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask",
-                 testStatus );
 
     /*************************************************************************/
@@ -296 +273 @@
     float realUQWithMask = MEAN;
     psS32 realN50WithMask = N / 4;
-    psS32 realNfitWithMask = N / 4;
     /*************************************************************************/
     /*  Allocate and initialize data structures                              */
     /*************************************************************************/
-    myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN |
-                            PS_STAT_ROBUST_MEDIAN |
-                            PS_STAT_ROBUST_MODE |
-                            PS_STAT_ROBUST_STDEV |
-                            PS_STAT_ROBUST_QUARTILE );
+    myStats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN |
+                           PS_STAT_ROBUST_STDEV |
+                           PS_STAT_ROBUST_QUARTILE |
+                           PS_STAT_FITTED_MEAN |
+                           PS_STAT_FITTED_STDEV);
 
     maskVector = psVectorAlloc( N, PS_TYPE_U8 );
@@ -334 +310 @@
     printf( "Called psVectorStats() on a vector with elements masked.\n" );
     printf( "The expected Mean was %.2f; the calculated Mean was %.2f\n",
-            realMeanWithMask, myStats->robustMean );
-    if ( fabs( myStats->robustMean - realMeanWithMask ) < ( ERROR_TOLERANCE * realMeanWithMask ) ) {
+            realMeanWithMask, myStats->fittedMean );
+    if ( fabs( myStats->fittedMean - realMeanWithMask ) < ( ERROR_TOLERANCE * realMeanWithMask ) ) {
         testStatus = true;
     } else {
@@ -389 +365 @@
 
     printf( "The expected Stdev was %.2f; the calculated Stdev was %.2f\n",
-            realStdevWithMask, myStats->robustStdev );
-    if ( fabs( myStats->robustStdev - realStdevWithMask ) < ( ERROR_TOLERANCE * realStdevWithMask ) ) {
+            realStdevWithMask, myStats->fittedStdev );
+    if ( fabs( myStats->fittedStdev - realStdevWithMask ) < ( ERROR_TOLERANCE * realStdevWithMask ) ) {
         testStatus = true;
     } else {
@@ -459 +435 @@
                  "PS_STAT_ROBUST_STATS: robust N50: with vector mask",
                  testStatus );
-
-
-
-    printPositiveTestHeader( stdout,
-                             "psStats functions",
-                             "PS_STAT_ROBUST_STATS: robust Nfit: with vector mask" );
-
-    printf( "The expected Nfit was %d; the calculated Nfit was %d\n",
-            realNfitWithMask, myStats->robustNfit );
-    /* XXX: fix
-        if ( fabs( myStats->robustNfit - realNfitWithMask ) < ( ERROR_TOLERANCE * realNfitWithMask ) ) {
-            testStatus = true;
-        } else {
-            testStatus = false;
-            globalTestStatus = false;
-        }
-    */
-    printFooter( stdout,
-                 "psVector functions",
-                 "PS_STAT_ROBUST_STATS: robust Nfit: with vector mask",
-                 testStatus );
-