Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 1129) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 1176) @@ -45,122 +45,4 @@ // 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; - - for (i=0;i<3;i++) { - printf("(expb_f(): gsl_vector_get(x, %d) is %.1f\n", i, - gsl_vector_get(x, i)); - } - - - 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; -} - - @@ -216,12 +98,4 @@ tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32); - for (i=0;i<3;i++) { - printf("(FART: gsl_vector_get(x, %d) is %.1f\n", i, - gsl_vector_get(x, i)); - } - - - - // The GSL routines will call this function with the masked parameters // removed. However, the user-supplied function (to be modified) does not @@ -231,5 +105,4 @@ if (mask != NULL) { - printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) { @@ -245,28 +118,13 @@ } } - for (i=0;in;i++) { - printf("(gsl_f(): inputParameterList->data.F32[%d] is %.1f\n", i, inputParameterList->data.F32[i]); - } for (i=0;inumRows;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]); } 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("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], - //inputParameterList->data.F32[0], - //inputParameterList->data.F32[1], - //inputParameterList->data.F32[2], - //(tmpf - data->data.F32[i])/errors->data.F32[i]); } @@ -301,5 +159,4 @@ if (mask != NULL) { - printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) { @@ -323,16 +180,5 @@ for (j=0;jn;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])); - printf("gsl_matrix_set(J, %d, %d, %.1f)\n", i, j, - gsl_matrix_get(J, i, j)); } } @@ -353,103 +199,4 @@ 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; -} - - - - /****************************************************************************** @@ -466,7 +213,4 @@ float *chiSq) { - printf("Calling psMinimizeChi2()\n"); - try03(); - int numData = domain->numRows; // Number of data points int status; // Return status for the GSL solver. @@ -484,4 +228,5 @@ gsl_multifit_fdfsolver *s; // GSL data structure. psModelData inputData; + float chiSqOld = 0.0; inputData.n = numData; @@ -511,9 +256,7 @@ xInit = (double *) psAlloc(inputData.paramCount * sizeof(double)); if (paramMask != NULL) { - printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) { if (paramMask->data.U8[i] == 0) { - printf("xInit[%d] is %f (masked loop)\n", j, initialGuess->data.F32[i]); xInit[j++] = initialGuess->data.F32[i]; } @@ -522,5 +265,4 @@ for (i=0;in;i++) { xInit[i] = initialGuess->data.F32[i]; - printf("xInit[%d] is %f\n", i, initialGuess->data.F32[i]); } } @@ -545,34 +287,35 @@ f.params = &inputData; - printf("inputData.paramCount is %d\n", inputData.paramCount); - printf("numData is %d\n", numData); - - gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount); T = gsl_multifit_fdfsolver_lmsder; s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount); gsl_multifit_fdfsolver_set(s, &f, &x.vector); + *chiSq = 0.0; + chiSqOld = 0.0; 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. - status = gsl_multifit_fdfsolver_iterate(s); - printf ("psMinimize() status = %s\n", gsl_strerror (status)); + // printf ("psMinimize() status = %s\n", gsl_strerror(status)); // If there was a problem, abort. if (status) { - printf("ABORT: status is %d\n", status); - // psAbort(__func__, "gsl_multifit_fdfsolver_iterate()\n"); + 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); + // NOTE: This wasn't working right when the parameters fit exactly. + // Figure out why. + // status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4); + + // We test for convergence if chiSquared changes by less than 1.0 + // as specified in the ADD. + *chiSq = gsl_blas_dnrm2(s->f); + if (fabs(*chiSq - chiSqOld) < 1.0) { + status = GSL_SUCCESS; + } else { + status = GSL_CONTINUE; + } + chiSqOld = *chiSq; + } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS); @@ -582,5 +325,4 @@ // into the solution. if (paramMask != NULL) { - printf("--------------- WHOA ---------------: MASK NOT NULL\n"); j = 0; for (i=0;in;i++) {