Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 1079)
+++ trunk/psLib/src/math/psMinimize.c	(revision 1086)
@@ -13,4 +13,5 @@
 #include "psMemory.h"
 #include "psVector.h"
+#include "psImage.h"
 #include "psTrace.h"
 #include "psError.h"
@@ -27,15 +28,21 @@
 typedef struct
 {
-    size_t n;
-    psVector *initialGuess;
+    size_t n;   // Number of data points points in domain.
+    int count;   // Number of non-masked parameters.
+    psVector *restrict initialGuess;
+    psImage  *domain;  //
     psVector *data;
     psVector *errors;
     psVector *paramMask;
+    // 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.
     float (*evalModel) (psVector *, psVector *);
     float (*d_evalModel) (psVector *, psVector *);
 }
-psmodeldata;
-
-
+psModelData;
 
 
@@ -53,101 +60,118 @@
 
 
+/******************************************************************************
+ ******************************************************************************
+gsl_function_f(x, params, f): This function serves as a standardized wrapper
+for the user supplied function which is to be minimized.  The GSL
+minimization routines have no knowledge of the user-supplied function.  They
+only clal this routine, which then calls the user-supplied function.  The
+arguments are:
+    x: These are the parameter which are to be varied by GSL in order to
+       minimized chi2 over the data set.
+    params:
+    f
+ 
+ 
+ ******************************************************************************
+ *****************************************************************************/
 int gsl_function_f(const gsl_vector *x,
                    void *params,
                    gsl_vector *f)
 {
-    size_t n  = ((psmodeldata *)params)->n;
-    psVector *data   = ((psmodeldata *)params)->data;
-    psVector *errors = ((psmodeldata *) params)->errors;
-    psVector *mask   = ((psmodeldata *) params)->paramMask;
-    psVector *initialguess = ((psmodeldata *)params)->initialguess;
-    int count = 0;
-    float (*evalModel) = ((psmodeldata *) params)->evalModel;
+    psImage  *domain   = ((psModelData *)params)->domain;
+    psVector *data     = ((psModelData *)params)->data;
+    psVector *errors   = ((psModelData *) params)->errors;
+    psVector *mask     = ((psModelData *) params)->paramMask;
+    psVector *initialGuess = ((psModelData *)params)->initialGuess;
+    psVector *inputParameterList = NULL;
+    //    float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel;
     size_t i;
-
-    for (i=0;i<initialguess->n;i++) {
-        if (mask->data.U8[i] != 0) {
-            count++;
-        }
-    }
+    int j;
+    float tmpf;
+
+    // The GSL routines will call this functions with the masked parameters
+    // removed.  However, the user-supplied function (to be modified) does not
+    // have those parameters removed.  Here will create a new parameter list
+    // with the masked parameters added (we expand initialGuess).
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
 
     if (mask != NULL) {
         j = 0;
-        inputparameterlist = psVectorAlloc(count, PS_TYPE_F32);
-        for (i=0;i<initialguess->n;i++) {
-            if (mask->data.U8[i] == 0) {
-                inputparameterlist->data.F32[j++] = initialguess->data.F32[i];
+        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(x, j++);
             }
         }
     } else {
-        for (i=0;i<initialguess->n;i++) {
-            inputparameterlist->data.F32[i] = initialguess->data.F32[i];
-        }
-    }
-    // PAUL: what is this?
-    evalModel(entiredomain, inputparameterlist);
-    ???
-
-    loop i over domain size n
-    {
-        gsl_vector_set (f, i, (modelvalue[i] - data[i])/errors[i]);
-    }
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+        }
+    }
+
+    for (i=0;i<domain->numRows;i++) {
+        //GUS        tmpf = evalModel(inputParameterList, domain->data.F32[i]);
+        gsl_vector_set(f, i, (tmpf - data->data.F32[i])/errors->data.F32[i]);
+    }
+
     return GSL_SUCCESS;
 }
 
-int gsl_function_df(const gsl_vector * x,
+int gsl_function_df(const gsl_vector *x,
                     void *params,
-                    gsl_matrix* J)
-{
-    size_t n  = ((psmodeldata *)params)->n;
-    psVector *data   = ((psmodeldata *)params)->data;
-    psVector *errors = ((psmodeldata *) params)->errors;
-    psVector *mask   = ((psmodeldata *) params)->paramMask;
-    psVector *initialguess = ((psmodeldata *)params)->initialguess;
-    int count = 0;
-    float (*d_evalModel) = ((psmodeldata *) params)->evalModel;
+                    gsl_matrix *J)
+{
+    //    size_t n  = ((psModelData *)params)->n;
+    psImage  *domain   = ((psModelData *)params)->domain;
+    //    psVector *data   = ((psModelData *)params)->data;
+    psVector *errors = ((psModelData *) params)->errors;
+    psVector *mask   = ((psModelData *) params)->paramMask;
+    psVector *initialGuess = ((psModelData *)params)->initialGuess;
+    psVector *inputParameterList = NULL;
+    //    int count = 0;
+    //    float (*d_evalModel) = ((psModelData *) params)->evalModel;
     size_t i;
-
-    for (i=0;i<initialguess->n;i++) {
-        if (mask->data.U8[i] != 0) {
-            count++;
-        }
-    }
-
-
+    int j;
+    float tmpf;
+
+    // The GSL routines will call this functions with the masked parameters
+    // removed.  However, the user-supplied function (to be modified) does not
+    // have those parameters removed.  Here will create a new parameter list
+    // with the masked parameters added (we expand initialGuess).
+    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
     if (mask != NULL) {
         j = 0;
-        inputparameterlist = psVectorAlloc(count, PS_TYPE_F32);
-        for (i=0;i<initialguess->n;i++) {
-            if (mask->data.U8[i] == 0) {
-                inputparameterlist->data.F32[j++] = initialguess->data.F32[i];
+        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(x, j++);
             }
         }
     } else {
-        for (i=0;i<initialguess->n;i++) {
-            inputparameterlist->data.F32[i] = initialguess->data.F32[i];
-        }
-    }
-
-    d_evalModel(entiredomain, inputparameterlist);
-    ???
-    loop i over domain size n
-    // Jacobian matrix J(i,j) = dfi / dxj, where fi = (Modeli - datai)/errors[i]
-    {
-        loop j over the number of coefficients in model (excluding masks)
-        {
-            gsl_matrix_set (J, i, j, (d_modelvalue[i][j])/errors[i]);
-        }
-    }
+        for (i=0;i<initialGuess->n;i++) {
+            inputParameterList->data.F32[i] = initialGuess->data.F32[i];
+        }
+    }
+
+    for (i=0;i<domain->numRows;i++) {
+        for (j=0;j<inputParameterList->n;j++) {
+            //GUS            tmpf = d_evalModel(inputParameterList, domain->data.F32[i], j);
+
+            gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i]));
+        }
+    }
+
     return GSL_SUCCESS;
 }
 
-int gsl_function_fdf(const gsl_vector * x,
-                     void *params,
-                     gsl_vector * f,
-                     gsl_matrix * J)
-{
-    gsl_function_f (x, params, f);
-    gsl_function_df (x, params, J);
+int gsl_my_function_fdf(const gsl_vector *x,
+                        void *params,
+                        gsl_vector *f,
+                        gsl_matrix *J)
+{
+    gsl_function_f(x, params, f);
+    gsl_function_df(x, params, J);
 
     return GSL_SUCCESS;
@@ -166,116 +190,136 @@
                const psVector *restrict errors,
                psVector *restrict initialGuess,
-               const psVector *restrict paramMask)
-{
-    size_t n = size(domain);
-    size_t p = size(initialGuess);
-    gsl_multifit_fdfsolver_type *T = NULL;
-    gsl_multifit_fdfsolver *s = NULL;
+               const psVector *restrict paramMask,
+               float *chiSq)
+{
+    size_t n = domain->n;
+    const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
+    gsl_multifit_fdfsolver *s      = NULL;
     int status;
     size_t i = 0;
+    int j = 0;
     size_t iter = 0;
-    psmodeldata inputdata = { n, initialguess, data, errors, paramMask, evalmodel, d_evalmodel};
+    psModelData *inputData = NULL;
     gsl_multifit_function_fdf f;
-    double *x_init = NULL;
-
-
-    //GUS: call psStatsGetnValues, allocate x_init.
-    /* compress initialguess in light of mask */
-    j = 0;
-    loop i over length of initialguess
-    {
-        if not masked[i]
-        {
-            x_init[j] = initialguess[i];
-            j++;
-        } else
-        { /* reduce p accordingly */
-            p--;
-        }
-}
-
-
-// Creates the vector for x which GSL uses.  Must deallocate.
-gsl_vector_view x = gsl_vector_view_array (x_init, p);
-
-    const gsl_rng_type * type;
+    double *xInit = NULL;
+
+    inputData = (psModelData *) psAlloc(sizeof(psModelData));
+    inputData->n = n;
+    inputData->count = 0;
+    inputData->initialGuess = initialGuess;
+    //    inputData->domain = domain;
+    //    inputData->data = data;
+    //    inputData->errors = errors;
+    //    inputData->paramMask = paramMask;
+    //    inputData->evalModel = evalModel;
+    //    inputData->d_evalModel = DevalModel;
+
+    inputData->count = 0;
+    if (paramMask != NULL) {
+        for (i=0;i<paramMask->n;i++) {
+            if (paramMask->data.U8[i] != 0) {
+                inputData->count++;
+            }
+        }
+    } else {
+        inputData->count= 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.  How can this possibly work?  The user-
+    // supplied function will require a fixed number of parameters.
+
+    xInit = (double *) psAlloc(inputData->count * 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];
+        }
+    }
+
+    // Creates the vector for x which GSL uses.  Must deallocate.
+    gsl_vector_view x = gsl_vector_view_array(xInit, inputData->count);
+
+    const gsl_rng_type *type;
     gsl_rng *r;
     gsl_rng_env_setup();
 
     type = gsl_rng_default;
-    r = gsl_rng_alloc (type);
+    r = gsl_rng_alloc(type);
 
     f.f = &gsl_function_f;
     f.df = &gsl_function_df;
-    f.fdf = &gsl_function_fdf;
+    f.fdf = &gsl_my_function_fdf;
     f.n = n;
-    f.p = p;
-    f.params = &inputdata;
-
-    T = gsl_multifit_fdfsolver_lmsder;
-    s = gsl_multifit_fdfsolver_alloc (T, n, p);
-    gsl_multifit_fdfsolver_set (s, &f, &x.vector);
-
-    print_state (iter, s);
-
+    f.p = inputData->count;
+    f.params = &inputData;
+
+    // This tells GSL to use the Levenberg-Marquardt algorithm for chi2
+    // minimization.
+    //    T = gsl_multifit_fdfsolver_lmsder;
+
+    // Creates an instance of the GSL solver that we will be iterating on.
+    s = gsl_multifit_fdfsolver_alloc(T, n, inputData->count);
+
+    // Initialize the GSL minimizer.
+    gsl_multifit_fdfsolver_set(s, &f, &x.vector);
+
+    //    print_state(iter, s);
     do {
         iter++;
-        status = gsl_multifit_fdfsolver_iterate (s);
-
-        printf ("status = %s\n", gsl_strerror (status));
-
-        print_state (iter, s);
-
-        // Figure this out
-        if (status)
-            break;
+        // Perform an iteration of the GSL solver.
+        status = gsl_multifit_fdfsolver_iterate(s);
+
+        //        printf("status = %s\n", gsl_strerror(status));
+        //        print_state(iter, s);
+
+        // If there was a problem, abort.
+        if (status) {
+            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
+        }
 
         // Checks whether the (L2-norm) computed derivative and the difference
-        // between the real/actual for that test x-vector.  If were close enough
-        // exit loop.
-        status = gsl_multifit_test_delta (s->dx, s->x,
-                                          1e-4, 1e-4);
+        // between the real/actual for that test x-vector.  If were close
+        // enough exit loop.
+
+        // Test if the parameters changed by a small enough amount.
+        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
     } while (status == GSL_CONTINUE && iter < 500);
-    gsl_matrix *covar = gsl_matrix_alloc (p, p);
-    gsl_multifit_covar (s->J, 0.0, covar);
-    gsl_matrix_fprintf (stdout, covar, "%g");
-
-    /* add in masked values */
-    j = 0;
-    loop i over the length of initialguess
-    {
-        if mask(i)
-    {
-        // do nothing, initialguess shouldn't have changed */
-    } else
-    {
-        initialguess(i)  = gsl_vector_get(s->x, j);  /* output */
-            uncertainity = sqrt(gsl_matrix_get(covar,j,j));
-            j++;
-        }
-    }
-    double chi = gsl_blas_dnrm2(s->f);
-    gsl_multifit_fdfsolver_free (s);
-    return 0;
-}
-
-// For 3-element vectors
-int print_state(size_t iter,
-                gsl_multifit_fdfsolver * s)
-{
-    printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
-            "|f(x)| = %g\n",
-            iter,
-            gsl_vector_get (s->x, 0),
-            gsl_vector_get (s->x, 1),
-            gsl_vector_get (s->x, 2),
-            gsl_blas_dnrm2 (s->f));
-}
-
-
-
-
-
-
+
+
+
+
+    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);
+        }
+    }
+
+    *chiSq = gsl_blas_dnrm2(s->f);
+    gsl_multifit_fdfsolver_free(s);
+    // Free all allocated memory
+    psFree(xInit);
+    psFree(inputData);
+
+    return(initialGuess);
+}
 
 /** @brief This procedure calculates various combinations of powers of x and y