Index: /trunk/psLib/src/dataManip/psMinimize.c
===================================================================
--- /trunk/psLib/src/dataManip/psMinimize.c	(revision 1078)
+++ /trunk/psLib/src/dataManip/psMinimize.c	(revision 1079)
@@ -5,4 +5,9 @@
 #include <float.h>
 #include <math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multifit_nlin.h>
 
 #include "psMemory.h"
@@ -19,4 +24,20 @@
 #define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
 int MyInfoLevel = 0;
+
+typedef struct
+{
+    size_t n;
+    psVector *initialGuess;
+    psVector *data;
+    psVector *errors;
+    psVector *paramMask;
+    float (*evalModel) (psVector *, psVector *);
+    float (*d_evalModel) (psVector *, psVector *);
+}
+psmodeldata;
+
+
+
+
 /******************************************************************************
     This routine must minimize an arbitrary function.
@@ -30,10 +51,115 @@
 }
 
+
+
+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;
+    size_t i;
+
+    for (i=0;i<initialguess->n;i++) {
+        if (mask->data.U8[i] != 0) {
+            count++;
+        }
+    }
+
+    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];
+            }
+        }
+    } 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]);
+    }
+    return GSL_SUCCESS;
+}
+
+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;
+    size_t i;
+
+    for (i=0;i<initialguess->n;i++) {
+        if (mask->data.U8[i] != 0) {
+            count++;
+        }
+    }
+
+
+    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];
+            }
+        }
+    } 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]);
+        }
+    }
+    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);
+
+    return GSL_SUCCESS;
+}
+
+
+
 /******************************************************************************
     This routine must minimize an arbitrary function.
  *****************************************************************************/
 psVector *
-psMinimizeChi2(float (*evalModel)(const psVector *restrict,
-                                  const psVector *restrict),
+psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict),
+               float (*DevalModel)(psVector *, psVector *),
                const psVector *restrict domain,
                const psVector *restrict data,
@@ -42,6 +168,114 @@
                const psVector *restrict paramMask)
 {
-    return(NULL);
-}
+    size_t n = size(domain);
+    size_t p = size(initialGuess);
+    gsl_multifit_fdfsolver_type *T = NULL;
+    gsl_multifit_fdfsolver *s = NULL;
+    int status;
+    size_t i = 0;
+    size_t iter = 0;
+    psmodeldata inputdata = { n, initialguess, data, errors, paramMask, evalmodel, d_evalmodel};
+    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;
+    gsl_rng *r;
+    gsl_rng_env_setup();
+
+    type = gsl_rng_default;
+    r = gsl_rng_alloc (type);
+
+    f.f = &gsl_function_f;
+    f.df = &gsl_function_df;
+    f.fdf = &gsl_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);
+
+    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;
+
+        // 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);
+    } 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));
+}
+
+
+
+
+
+
 
 /** @brief This procedure calculates various combinations of powers of x and y
Index: /trunk/psLib/src/image/psImageStats.h
===================================================================
--- /trunk/psLib/src/image/psImageStats.h	(revision 1078)
+++ /trunk/psLib/src/image/psImageStats.h	(revision 1079)
@@ -10,6 +10,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-14 19:40:14 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-24 02:39:34 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -48,3 +48,6 @@
                      );
 
+float psImagePixelInterpolation(psImage *input,
+                                float x,
+                                float y);
 #endif
Index: /trunk/psLib/src/imageops/psImageStats.h
===================================================================
--- /trunk/psLib/src/imageops/psImageStats.h	(revision 1078)
+++ /trunk/psLib/src/imageops/psImageStats.h	(revision 1079)
@@ -10,6 +10,6 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-14 19:40:14 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-24 02:39:34 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -48,3 +48,6 @@
                      );
 
+float psImagePixelInterpolation(psImage *input,
+                                float x,
+                                float y);
 #endif
Index: /trunk/psLib/src/math/psMinimize.c
===================================================================
--- /trunk/psLib/src/math/psMinimize.c	(revision 1078)
+++ /trunk/psLib/src/math/psMinimize.c	(revision 1079)
@@ -5,4 +5,9 @@
 #include <float.h>
 #include <math.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multifit_nlin.h>
 
 #include "psMemory.h"
@@ -19,4 +24,20 @@
 #define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
 int MyInfoLevel = 0;
+
+typedef struct
+{
+    size_t n;
+    psVector *initialGuess;
+    psVector *data;
+    psVector *errors;
+    psVector *paramMask;
+    float (*evalModel) (psVector *, psVector *);
+    float (*d_evalModel) (psVector *, psVector *);
+}
+psmodeldata;
+
+
+
+
 /******************************************************************************
     This routine must minimize an arbitrary function.
@@ -30,10 +51,115 @@
 }
 
+
+
+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;
+    size_t i;
+
+    for (i=0;i<initialguess->n;i++) {
+        if (mask->data.U8[i] != 0) {
+            count++;
+        }
+    }
+
+    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];
+            }
+        }
+    } 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]);
+    }
+    return GSL_SUCCESS;
+}
+
+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;
+    size_t i;
+
+    for (i=0;i<initialguess->n;i++) {
+        if (mask->data.U8[i] != 0) {
+            count++;
+        }
+    }
+
+
+    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];
+            }
+        }
+    } 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]);
+        }
+    }
+    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);
+
+    return GSL_SUCCESS;
+}
+
+
+
 /******************************************************************************
     This routine must minimize an arbitrary function.
  *****************************************************************************/
 psVector *
-psMinimizeChi2(float (*evalModel)(const psVector *restrict,
-                                  const psVector *restrict),
+psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict),
+               float (*DevalModel)(psVector *, psVector *),
                const psVector *restrict domain,
                const psVector *restrict data,
@@ -42,6 +168,114 @@
                const psVector *restrict paramMask)
 {
-    return(NULL);
-}
+    size_t n = size(domain);
+    size_t p = size(initialGuess);
+    gsl_multifit_fdfsolver_type *T = NULL;
+    gsl_multifit_fdfsolver *s = NULL;
+    int status;
+    size_t i = 0;
+    size_t iter = 0;
+    psmodeldata inputdata = { n, initialguess, data, errors, paramMask, evalmodel, d_evalmodel};
+    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;
+    gsl_rng *r;
+    gsl_rng_env_setup();
+
+    type = gsl_rng_default;
+    r = gsl_rng_alloc (type);
+
+    f.f = &gsl_function_f;
+    f.df = &gsl_function_df;
+    f.fdf = &gsl_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);
+
+    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;
+
+        // 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);
+    } 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));
+}
+
+
+
+
+
+
 
 /** @brief This procedure calculates various combinations of powers of x and y