Changeset 3349 for trunk/psLib/src/math/psStats.c

Timestamp:

Mar 1, 2005, 11:03:25 AM (21 years ago)

Author:

gusciora

Message:

Fixing bugs from 316. Not tested yet.

File:

• trunk/psLib/src/math/psStats.c (modified) (10 diffs)

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.114 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-03-01 21:02:21 $
+ *  @version $Revision: 1.115 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-03-01 21:03:25 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -1002 +1002 @@
     -1: error
     -2: warning
+ 
+XXX: Create a new psStats object for the sample mean and sample stdev.
+ 
+XXX: Do we really need to calculate median and mean?
+ 
+XXX: Use static vectors for tmpMask.
  *****************************************************************************/
 psS32 p_psVectorClippedStats(const psVector* myVector,
@@ -1013 +1019 @@
     psF32 clippedMean = 0.0;    // self-explanatory
     psF32 clippedStdev = 0.0;   // self-explanatory
-    psF32 oldStanMean = 0.0;    // Temporary variable
-    psF32 oldStanStdev = 0.0;   // Temporary variable
     psVector* tmpMask = NULL;   // Temporary vector
+    psStats *statsTmp = NULL;
 
     // Ensure that stats->clipIter is within the proper range.
@@ -1026 +1031 @@
                        PS_CLIPPED_SIGMA_LB,
                        PS_CLIPPED_SIGMA_UB, -1);
+
+    if (statsTmp == NULL) {
+        statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
+        p_psMemSetPersistent(statsTmp, true);
+    }
 
     // We allocate a temporary mask vector since during the iterative
@@ -1040 +1050 @@
     }
     // 1. Compute the sample median.
-    p_psVectorSampleMedian(myVector, maskVector, maskVal, stats);
-    if (isnan(stats->sampleMedian)) {
+    p_psVectorSampleMedian(myVector, maskVector, maskVal, statsTmp);
+    if (isnan(statsTmp->sampleMedian)) {
         psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleMedian returned NAN\n");
         stats->clippedMean = NAN;
@@ -1049 +1059 @@
 
     // 2. Compute the sample standard deviation.
-    p_psVectorSampleStdev(myVector, errors, maskVector, maskVal, stats);
-    if (isnan(stats->sampleStdev)) {
+    p_psVectorSampleStdev(myVector, errors, maskVector, maskVal, statsTmp);
+    if (isnan(statsTmp->sampleStdev)) {
         psLogMsg(__func__, PS_LOG_WARN, "Call to p_psVectorSampleStdev returned NAN\n");
         stats->clippedMean = NAN;
@@ -1058 +1068 @@
 
     // 3. Use the sample median as the first estimator of the mean X.
-    clippedMean = stats->sampleMean;
-
+    clippedMean = statsTmp->sampleMean;
     // 4. Use the sample stdev as the first estimator of the mean stdev.
-    clippedStdev = stats->sampleStdev;
-
-    // Must save the old sampleMean and sampleStdev since the following code
-    // block overwrites them.
-    oldStanMean = stats->sampleMean;
-    oldStanStdev = stats->sampleStdev;
+    clippedStdev = statsTmp->sampleStdev;
 
     // 5. Repeat N times:
@@ -1073 +1077 @@
         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;
+                if (fabs(myVector->data.F32[j] - clippedMean) >
+                        (stats->clipSigma * errors->data.F32[j])) {
+                    tmpMask->data.U8[j] = 0xff;
                 }
             }
         } else {
             for (j = 0; j < myVector->n; j++) {
-                if (fabs(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) {
-                    tmpMask->data.U8[i] = 0xff;
+                if (fabs(myVector->data.F32[j] - clippedMean) >
+                        (stats->clipSigma * clippedStdev)) {
+                    tmpMask->data.U8[j] = 0xff;
                 }
             }
@@ -1086 +1092 @@
 
         // b) compute new mean and stdev
-        p_psVectorSampleMean(myVector, errors, tmpMask, maskVal, stats);
-        p_psVectorSampleStdev(myVector, errors, tmpMask, maskVal, stats);
+        p_psVectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp);
+        p_psVectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp);
 
         // If the new mean and stdev are not NAN, then we use them.
         // Otherwise, keep the old ones, and exit the loop.
-        if (! (isnan(stats->sampleMean) || isnan(stats->sampleStdev))) {
-            clippedMean = stats->sampleMean;
-            clippedStdev = stats->sampleStdev;
-
+        if (! (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev))) {
+            clippedMean = statsTmp->sampleMean;
+            clippedStdev = statsTmp->sampleStdev;
+        } else {
             // Exit loop
+            // Should we throw an error or warning here?
             i = stats->clipIter;
         }
     }
-
-    stats->sampleMean = oldStanMean;
-    stats->sampleStdev = oldStanStdev;
 
     // 7. The last calcuated value of x is the cliped mean.
@@ -2240 +2244 @@
     // 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, errorsF32, mask, maskVal, stats)) {
+        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");
             stats->clippedMean = NAN;
             stats->clippedStdev = NAN;
-        } else if (-2 == p_psVectorClippedStats(inF32, errorsF32, mask, maskVal, stats)) {
+        } else if (-2 == rc) {
             psLogMsg(__func__, PS_LOG_WARN, "Failed to calculate clipped statistics for input psVector.");
             stats->clippedMean = NAN;