// BLMinimize.hh #define USE_GSL_MULTIMIN //#define USE_TMINUIT /* This source file is part of G4beamline, http://g4beamline.muonsinc.com Copyright (C) 2003,2004,2005,2006 by Tom Roberts, all rights reserved. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. http://www.gnu.org/copyleft/gpl.html */ #ifndef BLMINIMIZE_HH #define BLMINIMIZE_HH #include #include #include #include "globals.hh" #include "BLAssert.hh" /** class BLMinimizeFunction is the base class for a class defining a * function to minimize. **/ class BLMinimizeFunction { public: BLMinimizeFunction() { } virtual ~BLMinimizeFunction() { } /// operator() is the function to minimize. Note it is not const /// for flexibility in derived classes. /// Its argument is necessarily const. virtual double operator()(const std::vector &x) = 0; }; #ifdef USE_GSL_MULTIMIN #include /** class BLMinimize implements a general function minimizer using * the GSL multimin functions. * * Uses the the Simplex algorithm of Nelder and Mead. * * Compilation for simple programs: * g++ -I... test.cc -lgsl -lgslcblas -lm * * NOTE: the parameters to the minimizer must all have scales that are * comparable, or the minimizer will often converge to a wildly incorrect * value. This applies to track fitting, where positions have scales of * 1 (millimeters), but angles have scales of 0.001 (milliradians). * This class scales the parameters appropriately. it also permits the * function to have arguments that are held constant during the * minimization by setting their scale to 0.0. * * NOTE: only a single instance of BLMinimize can be active at a time * (i.e. inside the minimize() function). Any number of instances can * exist. This basically means your BLMinimizeFunction cannot itself * use a minimizer internally. **/ class BLMinimize { std::vector scale; std::vector arg; std::vector name; double val; double size; int iter; int print; double chisqMin; BLMinimizeFunction *func; gsl_multimin_function ff; gsl_multimin_fminimizer *m; private: // f() is the function the gsl routines call (they are in C) static double f(const gsl_vector *x, void *param) { BLMinimize *p = (BLMinimize*)param; return p->call(x); } // call() will call the user's function, scaling the argument vector. double call(const gsl_vector *x) { if(print >= 2) printf("BLMinimize call args:"); int j=0; for(unsigned i=0; i= 2) printf(" %.4g",arg[i]); } if(print >= 2) printf("\n"); double v = (*func)(arg); if(print >= 2) printf("BLMinimize call returns value=%.4g\n",v); return v; } public: /// Constructor. BLMinimize(int _print=2) : scale(), arg(), name() { val=GSL_NAN; size=GSL_NAN; iter=0; print=_print; chisqMin=GSL_NAN; func=0; ff.f=0; ff.n=0; ff.params=0; m=0; } virtual ~BLMinimize() { if(m) gsl_multimin_fminimizer_free(m); } /// setFunc() will set the function to minimize and the scale of its /// arguments. The number of non-zero entries in _scale[] determines /// the number of variables to minimize. names is a colon-separated /// list of parameter names. void setFunc(BLMinimizeFunction &f, const std::vector &_scale, const std::string &names="") { func = &f; scale.clear(); arg.clear(); name.clear(); std::vector n = splitString(names); int nx = 0; for(unsigned i=0; i<_scale.size(); ++i) { double s=_scale[i]; scale.push_back(s); arg.push_back(s); name.push_back(i &x, double tolerance=1e-3, int maxIter=100) { if(ff.n <= 0) return GSL_EINVAL; BLAssert(x.size()==scale.size()); gsl_vector *init = gsl_vector_alloc(ff.n); gsl_vector *istep = gsl_vector_alloc(ff.n); unsigned j=0; for(unsigned i=0; i= 1) printf("BLMinimize FAILED status=%d\n", status); gsl_vector_free(init); gsl_vector_free(istep); return status; } val = m->fval; size = gsl_multimin_fminimizer_size(m); if(print >= 1) printf("BLMinimize Iter %d: val=%.4g size=%.4g\n", iter,val,size); status = gsl_multimin_test_size(size,tolerance); if(status == GSL_SUCCESS || (!gsl_isnan(chisqMin) && val <= chisqMin)) { unsigned j=0; for(unsigned i=0; ix,j++) * scale[i]; } gsl_vector_free(init); gsl_vector_free(istep); if(print >= 1) printf("BLMinimize SUCCESS\n"); return 0; } } gsl_vector_free(init); gsl_vector_free(istep); if(print >= 1) printf("BLMinimize FAILED iter > %d\n",maxIter); return GSL_EMAXITER; } /// value() returns the function value from the previous minimize(). double value() const { return val; } /// getSize() returns the Simplex size from the previous minimize(). double getSize() { return size; } /// getIterations() returns the # iterations of the previous minimize(). int getIterations() { return iter; } private: std::vector splitString(const std::string &s) { std::string tmp(s); std::vector v; while(tmp.size() > 0) { unsigned i=tmp.find_first_of(":"); if(i > tmp.size()) i = tmp.size(); v.push_back(tmp.substr(0,i)); tmp.erase(0,i+1); } return v; } }; #endif // USE_GSL_MULTIMIN #ifdef USE_TMINUIT #include "TMinuit.h" /** class BLMinimize implements a general function minimizer using * the TMinuit class from Root. * * NOTE: only a single instance of BLMinimize can be active at a time * (i.e. inside the minimize() function). Any number of instances can * exist. This basically means your BLMinimizeFunction cannot itself * use a minimizer internally. **/ class BLMinimize { std::vector scale; std::vector name; std::vector arg; double val; int iter; int print; double chisqMin; BLMinimizeFunction *func; TMinuit *tm; static BLMinimize *current; private: static void fcn(Int_t &npar, Double_t *grad, Double_t &val, Double_t *x, Int_t iflag) { BLAssert(current != 0); val = current->call(npar,x); } double call(int npar, double *x) { // npar seems to be 1 too small assert(npar == (int)arg.size()); for(unsigned i=0; i &_scale, const std::string &names="") { func = &f; scale.clear(); name.clear(); arg.clear(); std::vector n = splitString(names); for(unsigned i=0; i<_scale.size(); ++i) { double s=_scale[i]; scale.push_back(s); name.push_back(iSetFCN(fcn); tm->SetPrintLevel(-1); tm->Command("SET PRINTOUT -1"); } /// setPrintDetail() sets the detail of printing per iteration: /// 0=no print, 1=print value and size only, 2=print value, size, and x. /// the default is 2. void setPrintDetail(int v=0) { print = v; } /// setChisqMin() sets an additional convergence criterium of chisq /// less than the given value. Can only be used when the function is /// positive definite, usually a chi-squared. Required for track- /// fitting of test tracks that have truly 0 chisq -- the normal /// convergence often fails for such tracks. void setChisqMin(double v) { chisqMin = v; } /// minimize() finds the minimum, returning status (0=success). /// x the initial and then the final argument values /// tolerance the convergence limit on the simplex size /// NOTE: this is after scaling! /// maxIter limits the number of iterations /// The initial step for x[i] is scale[i]. int minimize(std::vector &x, double tolerance=1e-3, int maxIter=100) { BLAssert(tm != 0); BLAssert(x.size() == scale.size()); for(unsigned i=0; iDefineParameter(i,name[i].c_str(),x[i],scale[i],0.0,0.0); } //@ What do I do with tolerance???? tm->SetMaxIterations(maxIter); current = this; int retval = tm->Migrad(); current = 0; for(unsigned i=0; iGetParameter(i,x[i],err); } val = (*func)(x); return retval; } /// value() returns the function value from the previous minimize(). double value() const { return val; } /// getSize() returns the Simplex size from the previous minimize(). double getSize() { return 0.0; } /// getIterations() returns the # iterations of the previous minimize(). int getIterations() { return iter; } private: std::vector splitString(const std::string &s) { std::string tmp(s); std::vector v; while(tmp.size() > 0) { unsigned i=tmp.find_first_of(":"); if(i > tmp.size()) i = tmp.size(); v.push_back(tmp.substr(0,i)); tmp.erase(0,i+1); } return v; } }; #endif // USE_TMINUIT #endif // BLMINIMIZE_HH