Changeset 1122 for trunk/psLib/src/dataManip/psMinimize.c

Timestamp:

Jun 28, 2004, 7:37:45 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/dataManip/psMinimize.c (modified) (14 diffs)

trunk/psLib/src/dataManip/psMinimize.c

@@ -5 +5 @@
 #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"
@@ -23 +18 @@
 #include <math.h>
 
+#include <gsl/gsl_multifit_nlin.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
 #define MAX_LMM_NUM_ITERATIONS 500
 typedef struct
 {
     size_t n;   // Number of data points points in domain.
-    int count;   // Number of non-masked parameters.
+    int paramCount;   // Number of non-masked parameters.
     psVector *restrict       initialGuess;
     const psImage *restrict  domain;
@@ -33 +34 @@
     const psVector *restrict  errors;
     const psVector *restrict paramMask;
-    float (*evalModel) (psVector *, psVector *);
-    float (*d_evalModel) (psVector *, psVector *, int);
+    float (*evalModel) (const psVector *, const psVector *);
+    float (*d_evalModel) (const psVector *, const psVector *, int);
 }
 psModelData;
@@ -66 +67 @@
 only call 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:
  
- 
+    x: These are the parameters which are to be varied by GSL in order to
+ minimized chi2 over the data set.
+    params: this data structure contains the input values over which the
+ function will be evaluated, the expected value of the function at
+ those points, the amount of error tolerable at those points, a mask
+ vector which specifies which parameters to the function are to be
+ constant, and an initial guess at the parameters.
+    outData: The function is evaluated at each point, then subtract the
+ expected value and divide by the error.
  ******************************************************************************
  *****************************************************************************/
 int gsl_function_f(const gsl_vector *x,
                    void *params,
-                   gsl_vector *f)
-{
-    return 0;
+                   gsl_vector *outData)
+{
     int i;    // Loop index variable.
     int j;    // Loop index variable.
@@ -87 +91 @@
     const psVector *restrict mask     = ((psModelData *) params)->paramMask;
     psVector *restrict initialGuess = ((psModelData *)params)->initialGuess;
-    float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel;
+    float (*evalModel)(const psVector *, const psVector *) = ((psModelData *) params)->evalModel;
     psVector *inputParameterList = NULL;
     psVector *tmpVecPtr = NULL;
-
 
     tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
@@ -114 +117 @@
         }
     }
+    for (i=0;i<inputParameterList->n;i++) {
+        printf("(gsl_f(): inputParameterList->data.F32[%d] is %.1f\n", i, inputParameterList->data.F32[i]);
+    }
 
     for (i=0;i<domain->numRows;i++) {
-        for (j=0;j<tmpVecPtr->n;j++) {
+        //        printf("Data item %d is ( ", i);
+        for (j=0;j<domain->numCols;j++) {
             tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
-        }
-        tmpf = evalModel(inputParameterList, tmpVecPtr);
-        gsl_vector_set(f, i, (tmpf - data->data.F32[i])/errors->data.F32[i]);
+            //            printf("%.1f ", tmpVecPtr->data.F32[j]);
+        }
+        tmpf = evalModel(tmpVecPtr, inputParameterList);
+        //        printf(" ).  Output is %.1f\n", tmpf);
+
+        gsl_vector_set(outData, i, (tmpf - data->data.F32[i])/
+                       errors->data.F32[i]);
+        //printf("--------- MYFUNC((%.1f) %.1f %.1f %.1f) is [%.1f] ---------\n",
+        //tmpVecPtr->data.F32[0],
+        //inputParameterList->data.F32[0],
+        //inputParameterList->data.F32[1],
+        //inputParameterList->data.F32[2],
+        //(tmpf - data->data.F32[i])/errors->data.F32[i]);
+
+    }
+
+    for (i=0;i<domain->numRows;i++) {
+        printf("gsl_vector_set(outData, %d, %.1f)\n", i,
+               gsl_vector_get(outData, i));
     }
 
@@ -138 +161 @@
     psVector *inputParameterList = NULL;
     psVector *tmpVecPtr = NULL;
-    float (*d_evalModel)(psVector *, psVector *, int) = ((psModelData *) params)->d_evalModel;
+    float (*d_evalModel)(const psVector *, const psVector *, int) = ((psModelData *) params)->d_evalModel;
 
     size_t i;
@@ -150 +173 @@
     // have those parameters removed.  Here will create a new parameter list
     // with the masked parameters added (we expand initialGuess).
-
-    psMemCheckCorruption(1);
-    printf("Calling psVectorAlloc(%d)\n", initialGuess->n);
-
     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
-    printf("Called psVectorAlloc(%d)\n", initialGuess->n);
 
     if (mask != NULL) {
@@ -178 +196 @@
 
         for (j=0;j<inputParameterList->n;j++) {
-            tmpf = d_evalModel(inputParameterList, tmpVecPtr, j);
+            tmpf = d_evalModel(tmpVecPtr, inputParameterList, j);
+
+            //printf("--------- MYFUNCDERIV((%.1f) %.1f %.1f %.1f, %d) is [%.1f] ---------\n",
+            //tmpVecPtr->data.F32[0],
+            //inputParameterList->data.F32[0],
+            //inputParameterList->data.F32[1],
+            //inputParameterList->data.F32[2],
+            //j,
+            //tmpf/errors->data.F32[i]);
+
             gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i]));
         }
     }
+
+
+    for (i=0;i<domain->numRows;i++) {
+        for (j=0;j<inputParameterList->n;j++) {
+            printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j,
+                   gsl_matrix_get(J, i, j));
+        }
+    }
+
+
 
     psFree(inputParameterList);
@@ -193 +230 @@
                         gsl_matrix *J)
 {
-    return 0;
     gsl_function_f(x, params, f);
     gsl_function_df(x, params, J);
@@ -215 +251 @@
                float *chiSq)
 {
-
-    int numData = domain->numCols; // Number of data points
+    printf("Calling psMinimizeChi2()\n");
+
+    int numData = domain->numRows; // Number of data points
     int status;    // Return status for the GSL solver.
     int i = 0;    // Loop index variable.
@@ -225 +262 @@
     double *xInit = NULL;     // The initial guess at the parameters
     // with masked parameters removed.
-    const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
+    const gsl_multifit_fdfsolver_type *T;
     // This tells GSL to use the Levenberg-
     // Marquardt algorithm for chi2
     // minimization.
     gsl_multifit_fdfsolver *s; // GSL data structure.
-
     psModelData inputData;
+
     inputData.n = numData;
-    inputData.count = 0;
+    inputData.paramCount = 0;
     inputData.initialGuess = initialGuess;
     inputData.domain = domain;
-    inputData.count = 0;
+    inputData.data = data;
+    inputData.errors = errors;
+    inputData.paramMask = paramMask;
+    inputData.evalModel = evalModel;
+    inputData.d_evalModel = DevalModel;
+
     if (paramMask != NULL) {
         for (i=0;i<paramMask->n;i++) {
             if (paramMask->data.U8[i] != 0) {
-                inputData.count++;
+                inputData.paramCount++;
             }
         }
     } else {
-        inputData.count= initialGuess->n;
+        inputData.paramCount= initialGuess->n;
     }
 
@@ -251 +293 @@
     // parameters are masked out.
 
-    xInit = (double *) psAlloc(inputData.count * sizeof(double));
+    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++] = (double) initialGuess->data.F32[i];
+                printf("xInit[%d] is %f (masked loop)\n", j, initialGuess->data.F32[i]);
+                xInit[j++] = initialGuess->data.F32[i];
             }
         }
     } else {
         for (i=0;i<initialGuess->n;i++) {
-            xInit[i] = (double) initialGuess->data.F32[i];
+            xInit[i] = initialGuess->data.F32[i];
+            printf("xInit[%d] is %f\n", i, initialGuess->data.F32[i]);
         }
     }
@@ -281 +325 @@
     f.fdf = &gsl_my_function_fdf;
     f.n = numData;
-    f.p = inputData.count;
+    f.p = inputData.paramCount;
     f.params = &inputData;
 
-    printf("inputData.count is %d\n", inputData.count);
+    printf("inputData.paramCount is %d\n", inputData.paramCount);
     printf("numData is %d\n", numData);
     // Creates the vector for x which GSL uses.  Must deallocate.
-    gsl_vector_view x = gsl_vector_view_array(xInit, inputData.count);
+    gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount);
+    T = gsl_multifit_fdfsolver_lmsder;
 
     // Create an instance of the GSL solver that we will be iterating on.
-    // It will have numData data points and inputData.count parameters.
-    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.count);
-
-    printf("HERE 04 (HMMM)\n");
+    // It will have numData data points and inputData.paramCount parameters.
+    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount);
+
     // Initialize the GSL minimizer to use function defined by the data
     // structure "f" and x.vector as an initial guess for the parameters.
-    psMemCheckCorruption(1);
-    printf("HERE 05 (HMMM)\n");
+
     gsl_multifit_fdfsolver_set(s, &f, &x.vector);
+
     // Each iteration of the following loop will perform one step in an
     // attempt to minimized chi-squared for the function.  The loop exits
     // either when the change in parameters is small enough, or when the
     // maximum number of iterations is reached.
-    printf("HERE 06.0\n");
     do {
         iter++;
+        printf("##################################################\n");
+        printf("##################################################\n");
+        printf("##################################################\n");
+        printf("################## Iteration %d ##################\n", iter);
+        printf("##################################################\n");
+        printf("##################################################\n");
+        printf("##################################################\n");
+
         // Perform an iteration of the GSL solver.
-        printf("HERE 06\n");
+
         status = gsl_multifit_fdfsolver_iterate(s);
-        printf("HERE 07\n");
         // If there was a problem, abort.
         if (status) {
-            exit(1);
-            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
-        }
-        printf("HERE 08\n");
+            printf("ABORT: status is %d\n", status);
+            //            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
+        }
 
         // Test if the parameters changed by a small enough amount.
         status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
-        printf("HERE 09\n");
     } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS);
-    printf("HERE 12\n");