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

Timestamp:

Dec 21, 2004, 10:42:07 AM (22 years ago)

Author:

gusciora

Message:

Added errors arguments to the VectorStats and VectorHistogram functions.
The arguments are ignored in the histogram and robust stats code.

File:

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

trunk/psLib/src/math/psStats.c

@@ -9 +9 @@
  *  @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 +146 @@
  *****************************************************************************/
 int p_psVectorSampleMean(const psVector* restrict myVector,
+                         const psVector* restrict errors,
                          const psVector* restrict maskVector,
                          psU32 maskVal,
@@ -156 +157 @@
     // 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 +409 @@
 /******************************************************************************
 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 +507 @@
 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 +584 @@
     }
 
-    // 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 +786 @@
  
  *****************************************************************************/
-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 +803 @@
     // 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 +865 @@
     } 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 +995 @@
  *****************************************************************************/
 int p_psVectorClippedStats(const psVector* restrict myVector,
+                           const psVector* restrict errors,
                            const psVector* restrict maskVector,
                            psU32 maskVal,
@@ -873 +1040 @@
 
     // 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 +1061 @@
     // 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 +1417 @@
 *****************************************************************************/
 int p_psVectorRobustStats(const psVector* restrict myVector,
+                          const psVector* restrict errors,
                           const psVector* restrict maskVector,
                           psU32 maskVal,
@@ -1274 +1450 @@
     // 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 +1510 @@
 
     // Populate the histogram array.
-    psVectorHistogram(robustHistogram, myVector, maskVector, maskVal);
+    psVectorHistogram(robustHistogram, myVector, errors, maskVector, maskVal);
 
     // Smooth the histogram, Gaussian-style.
@@ -1416 +1593 @@
         }
     }
-    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 +1851 @@
 
 /*****************************************************************************
-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 +1863 @@
 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 +1880 @@
         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 +2012 @@
 psStats* psVectorStats(psStats* stats,
                        psVector* in,
+                       psVector* errors,
                        psVector* mask,
                        psU32 maskVal)
@@ -1837 +2022 @@
         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 +2049 @@
     // ************************************************************************
     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 +2063 @@
     // ************************************************************************
     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 +2084 @@
             (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 +2091 @@
     }
 
+    // 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 +2123 @@
     }
 
-    if (mustFreeTmp == 1) {
+    if (mustFreeVectorIn == 1) {
         psFree(inF32);
     }
+    if (mustFreeVectorErrors == 1) {
+        psFree(errorsF32);
+    }
     return (stats);
 }