Index: /branches/eam_rel8_b1/psLib/src/math/psMinimize.c
===================================================================
--- /branches/eam_rel8_b1/psLib/src/math/psMinimize.c	(revision 5566)
+++ /branches/eam_rel8_b1/psLib/src/math/psMinimize.c	(revision 5567)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.142 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-10-12 21:02:20 $
+ *  @version $Revision: 1.142.6.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-11-22 18:06:10 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -1443,4 +1443,66 @@
 
 /******************************************************************************
+BuildSums4D(sums, x, y, z, t, nXterm, nYterm, nZterm, nTterm). equiv to 
+BuildSums2D(). The result is returned as a double ****
+ *****************************************************************************/
+static psImage *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*PSELEMTYPE_SIZEOF(psF64));
+        for (int i = 0; i < nXsum; i++) {
+            sums[i] = (psF64 ***) psAlloc (nYsum*PSELEMTYPE_SIZEOF(psF64));
+            for (int j = 0; j < nYsum; j++) {
+                sums[i][j] = (psF64 **) psAlloc (nZsum*PSELEMTYPE_SIZEOF(psF64));
+                for (int k = 0; k < nZsum; k++) {
+                    sums[i][j][k] = (psF64 *) psAlloc (nTsum*PSELEMTYPE_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);
+}
+
+/******************************************************************************
 Polynomial2DEvalVectorD(myPoly, x, y): This routine takes as input two
 psVectors x and y, and evaluates myPoly for each pair of (x, y), and stores it
@@ -1658,7 +1720,10 @@
     // 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,6 +1735,6 @@
             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,29 +1751,11 @@
     }
 
-    // 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,5 +1866,5 @@
 }
 
-
+// This function accepts F32 and F64 input vectors.
 psPolynomial1D *psVectorClipFitPolynomial1D(
     psPolynomial1D *poly,
@@ -1827,8 +1874,8 @@
     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,24 +1883,93 @@
     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);
+    } else {
+        x = xIn;
+    }
+
+    // clipping range defined by min and max and/or clipSigma
+    float minClipSigma;
+    float maxClipSigma;
+    if (finite(stats->max)) {
+        maxClipSigma = fabs(stats->clipSigma);
+    } else {
+        maxClipSigma = fabs(stats->max);
+    }
+    if (finite(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,6 +2045,8 @@
     // 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,8 +2054,6 @@
             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 : 1.0 / PS_SQR(fErr->data.F64[k]);
         }
 
@@ -1962,11 +2078,7 @@
 
     // 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++) {
@@ -2251,5 +2363,5 @@
     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,21 +2370,75 @@
     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 (finite(stats->max)) {
+        maxClipSigma = fabs(stats->clipSigma);
     } else {
-        // XXX: Use maskValue.
-        // XXX: Change argument order.
-        poly = RobustFit2D(poly, mask, x, y, f, fErr);
-    }
+        maxClipSigma = fabs(stats->max);
+    }
+    if (finite(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 = 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) {