From f011e48326aff652ff889ca2ac061355ce620534 Mon Sep 17 00:00:00 2001 From: pospelov <pospelov@fz-juelich.de> Date: Tue, 11 Dec 2012 09:51:08 +0100 Subject: [PATCH] New TestToyExperiment to test whole fitting chain --- App/App.pro | 8 ++- App/inc/IsGISAXSTools.h | 1 + App/inc/ROOTMinimizerFunction.h | 79 ++++++++++++++++++++ App/inc/TestFumiliLMA.h | 1 - App/inc/TestToyExperiment.h | 83 +++++++++++++++++++++ App/src/FunctionalTestFactory.cpp | 3 + App/src/IsGISAXSTools.cpp | 13 ++++ App/src/ROOTMinimizer.cpp | 3 + App/src/ROOTMinimizerFunction.cpp | 5 ++ App/src/TestFumiliLMA.cpp | 26 ++----- App/src/TestToyExperiment.cpp | 115 ++++++++++++++++++++++++++++++ Core/Tools/inc/FitSuite.h | 4 -- Core/Tools/inc/IMinimizer.h | 4 -- Core/Tools/src/FitSuite.cpp | 11 --- 14 files changed, 315 insertions(+), 41 deletions(-) create mode 100644 App/inc/ROOTMinimizerFunction.h create mode 100644 App/inc/TestToyExperiment.h create mode 100644 App/src/ROOTMinimizerFunction.cpp create mode 100644 App/src/TestToyExperiment.cpp diff --git a/App/App.pro b/App/App.pro index caf94fb33fc..14afb166f9c 100644 --- a/App/App.pro +++ b/App/App.pro @@ -51,7 +51,9 @@ SOURCES += \ src/TestPerformance.cpp \ src/TestRootTree.cpp \ src/TestRoughness.cpp \ - src/TreeEventStructure.cpp + src/TestToyExperiment.cpp \ + src/TreeEventStructure.cpp \ + src/ROOTMinimizerFunction.cpp HEADERS += \ inc/App.h \ @@ -93,7 +95,9 @@ HEADERS += \ inc/TestPerformance.h \ inc/TestRootTree.h \ inc/TestRoughness.h \ - inc/TreeEventStructure.h + inc/TestToyExperiment.h \ + inc/TreeEventStructure.h \ + inc/ROOTMinimizerFunction.h INCLUDEPATH += ./inc ../Core/Algorithms/inc ../Core/FormFactors/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc ../Core/PythonAPI/inc DEPENDPATH += ./inc ../Core/Algorithms/inc ../Core/FormFactors/inc ../Core/Geometry/inc ../Core/Samples/inc ../Core/Tools/inc ../Core/PythonAPI/inc diff --git a/App/inc/IsGISAXSTools.h b/App/inc/IsGISAXSTools.h index 621dd4cdf9f..0b0318b4117 100644 --- a/App/inc/IsGISAXSTools.h +++ b/App/inc/IsGISAXSTools.h @@ -106,6 +106,7 @@ public: //! create noisy data static OutputData<double > *createNoisyData(const OutputData<double> &exact_data, double noise_factor = 0.1); + static OutputData<double > *createDataWithGaussianNoise(const OutputData<double> &exact_data, double sigma); private: static double m_hist_min; // minimum value of y-axis (for 1D histograms), or z-axis (for 2D histograms) diff --git a/App/inc/ROOTMinimizerFunction.h b/App/inc/ROOTMinimizerFunction.h new file mode 100644 index 00000000000..09489e3095f --- /dev/null +++ b/App/inc/ROOTMinimizerFunction.h @@ -0,0 +1,79 @@ +#ifndef ROOTMINIMIZERFUNCTION_H +#define ROOTMINIMIZERFUNCTION_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 ROOTMinimizerFunction.h +//! @brief Definition of ROOTMinimizerFunction +//! @author Scientific Computing Group at FRM II +//! @date Dec 10, 2012 + + + +//#include "Math/Minimizer.h" +//#include "Math/Factory.h" +//#include "Math/FitMethodFunction.h" + + +////- ------------------------------------------------------------------- +////! @class ROOTMinimizerFunction +////! @brief +////- ------------------------------------------------------------------- + +//class ROOTMinimizerFunction : public ROOT::Math::FitMethodFunction +//{ +//public: +// typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t Type_t; +// typedef boost::function<double(const double *)> function_eval_t; +// typedef boost::function<double(const double *, unsigned int, double *)> function_element_t; + +// ROOTMinimizerFunction(); +// virtual ~ROOTMinimizerFunction(){} + +// Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; } +// ROOT::Math::IMultiGenFunction * Clone() const { return new MyChi2Function(m_test); } + +// double DoEval(const double * par) const { + + +//}; + + +//class MyChi2Function : public ROOT::Math::FitMethodFunction +//{ +//public: +// typedef ROOT::Math::BasicFitMethodFunction<ROOT::Math::IMultiGenFunction>::Type_t Type_t; + +// MyChi2Function(TestFumiliLMA *test) : ROOT::Math::FitMethodFunction(test->m_ndim, test->m_real_data->getAllocatedSize()), m_test(test) {} +// virtual ~MyChi2Function(){} + +// Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; } +// ROOT::Math::IMultiGenFunction * Clone() const { return new MyChi2Function(m_test); } + +// // evaluation of the all chi2 +// double DoEval(const double * par) const { +// int ndata = NPoints(); +// std::cout << "DoEval: " << ndata << std::endl; +// double chi2 = 0; +// for (int i = 0; i < ndata; ++i) { +// double res = DataElement( par, i); +// chi2 += res*res; +// } +// //std::cout << "DoEval: chi" << chi2/double(ndata) << std::endl; +// return chi2/double(ndata); +// } + +// double DataElement(const double *par, unsigned int i, double *g = 0) const; + +// TestFumiliLMA *m_test; +//}; + + + +#endif // ROOTMINIMIZERFUNCTION_H diff --git a/App/inc/TestFumiliLMA.h b/App/inc/TestFumiliLMA.h index 4d9b4228bd5..acc06904d94 100644 --- a/App/inc/TestFumiliLMA.h +++ b/App/inc/TestFumiliLMA.h @@ -47,7 +47,6 @@ public: double functionToMinimize(const double *pars); private: - OutputData<double > *createDataWithGaussianNoise(const OutputData<double> &exact_data, double sigma); void FillOutputDataFromFunction(OutputData<double> &data, TF2 *fun, int nbinsx=100, int nbinsy=100); ROOT::Math::Minimizer *m_root_minimizer; //! minimizer diff --git a/App/inc/TestToyExperiment.h b/App/inc/TestToyExperiment.h new file mode 100644 index 00000000000..2c7a4084c49 --- /dev/null +++ b/App/inc/TestToyExperiment.h @@ -0,0 +1,83 @@ +#ifndef TESTTOYEXPERIMENT_H +#define TESTTOYEXPERIMENT_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 TestToyExperiment.h +//! @brief Tests of fitting chain using toy experiment +//! @author Scientific Computing Group at FRM II +//! @date 10.12.2012 + +#include "IFunctionalTest.h" +#include "Experiment.h" +#include "OutputData.h" +#include "TestFumiliLMA.h" +#include "FitSuite.h" + +#include <vector> +#include "TF2.h" + + +//- ------------------------------------------------------------------- +//! @class ToyExperiment +//! @brief Toy experiment to test whole fitting chain with simple 2D functions +//- ------------------------------------------------------------------- +class ToyExperiment : public Experiment +{ +public: + ToyExperiment(TF2 *func) : m_func(func) { pars.resize(func->GetNpar(), 0.0); setName("ToyExperiment"); init_parameters(); } + virtual ~ToyExperiment() {} + virtual void runSimulation(); + virtual ToyExperiment *clone() const { return new ToyExperiment(*this); } + void setParameter(size_t i, double value) { pars[i] = value; } +protected: + virtual void init_parameters(); +private: + ToyExperiment(const ToyExperiment &other) : Experiment(other) + { + m_func = other.m_func; + pars = other.pars; + setName("ToyExperiment"); + init_parameters(); + } + TF2 *m_func; + std::vector<double > pars; +}; + + + + +//- ------------------------------------------------------------------- +//! @class TestToyExperiment +//! @brief Test of fitting chain using toy experiment +//- ------------------------------------------------------------------- +class TestToyExperiment : public IFunctionalTest +{ +public: + TestToyExperiment(); + virtual ~TestToyExperiment(); + virtual void execute(); + +private: + void initializeExperimentAndRealData(); + + IFunctionObject *m_func_object; //! simulation function + TF2 *m_func; //! ROOT representation of the simulation function with min, max defined + double m_sigma_noise; + ToyExperiment *m_experiment; + OutputData<double > *m_real_data; + FitSuite *m_fitSuite; + +}; + + + + + +#endif // TESTTOYEXPERIMENT_H diff --git a/App/src/FunctionalTestFactory.cpp b/App/src/FunctionalTestFactory.cpp index f87a1d06a54..dedf8e32104 100644 --- a/App/src/FunctionalTestFactory.cpp +++ b/App/src/FunctionalTestFactory.cpp @@ -26,6 +26,7 @@ #include "TestMiscellaneous.h" #include "TestFittingBenchmark.h" #include "TestFourier.h" +#include "TestToyExperiment.h" #include "TBenchmark.h" @@ -90,6 +91,8 @@ FunctionalTestFactory::FunctionalTestFactory() : m_benchmark(0) "functional test: test of Fourier transformation of OutputData maps"); registerItem("fumili", IFactoryCreateFunction<TestFumiliLMA, IFunctionalTest>, "functional test: test of ROOT's LMA-based minimizers Fumili and GSLMultiFit"); + registerItem("toyexp", IFactoryCreateFunction<TestToyExperiment, IFunctionalTest>, + "functional test: test fitting algorithms with toy experiment"); m_benchmark = new TBenchmark(); } diff --git a/App/src/IsGISAXSTools.cpp b/App/src/IsGISAXSTools.cpp index ad21caa47d9..606db559ab9 100644 --- a/App/src/IsGISAXSTools.cpp +++ b/App/src/IsGISAXSTools.cpp @@ -590,3 +590,16 @@ OutputData<double > *IsGISAXSTools::createNoisyData(const OutputData<double> &ex } +OutputData<double > *IsGISAXSTools::createDataWithGaussianNoise(const OutputData<double> &exact_data, double sigma) +{ + OutputData<double > *real_data = exact_data.clone(); + OutputData<double>::iterator it = real_data->begin(); + while (it != real_data->end()) { + double current = *it; + double random = MathFunctions::GenerateNormalRandom(0.0, sigma); + *it = current+random; + ++it; + } + return real_data; +} + diff --git a/App/src/ROOTMinimizer.cpp b/App/src/ROOTMinimizer.cpp index 7fa84a78677..5ae2f91574e 100644 --- a/App/src/ROOTMinimizer.cpp +++ b/App/src/ROOTMinimizer.cpp @@ -6,6 +6,9 @@ ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type) : m_fcn(0) { + if(minimizer_name == "Minuit2" && algo_type == "Fumili2") { + throw LogicErrorException("ROOTMinimizer::ROOTMinimizer() -> Error! Use word Fumili instead of 'Fumili2'"); + } m_root_minimizer = ROOT::Math::Factory::CreateMinimizer(minimizer_name.c_str(), algo_type.c_str() ); m_root_minimizer->SetMaxFunctionCalls(20000); m_root_minimizer->SetMaxIterations(20000); diff --git a/App/src/ROOTMinimizerFunction.cpp b/App/src/ROOTMinimizerFunction.cpp new file mode 100644 index 00000000000..5bad80951b4 --- /dev/null +++ b/App/src/ROOTMinimizerFunction.cpp @@ -0,0 +1,5 @@ +#include "ROOTMinimizerFunction.h" + +//ROOTMinimizerFunction::ROOTMinimizerFunction() +//{ +//} diff --git a/App/src/TestFumiliLMA.cpp b/App/src/TestFumiliLMA.cpp index b379321303c..cf6fc7be325 100644 --- a/App/src/TestFumiliLMA.cpp +++ b/App/src/TestFumiliLMA.cpp @@ -46,7 +46,7 @@ TestFumiliLMA::TestFumiliLMA() // "real" data OutputData<double> data; FillOutputDataFromFunction(data, m_func); - m_real_data = createDataWithGaussianNoise(data, m_sigma); + m_real_data = IsGISAXSTools::createDataWithGaussianNoise(data, m_sigma); } @@ -72,9 +72,9 @@ void TestFumiliLMA::execute() //m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad" ); //m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Fumili" ); // same as ("Fumili2") //m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("Fumili2"); // same as ("Minuit2", "Fumili" ) - //m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Fumili2" ); // same as minuit2/migrad - m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("Fumili"); // - //m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("GSLMultiFit"); + //m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Fumili2" ); // same as ("Minuit2", "Migrad" ), i.e. Fumili2 is wrong key + //m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("Fumili"); // + m_root_minimizer = ROOT::Math::Factory::CreateMinimizer("GSLMultiFit"); m_root_minimizer->SetVariable(0, "p0", 2.0, 0.01); m_root_minimizer->SetVariable(1, "p1", 0.0, 0.01); m_root_minimizer->SetVariable(2, "p2", 0.0, 0.01); @@ -157,20 +157,6 @@ double TestFumiliLMA::functionToMinimize(const double *pars) /* ************************************************************************* */ // /* ************************************************************************* */ -OutputData<double > *TestFumiliLMA::createDataWithGaussianNoise(const OutputData<double> &exact_data, double sigma) -{ - OutputData<double > *real_data = exact_data.clone(); - OutputData<double>::iterator it = real_data->begin(); - while (it != real_data->end()) { - double current = *it; - double random = MathFunctions::GenerateNormalRandom(0.0, sigma); - *it = current+random; - ++it; - } - return real_data; -} - - void TestFumiliLMA::FillOutputDataFromFunction(OutputData<double> &data, TF2 *fun, int nbinsx, int nbinsy) { data.clear(); @@ -197,7 +183,9 @@ void TestFumiliLMA::FillOutputDataFromFunction(OutputData<double> &data, TF2 *fu double MyChi2Function::DataElement(const double *pars, unsigned int i, double *g ) const { double xx[2]; - std::cout << " DataElement() -> " << g << std::endl; + std::cout << " DataElement() -> " << g << " " << i; + for(size_t i=0; i<m_test->m_func->GetNpar(); ++i) std::cout << " " << pars[i]; + std::cout << std::endl; const AxisDouble *xaxis = m_test->m_real_data->getAxis(0); const AxisDouble *yaxis = m_test->m_real_data->getAxis(1); size_t index_x = m_test->m_real_data->toCoordinates(i)[0]; diff --git a/App/src/TestToyExperiment.cpp b/App/src/TestToyExperiment.cpp new file mode 100644 index 00000000000..f32e2afa090 --- /dev/null +++ b/App/src/TestToyExperiment.cpp @@ -0,0 +1,115 @@ +#include "TestToyExperiment.h" +#include "Exceptions.h" +#include "IsGISAXSTools.h" +#include "FitSuite.h" +#include "FitSuiteHelper.h" +#include "ROOTMinimizer.h" +#include <iostream> + + +/* ************************************************************************* */ +// +/* ************************************************************************* */ +void ToyExperiment::runSimulation() +{ + if( !m_func ) throw NullPointerException("ToyExperiment::runSimulation() -> Error! No function is defined."); + + m_func->SetParameters(&pars[0]); + m_intensity_map.setAllTo(0.0); + OutputData<double>::iterator it = m_intensity_map.begin(); + while( it!= m_intensity_map.end() ) { + double phi_f = m_intensity_map.getValueOfAxis("phi_f", it.getIndex()); + double alpha_f = m_intensity_map.getValueOfAxis("alpha_f", it.getIndex()); + double value = m_func->Eval(phi_f, alpha_f); + (*it) = value; + ++it; + } +} + +void ToyExperiment::init_parameters() +{ + getParameterPool()->clear(); + for(size_t i=0; i<pars.size(); ++i) { + std::ostringstream ostr; + ostr << "par"<< i; + getParameterPool()->registerParameter(ostr.str(), &pars[i]); + } +} + + +/* ************************************************************************* */ +// +/* ************************************************************************* */ +TestToyExperiment::TestToyExperiment() + : m_func_object(0) + , m_func(0) + , m_sigma_noise(0) + , m_experiment(0) + , m_real_data(0) +{ + m_sigma_noise = 0.01; + m_func_object = new SincXSincYFunctionObject(); + m_func = new TF2("sincxy", m_func_object, -5.,5., -5.,5., 3, "SincXSincYFunctionObject"); + //m_func->SetParameters(1.0, 1.0, 0.5); // parameters we have to find +} + + +TestToyExperiment::~TestToyExperiment() +{ + delete m_func_object; + delete m_func; + delete m_experiment; + delete m_real_data; + delete m_fitSuite; +} + + +void TestToyExperiment::execute() +{ + std::cout << "TestToyExperiment()::execute() -> Hello World!" << std::endl; + + initializeExperimentAndRealData(); + + // setting up fitSuite + FitSuite *m_fitSuite = new FitSuite(); + m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") ); + m_fitSuite->attachObserver( new FitSuiteObserverPrint() ); + m_fitSuite->attachObserver( new FitSuiteObserverDraw() ); + + m_fitSuite->addFitParameter("*/par0", 1.0, 0.01); + m_fitSuite->addFitParameter("*/par1", 0.0, 0.01); + m_fitSuite->addFitParameter("*/par2", 0.0, 0.01); + + ChiSquaredModule chi_module; + chi_module.setChiSquaredFunction(SquaredFunctionWithGaussianError(m_sigma_noise) ); + m_fitSuite->addExperimentAndRealData(*m_experiment, *m_real_data, chi_module); +// m_fitSuite->addExperimentAndRealData(*m_experiment, *m_real_data); + m_fitSuite->runFit(); + +} + + + +/* ************************************************************************* */ +// +/* ************************************************************************* */ +void TestToyExperiment::initializeExperimentAndRealData() +{ + delete m_experiment; + m_experiment = new ToyExperiment(m_func); + + OutputData<double > tmp; + tmp.addAxis("phi_f", 100, m_func->GetXmin(), m_func->GetXmax()); + tmp.addAxis("alpha_f", 100, m_func->GetYmin(), m_func->GetYmax()); + m_experiment->setDetectorParameters(tmp); + + // generating real data + delete m_real_data; + m_experiment->setParameter(0, 2.0); + m_experiment->setParameter(1, 1.0); + m_experiment->setParameter(2, 0.5); + m_experiment->runSimulation(); + m_real_data = IsGISAXSTools::createDataWithGaussianNoise(*m_experiment->getOutputData(), m_sigma_noise); +} + + diff --git a/Core/Tools/inc/FitSuite.h b/Core/Tools/inc/FitSuite.h index a954400df1a..49833193e84 100644 --- a/Core/Tools/inc/FitSuite.h +++ b/Core/Tools/inc/FitSuite.h @@ -72,10 +72,6 @@ public: //! function to minimize double functionToMinimize(const double *pars_current_values); -// //! function gradient -// double functionGradient(const double *pars, int icoord); - - //! return reference to the kit with data FitSuiteObjects *getFitObjects() { return &m_fit_objects; } diff --git a/Core/Tools/inc/IMinimizer.h b/Core/Tools/inc/IMinimizer.h index bdef292947e..f22fd3b63ae 100644 --- a/Core/Tools/inc/IMinimizer.h +++ b/Core/Tools/inc/IMinimizer.h @@ -37,9 +37,6 @@ public: //! set function to minimize virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1) = 0; -// //! set function and its gradient -// virtual void setFunctionAndGradient(boost::function<double(const double *)> fcn, boost::function<double(const double *, int)> fcn_deriv, int ndim = 1) = 0; - //! run minimization virtual void minimize() = 0; @@ -80,7 +77,6 @@ public: //! set function to minimize virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1) { m_fcn = fcn, m_ndim = ndim; } -// virtual void setFunctionAndGradient(boost::function<double(const double *)> fcn, boost::function<double(const double *, int)> fcn_deriv, int ndim = 1) { (void)fcn; (void)fcn_deriv; (void)ndim; throw NotImplementedException("TestMinimizer::setFunctionGradient"); } //! run minimization virtual void minimize() diff --git a/Core/Tools/src/FitSuite.cpp b/Core/Tools/src/FitSuite.cpp index f508af25b2a..bde80ad6744 100644 --- a/Core/Tools/src/FitSuite.cpp +++ b/Core/Tools/src/FitSuite.cpp @@ -84,13 +84,8 @@ void FitSuite::link_fit_parameters() void FitSuite::minimize() { // initializing minimizer with fcn function belonging to given class - //m_minimizer->setFunction( std::bind1st(std::mem_fun(&FitSuite::functionToMinimize), this), (int)m_fit_parameters.size() ); m_minimizer->setFunction( boost::bind(&FitSuite::functionToMinimize, this, _1), (int)m_fit_parameters.size() ); -// m_minimizer->setFunctionAndGradient(boost::bind(&FitSuite::functionToMinimize, this, _1), -// boost::bind(&FitSuite::functionGradient, this, _1, _2), -// (int)m_fit_parameters.size() ); - // propagating local fit parameters to the minimizer's internal list of parameters for(size_t i_par = 0; i_par<m_fit_parameters.size(); i_par++) { m_minimizer->setVariable((int)i_par, m_fit_parameters[i_par] ); @@ -171,9 +166,3 @@ double FitSuite::functionToMinimize(const double *pars_current_values) } -//double FitSuite::functionGradient(const double *pars, int icoord) -//{ -// std::cout << " functionGradient " << icoord << std::endl; -// return 0.0; -//} - -- GitLab