Changeset 1128 for trunk/psLib/src/math/psMinimize.c

Timestamp:

Jun 29, 2004, 1:21:08 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/math/psMinimize.c (modified) (15 diffs)

trunk/psLib/src/math/psMinimize.c

@@ -5 +5 @@
 #include <float.h>
 #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>
 
 #include "psMemory.h"
@@ -18 +24 @@
 #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
@@ -45 +45 @@
 // The second argument must have the same size
 // as *initialGuess and *paramMask.
+
+
+
+#define NUM_DATA 40
+
+
+struct data
+{
+    size_t n;
+    double *y;
+    double *sigma;
+};
+
+/******************************************************************************
+expb_f(x, params, f): This function performs
+ 
+    x: these are the parameters of the function that must be determined.
+    params: the actual input data is here.
+ params->n the number of input data points
+ params->y the data range
+ params->sigma the errors associated with each data point.
+    f: the result is returned here.  Actually, it is the function evaluated
+ at the data point, minus the expected value of the function, and then
+ divided by the errors.
+ *****************************************************************************/
+int expb_f(const gsl_vector *x,
+           void *params,
+           gsl_vector *f)
+{
+    size_t n = ((struct data *)params)->n;
+    double *y = ((struct data *)params)->y;
+    double *sigma = ((struct data *) params)->sigma;
+    double A;
+    double lambda;
+    double b;
+    int i;
+    double Yi;
+
+    A      = gsl_vector_get(x, 0);
+    lambda = gsl_vector_get(x, 1);
+    b      = gsl_vector_get(x, 2);
+    for (i = 0; i < n; i++) {
+        Yi = A * exp(-lambda * ((double) i)) + b;
+        Yi = (Yi - y[i]) / sigma[i];
+        gsl_vector_set(f, i, Yi);
+        printf("gsl_vector_set(outData, %d, %.1f)\n", i,
+               gsl_vector_get(f, i));
+
+        //        printf("--------- expb_f((%.1f) %.1f %.1f %.1f) is %.1f [%.1f] ---------\n",
+        //               (double) i, A, lambda, b, A * exp(-lambda * ((double) i)) + b,
+        //               Yi);
+    }
+
+    return GSL_SUCCESS;
+}
+
+/******************************************************************************
+ *****************************************************************************/
+int expb_df(const gsl_vector *x,
+            void *params,
+            gsl_matrix *J)
+{
+    size_t n = ((struct data *)params)->n;
+    double *sigma = ((struct data *) params)->sigma;
+    double t, s, e;
+    size_t i, j;
+
+    double A = gsl_vector_get(x, 0);
+    double lambda = gsl_vector_get(x, 1);
+
+
+    for (i = 0; i < n; i++) {
+        // Jacobian matrix J(i,j) = dfi / dxj,
+        // where fi = (Yi - yi)/sigma[i],
+        //       Yi = A * exp(-lambda * i) + b
+        // and the xj are the parameters (A,lambda,b)
+        t = i;
+        s = sigma[i];
+        e = exp(-lambda *t);
+        gsl_matrix_set(J, i, 0, e/s);
+        gsl_matrix_set(J, i, 1, -t * A * e/s);
+        gsl_matrix_set(J, i, 2, 1/s);
+        //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
+        //              (double) i, A, lambda, gsl_vector_get(x, 2), 0, e, e/s);
+        //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
+        //              (double) i, A, lambda, gsl_vector_get(x, 2), 1, -t * A * e, -t * A * e/s);
+        //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
+        //              (double) i, A, lambda, gsl_vector_get(x, 2), 2, 1.0, 1.0/s);
+        for (j=0;j<3;j++) {
+            printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j,
+                   gsl_matrix_get(J, i, j));
+        }
+
+
+    }
+
+    return GSL_SUCCESS;
+}
+
+/******************************************************************************
+ *****************************************************************************/
+int expb_fdf(const gsl_vector *x,
+             void *params,
+             gsl_vector *f, gsl_matrix *J)
+{
+    expb_f (x, params, f);
+    expb_df (x, params, J);
+
+    return GSL_SUCCESS;
+}
+
+
+
 
 
@@ -104 +217 @@
 
     if (mask != NULL) {
+        printf("--------------- WHOA ---------------: MASK NOT NULL\n");
         j = 0;
         for (i=0;i<mask->n;i++) {
@@ -122 +236 @@
 
     for (i=0;i<domain->numRows;i++) {
-        printf("Data item %d is ( ", i);
+        //        printf("Data item %d is ( ", i);
         for (j=0;j<domain->numCols;j++) {
             tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
-            printf("%.1f ", tmpVecPtr->data.F32[j]);
+            //            printf("%.1f ", tmpVecPtr->data.F32[j]);
         }
         tmpf = evalModel(tmpVecPtr, inputParameterList);
-        printf(" ).  Output is %.1f\n", tmpf);
+        //        printf(" ).  Output is %.1f\n", tmpf);
 
         gsl_vector_set(outData, i, (tmpf - data->data.F32[i])/
                        errors->data.F32[i]);
+        printf("gsl_vector_set(outData, %d, %.1f)\n", i,
+               gsl_vector_get(outData, i));
+
         //printf("--------- MYFUNC((%.1f) %.1f %.1f %.1f) is [%.1f] ---------\n",
         //tmpVecPtr->data.F32[0],
@@ -138 +255 @@
         //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));
     }
 
@@ -176 +287 @@
 
     if (mask != NULL) {
+        printf("--------------- WHOA ---------------: MASK NOT NULL\n");
         j = 0;
         for (i=0;i<mask->n;i++) {
@@ -207 +319 @@
 
             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);
@@ -235 +339 @@
     return GSL_SUCCESS;
 }
+
+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));
+    return 0;
+}
+
+#define FIT(i) gsl_vector_get(s->x, i)
+#define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
+
+
+int try03()
+{
+    const gsl_multifit_fdfsolver_type *T;
+    gsl_multifit_fdfsolver *s;
+    int status;
+    size_t i, iter = 0;
+    const size_t n = NUM_DATA;
+    const size_t p = 3;
+    gsl_multifit_function_fdf f;
+    double y[NUM_DATA], sigma[NUM_DATA];
+    gsl_matrix *covar = gsl_matrix_alloc (p, p);
+    struct data d = {
+                        n, y, sigma
+                    };
+    double x_init[3] = { 1.0, 0.0, 0.0 };
+    gsl_vector_view x = gsl_vector_view_array (x_init, p);
+
+    printf("Calling try03()\n");
+    const gsl_rng_type *type;
+    gsl_rng *r;
+
+    gsl_rng_env_setup();
+
+    type = gsl_rng_default;
+    r = gsl_rng_alloc (type);
+
+    f.f = &expb_f;
+    f.df = &expb_df;
+    f.fdf = &expb_fdf;
+    f.n = n;   // Equals N, equals the number of data points.
+    f.p = p;   // Equals 3
+    f.params = &d;
+
+    /* This is the data to be fitted */
+
+    for (i = 0; i < n; i++) {
+        double t = i;
+        y[i] = 1.0 + 5 * exp (-0.1 * t)
+               + gsl_ran_gaussian (r, 0.1);
+        sigma[i] = 0.1;
+        //      printf ("data: %d %g %g\n", i, y[i], sigma[i]);
+        //      printf("data->data.F32[%d][0] = %f;\n", i, y[i]);
+    };
+
+    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);
+        if (status)
+            break;
+        status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
+    } while (status == GSL_CONTINUE && iter < 500);
+
+    gsl_multifit_covar (s->J, 0.0, covar);
+
+    gsl_matrix_fprintf (stdout, covar, "%g");
+
+    printf ("A      = %.5f +/- %.5f\n", FIT(0), ERR(0));
+    printf ("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
+    printf ("b      = %.5f +/- %.5f\n", FIT(2), ERR(2));
+
+    {
+        double chi = gsl_blas_dnrm2(s->f);
+        printf("chisq/dof = %g\n",  pow(chi, 2.0)/ (n - p));
+    }
+
+    printf ("status = %s\n", gsl_strerror (status));
+
+    gsl_multifit_fdfsolver_free (s);
+    printf("Called try03()\n");
+    return 0;
+}
+
+
 
 
@@ -252 +453 @@
 {
     printf("Calling psMinimizeChi2()\n");
+    try03();
 
     int numData = domain->numRows; // Number of data points
@@ -295 +497 @@
     xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
     if (paramMask != NULL) {
+        printf("--------------- WHOA ---------------: MASK NOT NULL\n");
         j = 0;
         for (i=0;i<initialGuess->n;i++) {
@@ -330 +533 @@
     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.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.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.
-
     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.
     do {
         iter++;
@@ -360 +552 @@
 
         status = gsl_multifit_fdfsolver_iterate(s);
+        printf ("psMinimize() status = %s\n", gsl_strerror (status));
         // If there was a problem, abort.
         if (status) {
@@ -375 +568 @@
     // into the solution.
     if (paramMask != NULL) {
+        printf("--------------- WHOA ---------------: MASK NOT NULL\n");
         j = 0;
         for (i=0;i<initialGuess->n;i++) {
@@ -583 +777 @@
  * success.
  */
+
+
+
 void polyOrderCheck(float **A,
                     int N,