diff --git a/App/App.pro b/App/App.pro index cb11d808964c63c884cab7291e48e41fb758d260..e0ca2ee5acfa01092bb13a861038427c90702976 100644 --- a/App/App.pro +++ b/App/App.pro @@ -56,7 +56,6 @@ SOURCES += \ src/TestRoughness.cpp \ src/TestToyExperiment.cpp \ src/TreeEventStructure.cpp \ - src/ROOTMinimizerFunction.cpp \ src/ROOTGSLNLSMinimizer.cpp \ $${FUNCTIONAL_TESTS}/IsGISAXS01/IsGISAXS01.cpp diff --git a/App/inc/ROOTMinimizer.h b/App/inc/ROOTMinimizer.h index 3c7fe5c0492827d3b08b32f01a2f45ce0d90da5e..c610b715fe0780384cb61e612a93635318728c15 100644 --- a/App/inc/ROOTMinimizer.h +++ b/App/inc/ROOTMinimizer.h @@ -37,44 +37,33 @@ public: ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type=std::string()); virtual ~ROOTMinimizer(); - virtual void setParameter(size_t index, const FitParameter *par); - virtual void setParameters(const FitSuiteParameters ¶meters); + virtual void minimize(); + virtual void setParameter(size_t index, const FitParameter *par); - virtual void setFunction(function_chi2_t fun_chi2, size_t nparameters, function_gradient_t fun_gradient, size_t ndatasize); + virtual void setParameters(const FitSuiteParameters ¶meters); - virtual void minimize(); + virtual void setChiSquaredFunction(function_chi2_t fun_chi2, size_t nparameters); - //! return created minimizer - ROOT::Math::Minimizer *getROOTMinimizer() { return m_root_minimizer; } + virtual void setGradientFunction(function_gradient_t fun_gradient, size_t nparameters, size_t ndatasize); - //! get number of variables to fit virtual size_t getNumberOfVariables() const { return m_root_minimizer->NDim(); } - //! return minimum function value virtual double getMinValue() const { return m_root_minimizer->MinValue(); } - //! return value of variable corresponding the minimum of the function - virtual double getValueOfVariableAtMinimum(size_t i) const { - if(i >= getNumberOfVariables() ) throw OutOfBoundsException("ROOTMinimizer::getValueOfVariableAtMinimum() -> Wrong number of the variable"); - return m_root_minimizer->X()[i]; - } + virtual double getValueOfVariableAtMinimum(size_t i) const {return m_root_minimizer->X()[check_index(i)]; } - //! return value of variable corresponding the minimum of the function - virtual double getErrorOfVariable(size_t i) const { - if(i >= getNumberOfVariables() ) throw OutOfBoundsException("ROOTMinimizer::getErrorOfVariable() -> Wrong number of the variable"); - return (m_root_minimizer->Errors() == 0? 0 : m_root_minimizer->Errors()[i]); - } + virtual double getErrorOfVariable(size_t i) const { return (m_root_minimizer->Errors() == 0? 0 : m_root_minimizer->Errors()[check_index(i)]); } - //! printing results virtual void printResults() const; - //! clear resources (parameters) for consecutives minimizations virtual void clear() { m_root_minimizer->Clear(); } - //! printing minimizer description virtual void printOptions() const; + //! return created minimizer + ROOT::Math::Minimizer *getROOTMinimizer() { return m_root_minimizer; } + //! checking validity of the combination minimizer_name and algo_type bool isValidNames(const std::string &minimizer_name, const std::string &algo_type); @@ -82,20 +71,16 @@ public: bool isGradientBasedAgorithm(); private: + ROOTMinimizer(const ROOTMinimizer &); + ROOTMinimizer &operator=(const ROOTMinimizer &); + + size_t check_index(size_t index) const { return index<getNumberOfVariables() ? index : throw OutOfBoundsException("ROOTMinimizer::getErrorOfVariable() -> Wrong number of the variable"); } + std::string m_minimizer_name; std::string m_algo_type; - ROOT::Math::Minimizer *m_root_minimizer; - ROOTMinimizerFunction * m_minfunc; - ROOTMinimizerElementFunction * m_minfunc_element; - -// IMinimizerFunction *m_minimizer_function; -// ROOT::Math::Functor *m_fcn; //! function to minimize -// ROOT::Math::GradFunctor *m_fcn_grad; //! gradient of function to minimize -// function_t m_fcn; -// element_function_t m_element_fcn; -// int m_ndims; -// int m_nelements; + ROOTMinimizerChiSquaredFunction *m_chi2_func; + ROOTMinimizerGradientFunction *m_gradient_func; }; #endif // ROOTMINIMIZER_H diff --git a/App/inc/ROOTMinimizerFunction.h b/App/inc/ROOTMinimizerFunction.h index 6e242c2ced6d29819145ffcb8bbd406d874a663c..407c7bbc5847a974990380306b7f959ac6aeb99e 100644 --- a/App/inc/ROOTMinimizerFunction.h +++ b/App/inc/ROOTMinimizerFunction.h @@ -14,72 +14,61 @@ //! @author Scientific Computing Group at FRM II //! @date Dec 10, 2012 - -#include "Numeric.h" #include "IMinimizer.h" #include "Math/Functor.h" #include "Math/FitMethodFunction.h" -#include <vector> -#include <cmath> + //- ------------------------------------------------------------------- -//! @class ROOTMinimizerFunction -//! @brief Basic minimizer function +//! @class ROOTMinimizerChiSquaredFunction +//! @brief minimizer chi2 function //- ------------------------------------------------------------------- -class ROOTMinimizerFunction : public ROOT::Math::Functor +class ROOTMinimizerChiSquaredFunction : public ROOT::Math::Functor { public: - ROOTMinimizerFunction(IMinimizer::function_chi2_t fcn, int ndims ) : ROOT::Math::Functor(fcn, ndims), m_fcn(fcn) {} - virtual ~ROOTMinimizerFunction(){} + ROOTMinimizerChiSquaredFunction(IMinimizer::function_chi2_t fcn, int ndims ) : ROOT::Math::Functor(fcn, ndims), m_fcn(fcn) {} + virtual ~ROOTMinimizerChiSquaredFunction(){} IMinimizer::function_chi2_t m_fcn; }; //- ------------------------------------------------------------------- -//! @class ROOTMinimizerElementFunction +//! @class ROOTMinimizerGradientFunction //! @brief Minimizer function with access to single data element residuals. //! Required by Fumili, Fumili2 and GSLMultiMin minimizers //- ------------------------------------------------------------------- -class ROOTMinimizerElementFunction : public ROOT::Math::FitMethodFunction +class ROOTMinimizerGradientFunction : public ROOT::Math::FitMethodFunction { public: typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t Type_t; - ROOTMinimizerElementFunction(IMinimizer::function_gradient_t fun_gradient, size_t npars, size_t ndatasize) + ROOTMinimizerGradientFunction(IMinimizer::function_gradient_t fun_gradient, size_t npars, size_t ndatasize) : ROOT::Math::FitMethodFunction(npars, ndatasize) , m_fun_gradient(fun_gradient) { } - virtual ~ROOTMinimizerElementFunction(){} + virtual ~ROOTMinimizerGradientFunction(){} Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; } - ROOT::Math::IMultiGenFunction * Clone() const { return new ROOTMinimizerElementFunction(m_fun_gradient, NDim(), NPoints()); } + ROOT::Math::IMultiGenFunction * Clone() const { return new ROOTMinimizerGradientFunction(m_fun_gradient, NDim(), NPoints()); } //! evaluation of single data element residual - double DataElement(const double *pars, unsigned int i_data, double *gradient = 0) const - { + double DataElement(const double *pars, unsigned int i_data, double *gradient = 0) const { return m_fun_gradient(pars, i_data, gradient); } private: //! evaluation of chi2 - double DoEval(const double * pars) const - { - //std::cout << "ROOTMinimizerFuncion::DoEval() -> 1.1 ndim:" <<NDim() <<" npoint:" << NPoints() << std::endl; + double DoEval(const double * pars) const { double chi2 = 0.0; for(size_t i_data=0; i_data<NPoints(); ++i_data) { double res = DataElement(pars, i_data); chi2 += res*res; } - //std::cout << "ROOTMinimizerFuncion::DoEval() -> 1.2 chi:" <<chi2 << " ndim:" << NDim() << " np:" << NPoints() << std::endl; return chi2/double(NPoints()); - //return chi2; } - IMinimizer::function_gradient_t m_fun_gradient; }; - - #endif // ROOTMINIMIZERFUNCTION_H diff --git a/App/src/IsGISAXSData.cpp b/App/src/IsGISAXSData.cpp index 419aec4471854562d72b6f8cc2976152bb669fb0..d209a2f9c5f966c61cc84cdce01bfe9180a59c44 100644 --- a/App/src/IsGISAXSData.cpp +++ b/App/src/IsGISAXSData.cpp @@ -162,7 +162,6 @@ OutputData<double> *IsGISAXSData::convert_isgi_scan(std::vector<IsgiData > &isgi for(size_t i_point=0; i_point<isgi_data.size(); ++i_point) { phi_axis.push_back(isgi_data[i_point].phif); } - } OutputData<double > * data = new OutputData<double >; data->addAxis(phi_axis); diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp index 95ea8cc788eb4104174ecc923494ce0a65ae5012..1fe11f34fac822aace4a996da59c62ae8dc0453d 100644 --- a/App/src/IsGISAXSTools.cpp +++ b/App/src/IsGISAXSTools.cpp @@ -20,6 +20,7 @@ #include <iomanip> #include <cmath> #include <algorithm> +#include <cassert> double IsGISAXSTools::m_hist_min = 0; double IsGISAXSTools::m_hist_max = 0; @@ -34,6 +35,7 @@ void IsGISAXSTools::drawLogOutputData(const OutputData<double>& output, const std::string& canvas_name, const std::string& canvas_title, const std::string& draw_options, const std::string &histogram_title) { + assert(&output); TCanvas *c1 = new TCanvas(canvas_name.c_str(), canvas_title.c_str(), 0, 0, 1024, 768); IsGISAXSTools::setMinimum(1.); c1->cd(); @@ -48,6 +50,7 @@ void IsGISAXSTools::drawOutputData(const OutputData<double>& output, const std::string& canvas_name, const std::string& canvas_title, const std::string& draw_options, const std::string &histogram_title) { + assert(&output); TCanvas *c1 = new TCanvas(canvas_name.c_str(), canvas_title.c_str(), 0, 0, 1024, 768); c1->cd(); drawOutputDataInPad(output, draw_options, histogram_title); @@ -59,6 +62,7 @@ void IsGISAXSTools::drawOutputData(const OutputData<double>& output, void IsGISAXSTools::drawOutputDataInPad(const OutputData<double>& output, const std::string& draw_options, const std::string &histogram_title) { + assert(&output); if(!gPad) { throw NullPointerException("IsGISAXSTools::drawOutputDataInPad() -> Error! No canvas exists."); } @@ -85,6 +89,7 @@ void IsGISAXSTools::drawOutputDataInPad(const OutputData<double>& output, /* ************************************************************************* */ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const std::string &histo_name) { + assert(&output); if (output.getNdimensions() !=2) throw( "IsGISAXSTools::getOutputDataTH2D() -> Warning! Expected number of dimensiobs is 2."); std::vector<AxisStructure > haxises; @@ -146,6 +151,7 @@ TH2D *IsGISAXSTools::getOutputDataTH2D(const OutputData<double>& output, const s /* ************************************************************************* */ TH1 *IsGISAXSTools::getOutputDataTH123D(const OutputData<double>& output, const std::string &histo_name) { + assert(&output); if (output.getNdimensions() >3) { std::cout << "IsGISAXSTools::getOutputDataTH123D() -> Warning! Expected number of dimensions should be not more than 3" << std::endl; return 0; @@ -270,6 +276,8 @@ void IsGISAXSTools::drawOutputDataDistribution1D(const OutputData<double> &outpu /* ************************************************************************* */ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, const OutputData<double> &right, const std::string &draw_options, const std::string &histogram_title) { + assert(&left); + assert(&right); if(!gPad) { throw NullPointerException("IsGISAXSTools::drawOutputDataDifference1D() -> Error! No canvas exists."); } @@ -320,6 +328,8 @@ void IsGISAXSTools::drawOutputDataDifference1D(const OutputData<double> &left, c /* ************************************************************************* */ void IsGISAXSTools::drawOutputDataRelativeDifference2D(const OutputData<double> &left, const OutputData<double> &right, const std::string &draw_options, const std::string &histogram_title) { + assert(&left); + assert(&right); if(!gPad) { throw NullPointerException("IsGISAXSTools::drawOutputDataDifference2D -> Error! No canvas exists."); } @@ -342,6 +352,8 @@ void IsGISAXSTools::drawOutputDataRelativeDifference2D(const OutputData<double> /* ************************************************************************* */ void IsGISAXSTools::drawOutputDataChi2Difference2D(const OutputData<double> &left, const OutputData<double> &right, const std::string &draw_options, const std::string &histogram_title) { + assert(&left); + assert(&right); if(!gPad) { throw NullPointerException("IsGISAXSTools::drawOutputDataDifference2D -> Error! No canvas exists."); } @@ -465,6 +477,7 @@ void IsGISAXSTools::exportOutputDataInVectors2D(const OutputData<double> &output , std::vector<std::vector<double > > &v_axis0 , std::vector<std::vector<double > > &v_axis1) { + assert(&output_data); if (output_data.getRank() != 2) return; const IAxis *p_axis0 = output_data.getAxis(0); @@ -511,6 +524,8 @@ void IsGISAXSTools::exportOutputDataInVectors2D(const OutputData<double> &output /* ************************************************************************* */ TLine *IsGISAXSTools::getOutputDataScanLine(const OutputData<double> &data) { + assert(&data); + if(data.getNdimensions() != 2) throw LogicErrorException("IsGISAXSTools::getOutputDataScanLine() -> Error! Number of dimensions should be 2"); double x1(0), x2(0), y1(0), y2(0); if( data.getAxis(NDetector2d::ALPHA_AXIS_NAME) && data.getAxis(NDetector2d::ALPHA_AXIS_NAME)->getSize() == 1) { @@ -540,6 +555,8 @@ TLine *IsGISAXSTools::getOutputDataScanLine(const OutputData<double> &data) /* ************************************************************************* */ TH1D *IsGISAXSTools::getOutputDataScanHist(const OutputData<double> &data, const std::string &hname) { + assert(&data); + if(data.getNdimensions() != 2) throw LogicErrorException("IsGISAXSTools::getOutputDataScanHist() -> Error! Number of dimensions should be 2"); // one of axis should have dimension 1 if( (data.getAxis(NDetector2d::ALPHA_AXIS_NAME) && data.getAxis(NDetector2d::ALPHA_AXIS_NAME)->getSize() != 1) && (data.getAxis(NDetector2d::PHI_AXIS_NAME) && data.getAxis(NDetector2d::PHI_AXIS_NAME)->getSize() != 1)) @@ -578,6 +595,7 @@ TH1D *IsGISAXSTools::getOutputDataScanHist(const OutputData<double> &data, const /* ************************************************************************* */ OutputData<double > *IsGISAXSTools::createNoisyData(const OutputData<double> &exact_data, double noise_factor) { + assert(&exact_data); OutputData<double > *real_data = exact_data.clone(); OutputData<double>::iterator it = real_data->begin(); while (it != real_data->end()) { @@ -594,6 +612,7 @@ OutputData<double > *IsGISAXSTools::createNoisyData(const OutputData<double> &ex OutputData<double > *IsGISAXSTools::createDataWithGaussianNoise(const OutputData<double> &exact_data, double sigma) { + assert(&exact_data); OutputData<double > *real_data = exact_data.clone(); OutputData<double>::iterator it = real_data->begin(); while (it != real_data->end()) { @@ -608,6 +627,8 @@ OutputData<double > *IsGISAXSTools::createDataWithGaussianNoise(const OutputData void IsGISAXSTools::drawOutputDataComparisonResults(const OutputData<double> &data, const OutputData<double> &reference, const std::string &name, const std::string &title) { + assert(&data); + assert(&reference); TCanvas *c1 = DrawHelper::instance().createAndRegisterCanvas(name, title); c1->Divide(2,2); diff --git a/App/src/ROOTMinimizer.cpp b/App/src/ROOTMinimizer.cpp index b24af7286c200c88814368f042587f61e141cea6..45f7efc9bfe5a46a4390c80eaf5e00e2846c17b2 100644 --- a/App/src/ROOTMinimizer.cpp +++ b/App/src/ROOTMinimizer.cpp @@ -10,21 +10,20 @@ /* ************************************************************************* */ +// ROOTMinimizer c-tor // +// some usefull staff +// see http://root.cern.ch/phpBB3/viewtopic.php?f=15&t=14230&p=61216&hilit=minimizer+precision#p61216 +// see http://root.cern.ch/phpBB3/viewtopic.php?f=3&t=9181&hilit=precision+tolerance /* ************************************************************************* */ ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type) : m_minimizer_name(minimizer_name) , m_algo_type(algo_type) - , m_minfunc(0) - , m_minfunc_element(0) + , m_chi2_func(0) + , m_gradient_func(0) { - if( !isValidNames(m_minimizer_name, m_algo_type) ) throw LogicErrorException("ROOTMinimizer::ROOTMinimizer() -> Error! Wrong minimizer initialization parameters."); - // see http://root.cern.ch/phpBB3/viewtopic.php?f=15&t=14230&p=61216&hilit=minimizer+precision#p61216 - // see http://root.cern.ch/phpBB3/viewtopic.php?f=3&t=9181&hilit=precision+tolerance - //ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.1); - if( m_minimizer_name == "GSLMultiFit") { // hacked version of ROOT's GSL Levenberg-Marquardt minimizer m_root_minimizer = new ROOT::Patch::GSLNLSMinimizer(2); @@ -32,20 +31,14 @@ ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::strin m_root_minimizer = ROOT::Math::Factory::CreateMinimizer(minimizer_name, algo_type ); } if( !m_root_minimizer ) throw NullPointerException("ROOTMinimizer::ROOTMinimizer() -> Error! Can't create minimizer."); - -// m_root_minimizer->SetMaxFunctionCalls(20000); -// m_root_minimizer->SetMaxIterations(20000); -// m_root_minimizer->SetPrintLevel(4); -// m_root_minimizer->SetTolerance(0.01); -// m_root_minimizer->SetPrecision(1e-6); } ROOTMinimizer::~ROOTMinimizer() { delete m_root_minimizer; - delete m_minfunc; - delete m_minfunc_element; + delete m_chi2_func; + delete m_gradient_func; } @@ -149,21 +142,18 @@ void ROOTMinimizer::minimize() /* ************************************************************************* */ // set fcn function for minimizer /* ************************************************************************* */ -// FIXME ROOTMinimizer::setFunction Implement Multiple inheretiance in ROOTMinimizerElementFunction ;) -void ROOTMinimizer::setFunction(function_chi2_t fun_chi2, size_t nparameters, function_gradient_t fun_gradient, size_t ndatasize) +void ROOTMinimizer::setChiSquaredFunction(function_chi2_t fun_chi2, size_t nparameters) { - if( isGradientBasedAgorithm() ) { - std::cout << " ROOTMinimizer::setFunction() -> XXX 1.1 making ROOTMinimizerElementFunction " << std::endl; - delete m_minfunc_element; - m_minfunc_element = new ROOTMinimizerElementFunction(fun_gradient, nparameters, ndatasize); - m_root_minimizer->SetFunction(*m_minfunc_element); + delete m_chi2_func; + m_chi2_func = new ROOTMinimizerChiSquaredFunction(fun_chi2, nparameters); + if( !isGradientBasedAgorithm() ) m_root_minimizer->SetFunction(*m_chi2_func); +} - } else { - std::cout << " ROOTMinimizer::setFunction() -> XXX 1.2 making ROOTMinimizerFunction" << std::endl; - delete m_minfunc; - m_minfunc = new ROOTMinimizerFunction(fun_chi2, nparameters); - m_root_minimizer->SetFunction(*m_minfunc); - } +void ROOTMinimizer::setGradientFunction(function_gradient_t fun_gradient, size_t nparameters, size_t ndatasize) +{ + delete m_gradient_func; + m_gradient_func = new ROOTMinimizerGradientFunction(fun_gradient, nparameters, ndatasize); + if( isGradientBasedAgorithm() ) m_root_minimizer->SetFunction(*m_gradient_func); } @@ -204,8 +194,8 @@ void ROOTMinimizer::printResults() const std::cout << std::setw(25) << std::left << " Status " << ": " << m_root_minimizer->Status() << " '" << minimizerStatus[m_root_minimizer->Status()] << "'" << std::endl; std::cout << std::setw(25) << std::left << " IsValidError " << ": " << m_root_minimizer->IsValidError() << " '" << validErrorStatus[m_root_minimizer->Status()] << "'" <<std::endl; std::cout << std::setw(25) << std::left << " NCalls" << ": " << m_root_minimizer->NCalls() << std::endl; - if(m_minfunc_element) { - std::cout << std::setw(25) << std::left << " NCallsElement " << ": " << m_minfunc_element->NCalls() << std::endl; + if(m_gradient_func) { + std::cout << std::setw(25) << std::left << " NCallsElement " << ": " << m_gradient_func->NCalls() << std::endl; } std::cout << std::setw(25) << std::left << " MinValue" << ": " << std::scientific << std::setprecision(8) << getMinValue() << std::endl; std::cout << std::setw(25) << std::left << " Edm" << ": " << std::scientific << std::setprecision(8) << m_root_minimizer->Edm() << std::endl; diff --git a/App/src/ROOTMinimizerFunction.cpp b/App/src/ROOTMinimizerFunction.cpp deleted file mode 100644 index cc33dc7511b476522f427ec9b2d56adf6b6dc3be..0000000000000000000000000000000000000000 --- a/App/src/ROOTMinimizerFunction.cpp +++ /dev/null @@ -1,129 +0,0 @@ -#include "ROOTMinimizerFunction.h" - -//ROOTMinimizerFunction::ROOTMinimizerFunction() -//{ -//} - - -///* *************************************************************************** */ -//// evaluation of single data element residual -///* *************************************************************************** */ -//double ROOTMinimizerElementFunction::DataElement(const double *pars, unsigned int i_selected_element, double *gradient) const -//{ -// bool parameters_changed = parametersHaveChanged(pars); -// for(size_t i=0; i<m_ndims; ++i) { -// std::cout << pars[i] << " "; -// } -// std::cout << " parameters_changed:" << parameters_changed << " selected_element:" << i_selected_element << " index_prev:" << m_prev_element << std::endl; - -// int diff = i_selected_element-m_prev_element; -// if( diff >1 && ( diff != -1.*(int)m_nelements) ) std::cout << "XXXXX 1.1 Warning! (i-prev>1) Non sequential access to elements!" << std::endl; -// if( parameters_changed && i_selected_element!= 0 )std::cout << "XXXXX 1.2 Warning! Parameters changed and i!= 0" << std::endl; -// if( parameters_changed && i_selected_element == m_prev_element )std::cout << "XXXXX 1.2 Warning! Parameters changed and i=prev_element" << std::endl; -// m_prev_element = i_selected_element; - -// if(parameters_changed) { -// m_fcn(pars); -// for(size_t i_element=0; i_element<m_nelements; ++i_element) { -// m_residuals[i_element] = m_element_fcn(pars, i_element, gradient); -// } -// } - -// if( gradient ) { -// std::cout << " gradient " << std::endl; -// if(i_selected_element == 0 || parameters_changed) { -// std::cout << "g!=0 parameters_changed " << parameters_changed << std::endl; -// } -// for(size_t i_par=0; i_par<m_ndims; ++i_par) { -// gradient[i_par] = m_gradients[i_par][i_selected_element]; -// std::cout << " gradient_result i_par:" << i_par << " " << gradient[i_par] << std::endl; -// } -// } - -// // ./root/hist/hist/src/TF1.cxx - how root estimate derivative - -// m_ncalls_element++; -// return m_residuals[i_selected_element]; -//} - - - -//class ROOTMinimizerElementFunction : public ROOT::Math::FitMethodFunction -//{ -//public: -// typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t Type_t; - -// ROOTMinimizerElementFunction(IMinimizer::function_chi2_t fcn, size_t ndims, IMinimizer::function_gradient_t element_fcn, size_t nelements) -// : ROOT::Math::FitMethodFunction(ndims, nelements) -// , m_fcn(fcn) -// , m_element_fcn(element_fcn) -// , m_ndims(ndims) -// , m_nelements(nelements) -// , m_ncalls(0) -// , m_ncalls_element(0) -// , m_prev_element(0) -// { -// m_parameter_values.resize(ndims, 0.0); -// m_residuals.resize(nelements, 0.0); -// m_gradients.resize(ndims); -// for(size_t i_dim=0; i_dim<ndims; ++i_dim) { -// m_gradients[i_dim].resize(nelements, 0.0); -// } -// } - -// virtual ~ROOTMinimizerElementFunction(){} - -// Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; } -// ROOT::Math::IMultiGenFunction * Clone() const { return new ROOTMinimizerElementFunction(m_fcn, m_ndims, m_element_fcn, m_nelements); } - -// //! evaluation of chi2 -// double DoEval(const double * pars) const -// { -// double chi2 = 0.0; -// for(size_t i_data=0; i_data<m_nelements; ++i_data) { -// chi2 += DataElement(pars, i_data); -// } -// ++m_ncalls; -// std::cout << " ROOTMinimizerElementFunction::DoEval() -> " << chi2 << std::endl; -// return chi2; -// } - -// //! evaluation of single data element residual -// double DataElement(const double *pars, unsigned int i_selected_element, double *gradient = 0) const ; - -// virtual unsigned int NCalls() const { return m_ncalls; } -// virtual unsigned int NCallsElement() const { return m_ncalls_element; } - -// bool isEqual(double a, double b) const -// { -// return std::abs(a-b) < 2.*Numeric::double_epsilon; -// } - -// bool parametersHaveChanged(const double *pars) const -// { -// for(size_t i=0; i<m_ndims; ++i) { -// if( !isEqual(pars[i], m_parameter_values[i])) { -// std::copy(pars, pars+m_ndims, m_parameter_values.begin()); -// return true; -// } -// } -// return false; -// } - - -//private: - -// IMinimizer::function_chi2_t m_fcn; -// IMinimizer::function_gradient_t m_element_fcn; -// size_t m_ndims; -// size_t m_nelements; -// mutable std::vector<double > m_parameter_values; -// mutable size_t m_ncalls; -// mutable size_t m_ncalls_element; -// mutable size_t m_prev_element; -// mutable std::vector<double > m_residuals; // [m_nelements] -// mutable std::vector<std::vector<double> > m_gradients; // [m_ndims][m_nelements] -//}; - - - diff --git a/App/src/TestIsGISAXS12.cpp b/App/src/TestIsGISAXS12.cpp index da0bd1f8f250c31ea8e5f3ae586b7dc6cb96efd9..30c6edbfa931f4bf0ea75915ceaa9453ff9c6800 100644 --- a/App/src/TestIsGISAXS12.cpp +++ b/App/src/TestIsGISAXS12.cpp @@ -23,6 +23,7 @@ #include "ResolutionFunction2DSimple.h" #include "MathFunctions.h" #include "ROOTMinimizer.h" +#include "MinimizerTest.h" #include "OutputDataFunctions.h" #include "ExperimentConstants.h" #include "OutputDataIOFactory.h" @@ -366,7 +367,7 @@ void TestIsGISAXS12::run_test_minimizer() // our simulation produces numerically same results m_fitSuite = new FitSuite(); - m_fitSuite->setMinimizer( new TestMinimizer() ); + m_fitSuite->setMinimizer( new MinimizerTest() ); m_fitSuite->addFitParameter("*Normalizer/scale", 1.31159E+05, 100, AttLimits::limited(1e4, 2e5)); m_fitSuite->addFitParameter("*Normalizer/shift", -8.10009E-02, 1, AttLimits::limited(-10., 20.)); diff --git a/App/src/TestIsGISAXS13.cpp b/App/src/TestIsGISAXS13.cpp index 115ccd002b285e80dfe968d34756cfaf7b3e79af..7919f3c86bb81814559517ba5e76e99b83b99d39 100644 --- a/App/src/TestIsGISAXS13.cpp +++ b/App/src/TestIsGISAXS13.cpp @@ -28,6 +28,7 @@ #include "ExperimentConstants.h" #include "OutputDataIOFactory.h" #include "IsGISAXSData.h" +#include "MinimizerScan.h" #include <iostream> #include <fstream> @@ -80,20 +81,24 @@ void TestIsGISAXS13::run_isgisaxs_fit() //m_fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiFit", "") ); //m_fitSuite->setMinimizer( new ROOTMinimizer("Fumili", "") ); - m_fitSuite->setMinimizer( new ROOTMinimizer("GSLSimAn") ); - - m_fitSuite->attachObserver( new FitSuiteObserverPrint(10) ); - m_fitSuite->attachObserver( new FitSuiteObserverDraw(10) ); - ROOT::Math::Minimizer *minim = (dynamic_cast<ROOTMinimizer *>(m_fitSuite->getMinimizer()))->getROOTMinimizer(); - minim->SetPrintLevel(4); - minim->SetMaxIterations(400); - minim->SetMaxFunctionCalls(1000); - minim->SetTolerance(0.1); +// m_fitSuite->setMinimizer( new ROOTMinimizer("GSLSimAn") ); +// ROOT::Math::Minimizer *minim = (dynamic_cast<ROOTMinimizer *>(m_fitSuite->getMinimizer()))->getROOTMinimizer(); +// minim->SetPrintLevel(4); +// minim->SetMaxIterations(400); +// minim->SetMaxFunctionCalls(1000); +// minim->SetTolerance(0.1); //minim->SetPrecision(0.004); // minim->SetStrategy(1); // minim->SetPrecision(1.); + m_fitSuite->setMinimizer( new MinimizerScan(20) ); + + +// m_fitSuite->attachObserver( new FitSuiteObserverPrint(10) ); +// m_fitSuite->attachObserver( new FitSuiteObserverDraw(10) ); + + // Migrad // 0 *Normalizer/scale 1.333150e+05 7.316920e+04 6.781978e-01 // 1 *Normalizer/shift 1.212894e+00 1.469154e+01 2.140932e-01 @@ -125,9 +130,10 @@ void TestIsGISAXS13::run_isgisaxs_fit() // m_fitSuite->addFitParameter("*Normalizer/shift", 10, 10, AttLimits::limited(0., 20.)); m_fitSuite->addFitParameter("*Normalizer/scale", 1.333150e+05, 1e5, AttLimits::fixed()); m_fitSuite->addFitParameter("*Normalizer/shift", 1.212894e+00, 10, AttLimits::fixed()); - m_fitSuite->addFitParameter("*SampleBuilder/particle_radius", 4.0e+00*Units::nanometer, 0.01*Units::nanometer, AttLimits::limited(1., 10.) ); + m_fitSuite->addFitParameter("*SampleBuilder/particle_radius", 4.0e+00*Units::nanometer, 0.01*Units::nanometer, AttLimits::limited(2.0, 8.0) ); m_fitSuite->addFitParameter("*SampleBuilder/dispersion_radius", 1.991729e-01, 0.2, AttLimits::fixed() ); - m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio", 9.937262e-01, 1.0, AttLimits::fixed() ); + //m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio", 9.937262e-01, 1.0, AttLimits::fixed() ); + m_fitSuite->addFitParameter("*SampleBuilder/height_aspect_ratio", 9.937262e-01, 1.0, AttLimits::limited(0.5, 1.5) ); m_fitSuite->addFitParameter("*SampleBuilder/interf_distance", 1.497518e+01*Units::nanometer, 15.*Units::nanometer, AttLimits::fixed() ); m_fitSuite->addFitParameter("*SampleBuilder/interf_width", 3.023607e+00*Units::nanometer, 3.0*Units::nanometer, AttLimits::fixed() ); @@ -145,51 +151,73 @@ void TestIsGISAXS13::run_isgisaxs_fit() m_fitSuite->runFit(); - // drawing results - TCanvas *c2 = new TCanvas("c2","GISASFW fit results",800,600); - c2->Divide(2,2); - TLegend *leg1 = new TLegend(0.5,0.6,0.85,0.85); - leg1->SetBorderSize(1); - leg1->SetFillStyle(0); - for(size_t i_set=0; i_set<m_fitSuite->getFitObjects()->size(); ++i_set) { - c2->cd((int)i_set+1); - const FitObject *obj = m_fitSuite->getFitObjects()->getObject(i_set); - TH1D *hreal = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getRealData(),"gisasfw_real"); - TH1D *hsimul = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getSimulationData(),"gisasfw_simul"); - hreal->SetLineColor(kBlue); - gPad->SetLogy(); - hreal->DrawCopy(); - hsimul->DrawCopy("same"); - if(i_set==0) leg1->AddEntry(hreal,"GISASFW data","lp"); - if(i_set==0) leg1->AddEntry(hsimul,"GISASFW simul","lp"); + TCanvas *c1 = new TCanvas("c1","c1",1024,768); + c1->cd(); + + MinimizerScan *minim = dynamic_cast<MinimizerScan *>(m_fitSuite->getMinimizer()); + std::cout << "Results " << minim->getMinValue() << std::endl; + for(size_t i=0; i<minim->getNumberOfVariables(); ++i) { + std::cout << i << " " << minim->getValueOfVariableAtMinimum(i) << std::endl; } - c2->cd(1); leg1->Draw(); - c2->cd(2); leg1->Draw(); - - // drawing ratio - TLegend *leg2 = new TLegend(0.5,0.6,0.85,0.85); - leg2->SetBorderSize(1); - leg2->SetFillStyle(0); - for(size_t i_set=0; i_set<m_fitSuite->getFitObjects()->size(); ++i_set) { - c2->cd(3+1); - const FitObject *obj = m_fitSuite->getFitObjects()->getObject(i_set); - OutputData<double > *real = obj->getChiSquaredModule()->getRealData()->clone(); - OutputData<double > *simul = obj->getChiSquaredModule()->getSimulationData()->clone(); - - c2->cd((int)(i_set+3)); - *simul /= *real; - TH1D *hratio = IsGISAXSTools::getOutputDataScanHist(*simul,"gisasfw_real_simul_ratio"); - hratio->DrawCopy(); - if(i_set==0) { - leg2->AddEntry(hratio,"GISASFW simul/real","lp"); - } - delete real; - delete simul; + + const OutputData<double> *data = minim->getOutputData(); + if(data->getNdimensions()==1) { + TH1 * h1 = IsGISAXSTools::getOutputDataTH123D(*data,"hist"); + h1->Draw(); + gPad->SetLogy(); + } else { + TH2D *hist = IsGISAXSTools::getOutputDataTH2D( *data, "hist"); + hist->Draw("colz"); + gPad->SetLogz(); } - c2->cd(3); leg2->Draw(); - c2->cd(4); leg2->Draw(); - c2->Update(); + + +// // drawing results +// TCanvas *c2 = new TCanvas("c2","GISASFW fit results",800,600); +// c2->Divide(2,2); +// TLegend *leg1 = new TLegend(0.5,0.6,0.85,0.85); +// leg1->SetBorderSize(1); +// leg1->SetFillStyle(0); +// for(size_t i_set=0; i_set<m_fitSuite->getFitObjects()->size(); ++i_set) { +// c2->cd((int)i_set+1); +// const FitObject *obj = m_fitSuite->getFitObjects()->getObject(i_set); +// TH1D *hreal = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getRealData(),"gisasfw_real"); +// TH1D *hsimul = IsGISAXSTools::getOutputDataScanHist(*obj->getChiSquaredModule()->getSimulationData(),"gisasfw_simul"); +// hreal->SetLineColor(kBlue); +// gPad->SetLogy(); +// hreal->DrawCopy(); +// hsimul->DrawCopy("same"); +// if(i_set==0) leg1->AddEntry(hreal,"GISASFW data","lp"); +// if(i_set==0) leg1->AddEntry(hsimul,"GISASFW simul","lp"); +// } +// c2->cd(1); leg1->Draw(); +// c2->cd(2); leg1->Draw(); + +// // drawing ratio +// TLegend *leg2 = new TLegend(0.5,0.6,0.85,0.85); +// leg2->SetBorderSize(1); +// leg2->SetFillStyle(0); +// for(size_t i_set=0; i_set<m_fitSuite->getFitObjects()->size(); ++i_set) { +// c2->cd(3+1); +// const FitObject *obj = m_fitSuite->getFitObjects()->getObject(i_set); +// OutputData<double > *real = obj->getChiSquaredModule()->getRealData()->clone(); +// OutputData<double > *simul = obj->getChiSquaredModule()->getSimulationData()->clone(); + +// c2->cd((int)(i_set+3)); +// *simul /= *real; +// TH1D *hratio = IsGISAXSTools::getOutputDataScanHist(*simul,"gisasfw_real_simul_ratio"); +// hratio->DrawCopy(); +// if(i_set==0) { +// leg2->AddEntry(hratio,"GISASFW simul/real","lp"); +// } +// delete real; +// delete simul; +// } +// c2->cd(3); leg2->Draw(); +// c2->cd(4); leg2->Draw(); + +// c2->Update(); } diff --git a/Core/Core.pro b/Core/Core.pro index aa706e4fc4543516ce439c4b4d22ff3eaded24e5..bc1368efd575c5213f35bae9cdc82607c2e25295 100644 --- a/Core/Core.pro +++ b/Core/Core.pro @@ -51,6 +51,8 @@ SOURCES += \ Fitting/src/FitSuite.cpp \ Fitting/src/FitSuiteFunctions.cpp \ Fitting/src/FitSuiteParameters.cpp \ + Fitting/src/MinimizerTest.cpp \ + Fitting/src/MinimizerScan.cpp \ \ FormFactors/src/FormFactorBox.cpp \ FormFactors/src/FormFactorCrystal.cpp \ @@ -175,6 +177,8 @@ HEADERS += \ Fitting/inc/FitSuiteParameters.h \ Fitting/inc/FitSuiteStrategy.h \ Fitting/inc/IMinimizer.h \ + Fitting/inc/MinimizerTest.h \ + Fitting/inc/MinimizerScan.h \ \ FormFactors/inc/FormFactorBox.h \ FormFactors/inc/FormFactorCrystal.h \ diff --git a/Core/Fitting/inc/FitParameterLinked.h b/Core/Fitting/inc/FitParameterLinked.h index 2482ed781f031bc3c3257a379ae50c919c7f2849..cf3758db22965723ced5d9cf75a72651a2f171c1 100644 --- a/Core/Fitting/inc/FitParameterLinked.h +++ b/Core/Fitting/inc/FitParameterLinked.h @@ -30,18 +30,17 @@ class FitParameterLinked : public FitParameter { public: - typedef std::vector<ParameterPool::parameter_t > PoolParameterColl_t; + typedef std::vector<ParameterPool::parameter_t > pool_parameters_t; FitParameterLinked(); FitParameterLinked(const std::string &name, double value, double step, const AttLimits &attlim=AttLimits::limitless(), double error=0.0); virtual ~FitParameterLinked(){} //! set given value for all binded parameters - virtual void setValue(double value) - { + virtual void setValue(double value) { FitParameter::setValue(value); - for(PoolParameterColl_t::iterator it=m_parametercoll.begin(); it!=m_parametercoll.end(); ++it) { - (*it).setValue(value); // setting value for all registered parameters + for(pool_parameters_t::iterator it=m_pool_parameters.begin(); it!=m_pool_parameters.end(); ++it) { + (*it).setValue(value); } } @@ -55,14 +54,14 @@ public: friend std::ostream &operator<<(std::ostream &ostr, const FitParameterLinked &m) { m.print(ostr); return ostr; } protected: - // disabled copy constructor and assignment operator - FitParameterLinked(const FitParameterLinked &other); - FitParameterLinked &operator=(const FitParameterLinked &other); - //! print class void print(std::ostream &ostr) const; - PoolParameterColl_t m_parametercoll; //! collection of parameters from parameter pools + pool_parameters_t m_pool_parameters; //! collection of parameters from parameter pools + +private: + FitParameterLinked(const FitParameterLinked &); + FitParameterLinked &operator=(const FitParameterLinked &); }; #endif // FITMULTIPARAMETER_H diff --git a/Core/Fitting/inc/FitSuite.h b/Core/Fitting/inc/FitSuite.h index 09376b26397672134c198978e654e1804e67ca88..fc08eaa812238b8dd5ceaea083869e5c1a9df56e 100644 --- a/Core/Fitting/inc/FitSuite.h +++ b/Core/Fitting/inc/FitSuite.h @@ -70,12 +70,6 @@ public: //! run fitting which may consist of several minimization rounds virtual void runFit(); - //! function to minimize - //double fittingChiSquaredFunction(const double *pars_current_values); - - //! provides minimizer with gradients wrt parameters for single data element - //double fittingGradientFunction(const double *pars_current_values, unsigned int index, double *deriv); - //! return reference to the kit with data FitSuiteObjects *getFitObjects() { return &m_fit_objects; } diff --git a/Core/Fitting/inc/FitSuiteParameters.h b/Core/Fitting/inc/FitSuiteParameters.h index d1bb7407912361534ea1b3a6167462012ceb0bf8..7099e46de481888e31fb3efd2d099c44978e5f13 100644 --- a/Core/Fitting/inc/FitSuiteParameters.h +++ b/Core/Fitting/inc/FitSuiteParameters.h @@ -15,28 +15,30 @@ //! @date 15.11.2012 #include "Exceptions.h" +#include "SafePointerVector.h" #include "FitParameterLinked.h" #include <vector> //class Experiment; class ParameterPool; + //- ------------------------------------------------------------------- //! @class FitSuiteParameters //! @brief Class holds vector of parameters for FitSuite //- ------------------------------------------------------------------- -class FitSuiteParameters +class FitSuiteParameters : public SafePointerVector<FitParameter > { public: - typedef std::vector<FitParameter *> parameters_t; - typedef parameters_t::iterator iterator; - typedef parameters_t::const_iterator const_iterator; +// typedef std::vector<FitParameter *> parameters_t; +// typedef parameters_t::iterator iterator; +// typedef parameters_t::const_iterator const_iterator; - FitSuiteParameters(); - virtual ~FitSuiteParameters(); + FitSuiteParameters() {} + virtual ~FitSuiteParameters(){} - //! clear all defined parameters - void clear(); +// //! clear all defined parameters +// void clear(); //! add fit parameter void addParameter(const std::string &name, double value, double step, const AttLimits &attlim, double error=0.0); @@ -52,22 +54,22 @@ public: //! get values of all defined parameters std::vector<double > getValues() const; - //! return number of parameters - size_t size() const { return m_parameters.size(); } +// //! return number of parameters +// size_t size() const { return m_parameters.size(); } - //! return begin of container - iterator begin() { return m_parameters.begin(); } - const_iterator begin() const { return m_parameters.begin(); } +// //! return begin of container +// iterator begin() { return m_parameters.begin(); } +// const_iterator begin() const { return m_parameters.begin(); } - //! return end of container - iterator end() { return m_parameters.end(); } - const_iterator end() const { return m_parameters.end(); } +// //! return end of container +// iterator end() { return m_parameters.end(); } +// const_iterator end() const { return m_parameters.end(); } - //! access to parameters - const FitParameter *operator[](size_t index) const { return m_parameters[check_index(index)]; } - FitParameter *operator[](size_t index) { return m_parameters[check_index(index)]; } - const FitParameter *operator[](std::string name) const { return getParameter(name); } - FitParameter *operator[](std::string name) { return getParameter(name); } +// //! access to parameters +// const FitParameter *operator[](size_t index) const { return m_parameters[check_index(index)]; } +// FitParameter *operator[](size_t index) { return m_parameters[check_index(index)]; } +// const FitParameter *operator[](std::string name) const { return getParameter(name); } +// FitParameter *operator[](std::string name) { return getParameter(name); } //! linking fit parameters with pool parameters void link_to_pool(const ParameterPool *pool); @@ -80,13 +82,13 @@ public: private: //! disabled copy constructor and assignment operator - FitSuiteParameters &operator=(const FitSuiteParameters &other); - FitSuiteParameters(const FitSuiteParameters &other); +// FitSuiteParameters &operator=(const FitSuiteParameters &); +// FitSuiteParameters(const FitSuiteParameters &); - inline size_t check_index(size_t index) const { return (index < m_parameters.size() ? index : throw OutOfBoundsException("FitSuiteParameters::check_index() -> Index out of bounds") ); } - parameters_t m_parameters; //! collection of fit parameters +// inline size_t check_index(size_t index) const { return (index < m_parameters.size() ? index : throw OutOfBoundsException("FitSuiteParameters::check_index() -> Index out of bounds") ); } +// parameters_t m_parameters; //! collection of fit parameters - static double m_default_parameter_error; +// static double m_default_parameter_error; }; #endif // FITSUITEPARAMETERS_H diff --git a/Core/Fitting/inc/IMinimizer.h b/Core/Fitting/inc/IMinimizer.h index 885a8ff8d810c14fd16cb2a25d41f0b1032e2f5a..4b151e72ff312cc938131bdde3f526e87d042b8a 100644 --- a/Core/Fitting/inc/IMinimizer.h +++ b/Core/Fitting/inc/IMinimizer.h @@ -14,13 +14,12 @@ //! @author Scientific Computing Group at FRM II //! @date 05.10.2012 - -#include "FitParameter.h" -#include "FitSuiteParameters.h" #include <boost/function.hpp> -#include <map> #include "Exceptions.h" +class FitParameter; +class FitSuiteParameters; + //- ------------------------------------------------------------------- //! @class IMinimizer @@ -29,107 +28,108 @@ class IMinimizer { public: - //! signature of function to minimize + //! signature of chi squared function to minimize typedef boost::function<double(const double *)> function_chi2_t; - //! signature of function to minimize with acess to single element residual and gradient - typedef boost::function<double(const double *, unsigned int, double *)> function_gradient_t; - IMinimizer(){} - virtual ~IMinimizer(){} - - //! set parameter - virtual void setParameter(size_t index, const FitParameter *par) = 0; - virtual void setParameters(const FitSuiteParameters ¶meters) = 0; + //! signature of gradient to minimize with acess to single data element residuals + typedef boost::function<double(const double *, unsigned int, double *)> function_gradient_t; - //! set function to minimize - virtual void setFunction(function_chi2_t fun_chi2, size_t nparameters, function_gradient_t fun_gradient = function_gradient_t(), size_t ndatasize = 0) = 0; + IMinimizer() { } + virtual ~IMinimizer() { } //! run minimization virtual void minimize() = 0; + //! set internal minimizer parameter + virtual void setParameter(size_t index, const FitParameter *par); + + //! set internal minimizer parameters using external parameter list + virtual void setParameters(const FitSuiteParameters ¶meters); + + //! set chi squared function to minimize + virtual void setChiSquaredFunction(function_chi2_t fun_chi2, size_t nparameters); + + //! set gradient function to minimize + virtual void setGradientFunction(function_gradient_t fun_gradient, size_t nparameters, size_t ndatasize); + //! get number of variables to fit - virtual size_t getNumberOfVariables() const = 0; + virtual size_t getNumberOfVariables() const; //! return minimum function value - virtual double getMinValue() const = 0; + virtual double getMinValue() const; //! return pointer to the parameters values at the minimum - virtual double getValueOfVariableAtMinimum(size_t i) const = 0; + virtual double getValueOfVariableAtMinimum(size_t i) const; //! return pointer to the parameters values at the minimum - virtual double getErrorOfVariable(size_t i) const = 0; + virtual double getErrorOfVariable(size_t i) const; //! clear resources (parameters) for consecutives minimizations - virtual void clear() = 0; + virtual void clear(); //! print fit results - virtual void printResults() const = 0; - + virtual void printResults() const; }; - -//- ------------------------------------------------------------------- -//! @class TestMinimizer -//! @brief Minimizer which calls minimization function once to test whole chain -//- ------------------------------------------------------------------- -class TestMinimizer : public IMinimizer +inline void IMinimizer::setParameter(size_t index, const FitParameter *par) { -public: - TestMinimizer(){} - virtual ~TestMinimizer(){} - - //! set variable - virtual void setParameter(size_t index, const FitParameter *par) { m_values[index] = par->getValue(); } - virtual void setParameters(const FitSuiteParameters &/*parameters */) { throw NotImplementedException("TestMinimizer::setParameters() -> Error! Not implemented."); } - - //! set function to minimize - virtual void setFunction(function_chi2_t fun_chi2, size_t /* nparameters */, function_gradient_t /* fun_gradient */, size_t /* ndatasize */ ) { m_fcn = fun_chi2; } - - //! run minimization - virtual void minimize() - { - std::vector<double > buffer; - buffer.resize(m_values.size(), 0.0); - for(std::map<int, double >::iterator it=m_values.begin(); it!= m_values.end(); ++it ) { - buffer[it->first] = it->second; - std::cout << " minimize(): " << it->first << " " << it->second << std::endl; - } - std::cout << "TestMinimizer::minimize() -> Info. Calling fcn" << std::endl; - m_fcn(&buffer[0]); - } + (void)index; + (void)par; + throw NotImplementedException("IMinimizer::setParameter() -> Not implemented."); +} - //! get number of variables to fit - virtual size_t getNumberOfVariables() const { return m_values.size(); } +inline void IMinimizer::setParameters(const FitSuiteParameters ¶meters) +{ + (void)parameters; + throw NotImplementedException("IMinimizer::setParameters() -> Not implemented."); +} - //! return minimum function value - virtual double getMinValue() const { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); return 0.0; } +inline void IMinimizer::setChiSquaredFunction(function_chi2_t fun_chi2, size_t nparameters) +{ + (void)fun_chi2; + (void)nparameters; + throw NotImplementedException("IMinimizer::setChiSquaredFunction() -> Not implemented."); +} - //! return pointer to the parameters values at the minimum - virtual double getValueOfVariableAtMinimum(size_t i) const - { - std::map<int, double >::const_iterator pos = m_values.find((int)i); - if(pos != m_values.end()){ - return pos->second; - } else { - throw LogicErrorException("TestMinimizer::getValueOfVariableAtMinimum() -> Not found!"); - } - } +inline void IMinimizer::setGradientFunction(function_gradient_t fun_gradient, size_t nparameters, size_t ndatasize) +{ + (void)fun_gradient; + (void)nparameters; + (void)ndatasize; + throw NotImplementedException("IMinimizer::setGradientFunction() -> Not implemented."); +} - //! return pointer to the parameters values at the minimum - virtual double getErrorOfVariable(size_t /* i*/) const { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); return 0.0; } +inline size_t IMinimizer::getNumberOfVariables() const +{ + throw NotImplementedException("IMinimizer::getNumberOfVariables() -> Not implemented."); +} - //! clear resources (parameters) for consecutives minimizations - virtual void clear() { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); } +inline double IMinimizer::getMinValue() const +{ + throw NotImplementedException("IMinimizer::getMinValue() -> Not implemented."); +} - //! print fit results - virtual void printResults() const { throw NotImplementedException("TestMinimizer::getMinValue() -> Not implemented. "); } +inline double IMinimizer::getValueOfVariableAtMinimum(size_t index) const +{ + (void)index; + throw NotImplementedException("IMinimizer::getValueOfVariableAtMinimum() -> Not implemented."); +} -private: - std::map<int, double > m_values; - function_chi2_t m_fcn; -}; +inline double IMinimizer::getErrorOfVariable(size_t index) const +{ + (void)index; + throw NotImplementedException("IMinimizer::getErrorOfVariable() -> Not implemented."); +} +inline void IMinimizer::clear() +{ + throw NotImplementedException("IMinimizer::clear() -> Not implemented."); +} +inline void IMinimizer::printResults() const +{ + throw NotImplementedException("IMinimizer::printResults() -> Not implemented."); +} #endif // IMINIMIZER_H diff --git a/Core/Fitting/inc/MinimizerScan.h b/Core/Fitting/inc/MinimizerScan.h new file mode 100644 index 0000000000000000000000000000000000000000..dcc455c462d49d7394bd5b71b519c90ccc28560c --- /dev/null +++ b/Core/Fitting/inc/MinimizerScan.h @@ -0,0 +1,68 @@ +#ifndef MINIMIZERSCAN_H +#define MINIMIZERSCAN_H +// ******************************************************************** +// * The BornAgain project * +// * Simulation of neutron and x-ray scattering at grazing incidence * +// * * +// * LICENSE AND DISCLAIMER * +// * Lorem ipsum dolor sit amet, consectetur adipiscing elit. Mauris * +// * eget quam orci. Quisque porta varius dui, quis posuere nibh * +// * mollis quis. Mauris commodo rhoncus porttitor. * +// ******************************************************************** +//! @file MinimizerScan.h +//! @brief Definition of MinimizerScan class +//! @author Scientific Computing Group at FRM II +//! @date 21.01.2013 + + +#include "IMinimizer.h" +#include "FitSuiteParameters.h" +#include "SafePointerVector.h" +#include "OutputData.h" + +//- ------------------------------------------------------------------- +//! @class MinimizerScan +//! @brief Simple scan minimizer +//- ------------------------------------------------------------------- +class MinimizerScan : public IMinimizer +{ +public: + MinimizerScan(int nbins = 10) : m_parameter_map(0), m_nbins(nbins) { } + virtual ~MinimizerScan() { delete m_parameter_map; } + + virtual void minimize(); + + virtual void setParameters(const FitSuiteParameters ¶meters); + + virtual void setChiSquaredFunction(function_chi2_t fun_chi2, size_t nparameters); + + virtual void setGradientFunction(function_gradient_t fun_gradient, size_t nparameters, size_t ndatasize); + + virtual size_t getNumberOfVariables() const { return m_fit_parameters.size(); } + + virtual double getMinValue() const; + + virtual double getValueOfVariableAtMinimum(size_t i) const; + + virtual void printResults() const; + + + void setNbins(int nbins) { m_nbins = nbins; } + size_t getNbins() const { return m_nbins; } + + const OutputData<double > *getOutputData() { return m_parameter_map; } + +private: + + void construct_parameter_map(); + void set_parvalues_to_minimum(); + + OutputData<double > *m_parameter_map; + size_t m_nbins; //! number of bins per one parameter + //SafePointerVector<FitParameter > m_fit_parameters; + FitSuiteParameters m_fit_parameters; + std::vector<double> m_parvalues_at_minimum; + function_chi2_t m_fcn; +}; + +#endif // MINIMIZERSCAN_H diff --git a/Core/Fitting/inc/MinimizerTest.h b/Core/Fitting/inc/MinimizerTest.h new file mode 100644 index 0000000000000000000000000000000000000000..8bfbe325d9bc6f0c1bb90e95c3af58817a94f438 --- /dev/null +++ b/Core/Fitting/inc/MinimizerTest.h @@ -0,0 +1,53 @@ +#ifndef MINIMIZERTEST_H +#define MINIMIZERTEST_H +// ******************************************************************** +// * The BornAgain project * +// * Simulation of neutron and x-ray scattering at grazing incidence * +// * * +// * LICENSE AND DISCLAIMER * +// * Lorem ipsum dolor sit amet, consectetur adipiscing elit. Mauris * +// * eget quam orci. Quisque porta varius dui, quis posuere nibh * +// * mollis quis. Mauris commodo rhoncus porttitor. * +// ******************************************************************** +//! @file MinimizerTest.h +//! @brief Definition of MinimizerTest class +//! @author Scientific Computing Group at FRM II +//! @date 21.01.2013 + +#include "IMinimizer.h" +#include "FitSuiteParameters.h" +#include <map> + + +//- ------------------------------------------------------------------- +//! @class MinimizerTest +//! @brief Minimizer which calls minimization function once to test whole chain +//- ------------------------------------------------------------------- +class MinimizerTest : public IMinimizer +{ +public: + MinimizerTest(){} + virtual ~MinimizerTest(){} + + //! run minimization + virtual void minimize(); + + //! set variable + virtual void setParameter(size_t index, const FitParameter *par) { m_values[index] = par->getValue(); } + + //! set function to minimize + virtual void setChiSquaredFunction(function_chi2_t fun_chi2, size_t nparameters) { (void)nparameters, m_fcn = fun_chi2; } + + //! get number of variables to fit + virtual size_t getNumberOfVariables() const { return m_values.size(); } + + //! return pointer to the parameters values at the minimum + virtual double getValueOfVariableAtMinimum(size_t index) const; + +private: + std::map<int, double > m_values; + function_chi2_t m_fcn; +}; + +#endif // MINIMIZERTEST_H + diff --git a/Core/Fitting/src/FitParameterLinked.cpp b/Core/Fitting/src/FitParameterLinked.cpp index 27d4382e781949cf97b8d812f482aacae1cba88e..ad91792ea0b539be94db1d6a473b798a7cbe86f6 100644 --- a/Core/Fitting/src/FitParameterLinked.cpp +++ b/Core/Fitting/src/FitParameterLinked.cpp @@ -20,7 +20,7 @@ FitParameterLinked::FitParameterLinked(const std::string &name, double value, do void FitParameterLinked::addParameter(ParameterPool::parameter_t par) { if( !par.isNull() ) { - m_parametercoll.push_back(par); + m_pool_parameters.push_back(par); } else { throw LogicErrorException("FitMultiParameter::addParameter() -> Attempt to add null parameter"); } @@ -35,8 +35,8 @@ void FitParameterLinked::addMatchedParametersFromPool(const ParameterPool *pool, std::string wildcard_to_use = getName(); if( !wildcard.empty()) wildcard_to_use = wildcard; - PoolParameterColl_t matched_pars = pool->getMatchedParameters(wildcard_to_use); - m_parametercoll.insert(m_parametercoll.end(), matched_pars.begin(), matched_pars.end()); + pool_parameters_t matched_pars = pool->getMatchedParameters(wildcard_to_use); + m_pool_parameters.insert(m_pool_parameters.end(), matched_pars.begin(), matched_pars.end()); if( matched_pars.empty() ) { throw LogicErrorException("FitMultiParameter::addMatchedParametersFromPool() -> Error! Failed to add anything from pool using wildcard '"+wildcard_to_use+"'"); @@ -48,7 +48,7 @@ void FitParameterLinked::addMatchedParametersFromPool(const ParameterPool *pool, void FitParameterLinked::print(std::ostream &ostr) const { FitParameter::print(ostr); - ostr << "FitParameterLinked '" << getName() << "'" << " value:" << m_value << " collsize:" << m_parametercoll.size(); + ostr << "FitParameterLinked '" << getName() << "'" << " value:" << m_value << " collsize:" << m_pool_parameters.size(); // if(m_parametercoll.size() ) { // ostr << " addresses: "; // for(parametercoll_t::const_iterator it=m_parametercoll.begin(); it!=m_parametercoll.end(); it++) { diff --git a/Core/Fitting/src/FitSuite.cpp b/Core/Fitting/src/FitSuite.cpp index f5913009fd14b9f85028b3dbcd5d4881d399dddc..2caf51c5c107e36bf31f93b2c61898696dc93bb0 100644 --- a/Core/Fitting/src/FitSuite.cpp +++ b/Core/Fitting/src/FitSuite.cpp @@ -86,13 +86,11 @@ void FitSuite::link_fit_parameters() void FitSuite::minimize() { // initializing minimizer with fitting functions - //IMinimizer::function_chi2_t fun_chi2 = boost::bind(&FitSuite::fittingChiSquaredFunction, this, _1); - //IMinimizer::function_gradient_t fun_gradient = boost::bind(&FitSuite::fittingGradientFunction, this, _1, _2, _3); - IMinimizer::function_chi2_t fun_chi2 = boost::bind(&FitSuiteChiSquaredFunction::evaluate, &m_function_chi2, _1); - IMinimizer::function_gradient_t fun_gradient = boost::bind(&FitSuiteGradientFunction::evaluate, &m_function_gradient, _1, _2, _3); + m_minimizer->setChiSquaredFunction( fun_chi2, m_fit_parameters.size()); - m_minimizer->setFunction( fun_chi2, m_fit_parameters.size(), fun_gradient, m_fit_objects.getSizeOfDataSet() ); + IMinimizer::function_gradient_t fun_gradient = boost::bind(&FitSuiteGradientFunction::evaluate, &m_function_gradient, _1, _2, _3); + m_minimizer->setGradientFunction( fun_gradient, m_fit_parameters.size(), m_fit_objects.getSizeOfDataSet() ); // initializing minimizer's parameters with the list of local fit parameters m_minimizer->setParameters(m_fit_parameters); @@ -108,7 +106,7 @@ void FitSuite::minimize() bool FitSuite::check_prerequisites() const { if( !m_minimizer ) throw LogicErrorException("FitSuite::check_prerequisites() -> Error! No minimizer found."); - if( !m_fit_objects.size() ) throw LogicErrorException("FitSuite::check_prerequisites() -> Error! No experiment defined"); + if( !m_fit_objects.size() ) throw LogicErrorException("FitSuite::check_prerequisites() -> Error! No experiment/data description defined"); if( !m_fit_parameters.size() ) throw LogicErrorException("FitSuite::check_prerequisites() -> Error! No fit parameters defined"); return true; } @@ -149,38 +147,3 @@ void FitSuite::runFit() } -/* ************************************************************************* */ -// function to minimize -/* ************************************************************************* */ -//double FitSuite::fittingChiSquaredFunction(const double *pars_current_values) -//{ -// std::cout << "FitSuite::functionToMinimize() -> Info" << std::endl; -// // set fitting parameters to values suggested by the minimizer -// m_fit_parameters.setValues(pars_current_values); - -// // run simulations -// m_fit_objects.runSimulation(); - -// // caclulate chi2 value -// double chi_squared = m_fit_objects.getChiSquaredValue((int)m_fit_parameters.getNfreeParameters()); - -// notifyObservers(); -// m_n_call++; -// return chi_squared; -//} - - -/* ************************************************************************* */ -// provides minimizer with gradients wrt parameters for single data element -/* ************************************************************************* */ -//double FitSuite::fittingGradientFunction(const double *pars_current_values, unsigned int index, double *deriv) -//{ -// throw 1; -// (void)pars_current_values; -// (void) deriv; -// double residual = m_fit_objects.getResidualValue(index); -// return residual; -//} - - - diff --git a/Core/Fitting/src/FitSuiteParameters.cpp b/Core/Fitting/src/FitSuiteParameters.cpp index ee5a764be166d8136a76604650c38c71c3be5418..ef5e90d49d743bf76c903d5878481d17e1daca92 100644 --- a/Core/Fitting/src/FitSuiteParameters.cpp +++ b/Core/Fitting/src/FitSuiteParameters.cpp @@ -2,83 +2,136 @@ #include "Experiment.h" -double FitSuiteParameters::m_default_parameter_error=0.001; +//double FitSuiteParameters::m_default_parameter_error=0.001; -FitSuiteParameters::FitSuiteParameters() -{ -} +//FitSuiteParameters::FitSuiteParameters() +//{ +//} -FitSuiteParameters::~FitSuiteParameters() -{ - clear(); -} +//FitSuiteParameters::~FitSuiteParameters() +//{ +// clear(); +//} // clear all defined parameters -void FitSuiteParameters::clear() -{ - for(parameters_t::iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) { - delete (*it); - } - m_parameters.clear(); -} +//void FitSuiteParameters::clear() +//{ +//// for(parameters_t::iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) { +// for(parameters_t::iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) { +// delete (*it); +// } +// m_parameters.clear(); +//} // add fit parameter +//void FitSuiteParameters::addParameter(const std::string &name, double value, double step, const AttLimits &attlim, double error) +//{ +// for(parameters_t::const_iterator it = m_parameters.begin(); it!=m_parameters.end(); ++it) { +// if( (*it)->getName() == name ) throw LogicErrorException("FitSuiteParameters:addtFitParameter() -> Error. Existing parameter '"+name+"'"); +// } +// // defining default parameter error +//// if(error == 0.0) error = value*m_default_parameter_error; +// m_parameters.push_back(new FitParameterLinked(name, value, step, attlim, error) ); +//} + void FitSuiteParameters::addParameter(const std::string &name, double value, double step, const AttLimits &attlim, double error) { - for(parameters_t::const_iterator it = m_parameters.begin(); it!=m_parameters.end(); ++it) { + for(const_iterator it = begin(); it!=end(); ++it) { if( (*it)->getName() == name ) throw LogicErrorException("FitSuiteParameters:addtFitParameter() -> Error. Existing parameter '"+name+"'"); } // defining default parameter error - if(error == 0.0) error = value*m_default_parameter_error; - m_parameters.push_back(new FitParameterLinked(name, value, step, attlim, error) ); +// if(error == 0.0) error = value*m_default_parameter_error; + push_back(new FitParameterLinked(name, value, step, attlim, error) ); } + // return const fit parameter with given name +//const FitParameter *FitSuiteParameters::getParameter(const std::string &name) const +//{ +// for(parameters_t::const_iterator it = m_parameters.begin(); it!=m_parameters.end(); ++it) { +// if( (*it)->getName() == name ) return (*it); +// } +// throw LogicErrorException("FitSuiteParameters::getFitParameter() -> Error. No parameter with name '"+name+"'"); +//} const FitParameter *FitSuiteParameters::getParameter(const std::string &name) const { - for(parameters_t::const_iterator it = m_parameters.begin(); it!=m_parameters.end(); ++it) { + for(const_iterator it = begin(); it!=end(); ++it) { if( (*it)->getName() == name ) return (*it); } throw LogicErrorException("FitSuiteParameters::getFitParameter() -> Error. No parameter with name '"+name+"'"); } + +//// return fit parameter with given name +//FitParameter *FitSuiteParameters::getParameter(const std::string &name) +//{ +// for(parameters_t::iterator it = m_parameters.begin(); it!=m_parameters.end(); ++it) { +// if( (*it)->getName() == name ) return (*it); +// } +// throw LogicErrorException("FitSuiteParameters::getFitParameter() -> Error. No parameter with name '"+name+"'"); +//} // return fit parameter with given name FitParameter *FitSuiteParameters::getParameter(const std::string &name) { - for(parameters_t::iterator it = m_parameters.begin(); it!=m_parameters.end(); ++it) { + for(iterator it = begin(); it!=end(); ++it) { if( (*it)->getName() == name ) return (*it); } throw LogicErrorException("FitSuiteParameters::getFitParameter() -> Error. No parameter with name '"+name+"'"); } + // set values for all defined parameters // FIXME FitSuiteParameters::setValues remove check for attempt to set parameter values +//void FitSuiteParameters::setValues(const double *pars_values) +//{ +// if( !valuesAreDifferrent(pars_values) ) { +// std::cout << "FitSuiteParameters::setValues() -> Warning! Small or absent change in parameter values." << std::endl; +// for(size_t i=0; i<m_parameters.size(); i++) std::cout << (m_parameters[i]->getValue() -pars_values[i]) << " " << Numeric::areAlmostEqual(m_parameters[i]->getValue(), pars_values[i]) << std::endl; +// } +// size_t index(0); +// for(parameters_t::iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) (*it)->setValue(pars_values[index++]); +//} void FitSuiteParameters::setValues(const double *pars_values) { if( !valuesAreDifferrent(pars_values) ) { std::cout << "FitSuiteParameters::setValues() -> Warning! Small or absent change in parameter values." << std::endl; - for(size_t i=0; i<m_parameters.size(); i++) std::cout << (m_parameters[i]->getValue() -pars_values[i]) << " " << Numeric::areAlmostEqual(m_parameters[i]->getValue(), pars_values[i]) << std::endl; + for(size_t i=0; i<size(); i++) std::cout << ((*this)[i]->getValue() -pars_values[i]) << " " << Numeric::areAlmostEqual( (*this)[i]->getValue(), pars_values[i]) << std::endl; } size_t index(0); - for(parameters_t::iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) (*it)->setValue(pars_values[index++]); + for(iterator it=begin(); it!=end(); ++it) (*it)->setValue(pars_values[index++]); } - +//void FitSuiteParameters::setValues(const std::vector<double> &pars_values) +//{ +// if(pars_values.size() != m_parameters.size() ) throw OutOfBoundsException("FitSuiteParameters::setValues() -> Wrong size of array with parameter values"); +// setValues(&pars_values[0]); +//} void FitSuiteParameters::setValues(const std::vector<double> &pars_values) { - if(pars_values.size() != m_parameters.size() ) throw OutOfBoundsException("FitSuiteParameters::setValues() -> Wrong size of array with parameter values"); + if(pars_values.size() != size() ) throw OutOfBoundsException("FitSuiteParameters::setValues() -> Wrong size of array with parameter values"); setValues(&pars_values[0]); } + +//std::vector<double > FitSuiteParameters::getValues() const +//{ +// std::vector<double > result; +// for(parameters_t::const_iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) +// { +// result.push_back((*it)->getValue()); +// } +// return result; +//} + std::vector<double > FitSuiteParameters::getValues() const { std::vector<double > result; - for(parameters_t::const_iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) + for(const_iterator it=begin(); it!=end(); ++it) { result.push_back((*it)->getValue()); } @@ -86,10 +139,19 @@ std::vector<double > FitSuiteParameters::getValues() const } +//size_t FitSuiteParameters::getNfreeParameters() const +//{ +// size_t result(0); +// for(parameters_t::const_iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) +// { +// if( !(*it)->isFixed() ) result++; +// } +// return result; +//} size_t FitSuiteParameters::getNfreeParameters() const { size_t result(0); - for(parameters_t::const_iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) + for(const_iterator it=begin(); it!=end(); ++it) { if( !(*it)->isFixed() ) result++; } @@ -101,11 +163,21 @@ size_t FitSuiteParameters::getNfreeParameters() const /* ************************************************************************* */ // linking fit parameters with pool parameters /* ************************************************************************* */ +//void FitSuiteParameters::link_to_pool(const ParameterPool *pool) +//{ +// // linking fit parameter with whose pool parameters which match name of fit parameter +// // going through all fit parameters defined +// for(parameters_t::iterator it = m_parameters.begin(); it!= m_parameters.end(); ++it) { +// FitParameterLinked *par = dynamic_cast<FitParameterLinked *>((*it)); +// if( !par ) throw LogicErrorException("FitSuiteParameters::link_to_experiment() -> Error! Can't cast to FitParameterLinked."); +// par->addMatchedParametersFromPool(pool); +// } +//} void FitSuiteParameters::link_to_pool(const ParameterPool *pool) { // linking fit parameter with whose pool parameters which match name of fit parameter // going through all fit parameters defined - for(parameters_t::iterator it = m_parameters.begin(); it!= m_parameters.end(); ++it) { + for(iterator it = begin(); it!= end(); ++it) { FitParameterLinked *par = dynamic_cast<FitParameterLinked *>((*it)); if( !par ) throw LogicErrorException("FitSuiteParameters::link_to_experiment() -> Error! Can't cast to FitParameterLinked."); par->addMatchedParametersFromPool(pool); @@ -113,10 +185,18 @@ void FitSuiteParameters::link_to_pool(const ParameterPool *pool) } +//bool FitSuiteParameters::valuesAreDifferrent(const double *pars_values, double tolerance_factor) const +//{ +// size_t index(0); +// for(parameters_t::const_iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) { +// if( !Numeric::areAlmostEqual(pars_values[index++], (*it)->getValue(), tolerance_factor )) return true; +// } +// return false; +//} bool FitSuiteParameters::valuesAreDifferrent(const double *pars_values, double tolerance_factor) const { size_t index(0); - for(parameters_t::const_iterator it=m_parameters.begin(); it!=m_parameters.end(); ++it) { + for(const_iterator it=begin(); it!=end(); ++it) { if( !Numeric::areAlmostEqual(pars_values[index++], (*it)->getValue(), tolerance_factor )) return true; } return false; diff --git a/Core/Fitting/src/MinimizerScan.cpp b/Core/Fitting/src/MinimizerScan.cpp new file mode 100644 index 0000000000000000000000000000000000000000..7d79340b160342e1862be81402dbbc8f22213ae7 --- /dev/null +++ b/Core/Fitting/src/MinimizerScan.cpp @@ -0,0 +1,101 @@ +#include "MinimizerScan.h" +#include <algorithm> + + +/* ************************************************************************* */ +// +/* ************************************************************************* */ +void MinimizerScan::minimize() +{ + construct_parameter_map(); + + // scanning values of fit parameters + for(OutputData<double>::iterator it = m_parameter_map->begin(); it!=m_parameter_map->end(); ++it) { + for(size_t i_axis=0; i_axis<m_parameter_map->getNdimensions(); ++i_axis) { + size_t xbin = m_parameter_map->toCoordinate(it.getIndex(), i_axis); + double value = (*m_parameter_map->getAxis(i_axis))[xbin]; + std::string parname = m_parameter_map->getAxis(i_axis)->getName(); + m_fit_parameters.getParameter(parname)->setValue(value); + } + std::vector<double> current_values=m_fit_parameters.getValues(); + (*it) = m_fcn(¤t_values[0]); + } + + set_parvalues_to_minimum(); +} + + +/* ************************************************************************* */ +// Construct N dimensional space over all fit parameters with lower and upper limits +// defined. +/* ************************************************************************* */ +void MinimizerScan::construct_parameter_map() +{ + delete m_parameter_map; + m_parameter_map = new OutputData<double>; + for(size_t i_par=0; i_par < m_fit_parameters.size(); i_par++ ) { + const FitParameter *par = m_fit_parameters[i_par]; + if( par->hasLowerLimit() && par->hasUpperLimit() ) { + AxisDouble axis(par->getName(), m_nbins, par->getLowerLimit(), par->getUpperLimit()); + m_parameter_map->addAxis(axis); + } + } + m_parameter_map->setAllTo(0.0); +} + + +void MinimizerScan::set_parvalues_to_minimum() +{ + assert(m_parameter_map); + OutputData<double>::iterator it = std::min_element(m_parameter_map->begin(), m_parameter_map->end()); + for(size_t i_axis=0; i_axis<m_parameter_map->getNdimensions(); ++i_axis) { + size_t xbin = m_parameter_map->toCoordinate(it.getIndex(), i_axis); + double value = (*m_parameter_map->getAxis(i_axis))[xbin]; + std::string parname = m_parameter_map->getAxis(i_axis)->getName(); + m_fit_parameters.getParameter(parname)->setValue(value); + } + +} + + +double MinimizerScan::getMinValue() const +{ + assert(m_parameter_map); + return *std::min_element(m_parameter_map->begin(), m_parameter_map->end()); +} + + +void MinimizerScan::setGradientFunction(function_gradient_t fun_gradient, size_t nparameters, size_t ndatasize) +{ + (void) fun_gradient; + (void) nparameters; + (void) ndatasize; +} + + +double MinimizerScan::getValueOfVariableAtMinimum(size_t index) const +{ + return m_fit_parameters[index]->getValue(); +} + + +void MinimizerScan::printResults() const +{ + std::cout << "MinimizerScan::printResults() -> Info" << std::endl; +} + + +void MinimizerScan::setChiSquaredFunction(function_chi2_t fun_chi2, size_t nparameters) +{ + (void)nparameters; + m_fcn = fun_chi2; +} + + +void MinimizerScan::setParameters(const FitSuiteParameters ¶meters) +{ + m_fit_parameters.clear(); + for(size_t i_par = 0; i_par<parameters.size(); ++i_par) { + m_fit_parameters.push_back(new FitParameter( *parameters[i_par] ) ); + } +} diff --git a/Core/Fitting/src/MinimizerTest.cpp b/Core/Fitting/src/MinimizerTest.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8bd15a9b75c8383c682c0516bc0f72481fdf8e24 --- /dev/null +++ b/Core/Fitting/src/MinimizerTest.cpp @@ -0,0 +1,27 @@ +#include "MinimizerTest.h" + + +// run minimization +void MinimizerTest::minimize() +{ + std::vector<double > buffer; + buffer.resize(m_values.size(), 0.0); + for(std::map<int, double >::iterator it=m_values.begin(); it!= m_values.end(); ++it ) { + buffer[it->first] = it->second; + std::cout << " minimize(): " << it->first << " " << it->second << std::endl; + } + std::cout << "TestMinimizer::minimize() -> Info. Calling fcn" << std::endl; + m_fcn(&buffer[0]); +} + + +// return pointer to the parameters values at the minimum +double MinimizerTest::getValueOfVariableAtMinimum(size_t index) const +{ + std::map<int, double >::const_iterator pos = m_values.find((int)index); + if(pos != m_values.end()){ + return pos->second; + } else { + throw LogicErrorException("TestMinimizer::getValueOfVariableAtMinimum() -> Not found!"); + } +} diff --git a/Core/Tools/inc/SafePointerVector.h b/Core/Tools/inc/SafePointerVector.h index 6ae4143d39bacacc4db0be9ffa7d092393b7b8d3..7110a591892ae2cc8385337c9a87ced652be6f9b 100644 --- a/Core/Tools/inc/SafePointerVector.h +++ b/Core/Tools/inc/SafePointerVector.h @@ -30,7 +30,7 @@ public: typedef typename std::vector<T *>::const_iterator const_iterator; SafePointerVector(); SafePointerVector(const SafePointerVector &other); - ~SafePointerVector(); + virtual ~SafePointerVector(); SafePointerVector &operator=(const SafePointerVector &right); size_t size() const; diff --git a/Tests/UnitTests/TestCore/TestCore b/Tests/UnitTests/TestCore/TestCore index f4d604aaee1ed3bc5d972dd439bf5b150c90cda4..24bf8fa1f1d558f78c018fd07a67d30b5a67689b 100755 Binary files a/Tests/UnitTests/TestCore/TestCore and b/Tests/UnitTests/TestCore/TestCore differ