diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile | 2 | ||||
| -rw-r--r-- | src/fit.c | 132 |
2 files changed, 52 insertions, 82 deletions
diff --git a/src/Makefile b/src/Makefile index 2359e66..ae8807b 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,5 +1,5 @@ CFLAGS=-O0 -Wall -g -DSWAP_BYTES -LDLIBS=-lm -lgsl -lgslcblas +LDLIBS=-lm -lgsl -lgslcblas -lnlopt_cxx -lstdc++ all: test test-vector test-likelihood fit test-charge @@ -4,39 +4,51 @@ #include "Record_Info.h" #include "event.h" #include "zdab_utils.h" -#include <gsl/gsl_vector.h> -#include <gsl/gsl_multimin.h> #include "scattering.h" #include "pmt.h" #include "sno_charge.h" #include "db.h" #include "dqxx.h" +#include <nlopt.h> +#include <math.h> /* for sin(), cos(), etc. */ #define EV_RECORD 0x45562020 // 'EV ' (as written to ZDAB file) -double nll(const gsl_vector *v, void *params) +double nll(unsigned int n, const double *x, double *grad, void *params) { double T, theta, phi, t0; double pos[3], dir[3]; + static size_t iter; + double fval; - T = gsl_vector_get(v,0); + T = x[0]; - /* FIXME: should use a minimization algorithm with bounds. */ - if (T < 0) return 1e9; + pos[0] = x[1]; + pos[1] = x[2]; + pos[2] = x[3]; - pos[0] = gsl_vector_get(v,1); - pos[1] = gsl_vector_get(v,2); - pos[2] = gsl_vector_get(v,3); - - theta = gsl_vector_get(v,4); - phi = gsl_vector_get(v,5); + theta = x[4]; + phi = x[5]; dir[0] = sin(theta)*cos(phi); dir[1] = sin(theta)*sin(phi); dir[2] = cos(theta); - t0 = gsl_vector_get(v,6); + t0 = x[6]; + + fval = nll_muon((event *) params, T, pos, dir, t0); - return nll_muon((event *) params, T, pos, dir, t0); + printf("%5zu %.2g %.2g %.2g %.2g %.2g %.2g %.2g f() = %7.3e\n", + iter++, + x[0], + x[1], + x[2], + x[3], + x[4], + x[5], + x[6], + fval); + + return fval; } int main(int argc, char **argv) @@ -50,6 +62,9 @@ int main(int argc, char **argv) event ev = {0}; int crate, card, channel; int id; + double x[7]; + double ss[7]; + double opt_f; if (argc < 2) { fprintf(stderr, "usage: fit [filename]\n"); @@ -72,23 +87,8 @@ int main(int argc, char **argv) ev.run = -1; - const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex2; - gsl_multimin_fminimizer *s = NULL; - gsl_vector *ss, *x; - gsl_multimin_function minex_func; - - size_t iter = 0; - int status; - double size; - - /* Starting point */ - x = gsl_vector_alloc(7); - - ss = gsl_vector_alloc(7); - - minex_func.n = 7; - minex_func.f = nll; - minex_func.params = &ev; + nlopt_opt opt = nlopt_create(NLOPT_LN_SBPLX, 7); + nlopt_set_min_objective(opt,nll,&ev); init_interpolation(); init_charge(); @@ -132,54 +132,25 @@ int main(int argc, char **argv) if (ev.run != -1) { if (ev.run != bev.run || ev.gtid != bev.gtr_id) { /* New event, so we need to fit the old event. */ - iter = 0; - - gsl_vector_set(x,0,500.0); - gsl_vector_set(x,1,100.0); - gsl_vector_set(x,2,100.0); - gsl_vector_set(x,3,100.0); - gsl_vector_set(x,4,1.57); - gsl_vector_set(x,5,0.0); - gsl_vector_set(x,6,120.0); - - gsl_vector_set(ss,0,10.0); - gsl_vector_set(ss,1,10.0); - gsl_vector_set(ss,2,10.0); - gsl_vector_set(ss,3,10.0); - gsl_vector_set(ss,4,0.1); - gsl_vector_set(ss,5,0.1); - gsl_vector_set(ss,6,1.0); - - s = gsl_multimin_fminimizer_alloc(T, 7); - gsl_multimin_fminimizer_set(s, &minex_func, x, ss); - - do { - iter++; - status = gsl_multimin_fminimizer_iterate(s); - - if (status) break; - - size = gsl_multimin_fminimizer_size(s); - status = gsl_multimin_test_size(size,1e-2); - - if (status == GSL_SUCCESS) { - printf("converged to minimum at\n"); - } - - printf("%5zu %.2g %.2g %.2g %.2g %.2g %.2g %.2g f() = %7.3e size = %.3f\n", - iter, - gsl_vector_get(s->x,0), - gsl_vector_get(s->x,1), - gsl_vector_get(s->x,2), - gsl_vector_get(s->x,3), - gsl_vector_get(s->x,4), - gsl_vector_get(s->x,5), - gsl_vector_get(s->x,6), - s->fval, - size); - } while (status == GSL_CONTINUE && iter < 1000); - - gsl_multimin_fminimizer_free(s); + x[0] = 500.0; + x[1] = 100.0; + x[2] = 100.0; + x[3] = 100.0; + x[4] = 1.57; + x[5] = 0.0; + x[6] = 120.0; + + ss[0] = 10.0; + ss[1] = 10.0; + ss[2] = 10.0; + ss[3] = 10.0; + ss[4] = 0.1; + ss[5] = 0.1; + ss[6] = 1.0; + + nlopt_set_initial_step(opt, ss); + + nlopt_optimize(opt,x,&opt_f); } } ev.run = bev.run; @@ -191,8 +162,7 @@ int main(int argc, char **argv) db_free(db); - gsl_vector_free(x); - gsl_vector_free(ss); + nlopt_destroy(opt); zebra_close(f); |