Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 1176)
+++ trunk/psLib/src/math/psMinimize.c	(revision 1185)
@@ -7,4 +7,5 @@
 
 #include <gsl/gsl_multifit_nlin.h>
+#include <gsl/gsl_multimin.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
@@ -38,5 +39,4 @@
 }
 psModelData;
-
 // The first argument to evalModel() and
 // d_evalModel() specifies the data point.
@@ -47,17 +47,241 @@
 
 
+typedef struct
+{
+    int paramCount;   // Number of non-masked parameters.
+    psVector *restrict       initialGuess;
+    const psVector *restrict  coord;
+    const psVector *restrict paramMask;
+    float (*evalModel) (const psVector *, const psVector *);
+    float (*d_evalModel) (const psVector *, const psVector *, int);
+}
+psModelData2;
+
 
 /******************************************************************************
-    This routine must minimize an arbitrary function.
+p_psMinFunc(*v, *params): This routine calls the user supplied function to
+be minimized.  The "v" argument of this function corresponds to the
+parameters of the function to be minimized.
+ 
+ 
+ *****************************************************************************/
+double p_psMinFunc(const gsl_vector *v,
+                   void *params)
+{
+    int i;    // Loop index variable.
+    int j;    // Loop index variable.
+    float tmpf;    // Temporary floating point variable.
+    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
+    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
+    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
+    float (*evalModel)(const psVector *, const psVector *) =
+        ((psModelData2 *) params)->evalModel;
+    psVector *inputParameterList = NULL;
+
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
+
+    if (mask != NULL) {
+        j = 0;
+        for (i=0;i<mask->n;i++) {
+            if (mask->data.U8[i] != 0) {
+                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+            } else {
+                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
+        }
+    }
+    //printf("HMMM: (%.1f %.1f %.1f %.1f)\n",
+    //       inputParameterList->data.F32[0],
+    //       inputParameterList->data.F32[1],
+    //       coord->data.F32[0],
+    //       coord->data.F32[1]);
+
+    tmpf = evalModel(inputParameterList, coord);
+    psFree(inputParameterList);
+    //printf("Called p_psMinFunc()\n");
+    return(tmpf);
+}
+
+void p_psMinFuncDeriv(const gsl_vector *v,
+                      void *params,
+                      gsl_vector *df)
+{
+    int i;    // Loop index variable.
+    int j;    // Loop index variable.
+    float tmpf;    // Temporary floating point variable.
+    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
+    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
+    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
+    float (*d_evalModel)(const psVector *, const psVector *, int) =
+        ((psModelData2 *) params)->d_evalModel;
+    psVector *inputParameterList = NULL;
+
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
+
+    if (mask != NULL) {
+        j = 0;
+        for (i=0;i<mask->n;i++) {
+            if (mask->data.U8[i] != 0) {
+                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+            } else {
+                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
+        }
+    }
+
+    for (i=0;i<initialGuess->n;i++) {
+        tmpf = d_evalModel(inputParameterList, coord, i);
+        gsl_vector_set(df, i, tmpf);
+    }
+    psFree(inputParameterList);
+}
+
+/* Compute both f and df together. */
+void p_psMinFuncFuncDeriv(const gsl_vector *x,
+                          void *params,
+                          double *f,
+                          gsl_vector *df)
+{
+    *f = p_psMinFunc(x, params);
+    p_psMinFuncDeriv(x, params, df);
+}
+
+
+/******************************************************************************
+psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
+ 
+This routine must minimize an arbitrary function; it must determine the set
+of parameters of that function such that the 
+ 
  *****************************************************************************/
 psVector *
-psMinimize(float (*myFunction)(const psVector *restrict),
-           psVector *restrict initialGuess,
-           psVector *restrict paramMask)
-{
-    return(NULL);
-}
-
-
+psMinimize(psVector *restrict initialGuess,
+           float (*myFunction)(const psVector *restrict, const psVector *restrict),
+           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
+           const psVector *restrict coord,
+           const psVector *restrict paramMask)
+{
+    int status;
+    int i = 0;
+    int j = 0;
+    int iter = 0;
+    gsl_multimin_function_fdf f;
+    double *xInit = NULL;
+    const gsl_multimin_fdfminimizer_type *T;
+    gsl_multimin_fdfminimizer *s;
+    psModelData2 inputData;
+    gsl_vector *x;
+
+    inputData.initialGuess = initialGuess;
+    inputData.coord = coord;
+    inputData.paramMask = paramMask;
+    inputData.evalModel = myFunction;
+    inputData.d_evalModel = myFunctionDeriv;
+    inputData.paramCount = 0;
+
+    printf("psMinimize(): initialGuess->n is %d\n", initialGuess->n);
+    printf("psMinimize(): coord->n is %d\n", coord->n);
+    printf("psMinimize(): coord is (%.1f, %.1f)\n", coord->data.F32[0],
+           coord->data.F32[1]);
+
+    if (paramMask != NULL) {
+        for (i=0;i<paramMask->n;i++) {
+            if (paramMask->data.U8[i] != 0) {
+                inputData.paramCount++;
+            }
+        }
+    } else {
+        inputData.paramCount= initialGuess->n;
+    }
+
+    // The initial guess at the parameters for the function are written into
+    // the vector inputParameterList.  If the paramMask is not NULL, then those
+    // parameters are masked out.
+
+    xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
+    if (paramMask != NULL) {
+        j = 0;
+        for (i=0;i<initialGuess->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                xInit[j++] = initialGuess->data.F32[i];
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            xInit[i] = initialGuess->data.F32[i];
+        }
+    }
+    f.f = &p_psMinFunc;
+    f.df = &p_psMinFuncDeriv;
+    f.fdf = &p_psMinFuncFuncDeriv;
+    f.n = inputData.paramCount;
+    // GUS: is this correct?
+    //f.params = xInit;
+    f.params = &inputData;
+
+
+    printf("inputData.paramCount is %d\n", inputData.paramCount);
+    x = gsl_vector_alloc(inputData.paramCount);
+    for (i=0;i<inputData.paramCount;i++) {
+        gsl_vector_set(x, i, xInit[i]);
+    }
+    T = gsl_multimin_fdfminimizer_conjugate_fr;
+    s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
+
+    gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
+
+    do {
+        iter++;
+        status = gsl_multimin_fdfminimizer_iterate(s);
+
+        if (status)
+            break;
+
+        status = gsl_multimin_test_gradient(s->gradient, 1e-3);
+
+        if (status == GSL_SUCCESS)
+            printf ("Minimum found at:\n");
+
+    } while (status == GSL_CONTINUE && iter < 100);
+
+    // In the above steps we had removed the masked elements from the
+    // the solver.  This next code blocks puts those masked elements
+    // into the solution.
+    if (paramMask != NULL) {
+        j = 0;
+        for (i=0;i<initialGuess->n;i++) {
+            if (paramMask->data.U8[i] == 0) {
+                initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
+            } else {
+                initialGuess->data.F32[i] = initialGuess->data.F32[i];
+            }
+        }
+    } else {
+        for (i=0;i<initialGuess->n;i++) {
+            initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
+        }
+    }
+    psFree(xInit);
+    return(initialGuess);
+}
+
+
+
+
+
+// The first argument to evalModel() and
+// d_evalModel() specifies the data point.
+// It must have the same size as the second
+// dimension of *domain.
+// The second argument must have the same size
+// as *initialGuess and *paramMask.
 
 /******************************************************************************