Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 4944)
+++ trunk/psLib/src/math/psMinimize.c	(revision 4958)
@@ -9,6 +9,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.133 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-02 21:32:06 $
+ *  @version $Revision: 1.134 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-07 21:35:50 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -29,4 +29,5 @@
 #include "psImageStructManip.h"
 #include "psBinaryOp.h"
+#include "psLogMsg.h"
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
@@ -112,9 +113,11 @@
 
 /******************************************************************************
-p_psBuildSums1D(x, polyOrder, sums): this routine calculates the powers of
-input parameter "x" between 0 and input parameter polyOrder.  The result is
-returned as a psVector sums.
+buildSums1D(x, polyOrder, sums): this routine calculates the powers of input
+parameter "x" between 0 and input parameter polyOrder.  The result is returned
+as a psVector sums.
  
 XXX: Use a static vector.
+ 
+XXX: should the argument be polyOrder or numTerms?
  *****************************************************************************/
 static void buildSums1D(psF64 x,
@@ -124,14 +127,14 @@
     psS32 i = 0;
     psF64 xSum = 0.0;
+    psS32 numTerms = polyOrder + 1;
 
     if (sums == NULL) {
-        sums = psVectorAlloc(polyOrder, PS_TYPE_F64);
-    }
-    if (polyOrder > sums->n) {
-        sums = psVectorRealloc(sums, polyOrder);
+        sums = psVectorAlloc(numTerms, PS_TYPE_F64);
+    } else if (numTerms > sums->n) {
+        sums = psVectorRealloc(sums, numTerms);
     }
 
     xSum = 1.0;
-    for (i = 0; i <= polyOrder; i++) {
+    for (i = 0; i < numTerms; i++) {
         sums->data.F64[i] = xSum;
         xSum *= x;
@@ -708,5 +711,5 @@
 those coords.
  
-// XXX EAM this is my re-implementation of MinLM
+XXX EAM this is my re-implementation of MinLM
  
 XXX: This must work for both F32 and F64.  F32 is currently implemented.
@@ -1306,4 +1309,6 @@
  
 XXX: type F32 is done via vector conversion only.
+ 
+XXX: Get rid of this.  Use new argument list.
  *****************************************************************************/
 psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* poly,
@@ -2143,8 +2148,13 @@
 
 /******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
 EAM Code:
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
  *****************************************************************************/
 
-// XXX EAM : my alternate EAMBuildSums1D
 static psVector *EAMBuildSums1D(
     psVector* sums,
@@ -2155,9 +2165,12 @@
     psF64 xSum = 0.0;
 
+    //
+    // XXX: Why do we multiply by 2 here?  Better to do it outside and have the
+    // definition of this function remain sensible.
+    //
     nSum = 2*nTerm;
     if (sums == NULL) {
         sums = psVectorAlloc(nSum, PS_TYPE_F64);
-    }
-    if (nSum > sums->n) {
+    } else if (nSum > sums->n) {
         sums = psVectorRealloc(sums, nSum);
     }
@@ -2177,6 +2190,5 @@
     const psVector *x,
     const psVector *y,
-    const psVector *yErr
-)
+    const psVector *yErr)
 {
     // I think this is 1 dimension down
@@ -2277,6 +2289,5 @@
     const psPolynomial2D *myPoly,
     const psVector *x,
-    const psVector *y
-)
+    const psVector *y)
 {
     PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
@@ -2312,6 +2323,5 @@
     psF64 y,
     psS32 nXterm,
-    psS32 nYterm
-)
+    psS32 nYterm)
 {
     psS32 nXsum = 0;
@@ -2349,6 +2359,5 @@
     const psVector* y,
     const psVector* z,
-    const psVector* zErr
-)
+    const psVector* zErr)
 {
     // I think this is 1 dimension down
@@ -2439,6 +2448,5 @@
     const psVector* y,
     const psVector* z,
-    const psVector* dz
-)
+    const psVector* dz)
 {
     psVector *X;
@@ -2487,4 +2495,15 @@
 
 // XXX EAM : be careful here with F32 vs F64 vectors
+/*
+ Basically, you repetitively fit a polynomial to a set of data points,
+ reject the points that did not fit well, then refit the polynomial.
+ 
+ Basically, simply fit the polynomial to the data.  They compare the
+ fit (by evaluating the x data with that polynomial and subtracting
+ from the original f data).  That's the residual.  Loop through all
+ data and if the ((residual - mean) > 3*stDev), mask that data point,
+ and fit the polynomial again.  Do this 3 times.
+*/
+
 psPolynomial2D* RobustFit2D(psPolynomial2D* poly,
                             psVector* mask,
@@ -2531,3 +2550,412 @@
 }
 
-
+/******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+NEW Code:
+ ******************************************************************************
+ ******************************************************************************
+ ******************************************************************************
+ *****************************************************************************/
+
+#define PS_VECTOR_GEN_INDEX_F32(VEC, SIZE) \
+VEC = psVectorAlloc(SIZE, PS_TYPE_F32); \
+for (psS32 i = 0 ; i < SIZE ; i++) { \
+    VEC->data.F32[i] = (psF32) i; \
+} \
+
+psPolynomial1D *psVectorFitPolynomial1D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *x32 = NULL;
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, 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);
+    } else {
+        PS_VECTOR_GEN_INDEX_F32(x32, f->n);
+    }
+
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        // XXX: Call EAM code
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        // XXX: Call my code
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+
+    // Free psVectors that were created for NULL arguments.
+    if (x == NULL) {
+        psFree(x32);
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+
+    return(NULL);
+}
+
+
+psPolynomial2D *psVectorFitPolynomial2D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    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_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial3D *psVectorFitPolynomial3D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    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_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial4D *psVectorFitPolynomial4D_NEW(
+    psPolynomial1D *poly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    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_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    if (poly->type == PS_POLYNOMIAL_ORD) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented.  Returning NULL.\n");
+    } else {
+        printf("XXX: ERROR: incorrect polynomial type.\n");
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+
+psPolynomial1D *psVectorClipFitPolynomial1D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *x32 = NULL;
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, 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 (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    // Free psVectors that were created for NULL arguments.
+    if (x == NULL) {
+        psFree(x32);
+    }
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial2D *psVectorClipFitPolynomial2D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    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_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial3D *psVectorClipFitPolynomial3D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    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_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+psPolynomial4D *psVectorClipFitPolynomial4D_NEW(
+    psPolynomial1D *poly,
+    psStats *stats,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x,
+    const psVector *y,
+    const psVector *z,
+    const psVector *t)
+{
+    // Internal pointers for possibly NULL vectors.
+    psVector *fErr32 = NULL;
+
+    PS_ASSERT_POLY_NON_NULL(poly, NULL);
+    PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(f, NULL);
+    PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL);
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+    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_NON_NULL(y, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(z, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL);
+    PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(t, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL);
+    PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL);
+    if (fErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL);
+    } else {
+        fErr32 = psVectorAlloc(f->n, PS_TYPE_F32);
+        PS_VECTOR_SET_F32(fErr32, 1.0);
+    }
+
+    psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented.  Returning NULL.\n");
+    if (fErr == NULL) {
+        psFree(fErr32);
+    }
+    return(NULL);
+}
+
+