Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 1123) +++ trunk/psLib/src/math/psMinimize.c (revision 1128) @@ -5,4 +5,10 @@ #include #include + +#include +#include +#include +#include +#include #include "psMemory.h" @@ -18,10 +24,4 @@ #include -#include -#include -#include -#include -#include - #define MAX_LMM_NUM_ITERATIONS 500 typedef struct @@ -45,4 +45,117 @@ // 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,4 +217,5 @@ if (mask != NULL) { + printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) { @@ -122,14 +236,17 @@ for (i=0;inumRows;i++) { - printf("Data item %d is ( ", i); + // printf("Data item %d is ( ", i); for (j=0;jnumCols;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,10 +255,4 @@ //inputParameterList->data.F32[2], //(tmpf - data->data.F32[i])/errors->data.F32[i]); - - } - - for (i=0;inumRows;i++) { - printf("gsl_vector_set(outData, %d, %.1f)\n", i, - gsl_vector_get(outData, i)); } @@ -176,4 +287,5 @@ if (mask != NULL) { + printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) { @@ -207,16 +319,8 @@ gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i])); - } - } - - - for (i=0;inumRows;i++) { - for (j=0;jn;j++) { printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j, gsl_matrix_get(J, i, j)); } } - - psFree(inputParameterList); @@ -235,4 +339,101 @@ 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,4 +453,5 @@ { printf("Calling psMinimizeChi2()\n"); + try03(); int numData = domain->numRows; // Number of data points @@ -295,4 +497,5 @@ xInit = (double *) psAlloc(inputData.paramCount * sizeof(double)); if (paramMask != NULL) { + printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) { @@ -330,21 +533,10 @@ 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,4 +552,5 @@ status = gsl_multifit_fdfsolver_iterate(s); + printf ("psMinimize() status = %s\n", gsl_strerror (status)); // If there was a problem, abort. if (status) { @@ -375,4 +568,5 @@ // into the solution. if (paramMask != NULL) { + printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) { @@ -583,4 +777,7 @@ * success. */ + + + void polyOrderCheck(float **A, int N,