diff --git a/App/App.pro b/App/App.pro
index cb555dc6d86055813a10c09922a19d3546dbcb84..caf94fb33fcf7460f994378ba9f8e3ba78b514d0 100644
--- a/App/App.pro
+++ b/App/App.pro
@@ -34,6 +34,7 @@ SOURCES += \
src/TestFormFactor.cpp \
src/TestFourier.cpp \
src/TestFresnelCoeff.cpp \
+ src/TestFumiliLMA.cpp \
src/TestIsGISAXS1.cpp \
src/TestIsGISAXS2.cpp \
src/TestIsGISAXS3.cpp \
@@ -75,6 +76,7 @@ HEADERS += \
inc/TestFormFactor.h \
inc/TestFourier.h \
inc/TestFresnelCoeff.h \
+ inc/TestFumiliLMA.h \
inc/TestIsGISAXS1.h \
inc/TestIsGISAXS2.h \
inc/TestIsGISAXS3.h \
diff --git a/App/inc/TestFumiliLMA.h b/App/inc/TestFumiliLMA.h
new file mode 100644
index 0000000000000000000000000000000000000000..4d9b4228bd54fadfb4b770b2d34293b9655f3df7
--- /dev/null
+++ b/App/inc/TestFumiliLMA.h
@@ -0,0 +1,135 @@
+#ifndef TESTFUMILILMA_H
+#define TESTFUMILILMA_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 TestFittingModule1.h
+//! @brief Definition of TestFittingModule class
+//! @author Scientific Computing Group at FRM II
+//! @date Dec 5, 2012
+
+#include "MathFunctions.h"
+#include "IFunctionalTest.h"
+#include "Math/Minimizer.h"
+#include "Math/Factory.h"
+#include "Math/FitMethodFunction.h"
+#include "OutputData.h"
+#include "Math/Functor.h"
+#include "gsl/gsl_sf_trig.h"
+#include "Fit/Fitter.h"
+#include <cmath>
+
+class TF2;
+class IFunctionObject;
+class IChiSquaredModule;
+class TCanvas;
+
+
+//- -------------------------------------------------------------------
+//! @class TestFumiliLMA
+//! @brief Test of ROOT's LMA-based minimizers Fumili and GSLMultiFit
+//- -------------------------------------------------------------------
+class TestFumiliLMA : public IFunctionalTest
+{
+public:
+ friend class MyChi2Function;
+
+ TestFumiliLMA();
+ virtual ~TestFumiliLMA();
+ virtual void execute();
+
+ 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
+ ROOT::Math::Functor *m_fcn; //! fit function
+ IFunctionObject *m_func_object; //! simulation function
+ TF2 *m_func; //! ROOT representation of the simulation function with min, max defined
+ IChiSquaredModule *m_chi_module; //! chi squared module
+ OutputData<double > *m_real_data; //! real data
+ size_t m_ndim; //! number of fit parametrs
+ double m_sigma; //! gaussian noise
+ TCanvas *m_c1;
+};
+
+
+class IFunctionObject
+{
+public:
+ virtual ~IFunctionObject(){}
+ virtual double operator()(const double *xx, const double *pars ) = 0;
+};
+
+
+class RosenBrockFunctionObject : public IFunctionObject
+{
+public:
+ double operator()(const double *xx, const double *pars )
+ {
+ const double x = xx[0];
+ const double y = xx[1];
+ const double tmp1 = y-x*x;
+ const double tmp2 = 1-x;
+ const double p0=pars[0];
+ const double p1=pars[1];
+ return p0*100*tmp1*tmp1 + p1*tmp2*tmp2;
+ }
+};
+
+
+class SincXSincYFunctionObject : public IFunctionObject
+{
+public:
+ double operator()(const double *xx, const double *pars )
+ {
+ const double x = xx[0];
+ const double y = xx[1];
+ const double p0=pars[0];
+ const double p1=pars[1];
+ const double p2=pars[2];
+ double value = p0 * MathFunctions::Sinc(x-p1) * MathFunctions::Sinc(y-p2);
+ return value;
+ }
+};
+
+
+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 // TESTFUMILILMA_H
diff --git a/App/src/FunctionalTestFactory.cpp b/App/src/FunctionalTestFactory.cpp
index 270d2025ec5ee8084e270f60f6db107cfd80bc87..f87a1d06a54f4f5e3a416ca8b888f69abc3e44f9 100644
--- a/App/src/FunctionalTestFactory.cpp
+++ b/App/src/FunctionalTestFactory.cpp
@@ -2,6 +2,7 @@
#include "TestRoughness.h"
#include "TestFresnelCoeff.h"
#include "TestFormFactor.h"
+#include "TestFumiliLMA.h"
#include "TestDiffuseReflection.h"
#include "TestIsGISAXS1.h"
#include "TestIsGISAXS2.h"
@@ -87,6 +88,8 @@ FunctionalTestFactory::FunctionalTestFactory() : m_benchmark(0)
"functional test: test of minimizers with hard-to-minimize test functions");
registerItem("fourier", IFactoryCreateFunction<TestFourier, IFunctionalTest>,
"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");
m_benchmark = new TBenchmark();
}
diff --git a/App/src/TestFumiliLMA.cpp b/App/src/TestFumiliLMA.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b379321303c2b884c1b7a9a72bbd87d480550346
--- /dev/null
+++ b/App/src/TestFumiliLMA.cpp
@@ -0,0 +1,228 @@
+#include "TestFumiliLMA.h"
+#include "OutputData.h"
+#include "IsGISAXSTools.h"
+#include "MathFunctions.h"
+#include "ChiSquaredModule.h"
+#include "ISquaredFunction.h"
+
+#include <iostream>
+#include <cmath>
+#include <boost/bind.hpp>
+#include <boost/function.hpp>
+
+#include "TCanvas.h"
+#include "TF1.h"
+#include "TF2.h"
+#include "TH2D.h"
+#include "TRandom.h"
+
+
+
+/* ************************************************************************* */
+//
+/* ************************************************************************* */
+TestFumiliLMA::TestFumiliLMA()
+ : m_root_minimizer(0)
+ , m_fcn(0)
+ , m_func_object(0)
+ , m_func(0)
+ , m_chi_module(0)
+ , m_real_data(0)
+ , m_ndim(0)
+ , m_sigma(0.01)
+ , m_c1(0)
+{
+ // "simulation"
+ m_ndim = 3;
+ m_func_object = new SincXSincYFunctionObject();
+ m_func = new TF2("sincxy", m_func_object, -50.,50., -50.,50., m_ndim, "SincXSincYFunctionObject");
+ m_func->SetParameters(1.0, 1.0, 0.5); // parameters we have to find
+
+
+ // chi module
+ m_chi_module = new ChiSquaredModule();
+ m_chi_module->setChiSquaredFunction(SquaredFunctionWithGaussianError(m_sigma) );
+
+ // "real" data
+ OutputData<double> data;
+ FillOutputDataFromFunction(data, m_func);
+ m_real_data = createDataWithGaussianNoise(data, m_sigma);
+
+
+}
+
+
+TestFumiliLMA::~TestFumiliLMA()
+{
+ delete m_root_minimizer;
+ delete m_fcn;
+ //delete m_func_object;
+ //delete m_func;
+ delete m_chi_module;
+ delete m_real_data;
+}
+
+
+/* ************************************************************************* */
+//
+/* ************************************************************************* */
+void TestFumiliLMA::execute()
+{
+ std::cout << "TestFumiliLMA::execute() -> Hello world" << std::endl;
+ //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->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);
+// m_root_minimizer->SetVariable(0, "p0", 5.0, 0.01);
+// m_root_minimizer->SetVariable(1, "p1", -10.0, 0.01);
+// m_root_minimizer->SetVariable(2, "p2", 10.0, 0.01);
+
+
+ m_c1 = new TCanvas("c1","c1", 1024, 768);
+ m_c1->Divide(2,2);
+
+ m_c1->cd(1);
+ m_func->Draw("SURF");
+
+ TH2D *h2_real = IsGISAXSTools::getOutputDataTH2D(*m_real_data, "real_data");
+ m_c1->cd(2);
+ h2_real->DrawCopy("SURF");
+
+ // fitting function
+// boost::function<double(const double *)> fcn = boost::bind(&TestFumiliLMA::functionToMinimize, this, _1);
+// m_fcn = new ROOT::Math::Functor(fcn, m_ndim);
+// m_root_minimizer->SetFunction(*m_fcn);
+
+ // fitting function with gradient
+ MyChi2Function *chi = new MyChi2Function(this);
+ m_root_minimizer->SetFunction(*chi);
+
+ m_root_minimizer->SetPrintLevel(3);
+ m_root_minimizer->Minimize();
+ m_root_minimizer->PrintResults();
+
+ std::cout << "MyResults: " << m_root_minimizer->MinValue() << std::endl;
+ for(size_t i=0; i<m_root_minimizer->NFree(); ++i) {
+ std::cout << " i " << m_root_minimizer->X()[i]<< std::endl;
+ }
+
+// ROOT::Fit::Fitter fitter;
+// const int Npar = 3;
+// double par0[Npar] = { 5.0, -10.0, 10.0};
+
+// // create before the parameter settings in order to fix or set range on them
+// fitter.Config().SetParamsSettings(3,par0);
+// fitter.Config().ParSettings(0).SetStepSize(0.01);
+// fitter.Config().ParSettings(1).SetStepSize(0.01);
+// fitter.Config().ParSettings(2).SetStepSize(0.01);
+
+// fitter.Config().MinimizerOptions().SetPrintLevel(3);
+// fitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit");
+
+// MyChi2Function *chi = new MyChi2Function(this);
+//// m_root_minimizer->SetFunction(*chi);
+// fitter.FitFCN(*chi);
+// ROOT::Fit::FitResult result = fitter.Result();
+// result.Print(std::cout);
+
+}
+
+
+/* ************************************************************************* */
+//
+/* ************************************************************************* */
+double TestFumiliLMA::functionToMinimize(const double *pars)
+{
+ m_func->SetParameters(pars);
+ OutputData<double> simulated_data;
+ FillOutputDataFromFunction(simulated_data, m_func);
+ m_chi_module->setRealAndSimulatedData(*m_real_data, simulated_data);
+ double chi = m_chi_module->calculateChiSquared();
+
+ TH2D *h2_simul = IsGISAXSTools::getOutputDataTH2D(simulated_data, "simulated_data");
+ m_c1->cd(3);
+ h2_simul->DrawCopy("SURF");
+ delete h2_simul;
+ m_c1->Update();
+
+ return chi;
+}
+
+
+/* ************************************************************************* */
+//
+/* ************************************************************************* */
+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();
+ data.addAxis("x", nbinsx, fun->GetXmin(), fun->GetXmax());
+ data.addAxis("y", nbinsy, fun->GetYmin(), fun->GetYmax());
+ data.setAllTo(0.0);
+
+ OutputData<double>::iterator it=data.begin();
+ const AxisDouble *xaxis = data.getAxis(0);
+ const AxisDouble *yaxis = data.getAxis(1);
+ while( it!= data.end() )
+ {
+ size_t index_x = data.toCoordinates(it.getIndex())[0];
+ size_t index_y = data.toCoordinates(it.getIndex())[1];
+ double x = (*xaxis)[index_x];
+ double y = (*yaxis)[index_y];
+ double value = fun->Eval(x,y);
+ (*it)=value;
+ ++it;
+ }
+}
+
+
+double MyChi2Function::DataElement(const double *pars, unsigned int i, double *g ) const
+{
+ double xx[2];
+ std::cout << " DataElement() -> " << g << 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];
+ size_t index_y = m_test->m_real_data->toCoordinates(i)[1];
+ double x = (*xaxis)[index_x];
+ double y = (*yaxis)[index_y];
+ m_test->m_func->SetParameters(pars);
+ double value = m_test->m_func->Eval(x,y);
+ double real_value = (*m_test->m_real_data)[i];
+// std::cout << " DataElement " << i << " ix:" << index_x << " iy:" << index_y << " x:" << x << " y:" << y << " value:" << value << " real_val:" << real_value<< std::endl;
+ double weight = 1./(m_test->m_sigma * m_test->m_sigma);
+ double residual = weight * ( real_value - value ) ;
+ if (g != 0) {
+// std::cout << " gradient ! " << std::endl;
+ m_test->m_func->SetParameters(pars);
+ xx[0] = x;
+ xx[1] = y;
+ m_test->m_func->GradientPar( xx, g,1.E-9);
+ // need to mutiply by -1*weight
+ for (int ipar = 0; ipar < m_test->m_func->GetNpar(); ++ipar) {
+ g[ipar] *= -1 * weight;
+ }
+ }
+ return residual;
+}
+
+
+
diff --git a/App/src/TestIsGISAXS12.cpp b/App/src/TestIsGISAXS12.cpp
index fdac6ae83ecba87b9a8637f80a29e229ae20dc3b..a46c265841fa0c9f9fbe86d9c0986afcfb2c9b3e 100644
--- a/App/src/TestIsGISAXS12.cpp
+++ b/App/src/TestIsGISAXS12.cpp
@@ -152,12 +152,12 @@ void TestIsGISAXS12::run_isgisaxs_fit()
m_fitSuite = new FitSuite();
//m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
//m_fitSuite->setMinimizer( new ROOTMinimizer("GSLMultiFit", "") );
- m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Fumili") );
+ m_fitSuite->setMinimizer( new ROOTMinimizer("Minuit2", "Migrad") );
m_fitSuite->attachObserver( new FitSuiteObserverPrint(10) );
m_fitSuite->attachObserver( new FitSuiteObserverDraw(50) );
- ROOT::Math::Minimizer *minim = (dynamic_cast<ROOTMinimizer *>(m_fitSuite->getMinimizer()))->getROOTMinimizer();
- minim->SetStrategy(1);
+// ROOT::Math::Minimizer *minim = (dynamic_cast<ROOTMinimizer *>(m_fitSuite->getMinimizer()))->getROOTMinimizer();
+// minim->SetStrategy(1);
// minim->SetPrecision(1.);
m_fitSuite->addFitParameter("*Normalizer/scale", 1e5, 1, AttLimits::limited(1e4, 2e5));
diff --git a/Core/Algorithms/inc/ISquaredFunction.h b/Core/Algorithms/inc/ISquaredFunction.h
index a3d8235efca456499d6b268dc97e7352d8dca79b..9c2c28bd37a7b3cfb7696742e0d920f233df0570 100644
--- a/Core/Algorithms/inc/ISquaredFunction.h
+++ b/Core/Algorithms/inc/ISquaredFunction.h
@@ -90,4 +90,23 @@ private:
};
+class SquaredFunctionWithGaussianError : public ISquaredFunction
+{
+public:
+ SquaredFunctionWithGaussianError(double sigma = 0.01) : m_sigma(sigma){}
+ virtual ~SquaredFunctionWithGaussianError() {}
+ virtual SquaredFunctionWithGaussianError *clone() const { return new SquaredFunctionWithGaussianError(*this); }
+
+ virtual inline double calculateSquaredDifference(double real_value, double simulated_value) const
+ {
+ double diff_squared = (simulated_value-real_value)*(simulated_value-real_value);
+ double sigma_squared = m_sigma*m_sigma;
+ return diff_squared/sigma_squared;
+ }
+private:
+ double m_sigma;
+};
+
+
+
#endif /* ISQUAREDFUNCTION_H_ */
diff --git a/Core/Tools/inc/ParameterPool.h b/Core/Tools/inc/ParameterPool.h
index f2a3cb6df7bcddda6ec04d04c3974d1c4e53e748..e57f04715f800cae046e25998489ad82b41de56a 100644
--- a/Core/Tools/inc/ParameterPool.h
+++ b/Core/Tools/inc/ParameterPool.h
@@ -38,7 +38,7 @@ public:
void checkNull() const { if(!m_data) throw NullPointerException("ParameterPool::RealPar::getValue() -> Attempt to access uninitialised pointer."); }
friend std::ostream &operator<<(std::ostream &ostr, const RealPar &p) { ostr << p.m_data; return ostr; }
private:
- double *m_data;
+ volatile double *m_data;
};
//! definition of parameter type and parameter container