IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1176


Ignore:
Timestamp:
Jul 1, 2004, 12:22:31 PM (22 years ago)
Author:
gusciora
Message:

Does it work? Does it not work? Nobody knows.

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psMinimize.c

    r1129 r1176  
    4545// The second argument must have the same size
    4646// as *initialGuess and *paramMask.
    47 
    48 
    49 
    50 #define NUM_DATA 40
    51 
    52 
    53 struct data
    54 {
    55     size_t n;
    56     double *y;
    57     double *sigma;
    58 };
    59 
    60 /******************************************************************************
    61 expb_f(x, params, f): This function performs
    62  
    63     x: these are the parameters of the function that must be determined.
    64     params: the actual input data is here.
    65  params->n the number of input data points
    66  params->y the data range
    67  params->sigma the errors associated with each data point.
    68     f: the result is returned here.  Actually, it is the function evaluated
    69  at the data point, minus the expected value of the function, and then
    70  divided by the errors.
    71  *****************************************************************************/
    72 int expb_f(const gsl_vector *x,
    73            void *params,
    74            gsl_vector *f)
    75 {
    76     size_t n = ((struct data *)params)->n;
    77     double *y = ((struct data *)params)->y;
    78     double *sigma = ((struct data *) params)->sigma;
    79     double A;
    80     double lambda;
    81     double b;
    82     int i;
    83     double Yi;
    84 
    85     for (i=0;i<3;i++) {
    86         printf("(expb_f(): gsl_vector_get(x, %d) is %.1f\n", i,
    87                gsl_vector_get(x, i));
    88     }
    89 
    90 
    91     A      = gsl_vector_get(x, 0);
    92     lambda = gsl_vector_get(x, 1);
    93     b      = gsl_vector_get(x, 2);
    94     for (i = 0; i < n; i++) {
    95         Yi = A * exp(-lambda * ((double) i)) + b;
    96         Yi = (Yi - y[i]) / sigma[i];
    97         gsl_vector_set(f, i, Yi);
    98         printf("gsl_vector_set(outData, %d, %.1f)\n", i,
    99                gsl_vector_get(f, i));
    100 
    101         //        printf("--------- expb_f((%.1f) %.1f %.1f %.1f) is %.1f [%.1f] ---------\n",
    102         //               (double) i, A, lambda, b, A * exp(-lambda * ((double) i)) + b,
    103         //               Yi);
    104     }
    105 
    106     return GSL_SUCCESS;
    107 }
    108 
    109 /******************************************************************************
    110  *****************************************************************************/
    111 int expb_df(const gsl_vector *x,
    112             void *params,
    113             gsl_matrix *J)
    114 {
    115     size_t n = ((struct data *)params)->n;
    116     double *sigma = ((struct data *) params)->sigma;
    117     double t, s, e;
    118     size_t i, j;
    119 
    120     double A = gsl_vector_get(x, 0);
    121     double lambda = gsl_vector_get(x, 1);
    122 
    123 
    124     for (i = 0; i < n; i++) {
    125         // Jacobian matrix J(i,j) = dfi / dxj,
    126         // where fi = (Yi - yi)/sigma[i],
    127         //       Yi = A * exp(-lambda * i) + b
    128         // and the xj are the parameters (A,lambda,b)
    129         t = i;
    130         s = sigma[i];
    131         e = exp(-lambda *t);
    132         gsl_matrix_set(J, i, 0, e/s);
    133         gsl_matrix_set(J, i, 1, -t * A * e/s);
    134         gsl_matrix_set(J, i, 2, 1/s);
    135         //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
    136         //              (double) i, A, lambda, gsl_vector_get(x, 2), 0, e, e/s);
    137         //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
    138         //              (double) i, A, lambda, gsl_vector_get(x, 2), 1, -t * A * e, -t * A * e/s);
    139         //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
    140         //              (double) i, A, lambda, gsl_vector_get(x, 2), 2, 1.0, 1.0/s);
    141         for (j=0;j<3;j++) {
    142             printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j,
    143                    gsl_matrix_get(J, i, j));
    144         }
    145 
    146 
    147     }
    148 
    149     return GSL_SUCCESS;
    150 }
    151 
    152 /******************************************************************************
    153  *****************************************************************************/
    154 int expb_fdf(const gsl_vector *x,
    155              void *params,
    156              gsl_vector *f, gsl_matrix *J)
    157 {
    158     expb_f (x, params, f);
    159     expb_df (x, params, J);
    160 
    161     return GSL_SUCCESS;
    162 }
    163 
    164 
    16547
    16648
    … …  
    21698    tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    21799
    218     for (i=0;i<3;i++) {
    219         printf("(FART: gsl_vector_get(x, %d) is %.1f\n", i,
    220                gsl_vector_get(x, i));
    221     }
    222 
    223 
    224 
    225 
    226100    // The GSL routines will call this function with the masked parameters
    227101    // removed.  However, the user-supplied function (to be modified) does not
    … …  
    231105
    232106    if (mask != NULL) {
    233         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    234107        j = 0;
    235108        for (i=0;i<mask->n;i++) {
    … …  
    245118        }
    246119    }
    247     for (i=0;i<inputParameterList->n;i++) {
    248         printf("(gsl_f(): inputParameterList->data.F32[%d] is %.1f\n", i, inputParameterList->data.F32[i]);
    249     }
    250120
    251121    for (i=0;i<domain->numRows;i++) {
    252         //        printf("Data item %d is ( ", i);
    253122        for (j=0;j<domain->numCols;j++) {
    254123            tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
    255             //            printf("%.1f ", tmpVecPtr->data.F32[j]);
    256124        }
    257125        tmpf = evalModel(tmpVecPtr, inputParameterList);
    258         //        printf(" ).  Output is %.1f\n", tmpf);
    259126
    260127        gsl_vector_set(outData, i, (tmpf - data->data.F32[i])/
    261128                       errors->data.F32[i]);
    262         printf("gsl_vector_set(outData, %d, %.1f)\n", i,
    263                gsl_vector_get(outData, i));
    264 
    265         //printf("--------- MYFUNC((%.1f) %.1f %.1f %.1f) is [%.1f] ---------\n",
    266         //tmpVecPtr->data.F32[0],
    267         //inputParameterList->data.F32[0],
    268         //inputParameterList->data.F32[1],
    269         //inputParameterList->data.F32[2],
    270         //(tmpf - data->data.F32[i])/errors->data.F32[i]);
    271129    }
    272130
    … …  
    301159
    302160    if (mask != NULL) {
    303         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    304161        j = 0;
    305162        for (i=0;i<mask->n;i++) {
    … …  
    323180        for (j=0;j<inputParameterList->n;j++) {
    324181            tmpf = d_evalModel(tmpVecPtr, inputParameterList, j);
    325 
    326             //printf("--------- MYFUNCDERIV((%.1f) %.1f %.1f %.1f, %d) is [%.1f] ---------\n",
    327             //tmpVecPtr->data.F32[0],
    328             //inputParameterList->data.F32[0],
    329             //inputParameterList->data.F32[1],
    330             //inputParameterList->data.F32[2],
    331             //j,
    332             //tmpf/errors->data.F32[i]);
    333 
    334182            gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i]));
    335             printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j,
    336                    gsl_matrix_get(J, i, j));
    337183        }
    338184    }
    … …  
    353199    return GSL_SUCCESS;
    354200}
    355 
    356 int print_state(size_t iter,
    357                 gsl_multifit_fdfsolver *s)
    358 {
    359     printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
    360             "|f(x)| = %g\n",
    361             iter,
    362             gsl_vector_get(s->x, 0),
    363             gsl_vector_get(s->x, 1),
    364             gsl_vector_get(s->x, 2),
    365             gsl_blas_dnrm2(s->f));
    366     return 0;
    367 }
    368 
    369 #define FIT(i) gsl_vector_get(s->x, i)
    370 #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
    371 
    372 
    373 int try03()
    374 {
    375     const gsl_multifit_fdfsolver_type *T;
    376     gsl_multifit_fdfsolver *s;
    377     int status;
    378     size_t i, iter = 0;
    379     const size_t n = NUM_DATA;
    380     const size_t p = 3;
    381     gsl_multifit_function_fdf f;
    382     double y[NUM_DATA], sigma[NUM_DATA];
    383     gsl_matrix *covar = gsl_matrix_alloc (p, p);
    384     struct data d = {
    385                         n, y, sigma
    386                     };
    387     double x_init[3] = { 1.0, 0.0, 0.0 };
    388     gsl_vector_view x = gsl_vector_view_array (x_init, p);
    389 
    390     printf("Calling try03()\n");
    391     const gsl_rng_type *type;
    392     gsl_rng *r;
    393 
    394     gsl_rng_env_setup();
    395 
    396     type = gsl_rng_default;
    397     r = gsl_rng_alloc (type);
    398 
    399     f.f = &expb_f;
    400     f.df = &expb_df;
    401     f.fdf = &expb_fdf;
    402     f.n = n;   // Equals N, equals the number of data points.
    403     f.p = p;   // Equals 3
    404     f.params = &d;
    405 
    406     /* This is the data to be fitted */
    407 
    408     for (i = 0; i < n; i++) {
    409         double t = i;
    410         y[i] = 1.0 + 5 * exp (-0.1 * t)
    411                + gsl_ran_gaussian (r, 0.1);
    412         sigma[i] = 0.1;
    413         //      printf ("data: %d %g %g\n", i, y[i], sigma[i]);
    414         //      printf("data->data.F32[%d][0] = %f;\n", i, y[i]);
    415     };
    416 
    417     T = gsl_multifit_fdfsolver_lmsder;
    418     s = gsl_multifit_fdfsolver_alloc (T, n, p);
    419     gsl_multifit_fdfsolver_set(s, &f, &x.vector);
    420     print_state(iter, s);
    421     do {
    422         iter++;
    423         status = gsl_multifit_fdfsolver_iterate (s);
    424         printf ("status = %s\n", gsl_strerror (status));
    425         print_state (iter, s);
    426         if (status)
    427             break;
    428         status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
    429     } while (status == GSL_CONTINUE && iter < 500);
    430 
    431     gsl_multifit_covar (s->J, 0.0, covar);
    432 
    433     gsl_matrix_fprintf (stdout, covar, "%g");
    434 
    435     printf ("A      = %.5f +/- %.5f\n", FIT(0), ERR(0));
    436     printf ("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
    437     printf ("b      = %.5f +/- %.5f\n", FIT(2), ERR(2));
    438 
    439     {
    440         double chi = gsl_blas_dnrm2(s->f);
    441         printf("chisq/dof = %g\n",  pow(chi, 2.0)/ (n - p));
    442     }
    443 
    444     printf ("status = %s\n", gsl_strerror (status));
    445 
    446     gsl_multifit_fdfsolver_free (s);
    447     printf("Called try03()\n");
    448     return 0;
    449 }
    450 
    451 
    452 
    453 
    454201
    455202/******************************************************************************
    … …  
    466213               float *chiSq)
    467214{
    468     printf("Calling psMinimizeChi2()\n");
    469     try03();
    470 
    471215    int numData = domain->numRows; // Number of data points
    472216    int status;    // Return status for the GSL solver.
    … …  
    484228    gsl_multifit_fdfsolver *s; // GSL data structure.
    485229    psModelData inputData;
     230    float chiSqOld = 0.0;
    486231
    487232    inputData.n = numData;
    … …  
    511256    xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
    512257    if (paramMask != NULL) {
    513         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    514258        j = 0;
    515259        for (i=0;i<initialGuess->n;i++) {
    516260            if (paramMask->data.U8[i] == 0) {
    517                 printf("xInit[%d] is %f (masked loop)\n", j, initialGuess->data.F32[i]);
    518261                xInit[j++] = initialGuess->data.F32[i];
    519262            }
    … …  
    522265        for (i=0;i<initialGuess->n;i++) {
    523266            xInit[i] = initialGuess->data.F32[i];
    524             printf("xInit[%d] is %f\n", i, initialGuess->data.F32[i]);
    525267        }
    526268    }
    … …  
    545287    f.params = &inputData;
    546288
    547     printf("inputData.paramCount is %d\n", inputData.paramCount);
    548     printf("numData is %d\n", numData);
    549 
    550 
    551289    gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount);
    552290    T = gsl_multifit_fdfsolver_lmsder;
    553291    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount);
    554292    gsl_multifit_fdfsolver_set(s, &f, &x.vector);
     293    *chiSq = 0.0;
     294    chiSqOld = 0.0;
    555295    do {
    556296        iter++;
    557         printf("##################################################\n");
    558         printf("##################################################\n");
    559         printf("##################################################\n");
    560         printf("################## Iteration %d ##################\n", iter);
    561         printf("##################################################\n");
    562         printf("##################################################\n");
    563         printf("##################################################\n");
    564 
    565297        // Perform an iteration of the GSL solver.
    566 
    567298        status = gsl_multifit_fdfsolver_iterate(s);
    568         printf ("psMinimize() status = %s\n", gsl_strerror (status));
     299        //        printf ("psMinimize() status = %s\n", gsl_strerror(status));
    569300        // If there was a problem, abort.
    570301        if (status) {
    571             printf("ABORT: status is %d\n", status);
    572             //            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
     302            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
    573303        }
    574304
    575305        // Test if the parameters changed by a small enough amount.
    576         status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
     306        // NOTE: This wasn't working right when the parameters fit exactly.
     307        // Figure out why.
     308        //        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
     309
     310        // We test for convergence if chiSquared changes by less than 1.0
     311        // as specified in the ADD.
     312        *chiSq = gsl_blas_dnrm2(s->f);
     313        if (fabs(*chiSq - chiSqOld) < 1.0) {
     314            status = GSL_SUCCESS;
     315        } else {
     316            status = GSL_CONTINUE;
     317        }
     318        chiSqOld = *chiSq;
     319
    577320    } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS);
    578321
    … …  
    582325    // into the solution.
    583326    if (paramMask != NULL) {
    584         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    585327        j = 0;
    586328        for (i=0;i<initialGuess->n;i++) {

  • trunk/psLib/src/math/psMinimize.c

    r1129 r1176  
    4545// The second argument must have the same size
    4646// as *initialGuess and *paramMask.
    47 
    48 
    49 
    50 #define NUM_DATA 40
    51 
    52 
    53 struct data
    54 {
    55     size_t n;
    56     double *y;
    57     double *sigma;
    58 };
    59 
    60 /******************************************************************************
    61 expb_f(x, params, f): This function performs
    62  
    63     x: these are the parameters of the function that must be determined.
    64     params: the actual input data is here.
    65  params->n the number of input data points
    66  params->y the data range
    67  params->sigma the errors associated with each data point.
    68     f: the result is returned here.  Actually, it is the function evaluated
    69  at the data point, minus the expected value of the function, and then
    70  divided by the errors.
    71  *****************************************************************************/
    72 int expb_f(const gsl_vector *x,
    73            void *params,
    74            gsl_vector *f)
    75 {
    76     size_t n = ((struct data *)params)->n;
    77     double *y = ((struct data *)params)->y;
    78     double *sigma = ((struct data *) params)->sigma;
    79     double A;
    80     double lambda;
    81     double b;
    82     int i;
    83     double Yi;
    84 
    85     for (i=0;i<3;i++) {
    86         printf("(expb_f(): gsl_vector_get(x, %d) is %.1f\n", i,
    87                gsl_vector_get(x, i));
    88     }
    89 
    90 
    91     A      = gsl_vector_get(x, 0);
    92     lambda = gsl_vector_get(x, 1);
    93     b      = gsl_vector_get(x, 2);
    94     for (i = 0; i < n; i++) {
    95         Yi = A * exp(-lambda * ((double) i)) + b;
    96         Yi = (Yi - y[i]) / sigma[i];
    97         gsl_vector_set(f, i, Yi);
    98         printf("gsl_vector_set(outData, %d, %.1f)\n", i,
    99                gsl_vector_get(f, i));
    100 
    101         //        printf("--------- expb_f((%.1f) %.1f %.1f %.1f) is %.1f [%.1f] ---------\n",
    102         //               (double) i, A, lambda, b, A * exp(-lambda * ((double) i)) + b,
    103         //               Yi);
    104     }
    105 
    106     return GSL_SUCCESS;
    107 }
    108 
    109 /******************************************************************************
    110  *****************************************************************************/
    111 int expb_df(const gsl_vector *x,
    112             void *params,
    113             gsl_matrix *J)
    114 {
    115     size_t n = ((struct data *)params)->n;
    116     double *sigma = ((struct data *) params)->sigma;
    117     double t, s, e;
    118     size_t i, j;
    119 
    120     double A = gsl_vector_get(x, 0);
    121     double lambda = gsl_vector_get(x, 1);
    122 
    123 
    124     for (i = 0; i < n; i++) {
    125         // Jacobian matrix J(i,j) = dfi / dxj,
    126         // where fi = (Yi - yi)/sigma[i],
    127         //       Yi = A * exp(-lambda * i) + b
    128         // and the xj are the parameters (A,lambda,b)
    129         t = i;
    130         s = sigma[i];
    131         e = exp(-lambda *t);
    132         gsl_matrix_set(J, i, 0, e/s);
    133         gsl_matrix_set(J, i, 1, -t * A * e/s);
    134         gsl_matrix_set(J, i, 2, 1/s);
    135         //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
    136         //              (double) i, A, lambda, gsl_vector_get(x, 2), 0, e, e/s);
    137         //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
    138         //              (double) i, A, lambda, gsl_vector_get(x, 2), 1, -t * A * e, -t * A * e/s);
    139         //      printf("--------- myFuncDeriv((%.1f) %.1f %.1f %.1f (%d)) is %.1f [%.1f] ---------\n",
    140         //              (double) i, A, lambda, gsl_vector_get(x, 2), 2, 1.0, 1.0/s);
    141         for (j=0;j<3;j++) {
    142             printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j,
    143                    gsl_matrix_get(J, i, j));
    144         }
    145 
    146 
    147     }
    148 
    149     return GSL_SUCCESS;
    150 }
    151 
    152 /******************************************************************************
    153  *****************************************************************************/
    154 int expb_fdf(const gsl_vector *x,
    155              void *params,
    156              gsl_vector *f, gsl_matrix *J)
    157 {
    158     expb_f (x, params, f);
    159     expb_df (x, params, J);
    160 
    161     return GSL_SUCCESS;
    162 }
    163 
    164 
    16547
    16648
    … …  
    21698    tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    21799
    218     for (i=0;i<3;i++) {
    219         printf("(FART: gsl_vector_get(x, %d) is %.1f\n", i,
    220                gsl_vector_get(x, i));
    221     }
    222 
    223 
    224 
    225 
    226100    // The GSL routines will call this function with the masked parameters
    227101    // removed.  However, the user-supplied function (to be modified) does not
    … …  
    231105
    232106    if (mask != NULL) {
    233         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    234107        j = 0;
    235108        for (i=0;i<mask->n;i++) {
    … …  
    245118        }
    246119    }
    247     for (i=0;i<inputParameterList->n;i++) {
    248         printf("(gsl_f(): inputParameterList->data.F32[%d] is %.1f\n", i, inputParameterList->data.F32[i]);
    249     }
    250120
    251121    for (i=0;i<domain->numRows;i++) {
    252         //        printf("Data item %d is ( ", i);
    253122        for (j=0;j<domain->numCols;j++) {
    254123            tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
    255             //            printf("%.1f ", tmpVecPtr->data.F32[j]);
    256124        }
    257125        tmpf = evalModel(tmpVecPtr, inputParameterList);
    258         //        printf(" ).  Output is %.1f\n", tmpf);
    259126
    260127        gsl_vector_set(outData, i, (tmpf - data->data.F32[i])/
    261128                       errors->data.F32[i]);
    262         printf("gsl_vector_set(outData, %d, %.1f)\n", i,
    263                gsl_vector_get(outData, i));
    264 
    265         //printf("--------- MYFUNC((%.1f) %.1f %.1f %.1f) is [%.1f] ---------\n",
    266         //tmpVecPtr->data.F32[0],
    267         //inputParameterList->data.F32[0],
    268         //inputParameterList->data.F32[1],
    269         //inputParameterList->data.F32[2],
    270         //(tmpf - data->data.F32[i])/errors->data.F32[i]);
    271129    }
    272130
    … …  
    301159
    302160    if (mask != NULL) {
    303         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    304161        j = 0;
    305162        for (i=0;i<mask->n;i++) {
    … …  
    323180        for (j=0;j<inputParameterList->n;j++) {
    324181            tmpf = d_evalModel(tmpVecPtr, inputParameterList, j);
    325 
    326             //printf("--------- MYFUNCDERIV((%.1f) %.1f %.1f %.1f, %d) is [%.1f] ---------\n",
    327             //tmpVecPtr->data.F32[0],
    328             //inputParameterList->data.F32[0],
    329             //inputParameterList->data.F32[1],
    330             //inputParameterList->data.F32[2],
    331             //j,
    332             //tmpf/errors->data.F32[i]);
    333 
    334182            gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i]));
    335             printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j,
    336                    gsl_matrix_get(J, i, j));
    337183        }
    338184    }
    … …  
    353199    return GSL_SUCCESS;
    354200}
    355 
    356 int print_state(size_t iter,
    357                 gsl_multifit_fdfsolver *s)
    358 {
    359     printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
    360             "|f(x)| = %g\n",
    361             iter,
    362             gsl_vector_get(s->x, 0),
    363             gsl_vector_get(s->x, 1),
    364             gsl_vector_get(s->x, 2),
    365             gsl_blas_dnrm2(s->f));
    366     return 0;
    367 }
    368 
    369 #define FIT(i) gsl_vector_get(s->x, i)
    370 #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
    371 
    372 
    373 int try03()
    374 {
    375     const gsl_multifit_fdfsolver_type *T;
    376     gsl_multifit_fdfsolver *s;
    377     int status;
    378     size_t i, iter = 0;
    379     const size_t n = NUM_DATA;
    380     const size_t p = 3;
    381     gsl_multifit_function_fdf f;
    382     double y[NUM_DATA], sigma[NUM_DATA];
    383     gsl_matrix *covar = gsl_matrix_alloc (p, p);
    384     struct data d = {
    385                         n, y, sigma
    386                     };
    387     double x_init[3] = { 1.0, 0.0, 0.0 };
    388     gsl_vector_view x = gsl_vector_view_array (x_init, p);
    389 
    390     printf("Calling try03()\n");
    391     const gsl_rng_type *type;
    392     gsl_rng *r;
    393 
    394     gsl_rng_env_setup();
    395 
    396     type = gsl_rng_default;
    397     r = gsl_rng_alloc (type);
    398 
    399     f.f = &expb_f;
    400     f.df = &expb_df;
    401     f.fdf = &expb_fdf;
    402     f.n = n;   // Equals N, equals the number of data points.
    403     f.p = p;   // Equals 3
    404     f.params = &d;
    405 
    406     /* This is the data to be fitted */
    407 
    408     for (i = 0; i < n; i++) {
    409         double t = i;
    410         y[i] = 1.0 + 5 * exp (-0.1 * t)
    411                + gsl_ran_gaussian (r, 0.1);
    412         sigma[i] = 0.1;
    413         //      printf ("data: %d %g %g\n", i, y[i], sigma[i]);
    414         //      printf("data->data.F32[%d][0] = %f;\n", i, y[i]);
    415     };
    416 
    417     T = gsl_multifit_fdfsolver_lmsder;
    418     s = gsl_multifit_fdfsolver_alloc (T, n, p);
    419     gsl_multifit_fdfsolver_set(s, &f, &x.vector);
    420     print_state(iter, s);
    421     do {
    422         iter++;
    423         status = gsl_multifit_fdfsolver_iterate (s);
    424         printf ("status = %s\n", gsl_strerror (status));
    425         print_state (iter, s);
    426         if (status)
    427             break;
    428         status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);
    429     } while (status == GSL_CONTINUE && iter < 500);
    430 
    431     gsl_multifit_covar (s->J, 0.0, covar);
    432 
    433     gsl_matrix_fprintf (stdout, covar, "%g");
    434 
    435     printf ("A      = %.5f +/- %.5f\n", FIT(0), ERR(0));
    436     printf ("lambda = %.5f +/- %.5f\n", FIT(1), ERR(1));
    437     printf ("b      = %.5f +/- %.5f\n", FIT(2), ERR(2));
    438 
    439     {
    440         double chi = gsl_blas_dnrm2(s->f);
    441         printf("chisq/dof = %g\n",  pow(chi, 2.0)/ (n - p));
    442     }
    443 
    444     printf ("status = %s\n", gsl_strerror (status));
    445 
    446     gsl_multifit_fdfsolver_free (s);
    447     printf("Called try03()\n");
    448     return 0;
    449 }
    450 
    451 
    452 
    453 
    454201
    455202/******************************************************************************
    … …  
    466213               float *chiSq)
    467214{
    468     printf("Calling psMinimizeChi2()\n");
    469     try03();
    470 
    471215    int numData = domain->numRows; // Number of data points
    472216    int status;    // Return status for the GSL solver.
    … …  
    484228    gsl_multifit_fdfsolver *s; // GSL data structure.
    485229    psModelData inputData;
     230    float chiSqOld = 0.0;
    486231
    487232    inputData.n = numData;
    … …  
    511256    xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
    512257    if (paramMask != NULL) {
    513         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    514258        j = 0;
    515259        for (i=0;i<initialGuess->n;i++) {
    516260            if (paramMask->data.U8[i] == 0) {
    517                 printf("xInit[%d] is %f (masked loop)\n", j, initialGuess->data.F32[i]);
    518261                xInit[j++] = initialGuess->data.F32[i];
    519262            }
    … …  
    522265        for (i=0;i<initialGuess->n;i++) {
    523266            xInit[i] = initialGuess->data.F32[i];
    524             printf("xInit[%d] is %f\n", i, initialGuess->data.F32[i]);
    525267        }
    526268    }
    … …  
    545287    f.params = &inputData;
    546288
    547     printf("inputData.paramCount is %d\n", inputData.paramCount);
    548     printf("numData is %d\n", numData);
    549 
    550 
    551289    gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount);
    552290    T = gsl_multifit_fdfsolver_lmsder;
    553291    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount);
    554292    gsl_multifit_fdfsolver_set(s, &f, &x.vector);
     293    *chiSq = 0.0;
     294    chiSqOld = 0.0;
    555295    do {
    556296        iter++;
    557         printf("##################################################\n");
    558         printf("##################################################\n");
    559         printf("##################################################\n");
    560         printf("################## Iteration %d ##################\n", iter);
    561         printf("##################################################\n");
    562         printf("##################################################\n");
    563         printf("##################################################\n");
    564 
    565297        // Perform an iteration of the GSL solver.
    566 
    567298        status = gsl_multifit_fdfsolver_iterate(s);
    568         printf ("psMinimize() status = %s\n", gsl_strerror (status));
     299        //        printf ("psMinimize() status = %s\n", gsl_strerror(status));
    569300        // If there was a problem, abort.
    570301        if (status) {
    571             printf("ABORT: status is %d\n", status);
    572             //            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
     302            psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n");
    573303        }
    574304
    575305        // Test if the parameters changed by a small enough amount.
    576         status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
     306        // NOTE: This wasn't working right when the parameters fit exactly.
     307        // Figure out why.
     308        //        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
     309
     310        // We test for convergence if chiSquared changes by less than 1.0
     311        // as specified in the ADD.
     312        *chiSq = gsl_blas_dnrm2(s->f);
     313        if (fabs(*chiSq - chiSqOld) < 1.0) {
     314            status = GSL_SUCCESS;
     315        } else {
     316            status = GSL_CONTINUE;
     317        }
     318        chiSqOld = *chiSq;
     319
    577320    } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS);
    578321
    … …  
    582325    // into the solution.
    583326    if (paramMask != NULL) {
    584         printf("--------------- WHOA ---------------: MASK NOT NULL\n");
    585327        j = 0;
    586328        for (i=0;i<initialGuess->n;i++) {
Note: See TracChangeset for help on using the changeset viewer.