Index: trunk/psLib/src/dataManip/psMinimize.c
===================================================================
--- trunk/psLib/src/dataManip/psMinimize.c	(revision 1629)
+++ trunk/psLib/src/dataManip/psMinimize.c	(revision 1719)
@@ -9,6 +9,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-08-26 20:10:25 $
+ *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-08 06:02:23 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -43,4 +43,5 @@
 #include "psMinimize.h"
 #include "psMatrix.h"
+#include "psTrace.h"
 
 /*****************************************************************************/
@@ -569,4 +570,7 @@
     gsl_vector *x;
 
+    psTrace(".psLib.dataManip.psMinimize", 4,
+            "---- psMinimize() begin ----\n");
+
     PS_CHECK_NULL_VECTOR(initialGuess);
     PS_CHECK_EMPTY_VECTOR(initialGuess);
@@ -655,4 +659,6 @@
         }
     }
+    psTrace(".psLib.dataManip.psMinimize", 4,
+            "---- psMinimize() end ----\n");
     return (initialGuess);
 }
@@ -691,4 +697,7 @@
     float chiSqOld = 0.0;
 
+    psTrace(".psLib.dataManip.psMinimizeChi2", 4,
+            "---- psMinimizeChi2() begin ----\n");
+
     PS_CHECK_NULL_IMAGE(domain);
     PS_CHECK_EMPTY_IMAGE(domain);
@@ -837,4 +846,6 @@
 
     // Bye bye.
+    psTrace(".psLib.dataManip.psMinimizeChi2", 4,
+            "---- psMinimizeChi2() end ----\n");
     return (initialGuess);
 }
@@ -863,4 +874,6 @@
     psVector* xSums = NULL;
 
+    psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
+            "---- psVectorFitPolynomial1D() begin ----\n");
     // printf("psVectorFitPolynomial1D()\n");
     // for (i=0;i<x->n;i++) {
@@ -936,12 +949,27 @@
     psFree(xSums);
 
+    psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4,
+            "---- psVectorFitPolynomial1D() begin ----\n");
     return (myPoly);
+}
+
+psMinimization *psMinimizationAlloc(int maxIter,
+                                    float tol)
+{
+    psMinimization *min = psAlloc(sizeof(psMinimization));
+    min->maxIter = maxIter;
+    min->tol = tol;
+    min->value = 0.0;
+    min->iter = 0;
+    min->lastDelta = 0.0;
+
+    return(min);
 }
 
 #define STEP_SIZE 0.10
 /******************************************************************************
-    This routine takes as input an arbitrary function of a single parameter.
-    This function produces as output a bracket [a, b, c] such that f(b) is
-    less than f(a) and f(b).
+    This routine takes as input an arbitrary function, and the parameter to
+    vary.  This function produces as output a bracket [a, b, c] such that
+    f(b) is less than f(a) and f(b).
  
     Algorithm: XXX completely ad hoc: start with the user-supplied starting
@@ -950,7 +978,8 @@
     found.
  *****************************************************************************/
-psVector *p_psDetermineBracket(psVector *params,
-                               const psArray *coords,
-                               psMinimizePowellFunc func)
+psVector *p_psDetermineBracketOld(psVector *params,
+                                  int dim,
+                                  const psArray *coords,
+                                  psMinimizePowellFunc func)
 {
     float a = 0.0;
@@ -966,16 +995,20 @@
     float new_cDir = 0.0;
     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
-    psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
-
+    float stepSize = params->data.F32[dim] * STEP_SIZE;
+    //    float initialParam = params->data.F32[dim];
+
+    if (stepSize == 0.0) {
+        stepSize = 1.0;
+    }
     a = b = c = params->data.F32[0];
-    a = a * (1.0 - STEP_SIZE);
-    c = c * (1.0 + STEP_SIZE);
-
-    dummyParams->data.F32[0] = a;
-    fa = func(dummyParams, coords);
-    dummyParams->data.F32[0] = b;
-    fb = func(dummyParams, coords);
-    dummyParams->data.F32[0] = c;
-    fc = func(dummyParams, coords);
+    a-= stepSize;
+    c+= stepSize;
+
+    params->data.F32[dim] = a;
+    fa = func(params, coords);
+    params->data.F32[dim] = b;
+    fb = func(params, coords);
+    params->data.F32[dim] = c;
+    fc = func(params, coords);
     if (fa < fb) {
         aDir = -1;
@@ -997,10 +1030,14 @@
             return(bracket);
         }
-        a = a * (1.0 - STEP_SIZE);
-        c = c * (1.0 + STEP_SIZE);
-        dummyParams->data.F32[0] = a;
-        fa = func(dummyParams, coords);
-        dummyParams->data.F32[0] = c;
-        fc = func(dummyParams, coords);
+        stepSize*= (1.0 + stepSize);
+        a = a - stepSize;
+        c = c + stepSize;
+
+        params->data.F32[dim] = a;
+        fa = func(params, coords);
+        params->data.F32[dim] = c;
+        fc = func(params, coords);
+
+        //printf("HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
 
         if (fa < fb) {
@@ -1017,57 +1054,251 @@
         if ((new_aDir == 1) && (aDir == -1)) {
             bracket->data.F32[0] = a;
-            bracket->data.F32[1] = bracket->data.F32[0] / (1.0 - STEP_SIZE);
-            bracket->data.F32[2] = bracket->data.F32[1] / (1.0 - STEP_SIZE);
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
             return(bracket);
         }
 
         if ((new_cDir == 1) && (cDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
             bracket->data.F32[2] = c;
-            bracket->data.F32[1] = bracket->data.F32[0] / (1.0 + STEP_SIZE);
-            bracket->data.F32[0] = bracket->data.F32[1] / (1.0 + STEP_SIZE);
             return(bracket);
         }
         aDir = new_aDir;
         cDir = new_cDir;
-    }
-    psFree(dummyParams);
+        iter--;
+    }
     psFree(bracket);
     return(NULL);
 }
 
-
-
 /******************************************************************************
-    This routine must minimize a function of exactly one parameter.
+    This routine takes as input an arbitrary function, and the parameter to
+    vary, and the line along which it must vary.  This function produces as
+    output a bracket [a, b, c] such that
+        f(param + b * line) is less than f(param + a * line) and
+                                         f(param + c * line).
+ 
+    Algorithm: XXX completely ad hoc: start with the user-supplied starting
+    parameter and call that b.  Calculate a/c as a fractional amount
+    smaller/larger than b.  Repeat this process until a local minimum is
+    found.
+ *****************************************************************************/
+psVector *p_psDetermineBracket(psVector *params,
+                               psVector *line,
+                               const psVector *paramMask,
+                               const psArray *coords,
+                               psMinimizePowellFunc func)
+{
+    float a = 0.0;
+    float b = 0.0;
+    float c = 0.0;
+    float fa = 0.0;
+    float fb = 0.0;
+    float fc = 0.0;
+    int iter = 100;
+    float aDir = 0.0;
+    float cDir = 0.0;
+    float new_aDir = 0.0;
+    float new_cDir = 0.0;
+    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
+    float stepSize = STEP_SIZE;
+    psVector *tmp = NULL;
+    int i = 0;
+    int null = 0;
+
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+            "---- p_psDetermineBracket() begin ----\n");
+
+    for (i=0;i<params->n;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            if (line->data.F32[i] >= FLT_EPSILON) {
+                null = 1;
+            }
+        }
+    }
+    if (null == 0) {
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 2,
+                "p_psDetermineBracket() called with zero line vector.\n");
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                "---- p_psDetermineBracket() end (NULL) ----\n");
+        return(NULL);
+    }
+
+    tmp = psVectorAlloc(params->n, PS_TYPE_F32);
+
+    b = 0;
+    a = -stepSize;
+    c = stepSize;
+
+    for (i=0;i<params->n;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
+        } else {
+            tmp->data.F32[i] = params->data.F32[i];
+        }
+    }
+    fa = func(tmp, coords);
+
+    for (i=0;i<params->n;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            tmp->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
+        } else {
+            tmp->data.F32[i] = params->data.F32[i];
+        }
+    }
+    fb = func(tmp, coords);
+
+    for (i=0;i<params->n;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
+        } else {
+            tmp->data.F32[i] = params->data.F32[i];
+        }
+    }
+    fc = func(tmp, coords);
+
+    if (fa < fb) {
+        aDir = -1;
+    } else {
+        aDir = 1;
+    }
+
+    if (fc < fb) {
+        cDir = -1;
+    } else {
+        cDir = 1;
+    }
+
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+            "HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
+    while (iter > 0) {
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                "psDetermineBracket(): iteration %d\n", iter);
+        if ((b < a) && (b < c)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+        stepSize*= (1.0 + stepSize);
+        a =- stepSize;
+        c =+ stepSize;
+
+        for (i=0;i<params->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
+            } else {
+                tmp->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fa = func(tmp, coords);
+
+        for (i=0;i<params->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
+            } else {
+                tmp->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fc = func(tmp, coords);
+
+        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
+                "HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
+
+        if (fa < fb) {
+            new_aDir = -1;
+        } else {
+            new_aDir = 1;
+        }
+
+        if (fc < fb) {
+            new_cDir = -1;
+        } else {
+            new_cDir = 1;
+        }
+        if ((new_aDir == 1) && (aDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+
+        if ((new_cDir == 1) && (cDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            psFree(tmp);
+            psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+                    "---- p_psDetermineBracket() end ----\n");
+            return(bracket);
+        }
+        aDir = new_aDir;
+        cDir = new_cDir;
+        iter--;
+    }
+    psFree(tmp);
+    psFree(bracket);
+    psTrace(".psLib.dataManip.p_psDetermineBracket", 4,
+            "---- p_psDetermineBracket() end (NULL) ----\n");
+    return(NULL);
+}
+
+bool p_psDummyFunc(int dimNum,
+                   psVector *params,
+                   const psArray *coords,
+                   psMinimizePowellFunc func)
+{
+    return(true);
+}
+
+
+/******************************************************************************
+    This routine must minimize a possibly multi-dimensional function
+    (several parameters) along a single dimension.
  *****************************************************************************/
 bool psMinimize1DFunc(psMinimization *min,
                       psVector *params,
+                      int dim,
                       const psArray *coords,
                       psMinimizePowellFunc func)
 {
     psVector *bracket;
-    psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
     float a = 0.0;
     float b = 0.0;
     float c = 0.0;
     float n = 0.0;
+    float fa = 0.0;
     float fb = 0.0;
+    float fc = 0.0;
     float fn = 0.0;
-
-    if (params->n != 1) {
-        psAbort(__func__, "params is larger than a single dimension.");
-    }
-    bracket = p_psDetermineBracket(params, coords, func);
-
+    //    float initialParam = params->data.F32[dim];
+
+    bracket = p_psDetermineBracketOld(params, dim, coords, func);
+    if (bracket == NULL) {
+        psAbort(__func__, "Could not bracket minimum.");
+    }
 
     min->iter = 0;
-    while (1) {
+    while (min->iter < min->maxIter) {
         min->iter++;
-
+        //printf("psMinimize1DFunc(): iteration %d\n", min->iter);
         a = bracket->data.F32[0];
         b = bracket->data.F32[1];
         c = bracket->data.F32[2];
-        dummyParams->data.F32[0] = b;
-        fb = func(dummyParams, coords);
+
+        params->data.F32[dim] = a;
+        fa = func(params, coords);
+        params->data.F32[dim] = b;
+        fb = func(params, coords);
+        params->data.F32[dim] = c;
+        fc = func(params, coords);
+        //printf("Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
 
         // We determine which is the biggest segment in [a,b,c] then split
@@ -1075,6 +1306,6 @@
         if ((b-a) > (c-b)) {
             // This is the golden section formula
-            dummyParams->data.F32[0] = n = a + (0.69 * (b-a));
-            fn = func(dummyParams, coords);
+            params->data.F32[dim] = n = a + (0.69 * (b-a));
+            fn = func(params, coords);
             if (fn > fb) {
                 // a = n, b = b, c = c
@@ -1086,6 +1317,6 @@
             }
         } else {
-            dummyParams->data.F32[0] = b + (0.69 * (c-b));
-            fn = func(dummyParams, coords);
+            params->data.F32[dim] = n = b + (0.69 * (c-b));
+            fn = func(params, coords);
             if (fn > fb) {
                 // a = a, b = b, c = n
@@ -1097,9 +1328,183 @@
             }
         }
-    }
-
-    psFree(dummyParams);
-}
-
+
+        if ((fabs(a-b) < min->tol) &&
+                (fabs(b-c) < min->tol)) {
+            //            psFree(bracket);
+            //  XXX: is this line correct?
+            params->data.F32[dim] = bracket->data.F32[1];
+            min->value = func(params, coords);
+            return(true);
+        }
+    }
+
+    //    psFree(bracket);
+    return(false);
+}
+
+/******************************************************************************
+    This routine must minimize a possibly multi-dimensional function
+    along a vector defined by line.
+ *****************************************************************************/
+float psLineMin(psMinimization *min,
+                psVector *params,
+                psVector *line,
+                const psVector *paramMask,
+                const psArray *coords,
+                psMinimizePowellFunc func)
+{
+    psVector *bracket;
+    float a = 0.0;
+    float b = 0.0;
+    float c = 0.0;
+    float n = 0.0;
+    float fa = 0.0;
+    float fb = 0.0;
+    float fc = 0.0;
+    float fn = 0.0;
+    float mul = 0.0;
+    psVector *tmpa = NULL;
+    psVector *tmpb = NULL;
+    psVector *tmpc = NULL;
+    psVector *tmpn = NULL;
+    int i = 0;
+    int null = 0;
+
+    psTrace(".psLib.dataManip.psLineMin", 4,
+            "---- psLineMin() begin ----\n");
+
+    for (i=0;i<params->n;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            if (line->data.F32[i] >= FLT_EPSILON) {
+                null = 1;
+                psTrace(".psLib.dataManip.psLineMin", 4,
+                        "line->data.F32[%d] is %f\n", i, line->data.F32[i]);
+            }
+        }
+    }
+    if (null == 0) {
+        psTrace(".psLib.dataManip.psLineMin", 2,
+                "psLineMin() called with zero line vector.\n");
+        return(0.0);
+    }
+
+    tmpa = psVectorAlloc(params->n, PS_TYPE_F32);
+    tmpb = psVectorAlloc(params->n, PS_TYPE_F32);
+    tmpc = psVectorAlloc(params->n, PS_TYPE_F32);
+    tmpn = psVectorAlloc(params->n, PS_TYPE_F32);
+
+
+    for (i=0;i<params->n;i++) {
+        psTrace(".psLib.dataManip.psLineMin", 6,
+                "params, paramMask, line [%d] is (%f %d %f)\n", i,
+                params->data.F32[i],
+                paramMask->data.U8[i],
+                line->data.F32[i]);
+    }
+
+    bracket = p_psDetermineBracket(params, line, paramMask, coords, func);
+    if (bracket == NULL) {
+        psAbort(__func__, "(1) Could not bracket minimum.");
+    }
+
+    min->iter = 0;
+    while (min->iter < min->maxIter) {
+        min->iter++;
+        psTrace(".psLib.dataManip.psLineMin", 6,
+                "psLineMin(): iteration %d\n", min->iter);
+
+        a = bracket->data.F32[0];
+        b = bracket->data.F32[1];
+        c = bracket->data.F32[2];
+
+        for (i=0;i<params->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                tmpa->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]);
+                tmpb->data.F32[i] = params->data.F32[i] + (b * line->data.F32[i]);
+                tmpc->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]);
+            } else {
+                tmpa->data.F32[i] = params->data.F32[i];
+                tmpb->data.F32[i] = params->data.F32[i];
+                tmpc->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fa = func(tmpa, coords);
+        fb = func(tmpb, coords);
+        fc = func(tmpc, coords);
+        psTrace(".psLib.dataManip.psLineMin", 6,
+                "Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
+
+        // We determine which is the biggest segment in [a,b,c] then split
+        // that with the point n.
+        if ((b-a) > (c-b)) {
+            // This is the golden section formula
+            n = a + (0.69 * (b-a));
+            for (i=0;i<params->n;i++) {
+                tmpn->data.F32[i] = params->data.F32[i] + (n * line->data.F32[i]);
+            }
+            fn = func(tmpn, coords);
+
+            if (fn > fb) {
+                // a = n, b = b, c = c
+                bracket->data.F32[0] = n;
+            } else {
+                // a = a, b = n, c = b
+                bracket->data.F32[1] = n;
+                bracket->data.F32[2] = b;
+            }
+        } else {
+            n = b + (0.69 * (c-b));
+            for (i=0;i<params->n;i++) {
+                tmpn->data.F32[i] = params->data.F32[i] + (n * line->data.F32[i]);
+            }
+            fn = func(tmpn, coords);
+
+            if (fn > fb) {
+                // a = a, b = b, c = n
+                bracket->data.F32[2] = n;
+            } else {
+                // a = b, b = n, c = c
+                bracket->data.F32[0] = b;
+                bracket->data.F32[1] = n;
+            }
+        }
+        psTrace(".psLib.dataManip.psLineMin", 6,
+                "Iteration %d: new bracket is (%f %f %f)\n", min->iter, bracket->data.F32[0], bracket->data.F32[1], bracket->data.F32[2]);
+
+        mul = bracket->data.F32[1];
+        if ((fabs(a-b) < min->tol) &&
+                (fabs(b-c) < min->tol)) {
+            for (i=0;i<params->n;i++) {
+                if (paramMask->data.U8[i] == 0) {
+                    params->data.F32[i]+= mul * line->data.F32[i];
+                }
+            }
+            min->value = func(params, coords);
+            psFree(bracket);
+            psFree(tmpa);
+            psFree(tmpb);
+            psFree(tmpc);
+            psFree(tmpn);
+            psTrace(".psLib.dataManip.psLineMin", 4,
+                    "---- psLineMin() end (%f) ----\n", mul);
+            return(mul);
+        }
+    }
+
+    psFree(bracket);
+    psFree(tmpa);
+    psFree(tmpb);
+    psFree(tmpc);
+    psFree(tmpn);
+
+    psTrace(".psLib.dataManip.psLineMin", 4,
+            "---- psLineMin() end (0.0) ----\n");
+    return(0.0);
+}
+
+
+/******************************************************************************
+    This routine must minimize a possibly multi-dimensional function.
+ *****************************************************************************/
 bool psMinimizePowell(psMinimization *min,
                       psVector *params,
@@ -1108,12 +1513,161 @@
                       psMinimizePowellFunc func)
 {
-
-
-
-
-
-
-
-
-    return(true);
-}
+    int numDims = params->n;
+    int i = 0;
+    int j = 0;
+    psVector **v = NULL;
+    psVector *pQP = NULL;
+    psVector *Q = NULL;
+    psVector *u = NULL;
+    psMinimization dummyMin;
+    float mul = 0.0;
+    float baseFuncVal = 0.0;
+    float currFuncVal = 0.0;
+    int biggestIter = 0;
+    float biggestDiff = 0.0;
+    float term1 = 0.0;
+    float term2 = 0.0;
+    float fqp = 0.0;
+    float diff = 0.0;
+
+    psTrace(".psLib.dataManip.psMinimizePowell", 4,
+            "---- psMinimizePowell() begin ----\n");
+
+    pQP = psVectorAlloc(numDims, PS_TYPE_F32);
+    u = psVectorAlloc(numDims, PS_TYPE_F32);
+    Q = psVectorAlloc(numDims, PS_TYPE_F32);
+
+    // 1: Set v[i] to be the unit vectors for each dimension in params
+    v = (psVector **) psAlloc(params->n * sizeof(psVector **));
+    for (i=0;i<numDims;i++) {
+        v[i] = psVectorAlloc(numDims, PS_TYPE_F32);
+        for (j=0;j<numDims;j++) {
+            if (i == j) {
+                v[i]->data.F32[j] = 1.0;
+            } else {
+                v[i]->data.F32[j] = 0.0;
+            }
+        }
+    }
+
+    // 2: Set Q to be the initial params (P in the ADD)
+    for (i=0;i<numDims;i++) {
+        if (paramMask->data.U8[i] == 0) {
+            Q->data.F32[i] = params->data.F32[i];
+        }
+    }
+
+    min->iter = 0;
+    while (min->iter < min->maxIter) {
+        min->iter++;
+        psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                "psMinimizePowell() iteration %d\n", min->iter);
+
+        // 3: For each dimension in params, move Q only in the vector v[i] to
+        //    minimize the function.
+
+        baseFuncVal = func(Q, coords);
+        currFuncVal = baseFuncVal;
+        diff = 0.0;
+        biggestDiff = 0;
+        biggestIter = 0;
+        for (i=0;i<numDims;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                mul = psLineMin(&dummyMin, Q, v[i], paramMask, coords, func);
+                if (fabs(dummyMin.value - currFuncVal) > biggestDiff) {
+                    biggestDiff = fabs(dummyMin.value - currFuncVal);
+                    biggestIter = i;
+                }
+                currFuncVal = dummyMin.value;
+            }
+        }
+
+        // 4: Set the vector u = Q - P
+        for (i=0;i<numDims;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                u->data.F32[i] = Q->data.F32[i] - params->data.F32[i];
+
+                psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                        "WOW (%f %f %f)\n", u->data.F32[i], Q->data.F32[i], params->data.F32[i]);
+
+            } else {
+                u->data.F32[i] = 0.0;
+            }
+        }
+
+        // 5: Move Q only in the direction u, and minimize the function.
+        psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
+        psTrace(".psLib.dataManip.psMinimizePowell", 6,
+                "u[0] is %f\n", u->data.F32[0]);
+
+        mul = psLineMin(min, params, u, paramMask, coords, func);
+
+        // 6:
+        if (min->value > currFuncVal) {
+            for (i=0;i<numDims;i++) {
+                psFree(v[i]);
+            }
+            psFree(v);
+            psFree(u);
+            psFree(pQP);
+            psFree(Q);
+            psTrace(".psLib.dataManip.psMinimizePowell", 4,
+                    "---- psMinimizePowell() end (true) ----\n");
+            return(true);
+        }
+
+        for (i=0;i<numDims;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                pQP->data.F32[i] = (2 * Q->data.F32[i]) - params->data.F32[i];
+            } else {
+                pQP->data.F32[i] = params->data.F32[i];
+            }
+        }
+        fqp = func(pQP, coords);
+
+        term1 = (baseFuncVal - currFuncVal) - biggestDiff;
+        term1*= term1;
+        term1*= 2.0 * (baseFuncVal - (2.0 * currFuncVal) + fqp);
+        term2 = baseFuncVal - fqp;
+        term2*= term2 * biggestDiff;
+        if (term1 < term2) {
+            for (i=0;i<numDims;i++) {
+                if (paramMask->data.U8[i] == 0) {
+                    v[biggestIter]->data.F32[i] = u->data.F32[i];
+                }
+            }
+        }
+
+        // 7: Set P to Q
+        for (i=0;i<numDims;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                params->data.F32[i] = Q->data.F32[i];
+            }
+        }
+
+        // 8: Go to step 3 until the change is less than some tolerance.
+        if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
+            for (i=0;i<numDims;i++) {
+                psFree(v[i]);
+            }
+            psFree(v);
+            psFree(u);
+            psFree(pQP);
+            psFree(Q);
+            psTrace(".psLib.dataManip.psMinimizePowell", 4,
+                    "---- psMinimizePowell() end (true) ----\n");
+            return(true);
+        }
+    }
+
+    for (i=0;i<numDims;i++) {
+        psFree(v[i]);
+    }
+    psFree(v);
+    psFree(u);
+    psFree(pQP);
+    psFree(Q);
+    psTrace(".psLib.dataManip.psMinimizePowell", 4,
+            "---- psMinimizePowell() end (false) ----\n");
+    return(false);
+}