From f8bb6d63f708fb612610f0c3af5c47700d8bdaa5 Mon Sep 17 00:00:00 2001 From: Gennady Pospelov <g.pospelov@fz-juelich.de> Date: Tue, 27 Aug 2013 11:29:56 +0200 Subject: [PATCH] libMathMore is merged into libRootMinimizers, functional TestFit now are working under windows --- BornAgain.pro | 16 +- Core/Tools/inc/IObserver.h | 4 +- Core/Tools/inc/MessageService.h | 3 +- Core/Tools/inc/OutputDataFunctions.h | 2 +- Core/Tools/inc/Utils.h | 2 +- Fit/Factory/inc/FitSuite.h | 2 +- Fit/Factory/inc/FitSuiteParameters.h | 2 +- Fit/Factory/inc/IMinimizer.h | 9 +- Fit/Factory/inc/MinimizerFactory.h | 2 +- Fit/Factory/src/FitSuite.cpp | 5 + Fit/Factory/src/MinimizerFactory.cpp | 4 + Fit/Factory/src/ROOTMinimizer.cpp | 7 + Fit/Factory/src/ROOTMinimizerHelper.cpp | 5 +- Fit/Fit.pro | 7 +- Tests/FunctionalTests/TestCore/TestCore.pro | 22 +- .../TestFit/TestFit01/TestFit01.pro | 2 +- .../TestFit/TestFit02/TestFit02.pro | 2 +- ThirdParty/RootMinimizers/RootMinimizers.pro | 38 ++ .../RootMinimizers/inc/Math/Derivator.h | 251 ++++++++++ .../inc/Math/GSLFunctionAdapter.h | 97 ++++ .../inc/Math/GSLFunctionWrapper.h | 150 ++++++ .../RootMinimizers/inc/Math/GSLMinimizer.h | 231 +++++++++ .../RootMinimizers/inc/Math/GSLMinimizer1D.h | 224 +++++++++ .../RootMinimizers/inc/Math/GSLNLSMinimizer.h | 297 +++++++++++ .../RootMinimizers/inc/Math/GSLRndmEngines.h | 437 +++++++++++++++++ .../inc/Math/GSLSimAnMinimizer.h | 204 ++++++++ .../RootMinimizers/inc/Math/GSLSimAnnealing.h | 257 ++++++++++ .../inc/Math/MultiNumGradFunction.h | 143 ++++++ .../RootMinimizers/src/Math/Derivator.cxx | 161 ++++++ .../src/Math/GSL1DMinimizerWrapper.h | 76 +++ .../RootMinimizers/src/Math/GSLDerivator.cxx | 128 +++++ .../RootMinimizers/src/Math/GSLDerivator.h | 176 +++++++ .../src/Math/GSLFunctionWrapper.h | 150 ++++++ .../RootMinimizers/src/Math/GSLMinimizer.cxx | 400 +++++++++++++++ .../src/Math/GSLMinimizer1D.cxx | 222 +++++++++ .../RootMinimizers/src/Math/GSLMultiFit.h | 209 ++++++++ .../src/Math/GSLMultiFitFunctionAdapter.h | 131 +++++ .../src/Math/GSLMultiFitFunctionWrapper.h | 100 ++++ .../src/Math/GSLMultiMinFunctionAdapter.h | 99 ++++ .../src/Math/GSLMultiMinFunctionWrapper.h | 163 ++++++ .../src/Math/GSLMultiMinimizer.h | 213 ++++++++ .../src/Math/GSLNLSMinimizer.cxx | 463 ++++++++++++++++++ .../src/Math/GSLRndmEngines.cxx | 380 ++++++++++++++ .../RootMinimizers/src/Math/GSLRngWrapper.h | 140 ++++++ .../src/Math/GSLSimAnMinimizer.cxx | 230 +++++++++ .../src/Math/GSLSimAnnealing.cxx | 240 +++++++++ .../src/Math/MultiNumGradFunction.cxx | 61 +++ shared.pri | 2 +- 48 files changed, 6136 insertions(+), 33 deletions(-) create mode 100644 ThirdParty/RootMinimizers/inc/Math/Derivator.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLFunctionAdapter.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLFunctionWrapper.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLMinimizer.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLMinimizer1D.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLNLSMinimizer.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLSimAnMinimizer.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLSimAnnealing.h create mode 100644 ThirdParty/RootMinimizers/inc/Math/MultiNumGradFunction.h create mode 100644 ThirdParty/RootMinimizers/src/Math/Derivator.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/GSL1DMinimizerWrapper.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLDerivator.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLDerivator.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLFunctionWrapper.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMinimizer.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMinimizer1D.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionAdapter.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionWrapper.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionAdapter.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionWrapper.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLMultiMinimizer.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLNLSMinimizer.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLRngWrapper.h create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLSimAnMinimizer.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/GSLSimAnnealing.cxx create mode 100644 ThirdParty/RootMinimizers/src/Math/MultiNumGradFunction.cxx diff --git a/BornAgain.pro b/BornAgain.pro index e259719b7f1..c5394641d1e 100644 --- a/BornAgain.pro +++ b/BornAgain.pro @@ -5,22 +5,24 @@ include($$PWD/shared.pri) SUBDIRS += Core SUBDIRS += ThirdParty/gtest SUBDIRS += Tests/UnitTests/TestCore -TestCore.depends = ThirdParty/gtest +#TestCore.depends = ThirdParty/gtest SUBDIRS += Tests/FunctionalTests/TestCore -#isEmpty(ROOT_FRAMEWORK) { -# message("No ROOT installation found. Additional library libRootMinimizers.so will be compiled.") -# SUBDIRS += ThirdParty/RootMinimizers -#} +isEmpty(ROOT_FRAMEWORK) { + message("No ROOT installation found. Additional library libRootMinimizers.so will be compiled.") + SUBDIRS += ThirdParty/RootMinimizers +} #SUBDIRS += ThirdParty/RootMathMore -#SUBDIRS += Fit + +SUBDIRS += Fit +SUBDIRS += Tests/FunctionalTests/TestFit + #!isEmpty(ROOT_FRAMEWORK) { # SUBDIRS += App #} -#SUBDIRS += Tests/FunctionalTests/TestFit # compilation in lister order CONFIG += ordered diff --git a/Core/Tools/inc/IObserver.h b/Core/Tools/inc/IObserver.h index d908e6498e2..bc4bf7d048d 100644 --- a/Core/Tools/inc/IObserver.h +++ b/Core/Tools/inc/IObserver.h @@ -28,7 +28,7 @@ class IObservable; //! Observer interface from Observer pattern, for 1:n object dependencies. -class IObserver { +class BA_CORE_API_ IObserver { public: // IObserver() : m_observed_subject(0) {} virtual ~IObserver() {} @@ -44,7 +44,7 @@ class IObserver { //! Observable interface from Observer pattern, for 1:n object dependencies. -class IObservable { +class BA_CORE_API_ IObservable { public: typedef boost::shared_ptr<IObserver > observer_t; typedef std::list<observer_t > observerlist_t; diff --git a/Core/Tools/inc/MessageService.h b/Core/Tools/inc/MessageService.h index 7330f1be9ee..ce89b451ba3 100644 --- a/Core/Tools/inc/MessageService.h +++ b/Core/Tools/inc/MessageService.h @@ -1,6 +1,7 @@ #ifndef MESSAGESVC_H #define MESSAGESVC_H +#include "WinDllMacros.h" #include <iostream> #include <sstream> #include <string> @@ -15,7 +16,7 @@ namespace MSG enum MessageLevel { VERBOSE, DEBUG, INFO, WARNING, ERROR, FATAL }; -class Logger +class BA_CORE_API_ Logger { public: Logger(MessageLevel level) { diff --git a/Core/Tools/inc/OutputDataFunctions.h b/Core/Tools/inc/OutputDataFunctions.h index e74f358d769..4940312f9d0 100644 --- a/Core/Tools/inc/OutputDataFunctions.h +++ b/Core/Tools/inc/OutputDataFunctions.h @@ -27,7 +27,7 @@ namespace OutputDataFunctions { //! double the bin size for each dimension - OutputData<double> *doubleBinSize(const OutputData<double>& source); + BA_CORE_API_ OutputData<double> *doubleBinSize(const OutputData<double>& source); //! unnormalized Fourier transformation for real data void FourierTransform( diff --git a/Core/Tools/inc/Utils.h b/Core/Tools/inc/Utils.h index cd9c6463382..11c1eb02ad4 100644 --- a/Core/Tools/inc/Utils.h +++ b/Core/Tools/inc/Utils.h @@ -31,7 +31,7 @@ namespace Utils { //! Collection of utilities for std::string. -class String +class BA_CORE_API_ String { public: //! Parse double values from string to vector of double. diff --git a/Fit/Factory/inc/FitSuite.h b/Fit/Factory/inc/FitSuite.h index 2d08f1b3146..4b5451326e5 100644 --- a/Fit/Factory/inc/FitSuite.h +++ b/Fit/Factory/inc/FitSuite.h @@ -30,7 +30,7 @@ class ParameterPool; //! Main class to perform fitting -class FitSuite : public IObservable +class BA_CORE_API_ FitSuite : public IObservable { public: FitSuite(); diff --git a/Fit/Factory/inc/FitSuiteParameters.h b/Fit/Factory/inc/FitSuiteParameters.h index f00fbe9d980..688aaf2b087 100644 --- a/Fit/Factory/inc/FitSuiteParameters.h +++ b/Fit/Factory/inc/FitSuiteParameters.h @@ -25,7 +25,7 @@ class ParameterPool; //! Holds vector of parameters for FitSuite -class FitSuiteParameters +class BA_CORE_API_ FitSuiteParameters { public: //typedef SafePointerVector<FitParameter > parameters_t; diff --git a/Fit/Factory/inc/IMinimizer.h b/Fit/Factory/inc/IMinimizer.h index 937d4f71c60..2da417926c8 100644 --- a/Fit/Factory/inc/IMinimizer.h +++ b/Fit/Factory/inc/IMinimizer.h @@ -16,7 +16,14 @@ #ifndef IMINIMIZER_H #define IMINIMIZER_H +#include "WinDllMacros.h" + +#include "Macros.h" +GCC_DIAG_OFF(unused-local-typedefs); #include <boost/function.hpp> +GCC_DIAG_ON(unused-local-typedefs); + + #include "Exceptions.h" #include <vector> @@ -25,7 +32,7 @@ class FitSuiteParameters; //! Common interface for all kind minimizer's -class IMinimizer +class BA_CORE_API_ IMinimizer { public: //! signature of chi squared function to minimize diff --git a/Fit/Factory/inc/MinimizerFactory.h b/Fit/Factory/inc/MinimizerFactory.h index b11d9b70643..9ff3e61ca6d 100644 --- a/Fit/Factory/inc/MinimizerFactory.h +++ b/Fit/Factory/inc/MinimizerFactory.h @@ -24,7 +24,7 @@ //! @class MinimizerFactory //! @brief Factory to create minimizers //- ------------------------------------------------------------------- -class MinimizerFactory +class BA_CORE_API_ MinimizerFactory { public: static IMinimizer *createMinimizer(const std::string& minimizer, const std::string& algorithm = std::string(), const std::string& options=std::string() ); diff --git a/Fit/Factory/src/FitSuite.cpp b/Fit/Factory/src/FitSuite.cpp index 80d361d97e5..2de7c91aba7 100644 --- a/Fit/Factory/src/FitSuite.cpp +++ b/Fit/Factory/src/FitSuite.cpp @@ -19,7 +19,12 @@ #include "MessageService.h" #include "FitSuitePrintObserver.h" + +#include "Macros.h" +GCC_DIAG_OFF(unused-local-typedefs); #include <boost/bind.hpp> +GCC_DIAG_ON(unused-local-typedefs); + FitSuite::FitSuite() : m_minimizer(0), m_is_last_iteration(false) { diff --git a/Fit/Factory/src/MinimizerFactory.cpp b/Fit/Factory/src/MinimizerFactory.cpp index 1e9160fea9b..f82e9b72bfa 100644 --- a/Fit/Factory/src/MinimizerFactory.cpp +++ b/Fit/Factory/src/MinimizerFactory.cpp @@ -5,7 +5,11 @@ #include "MinimizerScan.h" #include "ROOTMinimizer.h" +#include "Macros.h" +GCC_DIAG_OFF(unused-local-typedefs); #include <boost/assign/list_of.hpp> +GCC_DIAG_ON(unused-local-typedefs); + #include <iomanip> MinimizerFactory::Catalogue MinimizerFactory::m_catalogue = MinimizerFactory::Catalogue(); diff --git a/Fit/Factory/src/ROOTMinimizer.cpp b/Fit/Factory/src/ROOTMinimizer.cpp index bc709dce398..efb2ded2b75 100644 --- a/Fit/Factory/src/ROOTMinimizer.cpp +++ b/Fit/Factory/src/ROOTMinimizer.cpp @@ -3,8 +3,15 @@ #include "FitSuiteParameters.h" #include "Utils.h" #include "ROOTMinimizerFunction.h" + #include <iomanip> #include <sstream> + +#include "Macros.h" +GCC_DIAG_OFF(unused-local-typedefs); +#include <boost/assign/list_of.hpp> +GCC_DIAG_ON(unused-local-typedefs); + #include <boost/assign/list_of.hpp> #include "ROOTGSLNLSMinimizer.h" #include "ROOTGSLSimAnMinimizer.h" diff --git a/Fit/Factory/src/ROOTMinimizerHelper.cpp b/Fit/Factory/src/ROOTMinimizerHelper.cpp index fac5011e0e4..e8eb806bfcb 100644 --- a/Fit/Factory/src/ROOTMinimizerHelper.cpp +++ b/Fit/Factory/src/ROOTMinimizerHelper.cpp @@ -1,8 +1,11 @@ #include "ROOTMinimizerHelper.h" #include "ROOTGSLSimAnMinimizer.h" #include "Utils.h" -#include <boost/lexical_cast.hpp> +#include "Macros.h" +GCC_DIAG_OFF(unused-local-typedefs); +#include <boost/lexical_cast.hpp> +GCC_DIAG_ON(unused-local-typedefs); // ---------------------------------------------------------------------------- // translate text with options into appropriate calls of minimizer's set method diff --git a/Fit/Fit.pro b/Fit/Fit.pro index c45f3a2f0e1..8036262a78b 100644 --- a/Fit/Fit.pro +++ b/Fit/Fit.pro @@ -59,6 +59,9 @@ contains(CONFIG, BORNAGAIN_PYTHON) { include($$PWD/python_module.pri) } +win32 { + DEFINES += BA_CORE_BUILD_DLL +} # ----------------------------------------------------------------------------- # includes @@ -66,7 +69,7 @@ contains(CONFIG, BORNAGAIN_PYTHON) { INCLUDEPATH += $$PWD/Factory/inc DEPENDPATH += $$PWD/Factory/inc -INCLUDEPATH += $${RootMathMore_INCLUDEPATH} +#INCLUDEPATH += $${RootMathMore_INCLUDEPATH} isEmpty(ROOT_FRAMEWORK) { INCLUDEPATH += $${RootMinimizers_INCLUDEPATH} } else { @@ -78,7 +81,7 @@ isEmpty(ROOT_FRAMEWORK) { # additional libraries # ----------------------------------------------------------------------------- LIBS += $$PWD/../lib/libBornAgainCore.$${SONAME} -LIBS += $${RootMathMore_LIB} +#LIBS += $${RootMathMore_LIB} isEmpty(ROOT_FRAMEWORK) { LIBS += $${RootMinimizers_LIB} } else { diff --git a/Tests/FunctionalTests/TestCore/TestCore.pro b/Tests/FunctionalTests/TestCore/TestCore.pro index dee3bb7551f..2d3a3e29f2b 100644 --- a/Tests/FunctionalTests/TestCore/TestCore.pro +++ b/Tests/FunctionalTests/TestCore/TestCore.pro @@ -2,16 +2,16 @@ TEMPLATE = subdirs SUBDIRS += \ IsGISAXS01 \ - IsGISAXS02 \ - IsGISAXS03 \ - IsGISAXS04 \ - IsGISAXS06 \ - IsGISAXS07 \ - IsGISAXS08 \ - IsGISAXS09 \ - IsGISAXS10 \ - IsGISAXS11 \ - IsGISAXS15 \ - MesoCrystal1 +# IsGISAXS02 \ +# IsGISAXS03 \ +# IsGISAXS04 \ +# IsGISAXS06 \ +# IsGISAXS07 \ +# IsGISAXS08 \ +# IsGISAXS09 \ +# IsGISAXS10 \ +# IsGISAXS11 \ +# IsGISAXS15 \ +# MesoCrystal1 CONFIG += ordered diff --git a/Tests/FunctionalTests/TestFit/TestFit01/TestFit01.pro b/Tests/FunctionalTests/TestFit/TestFit01/TestFit01.pro index d522a670f38..f9b8e049d60 100644 --- a/Tests/FunctionalTests/TestFit/TestFit01/TestFit01.pro +++ b/Tests/FunctionalTests/TestFit/TestFit01/TestFit01.pro @@ -11,7 +11,7 @@ INCLUDEPATH += $$myIncludes DEPENDPATH += $$myIncludes DEFINES += STANDALONE -LIBS += $$BornAgainCore_LIB $$BornAgainFit_LIB $$RootMathMore_LIB +LIBS += $$BornAgainCore_LIB $$BornAgainFit_LIB isEmpty(ROOT_FRAMEWORK) { LIBS += $$RootMinimizers_LIB } diff --git a/Tests/FunctionalTests/TestFit/TestFit02/TestFit02.pro b/Tests/FunctionalTests/TestFit/TestFit02/TestFit02.pro index 0d58e65ccaa..177f9d9e6c7 100644 --- a/Tests/FunctionalTests/TestFit/TestFit02/TestFit02.pro +++ b/Tests/FunctionalTests/TestFit/TestFit02/TestFit02.pro @@ -11,7 +11,7 @@ INCLUDEPATH += $$myIncludes DEPENDPATH += $$myIncludes DEFINES += STANDALONE -LIBS += $$BornAgainCore_LIB $$BornAgainFit_LIB $$RootMathMore_LIB +LIBS += $$BornAgainCore_LIB $$BornAgainFit_LIB isEmpty(ROOT_FRAMEWORK) { LIBS += $$RootMinimizers_LIB } diff --git a/ThirdParty/RootMinimizers/RootMinimizers.pro b/ThirdParty/RootMinimizers/RootMinimizers.pro index 08d5d711c5a..0fda8d11c21 100644 --- a/ThirdParty/RootMinimizers/RootMinimizers.pro +++ b/ThirdParty/RootMinimizers/RootMinimizers.pro @@ -84,6 +84,17 @@ SOURCES += \ src/Minuit2/SqrtUpParameterTransformation.cxx \ src/Minuit2/VariableMetricBuilder.cxx \ src/Minuit2/VariableMetricEDMEstimator.cxx \ + \ # from MathMore + src/Math/Derivator.cxx \ + src/Math/GSLDerivator.cxx \ + src/Math/GSLMinimizer1D.cxx \ + src/Math/GSLMinimizer.cxx \ + src/Math/GSLRndmEngines.cxx \ + src/Math/GSLSimAnnealing.cxx \ + src/Math/MultiNumGradFunction.cxx \ + src/Math/GSLNLSMinimizer.cxx \ + src/Math/GSLSimAnMinimizer.cxx \ + HEADERS += \ inc/Fit/BinData.h \ @@ -235,12 +246,39 @@ HEADERS += \ inc/Minuit2/VariableMetricEDMEstimator.h \ inc/Minuit2/VariableMetricMinimizer.h \ inc/Minuit2/VectorOuterProduct.h \ + \ # from MathMore + inc/Math/Derivator.h \ + inc/Math/GSLFunctionAdapter.h \ + inc/Math/GSLFunctionWrapper.h \ + inc/Math/GSLMinimizer1D.h \ + inc/Math/GSLMinimizer.h \ + inc/Math/GSLRndmEngines.h \ + inc/Math/GSLSimAnnealing.h \ + inc/Math/MultiNumGradFunction.h \ + inc/Math/GSLNLSMinimizer.h \ + inc/Math/GSLSimAnMinimizer.h \ + \ + src/Math/GSL1DMinimizerWrapper.h \ + src/Math/GSLDerivator.h \ + src/Math/GSLFunctionWrapper.h \ + src/Math/GSLMultiFitFunctionAdapter.h \ + src/Math/GSLMultiFitFunctionWrapper.h \ + src/Math/GSLMultiFit.h \ + src/Math/GSLMultiMinFunctionAdapter.h \ + src/Math/GSLMultiMinFunctionWrapper.h \ + src/Math/GSLMultiMinimizer.h \ + src/Math/GSLRngWrapper.h \ + INCLUDEPATH += inc $${RootMathMore_INCLUDEPATH} DEPENDPATH += inc $${RootMathMore_INCLUDEPATH} QMAKE_CXXFLAGS += -DMATH_NO_PLUGIN_MANAGER -DHAS_MINUIT2 -DR__HAS_MATHMORE +win32 { + DEFINES += BA_CORE_BUILD_DLL +} + # ----------------------------------------------------------------------------- # Installing library into dedicated directory at the end of compilation # ----------------------------------------------------------------------------- diff --git a/ThirdParty/RootMinimizers/inc/Math/Derivator.h b/ThirdParty/RootMinimizers/inc/Math/Derivator.h new file mode 100644 index 00000000000..0b5a0c65cca --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/Derivator.h @@ -0,0 +1,251 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class Derivator +// +// class for calculating Derivative of functions +// +// Created by: moneta at Sat Nov 13 14:46:00 2004 +// +// Last update: Sat Nov 13 14:46:00 2004 +// +#ifndef ROOT_Math_Derivator +#define ROOT_Math_Derivator + +/** + @defgroup Deriv Numerical Differentiation + @ingroup NumAlgo +*/ + +#ifndef ROOT_Math_IFunctionfwd +#include "Math/IFunctionfwd.h" +#endif + +#ifndef ROOT_Math_IParamFunctionfwd +#include "Math/IParamFunctionfwd.h" +#endif + + +namespace ROOT { +namespace Math { + + + +class GSLDerivator; + +//_______________________________________________________________________ +/** + Class for computing numerical derivative of a function. + Presently this class is implemented only using the numerical derivatives + algorithms provided by GSL + using the implementation class ROOT::Math::GSLDerivator + + This class does not support copying + + @ingroup Deriv +*/ + +class Derivator { + +public: + + /** + signature for function pointers used by GSL + */ + typedef double ( * GSLFuncPointer ) ( double, void * ); + + /** + Empty Construct for a Derivator class + Need to set the function afterwards with Derivator::SetFunction + */ + Derivator(); + /** + Construct using a ROOT::Math::IGenFunction interface + */ + explicit Derivator(const IGenFunction &f); + /** + Construct using a GSL function pointer type + @param f : free function pointer of the GSL required type + @param p : pointer to the object carrying the function state + (for example the function object itself) + */ + explicit Derivator(const GSLFuncPointer &f, void * p = 0); + + /// destructor + virtual ~Derivator(); + + // disable copying +private: + + Derivator(const Derivator &); + Derivator & operator = (const Derivator &); + +public: + + +#ifdef LATER + /** + Template methods for generic functions + Set the function f for evaluating the derivative. + The function type must implement the assigment operator, + <em> double operator() ( double x ) </em> + */ + template <class UserFunc> + inline void SetFunction(const UserFunc &f) { + const void * p = &f; + SetFunction( &GSLFunctionAdapter<UserFunc>::F, const_cast<void *>(p) ); + } +#endif + + /** + Set the function for calculating the derivatives. + The function must implement the ROOT::Math::IGenFunction signature + */ + void SetFunction(const IGenFunction &f); + + + /** + Set the function f for evaluating the derivative using a GSL function pointer type + @param f : free function pointer of the GSL required type + @param p : pointer to the object carrying the function state + (for example the function object itself) + */ + void SetFunction( const GSLFuncPointer &f, void * p = 0); + + + + /** + Computes the numerical derivative of a function f at a point x. + It uses Derivator::EvalCentral to compute the derivative using an + adaptive central difference algorithm with a step size h + */ + + double Eval(double x, double h = 1E-8) const; + + + + /** + Computes the numerical derivative at a point x using an adaptive central + difference algorithm with a step size h. + */ + double EvalCentral( double x, double h = 1E-8) const; + + /** + Computes the numerical derivative at a point x using an adaptive forward + difference algorithm with a step size h. + The function is evaluated only at points greater than x and at x itself. + */ + double EvalForward( double x, double h = 1E-8) const; + + /** + Computes the numerical derivative at a point x using an adaptive backward + difference algorithm with a step size h. + The function is evaluated only at points less than x and at x itself. + */ + double EvalBackward( double x, double h = 1E-8) const; + + /** @name --- Static methods --- + This methods don't require to use a Derivator object, and are designed to be used in + fast calculation. Error and status code cannot be retrieved in this case + */ + + /** + Computes the numerical derivative of a function f at a point x. + It uses Derivator::EvalCentral to compute the derivative using an + adaptive central difference algorithm with a step size h + */ + static double Eval(const IGenFunction & f, double x, double h = 1E-8); + + /** + Computes the numerical derivative of a function f at a point x using an adaptive central + difference algorithm with a step size h + */ + static double EvalCentral(const IGenFunction & f, double x, double h = 1E-8); + + + /** + Computes the numerical derivative of a function f at a point x using an adaptive forward + difference algorithm with a step size h. + The function is evaluated only at points greater than x and at x itself + */ + static double EvalForward(const IGenFunction & f, double x, double h = 1E-8); + + /** + Computes the numerical derivative of a function f at a point x using an adaptive backward + difference algorithm with a step size h. + The function is evaluated only at points less than x and at x itself + */ + static double EvalBackward(const IGenFunction & f, double x, double h = 1E-8); + + // Derivatives for multi-dimension functions + /** + Evaluate the partial derivative of a multi-dim function + with respect coordinate x_icoord at the point x[] + */ + static double Eval(const IMultiGenFunction & f, const double * x, unsigned int icoord = 0, double h = 1E-8); + + /** + Evaluate the derivative with respect a parameter for one-dim parameteric function + at the point ( x,p[]) with respect the parameter p_ipar + */ + static double Eval(IParamFunction & f, double x, const double * p, unsigned int ipar = 0, double h = 1E-8); + + /** + Evaluate the derivative with respect a parameter for a multi-dim parameteric function + at the point ( x[],p[]) with respect the parameter p_ipar + */ + static double Eval(IParamMultiFunction & f, const double * x, const double * p, unsigned int ipar = 0, double h = 1E-8); + + + /** + return the error status of the last derivative calculation + */ + int Status() const; + + /** + return the result of the last derivative calculation + */ + double Result() const; + + /** + return the estimate of the absolute error of the last derivative calculation + */ + double Error() const; + + +private: + + + mutable GSLDerivator * fDerivator; + +}; + + + + +} // namespace Math +} // namespace ROOT + + +#endif /* ROOT_Math_Derivator */ diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLFunctionAdapter.h b/ThirdParty/RootMinimizers/inc/Math/GSLFunctionAdapter.h new file mode 100644 index 00000000000..0e1b85e4acf --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLFunctionAdapter.h @@ -0,0 +1,97 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLFunctionAdapter +// +// Generic adapter for gsl_function signature +// usable for any c++ class which defines operator( ) +// +// Created by: Lorenzo Moneta at Fri Nov 12 16:58:51 2004 +// +// Last update: Fri Nov 12 16:58:51 2004 +// +#ifndef ROOT_Math_GSLFunctionAdapter +#define ROOT_Math_GSLFunctionAdapter + + +namespace ROOT { +namespace Math { + + /** + Function pointer corresponding to gsl_function signature + */ + + typedef double ( * GSLFuncPointer ) ( double, void *); + + + /** + Class for adapting any C++ functor class to C function pointers used by GSL. + The templated C++ function class must implement: + + <em> double operator( double x)</em> + and if the derivatives are required: + <em> double Gradient( double x)</em> + + This class defines static methods with will be used to fill the + \a gsl_function and \a gsl_function_fdf structs used by GSL. + See for examples the <A HREF="http://www.gnu.org/software/gsl/manual/gsl-ref_32.html#SEC432">GSL online manual</A> + */ + + + template<class UserFunc> + class GSLFunctionAdapter { + + public: + + GSLFunctionAdapter() {} + virtual ~GSLFunctionAdapter() {} + + static double F( double x, void * p) { + + UserFunc * function = reinterpret_cast< UserFunc *> (p); + return (*function)( x ); + } + + + static double Df( double x, void * p) { + + UserFunc * function = reinterpret_cast< UserFunc *> (p); + return (*function).Derivative( x ); + } + + static void Fdf( double x, void * p, double *f, double *df ) { + + UserFunc * function = reinterpret_cast< UserFunc *> (p); + *f = (*function) ( x ); + *df = (*function).Derivative( x ); + } + + }; + + +} // namespace Math +} // namespace ROOT + + +#endif /* ROOT_Math_GSLFunctionAdapter */ diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLFunctionWrapper.h b/ThirdParty/RootMinimizers/inc/Math/GSLFunctionWrapper.h new file mode 100644 index 00000000000..0dc09054801 --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLFunctionWrapper.h @@ -0,0 +1,150 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLFunctionWrapper +// +// Created by: moneta at Sat Nov 13 14:54:41 2004 +// +// Last update: Sat Nov 13 14:54:41 2004 +// +#ifndef ROOT_Math_GSLFunctionWrapper +#define ROOT_Math_GSLFunctionWrapper + +#include "gsl/gsl_math.h" + +#include "Math/GSLFunctionAdapter.h" + +#include <cassert> + +namespace ROOT { +namespace Math { + + + +typedef double ( * GSLFuncPointer ) ( double, void *); +typedef void ( * GSLFdfPointer ) ( double, void *, double *, double *); + + +/** + Wrapper class to the gsl_function C structure. + This class to fill the GSL C structure gsl_function with + the C++ function objcet. + Use the class ROOT::Math::GSLFunctionAdapter to adapt the + C++ function object to the right signature (function pointer type) + requested by GSL +*/ +class GSLFunctionWrapper { + +public: + + GSLFunctionWrapper() + { + fFunc.function = 0; + fFunc.params = 0; + } + + /// set in the GSL C struct the pointer to the function evaluation + void SetFuncPointer( GSLFuncPointer f) { fFunc.function = f; } + + /// set in the GSL C struct the extra-object pointer + void SetParams ( void * p) { fFunc.params = p; } + + /// fill the GSL C struct from a generic C++ callable object + /// implementing operator() + template<class FuncType> + void SetFunction(const FuncType &f) { + const void * p = &f; + assert (p != 0); + SetFuncPointer(&GSLFunctionAdapter<FuncType >::F); + SetParams(const_cast<void *>(p)); + } + + gsl_function * GetFunc() { return &fFunc; } + + GSLFuncPointer FunctionPtr() { return fFunc.function; } + + // evaluate the function + double operator() (double x) { return GSL_FN_EVAL(&fFunc, x); } + + /// check if function is valid (has been set) + bool IsValid() { + return (fFunc.function != 0) ? true : false; + } + +private: + gsl_function fFunc; + + +}; + + + /** + class to wrap a gsl_function_fdf (with derivatives) + */ + class GSLFunctionDerivWrapper { + + public: + + GSLFunctionDerivWrapper() + { + fFunc.f = 0; + fFunc.df = 0; + fFunc.fdf = 0; + fFunc.params = 0; + } + + + void SetFuncPointer( GSLFuncPointer f) { fFunc.f = f; } + void SetDerivPointer( GSLFuncPointer f) { fFunc.df = f; } + void SetFdfPointer( GSLFdfPointer f) { fFunc.fdf = f; } + void SetParams ( void * p) { fFunc.params = p; } + + + gsl_function_fdf * GetFunc() { return &fFunc; } + + // evaluate the function and derivatives + double operator() (double x) { return GSL_FN_FDF_EVAL_F(&fFunc, x); } + + double Derivative (double x) { return GSL_FN_FDF_EVAL_DF(&fFunc, x); } + + void Fdf(double x, double & f, double & df) { + return GSL_FN_FDF_EVAL_F_DF(&fFunc, x, &f, &df); + } + + /// check if function is valid (has been set) + bool IsValid() { + return (fFunc.f != 0 ) ? true : false; + } + + private: + gsl_function_fdf fFunc; + + }; + + + +} // namespace Math +} // namespace ROOT + +#endif /* ROOT_Math_GSLFunctionWrapper */ diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLMinimizer.h b/ThirdParty/RootMinimizers/inc/Math/GSLMinimizer.h new file mode 100644 index 00000000000..b0aee8051ff --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLMinimizer.h @@ -0,0 +1,231 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Oct 18 11:48:00 2006 + + /********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + + +// Header file for class GSLMinimizer + +#ifndef ROOT_Math_GSLMinimizer +#define ROOT_Math_GSLMinimizer + +#ifndef ROOT_Math_Minimizer +#include "Math/Minimizer.h" +#endif + + +#ifndef ROOT_Math_IFunctionfwd +#include "Math/IFunctionfwd.h" +#endif + +#ifndef ROOT_Math_IParamFunctionfwd +#include "Math/IParamFunctionfwd.h" +#endif + +#ifndef ROOT_Math_MinimizerVariable +#include "Math/MinimizerVariable.h" +#endif + + +#include <vector> +#include <map> +#include <string> + + + +namespace ROOT { + +namespace Math { + + + /** + enumeration specifying the types of GSL minimizers + @ingroup MultiMin + */ + enum EGSLMinimizerType { + kConjugateFR, + kConjugatePR, + kVectorBFGS, + kVectorBFGS2, + kSteepestDescent + }; + + + class GSLMultiMinimizer; + + class MinimTransformFunction; + + +//_____________________________________________________________________________________ +/** + GSLMinimizer class. + Implementation of the ROOT::Math::Minimizer interface using the GSL multi-dimensional + minimization algorithms. + + See <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Multidimensional-Minimization.html">GSL doc</A> + from more info on the GSL minimization algorithms. + + The class implements the ROOT::Math::Minimizer interface and can be instantiated using the + ROOT plugin manager (plugin name is "GSLMultiMin"). The varius minimization algorithms + (conjugatefr, conjugatepr, bfgs, etc..) can be passed as enumerations and also as a string. + The default algorithm is conjugatefr (Fletcher-Reeves conjugate gradient algorithm). + + @ingroup MultiMin +*/ +class GSLMinimizer : public ROOT::Math::Minimizer { + +public: + + /** + Default constructor + */ + GSLMinimizer (ROOT::Math::EGSLMinimizerType type = ROOT::Math::kConjugateFR ); + + /** + Constructor with a string giving name of algorithm + */ + GSLMinimizer (const char * type ); + + /** + Destructor + */ + virtual ~GSLMinimizer (); + +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + GSLMinimizer(const GSLMinimizer &) : Minimizer() {} + + /** + Assignment operator + */ + GSLMinimizer & operator = (const GSLMinimizer & rhs) { + if (this == &rhs) return *this; // time saving self-test + return *this; + } + +public: + + /// set the function to minimize + virtual void SetFunction(const ROOT::Math::IMultiGenFunction & func); + + /// set gradient the function to minimize + virtual void SetFunction(const ROOT::Math::IMultiGradFunction & func); + + /// set free variable + virtual bool SetVariable(unsigned int ivar, const std::string & name, double val, double step); + + + /// set lower limit variable (override if minimizer supports them ) + virtual bool SetLowerLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower ); + /// set upper limit variable (override if minimizer supports them ) + virtual bool SetUpperLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double upper ); + /// set upper/lower limited variable (override if minimizer supports them ) + virtual bool SetLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double /* lower */, double /* upper */); + /// set fixed variable (override if minimizer supports them ) + virtual bool SetFixedVariable(unsigned int /* ivar */, const std::string & /* name */, double /* val */); + /// set the value of an existing variable + virtual bool SetVariableValue(unsigned int ivar, double val ); + /// set the values of all existing variables (array must be dimensioned to the size of existing parameters) + virtual bool SetVariableValues(const double * x); + + /// method to perform the minimization + virtual bool Minimize(); + + /// return minimum function value + virtual double MinValue() const { return fMinVal; } + + /// return expected distance reached from the minimum + virtual double Edm() const { return 0; } // not impl. } + + /// return pointer to X values at the minimum + virtual const double * X() const { return &fValues.front(); } + + /// return pointer to gradient values at the minimum + virtual const double * MinGradient() const; + + /// number of function calls to reach the minimum + virtual unsigned int NCalls() const; + + /// this is <= Function().NDim() which is the total + /// number of variables (free+ constrained ones) + virtual unsigned int NDim() const { return fValues.size(); } + + /// number of free variables (real dimension of the problem) + /// this is <= Function().NDim() which is the total number of free parameters + virtual unsigned int NFree() const { return fObjFunc->NDim(); } + + /// minimizer provides error and error matrix + virtual bool ProvidesError() const { return false; } + + /// return errors at the minimum + virtual const double * Errors() const { + return 0; + } + + /** return covariance matrices elements + if the variable is fixed the matrix is zero + The ordering of the variables is the same as in errors + */ + virtual double CovMatrix(unsigned int , unsigned int ) const { return 0; } + + + // method of only GSL minimizer (not inherited from interface) + + /// return pointer to used objective gradient function + const ROOT::Math::IMultiGradFunction * ObjFunction() const { return fObjFunc; } + + /// return transformation function (NULL if not having a transformation) + const ROOT::Math::MinimTransformFunction * TransformFunction() const; + + +protected: + +private: + + // dimension of the function to be minimized + unsigned int fDim; + + ROOT::Math::GSLMultiMinimizer * fGSLMultiMin; + const ROOT::Math::IMultiGradFunction * fObjFunc; + + double fMinVal; + double fLSTolerance; // Line Search Tolerance + std::vector<double> fValues; + //mutable std::vector<double> fErrors; + std::vector<double> fSteps; + std::vector<std::string> fNames; + std::vector<ROOT::Math::EMinimVariableType> fVarTypes; // vector specifyng the type of variables + std::map< unsigned int, std::pair<double, double> > fBounds; // map specifying the bound using as key the parameter index + +}; + + } // end namespace Fit + +} // end namespace ROOT + + + +#endif /* ROOT_Math_GSLMinimizer */ diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLMinimizer1D.h b/ThirdParty/RootMinimizers/inc/Math/GSLMinimizer1D.h new file mode 100644 index 00000000000..3aa57443541 --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLMinimizer1D.h @@ -0,0 +1,224 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta, A. Zsenei 08/2005 + /********************************************************************** + * * + * Copyright (c) 2004 moneta, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLMinimizer1D +// +// Created by: moneta at Wed Dec 1 15:04:51 2004 +// +// Last update: Wed Dec 1 15:04:51 2004 +// + +#ifndef ROOT_Math_GSLMinimizer1D +#define ROOT_Math_GSLMinimizer1D + +#include "Math/IMinimizer1D.h" +#include "Math/GSLFunctionAdapter.h" + + +namespace ROOT { +namespace Math { + + namespace Minim1D { + + /** + Enumeration with One Dimensional Minimizer Algorithms. + The algorithms are implemented using GSL, see the + <A HREF="http://www.gnu.org/software/gsl/manual/gsl-ref_33.html#SEC447">GSL manual</A>. + + The algorithms available are: + <ul> + <li><em>Golden Section Algorithm</em>, simplest method of bracketing the minimum of a function + <li><em>Brent Algorithm</em>, which combines a parabolic interpolation with the golden section algorithm + </ul> + @ingroup Min1D + */ + + enum Type {kGOLDENSECTION, + kBRENT + }; + } + + class GSL1DMinimizerWrapper; + class GSLFunctionWrapper; + +//______________________________________________________________________________________ +/** + +Minimizer for arbitrary one dimensional functions. + +Implemented using GSL, for detailed description see: +<A HREF="http://www.gnu.org/software/gsl/manual/html_node/One-dimensional-Minimization.html">GSL online doc</A> + +The algorithms uspported are only bracketing algorithm which do not use derivatives information. +The algorithms which can be chosen at construction time are GOLDENSECTION, whic is the simplest method +but the slowest and BRENT (the default one) which combines the golden section with a parabolic interpolation. + + +This class does not support copying +@ingroup Min1D +*/ + + class GSLMinimizer1D: public IMinimizer1D { + + public: + + /** + Construct the minimizer passing the minimizer type using the Minim1D::Algorithm enumeration + */ + + explicit GSLMinimizer1D(Minim1D::Type type=Minim1D::kBRENT); + + /** + Destructor: free allocated resources + */ + virtual ~GSLMinimizer1D(); + + private: + // usually copying is non trivial, so we make this unaccessible + GSLMinimizer1D(const GSLMinimizer1D &); + GSLMinimizer1D & operator = (const GSLMinimizer1D &); + + public: + + + /** + Set, or reset, minimizer to use the function f and the initial search interval [xlow, xup], with a guess for the location of the minimum xmin. + The condition : \f$ f(xlow) > f(xmin) < f(xup)\f$ must be satisfied + */ + template <class UserFunc> + void SetFunction( const UserFunc & f, double xmin, double xlow, double xup) { + const void * p = &f; + SetFunction( &GSLFunctionAdapter<UserFunc>::F, const_cast<void *>(p), xmin, xlow, xup ); + } + + /** + Set, or reset, minimizer to use the function f and the initial search interval [xlow, xup], with a guess for the location of the minimum xmin. + The condition : \f$ f(xlow) > f(xmin) < f(xup) \f$ must be satisfied + + Method specialized on the GSL function type + */ + void SetFunction( GSLFuncPointer f, void * params, double xmin, double xlow, double xup); + + /** + Perform a minimizer iteration and + if an unexepcted problem occurr then an error code will be returned + */ + int Iterate(); + + + /** + Return current estimate of the position of the minimum + */ + double XMinimum() const; + + /** + Return current lower bound of the minimization interval + */ + double XLower() const; + + /** + Return current upper bound of the minimization interval + */ + double XUpper() const; + + /** + Return function value at current estimate of the minimum + */ + double FValMinimum() const; + + /** + Return function value at current lower bound of the minimization interval + */ + double FValLower() const; + + /** + Return function value at current upper bound of the minimization interval + */ + double FValUpper() const; + + + /** + Find minimum position iterating until convergence specified by the absolute and relative tolerance or + the maximum number of iteration is reached + Return true is result is successfull + \@param maxIter maximum number of iteration + \@param absTol desired absolute error in the minimum position + \@param absTol desired relative error in the minimum position + */ + bool Minimize( int maxIter, double absTol, double relTol); + + + /** + Return number of iteration used to find minimum + */ + int Iterations() const { + return fIter; + } + + /** + Return status of last minimization + */ + int Status() const { return fStatus; } + + /** + Return name of minimization algorithm + */ + const char * Name() const; + + /** + Test convergence of the interval. + The test returns success if + \f[ + |x_{min}-x_{truemin}| < epsAbs + epsRel *x_{truemin} + \f] + */ + static int TestInterval( double xlow, double xup, double epsAbs, double epsRel); + + + protected: + + + private: + + double fXmin; + double fXlow; + double fXup; + double fMin; + double fLow; + double fUp; + int fIter; + int fStatus; // status of last minimization (==0 ok =1 failed) + bool fIsSet; + + + GSL1DMinimizerWrapper * fMinimizer; + GSLFunctionWrapper * fFunction; + + }; + +} // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_GSLMinimizer1D */ diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLNLSMinimizer.h b/ThirdParty/RootMinimizers/inc/Math/GSLNLSMinimizer.h new file mode 100644 index 00000000000..5094e03b619 --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLNLSMinimizer.h @@ -0,0 +1,297 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Dec 20 17:16:32 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLNLSMinimizer + +#ifndef ROOT_Math_GSLNLSMinimizer +#define ROOT_Math_GSLNLSMinimizer + + + +#ifndef ROOT_Math_Minimizer +#include "Math/Minimizer.h" +#endif + + +#ifndef ROOT_Math_IFunctionfwd +#include "Math/IFunctionfwd.h" +#endif + +#ifndef ROOT_Math_IParamFunctionfwd +#include "Math/IParamFunctionfwd.h" +#endif + +#ifndef ROOT_Math_FitMethodFunction +#include "Math/FitMethodFunction.h" +#endif + +#ifndef ROOT_Math_MinimizerVariable +#include "Math/MinimizerVariable.h" +#endif + + +#include <vector> +#include <map> +#include <string> + +namespace ROOT { + + namespace Math { + + class GSLMultiFit; + + +//________________________________________________________________________________ +/** + LSResidualFunc class description. + Internal class used for accessing the residuals of the Least Square function + and their derivates which are estimated numerically using GSL numerical derivation. + The class contains a pointer to the fit method function and an index specifying + the i-th residual and wraps it in a multi-dim gradient function interface + ROOT::Math::IGradientFunctionMultiDim. + The class is used by ROOT::Math::GSLNLSMinimizer (GSL non linear least square fitter) + + @ingroup MultiMin +*/ +class LSResidualFunc : public IMultiGradFunction { +public: + + //default ctor (required by CINT) + LSResidualFunc() : fIndex(0), fChi2(0) + {} + + + LSResidualFunc(const ROOT::Math::FitMethodFunction & func, unsigned int i) : + fIndex(i), + fChi2(&func), + fX2(std::vector<double>(func.NDim() ) ) + {} + + + // copy ctor + LSResidualFunc(const LSResidualFunc & rhs) : + IMultiGenFunction(), + IMultiGradFunction() + { + operator=(rhs); + } + + // assignment + LSResidualFunc & operator= (const LSResidualFunc & rhs) + { + fIndex = rhs.fIndex; + fChi2 = rhs.fChi2; + fX2 = rhs.fX2; + return *this; + } + + IMultiGenFunction * Clone() const { + return new LSResidualFunc(*fChi2,fIndex); + } + + unsigned int NDim() const { return fChi2->NDim(); } + + void Gradient( const double * x, double * g) const { + double f0 = 0; + FdF(x,f0,g); + } + + void FdF (const double * x, double & f, double * g) const { + unsigned int n = NDim(); + std::copy(x,x+n,fX2.begin()); + const double kEps = 1.0E-4; + f = DoEval(x); + for (unsigned int i = 0; i < n; ++i) { + fX2[i] += kEps; + g[i] = ( DoEval(&fX2.front()) - f )/kEps; + fX2[i] = x[i]; + } + } + + +private: + + double DoEval (const double * x) const { + return fChi2->DataElement(x, fIndex); + } + + double DoDerivative(const double * x, unsigned int icoord) const { + //return ROOT::Math::Derivator::Eval(*this, x, icoord, 1E-8); + std::copy(x,x+NDim(),fX2.begin()); + const double kEps = 1.0E-4; + fX2[icoord] += kEps; + return ( DoEval(&fX2.front()) - DoEval(x) )/kEps; + } + + unsigned int fIndex; + const ROOT::Math::FitMethodFunction * fChi2; + mutable std::vector<double> fX2; // cached vector +}; + + +//_____________________________________________________________________________________________________ +/** + GSLNLSMinimizer class for Non Linear Least Square fitting + It Uses the Levemberg-Marquardt algorithm from + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Nonlinear-Least_002dSquares-Fitting.html"> + GSL Non Linear Least Square fitting</A>. + + @ingroup MultiMin +*/ +class GSLNLSMinimizer : public ROOT::Math::Minimizer { + +public: + + /** + Default constructor + */ + GSLNLSMinimizer (int type = 0); + + /** + Destructor (no operations) + */ + ~GSLNLSMinimizer (); + +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + GSLNLSMinimizer(const GSLNLSMinimizer &) : ROOT::Math::Minimizer() {} + + /** + Assignment operator + */ + GSLNLSMinimizer & operator = (const GSLNLSMinimizer & rhs) { + if (this == &rhs) return *this; // time saving self-test + return *this; + } + +public: + + /// set the function to minimize + virtual void SetFunction(const ROOT::Math::IMultiGenFunction & func); + + /// set gradient the function to minimize + virtual void SetFunction(const ROOT::Math::IMultiGradFunction & func); + + /// set free variable + virtual bool SetVariable(unsigned int ivar, const std::string & name, double val, double step); + + /// set lower limited variable + virtual bool SetLowerLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower ); + /// set upper limited variable + virtual bool SetUpperLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double upper ); + /// set upper/lower limited variable + virtual bool SetLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower , double upper ); + /// set fixed variable + virtual bool SetFixedVariable(unsigned int ivar , const std::string & name , double val ); + + /// set the value of an existing variable + virtual bool SetVariableValue(unsigned int ivar, double val ); + /// set the values of all existing variables (array must be dimensioned to the size of existing parameters) + virtual bool SetVariableValues(const double * x); + + /// method to perform the minimization + virtual bool Minimize(); + + /// return minimum function value + virtual double MinValue() const { return fMinVal; } + + /// return expected distance reached from the minimum + virtual double Edm() const { return fEdm; } // not impl. } + + /// return pointer to X values at the minimum + virtual const double * X() const { return &fValues.front(); } + + /// return pointer to gradient values at the minimum + virtual const double * MinGradient() const; + + /// number of function calls to reach the minimum + virtual unsigned int NCalls() const { return (fObjFunc) ? fObjFunc->NCalls() : 0; } + + /// this is <= Function().NDim() which is the total + /// number of variables (free+ constrained ones) + virtual unsigned int NDim() const { return fDim; } + + /// number of free variables (real dimension of the problem) + /// this is <= Function().NDim() which is the total + virtual unsigned int NFree() const { return fNFree; } + + /// minimizer provides error and error matrix + virtual bool ProvidesError() const { return true; } + + /// return errors at the minimum + virtual const double * Errors() const { return (fErrors.size() > 0) ? &fErrors.front() : 0; } +// { +// static std::vector<double> err; +// err.resize(fDim); +// return &err.front(); +// } + + /** return covariance matrices elements + if the variable is fixed the matrix is zero + The ordering of the variables is the same as in errors + */ + virtual double CovMatrix(unsigned int , unsigned int ) const; + + /// return covariance matrix status + virtual int CovMatrixStatus() const; + +protected: + +private: + + + unsigned int fDim; // dimension of the function to be minimized + unsigned int fNFree; // dimension of the internal function to be minimized + unsigned int fSize; // number of fit points (residuals) + + + ROOT::Math::GSLMultiFit * fGSLMultiFit; // pointer to GSL multi fit solver + const ROOT::Math::FitMethodFunction * fObjFunc; // pointer to Least square function + + double fMinVal; // minimum function value + double fEdm; // edm value + double fLSTolerance; // Line Search Tolerance + std::vector<double> fValues; + std::vector<double> fErrors; + std::vector<double> fCovMatrix; // cov matrix (stored as cov[ i * dim + j] + std::vector<double> fSteps; + std::vector<std::string> fNames; + std::vector<LSResidualFunc> fResiduals; //! transient Vector of the residual functions + + std::vector<ROOT::Math::EMinimVariableType> fVarTypes; // vector specifyng the type of variables + std::map< unsigned int, std::pair<double, double> > fBounds; // map specifying the bound using as key the parameter index + + +}; + + } // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_GSLNLSMinimizer */ diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h b/ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h new file mode 100644 index 00000000000..67fd5f21785 --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h @@ -0,0 +1,437 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLRandom +// +// Created by: moneta at Sun Nov 21 16:26:03 2004 +// +// Last update: Sun Nov 21 16:26:03 2004 +// +#ifndef ROOT_Math_GSLRndmEngines +#define ROOT_Math_GSLRndmEngines + +#include <string> +#include <vector> + +namespace ROOT { +namespace Math { + + + class GSLRngWrapper; + + //_________________________________________________________________ + /** + GSLRandomEngine + Base class for all GSL random engines, + normally user instantiate the derived classes + which creates internally the generator. + + The main GSL generators (see + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html"> + here</A>) are available as derived classes + In addition to generate uniform numbers it provides method for + generating numbers according to pre-defined distributions + using the GSL functions from + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distributions.html"> + GSL random number distributions</A>. + + + + @ingroup Random + */ + class GSLRandomEngine { + + public: + + /** + default constructor. No creation of rng is done. + If then Initialize() is called an engine is created + based on default GSL type (MT) + */ + GSLRandomEngine(); + + /** + create from an existing rng. + User manage the rng pointer which is then deleted olny by calling Terminate() + */ + GSLRandomEngine( GSLRngWrapper * rng); + + /** + initialize the generator + If no rng is present the default one based on Mersenne and Twister is created + */ + void Initialize(); + + /** + delete pointer to contained rng + */ + void Terminate(); + + /** + call Terminate() + */ + virtual ~GSLRandomEngine(); + + /** + Generate a random number between ]0,1] + 0 is excluded and 1 is included + */ + double operator() () const; + + /** + Generate an integer number between [0,max-1] (including 0 and max-1) + if max is larger than available range of algorithm + an error message is printed and zero is returned + */ + unsigned int RndmInt(unsigned int max) const; + + /** + Generate an array of random numbers. + The iterators points to the random numbers + */ + template<class Iterator> + void RandomArray(Iterator begin, Iterator end) const { + for ( Iterator itr = begin; itr != end; ++itr ) { + *itr = this->operator()(); + } + } + + /** + Generate an array of random numbers + The iterators points to the random numbers + */ + void RandomArray(double * begin, double * end) const; + + /** + return name of generator + */ + std::string Name() const; + + /** + return the state size of generator + */ + unsigned int Size() const; + + /** + set the random generator seed + */ + void SetSeed(unsigned int seed) const; + + + /** @name Random Distributions + Implemented using the + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distributions.html"> + GSL Random number Distributions</A> + **/ + //@{ + /** + Gaussian distribution - default method is Box-Muller (polar method) + */ + double Gaussian(double sigma) const; + + /** + Gaussian distribution - Ziggurat method + */ + double GaussianZig(double sigma) const; + + /** + Gaussian distribution - Ratio method + */ + double GaussianRatio(double sigma) const; + /** + Gaussian Tail distribution + */ + double GaussianTail(double a, double sigma) const; + + /** + Bivariate Gaussian distribution with correlation + */ + void Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const; + + /** + Exponential distribution + */ + double Exponential(double mu) const; + + /** + Cauchy distribution + */ + double Cauchy(double a) const; + + /** + Landau distribution + */ + double Landau() const; + + /** + Gamma distribution + */ + double Gamma(double a, double b) const; + + /** + Log Normal distribution + */ + double LogNormal(double zeta, double sigma) const; + + /** + Chi square distribution + */ + double ChiSquare(double nu) const; + + /** + F distrbution + */ + double FDist(double nu1, double nu2) const; + + /** + t student distribution + */ + double tDist(double nu) const; + + /** + generate random numbers in a 2D circle of radious 1 + */ + void Dir2D(double &x, double &y) const; + + /** + generate random numbers in a 3D sphere of radious 1 + */ + void Dir3D(double &x, double &y, double &z) const; + + /** + Poisson distribution + */ + unsigned int Poisson(double mu) const; + + /** + Binomial distribution + */ + unsigned int Binomial(double p, unsigned int n) const; + + /** + Negative Binomial distribution + */ + unsigned int NegativeBinomial(double p, double n) const; + + /** + Multinomial distribution + */ + std::vector<unsigned int> Multinomial( unsigned int ntot, const std::vector<double> & p ) const; + + //@} + + + + protected: + + /// internal method used by the derived class to set the type of generators + void SetType(GSLRngWrapper * r) { + fRng = r; + } + + private: + + GSLRngWrapper * fRng; // pointer to GSL generator wrapper (managed by the class) + mutable unsigned int fCurTime; // current time used to seed the generator + + + }; + + //_____________________________________________________________________________________ + /** + Mersenne-Twister generator + gsl_rng_mt19937 from + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + + @ingroup Random + */ + class GSLRngMT : public GSLRandomEngine { + public: + GSLRngMT(); + }; + + //_____________________________________________________________________________________ + /** + Old Ranlux generator (James, Luscher) (default luxury level, p = 223) + (This is eequivalent to TRandom1 with default luxury level) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngRanLux : public GSLRandomEngine { + public: + GSLRngRanLux(); + }; + + //_____________________________________________________________________________________ + /** + Second generation of Ranlux generator for single precision with luxury level of 1 + (It throws away 202 values for every 12 used) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngRanLuxS1 : public GSLRandomEngine { + public: + GSLRngRanLuxS1(); + }; + typedef GSLRngRanLuxS1 GSLRngRanLux1; // for backward compatibility + + //_____________________________________________________________________________________ + /** + Second generation of Ranlux generator for Single precision with luxury level of 2 + (It throws away 397 value for every 12 used) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngRanLuxS2 : public GSLRandomEngine { + public: + GSLRngRanLuxS2(); + }; + typedef GSLRngRanLuxS2 GSLRngRanLux2; // for backward compatibility + + //_____________________________________________________________________________________ + /** + Double precision (48 bits) version of Second generation of Ranlux generator with luxury level of 1 + (It throws away 202 value for every 12 used) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngRanLuxD1 : public GSLRandomEngine { + public: + GSLRngRanLuxD1(); + }; + + //_____________________________________________________________________________________ + /** + Double precision (48 bits) version of Second generation of Ranlux generator with luxury level of 2 + (It throws away 397 value for every 12 used) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngRanLuxD2 : public GSLRandomEngine { + public: + GSLRngRanLuxD2(); + }; + typedef GSLRngRanLuxD2 GSLRngRanLux48; // for backward compatibility + + + //_____________________________________________________________________________________ + /** + Tausworthe generator by L'Ecuyer + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngTaus : public GSLRandomEngine { + public: + GSLRngTaus(); + }; + + //_____________________________________________________________________________________ + /** + Lagged Fibonacci generator by Ziff + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngGFSR4 : public GSLRandomEngine { + public: + GSLRngGFSR4(); + }; + + //_____________________________________________________________________________________ + /** + Combined multiple recursive generator (L'Ecuyer) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngCMRG : public GSLRandomEngine { + public: + GSLRngCMRG(); + }; + + //_____________________________________________________________________________________ + /** + 5-th order multiple recursive generator (L'Ecuyer, Blouin and Coutre) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html">here</A> + + @ingroup Random + */ + class GSLRngMRG : public GSLRandomEngine { + public: + GSLRngMRG(); + }; + + //_____________________________________________________________________________________ + /** + BSD rand() generator + gsl_rmg_rand from + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Unix-random-number-generators.html">here</A> + + @ingroup Random + */ + class GSLRngRand : public GSLRandomEngine { + public: + GSLRngRand(); + }; + + //_____________________________________________________________________________________ + /** + RANMAR generator + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Unix-random-number-generators.html">here</A> + + @ingroup Random + */ + class GSLRngRanMar : public GSLRandomEngine { + public: + GSLRngRanMar(); + }; + + //_____________________________________________________________________________________ + /** + MINSTD generator (Park and Miller) + see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Unix-random-number-generators.html">here</A> + + @ingroup Random + */ + class GSLRngMinStd : public GSLRandomEngine { + public: + GSLRngMinStd(); + }; + + + + +} // namespace Math +} // namespace ROOT + + +#endif /* ROOT_Math_GSLRndmEngines */ + diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLSimAnMinimizer.h b/ThirdParty/RootMinimizers/inc/Math/GSLSimAnMinimizer.h new file mode 100644 index 00000000000..80e7dace8bd --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLSimAnMinimizer.h @@ -0,0 +1,204 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Dec 20 17:16:32 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLSimAnMinimizer + +#ifndef ROOT_Math_GSLSimAnMinimizer +#define ROOT_Math_GSLSimAnMinimizer + + + +#ifndef ROOT_Math_Minimizer +#include "Math/Minimizer.h" +#endif + + +#ifndef ROOT_Math_IFunctionfwd +#include "Math/IFunctionfwd.h" +#endif + +#ifndef ROOT_Math_IParamFunctionfwd +#include "Math/IParamFunctionfwd.h" +#endif + + + +#ifndef ROOT_Math_GSLSimAnnealing +#include "Math/GSLSimAnnealing.h" +#endif + +#ifndef ROOT_Math_MinimizerVariable +#include "Math/MinimizerVariable.h" +#endif + +#include <vector> +#include <map> + + + +namespace ROOT { + + namespace Math { + + class MinimTransformFunction; + + +//_____________________________________________________________________________________ +/** + GSLSimAnMinimizer class for minimization using simulated annealing + using the algorithm from + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Simulated-Annealing.html"> + GSL</A>. + It implements the ROOT::Minimizer interface and + a plug-in (name "GSLSimAn") exists to instantiate this class via the plug-in manager + + @ingroup MultiMin +*/ +class GSLSimAnMinimizer : public ROOT::Math::Minimizer { + +public: + + /** + Default constructor + */ + GSLSimAnMinimizer (int type = 0); + + /** + Destructor (no operations) + */ + ~GSLSimAnMinimizer (); + +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + GSLSimAnMinimizer(const GSLSimAnMinimizer &) : ROOT::Math::Minimizer() {} + + /** + Assignment operator + */ + GSLSimAnMinimizer & operator = (const GSLSimAnMinimizer & rhs) { + if (this == &rhs) return *this; // time saving self-test + return *this; + } + +public: + + /// set the function to minimize + virtual void SetFunction(const ROOT::Math::IMultiGenFunction & func); + + /// set gradient the function to minimize + virtual void SetFunction(const ROOT::Math::IMultiGradFunction & func); + + /// set free variable + virtual bool SetVariable(unsigned int ivar, const std::string & name, double val, double step); + + /// set fixed variable (override if minimizer supports them ) + virtual bool SetFixedVariable(unsigned int /* ivar */, const std::string & /* name */, double /* val */); + + /// set lower limit variable (override if minimizer supports them ) + virtual bool SetLowerLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower ); + /// set upper limit variable (override if minimizer supports them ) + virtual bool SetUpperLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double upper ); + /// set upper/lower limited variable (override if minimizer supports them ) + virtual bool SetLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double /* lower */, double /* upper */); + + /// set the value of an existing variable + virtual bool SetVariableValue(unsigned int ivar, double val ); + /// set the values of all existing variables (array must be dimensioned to the size of the existing parameters) + virtual bool SetVariableValues(const double * x); + + + /// method to perform the minimization + virtual bool Minimize(); + + /// return minimum function value + virtual double MinValue() const { return fMinVal; } + + /// return expected distance reached from the minimum + virtual double Edm() const { return 0; } // not impl. } + + /// return pointer to X values at the minimum + virtual const double * X() const { return &fValues.front(); } + + /// return pointer to gradient values at the minimum + virtual const double * MinGradient() const { return 0; } // not impl. + + /// number of function calls to reach the minimum + virtual unsigned int NCalls() const { return 0; } // not yet ipl. + + /// this is <= Function().NDim() which is the total + /// number of variables (free+ constrained ones) + virtual unsigned int NDim() const { return fDim; } + + /// number of free variables (real dimension of the problem) + /// this is <= Function().NDim() which is the total + virtual unsigned int NFree() const { return fDim; } + + /// minimizer provides error and error matrix + virtual bool ProvidesError() const { return false; } + + /// return errors at the minimum + virtual const double * Errors() const { return 0; } + + /** return covariance matrices elements + if the variable is fixed the matrix is zero + The ordering of the variables is the same as in errors + */ + virtual double CovMatrix(unsigned int , unsigned int ) const { return 0; } + + + /// return reference to the objective function + ///virtual const ROOT::Math::IGenFunction & Function() const; + + +protected: + +private: + + unsigned int fDim; // dimension of the function to be minimized + bool fOwnFunc; // flag to indicate if objective function is managed + + ROOT::Math::GSLSimAnnealing fSolver; + const ROOT::Math::IMultiGenFunction * fObjFunc; + + double fMinVal; // minimum values + + mutable std::vector<double> fValues; + + std::vector<double> fSteps; + std::vector<std::string> fNames; + std::vector<ROOT::Math::EMinimVariableType> fVarTypes; // vector specifyng the type of variables + std::map< unsigned int, std::pair<double, double> > fBounds; // map specifying the bound using as key the parameter index + +}; + + } // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_GSLSimAnMinimizer */ diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLSimAnnealing.h b/ThirdParty/RootMinimizers/inc/Math/GSLSimAnnealing.h new file mode 100644 index 00000000000..eb5705ee02d --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/GSLSimAnnealing.h @@ -0,0 +1,257 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Thu Jan 25 11:13:48 2007 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLSimAnnealing + +#ifndef ROOT_Math_GSLSimAnnealing +#define ROOT_Math_GSLSimAnnealing + +#include "Math/IFunctionfwd.h" + +#include <vector> + +namespace ROOT { + + namespace Math { + + class GSLRandomEngine; + +//_____________________________________________________________________________ +/** + GSLSimAnFunc class description. + Interface class for the objetive function to be used in simulated annealing + If user wants to re-implement some of the methods (like the one defining the metric) which are used by the + the simulated annealing algorithm must build a user derived class. + NOTE: Derived classes must re-implement the assignment and copy constructor to call them of the parent class + + @ingroup MultiMin + */ +class GSLSimAnFunc { +public: + + /** + construct from an interface of a multi-dimensional function + */ + GSLSimAnFunc(const ROOT::Math::IMultiGenFunction & func, const double * x); + + /** + construct from an interface of a multi-dimensional function + Use optionally a scale factor (for each coordinate) which can be used to scale the step sizes + (this is used for example by the minimization algorithm) + */ + GSLSimAnFunc(const ROOT::Math::IMultiGenFunction & func, const double * x, const double * scale); + +protected: + + /** + derived classes might need to re-define completely the class + */ + GSLSimAnFunc() : + fFunc(0) + {} + +public: + + + /// virtual distructor (no operations) + virtual ~GSLSimAnFunc() { } // + + + /** + fast copy method called by GSL simuated annealing internally + copy only the things which have been changed + must be re-implemented by derived classes if needed + */ + virtual GSLSimAnFunc & FastCopy(const GSLSimAnFunc & f); + + + /** + clone method. Needs to be re-implemented by the derived classes for deep copying + */ + virtual GSLSimAnFunc * Clone() const { + return new GSLSimAnFunc(*this); + } + + /** + evaluate the energy ( objective function value) + re-implement by derived classes if needed to be modified + */ + virtual double Energy() const; + + /** + change the x[i] value using a random value urndm generated between [0,1] + up to a maximum value maxstep + re-implement by derived classes if needed to be modified + */ + virtual void Step(const GSLRandomEngine & r, double maxstep); + + /** + calculate the distance (metric) between this one and another configuration + Presently a cartesian metric is used. + re-implement by derived classes if needed to be modified + */ + virtual double Distance(const GSLSimAnFunc & func) const; + + /** + print the position in the standard output ostream + GSL prints in addition n iteration, n function calls, temperature and energy + re-implement by derived classes if necessary + */ + virtual void Print(); + + /** + change the x values (used by sim annealing to take a step) + */ + void SetX(const double * x) { + std::copy(x, x+ fX.size(), fX.begin() ); + } + + template <class IT> + void SetX(IT begin, IT end) { + std::copy(begin, end, fX.begin() ); + } + + unsigned int NDim() const { return fX.size(); } + + double X(unsigned int i) const { return fX[i]; } + + const std::vector<double> & X() const { return fX; } + + double Scale(unsigned int i) const { return fScale[i]; } + + void SetX(unsigned int i, double x) { fX[i] = x; } + + // use compiler generated copy ctror and assignment operators + +private: + + std::vector<double> fX; + std::vector<double> fScale; + const ROOT::Math::IMultiGenFunction * fFunc; + +}; + +//_____________________________________________________ +/** + structure holding the simulated annealing parameters + + @ingroup MultiMin +*/ +struct GSLSimAnParams { + + // constructor with some default values + GSLSimAnParams() { + n_tries = 200; + iters_fixed_T = 10; + step_size = 10; + // the following parameters are for the Boltzmann distribution */ + k = 1.0; + t_initial = 0.002; + mu = 1.005; + t_min = 2.0E-6; + } + + + int n_tries; // number of points to try for each step + int iters_fixed_T; // number of iterations at each temperature + double step_size; // max step size used in random walk + /// parameters for the Boltzman distribution + double k; + double t_initial; + double mu; + double t_min; +}; + +//___________________________________________________________________________ +/** + GSLSimAnnealing class for performing a simulated annealing search of + a multidimensional function + + @ingroup MultiMin +*/ +class GSLSimAnnealing { + +public: + + /** + Default constructor + */ + GSLSimAnnealing (); + + /** + Destructor (no operations) + */ + ~GSLSimAnnealing () {} + +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + GSLSimAnnealing(const GSLSimAnnealing &) {} + + /** + Assignment operator + */ + GSLSimAnnealing & operator = (const GSLSimAnnealing & rhs) { + if (this == &rhs) return *this; // time saving self-test + return *this; + } + +public: + + + /** + solve the simulated annealing given a multi-dim function, the initial vector parameters + and a vector containing the scaling factors for the parameters + */ + int Solve(const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * scale, double * xmin, bool debug = false); + + /** + solve the simulated annealing given a GSLSimAnFunc object + The object will contain the initial state at the beginning and the final minimum state at the end + */ + int Solve(GSLSimAnFunc & func, bool debug = false); + + + GSLSimAnParams & Params() { return fParams; } + const GSLSimAnParams & Params() const { return fParams; } + + +protected: + + +private: + + GSLSimAnParams fParams; // parameters for GSLSimAnnealig + +}; + + } // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_GSLSimAnnealing */ diff --git a/ThirdParty/RootMinimizers/inc/Math/MultiNumGradFunction.h b/ThirdParty/RootMinimizers/inc/Math/MultiNumGradFunction.h new file mode 100644 index 00000000000..ddebcc28477 --- /dev/null +++ b/ThirdParty/RootMinimizers/inc/Math/MultiNumGradFunction.h @@ -0,0 +1,143 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Dec 20 14:36:31 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class MultiNumGradFunction + +#ifndef ROOT_Math_MultiNumGradFunction +#define ROOT_Math_MultiNumGradFunction + + +#ifndef ROOT_Math_IFunction +#include "Math/IFunction.h" +#endif + +#ifndef ROOT_Math_WrappedFunction +#include "Math/WrappedFunction.h" +#endif + + +namespace ROOT { + + namespace Math { + + +/** + MultiNumGradFunction class to wrap a normal function in a + gradient function using numerical gradient calculation + provided by the class Derivator (based on GSL numerical derivation) + + + @ingroup MultiMin +*/ +class MultiNumGradFunction : public IMultiGradFunction { + +public: + + + /** + Constructor from a IMultiGenFunction interface + */ + MultiNumGradFunction (const IMultiGenFunction & f) : + fFunc(&f), + fDim(f.NDim() ), + fNCalls(0), + fOwner(false) + {} + + /** + Constructor from a generic function (pointer or reference) and number of dimension + implementiong operator () (double * x) + */ + + template<class FuncType> + MultiNumGradFunction (FuncType f, int n) : + fDim( n ), + fNCalls(0), + fOwner(true) + { + // create a wrapped function + fFunc = new ROOT::Math::WrappedMultiFunction<FuncType> (f, n); + } + + /** + Destructor (no operations) + */ + ~MultiNumGradFunction () { + if (fOwner) delete fFunc; + } + + + // method inheritaed from IFunction interface + + unsigned int NDim() const { return fDim; } + + unsigned int NCalls() const { return fNCalls; } + + IMultiGenFunction * Clone() const { + if (!fOwner) + return new MultiNumGradFunction(*fFunc); + else { + // we need to copy the pointer to the wrapped function + MultiNumGradFunction * f = new MultiNumGradFunction(*(fFunc->Clone()) ); + f->fOwner = true; + return f; + } + } + + + /// precision value used for calculating the derivative step-size + /// h = eps * |x|. The default is 0.001, give a smaller in case function chanes rapidly + static void SetDerivPrecision(double eps); + + /// get precision value used for calculating the derivative step-size + static double GetDerivPrecision(); + + +private: + + + double DoEval(const double * x) const { + fNCalls++; + return (*fFunc)(x); + } + + // calculate derivative using mathcore derivator + double DoDerivative (const double * x, unsigned int icoord ) const; + + // adapat internal function type to IMultiGenFunction needed by derivative calculation + const IMultiGenFunction * fFunc; + unsigned int fDim; + mutable unsigned int fNCalls; + bool fOwner; + + static double fgEps; // epsilon used in derivative calculation h ~ eps |x| + +}; + + } // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_NumGradFunction */ diff --git a/ThirdParty/RootMinimizers/src/Math/Derivator.cxx b/ThirdParty/RootMinimizers/src/Math/Derivator.cxx new file mode 100644 index 00000000000..17158cdce3d --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/Derivator.cxx @@ -0,0 +1,161 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Implementation file for class GSLDerivator +// +// Created by: moneta at Sat Nov 13 14:46:00 2004 +// +// Last update: Sat Nov 13 14:46:00 2004 +// + +#include "Math/IFunction.h" +#include "Math/IParamFunction.h" +#include "Math/Derivator.h" +#include "GSLDerivator.h" + +#include "Math/OneDimFunctionAdapter.h" + +// for GSL greater then 1.5 +#include "gsl/gsl_deriv.h" +// for OLD GSL versions +//#include "gsl/gsl_diff.h" + +namespace ROOT { +namespace Math { + +Derivator::Derivator() { + fDerivator = new GSLDerivator(); +} + +Derivator::Derivator(const IGenFunction &f) +{ + // allocate a GSLDerivator + fDerivator = new GSLDerivator(); + fDerivator->SetFunction(f); +} + +Derivator::Derivator(const GSLFuncPointer &f, void * p) +{ + // allocate a GSLDerivator + fDerivator = new GSLDerivator(); + fDerivator->SetFunction(f,p); + +} + +Derivator::~Derivator() +{ + if (fDerivator) delete fDerivator; +} + + +Derivator::Derivator(const Derivator &) +{ +} + +Derivator & Derivator::operator = (const Derivator &rhs) +{ + if (this == &rhs) return *this; // time saving self-test + + return *this; +} + + +void Derivator::SetFunction(const IGenFunction &f) { + fDerivator->SetFunction(f); +} + +void Derivator::SetFunction( const GSLFuncPointer &f, void * p) { + fDerivator->SetFunction(f,p); +} + + +double Derivator::Eval( double x, double h) const { + return fDerivator->EvalCentral(x, h); +} + +double Derivator::EvalCentral( double x, double h) const { + return fDerivator->EvalCentral(x, h); +} + +double Derivator::EvalForward( double x, double h) const { + return fDerivator->EvalForward(x, h); +} + +double Derivator::EvalBackward( double x, double h) const { + return fDerivator->EvalBackward(x, h); +} + +// static methods +double Derivator::Eval(const IGenFunction & f, double x, double h ) { + return GSLDerivator::EvalCentral(f, x, h ); +} + +double Derivator::EvalCentral(const IGenFunction & f, double x, double h) { + return GSLDerivator::EvalCentral(f,x,h); +} + +double Derivator::EvalForward(const IGenFunction & f, double x, double h) { + return GSLDerivator::EvalForward(f, x, h); +} + +double Derivator::EvalBackward(const IGenFunction & f, double x, double h) { + return GSLDerivator::EvalBackward(f, x, h); +} + +double Derivator::Eval(const IMultiGenFunction & f, const double * x, unsigned int icoord, double h ) { + // partial derivative for a multi-dim function + GSLDerivator d; + OneDimMultiFunctionAdapter<> adapter(f,x,icoord); + d.SetFunction( &GSLFunctionAdapter<OneDimMultiFunctionAdapter<> >::F,static_cast<void *>(&adapter) ); + return d.EvalCentral(x[icoord],h); +} + +double Derivator::Eval(IParamFunction & f, double x, const double * p, unsigned int ipar, double h ) { + // derivative w.r.t parameter for a one-dim param function + GSLDerivator d; + const double xx = x; + OneDimParamFunctionAdapter<IParamFunction &> adapter(f,&xx,p,ipar); + d.SetFunction( &GSLFunctionAdapter<OneDimParamFunctionAdapter<IParamFunction &> >::F,static_cast<void *>(&adapter) ); + return d.EvalCentral(p[ipar],h); +} + +double Derivator::Eval(IParamMultiFunction & f, const double * x, const double * p, unsigned int ipar, double h ) { + // derivative w.r.t parameter for a multi-dim param function + GSLDerivator d; + OneDimParamFunctionAdapter<IParamMultiFunction &> adapter(f,x,p,ipar); + d.SetFunction( &GSLFunctionAdapter<OneDimParamFunctionAdapter<IParamMultiFunction &> >::F,static_cast<void *>(&adapter) ); + return d.EvalCentral(p[ipar],h); +} + + +double Derivator::Result() const { return fDerivator->Result(); } + +double Derivator::Error() const { return fDerivator->Error(); } + +int Derivator::Status() const { return fDerivator->Status(); } + + + +} // namespace Math +} // namespace ROOT diff --git a/ThirdParty/RootMinimizers/src/Math/GSL1DMinimizerWrapper.h b/ThirdParty/RootMinimizers/src/Math/GSL1DMinimizerWrapper.h new file mode 100644 index 00000000000..114a328d43a --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSL1DMinimizerWrapper.h @@ -0,0 +1,76 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + /********************************************************************** + * * + * Copyright (c) 2004 moneta, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSL1DMinimizerWrapper +// +// Created by: moneta at Wed Dec 1 17:25:44 2004 +// +// Last update: Wed Dec 1 17:25:44 2004 +// +#ifndef ROOT_Math_GSL1DMinimizerWrapper +#define ROOT_Math_GSL1DMinimizerWrapper + +#include "gsl/gsl_min.h" + + +namespace ROOT { + +namespace Math { + +/** + wrapper class for gsl_min_fminimizer structure + @ingroup Min1D +*/ +class GSL1DMinimizerWrapper { + +public: + GSL1DMinimizerWrapper( const gsl_min_fminimizer_type * T) + { + fMinimizer = gsl_min_fminimizer_alloc(T); + } + virtual ~GSL1DMinimizerWrapper() { + gsl_min_fminimizer_free(fMinimizer); + } + +private: +// usually copying is non trivial, so we make this unaccessible + GSL1DMinimizerWrapper(const GSL1DMinimizerWrapper &); + GSL1DMinimizerWrapper & operator = (const GSL1DMinimizerWrapper &); + +public: + + gsl_min_fminimizer * Get() const { + return fMinimizer; + } + + +private: + + gsl_min_fminimizer * fMinimizer; + +}; + +} // end namespace Math +} // end namespace ROOT + +#endif /* ROOT_Math_GSL1DMinimizerWrapper */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLDerivator.cxx b/ThirdParty/RootMinimizers/src/Math/GSLDerivator.cxx new file mode 100644 index 00000000000..ff9ae35ded1 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLDerivator.cxx @@ -0,0 +1,128 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Implementation file for class GSLDerivator +// +// Created by: moneta at Sat Nov 13 14:46:00 2004 +// +// Last update: Sat Nov 13 14:46:00 2004 +// + +#include "GSLDerivator.h" + +#include "GSLFunctionWrapper.h" +// for GSL greater then 1.5 +#include "gsl/gsl_deriv.h" +// for OLD GSL versions +//#include "gsl/gsl_diff.h" + +#include <iostream> + +namespace ROOT { +namespace Math { + + + +double GSLDerivator::EvalCentral( double x, double h) { + // Central evaluation using previously set function + if ( !fFunction.IsValid() ) { + std::cerr << "GSLDerivator: Error : The function has not been specified" << std::endl; + fStatus = -1; + return 0; + } + fStatus = gsl_deriv_central( fFunction.GetFunc(), x, h, &fResult, &fError); + return fResult; +} + +double GSLDerivator::EvalForward( double x, double h) { + // Forward evaluation using previously set function + if ( !fFunction.IsValid() ) { + std::cerr << "GSLDerivator: Error : The function has not been specified" << std::endl; + fStatus = -1; + return 0; + } + fStatus = gsl_deriv_forward( fFunction.GetFunc(), x, h, &fResult, &fError); + return fResult; +} + +double GSLDerivator::EvalBackward( double x, double h) { + // Backward evaluation using previously set function + if ( !fFunction.IsValid() ) { + std::cerr << "GSLDerivator: Error : The function has not been specified" << std::endl; + fStatus = -1; + return 0; + } + fStatus = gsl_deriv_backward( fFunction.GetFunc(), x, h, &fResult, &fError); + return fResult; +} + +// static methods not requiring the function +double GSLDerivator::EvalCentral(const IGenFunction & f, double x, double h) { + // Central evaluation using given function + GSLFunctionWrapper gslfw; + double result, error = 0; + gslfw.SetFunction(f); + gsl_deriv_central( gslfw.GetFunc(), x, h, &result, &error); + return result; +} + +double GSLDerivator::EvalForward(const IGenFunction & f, double x, double h) { + // Forward evaluation using given function + GSLFunctionWrapper gslfw; + double result, error = 0; + gslfw.SetFunction(f); + gsl_deriv_forward( gslfw.GetFunc(), x, h, &result, &error); + return result; +} + +double GSLDerivator::EvalBackward(const IGenFunction & f, double x, double h) { + // Backward evaluation using given function + GSLFunctionWrapper gslfw; + double result, error = 0; + gslfw.SetFunction(f); + gsl_deriv_backward( gslfw.GetFunc(), x, h, &result, &error); + return result; +} + + +double GSLDerivator::Result() const { return fResult; } + +double GSLDerivator::Error() const { return fError; } + +int GSLDerivator::Status() const { return fStatus; } + +// fill GSLFunctionWrapper with the pointer to the function + +void GSLDerivator::SetFunction( GSLFuncPointer fp, void * p) { + fFunction.SetFuncPointer( fp ); + fFunction.SetParams ( p ); +} + + +void GSLDerivator::SetFunction(const IGenFunction &f) { + fFunction.SetFunction(f); +} + +} // namespace Math +} // namespace ROOT diff --git a/ThirdParty/RootMinimizers/src/Math/GSLDerivator.h b/ThirdParty/RootMinimizers/src/Math/GSLDerivator.h new file mode 100644 index 00000000000..5c99a4a2700 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLDerivator.h @@ -0,0 +1,176 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class Derivator +// +// class for calculating Derivative of functions +// +// Created by: moneta at Sat Nov 13 14:46:00 2004 +// +// Last update: Sat Nov 13 14:46:00 2004 +// +#ifndef ROOT_Math_GSLDerivator +#define ROOT_Math_GSLDerivator + +/** +@defgroup Deriv Numerical Differentiation +*/ + +#include "Math/GSLFunctionAdapter.h" +#include "GSLFunctionWrapper.h" + + +#include "Math/IFunctionfwd.h" +#include "Math/IFunction.h" + +namespace ROOT { +namespace Math { + + +class GSLFunctionWrapper; + + +/** + Class for computing numerical derivative of a function based on the GSL numerical algorithm + This class is implemented using the numerical derivatives algorithms provided by GSL + (see <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Numerical-Differentiation.html">GSL Online Manual</A> ). + + @ingroup Deriv +*/ + +class GSLDerivator { + +public: + /** + Default Constructor of a GSLDerivator class based on GSL numerical differentiation algorithms + */ + GSLDerivator() : fStatus(0), fResult(0), fError(0) {} + + /// destructor (no operations) + virtual ~GSLDerivator() {} + +// // disable copying +// private: + +// GSLDerivator(const GSLDerivator &); +// GSLDerivator & operator = (const GSLDerivator &); + +// public: + + + + /** + Set the function for calculating the derivatives. + The function must implement the ROOT::Math::IGenFunction signature + */ + void SetFunction(const IGenFunction &f); + + /** + Set the function f for evaluating the derivative using a GSL function pointer type + @param f : free function pointer of the GSL required type + @param p : pointer to the object carrying the function state + (for example the function object itself) + */ + void SetFunction( GSLFuncPointer f, void * p = 0); + + /** + Computes the numerical derivative at a point x using an adaptive central + difference algorithm with a step size h. + */ + double EvalCentral( double x, double h); + + /** + Computes the numerical derivative at a point x using an adaptive forward + difference algorithm with a step size h. + The function is evaluated only at points greater than x and at x itself. + */ + double EvalForward( double x, double h); + + /** + Computes the numerical derivative at a point x using an adaptive backward + difference algorithm with a step size h. + The function is evaluated only at points less than x and at x itself. + */ + double EvalBackward( double x, double h); + + /** @name --- Static methods --- **/ + + /** + Computes the numerical derivative of a function f at a point x using an adaptive central + difference algorithm with a step size h + */ + static double EvalCentral(const IGenFunction & f, double x, double h); + + + /** + Computes the numerical derivative of a function f at a point x using an adaptive forward + difference algorithm with a step size h. + The function is evaluated only at points greater than x and at x itself + */ + static double EvalForward(const IGenFunction & f, double x, double h); + + /** + Computes the numerical derivative of a function f at a point x using an adaptive backward + difference algorithm with a step size h. + The function is evaluated only at points less than x and at x itself + */ + + static double EvalBackward(const IGenFunction & f, double x, double h); + + + + /** + return the error status of the last integral calculation + */ + int Status() const; + + /** + return the result of the last derivative calculation + */ + double Result() const; + + /** + return the estimate of the absolute error of the last derivative calculation + */ + double Error() const; + + +private: + + int fStatus; + double fResult; + double fError; + + GSLFunctionWrapper fFunction; + +}; + + + + +} // namespace Math +} // namespace ROOT + + +#endif /* ROOT_Math_GSLDerivator */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLFunctionWrapper.h b/ThirdParty/RootMinimizers/src/Math/GSLFunctionWrapper.h new file mode 100644 index 00000000000..0dc09054801 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLFunctionWrapper.h @@ -0,0 +1,150 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLFunctionWrapper +// +// Created by: moneta at Sat Nov 13 14:54:41 2004 +// +// Last update: Sat Nov 13 14:54:41 2004 +// +#ifndef ROOT_Math_GSLFunctionWrapper +#define ROOT_Math_GSLFunctionWrapper + +#include "gsl/gsl_math.h" + +#include "Math/GSLFunctionAdapter.h" + +#include <cassert> + +namespace ROOT { +namespace Math { + + + +typedef double ( * GSLFuncPointer ) ( double, void *); +typedef void ( * GSLFdfPointer ) ( double, void *, double *, double *); + + +/** + Wrapper class to the gsl_function C structure. + This class to fill the GSL C structure gsl_function with + the C++ function objcet. + Use the class ROOT::Math::GSLFunctionAdapter to adapt the + C++ function object to the right signature (function pointer type) + requested by GSL +*/ +class GSLFunctionWrapper { + +public: + + GSLFunctionWrapper() + { + fFunc.function = 0; + fFunc.params = 0; + } + + /// set in the GSL C struct the pointer to the function evaluation + void SetFuncPointer( GSLFuncPointer f) { fFunc.function = f; } + + /// set in the GSL C struct the extra-object pointer + void SetParams ( void * p) { fFunc.params = p; } + + /// fill the GSL C struct from a generic C++ callable object + /// implementing operator() + template<class FuncType> + void SetFunction(const FuncType &f) { + const void * p = &f; + assert (p != 0); + SetFuncPointer(&GSLFunctionAdapter<FuncType >::F); + SetParams(const_cast<void *>(p)); + } + + gsl_function * GetFunc() { return &fFunc; } + + GSLFuncPointer FunctionPtr() { return fFunc.function; } + + // evaluate the function + double operator() (double x) { return GSL_FN_EVAL(&fFunc, x); } + + /// check if function is valid (has been set) + bool IsValid() { + return (fFunc.function != 0) ? true : false; + } + +private: + gsl_function fFunc; + + +}; + + + /** + class to wrap a gsl_function_fdf (with derivatives) + */ + class GSLFunctionDerivWrapper { + + public: + + GSLFunctionDerivWrapper() + { + fFunc.f = 0; + fFunc.df = 0; + fFunc.fdf = 0; + fFunc.params = 0; + } + + + void SetFuncPointer( GSLFuncPointer f) { fFunc.f = f; } + void SetDerivPointer( GSLFuncPointer f) { fFunc.df = f; } + void SetFdfPointer( GSLFdfPointer f) { fFunc.fdf = f; } + void SetParams ( void * p) { fFunc.params = p; } + + + gsl_function_fdf * GetFunc() { return &fFunc; } + + // evaluate the function and derivatives + double operator() (double x) { return GSL_FN_FDF_EVAL_F(&fFunc, x); } + + double Derivative (double x) { return GSL_FN_FDF_EVAL_DF(&fFunc, x); } + + void Fdf(double x, double & f, double & df) { + return GSL_FN_FDF_EVAL_F_DF(&fFunc, x, &f, &df); + } + + /// check if function is valid (has been set) + bool IsValid() { + return (fFunc.f != 0 ) ? true : false; + } + + private: + gsl_function_fdf fFunc; + + }; + + + +} // namespace Math +} // namespace ROOT + +#endif /* ROOT_Math_GSLFunctionWrapper */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMinimizer.cxx b/ThirdParty/RootMinimizers/src/Math/GSLMinimizer.cxx new file mode 100644 index 00000000000..6404e7b5d74 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMinimizer.cxx @@ -0,0 +1,400 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Tue Dec 19 15:41:39 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Implementation file for class GSLMinimizer + +#include "Math/GSLMinimizer.h" + +#include "GSLMultiMinimizer.h" + +#include "Math/MultiNumGradFunction.h" +#include "Math/FitMethodFunction.h" + +#include "Math/MinimTransformFunction.h" + +#include <cassert> + +#include <iostream> +#include <iomanip> +#include <cmath> +#include <algorithm> +#include <functional> +#include <ctype.h> // need to use c version of tolower defined here +#include <limits> + +namespace ROOT { + + namespace Math { + + +GSLMinimizer::GSLMinimizer( ROOT::Math::EGSLMinimizerType type) : + fDim(0), + fObjFunc(0), + fMinVal(0) +{ + // Constructor implementation : create GSLMultiMin wrapper object + //std::cout << "create GSL Minimizer of type " << type << std::endl; + + fGSLMultiMin = new GSLMultiMinimizer((ROOT::Math::EGSLMinimizerType) type); + fValues.reserve(10); + fNames.reserve(10); + fSteps.reserve(10); + + fLSTolerance = 0.1; // line search tolerance (use fixed) + int niter = ROOT::Math::MinimizerOptions::DefaultMaxIterations(); + if (niter <=0 ) niter = 1000; + SetMaxIterations(niter); + SetPrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel()); +} + +GSLMinimizer::GSLMinimizer( const char * type) : + fDim(0), + fObjFunc(0), + fMinVal(0) +{ + // Constructor implementation from a string + std::string algoname(type); + std::transform(algoname.begin(), algoname.end(), algoname.begin(), (int(*)(int)) tolower ); + + ROOT::Math::EGSLMinimizerType algo = kVectorBFGS2; // default value + + if (algoname == "conjugatefr") algo = kConjugateFR; + if (algoname == "conjugatepr") algo = kConjugatePR; + if (algoname == "bfgs") algo = kVectorBFGS; + if (algoname == "bfgs2") algo = kVectorBFGS2; + if (algoname == "steepestdescent") algo = kSteepestDescent; + + + //std::cout << "create GSL Minimizer of type " << algo << std::endl; + + fGSLMultiMin = new GSLMultiMinimizer(algo); + fValues.reserve(10); + fNames.reserve(10); + fSteps.reserve(10); + + fLSTolerance = 0.1; // use 10**-4 + int niter = ROOT::Math::MinimizerOptions::DefaultMaxIterations(); + if (niter <=0 ) niter = 1000; + SetMaxIterations(niter); + SetPrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel()); +} + + +GSLMinimizer::~GSLMinimizer () { + assert(fGSLMultiMin != 0); + delete fGSLMultiMin; + if (fObjFunc) delete fObjFunc; +} + +bool GSLMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) { + // set variable in minimizer - support only free variables + // no transformation implemented - so far + if (ivar > fValues.size() ) return false; + if (ivar == fValues.size() ) { + fValues.push_back(val); + fNames.push_back(name); + fSteps.push_back(step); + fVarTypes.push_back(kDefault); + } + else { + fValues[ivar] = val; + fNames[ivar] = name; + fSteps[ivar] = step; + fVarTypes[ivar] = kDefault; + + // remove bounds if needed + std::map<unsigned int, std::pair<double, double> >::iterator iter = fBounds.find(ivar); + if ( iter != fBounds.end() ) fBounds.erase (iter); + + } + + return true; +} + +bool GSLMinimizer::SetLowerLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower) { + //MATH_WARN_MSGVAL("GSLMinimizer::SetLowerLimitedVariable","Ignore lower limit on variable ",ivar); + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( lower, lower); + fVarTypes[ivar] = kLowBound; + return true; +} +bool GSLMinimizer::SetUpperLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double upper ) { + //MATH_WARN_MSGVAL("GSLMinimizer::SetUpperLimitedVariable","Ignore upper limit on variable ",ivar); + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( upper, upper); + fVarTypes[ivar] = kUpBound; + return true; +} + +bool GSLMinimizer::SetLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower, double upper) { + //MATH_WARN_MSGVAL("GSLMinimizer::SetLimitedVariable","Ignore bounds on variable ",ivar); + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( lower, upper); + fVarTypes[ivar] = kBounds; + return true; +} + +bool GSLMinimizer::SetFixedVariable(unsigned int ivar , const std::string & name , double val ) { + /// set fixed variable (override if minimizer supports them ) + bool ret = SetVariable(ivar, name, val, 0.); + if (!ret) return false; + fVarTypes[ivar] = kFix; + return true; +} + + +bool GSLMinimizer::SetVariableValue(unsigned int ivar, double val) { + // set variable value in minimizer + // no change to transformation or variable status + if (ivar > fValues.size() ) return false; + fValues[ivar] = val; + return true; +} + +bool GSLMinimizer::SetVariableValues( const double * x) { + // set all variable values in minimizer + if (x == 0) return false; + std::copy(x,x+fValues.size(), fValues.begin() ); + return true; +} + + +void GSLMinimizer::SetFunction(const ROOT::Math::IMultiGenFunction & func) { + // set the function to minimizer + // need to calculate numerically the derivatives: do via class MultiNumGradFunction + fObjFunc = new MultiNumGradFunction( func); + fDim = fObjFunc->NDim(); +} + +void GSLMinimizer::SetFunction(const ROOT::Math::IMultiGradFunction & func) { + // set the function to minimizer + fObjFunc = dynamic_cast<const ROOT::Math::IMultiGradFunction *>( func.Clone()); + assert(fObjFunc != 0); + fDim = fObjFunc->NDim(); +} + +unsigned int GSLMinimizer::NCalls() const { + // return numbr of function calls + // if method support + const ROOT::Math::MultiNumGradFunction * fnumgrad = dynamic_cast<const ROOT::Math::MultiNumGradFunction *>(fObjFunc); + if (fnumgrad) return fnumgrad->NCalls(); + const ROOT::Math::FitMethodGradFunction * ffitmethod = dynamic_cast<const ROOT::Math::FitMethodGradFunction *>(fObjFunc); + if (ffitmethod) return ffitmethod->NCalls(); + // not supported in the other case + return 0; +} + +bool GSLMinimizer::Minimize() { + // set initial parameters of the minimizer + + if (fGSLMultiMin == 0) return false; + if (fObjFunc == 0) { + MATH_ERROR_MSG("GSLMinimizer::Minimize","Function has not been set"); + return false; + } + + unsigned int npar = fValues.size(); + if (npar == 0 || npar < fObjFunc->NDim() ) { + MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Wrong number of parameters",npar); + return false; + } + + // use a global step size = modules of step vectors + double stepSize = 0; + for (unsigned int i = 0; i < fSteps.size(); ++i) + stepSize += fSteps[i]*fSteps[i]; + stepSize = std::sqrt(stepSize); + + const double eps = std::numeric_limits<double>::epsilon(); + if (stepSize < eps) { + MATH_ERROR_MSGVAL("GSLMinimizer::Minimize","Step size is too small",stepSize); + return false; + } + + // check if a transformation is needed + bool doTransform = (fBounds.size() > 0); + unsigned int ivar = 0; + while (!doTransform && ivar < fVarTypes.size() ) { + doTransform = (fVarTypes[ivar++] != kDefault ); + } + + + std::vector<double> startValues(fValues.begin(), fValues.end() ); + + MinimTransformFunction * trFunc = 0; + + // in case of transformation wrap objective function in a new transformation function + // and transform from external variables to internals one + if (doTransform) { + trFunc = new MinimTransformFunction ( fObjFunc, fVarTypes, fValues, fBounds ); + trFunc->InvTransformation(&fValues.front(), &startValues[0]); + startValues.resize( trFunc->NDim() ); + fObjFunc = trFunc; + } + +// std::cout << " f has transform " << doTransform << " " << fBounds.size() << " " << startValues.size() << " ndim " << fObjFunc->NDim() << std::endl; std::cout << "InitialValues external : "; +// for (int i = 0; i < fValues.size(); ++i) std::cout << fValues[i] << " "; +// std::cout << "\n"; +// std::cout << "InitialValues internal : "; +// for (int i = 0; i < startValues.size(); ++i) std::cout << startValues[i] << " "; +// std::cout << "\n"; + + + // set parameters in internal GSL minimization class + fGSLMultiMin->Set(*fObjFunc, &startValues.front(), stepSize, fLSTolerance ); + + + int debugLevel = PrintLevel(); + + if (debugLevel >=1 ) std::cout <<"Minimize using GSLMinimizer " << fGSLMultiMin->Name() << std::endl; + + + //std::cout <<"print Level " << debugLevel << std::endl; + //debugLevel = 3; + + // start iteration + unsigned int iter = 0; + int status; + bool minFound = false; + bool iterFailed = false; + do { + status = fGSLMultiMin->Iterate(); + if (status) { + iterFailed = true; + break; + } + + status = fGSLMultiMin->TestGradient( Tolerance() ); + if (status == GSL_SUCCESS) { + minFound = true; + } + + if (debugLevel >=2) { + std::cout << "----------> Iteration " << std::setw(4) << iter; + int pr = std::cout.precision(18); + std::cout << " FVAL = " << fGSLMultiMin->Minimum() << std::endl; + std::cout.precision(pr); + if (debugLevel >=3) { + std::cout << " Parameter Values : "; + const double * xtmp = fGSLMultiMin->X(); + std::cout << std::endl; + if (trFunc != 0 ) { + xtmp = trFunc->Transformation(xtmp); + } + for (unsigned int i = 0; i < NDim(); ++i) { + std::cout << " " << fNames[i] << " = " << xtmp[i]; + // avoid nan + // if (std::isnan(xtmp[i])) status = -11; + } + std::cout << std::endl; + } + } + + + iter++; + + + } + while (status == GSL_CONTINUE && iter < MaxIterations() ); + + // save state with values and function value + double * x = fGSLMultiMin->X(); + if (x == 0) return false; + + // check to see if a transformation need to be applied + if (trFunc != 0) { + const double * xtrans = trFunc->Transformation(x); + assert(fValues.size() == trFunc->NTot() ); + assert( trFunc->NTot() == NDim() ); + std::copy(xtrans, xtrans + trFunc->NTot(), fValues.begin() ); + } + else { + // case of no transformation applied + assert( fValues.size() == NDim() ); + std::copy(x, x + NDim(), fValues.begin() ); + } + + fMinVal = fGSLMultiMin->Minimum(); + + fStatus = status; + + + if (minFound) { + if (debugLevel >=1 ) { + std::cout << "GSLMinimizer: Minimum Found" << std::endl; + int pr = std::cout.precision(18); + std::cout << "FVAL = " << fMinVal << std::endl; + std::cout.precision(pr); +// std::cout << "Edm = " << fState.Edm() << std::endl; + std::cout << "Niterations = " << iter << std::endl; + unsigned int ncalls = NCalls(); + if (ncalls) std::cout << "NCalls = " << ncalls << std::endl; + for (unsigned int i = 0; i < fDim; ++i) + std::cout << fNames[i] << "\t = " << fValues[i] << std::endl; + } + return true; + } + else { + if (debugLevel >= -1 ) { + std::cout << "GSLMinimizer: Minimization did not converge" << std::endl; + if (iterFailed) { + if (status == GSL_ENOPROG) // case status 27 + std::cout << "\t Iteration is not making progress towards solution" << std::endl; + else + std::cout << "\t Iteration failed with status " << status << std::endl; + + if (debugLevel >= 1) { + double * g = fGSLMultiMin->Gradient(); + double dg2 = 0; + for (unsigned int i = 0; i < fDim; ++i) dg2 += g[i] * g[1]; + std::cout << "Grad module is " << std::sqrt(dg2) << std::endl; + for (unsigned int i = 0; i < fDim; ++i) + std::cout << fNames[i] << "\t = " << fValues[i] << std::endl; + std::cout << "FVAL = " << fMinVal << std::endl; +// std::cout << "Edm = " << fState.Edm() << std::endl; + std::cout << "Niterations = " << iter << std::endl; + } + } + } + return false; + } + return false; +} + +const double * GSLMinimizer::MinGradient() const { + // return gradient (internal values) + return fGSLMultiMin->Gradient(); +} + +const MinimTransformFunction * GSLMinimizer::TransformFunction() const { + return dynamic_cast<const MinimTransformFunction *>(fObjFunc); +} + + } // end namespace Math + +} // end namespace ROOT + diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMinimizer1D.cxx b/ThirdParty/RootMinimizers/src/Math/GSLMinimizer1D.cxx new file mode 100644 index 00000000000..ede23fcc7c7 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMinimizer1D.cxx @@ -0,0 +1,222 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + /********************************************************************** + * * + * Copyright (c) 2004 moneta, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Implementation file for class GSLMinimizer1D +// +// Created by: moneta at Wed Dec 1 15:04:51 2004 +// +// Last update: Wed Dec 1 15:04:51 2004 +// + +#include <assert.h> + +#include "Math/GSLMinimizer1D.h" +#include "Math/Error.h" + +#include "GSLFunctionWrapper.h" +#include "GSL1DMinimizerWrapper.h" + + +#include "gsl/gsl_min.h" +#include "gsl/gsl_errno.h" + +#include <iostream> +#include <cmath> + +namespace ROOT { + +namespace Math { + + +GSLMinimizer1D::GSLMinimizer1D(Minim1D::Type type) : + fXmin(0), fXlow(0), fXup(0), fMin(0), fLow(0), fUp(0), + fIter(0), fStatus(-1), fIsSet(false), + fMinimizer(0), fFunction(0) +{ + // construct a minimizer passing the algorithm type as an enumeration + + const gsl_min_fminimizer_type* T = 0 ; + switch ( type ) + { + case Minim1D::kGOLDENSECTION : + T = gsl_min_fminimizer_goldensection; + break ; + case Minim1D::kBRENT : + T = gsl_min_fminimizer_brent; + break ; + default : + // default case is brent + T = gsl_min_fminimizer_brent; + break ; + } + + fMinimizer = new GSL1DMinimizerWrapper(T); + fFunction = new GSLFunctionWrapper(); + +} + +GSLMinimizer1D::~GSLMinimizer1D() +{ + // destructor: clean up minimizer and function pointers + + if (fMinimizer) delete fMinimizer; + if (fFunction) delete fFunction; +} + +GSLMinimizer1D::GSLMinimizer1D(const GSLMinimizer1D &): IMinimizer1D() +{ + // dummy copy ctr +} + +GSLMinimizer1D & GSLMinimizer1D::operator = (const GSLMinimizer1D &rhs) +{ + // dummy operator = + if (this == &rhs) return *this; // time saving self-test + return *this; +} + +void GSLMinimizer1D::SetFunction( GSLFuncPointer f, void * p, double xmin, double xlow, double xup) { + // set the function to be minimized + assert(fFunction); + assert(fMinimizer); + fXlow = xlow; + fXup = xup; + fXmin = xmin; + fFunction->SetFuncPointer( f ); + fFunction->SetParams( p ); + +#ifdef DEBUG + std::cout << " [ "<< xlow << " , " << xup << " ]" << std::endl; +#endif + + int status = gsl_min_fminimizer_set( fMinimizer->Get(), fFunction->GetFunc(), xmin, xlow, xup); + if (status != GSL_SUCCESS) + std::cerr <<"GSLMinimizer1D: Error: Interval [ "<< xlow << " , " << xup << " ] does not contain a minimum" << std::endl; + + + fIsSet = true; + fStatus = -1; + return; +} + +int GSLMinimizer1D::Iterate() { + // perform an iteration and update values + if (!fIsSet) { + std::cerr << "GSLMinimizer1D- Error: Function has not been set in Minimizer" << std::endl; + return -1; + } + + int status = gsl_min_fminimizer_iterate(fMinimizer->Get()); + // update values + fXmin = gsl_min_fminimizer_x_minimum(fMinimizer->Get() ); + fMin = gsl_min_fminimizer_f_minimum(fMinimizer->Get() ); + // update interval values + fXlow = gsl_min_fminimizer_x_lower(fMinimizer->Get() ); + fXup = gsl_min_fminimizer_x_upper(fMinimizer->Get() ); + fLow = gsl_min_fminimizer_f_lower(fMinimizer->Get() ); + fUp = gsl_min_fminimizer_f_upper(fMinimizer->Get() ); + return status; +} + +double GSLMinimizer1D::XMinimum() const { + // return x value at function minimum + return fXmin; +} + +double GSLMinimizer1D::XLower() const { + // return lower x value of bracketing interval + return fXlow; +} + +double GSLMinimizer1D::XUpper() const { + // return upper x value of bracketing interval + return fXup; +} + +double GSLMinimizer1D::FValMinimum() const { + // return function value at minimum + return fMin; +} + +double GSLMinimizer1D::FValLower() const { + // return function value at x lower + return fLow; +} + +double GSLMinimizer1D::FValUpper() const { + // return function value at x upper + return fUp; +} + +const char * GSLMinimizer1D::Name() const { + // return name of minimization algorithm + return gsl_min_fminimizer_name(fMinimizer->Get() ); +} + +bool GSLMinimizer1D::Minimize (int maxIter, double absTol, double relTol) +{ + // find the minimum via multiple iterations + fStatus = -1; + int iter = 0; + int status = 0; + do { + iter++; + status = Iterate(); + if (status != GSL_SUCCESS) { + MATH_ERROR_MSG("GSLMinimizer1D::Minimize","error returned when performing an iteration"); + fStatus = status; + return false; + } + +#ifdef DEBUG + std::cout << "Min1D - iteration " << iter << " interval : [ " << fXlow << " , " << fXup << " ] min = " << fXmin + << " fmin " << fMin << " f(a) " << fLow << " f(b) " << fUp << std::endl; +#endif + + + status = TestInterval(fXlow, fXup, absTol, relTol); + if (status == GSL_SUCCESS) { + fIter = iter; + fStatus = status; + return true; + } + } + while (status == GSL_CONTINUE && iter < maxIter); + if (status == GSL_CONTINUE) { + double tol = std::abs(fXup-fXlow); + MATH_INFO_MSGVAL("GSLMinimizer1D::Minimize","exceeded max iterations, reached tolerance is not sufficient",tol); + } + fStatus = status; + return false; +} + + +int GSLMinimizer1D::TestInterval( double xlow, double xup, double epsAbs, double epsRel) { +// static function to test interval + return gsl_min_test_interval(xlow, xup, epsAbs, epsRel); +} + +} // end namespace Math + +} // end namespace ROOT + diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h b/ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h new file mode 100644 index 00000000000..28b029482af --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h @@ -0,0 +1,209 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Dec 20 17:26:06 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLMultiFit + +#ifndef ROOT_Math_GSLMultiFit +#define ROOT_Math_GSLMultiFit + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_multifit_nlin.h" +#include "gsl/gsl_blas.h" +#include "GSLMultiFitFunctionWrapper.h" + +#include "Math/IFunction.h" +#include <string> +#include <cassert> + + +namespace ROOT { + + namespace Math { + + +/** + GSLMultiFit, internal class for implementing GSL non linear least square GSL fitting + + @ingroup MultiMin +*/ +class GSLMultiFit { + +public: + + /** + Default constructor + No need to specify the type so far since only one solver exists so far + */ + GSLMultiFit (const gsl_multifit_fdfsolver_type * type = 0) : + fSolver(0), + fVec(0), + fCov(0), + fType(type) + { + if (fType == 0) fType = gsl_multifit_fdfsolver_lmsder; // default value + } + + /** + Destructor (no operations) + */ + ~GSLMultiFit () { + if (fSolver) gsl_multifit_fdfsolver_free(fSolver); + if (fVec != 0) gsl_vector_free(fVec); + if (fCov != 0) gsl_matrix_free(fCov); + } + +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + GSLMultiFit(const GSLMultiFit &) {} + + /** + Assignment operator + */ + GSLMultiFit & operator = (const GSLMultiFit & rhs) { + if (this == &rhs) return *this; // time saving self-test + return *this; + } + + +public: + + /// create the minimizer from the type and size of number of fitting points and number of parameters + void CreateSolver(unsigned int npoints, unsigned int npar) { + if (fSolver) gsl_multifit_fdfsolver_free(fSolver); + fSolver = gsl_multifit_fdfsolver_alloc(fType, npoints, npar); + } + + /// set the solver parameters + template<class Func> + int Set(const std::vector<Func> & funcVec, const double * x) { + // create a vector of the fit contributions + // create function wrapper from an iterator of functions + unsigned int npts = funcVec.size(); + if (npts == 0) return -1; + + unsigned int npar = funcVec[0].NDim(); + //typedef typename std::vector<Func> FuncVec; + //FuncIt funcIter = funcVec.begin(); + fFunc.SetFunction(funcVec, npts, npar); + // create solver object + CreateSolver( npts, npar ); + // set initial values + if (fVec != 0) gsl_vector_free(fVec); + fVec = gsl_vector_alloc( npar ); + std::copy(x,x+npar, fVec->data); + assert(fSolver != 0); + return gsl_multifit_fdfsolver_set(fSolver, fFunc.GetFunc(), fVec); + } + + std::string Name() const { + if (fSolver == 0) return "undefined"; + return std::string(gsl_multifit_fdfsolver_name(fSolver) ); + } + + int Iterate() { + if (fSolver == 0) return -1; + return gsl_multifit_fdfsolver_iterate(fSolver); + } + + /// parameter values at the minimum + const double * X() const { + if (fSolver == 0) return 0; + gsl_vector * x = gsl_multifit_fdfsolver_position(fSolver); + return x->data; + } + + /// gradient value at the minimum + const double * Gradient() const { + if (fSolver == 0) return 0; + gsl_multifit_gradient(fSolver->J, fSolver->f,fVec); + return fVec->data; + } + + /// return covariance matrix of the parameters + const double * CovarMatrix() const { + if (fSolver == 0) return 0; + if (fCov != 0) gsl_matrix_free(fCov); + unsigned int npar = fSolver->fdf->p; + fCov = gsl_matrix_alloc( npar, npar ); + static double kEpsrel = 0.0001; + int ret = gsl_multifit_covar(fSolver->J, kEpsrel, fCov); + if (ret != GSL_SUCCESS) return 0; + return fCov->data; + } + + /// test gradient (ask from solver gradient vector) + int TestGradient(double absTol) const { + if (fSolver == 0) return -1; + Gradient(); + return gsl_multifit_test_gradient( fVec, absTol); + } + + /// test using abs and relative tolerance + /// |dx| < absTol + relTol*|x| for every component + int TestDelta(double absTol, double relTol) const { + if (fSolver == 0) return -1; + return gsl_multifit_test_delta(fSolver->dx, fSolver->x, absTol, relTol); + } + + // calculate edm 1/2 * ( grad^T * Cov * grad ) + double Edm() const { + // product C * g + double edm = -1; + const double * g = Gradient(); + if (g == 0) return edm; + const double * c = CovarMatrix(); + if (c == 0) return edm; + gsl_vector * tmp = gsl_vector_alloc( fSolver->fdf->p ); + int status = gsl_blas_dgemv(CblasNoTrans, 1.0, fCov, fVec, 0.,tmp); + if (status == 0) status |= gsl_blas_ddot(fVec, tmp, &edm); + gsl_vector_free(tmp); + if (status != 0) return -1; + // need to divide by 2 ?? + return 0.5*edm; + + } + + +private: + + GSLMultiFitFunctionWrapper fFunc; + gsl_multifit_fdfsolver * fSolver; + // cached vector to avoid re-allocating every time a new one + mutable gsl_vector * fVec; + mutable gsl_matrix * fCov; + const gsl_multifit_fdfsolver_type * fType; + +}; + + } // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_GSLMultiFit */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionAdapter.h b/ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionAdapter.h new file mode 100644 index 00000000000..ee6b88a321a --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionAdapter.h @@ -0,0 +1,131 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, Dec 2006 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLMultiMinFunctionAdapter +// +// Generic adapter for gsl_multimin_function signature +// usable for any c++ class which defines operator( ) +// +// Created by: Lorenzo Moneta at Fri Nov 12 16:58:51 2004 +// +// Last update: Fri Nov 12 16:58:51 2004 +// +#ifndef ROOT_Math_GSLMultiFitFunctionAdapter +#define ROOT_Math_GSLMultiFitFunctionAdapter + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" + +#include <cassert> + +#include <iostream> + +namespace ROOT { +namespace Math { + + + + /** + Class for adapting a C++ functor class to C function pointers used by GSL MultiFit + Algorithm + The templated C++ function class must implement: + + <em> double operator( const double * x)</em> + and if the derivatives are required: + <em> void Gradient( const double * x, double * g)</em> + and + <em> void FdF( const double * x, double &f, double * g)</em> + + This class defines static methods with will be used to fill the + \a gsl_multimin_function and + \a gsl_multimin_function_fdf structs used by GSL. + See for examples the + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Providing-a-function-to-minimize.html#Providing-a-function-to-minimize">GSL online manual</A> + + @ingroup MultiMin + */ + + +template<class FuncVector> +class GSLMultiFitFunctionAdapter { + +public: + + static int F( const gsl_vector * x, void * p, gsl_vector * f ) { + // p is a pointer to an iterator of functions + unsigned int n = f->size; + // need to copy iterator otherwise next time the function is called it wont work + FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) ); + if (n == 0) return -1; + for (unsigned int i = 0; i < n ; ++i) { + gsl_vector_set(f, i, (funcVec[i])(x->data) ); + } + return 0; + } + + + static int Df( const gsl_vector * x, void * p, gsl_matrix * h) { + + // p is a pointer to an iterator of functions + unsigned int n = h->size1; + unsigned int npar = h->size2; + if (n == 0) return -1; + if (npar == 0) return -2; + FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) ); + for (unsigned int i = 0; i < n ; ++i) { + double * g = (h->data)+i*npar; //pointer to start of i-th row + assert ( npar == (funcVec[i]).NDim() ); + (funcVec[i]).Gradient(x->data, g); + } + return 0; + } + + /// evaluate derivative and function at the same time + static int FDf( const gsl_vector * x, void * p, gsl_vector * f, gsl_matrix * h) { + // should be implemented in the function + // p is a pointer to an iterator of functions + unsigned int n = h->size1; + unsigned int npar = h->size2; + if (n == 0) return -1; + if (npar == 0) return -2; + FuncVector & funcVec = *( reinterpret_cast< FuncVector *> (p) ); + assert ( f->size == n); + for (unsigned int i = 0; i < n ; ++i) { + assert ( npar == (funcVec[i]).NDim() ); + double fval = 0; + double * g = (h->data)+i*npar; //pointer to start of i-th row + (funcVec[i]).FdF(x->data, fval, g); + gsl_vector_set(f, i, fval ); + } + return 0; + } + +}; + + +} // namespace Math +} // namespace ROOT + + +#endif /* ROOT_Math_GSLMultiMinFunctionAdapter */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionWrapper.h b/ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionWrapper.h new file mode 100644 index 00000000000..4eadcdcc1a5 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMultiFitFunctionWrapper.h @@ -0,0 +1,100 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta Dec 2006 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLMultiMinFunctionWrapper +// +// Created by: moneta at Sat Nov 13 14:54:41 2004 +// +// Last update: Sat Nov 13 14:54:41 2004 +// +#ifndef ROOT_Math_GSLMultiFitFunctionWrapper +#define ROOT_Math_GSLMultiFitFunctionWrapper + +#include "gsl/gsl_multifit.h" + +#include "GSLMultiFitFunctionAdapter.h" + + +#include <cassert> + +namespace ROOT { +namespace Math { + + + + typedef double ( * GSLMultiFitFPointer ) ( const gsl_vector *, void *, gsl_vector *); + typedef void ( * GSLMultiFitDfPointer ) ( const gsl_vector *, void *, gsl_matrix *); + typedef void ( * GSLMultiFitFdfPointer ) ( const gsl_vector *, void *, gsl_vector *, gsl_matrix *); + + +/** + wrapper to a multi-dim function withtout derivatives for multi-dimensional + minimization algorithm + + @ingroup MultiMin +*/ + +class GSLMultiFitFunctionWrapper { + +public: + + GSLMultiFitFunctionWrapper() + { + fFunc.f = 0; + fFunc.df = 0; + fFunc.fdf = 0; + fFunc.n = 0; + fFunc.p = 0; + fFunc.params = 0; + } + + + /// Fill gsl function structure from a C++ function iterator and size and number of residuals + template<class FuncVector> + void SetFunction(const FuncVector & f, unsigned int nres, unsigned int npar ) { + const void * p = &f; + assert (p != 0); + fFunc.f = &GSLMultiFitFunctionAdapter<FuncVector >::F; + fFunc.df = &GSLMultiFitFunctionAdapter<FuncVector >::Df; + fFunc.fdf = &GSLMultiFitFunctionAdapter<FuncVector >::FDf; + fFunc.n = nres; + fFunc.p = npar; + fFunc.params = const_cast<void *>(p); + } + + gsl_multifit_function_fdf * GetFunc() { return &fFunc; } + + + private: + + gsl_multifit_function_fdf fFunc; + +}; + + + +} // namespace Math +} // namespace ROOT + +#endif /* ROOT_Math_GSLMultiMinFunctionWrapper */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionAdapter.h b/ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionAdapter.h new file mode 100644 index 00000000000..0be4d687984 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionAdapter.h @@ -0,0 +1,99 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, 12/2006 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLMultiMinFunctionAdapter +// +// Generic adapter for gsl_multimin_function signature +// usable for any c++ class which defines operator( ) +// +// Created by: Lorenzo Moneta at Fri Nov 12 16:58:51 2004 +// +// Last update: Fri Nov 12 16:58:51 2004 +// +#ifndef ROOT_Math_GSLMultiMinFunctionAdapter +#define ROOT_Math_GSLMultiMinFunctionAdapter + +#include "gsl/gsl_vector.h" + +namespace ROOT { +namespace Math { + + + + + /** + Class for adapting any multi-dimension C++ functor class to C function pointers used by + GSL MultiMin algorithms. + The templated C++ function class must implement: + + <em> double operator( const double * x)</em> + and if the derivatives are required: + <em> void Gradient( const double * x, double * g)</em> + + This class defines static methods with will be used to fill the + \a gsl_multimin_function and + \a gsl_multimin_function_fdf structs used by GSL. + See for examples the + <A HREF="http://www.gnu.org/software/gsl/manual/html_node/Providing-a-function-to-minimize.html#Providing-a-function-to-minimize">GSL online manual</A> + + @ingroup MultiMin + + */ + + + template<class UserFunc> + struct GSLMultiMinFunctionAdapter { + + static double F( const gsl_vector * x, void * p) { + + UserFunc * function = reinterpret_cast< UserFunc *> (p); + // get pointer to data from gsl_vector + return (*function)( x->data ); + } + + + static void Df( const gsl_vector * x, void * p, gsl_vector * g) { + + UserFunc * function = reinterpret_cast< UserFunc *> (p); + (*function).Gradient( x->data, g->data ); + + } + + static void Fdf( const gsl_vector * x, void * p, double *f, gsl_vector * g ) { + + UserFunc * function = reinterpret_cast< UserFunc *> (p); +// *f = (*function) ( x ); +// *df = (*function).Gradient( x ); + + (*function).FdF( x->data, *f, g->data); + } + + }; + + +} // namespace Math +} // namespace ROOT + + +#endif /* ROOT_Math_GSLMultiMinFunctionAdapter */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionWrapper.h b/ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionWrapper.h new file mode 100644 index 00000000000..e70e5b83d6e --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMultiMinFunctionWrapper.h @@ -0,0 +1,163 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, 12/2006 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLMultiMinFunctionWrapper +// +// Created by: moneta at Sat Nov 13 14:54:41 2004 +// +// Last update: Sat Nov 13 14:54:41 2004 +// +#ifndef ROOT_Math_GSLMultiMinFunctionWrapper +#define ROOT_Math_GSLMultiMinFunctionWrapper + +#include "gsl/gsl_multimin.h" + +#include "GSLMultiMinFunctionAdapter.h" + + +#include <cassert> + +namespace ROOT { +namespace Math { + + + + typedef double ( * GSLMultiMinFuncPointer ) ( const gsl_vector *, void *); + typedef void ( * GSLMultiMinDfPointer ) ( const gsl_vector *, void *, gsl_vector *); + typedef void ( * GSLMultiMinFdfPointer ) ( const gsl_vector *, void *, double *, gsl_vector *); + + +/** + wrapper to a multi-dim function withtout derivatives for multi-dimensional + minimization algorithm + + @ingroup MultiMin +*/ + +class GSLMultiMinFunctionWrapper { + +public: + + GSLMultiMinFunctionWrapper() + { + fFunc.f = 0; + fFunc.n = 0; + fFunc.params = 0; + } + + void SetFuncPointer( GSLMultiMinFuncPointer f) { fFunc.f = f; } + void SetDim ( unsigned int n ) { fFunc.n = n; } + void SetParams ( void * p) { fFunc.params = p; } + + /// Fill gsl function structure from a C++ Function class + template<class FuncType> + void SetFunction(const FuncType &f) { + const void * p = &f; + assert (p != 0); + SetFuncPointer(&GSLMultiMinFunctionAdapter<FuncType >::F); + SetDim( f.NDim() ); + SetParams(const_cast<void *>(p)); + } + + gsl_multimin_function * GetFunc() { return &fFunc; } + + bool IsValid() { + return (fFunc.f != 0) ? true : false; + } + + + private: + + gsl_multimin_function fFunc; + + }; + + +/** + Wrapper for a multi-dimensional function with derivatives used in GSL multidim + minimization algorithm + + @ingroup MultiMin + + */ + class GSLMultiMinDerivFunctionWrapper { + + public: + + GSLMultiMinDerivFunctionWrapper() + { + fFunc.f = 0; + fFunc.df = 0; + fFunc.fdf = 0; + fFunc.n = 0; + fFunc.params = 0; + } + + + void SetFuncPointer( GSLMultiMinFuncPointer f) { fFunc.f = f; } + void SetDerivPointer( GSLMultiMinDfPointer f) { fFunc.df = f; } + void SetFdfPointer( GSLMultiMinFdfPointer f) { fFunc.fdf = f; } + void SetDim ( unsigned int n ) { fFunc.n = n; } + void SetParams ( void * p) { fFunc.params = p; } + + /// Fill gsl function structure from a C++ Function class + template<class FuncType> + void SetFunction(const FuncType &f) { + const void * p = &f; + assert (p != 0); + SetFuncPointer(&GSLMultiMinFunctionAdapter<FuncType >::F); + SetDerivPointer(&GSLMultiMinFunctionAdapter<FuncType >::Df); + SetFdfPointer(&GSLMultiMinFunctionAdapter<FuncType >::Fdf); + SetDim( f.NDim() ); + SetParams(const_cast<void *>(p)); + } + + gsl_multimin_function_fdf * GetFunc() { return &fFunc; } + +#ifdef NEEDED_LATER + // evaluate the function + double operator() (const double * x) { + // vx must be a gsl_vector + return GSL_MULTIMIN_FN_EVAL(&fFunc, vx); + } +#endif + + /// check if function is valid (has been set) + bool IsValid() { + return (fFunc.f != 0) ? true : false; + } + + private: + + gsl_multimin_function_fdf fFunc; + + }; + + + + +} // namespace Math +} // namespace ROOT + +#endif /* ROOT_Math_GSLMultiMinFunctionWrapper */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMultiMinimizer.h b/ThirdParty/RootMinimizers/src/Math/GSLMultiMinimizer.h new file mode 100644 index 00000000000..e8e83bcf9ed --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLMultiMinimizer.h @@ -0,0 +1,213 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Tue Dec 19 14:09:15 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLMultiMinimizer + +#ifndef ROOT_Math_GSLMultiMinimizer +#define ROOT_Math_GSLMultiMinimizer + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_multimin.h" +#include "GSLMultiMinFunctionWrapper.h" + +#include "Math/Error.h" + +#include "Math/IFunction.h" + +#include <cassert> + +namespace ROOT { + + namespace Math { + + +/** + GSLMultiMinimizer class , for minimizing multi-dimensional function + using derivatives + + @ingroup MultiMin + +*/ +class GSLMultiMinimizer { + +public: + + /** + Default constructor + */ + GSLMultiMinimizer (ROOT::Math::EGSLMinimizerType type) : + fMinimizer(0), + fType(0), + fVec(0) + { + switch(type) + { + case ROOT::Math::kConjugateFR : + fType = gsl_multimin_fdfminimizer_conjugate_fr; + break; + case ROOT::Math::kConjugatePR : + fType = gsl_multimin_fdfminimizer_conjugate_pr; + break; + case ROOT::Math::kVectorBFGS : + fType = gsl_multimin_fdfminimizer_vector_bfgs; + break; + case ROOT::Math::kVectorBFGS2 : +#if defined (GSL_VERSION_NUM) && GSL_VERSION_NUM >= 1009 + // bfgs2 is available only for v>= 1.9 + fType = gsl_multimin_fdfminimizer_vector_bfgs2; +#else + MATH_INFO_MSG("GSLMultiMinimizer","minimizer BFSG2 does not exist with this GSL version , use BFGS"); + fType = gsl_multimin_fdfminimizer_vector_bfgs; +#endif + break; + case ROOT::Math::kSteepestDescent: + fType = gsl_multimin_fdfminimizer_steepest_descent; + break; + default: + fType = gsl_multimin_fdfminimizer_conjugate_fr; + break; + } + + } + + /** + Destructor + */ + ~GSLMultiMinimizer () { + if (fMinimizer != 0 ) gsl_multimin_fdfminimizer_free(fMinimizer); + // can free vector (is copied inside) + if (fVec != 0) gsl_vector_free(fVec); + } + +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + GSLMultiMinimizer(const GSLMultiMinimizer &) {} + + /** + Assignment operator + */ + GSLMultiMinimizer & operator = (const GSLMultiMinimizer & rhs) { + if (this == &rhs) return *this; // time saving self-test + return *this; + } + +public: + + /** + set the function to be minimize the initial minimizer parameters, + step size and tolerance in the line search + */ + int Set(const ROOT::Math::IMultiGradFunction & func, const double * x, double stepSize, double tol) { + // create function wrapper + fFunc.SetFunction(func); + // create minimizer object (free previous one if already existing) + unsigned int ndim = func.NDim(); + CreateMinimizer( ndim ); + // set initial values + if (fVec != 0) gsl_vector_free(fVec); + fVec = gsl_vector_alloc( ndim ); + std::copy(x,x+ndim, fVec->data); + assert(fMinimizer != 0); + return gsl_multimin_fdfminimizer_set(fMinimizer, fFunc.GetFunc(), fVec, stepSize, tol); + } + + /// create the minimizer from the type and size + void CreateMinimizer(unsigned int n) { + if (fMinimizer) gsl_multimin_fdfminimizer_free(fMinimizer); + fMinimizer = gsl_multimin_fdfminimizer_alloc(fType, n); + } + + std::string Name() const { + if (fMinimizer == 0) return "undefined"; + return std::string(gsl_multimin_fdfminimizer_name(fMinimizer) ); + } + + int Iterate() { + if (fMinimizer == 0) return -1; + return gsl_multimin_fdfminimizer_iterate(fMinimizer); + } + + /// x values at the minimum + double * X() const { + if (fMinimizer == 0) return 0; + gsl_vector * x = gsl_multimin_fdfminimizer_x(fMinimizer); + return x->data; + } + + /// function value at the minimum + double Minimum() const { + if (fMinimizer == 0) return 0; + return gsl_multimin_fdfminimizer_minimum(fMinimizer); + } + + /// gradient value at the minimum + double * Gradient() const { + if (fMinimizer == 0) return 0; + gsl_vector * g = gsl_multimin_fdfminimizer_gradient(fMinimizer); + return g->data; + } + + /// restart minimization from current point + int Restart() { + if (fMinimizer == 0) return -1; + return gsl_multimin_fdfminimizer_restart(fMinimizer); + } + + /// test gradient (ask from minimizer gradient vector) + int TestGradient(double absTol) const { + if (fMinimizer == 0) return -1; + gsl_vector * g = gsl_multimin_fdfminimizer_gradient(fMinimizer); + return gsl_multimin_test_gradient( g, absTol); + } + + /// test gradient (require a vector gradient) + int TestGradient(const double * g, double absTol) const { + if (fVec == 0 ) return -1; + unsigned int n = fVec->size; + if (n == 0 ) return -1; + std::copy(g,g+n, fVec->data); + return gsl_multimin_test_gradient( fVec, absTol); + } + + +private: + + gsl_multimin_fdfminimizer * fMinimizer; + GSLMultiMinDerivFunctionWrapper fFunc; + const gsl_multimin_fdfminimizer_type * fType; + // cached vector to avoid re-allocating every time a new one + mutable gsl_vector * fVec; + +}; + + } // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_GSLMultiMinimizer */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLNLSMinimizer.cxx b/ThirdParty/RootMinimizers/src/Math/GSLNLSMinimizer.cxx new file mode 100644 index 00000000000..f22b9db4207 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLNLSMinimizer.cxx @@ -0,0 +1,463 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Dec 20 17:16:32 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * * + **********************************************************************/ + +// Implementation file for class GSLNLSMinimizer + +#include "Math/GSLNLSMinimizer.h" + +#include "Math/MinimTransformFunction.h" +#include "Math/MultiNumGradFunction.h" + +#include "Math/Error.h" +#include "GSLMultiFit.h" +#include "gsl/gsl_errno.h" + + +#include "Math/FitMethodFunction.h" +//#include "Math/Derivator.h" + +#include <iostream> +#include <iomanip> +#include <cassert> +#include <memory> + +namespace ROOT { + + namespace Math { + + +// class to implement transformation of chi2 function +// in general could make template on the fit method function type + +class FitTransformFunction : public FitMethodFunction { + +public: + + FitTransformFunction(const FitMethodFunction & f, const std::vector<EMinimVariableType> & types, const std::vector<double> & values, + const std::map<unsigned int, std::pair<double, double> > & bounds) : + FitMethodFunction( f.NDim(), f.NPoints() ), + fFunc(f), + fTransform(new MinimTransformFunction( new MultiNumGradFunction(f), types, values, bounds) ), + fGrad( std::vector<double>(f.NDim() ) ) + { + // constructor + // need to pass to MinimTransformFunction a new pointer which will be managed by the class itself + // pass a gradient pointer although it will not be used byb the class + } + + ~FitTransformFunction() { + assert(fTransform); + delete fTransform; + } + + // re-implement data element + virtual double DataElement(const double * x, unsigned i, double * g = 0) const { + // transform from x internal to x external + const double * xExt = fTransform->Transformation(x); + if ( g == 0) return fFunc.DataElement( xExt, i ); + // use gradient + double val = fFunc.DataElement( xExt, i, &fGrad[0]); + // transform gradient + fTransform->GradientTransformation( x, &fGrad.front(), g); + return val; + } + + + IMultiGenFunction * Clone() const { + // not supported + return 0; + } + + // dimension (this is number of free dimensions) + unsigned int NDim() const { + return fTransform->NDim(); + } + + unsigned int NTot() const { + return fTransform->NTot(); + } + + // forward of transformation functions + const double * Transformation( const double * x) const { return fTransform->Transformation(x); } + + + /// inverse transformation (external -> internal) + void InvTransformation(const double * xext, double * xint) const { fTransform->InvTransformation(xext,xint); } + + /// inverse transformation for steps (external -> internal) at external point x + void InvStepTransformation(const double * x, const double * sext, double * sint) const { fTransform->InvStepTransformation(x,sext,sint); } + + ///transform gradient vector (external -> internal) at internal point x + void GradientTransformation(const double * x, const double *gext, double * gint) const { fTransform->GradientTransformation(x,gext,gint); } + + void MatrixTransformation(const double * x, const double *cint, double * cext) const { fTransform->MatrixTransformation(x,cint,cext); } + +private: + + double DoEval(const double * x) const { + return fFunc( fTransform->Transformation(x) ); + } + + const FitMethodFunction & fFunc; // pointer to original fit method function + MinimTransformFunction * fTransform; // pointer to transformation function + mutable std::vector<double> fGrad; // cached vector of gradient values + +}; + + + + +// GSLNLSMinimizer implementation + +GSLNLSMinimizer::GSLNLSMinimizer( int type ) : + fDim(0), + fNFree(0), + fSize(0), + fObjFunc(0), + fMinVal(0) +{ + // Constructor implementation : create GSLMultiFit wrapper object + const gsl_multifit_fdfsolver_type * gsl_type = 0; // use default type defined in GSLMultiFit + if (type == 1) gsl_type = gsl_multifit_fdfsolver_lmsder; // scaled lmder version + if (type == 2) gsl_type = gsl_multifit_fdfsolver_lmder; // unscaled version + + fGSLMultiFit = new GSLMultiFit( gsl_type ); + fValues.reserve(10); + fNames.reserve(10); + fSteps.reserve(10); + + fEdm = -1; + + // defult tolerance and max iterations + int niter = ROOT::Math::MinimizerOptions::DefaultMaxIterations(); + if (niter <= 0) niter = 100; + SetMaxIterations(niter); + + fLSTolerance = ROOT::Math::MinimizerOptions::DefaultTolerance(); + if (fLSTolerance <=0) fLSTolerance = 0.0001; // default internal value + + SetPrintLevel(ROOT::Math::MinimizerOptions::DefaultPrintLevel()); +} + +GSLNLSMinimizer::~GSLNLSMinimizer () { + assert(fGSLMultiFit != 0); + delete fGSLMultiFit; +} + +bool GSLNLSMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) { + // set variable in minimizer - support only free variables + // no transformation implemented - so far + if (ivar > fValues.size() ) return false; + if (ivar == fValues.size() ) { + fValues.push_back(val); + fNames.push_back(name); + fSteps.push_back(step); + fVarTypes.push_back(kDefault); + } + else { + fValues[ivar] = val; + fNames[ivar] = name; + fSteps[ivar] = step; + fVarTypes[ivar] = kDefault; + + // remove bounds if needed + std::map<unsigned int, std::pair<double, double> >::iterator iter = fBounds.find(ivar); + if ( iter != fBounds.end() ) fBounds.erase (iter); + + } + return true; +} + +bool GSLNLSMinimizer::SetLowerLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower) { + //MATH_WARN_MSGVAL("GSLNLSMinimizer::SetLowerLimitedVariable","Ignore lower limit on variable ",ivar); + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( lower, lower); + fVarTypes[ivar] = kLowBound; + return true; +} +bool GSLNLSMinimizer::SetUpperLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double upper) { + //MATH_WARN_MSGVAL("GSLNLSMinimizer::SetUpperLimitedVariable","Ignore upper limit on variable ",ivar); + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( upper, upper); + fVarTypes[ivar] = kUpBound; + return true; +} +bool GSLNLSMinimizer::SetLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower, double upper ) { + //MATH_WARN_MSGVAL("GSLNLSMinimizer::SetLimitedVariable","Ignore bounds on variable ",ivar); + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( lower, upper); + fVarTypes[ivar] = kBounds; + return true; +} + +bool GSLNLSMinimizer::SetFixedVariable(unsigned int ivar , const std::string & name , double val ) { + /// set fixed variable (override if minimizer supports them ) + bool ret = SetVariable(ivar, name, val, 0.); + if (!ret) return false; + fVarTypes[ivar] = kFix; + return true; +} + +bool GSLNLSMinimizer::SetVariableValue(unsigned int ivar, double val) { + // set variable value in minimizer + // no transformation implemented - so far + if (ivar > fValues.size() ) return false; + fValues[ivar] = val; + return true; +} + +bool GSLNLSMinimizer::SetVariableValues( const double * x) { + // set all variable values in minimizer + if (x == 0) return false; + std::copy(x,x+fValues.size(), fValues.begin() ); + return true; +} + + + +void GSLNLSMinimizer::SetFunction(const ROOT::Math::IMultiGenFunction & func) { + // set the function to minimizer + // need to create vector of functions to be passed to GSL multifit + // support now only CHi2 implementation + + const ROOT::Math::FitMethodFunction * chi2Func = dynamic_cast<const ROOT::Math::FitMethodFunction *>(&func); + if (chi2Func == 0) { + if (PrintLevel() > 0) std::cout << "GSLNLSMinimizer: Invalid function set - only Chi2Func supported" << std::endl; + return; + } + fSize = chi2Func->NPoints(); + fDim = chi2Func->NDim(); + fNFree = fDim; + + // use vector by value + fResiduals.reserve(fSize); + for (unsigned int i = 0; i < fSize; ++i) { + fResiduals.push_back( LSResidualFunc(*chi2Func, i) ); + } + // keep pointers to the chi2 function + fObjFunc = chi2Func; + } + +void GSLNLSMinimizer::SetFunction(const ROOT::Math::IMultiGradFunction & func ) { + // set the function to minimizer using gradient interface + // not supported yet, implemented using the other SetFunction + return SetFunction(static_cast<const ROOT::Math::IMultiGenFunction &>(func) ); +} + + +bool GSLNLSMinimizer::Minimize() { + // set initial parameters of the minimizer + int debugLevel = PrintLevel(); + + + assert (fGSLMultiFit != 0); + if (fResiduals.size() != fSize || fObjFunc == 0) { + MATH_ERROR_MSG("GSLNLSMinimizer::Minimize","Function has not been set"); + return false; + } + + unsigned int npar = fValues.size(); + if (npar == 0 || npar < fDim) { + MATH_ERROR_MSGVAL("GSLNLSMinimizer::Minimize","Wrong number of parameters",npar); + return false; + } + + // set residual functions and check if a transformation is needed + + bool doTransform = (fBounds.size() > 0); + unsigned int ivar = 0; + while (!doTransform && ivar < fVarTypes.size() ) { + doTransform = (fVarTypes[ivar++] != kDefault ); + } + std::vector<double> startValues(fValues.begin(), fValues.end() ); + + std::auto_ptr<FitTransformFunction> trFunc; + + // in case of transformation wrap residual functions in new transformation functions + // and transform from external variables to internals one + if (doTransform) { + trFunc.reset(new FitTransformFunction ( *fObjFunc, fVarTypes, fValues, fBounds ) ); + for (unsigned int ires = 0; ires < fResiduals.size(); ++ires) { + fResiduals[ires] = LSResidualFunc(*trFunc, ires); + } + + trFunc->InvTransformation(&fValues.front(), &startValues[0]); + fNFree = trFunc->NDim(); // actual dimension + assert(fValues.size() == trFunc->NTot() ); + startValues.resize( fNFree ); + } + + if (debugLevel >=1 ) std::cout <<"Minimize using GSLNLSMinimizer " << std::endl; + +// // use a global step size = min (step vectors) +// double stepSize = 1; +// for (unsigned int i = 0; i < fSteps.size(); ++i) +// //stepSize += fSteps[i]; +// if (fSteps[i] < stepSize) stepSize = fSteps[i]; + + int iret = fGSLMultiFit->Set( fResiduals, &startValues.front() ); + if (iret) { + MATH_ERROR_MSGVAL("GSLNLSMinimizer::Minimize","Error setting the residual functions ",iret); + return false; + } + + if (debugLevel >=1 ) std::cout <<"GSLNLSMinimizer: " << fGSLMultiFit->Name() << " - start iterating......... " << std::endl; + + // start iteration + unsigned int iter = 0; + int status; + bool minFound = false; + do { + status = fGSLMultiFit->Iterate(); + + if (debugLevel >=1) { + std::cout << "----------> Iteration " << iter << " / " << MaxIterations() << " status " << gsl_strerror(status) << std::endl; + const double * x = fGSLMultiFit->X(); + if (trFunc.get()) x = trFunc->Transformation(x); + int pr = std::cout.precision(18); + std::cout << " FVAL = " << (*fObjFunc)(x) << std::endl; + std::cout.precision(pr); + std::cout << " X Values : "; + for (unsigned int i = 0; i < fDim; ++i) + std::cout << " " << fNames[i] << " = " << x[i]; + std::cout << std::endl; + } + + if (status) break; + + // check also the delta in X() + status = fGSLMultiFit->TestDelta( Tolerance(), Tolerance() ); + if (status == GSL_SUCCESS) { + minFound = true; + } + + // double-check with the gradient + int status2 = fGSLMultiFit->TestGradient( Tolerance() ); + if ( minFound && status2 != GSL_SUCCESS) { + // check now edm + fEdm = fGSLMultiFit->Edm(); + if (fEdm > Tolerance() ) { + // continue the iteration + status = status2; + minFound = false; + } + } + + if (debugLevel >=1) { + std::cout << " after Gradient and Delta tests: " << gsl_strerror(status); + if (fEdm > 0) std::cout << ", edm is: " << fEdm; + std::cout << std::endl; + } + + iter++; + + } + while (status == GSL_CONTINUE && iter < MaxIterations() ); + + // check edm + fEdm = fGSLMultiFit->Edm(); + if ( fEdm < Tolerance() ) { + minFound = true; + } + + // save state with values and function value + const double * x = fGSLMultiFit->X(); + if (x == 0) return false; + + // check to see if a transformation need to be applied + if (trFunc.get() != 0) { + const double * xtrans = trFunc->Transformation(x); + std::copy(xtrans, xtrans + trFunc->NTot(), fValues.begin() ); + } + else { + std::copy(x, x +fDim, fValues.begin() ); + } + + fMinVal = (*fObjFunc)(&fValues.front() ); + fStatus = status; + + fErrors.resize(fDim); + + // get errors from cov matrix + const double * cov = fGSLMultiFit->CovarMatrix(); + if (cov) { + + unsigned int ndim = fDim; + fCovMatrix.resize(ndim*ndim); + + if (trFunc.get() ) { + trFunc->MatrixTransformation(x, fGSLMultiFit->CovarMatrix(), &fCovMatrix[0] ); + } + else { + std::copy(cov, cov + ndim*ndim, fCovMatrix.begin() ); + } + + for (unsigned int i = 0; i < fDim; ++i) + fErrors[i] = std::sqrt(fCovMatrix[i*fDim + i]); + } + + if (minFound) { + + if (debugLevel >=1 ) { + std::cout << "GSLNLSMinimizer: Minimum Found" << std::endl; + int pr = std::cout.precision(18); + std::cout << "FVAL = " << fMinVal << std::endl; + std::cout << "Edm = " << fEdm << std::endl; + std::cout.precision(pr); + std::cout << "NIterations = " << iter << std::endl; + std::cout << "NFuncCalls = " << fObjFunc->NCalls() << std::endl; + for (unsigned int i = 0; i < fDim; ++i) + std::cout << std::setw(12) << fNames[i] << " = " << std::setw(12) << fValues[i] << " +/- " << std::setw(12) << fErrors[i] << std::endl; + } + + return true; + } + else { + if (debugLevel >=1 ) { + std::cout << "GSLNLSMinimizer: Minimization did not converge" << std::endl; + std::cout << "FVAL = " << fMinVal << std::endl; + std::cout << "Edm = " << fGSLMultiFit->Edm() << std::endl; + std::cout << "Niterations = " << iter << std::endl; + } + return false; + } + return false; +} + +const double * GSLNLSMinimizer::MinGradient() const { + // return gradient (internal values) + return fGSLMultiFit->Gradient(); +} + + +double GSLNLSMinimizer::CovMatrix(unsigned int i , unsigned int j ) const { + // return covariance matrix element + if ( fCovMatrix.size() == 0) return 0; + if (i > fDim || j > fDim) return 0; + return fCovMatrix[i*fDim + j]; +} + +int GSLNLSMinimizer::CovMatrixStatus( ) const { + // return covariance matrix status = 0 not computed, + // 1 computed but is approximate because minimum is not valid, 3 is fine + if ( fCovMatrix.size() == 0) return 0; + // case minimization did not finished correctly + if (fStatus != GSL_SUCCESS) return 1; + return 3; +} + + + } // end namespace Math + +} // end namespace ROOT + diff --git a/ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx b/ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx new file mode 100644 index 00000000000..89c92e9978a --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx @@ -0,0 +1,380 @@ +// @(#)root/mathmore:$Id$ +// Authors: L. Moneta, A. Zsenei 08/2005 + + /********************************************************************** + * * + * Copyright (c) 2004 ROOT Foundation, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// Header file for class GSLRandom +// +// Created by: moneta at Sun Nov 21 16:26:03 2004 +// +// Last update: Sun Nov 21 16:26:03 2004 +// + + + +// need to be included later +#include <time.h> +#include <stdlib.h> +#include <cassert> + +#include "gsl/gsl_rng.h" +#include "gsl/gsl_randist.h" + + +#include "Math/GSLRndmEngines.h" +#include "GSLRngWrapper.h" + +extern double gsl_ran_gaussian_acr( const gsl_rng * r, const double sigma); + +//#include <iostream> + +namespace ROOT { +namespace Math { + + + + + + // default constructor (need to call set type later) + GSLRandomEngine::GSLRandomEngine() : + fRng(0 ), + fCurTime(0) + { } + + // constructor from external rng + // internal generator will be managed or not depending on + // how the GSLRngWrapper is created + GSLRandomEngine::GSLRandomEngine( GSLRngWrapper * rng) : + fRng(new GSLRngWrapper(*rng) ), + fCurTime(0) + {} + +// // constructor from external rng +// GSLRandomEngine( GSLRngWrapper & rng) : +// fRng(new GSLRngWrapper(rng) ), +// fCurTime(0) +// {} + + GSLRandomEngine::~GSLRandomEngine() { + // destructor : call terminate if not yet called + if (fRng) Terminate(); + } + + + void GSLRandomEngine::Initialize() { + // initialize the generator by allocating the GSL object + // if type was not passed create with default generator + if (!fRng) fRng = new GSLRngWrapper(); + fRng->Allocate(); + } + + void GSLRandomEngine::Terminate() { + // terminate the generator by freeing the GSL object + if (!fRng) return; + fRng->Free(); + delete fRng; + fRng = 0; + } + + + double GSLRandomEngine::operator() () const { + // generate random between 0 and 1. + // 0 is excluded + return gsl_rng_uniform_pos(fRng->Rng() ); + } + + + unsigned int GSLRandomEngine::RndmInt(unsigned int max) const { + // generate a random integer number between 0 and MAX + return gsl_rng_uniform_int( fRng->Rng(), max ); + } + +// int GSLRandomEngine::GetMin() { +// // return minimum integer value used in RndmInt +// return gsl_rng_min( fRng->Rng() ); +// } + +// int GSLRandomEngine::GetMax() { +// // return maximum integr value used in RndmInt +// return gsl_rng_max( fRng->Rng() ); +// } + + void GSLRandomEngine::RandomArray(double * begin, double * end ) const { + // generate array of randoms betweeen 0 and 1. 0 is excluded + // specialization for double * (to be faster) + for ( double * itr = begin; itr != end; ++itr ) { + *itr = gsl_rng_uniform_pos(fRng->Rng() ); + } + } + + void GSLRandomEngine::SetSeed(unsigned int seed) const { + // set the seed, if = 0then the seed is set randomly using an std::rand() + // seeded with the current time. Be carefuk in case the current time is + // the same in consecutive calls + if (seed == 0) { + // use like in root (use time) + time_t curtime; + time(&curtime); + unsigned int ct = static_cast<unsigned int>(curtime); + if (ct != fCurTime) { + fCurTime = ct; + // set the seed for rand + srand(ct); + } + seed = rand(); + } + + assert(fRng); + gsl_rng_set(fRng->Rng(), seed ); + } + + std::string GSLRandomEngine::Name() const { + //---------------------------------------------------- + assert ( fRng != 0); + assert ( fRng->Rng() != 0 ); + return std::string( gsl_rng_name( fRng->Rng() ) ); + } + + unsigned int GSLRandomEngine::Size() const { + //---------------------------------------------------- + assert (fRng != 0); + return gsl_rng_size( fRng->Rng() ); + } + + + // Random distributions + + double GSLRandomEngine::GaussianZig(double sigma) const + { + // Gaussian distribution. Use fast ziggurat algorithm implemented since GSL 1.8 + return gsl_ran_gaussian_ziggurat( fRng->Rng(), sigma); + } + + double GSLRandomEngine::Gaussian(double sigma) const + { + // Gaussian distribution. Use default Box-Muller method + return gsl_ran_gaussian( fRng->Rng(), sigma); + } + + double GSLRandomEngine::GaussianRatio(double sigma) const + { + // Gaussian distribution. Use ratio method + return gsl_ran_gaussian_ratio_method( fRng->Rng(), sigma); + } + + + double GSLRandomEngine::GaussianTail(double a , double sigma) const + { + // Gaussian Tail distribution: eeturn values larger than a distributed + // according to the gaussian + return gsl_ran_gaussian_tail( fRng->Rng(), a, sigma); + } + + void GSLRandomEngine::Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) const + { + // Gaussian Bivariate distribution, with correlation coefficient rho + gsl_ran_bivariate_gaussian( fRng->Rng(), sigmaX, sigmaY, rho, &x, &y); + } + + double GSLRandomEngine::Exponential(double mu) const + { + // Exponential distribution + return gsl_ran_exponential( fRng->Rng(), mu); + } + + double GSLRandomEngine::Cauchy(double a) const + { + // Cauchy distribution + return gsl_ran_cauchy( fRng->Rng(), a); + } + + double GSLRandomEngine::Landau() const + { + // Landau distribution + return gsl_ran_landau( fRng->Rng()); + } + + double GSLRandomEngine::Gamma(double a, double b) const + { + // Gamma distribution + return gsl_ran_gamma( fRng->Rng(), a, b); + } + + double GSLRandomEngine::LogNormal(double zeta, double sigma) const + { + // Log normal distribution + return gsl_ran_lognormal( fRng->Rng(), zeta, sigma); + } + + double GSLRandomEngine::ChiSquare(double nu) const + { + // Chi square distribution + return gsl_ran_chisq( fRng->Rng(), nu); + } + + + double GSLRandomEngine::FDist(double nu1, double nu2) const + { + // F distribution + return gsl_ran_fdist( fRng->Rng(), nu1, nu2); + } + + double GSLRandomEngine::tDist(double nu) const + { + // t distribution + return gsl_ran_tdist( fRng->Rng(), nu); + } + + void GSLRandomEngine::Dir2D(double &x, double &y) const + { + // generate random numbers in a 2D circle of radious 1 + gsl_ran_dir_2d( fRng->Rng(), &x, &y); + } + + void GSLRandomEngine::Dir3D(double &x, double &y, double &z) const + { + // generate random numbers in a 3D sphere of radious 1 + gsl_ran_dir_3d( fRng->Rng(), &x, &y, &z); + } + + unsigned int GSLRandomEngine::Poisson(double mu) const + { + // Poisson distribution + return gsl_ran_poisson( fRng->Rng(), mu); + } + + unsigned int GSLRandomEngine::Binomial(double p, unsigned int n) const + { + // Binomial distribution + return gsl_ran_binomial( fRng->Rng(), p, n); + } + + unsigned int GSLRandomEngine::NegativeBinomial(double p, double n) const + { + // Negative Binomial distribution + return gsl_ran_negative_binomial( fRng->Rng(), p, n); + } + + + std::vector<unsigned int> GSLRandomEngine::Multinomial( unsigned int ntot, const std::vector<double> & p ) const + { + // Multinomial distribution return vector of integers which sum is ntot + std::vector<unsigned int> ival( p.size()); + gsl_ran_multinomial( fRng->Rng(), p.size(), ntot, &p.front(), &ival[0]); + return ival; + } + + + + //---------------------------------------------------- + // generators + //---------------------------------------------------- + + //---------------------------------------------------- + GSLRngMT::GSLRngMT() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_mt19937)); + } + + + // old ranlux - equivalent to TRandom1 + GSLRngRanLux::GSLRngRanLux() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_ranlux) ); + } + + // second generation of Ranlux (single precision version - luxury 1) + GSLRngRanLuxS1::GSLRngRanLuxS1() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_ranlxs1) ); + } + + // second generation of Ranlux (single precision version - luxury 2) + GSLRngRanLuxS2::GSLRngRanLuxS2() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_ranlxs2) ); + } + + // double precision version - luxury 1 + GSLRngRanLuxD1::GSLRngRanLuxD1() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_ranlxd1) ); + } + + // double precision version - luxury 2 + GSLRngRanLuxD2::GSLRngRanLuxD2() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_ranlxd2) ); + } + + //---------------------------------------------------- + GSLRngTaus::GSLRngTaus() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_taus2) ); + } + + //---------------------------------------------------- + GSLRngGFSR4::GSLRngGFSR4() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_gfsr4) ); + } + + //---------------------------------------------------- + GSLRngCMRG::GSLRngCMRG() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_cmrg) ); + } + + //---------------------------------------------------- + GSLRngMRG::GSLRngMRG() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_mrg) ); + } + + + //---------------------------------------------------- + GSLRngRand::GSLRngRand() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_rand) ); + } + + //---------------------------------------------------- + GSLRngRanMar::GSLRngRanMar() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_ranmar) ); + } + + //---------------------------------------------------- + GSLRngMinStd::GSLRngMinStd() : GSLRandomEngine() + { + SetType(new GSLRngWrapper(gsl_rng_minstd) ); + } + + + + + +} // namespace Math +} // namespace ROOT + + + diff --git a/ThirdParty/RootMinimizers/src/Math/GSLRngWrapper.h b/ThirdParty/RootMinimizers/src/Math/GSLRngWrapper.h new file mode 100644 index 00000000000..55a34ea197a --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLRngWrapper.h @@ -0,0 +1,140 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Fri Aug 24 17:20:45 2007 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * * + **********************************************************************/ + +// Header file for class GSLRngWrapper + +#ifndef ROOT_Math_GSLRngWrapper +#define ROOT_Math_GSLRngWrapper + + +namespace ROOT { + + namespace Math { + + +/** + GSLRngWrapper class to wrap gsl_rng structure +*/ +class GSLRngWrapper { + +public: + + + /** + Default constructor + */ + GSLRngWrapper () : + fOwn(0), + fRng(0), + fRngType(0) + { + } + + /** + Constructor with type + */ + GSLRngWrapper(const gsl_rng_type * type) : + fOwn(1), + fRng(0), + fRngType(type) + { + } + + /** + construct from an existing gsl_rng + it is managed externally - so will not be deleted at the end + */ + GSLRngWrapper(const gsl_rng * r ) : + fOwn(0), + fRngType(0) + { + fRng = const_cast<gsl_rng *>(r); + } + + /** + Copy constructor - pass ownership (need not to be const) + Just copy the pointer and do not manage it + */ + GSLRngWrapper(GSLRngWrapper & r) : + fOwn(r.fOwn), + fRng(r.fRng), + fRngType(r.fRngType) + { + // in case an rng exists must release it + if (fRng && fOwn) r.fOwn = false; + } + + + /** + Destructor (free the rng if not done before) + */ + ~GSLRngWrapper() { + if (fOwn) Free(); + } + + void Allocate() { + if (fRngType == 0) SetDefaultType(); + if (fRng != 0 && fOwn) Free(); + fRng = gsl_rng_alloc( fRngType ); + //std::cout << " allocate " << fRng << std::endl; + } + + void Free() { + if (!fOwn) return; // no operation if pointer is not own + //std::cout << "free gslrng " << fRngType << " " << fRng << std::endl; + if (fRng != 0) gsl_rng_free(fRng); + fRng = 0; + } + + + void SetType(const gsl_rng_type * type) { + fRngType = type; + } + + void SetDefaultType() { + // construct default engine + gsl_rng_env_setup(); + fRngType = gsl_rng_default; + } + + + + inline gsl_rng * Rng() { return fRng; } + + inline const gsl_rng * Rng() const { return fRng; } + +private: + // usually copying is non trivial, so we make this unaccessible + + + /** + Assignment operator + Disable since if don't want to change an already created wrapper + */ + GSLRngWrapper & operator = (const GSLRngWrapper & rhs) { + if (this == &rhs) return *this; // time saving self-test + return *this; + } + + +private: + + bool fOwn; // ownership of contained pointer + gsl_rng * fRng; + const gsl_rng_type * fRngType; +}; + + + } // end namespace Math + +} // end namespace ROOT + + +#endif /* ROOT_Math_GSLRngWrapper */ diff --git a/ThirdParty/RootMinimizers/src/Math/GSLSimAnMinimizer.cxx b/ThirdParty/RootMinimizers/src/Math/GSLSimAnMinimizer.cxx new file mode 100644 index 00000000000..39180d75098 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLSimAnMinimizer.cxx @@ -0,0 +1,230 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Dec 20 17:16:32 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * * + **********************************************************************/ + +// Implementation file for class GSLSimAnMinimizer + +#include "Math/GSLSimAnMinimizer.h" +#include "Math/WrappedParamFunction.h" +#include "Math/Error.h" + +#include "Math/MinimTransformFunction.h" +#include "Math/MultiNumGradFunction.h" // needed to use transformation function + +#include <iostream> +#include <cassert> + +namespace ROOT { + + namespace Math { + + + +// GSLSimAnMinimizer implementation + +GSLSimAnMinimizer::GSLSimAnMinimizer( int /* ROOT::Math::EGSLSimAnMinimizerType type */ ) : + fDim(0), + fOwnFunc(false), + fObjFunc(0), + fMinVal(0) +{ + // Constructor implementation : create GSLMultiFit wrapper object + + fValues.reserve(10); + fNames.reserve(10); + fSteps.reserve(10); + + SetMaxIterations(100); + SetPrintLevel(0); +} + +GSLSimAnMinimizer::~GSLSimAnMinimizer () { + if ( fOwnFunc && fObjFunc) delete fObjFunc; +} + +bool GSLSimAnMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) { + // set variable in minimizer - support only free variables + // no transformation implemented - so far + if (ivar > fValues.size() ) return false; + if (ivar == fValues.size() ) { + fValues.push_back(val); + fNames.push_back(name); + // step is the simmulated annealing scale + fSteps.push_back( step ); + fVarTypes.push_back(kDefault); + } + else { + fValues[ivar] = val; + fNames[ivar] = name; + fSteps[ivar] = step; + fVarTypes[ivar] = kDefault; + + // remove bounds if needed + std::map<unsigned int, std::pair<double, double> >::iterator iter = fBounds.find(ivar); + if ( iter != fBounds.end() ) fBounds.erase (iter); + } + return true; + +} + +bool GSLSimAnMinimizer::SetLowerLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower ) { + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( lower, lower); + fVarTypes[ivar] = kLowBound; + return true; +} +bool GSLSimAnMinimizer::SetUpperLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double upper) { + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( upper, upper); + fVarTypes[ivar] = kUpBound; + return true; +} +bool GSLSimAnMinimizer::SetLimitedVariable(unsigned int ivar, const std::string & name, double val, double step, double lower, double upper ) { + bool ret = SetVariable(ivar, name, val, step); + if (!ret) return false; + fBounds[ivar] = std::make_pair( lower, upper); + fVarTypes[ivar] = kBounds; + return true; +} + + + +bool GSLSimAnMinimizer::SetFixedVariable(unsigned int ivar, const std::string & name, double val) { + /// set fixed variable (override if minimizer supports them ) + // use zero step size + bool ret = SetVariable(ivar, name, val, 0.); + if (!ret) return false; + fVarTypes[ivar] = kFix; + return true; +} + +bool GSLSimAnMinimizer::SetVariableValue(unsigned int ivar, double val) { + // set variable value in minimizer + // no transformation implemented - so far + if (ivar > fValues.size() ) return false; + fValues[ivar] = val; + return true; +} + +bool GSLSimAnMinimizer::SetVariableValues( const double * x) { + // set all variable values in minimizer + if (x == 0) return false; + std::copy(x,x+fValues.size(), fValues.begin() ); + return true; +} + + +void GSLSimAnMinimizer::SetFunction(const ROOT::Math::IMultiGenFunction & func) { + // set the function to minimize + + // keep pointers to the chi2 function + fObjFunc = &func; + fDim = func.NDim(); +} + +void GSLSimAnMinimizer::SetFunction(const ROOT::Math::IMultiGradFunction & func ) { + // set the function to minimize + // use the other methods + SetFunction( static_cast<const ROOT::Math::IMultiGenFunction &>(func) ); +} + + +bool GSLSimAnMinimizer::Minimize() { + // set initial parameters of the minimizer + int debugLevel = PrintLevel(); + + if (debugLevel >=1 ) std::cout <<"Minimize using GSLSimAnMinimizer " << std::endl; + + + // adapt the steps (use largers) + for (unsigned int i = 0; i < fSteps.size() ; ++i) + fSteps[i] *= 10; + + // vector of internal values (copied by default) + std::vector<double> xvar (fValues ); + std::vector<double> steps (fSteps); + + + // check if a transformation is needed + bool doTransform = (fBounds.size() > 0); + unsigned int ivar = 0; + while (!doTransform && ivar < fVarTypes.size() ) { + doTransform = (fVarTypes[ivar++] != kDefault ); + } + + // if needed do transformation and wrap objective function in a new transformation function + // and transform from external variables (and steps) to internals one + MinimTransformFunction * trFunc = 0; + if (doTransform) { + // since objective function is gradient build a gradient function for the transformation + // although gradient is not needed + + trFunc = new MinimTransformFunction ( new MultiNumGradFunction( *fObjFunc), fVarTypes, fValues, fBounds ); + + trFunc->InvTransformation(&fValues.front(), &xvar[0]); + + // need to transform also the steps + trFunc->InvStepTransformation(&fValues.front(), &fSteps.front(), &steps[0]); + + xvar.resize( trFunc->NDim() ); + steps.resize( trFunc->NDim() ); + + fObjFunc = trFunc; + fOwnFunc = true; // flag to indicate we need to delete the function + } + + assert (xvar.size() == steps.size() ); + + // output vector + std::vector<double> xmin(xvar.size() ); + + int iret = fSolver.Solve(*fObjFunc, &xvar.front(), &steps.front(), &xmin[0], (debugLevel > 1) ); + + fMinVal = (*fObjFunc)(&xmin.front() ); + + // get the result (transform if needed) + if (trFunc != 0) { + const double * xtrans = trFunc->Transformation(&xmin.front()); + assert(fValues.size() == trFunc->NTot() ); + assert( trFunc->NTot() == NDim() ); + std::copy(xtrans, xtrans + trFunc->NTot(), fValues.begin() ); + } + else { + // case of no transformation applied + assert( fValues.size() == xmin.size() ); + std::copy(xmin.begin(), xmin.end(), fValues.begin() ); + } + + + + if (debugLevel >=1 ) { + if (iret == 0) + std::cout << "GSLSimAnMinimizer: Minimum Found" << std::endl; + else + std::cout << "GSLSimAnMinimizer: Error in solving" << std::endl; + + int pr = std::cout.precision(18); + std::cout << "FVAL = " << fMinVal << std::endl; + std::cout.precision(pr); + for (unsigned int i = 0; i < fDim; ++i) + std::cout << fNames[i] << "\t = " << fValues[i] << std::endl; + } + + + return ( iret == 0) ? true : false; +} + + + + } // end namespace Math + +} // end namespace ROOT + diff --git a/ThirdParty/RootMinimizers/src/Math/GSLSimAnnealing.cxx b/ThirdParty/RootMinimizers/src/Math/GSLSimAnnealing.cxx new file mode 100644 index 00000000000..92091ac967b --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/GSLSimAnnealing.cxx @@ -0,0 +1,240 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Thu Jan 25 11:13:48 2007 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * * + **********************************************************************/ + +// Implementation file for class GSLSimAnnealing + +#include "Math/GSLSimAnnealing.h" + +#include "gsl/gsl_siman.h" + +#include "Math/IFunction.h" +#include "Math/GSLRndmEngines.h" +#include "GSLRngWrapper.h" + + +#include <cassert> +#include <iostream> +#include <cmath> +#include <vector> + +namespace ROOT { + + namespace Math { + + +// implementation of GSLSimAnFunc + +GSLSimAnFunc::GSLSimAnFunc(const ROOT::Math::IMultiGenFunction & func, const double * x) : + fX( std::vector<double>(x, x + func.NDim() ) ), + fScale( std::vector<double>(func.NDim() )), + fFunc(&func) +{ + // set scale factors to 1 + fScale.assign(fScale.size(), 1.); +} + +GSLSimAnFunc::GSLSimAnFunc(const ROOT::Math::IMultiGenFunction & func, const double * x, const double * scale) : + fX( std::vector<double>(x, x + func.NDim() ) ), + fScale( std::vector<double>(scale, scale + func.NDim() ) ), + fFunc(&func) +{} + + +double GSLSimAnFunc::Energy() const { + // evaluate the energy + return (*fFunc)(&fX.front() ); +} + +void GSLSimAnFunc::Step(const GSLRandomEngine & random, double maxstep) { + // x -> x + Random[-step,step] for each coordinate + unsigned int ndim = NDim(); + for (unsigned int i = 0; i < ndim; ++i) { + double urndm = random(); + double sstep = maxstep * fScale[i]; + fX[i] += 2 * sstep * urndm - sstep; + } +} + + +double GSLSimAnFunc::Distance(const GSLSimAnFunc & f) const { + // calculate the distance with respect onother configuration + const std::vector<double> & x = fX; + const std::vector<double> & y = f.X(); + unsigned int n = x.size(); + assert (n == y.size()); + if (n > 1) { + double d2 = 0; + for (unsigned int i = 0; i < n; ++i) + d2 += ( x[i] - y[i] ) * ( x[i] - y[i] ); + return std::sqrt(d2); + } + else + // avoid doing a sqrt for 1 dim + return std::abs( x[0] - y[0] ); +} + +void GSLSimAnFunc::Print() { + // print the position x in standard ostream + // GSL prints also niter- ntrials - temperature and then the energy and energy min value (from 1.10) + std::cout << "\tx = ( "; + unsigned n = NDim(); + for (unsigned int i = 0; i < n-1; ++i) { + std::cout << fX[i] << " , "; + } + std::cout << fX.back() << " )\t"; + // energy us printed by GSL (and also end-line) + std::cout << "E / E_best = "; // GSL print then E and E best +} + +GSLSimAnFunc & GSLSimAnFunc::FastCopy(const GSLSimAnFunc & rhs) { + // copy only the information which is changed during the process + // in this case only the x values + std::copy(rhs.fX.begin(), rhs.fX.end(), fX.begin() ); + return *this; +} + + + +// definition and implementations of the static functions required by GSL + +namespace GSLSimAn { + + + double E( void * xp) { + // evaluate the energy given a state xp + GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp); + assert (fx != 0); + return fx->Energy(); + } + + void Step( const gsl_rng * r, void * xp, double step_size) { + // change xp according to the random step + GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp); + assert (fx != 0); + // create GSLRandomEngine class + // cast away constness (we make sure we don't delete (no call to Terminate() ) + GSLRngWrapper rng(const_cast<gsl_rng *>(r)); + GSLRandomEngine random(&rng); + // wrapper classes random and rng must exist during call to Step() + fx->Step(random, step_size); + } + + double Dist( void * xp, void * yp) { + // calculate the distance between two configuration + GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp); + GSLSimAnFunc * fy = reinterpret_cast<GSLSimAnFunc *> (yp); + + assert (fx != 0); + assert (fy != 0); + return fx->Distance(*fy); + } + + void Print(void * xp) { + // print the position xp + // GSL prints also first niter- ntrials - temperature and then the energy + GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp); + assert (fx != 0); + fx->Print(); + } + +// static function to pass to GSL copy - create and destroy the object + + void Copy( void * source, void * dest) { + GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (source); + assert (fx != 0); + GSLSimAnFunc * gx = reinterpret_cast<GSLSimAnFunc *> (dest); + assert (gx != 0); + gx->FastCopy(*fx); + } + + void * CopyCtor( void * xp) { + GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp); + assert (fx != 0); + return static_cast<void *> ( fx->Clone() ); + } + + void Destroy( void * xp) { + GSLSimAnFunc * fx = reinterpret_cast<GSLSimAnFunc *> (xp); + assert (fx != 0); + delete fx; + } + +} + +// implementation of GSLSimAnnealing class + + +GSLSimAnnealing::GSLSimAnnealing() +{ + // Default constructor implementation. +} + + + +// function for solving (from a Genfunction interface) + +int GSLSimAnnealing::Solve(const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * scale, double * xmin, bool debug) { + // solve the simulated annealing problem given starting point and objective function interface + + + // initial conditions + GSLSimAnFunc fx(func, x0, scale); + + int iret = Solve(fx, debug); + + if (iret == 0) { + // copy value of the minimum in xmin + std::copy(fx.X().begin(), fx.X().end(), xmin); + } + return iret; + +} + +int GSLSimAnnealing::Solve(GSLSimAnFunc & fx, bool debug) { + // solve the simulated annealing problem given starting point and GSLSimAnfunc object + + gsl_rng * r = gsl_rng_alloc(gsl_rng_mt19937); + + + + gsl_siman_params_t simanParams; + + // parameters for the simulated annealing + // copy them in GSL structure + + simanParams.n_tries = fParams.n_tries; /* how many points to try for each step */ + simanParams.iters_fixed_T = fParams.iters_fixed_T; /* how many iterations at each temperature? */ + simanParams.step_size = fParams.step_size; /* max step size in the random walk */ + // the following parameters are for the Boltzmann distribution */ + simanParams.k = fParams.k; + simanParams.t_initial = fParams.t_initial; + simanParams.mu_t = fParams.mu; + simanParams.t_min = fParams.t_min; + + + if (debug) + gsl_siman_solve(r, &fx, &GSLSimAn::E, &GSLSimAn::Step, &GSLSimAn::Dist, + &GSLSimAn::Print, &GSLSimAn::Copy, &GSLSimAn::CopyCtor , &GSLSimAn::Destroy, 0, simanParams ); + + else + gsl_siman_solve(r, &fx, &GSLSimAn::E, &GSLSimAn::Step, &GSLSimAn::Dist, + 0, &GSLSimAn::Copy, &GSLSimAn::CopyCtor , &GSLSimAn::Destroy, 0, simanParams ); + + return 0; + +} + + + + + } // end namespace Math + +} // end namespace ROOT + diff --git a/ThirdParty/RootMinimizers/src/Math/MultiNumGradFunction.cxx b/ThirdParty/RootMinimizers/src/Math/MultiNumGradFunction.cxx new file mode 100644 index 00000000000..1cecdca1069 --- /dev/null +++ b/ThirdParty/RootMinimizers/src/Math/MultiNumGradFunction.cxx @@ -0,0 +1,61 @@ +// @(#)root/mathmore:$Id$ +// Author: L. Moneta Wed Dec 20 14:36:31 2006 + +/********************************************************************** + * * + * Copyright (c) 2006 LCG ROOT Math Team, CERN/PH-SFT * + * * + * This library is free software; you can redistribute it and/or * + * modify it under the terms of the GNU General Public License * + * as published by the Free Software Foundation; either version 2 * + * of the License, or (at your option) any later version. * + * * + * This library is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * + * General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this library (see file COPYING); if not, write * + * to the Free Software Foundation, Inc., 59 Temple Place, Suite * + * 330, Boston, MA 02111-1307 USA, or contact the author. * + * * + **********************************************************************/ + +// implementation file for class MultiNumGradFunction + +#include "Math/MultiNumGradFunction.h" +#include <limits> +#include <cmath> +#include <algorithm> // needed for std::max on Solaris + +#ifndef ROOT_Math_Derivator +#include "Math/Derivator.h" +#endif + + +namespace ROOT { + + namespace Math { + + +double MultiNumGradFunction::fgEps = 0.001; + +double MultiNumGradFunction::DoDerivative (const double * x, unsigned int icoord ) const { + // calculate derivative using mathcore derivator class + // step size can be changes using SetDerivPrecision() + + static double kPrecision = std::sqrt ( std::numeric_limits<double>::epsilon() ); + double x0 = x[icoord]; + double step = std::max( fgEps* std::abs(x0), 8.0*kPrecision*(std::abs(x0) + kPrecision) ); + return ROOT::Math::Derivator::Eval(*fFunc, x, icoord, step); +} + +void MultiNumGradFunction::SetDerivPrecision(double eps) { fgEps = eps; } + +double MultiNumGradFunction::GetDerivPrecision( ) { return fgEps; } + + + } // end namespace Math + +} // end namespace ROOT diff --git a/shared.pri b/shared.pri index cde8aa49c46..0ea0548bd00 100644 --- a/shared.pri +++ b/shared.pri @@ -215,7 +215,7 @@ macx|unix { ROOT_FRAMEWORK = $$system(root-config --prefix) } win32 { - ROOT_FRAMEWORK = "C:/root" +# ROOT_FRAMEWORK = "C:/root" } # ----------------------------------------------------------------------------- -- GitLab