diff options
Diffstat (limited to 'src/fit.c')
| -rw-r--r-- | src/fit.c | 72 |
1 files changed, 63 insertions, 9 deletions
@@ -27,6 +27,26 @@ typedef struct fitParams { double epsrel; } fitParams; +/* In order to start the fitter close to the minimum, we first do a series of + * "quick" minimizations starting at the following points. We keep track of the + * parameters with the best likelihood value and then start the "real" + * minimization from those parameters. */ + +static struct startingParameters { + double x; + double y; + double z; + double theta; + double phi; +} startingParameters[] = { + {0.0, 0.0, 0.0, 0.0, 0.0}, + {0.0, 0.0, 0.0, M_PI/2, 0.0}, + {0.0, 0.0, 0.0, M_PI/2, M_PI/2}, + {0.0, 0.0, 0.0, M_PI/2, M_PI}, + {0.0, 0.0, 0.0, M_PI/2, 3*M_PI/2}, + {0.0, 0.0, 0.0, M_PI, 0.0}, +}; + double nll(unsigned int n, const double *x, double *grad, void *params) { fitParams *fpars = (fitParams *) params; @@ -59,7 +79,7 @@ double nll(unsigned int n, const double *x, double *grad, void *params) long long elapsed = (tv_stop.tv_sec - tv_start.tv_sec)*1000 + (tv_stop.tv_usec - tv_start.tv_usec)/1000; - printf("%5zu %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f %5.2f %5.2f f() = %7.3e took %lld ms\n", + printf("%5zu %10.2f %7.2f %7.2f %7.2f %5.2f %5.2f %5.2f %5.2f %5.2f f() = %7.3e took %lld ms\n", iter++, x[0], x[1], @@ -78,13 +98,13 @@ double nll(unsigned int n, const double *x, double *grad, void *params) int fit_event(event *ev, double *xopt, double *fmin) { + size_t i; fitParams fpars; double x[9], ss[9], lb[9], ub[9], fval, n; int rv; nlopt_opt opt = nlopt_create(NLOPT_LN_BOBYQA, 9); nlopt_set_min_objective(opt,nll,&fpars); - nlopt_set_ftol_abs(opt, 1e-5); x[0] = 1200.0; x[1] = 0.0; @@ -112,8 +132,8 @@ int fit_event(event *ev, double *xopt, double *fmin) lb[1] = -1000.0; lb[2] = -1000.0; lb[3] = -1000.0; - lb[4] = 0.0; - lb[5] = 0.0; + lb[4] = -INFINITY; + lb[5] = -INFINITY; lb[6] = 0.0; lb[7] = -10.0; lb[8] = -10.0; @@ -122,8 +142,8 @@ int fit_event(event *ev, double *xopt, double *fmin) ub[1] = 1000.0; ub[2] = 1000.0; ub[3] = 1000.0; - ub[4] = M_PI; - ub[5] = 2*M_PI; + ub[4] = INFINITY; + ub[5] = INFINITY; ub[6] = 400.0; ub[7] = 10.0; ub[8] = 10.0; @@ -134,10 +154,44 @@ int fit_event(event *ev, double *xopt, double *fmin) nlopt_set_initial_step(opt, ss); fpars.ev = ev; - fpars.epsrel = 1e-1; iter = 0; - rv = nlopt_optimize(opt,x,&fval); + + /* First we do a set of "quick" minimizations to try and start the + * minimizer close to the minimum. To make these function evaluations + * faster, we set the absolute tolerance on the likelihood to 1.0, the + * maximum number of function evaluations to 100, and the relative + * tolerance on the numerical integration to 10%. */ + fpars.epsrel = 1e-1; + nlopt_set_ftol_abs(opt, 1.0); + nlopt_set_maxeval(opt, 100); + + for (i = 0; i < sizeof(startingParameters)/sizeof(startingParameters[0]); i++) { + x[0] = 1200.0; + x[1] = startingParameters[i].x; + x[2] = startingParameters[i].y; + x[3] = startingParameters[i].z; + x[4] = startingParameters[i].theta; + x[5] = startingParameters[i].phi; + x[6] = 130.0; + x[7] = 0.0; + x[8] = 0.0; + + rv = nlopt_optimize(opt,x,&fval); + + if (fval < *fmin || i == 0) { + *fmin = fval; + memcpy(xopt,x,sizeof(x)); + } + } + + memcpy(x,xopt,sizeof(x)); + + /* Now, we do the "real" minimization. */ + fpars.epsrel = 1e-4; + nlopt_set_ftol_abs(opt, 1e-5); + nlopt_set_maxeval(opt, 1000); + do { *fmin = fval; @@ -270,7 +324,7 @@ int main(int argc, char **argv) fit_event(&ev,xopt,&fmin); if (fout) { - fprintf(fout, "%" PRIu32 " %10.2f %7.2f %7.2f %7.2f %4.2f %4.2f %5.2f %5.2f %5.2f f() = %7.3e\n", + fprintf(fout, "%" PRIu32 " %10.2f %7.2f %7.2f %7.2f %5.2f %5.2f %5.2f %5.2f %5.2f f() = %7.3e\n", ev.gtid, xopt[0], xopt[1], |