Changeset 5624 for trunk/psLib/src/math/psMinimize.c

Timestamp:

Nov 29, 2005, 4:00:11 PM (21 years ago)

Author:

desonia

Message:

merged eam_r8.0_b2 into CVS trunk.

File:

• trunk/psLib/src/math/psMinimize.c (modified) (15 diffs)

trunk/psLib/src/math/psMinimize.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.145 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-11-16 23:06:19 $
+ *  @version $Revision: 1.146 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-30 02:00:09 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -1441 +1441 @@
 }
 
+// here is the definition for BuildSums4D.  ifdef'ed away until it is used
+// by psPolynomial4DFit..
+# if (0)
+    /******************************************************************************
+    BuildSums4D(sums, x, y, z, t, nXterm, nYterm, nZterm, nTterm). equiv to
+    BuildSums2D(). The result is returned as a double ****
+     *****************************************************************************/
+    static double ****BuildSums4D(
+        psF64 ****sums,
+        psF64 x,
+        psF64 y,
+        psF64 z,
+        psF64 t,
+        psS32 nXterm,
+        psS32 nYterm,
+        psS32 nZterm,
+        psS32 nTterm)
+{
+    psS32 nXsum = 0;
+    psS32 nYsum = 0;
+    psS32 nZsum = 0;
+    psS32 nTsum = 0;
+    psF64 xSum = 1.0;
+    psF64 ySum = 1.0;
+    psF64 zSum = 1.0;
+    psF64 tSum = 1.0;
+
+    nXsum = 2*nXterm;
+    nYsum = 2*nYterm;
+    nZsum = 2*nZterm;
+    nTsum = 2*nTterm;
+    if (sums == NULL) {
+        sums = (psF64 ****) psAlloc (nXsum*sizeof(psF64));
+        for (int i = 0; i < nXsum; i++) {
+            sums[i] = (psF64 ***) psAlloc (nYsum*sizeof(psF64));
+            for (int j = 0; j < nYsum; j++) {
+                sums[i][j] = (psF64 **) psAlloc (nZsum*sizeof(psF64));
+                for (int k = 0; k < nZsum; k++) {
+                    sums[i][j][k] = (psF64 *) psAlloc (nTsum*sizeof(psF64));
+                }
+            }
+        }
+    }
+    // careful with this function: there is no size checking and realloc for reuse
+
+    tSum = 1.0;
+    for (int m = 0; m < nTsum; m++) {
+        zSum = tSum;
+        for (int k = 0; k < nZsum; k++) {
+            ySum = zSum;
+            for (int j = 0; j < nYsum; j++) {
+                xSum = ySum;
+                for (int i = 0; i < nXsum; i++) {
+                    sums[i][j][k][m] = xSum;
+                    xSum *= x;
+                }
+                ySum *= y;
+            }
+            zSum *= z;
+        }
+        tSum *= t;
+    }
+    return (sums);
+}
+# endif /* BuildSums4D */
 
 /******************************************************************************
@@ -1447 +1512 @@
 in the output vector.  This routine works on single-precision polynomials with
 double precision data.
+ 
+XXX EAM : this function is now deprecated: psPolynomial2DEvalVector handles F32 and F64
  *****************************************************************************/
 psVector *Polynomial2DEvalVectorD(
@@ -1658 +1725 @@
     // xSums look like: 1, x, x^2, ... x^(2n+1)
     // Build the B and A data structs.
+    // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects
+    // XXX EAM : this function is only valid for data vectors of F64
     for (int k = 0; k < f->n; k++) {
-        if ((mask != NULL) && mask->data.U8[k])
+        if ((mask != NULL) && (mask->data.U8[k] && maskValue)) {
             continue;
+        }
         if (x != NULL) {
             xSums = BuildSums1D(xSums, x->data.F64[k], nTerm);
@@ -1670 +1740 @@
             wt = 1.0;
         } else {
-            // this should filter fErr == 0 values
-            wt = 1.0 / PS_SQR(fErr->data.F64[k]);
+            // this filters fErr == 0 values
+            wt = (fErr->data.F64[k] == 0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]);
         }
         for (int i = 0; i < nTerm; i++) {
@@ -1686 +1756 @@
     }
 
-    // GaussJordan version
-    if (0) {
-        // does the solution in place
-        psGaussJordan (A, B);
-
-        // the first nTerm entries in B correspond directly to the desired
-        // polynomial coefficients.  this is only true for the 1D case
-        for (int k = 0; k < nTerm; k++) {
-            myPoly->coeff[k] = B->data.F64[k];
-        }
-    } else {
-        // LUD version of the fit
-        psImage *ALUD = NULL;
-        psVector* outPerm = NULL;
-        psVector* coeffs = NULL;
-
-        ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
-        ALUD = psMatrixLUD(ALUD, &outPerm, A);
-        coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-        for (int k = 0; k < nTerm; k++) {
-            myPoly->coeff[k] = coeffs->data.F64[k];
-        }
-        psFree(ALUD);
-        psFree(coeffs);
-        psFree(outPerm);
+    // does the solution in place
+    psGaussJordan (A, B);
+
+    // the first nTerm entries in B correspond directly to the desired
+    // polynomial coefficients.  this is only true for the 1D case
+    for (int k = 0; k < nTerm; k++) {
+        myPoly->coeff[k] = B->data.F64[k];
     }
 
@@ -1819 +1871 @@
 }
 
-
+// This function accepts F32 and F64 input vectors.
 psPolynomial1D *psVectorClipFitPolynomial1D(
     psPolynomial1D *poly,
@@ -1827 +1879 @@
     const psVector *f,
     const psVector *fErr,
-    const psVector *x)
+    const psVector *xIn)
 {
     // Internal pointers for possibly NULL vectors.
-    psVector *x32 = NULL;
+    psVector *x = NULL;
 
     PS_ASSERT_POLY_NON_NULL(poly, NULL);
@@ -1836 +1888 @@
     PS_ASSERT_PTR_NON_NULL(stats, NULL);
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
     if (mask != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
         PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
     }
-    if (x != NULL) {
-        PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL);
-        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    if (xIn != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, xIn, NULL);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(xIn, NULL);
     }
     if (fErr != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
-    }
-
-    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+    }
+
+    // assign sequence vector if xIn is NULL
+    if (xIn == NULL) {
+        x = psVectorCreate (NULL, 0, f->n, 1, f->type.type);
+    } else {
+        x = (psVector *) xIn;
+    }
+
+    // clipping range defined by min and max and/or clipSigma
+    float minClipSigma;
+    float maxClipSigma;
+    if (isfinite(stats->max)) {
+        maxClipSigma = fabs(stats->clipSigma);
+    } else {
+        maxClipSigma = fabs(stats->max);
+    }
+    if (isfinite(stats->min)) {
+        minClipSigma = fabs(stats->clipSigma);
+    } else {
+        minClipSigma = fabs(stats->min);
+    }
+    psVector *fit   = NULL;
+    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    // eventual expansion: user supplies one of various stats option pairs,
+    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
+    // evaluate the clipping sigma
+    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < stats->clipIter; N++) {
+        int Nkeep = 0;
+
+        poly = psVectorFitPolynomial1D (poly, mask, maskValue, f, fErr, x);
+        fit = psPolynomial1DEvalVector (poly, x);
+        resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit);
+
+        stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
+        float minClipValue = -minClipSigma*stats->sampleStdev;
+        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+
+        // set mask if pts are not valid
+        // we are masking out any point which is out of range
+        // recovery is not allowed with this scheme
+        for (int i = 0; i < resid->n; i++) {
+            if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
+                continue;
+            }
+            if (resid->data.F64[i] - stats->sampleMedian > maxClipValue) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            if (resid->data.F64[i] - stats->sampleMedian < minClipValue) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            Nkeep ++;
+        }
+
+        psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n",
+                 Nkeep, x->n);
+
+        psFree (fit);
+    }
     // Free psVectors that were created for NULL arguments.
-    if (x == NULL) {
-        psFree(x32);
-    }
-    return(NULL);
+    if (xIn == NULL) {
+        psFree(x);
+    }
+    // Free other local temporary variables
+    psFree (resid);
+
+    return (poly);
 }
 
@@ -1929 +2050 @@
     // Build the B and A data structs.
     for (int k  = 0; k < x->n; k++) {
-        if ((mask != NULL) && mask->data.U8[k])
+        if ((mask != NULL) && (mask->data.U8[k] & maskValue)) {
             continue;
+        }
+
         Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm);
 
@@ -1936 +2059 @@
             wt = 1.0;
         } else {
-            // XXX: this should probably by fErr^2 !!
-            // this should filter fErr == 0 values
-            // XXX: Why isn't this fErr^2?
-            wt = 1.0 / fErr->data.F64[k];
+            // this filters fErr == 0 values
+            wt = (fErr->data.F64[k] == 0.0) ? 0.0 : 1.0 / PS_SQR(fErr->data.F64[k]);
         }
 
@@ -1962 +2083 @@
 
     // does the solution in place
-    // XXX: Check return codes!
     psGaussJordan (A, B);
 
-    // XXX: Check return codes!
-    // ALUD = psMatrixLUD(ALUD, &outPerm, A);
-    // coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-
+    // select the appropriate solution entries
     for (int n = 0; n < nXterm; n++) {
         for (int m = 0; m < nYterm; m++) {
@@ -1985 +2102 @@
 
 
+// XXX EAM : I have implemented a single function to handle the mask/nomask cases
+//           this function can be deprecated
 psPolynomial2D* RobustFit2D_nomask(
     psPolynomial2D* poly,
@@ -2251 +2370 @@
     PS_ASSERT_PTR_NON_NULL(stats, NULL);
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
-    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL);
     if (mask != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
@@ -2258 +2377 @@
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL);
     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
     PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
-    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
     if (fErr != NULL) {
         PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
-        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
-    }
-
-    if (mask == NULL) {
-        // XXX: Change argument order.
-        poly = RobustFit2D_nomask(poly, x, y, f, fErr);
+        PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL);
+    }
+
+    // clipping range defined by min and max and/or clipSigma
+    float minClipSigma;
+    float maxClipSigma;
+    if (isfinite(stats->max)) {
+        maxClipSigma = fabs(stats->max);
     } else {
-        // XXX: Use maskValue.
-        // XXX: Change argument order.
-        poly = RobustFit2D(poly, mask, x, y, f, fErr);
-    }
+        maxClipSigma = fabs(stats->clipSigma);
+    }
+    if (isfinite(stats->min)) {
+        minClipSigma = fabs(stats->min);
+    } else {
+        minClipSigma = fabs(stats->clipSigma);
+    }
+    psVector *fit   = NULL;
+    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    // eventual expansion: user supplies one of various stats option pairs,
+    // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to
+    // evaluate the clipping sigma
+    // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used
+    stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    for (int N = 0; N < stats->clipIter; N++) {
+        int Nkeep = 0;
+
+        poly = psVectorFitPolynomial2D (poly, mask, maskValue, f, fErr, x, y);
+        fit = psPolynomial2DEvalVector (poly, x, y);
+        resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit);
+
+        stats  = psVectorStats (stats, resid, NULL, mask, maskValue);
+        float minClipValue = -minClipSigma*stats->sampleStdev;
+        float maxClipValue = +maxClipSigma*stats->sampleStdev;
+
+        // set mask if pts are not valid
+        // we are masking out any point which is out of range
+        // recovery is not allowed with this scheme
+        for (int i = 0; i < resid->n; i++) {
+            if ((mask != NULL) && (mask->data.U8[i] & maskValue)) {
+                continue;
+            }
+            if (resid->data.F64[i] - stats->sampleMedian > maxClipValue) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            if (resid->data.F64[i] - stats->sampleMedian < minClipValue) {
+                if (mask != NULL) {
+                    mask->data.U8[i] |= 0x01;
+                }
+                continue;
+            }
+            Nkeep ++;
+        }
+
+        psTrace (".psphot.VectorClipFit", 4, "keeping %d of %d pts for fit\n",
+                 Nkeep, x->n);
+
+        psFree (fit);
+    }
+    // Free local temporary variables
+    psFree (resid);
 
     if (poly == NULL) {