From 741b8433d5a5b19fd00f6eaee6a813b763960a2b Mon Sep 17 00:00:00 2001 From: pospelov <pospelov@fz-juelich.de> Date: Fri, 5 Oct 2012 18:44:19 +0200 Subject: [PATCH] FitSuite, part II --- App/inc/ROOTMinimizer.h | 13 +- App/src/ROOTMinimizer.cpp | 47 ++++++- App/src/TestMesoCrystal2.cpp | 65 +++++++-- Core/Algorithms/inc/Experiment.h | 7 + Core/Algorithms/src/Experiment.cpp | 11 ++ Core/Core.pro | 4 +- Core/Samples/inc/ICompositeSample.h | 4 + Core/Samples/inc/ISample.h | 5 +- Core/Samples/inc/ParameterPool.h | 14 +- Core/Samples/src/ISample.cpp | 10 +- Core/Samples/src/ParameterPool.cpp | 17 +-- Core/Tools/inc/FitMultiParameter.h | 6 + Core/Tools/inc/FitParameter.h | 37 +++++- Core/Tools/inc/FitSuite.h | 34 ++++- Core/Tools/inc/{Minimizer.h => IMinimizer.h} | 27 ++-- Core/Tools/src/FitMultiParameter.cpp | 8 ++ Core/Tools/src/FitParameter.cpp | 15 ++- Core/Tools/src/FitSuite.cpp | 131 ++++++++++++++++++- Core/Tools/src/IMinimizer.cpp | 5 + Core/Tools/src/Minimizer.cpp | 5 - 20 files changed, 402 insertions(+), 63 deletions(-) rename Core/Tools/inc/{Minimizer.h => IMinimizer.h} (61%) create mode 100644 Core/Tools/src/IMinimizer.cpp delete mode 100644 Core/Tools/src/Minimizer.cpp diff --git a/App/inc/ROOTMinimizer.h b/App/inc/ROOTMinimizer.h index 36e4115236c..bfdaa229725 100644 --- a/App/inc/ROOTMinimizer.h +++ b/App/inc/ROOTMinimizer.h @@ -14,25 +14,32 @@ //! @author Scientific Computing Group at FRM II //! @date 05.10.2012 -#include "Minimizer.h" +#include "IMinimizer.h" #include "OutputData.h" +#include <string> // from ROOT #include "Math/Minimizer.h" #include "Math/Factory.h" -#include <string> +#include "Math/Functor.h" + //- ------------------------------------------------------------------- //! @class FitSuite //! @brief Wrapper for ROOT minimizers to interface with our FitSuite //- ------------------------------------------------------------------- -class ROOTMinimizer : public Minimizer +class ROOTMinimizer : public IMinimizer { public: ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type); virtual ~ROOTMinimizer(); + virtual void setVariable(int i, const FitParameter *par) ; + virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1); + virtual void minimize(); + private: ROOT::Math::Minimizer *m_root_minimizer; + ROOT::Math::Functor *m_fcn; }; #endif // ROOTMINIMIZER_H diff --git a/App/src/ROOTMinimizer.cpp b/App/src/ROOTMinimizer.cpp index c6cc4cdf01b..9a04fbfc3c3 100644 --- a/App/src/ROOTMinimizer.cpp +++ b/App/src/ROOTMinimizer.cpp @@ -1,12 +1,57 @@ #include "ROOTMinimizer.h" +#include "Exceptions.h" -ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type) + +ROOTMinimizer::ROOTMinimizer(const std::string &minimizer_name, const std::string &algo_type) : m_fcn(0) { m_root_minimizer = ROOT::Math::Factory::CreateMinimizer(minimizer_name.c_str(), algo_type.c_str() ); + if( m_root_minimizer == 0 ) { + throw NullPointerException("ROOTMinimizer::ROOTMinimizer() -> Error! Can't build minimizer"); + } } ROOTMinimizer::~ROOTMinimizer() { delete m_root_minimizer; + delete m_fcn; +} + + +//void ROOTMinimizer::setVariable(int i, const std::string &name, double value, double step) +//{ +// m_root_minimizer->SetVariable(i, name.c_str(), value, step); +//} + +void ROOTMinimizer::setVariable(int i, const FitParameter *par) +{ + if(par->hasDoubleBound() ) { + m_root_minimizer->SetLimitedVariable(i, par->getName().c_str(), par->getValue(), par->getStep(), par->getLowerLimit(), par->getUpperLimit()); + } else if(par->hasLowerLimit() && !par->hasUpperLimit() ) { + m_root_minimizer->SetLowerLimitedVariable(i, par->getName().c_str(), par->getValue(), par->getStep(), par->getLowerLimit()); + } else if( par->hasUpperLimit() && !par->hasLowerLimit() ) { + m_root_minimizer->SetUpperLimitedVariable(i, par->getName().c_str(), par->getValue(), par->getStep(), par->getUpperLimit()); + } else if( !par->hasUpperLimit() && !par->hasLowerLimit() ) { + m_root_minimizer->SetVariable(i, par->getName().c_str(), par->getValue(), par->getStep()); + } else { + throw LogicErrorException("ROOTMinimizer::setVariable() -> I wish I knew how I got there..."); + } +} + + + +void ROOTMinimizer::minimize() +{ + m_root_minimizer->Minimize(); +} + + +/* ************************************************************************* */ +// set fcn function for minimizer +/* ************************************************************************* */ +void ROOTMinimizer::setFunction(boost::function<double(const double *)> fcn, int ndim) +{ + m_fcn = new ROOT::Math::Functor(fcn, ndim); + m_root_minimizer->SetFunction(*m_fcn); } + diff --git a/App/src/TestMesoCrystal2.cpp b/App/src/TestMesoCrystal2.cpp index f9420cbf9ae..41b65b76a00 100644 --- a/App/src/TestMesoCrystal2.cpp +++ b/App/src/TestMesoCrystal2.cpp @@ -21,6 +21,9 @@ #include "OutputDataIOFactory.h" #include "FitSuite.h" #include "ROOTMinimizer.h" +#include "SampleFactory.h" +#include "FitMultiParameter.h" +#include "TRange.h" #include "TCanvas.h" #include "TH2D.h" @@ -46,16 +49,30 @@ void TestMesoCrystal2::execute() { initializeExperiment(); +// m_sample->walk_and_print(); +// ParameterPool *p_param_pool = m_sample->createParameterTree(); +// std::cout << *p_param_pool; +// return; + m_experiment->runSimulation(); - IsGISAXSTools::drawLogOutputData(*(m_experiment->getOutputDataClone()), "c1_test_meso_crystal", "mesocrystal","CONT4 Z", "mesocrystal"); + m_experiment->normalize(); + + OutputData<double > *m_real_data = m_experiment->getOutputDataClone(); + +// IsGISAXSTools::setMinimum(1.); +// IsGISAXSTools::drawLogOutputData(*(m_experiment->getOutputDataClone()), "c1_test_meso_crystal", "mesocrystal","CONT4 Z", "mesocrystal"); //mp_exact_data = experiment.getOutputDataClone(); - return; +// return; // setting fitSuite FitSuite *fitSuite = new FitSuite(); - fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") ); fitSuite->setExperiment(m_experiment); + fitSuite->setRealData(*m_real_data); + //fitSuite->addFitParameter("*alpha",0.5, 0.1); + fitSuite->addFitParameter("*height",2.0, 0.1, TRange<double>(1.0, 10.0) ); + fitSuite->runFit(); + } @@ -66,15 +83,19 @@ void TestMesoCrystal2::execute() void TestMesoCrystal2::initializeExperiment() { delete m_experiment; - initializeSample(); + //initializeSample(); + + m_sample = SampleFactory::instance().createItem("IsGISAXS9_Pyramid"); + m_experiment = new GISASExperiment(); m_experiment->setSample( m_sample ); - m_experiment->setDetectorParameters(100, 0.3*Units::degree, 0.073, 100 , -0.4*Units::degree, 0.066); + m_experiment->setDetectorParameters(200, 0.3*Units::degree, 0.073, 200, -0.4*Units::degree, 0.066); m_experiment->setBeamParameters(1.77*Units::angstrom, -0.4*Units::degree, 0.0*Units::degree); - m_experiment->setBeamIntensity(1e7); + m_experiment->setBeamIntensity(1e10); } + /* ************************************************************************* */ // /* ************************************************************************* */ @@ -91,8 +112,8 @@ void TestMesoCrystal2::initializeSample() complex_t n_avg = std::sqrt(surface_filling_ratio*avg_n_squared_meso + 1.0 - surface_filling_ratio); complex_t n_particle_adapted = std::sqrt(n_avg*n_avg + n_particle*n_particle - 1.0); FormFactorCylinder ff_cyl(0.5*Units::micrometer, meso_width); - double sigma_h = 40*Units::nanometer; - double sigma_r = 120*Units::nanometer; + double sigma_h = 4*Units::nanometer; + double sigma_r = 50*Units::nanometer; FormFactorDecoratorDebyeWaller ff_meso(ff_cyl.clone(), sigma_h*sigma_h/2.0, sigma_r*sigma_r/2.0); // Create multilayer @@ -112,12 +133,26 @@ void TestMesoCrystal2::initializeSample() substrate_layer.setMaterial(p_substrate_material); IInterferenceFunction *p_interference_funtion = new InterferenceFunctionNone(); ParticleDecoration particle_decoration; - size_t n_max_phi_rotation_steps = 140; + size_t n_max_phi_rotation_steps = 180; size_t n_alpha_rotation_steps = 1; double alpha_step = 5.0*Units::degree/n_alpha_rotation_steps; double alpha_start = - (n_alpha_rotation_steps/2.0)*alpha_step; + // with optimization +// TRange<double> phi_range(0.0, 2.0*M_PI/3.0); +// TRange<double> z_range(0.0, 2.6); +// double max_rho = 2.6; +// std::vector<double> phi_angles = createLattice(6.1*Units::nanometer)->collectBraggAngles(n_max_phi_rotation_steps, max_rho, phi_range, z_range); +// for (size_t i=0; i<phi_angles.size(); ++i) { +// for (size_t j=0; j<n_alpha_rotation_steps; ++j) { +// Geometry::RotateZ3D transform1(phi_angles[i]); +// Geometry::RotateY3D transform2(alpha_start + j*alpha_step); +// Geometry::Transform3D *p_total_transform = new Geometry::Transform3D(transform1*transform2); +// particle_decoration.addParticle(createMesoCrystal(6.1*Units::nanometer, +// n_particle_adapted, &ff_meso), p_total_transform, 0.2*Units::micrometer); +// } +// } // without double phi_step = 2.0*M_PI/3.0/n_max_phi_rotation_steps; double phi_start = 0.0; @@ -126,7 +161,7 @@ void TestMesoCrystal2::initializeSample() Geometry::RotateZ3D transform1(phi_start + (double)i*phi_step); Geometry::RotateY3D transform2(alpha_start + j*alpha_step); Geometry::Transform3D *p_total_transform = new Geometry::Transform3D(transform1); - particle_decoration.addParticle(createMesoCrystal(6.1*Units::nanometer, + particle_decoration.addParticle(createMesoCrystal(4.7*Units::nanometer, n_particle_adapted, &ff_meso), p_total_transform, 0.5*Units::micrometer); } } @@ -135,7 +170,7 @@ void TestMesoCrystal2::initializeSample() particle_decoration.addInterferenceFunction(p_interference_funtion); LayerDecorator avg_layer_decorator(avg_layer, particle_decoration); - LayerRoughness roughness(0.5*Units::nanometer, 0.3, 500.0*Units::nanometer); + LayerRoughness roughness(2.0*Units::nanometer, 0.3, 500.0*Units::nanometer); p_multi_layer->addLayer(air_layer); p_multi_layer->addLayer(avg_layer_decorator); @@ -147,6 +182,10 @@ void TestMesoCrystal2::initializeSample() m_sample = p_multi_layer; + + std::cout << "Average layer index: " << n_avg << std::endl; + std::cout << "Adapted particle index: " << n_particle_adapted << std::endl; + std::cout << "Substrate layer index: " << n_substrate << std::endl; } MesoCrystal* TestMesoCrystal2::createMesoCrystal(double stacking_radius, complex_t n_particle, const IFormFactor* p_meso_form_factor) @@ -155,8 +194,8 @@ MesoCrystal* TestMesoCrystal2::createMesoCrystal(double stacking_radius, complex kvector_t bas_a = p_lat->getBasisVectorA(); kvector_t bas_b = p_lat->getBasisVectorB(); kvector_t bas_c = p_lat->getBasisVectorC(); - double sigma = 1.0*Units::nanometer; - Particle particle(n_particle, new FormFactorSphereGaussianRadius(stacking_radius-1.0*Units::nanometer, sigma)); + double sigma = 0.2*Units::nanometer; + Particle particle(n_particle, new FormFactorSphereGaussianRadius(stacking_radius-0.5*Units::nanometer, sigma)); kvector_t position_0 = kvector_t(0.0, 0.0, 0.0); kvector_t position_1 = 1.0/3.0*(2.0*bas_a + bas_b + bas_c); kvector_t position_2 = 1.0/3.0*(bas_a + 2.0*bas_b + 2.0*bas_c); diff --git a/Core/Algorithms/inc/Experiment.h b/Core/Algorithms/inc/Experiment.h index 2a14e34e473..0d174cb3fef 100644 --- a/Core/Algorithms/inc/Experiment.h +++ b/Core/Algorithms/inc/Experiment.h @@ -35,6 +35,10 @@ public: /// Set the sample to be tested void setSample(ISample *p_sample); + /// get the sample + ISample *getSample() { return mp_sample; } + const ISample *getSample() const { return mp_sample; } + /// Get data structure that contains the intensity map on the detector for all scan parameters OutputData<double>* getOutputDataClone() const; @@ -63,6 +67,9 @@ public: const OutputData<double>* getOutputDataMask() const; + //! create combined parameter pool of the whole sample and experiment (if any) + virtual ParameterPool *createParameterTree() const; + protected: /// Default implementation only adds the detector axes virtual void updateIntensityMapAxes(); diff --git a/Core/Algorithms/src/Experiment.cpp b/Core/Algorithms/src/Experiment.cpp index f7c7600cbde..46776a222dc 100644 --- a/Core/Algorithms/src/Experiment.cpp +++ b/Core/Algorithms/src/Experiment.cpp @@ -97,3 +97,14 @@ void Experiment::setOutputDataMask(size_t n_chunks_total, size_t n_chunk ) } } + + +/* ************************************************************************* */ +// create combined parameter pool of the whole sample and experiment (if any) +/* ************************************************************************* */ +ParameterPool *Experiment::createParameterTree() const +{ + if(mp_sample == 0 ) throw NullPointerException("Experiment::createParameterTree() -> Error! Sample is absent"); + return mp_sample->createParameterTree(); +} + diff --git a/Core/Core.pro b/Core/Core.pro index 15c37992dc8..d9a009292f2 100644 --- a/Core/Core.pro +++ b/Core/Core.pro @@ -111,7 +111,7 @@ SOURCES += \ PythonAPI/src/PythonModule.cpp \ PythonAPI/src/PythonPlusplusHelper.cpp \ PythonAPI/src/PythonOutputData.cpp \ - Tools/src/Minimizer.cpp + Tools/src/IMinimizer.cpp HEADERS += \ Algorithms/inc/Beam.h \ @@ -235,7 +235,7 @@ HEADERS += \ PythonAPI/inc/PythonModule.h \ PythonAPI/inc/PythonOutputData.h \ PythonAPI/inc/PythonPlusplusHelper.h \ - Tools/inc/Minimizer.h + Tools/inc/IMinimizer.h INCLUDEPATH += ./Algorithms/inc ./FormFactors/inc/ ./Geometry/inc ./Samples/inc ./Tools/inc ./PythonAPI/inc DEPENDPATH += ./Algorithms/inc ./FormFactors/inc/ ./Geometry/inc ./Samples/inc ./Tools/inc ./PythonAPI/inc diff --git a/Core/Samples/inc/ICompositeSample.h b/Core/Samples/inc/ICompositeSample.h index 80cec442393..592f3f19df8 100644 --- a/Core/Samples/inc/ICompositeSample.h +++ b/Core/Samples/inc/ICompositeSample.h @@ -31,11 +31,13 @@ public: //! definition of container for registered children typedef std::list<ISample *> samples_t; typedef samples_t::iterator iterator_t; + typedef samples_t::const_iterator const_iterator_t; ICompositeSample(); //! to confirm compound nature of given class virtual ICompositeSample *getCompositeSample() { return this; } + virtual const ICompositeSample *getCompositeSample() const { return this; } //! register child in the container virtual void registerChild(ISample *sample); @@ -46,6 +48,8 @@ public: //! iteration over local registered children iterator_t begin_shallow() { return m_samples.begin(); } iterator_t end_shallow() { return m_samples.end(); } + const_iterator_t begin_shallow() const { return m_samples.begin(); } + const_iterator_t end_shallow() const { return m_samples.end(); } //! size of children virtual size_t size() const { return m_samples.size(); } diff --git a/Core/Samples/inc/ISample.h b/Core/Samples/inc/ISample.h index b36aeba266a..780a2672412 100644 --- a/Core/Samples/inc/ISample.h +++ b/Core/Samples/inc/ISample.h @@ -35,6 +35,7 @@ public: //! return pointer to "this", if it is composite sample (to overload) virtual ICompositeSample *getCompositeSample() { return 0; } + virtual const ICompositeSample *getCompositeSample() const { return 0; } //! clone sample (to overload) virtual ISample *clone() const; @@ -43,7 +44,7 @@ public: ParameterPool *getParameterPool() { return &m_parameters; } //! create new parameter pool which contains all local parameter and parameters of children - virtual ParameterPool *createParameterTree(); + virtual ParameterPool *createParameterTree() const; //! check if this sample (or one of its subsamples) contains elements requiring DWBA corrections and return an ISimulation to calculate this virtual DWBASimulation *createDWBASimulation() const { return 0; } @@ -62,7 +63,7 @@ protected: virtual void init_parameters(); //! add parameters from local pool to external pool and call recursion over direct children - virtual void addParametersToExternalPool(std::string path, ParameterPool *external_pool, int copy_number=-1); + virtual void addParametersToExternalPool(std::string path, ParameterPool *external_pool, int copy_number=-1) const; ParameterPool m_parameters; //! parameter pool }; diff --git a/Core/Samples/inc/ParameterPool.h b/Core/Samples/inc/ParameterPool.h index b3a8c07304c..7beb0bf7135 100644 --- a/Core/Samples/inc/ParameterPool.h +++ b/Core/Samples/inc/ParameterPool.h @@ -52,10 +52,10 @@ public: ParameterPool *clone(); //! clone with adding preffix to every parameter key - ParameterPool *cloneWithPrefix(std::string prefix); + ParameterPool *cloneWithPrefix(const std::string &prefix); //! copy parameters of given pool to the external pool while adding prefix to local parameter keys - void copyToExternalPool(std::string prefix, ParameterPool *external_pool); + void copyToExternalPool(const std::string &prefix, ParameterPool *external_pool) const; //! clear and delete parameter map void clear(); @@ -64,10 +64,10 @@ public: size_t size() const { return m_map.size(); } //! main method to register data address in the pool - bool registerParameter(std::string name, double *parpointer); + bool registerParameter(const std::string &name, double *parpointer); //! add parameter to the pool - bool addParameter(std::string name, RealPar par); + bool addParameter(const std::string &name, RealPar par); //! access to parameter container iterator_t begin() { return m_map.begin(); } @@ -78,13 +78,13 @@ public: const_iterator_t end() const { return m_map.end(); } //! return parameter with given name - RealPar getParameter(std::string name) const; + RealPar getParameter(const std::string &name) const; //! set parameter value, return true in the case of success - bool setParameterValue(std::string name, double value); + bool setParameterValue(const std::string &name, double value); //! set parameter value, return number of changed parameters - int setMatchedParametersValue(std::string wildcards, double value); + int setMatchedParametersValue(const std::string &wildcards, double value); //! print parameter pool friend std::ostream &operator<<(std::ostream &ostr, const ParameterPool &obj) diff --git a/Core/Samples/src/ISample.cpp b/Core/Samples/src/ISample.cpp index 8a243c0904a..dc6a370236a 100644 --- a/Core/Samples/src/ISample.cpp +++ b/Core/Samples/src/ISample.cpp @@ -56,7 +56,7 @@ void ISample::init_parameters() // create new parameter pool which contains all local parameter and parameters of children // user have to delete it /* ************************************************************************* */ -ParameterPool *ISample::createParameterTree() +ParameterPool *ISample::createParameterTree() const { ParameterPool *newpool = new ParameterPool; std::string path("/"); @@ -68,7 +68,7 @@ ParameterPool *ISample::createParameterTree() /* ************************************************************************* */ // add parameters from local pool to external pool and call recursion over direct children /* ************************************************************************* */ -void ISample::addParametersToExternalPool(std::string path, ParameterPool *external_pool, int copy_number) +void ISample::addParametersToExternalPool(std::string path, ParameterPool *external_pool, int copy_number) const { // adding trailing slash, if it is not already there if( path[path.length()-1] != '/' ) path += "/"; @@ -82,19 +82,19 @@ void ISample::addParametersToExternalPool(std::string path, ParameterPool *exter m_parameters.copyToExternalPool(path, external_pool); // going through direct children of given sample and copy they parameters recursively - ICompositeSample *sample = getCompositeSample(); + const ICompositeSample *sample = getCompositeSample(); if( sample ) { // Here we need some default mechanism to handle cases with many children with same name. // Lets run through all direct children and save they names Utils::StringUsageMap strUsageMap; - for(ICompositeSample::iterator_t it=sample->begin_shallow(); it!=sample->end_shallow(); ++it) { + for(ICompositeSample::const_iterator_t it=sample->begin_shallow(); it!=sample->end_shallow(); ++it) { strUsageMap.add( path +(*it)->getName() ); // saving children name } // Now we run through direct children again, and assign copy number for all children with same name Utils::StringUsageMap strUsageMap2; - for(ICompositeSample::iterator_t it=sample->begin_shallow(); it!=sample->end_shallow(); ++it) { + for(ICompositeSample::const_iterator_t it=sample->begin_shallow(); it!=sample->end_shallow(); ++it) { std::string children_name = path +(*it)->getName(); strUsageMap2.add(children_name); int ncopy = strUsageMap2[children_name]-1; // staring from 0 diff --git a/Core/Samples/src/ParameterPool.cpp b/Core/Samples/src/ParameterPool.cpp index b11111fba6f..bd8736dbe4a 100644 --- a/Core/Samples/src/ParameterPool.cpp +++ b/Core/Samples/src/ParameterPool.cpp @@ -48,13 +48,13 @@ void ParameterPool::clear() // Registering parameter with given name. Name should be different for each // register. /* ************************************************************************* */ -bool ParameterPool::registerParameter(std::string name, double *parameter_address) +bool ParameterPool::registerParameter(const std::string &name, double *parameter_address) { RealPar par(parameter_address); return addParameter(name, par); } -bool ParameterPool::addParameter(std::string name, RealPar par) +bool ParameterPool::addParameter(const std::string &name, RealPar par) { parametermap_t::iterator it = m_map.find(name); if( it!=m_map.end() ) { @@ -79,7 +79,7 @@ ParameterPool *ParameterPool::clone() /* ************************************************************************* */ // Clone method that adds a prefix to parameter's keys /* ************************************************************************* */ -ParameterPool *ParameterPool::cloneWithPrefix(std::string prefix) +ParameterPool *ParameterPool::cloneWithPrefix(const std::string &prefix) { ParameterPool *new_pool = new ParameterPool; for(parametermap_t::iterator it=m_map.begin(); it!= m_map.end(); ++it) @@ -93,20 +93,21 @@ ParameterPool *ParameterPool::cloneWithPrefix(std::string prefix) // copy parameters of given pool to the external pool while adding prefix to // local parameter keys /* ************************************************************************* */ -void ParameterPool::copyToExternalPool(std::string prefix, ParameterPool *external_pool) +void ParameterPool::copyToExternalPool(const std::string &prefix, ParameterPool *external_pool) const { - for(parametermap_t::iterator it=m_map.begin(); it!= m_map.end(); ++it) { + for(parametermap_t::const_iterator it=m_map.begin(); it!= m_map.end(); ++it) { external_pool->addParameter(prefix+it->first, it->second); } } /* ************************************************************************* */ -ParameterPool::RealPar ParameterPool::getParameter(std::string name) const +ParameterPool::RealPar ParameterPool::getParameter(const std::string &name) const { parametermap_t::const_iterator it = m_map.find(name); if( it!=m_map.end() ) { return (*it).second; } else { + std::cout << "ParameterPool::getParameter() -> Warning. No parameter with name '" << name << std::endl; return RealPar(0); } } @@ -114,7 +115,7 @@ ParameterPool::RealPar ParameterPool::getParameter(std::string name) const /* ************************************************************************* */ // set parameter value /* ************************************************************************* */ -bool ParameterPool::setParameterValue(std::string name, double value) +bool ParameterPool::setParameterValue(const std::string &name, double value) { RealPar x = getParameter(name); if( x.isNull() ) { @@ -129,7 +130,7 @@ bool ParameterPool::setParameterValue(std::string name, double value) /* ************************************************************************* */ // set parameter value /* ************************************************************************* */ -int ParameterPool::setMatchedParametersValue(std::string wildcards, double value) +int ParameterPool::setMatchedParametersValue(const std::string &wildcards, double value) { int npars(0); for(iterator_t it=begin(); it!= end(); it++) { diff --git a/Core/Tools/inc/FitMultiParameter.h b/Core/Tools/inc/FitMultiParameter.h index 4a2bf284c56..75ad390d70b 100644 --- a/Core/Tools/inc/FitMultiParameter.h +++ b/Core/Tools/inc/FitMultiParameter.h @@ -33,6 +33,7 @@ public: typedef std::vector<parameter_t > parametercoll_t; FitMultiParameter(); + FitMultiParameter(const std::string &name, double value, double step, double error=0.0); virtual void setValue(double value) { @@ -55,6 +56,11 @@ public: friend std::ostream &operator<<(std::ostream &ostr, const FitMultiParameter &m) { m.print(ostr); return ostr; } protected: + // disabled copy constructor + FitMultiParameter(const FitMultiParameter &other); + // disabled assignment operator + FitMultiParameter &operator=(const FitMultiParameter &other); + //! print class void print(std::ostream &ostr) const; diff --git a/Core/Tools/inc/FitParameter.h b/Core/Tools/inc/FitParameter.h index 7bbc74951ea..bc2581b91b4 100644 --- a/Core/Tools/inc/FitParameter.h +++ b/Core/Tools/inc/FitParameter.h @@ -27,15 +27,42 @@ class FitParameter : public INamed { public: FitParameter(); + FitParameter(const std::string &name, double value, double step, double error=0.0); FitParameter(const FitParameter &other); FitParameter &operator=(const FitParameter &other); virtual ~FitParameter(){} - //! return value of parameter - virtual double getValue() const { return m_value;} - //! set value of parameter virtual void setValue(double value) { m_value = value;} + //! get value of parameter + virtual double getValue() const { return m_value;} + + //! set parameter step for minimizer + virtual void setStep(double value) { m_step = value; } + //! get parameter step for minimizer + virtual double getStep() const { return m_step;} + + //! set lower and upper limits + virtual void setLimits(double xmin, double xmax) { setLowerLimit(xmin); setUpperLimit(xmax); } + + //! set lower limit + virtual void setLowerLimit(double value) { m_lower_limit = value; m_has_lower_limit = true; } + //! get lower limit + double getLowerLimit() const { return m_lower_limit; } + + //! set upper limit + virtual void setUpperLimit(double value) { m_upper_limit = value; m_has_upper_limit = true; } + //! get upper limit + double getUpperLimit() const { return m_upper_limit; } + + //! if has lower limit + virtual bool hasLowerLimit() const { return m_has_lower_limit; } + + //! if has upper limit + virtual bool hasUpperLimit() const { return m_has_upper_limit; } + + //! if has upper limit + virtual bool hasDoubleBound() const { return (m_has_lower_limit && m_has_upper_limit); } //! print class friend std::ostream &operator<<(std::ostream &ostr, const FitParameter &m) { m.print(ostr); return ostr; } @@ -45,9 +72,13 @@ protected: void print(std::ostream &ostr) const; double m_value; //! parameter value + double m_step; //! parameter step size double m_error; //! parameter error + bool m_has_lower_limit; //! paramter has lower bound + bool m_has_upper_limit; //! paramter has upper bound double m_lower_limit; //! minimum allowed value double m_upper_limit; //! maximum allowed value }; + #endif // FITPARAMETER_H diff --git a/Core/Tools/inc/FitSuite.h b/Core/Tools/inc/FitSuite.h index 51bc7401643..a1c9f1ec682 100644 --- a/Core/Tools/inc/FitSuite.h +++ b/Core/Tools/inc/FitSuite.h @@ -16,8 +16,15 @@ #include "OutputData.h" +#include "TRange.h" +#include "FitMultiParameter.h" +#include <string> + class Experiment; -class Minimizer; +class IMinimizer; +class ParameterPool; +class ChiSquaredModule; + //- ------------------------------------------------------------------- //! @class FitSuite @@ -26,13 +33,34 @@ class Minimizer; class FitSuite { public: + typedef std::vector<FitMultiParameter *> fitmultiparameters_t; FitSuite(); + virtual ~FitSuite(); void setExperiment(Experiment *experiment) { m_experiment = experiment; } - void setMinimizer(Minimizer *minimizer) { m_minimizer = minimizer; } + void setMinimizer(IMinimizer *minimizer) { m_minimizer = minimizer; } + + //! add fit parameter + FitMultiParameter *addFitParameter(const std::string &name, double value, double step, double error=0.0); + FitMultiParameter *addFitParameter(const std::string &name, double value, double step, const TRange<double> &range); + + //! initialize fitting parameters + virtual void init_fit_parameters(); + + //! set real data + void setRealData(const OutputData<double> &data); + + //! run fit + virtual void runFit(); + + //! function to minimize + double functionToMinimize(const double *pars_current_values); + private: Experiment *m_experiment; - Minimizer *m_minimizer; + IMinimizer *m_minimizer; + fitmultiparameters_t m_fit_params; + ChiSquaredModule *m_chi2_module; }; #endif // FITSUITE_H diff --git a/Core/Tools/inc/Minimizer.h b/Core/Tools/inc/IMinimizer.h similarity index 61% rename from Core/Tools/inc/Minimizer.h rename to Core/Tools/inc/IMinimizer.h index 7d3f0072794..15e40512d72 100644 --- a/Core/Tools/inc/Minimizer.h +++ b/Core/Tools/inc/IMinimizer.h @@ -1,5 +1,5 @@ -#ifndef MINIMIZER_H -#define MINIMIZER_H +#ifndef IMINIMIZER_H +#define IMINIMIZER_H // ******************************************************************** // * The BornAgain project * // * Simulation of neutron and x-ray scattering at grazing incidence * @@ -9,21 +9,30 @@ // * eget quam orci. Quisque porta varius dui, quis posuere nibh * // * mollis quis. Mauris commodo rhoncus porttitor. * // ******************************************************************** -//! @file Minimizer.h +//! @file IMinimizer.h //! @brief Definition of Minimizer class //! @author Scientific Computing Group at FRM II //! @date 05.10.2012 +#include "FitParameter.h" +#include <boost/function.hpp> + + //- ------------------------------------------------------------------- -//! @class Minimizer -//! @brief Common wrapper for all kind minimizer's +//! @class IMinimizer +//! @brief Common interface for all kind minimizer's //- ------------------------------------------------------------------- -class Minimizer +class IMinimizer { public: - Minimizer(); - virtual ~Minimizer(){} + IMinimizer(); + virtual ~IMinimizer(){} + + virtual void setVariable(int i, const FitParameter *par) = 0; + //virtual void setVariable(int i, const std::string &name, double value, double step) = 0; + virtual void minimize() = 0; + virtual void setFunction(boost::function<double(const double *)> fcn, int ndim=1) = 0; }; -#endif // MINIMIZER_H +#endif // IMINIMIZER_H diff --git a/Core/Tools/src/FitMultiParameter.cpp b/Core/Tools/src/FitMultiParameter.cpp index 709ca7f07ec..b79ee6e4f07 100644 --- a/Core/Tools/src/FitMultiParameter.cpp +++ b/Core/Tools/src/FitMultiParameter.cpp @@ -1,11 +1,19 @@ #include "FitMultiParameter.h" #include "Utils.h" + FitMultiParameter::FitMultiParameter() { setName("FitMultiParameter"); } + +FitMultiParameter::FitMultiParameter(const std::string &name, double value, double step, double error) : FitParameter(name, value, step, error) +{ + +} + + /* ************************************************************************* */ //! add real parameter to the collection /* ************************************************************************* */ diff --git a/Core/Tools/src/FitParameter.cpp b/Core/Tools/src/FitParameter.cpp index c2aa5581af8..426dc299222 100644 --- a/Core/Tools/src/FitParameter.cpp +++ b/Core/Tools/src/FitParameter.cpp @@ -2,17 +2,27 @@ #include <iostream> -FitParameter::FitParameter() : m_value(0), m_error(0), m_lower_limit(0), m_upper_limit(0) +FitParameter::FitParameter() : m_value(0), m_step(0), m_error(0), + m_has_lower_limit(false), m_has_upper_limit(false), m_lower_limit(0), m_upper_limit(0) { setName("FitParameter"); } +FitParameter::FitParameter(const std::string &name, double value, double step, double error) : INamed(name), m_value(value), m_step(step), m_error(error), + m_has_lower_limit(false), m_has_upper_limit(false), m_lower_limit(0), m_upper_limit(0) +{ + +} + FitParameter::FitParameter(const FitParameter &other) : INamed(other) { setName("FitParameter"); m_value = other.m_value; + m_step = other.m_step; m_error = other.m_error; + m_has_lower_limit = other.m_has_lower_limit; + m_has_upper_limit = other.m_has_upper_limit; m_lower_limit = other.m_lower_limit; m_upper_limit = other.m_upper_limit; } @@ -24,7 +34,10 @@ FitParameter &FitParameter::operator=(const FitParameter &other) { INamed::operator=(other); m_value = other.m_value; + m_step = other.m_step; m_error = other.m_error; + m_has_lower_limit = other.m_has_lower_limit; + m_has_upper_limit = other.m_has_upper_limit; m_lower_limit = other.m_lower_limit; m_upper_limit = other.m_upper_limit; } diff --git a/Core/Tools/src/FitSuite.cpp b/Core/Tools/src/FitSuite.cpp index 6330d857923..d93d8d8a84c 100644 --- a/Core/Tools/src/FitSuite.cpp +++ b/Core/Tools/src/FitSuite.cpp @@ -1,5 +1,134 @@ #include "FitSuite.h" +#include "Exceptions.h" +#include "FitMultiParameter.h" +#include "ParameterPool.h" +#include "Experiment.h" +#include "IMinimizer.h" +#include "ChiSquaredModule.h" -FitSuite::FitSuite() : m_experiment(0), m_minimizer(0) + +FitSuite::FitSuite() : m_experiment(0), m_minimizer(0), m_chi2_module(0) +{ +} + + +FitSuite::~FitSuite() +{ + for(fitmultiparameters_t::iterator it = m_fit_params.begin(); it!= m_fit_params.end(); ++it) { + delete (*it); + } + delete m_minimizer; + delete m_chi2_module; +} + + +/* ************************************************************************* */ +// add fit parameter +/* ************************************************************************* */ +FitMultiParameter * FitSuite::addFitParameter(const std::string &name, double value, double step, double error) +{ + // making MultiFitParameter for later access to sample's parameters + FitMultiParameter *par = new FitMultiParameter(name, value, step, error); + m_fit_params.push_back( par ); + return par; +} + +FitMultiParameter * FitSuite::addFitParameter(const std::string &name, double value, double step, const TRange<double> &range) { + // making MultiFitParameter for later access to sample's parameters + FitMultiParameter *par = new FitMultiParameter(name, value, step); + par->setLimits(range.getMin(), range.getMax() ); + m_fit_params.push_back( par ); + return par; } + + + + +/* ************************************************************************* */ +// linking defined FitMultiParameters with sample parameters from pool +// initialising minimiser's fit parameters +/* ************************************************************************* */ +void FitSuite::init_fit_parameters() +{ + if( !m_experiment ) throw NullPointerException("FitSuite::init_fit_parameters() -> Error! Experiment is absent."); + if( !m_minimizer ) throw NullPointerException("FitSuite::init_fit_parameters() -> Error! Minimizer is absent."); + if( m_fit_params.empty() ) throw NullPointerException("FitSuite::init_fit_parameters() -> Error! No fit parameters defined"); + + // Linking FitMultiParameters with sample parameters from parameter pool. + ParameterPool *pool = m_experiment->createParameterTree(); + for(fitmultiparameters_t::iterator it = m_fit_params.begin(); it!= m_fit_params.end(); ++it) { + FitMultiParameter *par = (*it); + // name of FitMultiParameter is used to find in the pool links to sample's parameters + par->addMatchedParametersFromPool(par->getName(), pool); + } + + // propagating fit parameters in the minimizer + int npar(0); + for(fitmultiparameters_t::iterator it = m_fit_params.begin(); it!= m_fit_params.end(); ++it) { + m_minimizer->setVariable(npar++, (*it) ); + } + + delete pool; +} + + +/* ************************************************************************* */ +// linking defined FitMultiParameters with sample parameters from pool +// initialising minimiser's fit parameters +/* ************************************************************************* */ +void FitSuite::setRealData(const OutputData<double> &data) +{ + delete m_chi2_module; + m_chi2_module = new ChiSquaredModule(data); +} + + +/* ************************************************************************* */ +// run fit +/* ************************************************************************* */ +void FitSuite::runFit() +{ + if( !m_experiment || !m_minimizer || !m_chi2_module || m_fit_params.empty()) throw NullPointerException("FitSuite::runFit() -> Error! Something is wrong."); + + m_experiment->getSample()->walk_and_print(); + + // initializing fit parameters + init_fit_parameters(); + + std::cout << m_fit_params[0]->getName() << std::endl; + m_fit_params[0]->setValue(999); + m_experiment->getSample()->walk_and_print(); + + + // initializing minimizer with fcn function belonging to given class + m_minimizer->setFunction( std::bind1st(std::mem_fun(&FitSuite::functionToMinimize), this), m_fit_params.size() ); + + // run minimizer + m_minimizer->minimize(); + +} + + +/* ************************************************************************* */ +// function to minimize +/* ************************************************************************* */ +double FitSuite::functionToMinimize(const double *pars_current_values) +{ + std::cout << " FitSuit::functionToMinimize() -> Info. We are here." << std::endl; + + for(size_t i=0; i<m_fit_params.size(); ++i) { + std::cout << i << " " << pars_current_values[i] << std::endl; + m_fit_params[i]->setValue(pars_current_values[i]); + } + + m_experiment->runSimulation(); + m_experiment->normalize(); + + const OutputData<double> *p_simulated_data = m_experiment->getOutputData(); + double chi_squared = m_chi2_module->calculateChiSquared(p_simulated_data); + + std::cout << "chi squared = " << chi_squared << std::endl; + return chi_squared; +} + diff --git a/Core/Tools/src/IMinimizer.cpp b/Core/Tools/src/IMinimizer.cpp new file mode 100644 index 00000000000..223391c988b --- /dev/null +++ b/Core/Tools/src/IMinimizer.cpp @@ -0,0 +1,5 @@ +#include "IMinimizer.h" + +IMinimizer::IMinimizer() +{ +} diff --git a/Core/Tools/src/Minimizer.cpp b/Core/Tools/src/Minimizer.cpp deleted file mode 100644 index 12b9bd97d94..00000000000 --- a/Core/Tools/src/Minimizer.cpp +++ /dev/null @@ -1,5 +0,0 @@ -#include "Minimizer.h" - -Minimizer::Minimizer() -{ -} -- GitLab