From 0d7fb8011f5e2176b0aa20b701ab5ccb390b72a4 Mon Sep 17 00:00:00 2001
From: Gennady Pospelov <g.pospelov@fz-juelich.de>
Date: Thu, 21 Apr 2016 15:42:02 +0200
Subject: [PATCH] Updates in ROOT minimizers
---
.../RootMinimizers/inc/Math/Derivator.h | 4 -
.../inc/Math/GSLRandomFunctions.h | 212 ++++++++++
.../RootMinimizers/inc/Math/GSLRndmEngines.h | 24 ++
.../RootMinimizers/inc/Math/RandomFunctions.h | 303 +++++++++++++++
.../RootMinimizers/inc/Math/TRandomEngine.h | 74 ++++
.../inc/Minuit2/BasicMinimumParameters.h | 2 +-
.../inc/Minuit2/BasicMinimumState.h | 19 +-
.../RootMinimizers/inc/Minuit2/FCNBase.h | 10 +-
.../inc/Minuit2/MinimumParameters.h | 9 +-
.../RootMinimizers/inc/Minuit2/MinimumState.h | 13 +-
.../inc/Minuit2/VariableMetricBuilder.h | 3 +-
.../RootMinimizers/src/Math/GSLMultiFit.h | 14 +-
.../src/Math/GSLRndmEngines.cxx | 43 ++-
.../src/Math/RandomFunctions.cxx | 365 ++++++++++++++++++
.../src/Minuit2/VariableMetricBuilder.cxx | 67 ++--
cmake/modules/SearchInstalledSoftware.cmake | 5 +-
dev-tools/code-tools/update_minimizers.py | 2 +-
17 files changed, 1099 insertions(+), 70 deletions(-)
create mode 100644 ThirdParty/RootMinimizers/inc/Math/GSLRandomFunctions.h
create mode 100644 ThirdParty/RootMinimizers/inc/Math/RandomFunctions.h
create mode 100644 ThirdParty/RootMinimizers/inc/Math/TRandomEngine.h
create mode 100644 ThirdParty/RootMinimizers/src/Math/RandomFunctions.cxx
diff --git a/ThirdParty/RootMinimizers/inc/Math/Derivator.h b/ThirdParty/RootMinimizers/inc/Math/Derivator.h
index fffe6882dbf..b56445d3c9f 100644
--- a/ThirdParty/RootMinimizers/inc/Math/Derivator.h
+++ b/ThirdParty/RootMinimizers/inc/Math/Derivator.h
@@ -33,10 +33,6 @@
#ifndef ROOT_Math_Derivator
#define ROOT_Math_Derivator
-/**
- @defgroup Deriv Numerical Differentiation
- @ingroup NumAlgo
-*/
#ifndef ROOT_Math_IFunctionfwd
#include "Math/IFunctionfwd.h"
diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLRandomFunctions.h b/ThirdParty/RootMinimizers/inc/Math/GSLRandomFunctions.h
new file mode 100644
index 00000000000..2d7d380493b
--- /dev/null
+++ b/ThirdParty/RootMinimizers/inc/Math/GSLRandomFunctions.h
@@ -0,0 +1,212 @@
+// @(#)root/mathmore:$Id$
+// Authors: L. Moneta 8/2015
+
+/**********************************************************************
+ * *
+ * Copyright (c) 2015 , ROOT MathLib Team *
+ * *
+ * *
+ **********************************************************************/
+
+// Header file for random class
+//
+//
+// Created by: Lorenzo Moneta : Tue 4 Aug 2015
+//
+//
+#ifndef ROOT_Math_GSLRandomFunctions
+#define ROOT_Math_GSLRandomFunctions
+
+
+//#include <type_traits>
+
+#include "Math/RandomFunctions.h"
+
+#include "Math/GSLRndmEngines.h"
+
+namespace BA_ROOT {
+namespace Math {
+
+
+//___________________________________________________________________________________
+ /**
+ Specialized implementation of the Random functions based on the GSL library.
+ These will work onlmy with a GSLRandomEngine type
+
+ @ingroup Random
+ */
+
+
+ template <class EngineType >
+ class RandomFunctions<EngineType, BA_ROOT::Math::GSLRandomEngine> : public RandomFunctions<EngineType, DefaultEngineType> {
+ //class RandomFunctions<Engine, ROOT::Math::GSLRandomEngine> {
+
+ //typdef TRandomEngine DefaulEngineType;
+
+ public:
+
+ RandomFunctions() {}
+
+ RandomFunctions(EngineType & rng) : RandomFunctions<EngineType, DefaultEngineType>(rng) {}
+
+
+ inline EngineType & Engine() { return RandomFunctions<EngineType,DefaultEngineType>::Rng(); }
+
+ double GausZig(double mean, double sigma) {
+ return Engine().GaussianZig(sigma) + mean;
+ }
+ // double GausRatio(double mean, double sigma) {
+ // auto & r = RandomFunctions<Engine,DefaultEngineType>::Rng();
+ // return r.GaussianRatio(sigma) + mean;
+ // }
+
+ /**
+ Gaussian distribution. Default method (use Ziggurat)
+ */
+ double Gaus(double mean = 0, double sigma = 1) {
+ return mean + Engine().GaussianZig(sigma);
+ }
+
+ /**
+ Gaussian distribution (Box-Muller method)
+ */
+ double GausBM(double mean = 0, double sigma = 1) {
+ return mean + Engine().Gaussian(sigma);
+ }
+
+ /**
+ Gaussian distribution (Ratio Method)
+ */
+ double GausR(double mean = 0, double sigma = 1) {
+ return mean + Engine().GaussianRatio(sigma);
+ }
+
+ /**
+ Gaussian Tail distribution
+ */
+ double GaussianTail(double a, double sigma = 1) {
+ return Engine().GaussianTail(a,sigma);
+ }
+
+ /**
+ Bivariate Gaussian distribution with correlation
+ */
+ void Gaussian2D(double sigmaX, double sigmaY, double rho, double &x, double &y) {
+ Engine().Gaussian2D(sigmaX, sigmaY, rho, x, y);
+ }
+
+ /**
+ Exponential distribution
+ */
+ double Exp(double tau) {
+ return Engine().Exponential(tau);
+ }
+ /**
+ Breit Wigner distribution
+ */
+ double BreitWigner(double mean = 0., double gamma = 1) {
+ return mean + Engine().Cauchy( gamma/2.0 );
+ }
+
+ /**
+ Landau distribution
+ */
+ double Landau(double mean = 0, double sigma = 1) {
+ return mean + sigma*Engine().Landau();
+ }
+
+ /**
+ Gamma distribution
+ */
+ double Gamma(double a, double b) {
+ return Engine().Gamma(a,b);
+ }
+
+ /**
+ Log Normal distribution
+ */
+ double LogNormal(double zeta, double sigma) {
+ return Engine().LogNormal(zeta,sigma);
+ }
+
+ /**
+ Chi square distribution
+ */
+ double ChiSquare(double nu) {
+ return Engine().ChiSquare(nu);
+ }
+
+ /**
+ F distrbution
+ */
+ double FDist(double nu1, double nu2) {
+ return Engine().FDist(nu1,nu2);
+ }
+
+ /**
+ t student distribution
+ */
+ double tDist(double nu) {
+ return Engine().tDist(nu);
+ }
+
+ /**
+ generate random numbers in a 2D circle of radious 1
+ */
+ void Circle(double &x, double &y, double r = 1) {
+ Engine().Dir2D(x,y);
+ x *= r;
+ y *= r;
+ }
+
+ /**
+ generate random numbers in a 3D sphere of radious 1
+ */
+ void Sphere(double &x, double &y, double &z,double r = 1) {
+ Engine().Dir3D(x,y,z);
+ x *= r;
+ y *= r;
+ z *= r;
+ }
+
+ /**
+ Poisson distribution
+ */
+ unsigned int Poisson(double mu) {
+ return Engine().Poisson(mu);
+ }
+
+ /**
+ Binomial distribution
+ */
+ unsigned int Binomial(unsigned int ntot, double prob) {
+ return Engine().Binomial(prob,ntot);
+ }
+
+ /**
+ Negative Binomial distribution
+ First parameter is n, second is probability
+ To be consistent with Random::Binomial
+ */
+ unsigned int NegativeBinomial(double n, double prob) {
+ return Engine().NegativeBinomial(prob,n);
+ }
+
+ /**
+ Multinomial distribution
+ */
+ std::vector<unsigned int> Multinomial( unsigned int ntot, const std::vector<double> & p ) {
+ return Engine().Multinomial(ntot,p);
+ }
+
+
+
+ };
+
+
+
+
+} // namespace Math
+} // namespace ROOT
+
+#endif /* ROOT_Math_GSLRandomFunctions */
diff --git a/ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h b/ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h
index 20df9962bea..0772437fbd8 100644
--- a/ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h
+++ b/ThirdParty/RootMinimizers/inc/Math/GSLRndmEngines.h
@@ -34,6 +34,7 @@
#include <string>
#include <vector>
+
namespace BA_ROOT {
namespace Math {
@@ -108,6 +109,12 @@ namespace Math {
0 is excluded and 1 is included
*/
double operator() () const;
+
+ /**
+ Generate a random number between ]0,1]
+ 0 is excluded and 1 is included
+ */
+ double Rndm() const { return (*this)(); }
/**
Generate an integer number between [0,max-1] (including 0 and max-1)
@@ -279,6 +286,7 @@ namespace Math {
*/
class GSLRngMT : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngMT();
};
@@ -292,6 +300,7 @@ namespace Math {
*/
class GSLRngRanLux : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngRanLux();
};
@@ -305,6 +314,7 @@ namespace Math {
*/
class GSLRngRanLuxS1 : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngRanLuxS1();
};
typedef GSLRngRanLuxS1 GSLRngRanLux1; // for backward compatibility
@@ -319,6 +329,7 @@ namespace Math {
*/
class GSLRngRanLuxS2 : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngRanLuxS2();
};
typedef GSLRngRanLuxS2 GSLRngRanLux2; // for backward compatibility
@@ -333,6 +344,7 @@ namespace Math {
*/
class GSLRngRanLuxD1 : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngRanLuxD1();
};
@@ -346,6 +358,7 @@ namespace Math {
*/
class GSLRngRanLuxD2 : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngRanLuxD2();
};
typedef GSLRngRanLuxD2 GSLRngRanLux48; // for backward compatibility
@@ -360,6 +373,7 @@ namespace Math {
*/
class GSLRngTaus : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngTaus();
};
@@ -372,6 +386,7 @@ namespace Math {
*/
class GSLRngGFSR4 : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngGFSR4();
};
@@ -384,6 +399,7 @@ namespace Math {
*/
class GSLRngCMRG : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngCMRG();
};
@@ -396,6 +412,7 @@ namespace Math {
*/
class GSLRngMRG : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngMRG();
};
@@ -409,6 +426,7 @@ namespace Math {
*/
class GSLRngRand : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngRand();
};
@@ -421,6 +439,7 @@ namespace Math {
*/
class GSLRngRanMar : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngRanMar();
};
@@ -433,6 +452,7 @@ namespace Math {
*/
class GSLRngMinStd : public GSLRandomEngine {
public:
+ typedef GSLRandomEngine BaseType;
GSLRngMinStd();
};
@@ -442,6 +462,10 @@ namespace Math {
} // namespace Math
} // namespace ROOT
+// random functions specialization for GSL
+// needs to be defined after defining GSLRandomEngine class
+
+#include "Math/GSLRandomFunctions.h"
#endif /* ROOT_Math_GSLRndmEngines */
diff --git a/ThirdParty/RootMinimizers/inc/Math/RandomFunctions.h b/ThirdParty/RootMinimizers/inc/Math/RandomFunctions.h
new file mode 100644
index 00000000000..3348b7a9425
--- /dev/null
+++ b/ThirdParty/RootMinimizers/inc/Math/RandomFunctions.h
@@ -0,0 +1,303 @@
+// @(#)root/mathcore:$Id$
+// Authors: L. Moneta 8/2015
+
+/**********************************************************************
+ * *
+ * Copyright (c) 2015 , ROOT MathLib Team *
+ * *
+ * *
+ **********************************************************************/
+
+// Header file for random class
+//
+//
+// Created by: Lorenzo Moneta : Tue 4 Aug 2015
+//
+//
+#ifndef ROOT_Math_RandomFunctions
+#define ROOT_Math_RandomFunctions
+
+
+#include <type_traits>
+#include <cmath>
+//#include "Rtypes.h"
+//#include "TMath.h"
+#include "TMVA/Types.h"
+#include <cassert>
+#include <vector>
+
+#include "TRandomEngine.h"
+
+
+namespace BA_ROOT {
+namespace Math {
+
+
+//___________________________________________________________________________________
+
+
+ // class DefaultEngineType {};
+
+
+ /**
+ Documentation for the RandomFunction class
+
+ @ingroup Random
+ */
+
+
+ typedef TRandomEngine DefaultEngineType;
+ //class DefaultEngineType {}; // for generic types
+
+
+
+ /**
+ Definition of the generic impelmentation class for the RandomFunctions.
+ Needs to have specialized implementations on the different type of engines
+ */
+ template <class EngineBaseType>
+ class RandomFunctionsImpl {
+ public:
+ void SetEngine(void *) {}
+ };
+
+ /**
+ Implementation class for the RandomFunction for all the engined that derives from
+ TRandomEngine class, which defines an interface which has TRandomEngine::Rndm()
+ In this way we can have a common implementation for the RandomFunctions
+ */
+
+ template<>
+ class RandomFunctionsImpl<TRandomEngine> {
+
+ public:
+
+ /// class constructor
+ RandomFunctionsImpl() {}
+
+ void SetEngine(void *r) {
+ fBaseEngine = static_cast<TRandomEngine*>(r);
+ assert(fBaseEngine); // to be sure the static cast works
+ }
+
+
+ ///Generate binomial numbers
+ int Binomial(int ntot, double prob);
+
+ /// Return a number distributed following a BreitWigner function with mean and gamma.
+ double BreitWigner(double mean, double gamma);
+
+ /// Generates random vectors, uniformly distributed over a circle of given radius.
+ /// Input : r = circle radius
+ /// Output: x,y a random 2-d vector of length r
+ void Circle(double &x, double &y, double r);
+
+ /// Returns an exponential deviate.
+ /// exp( -t/tau )
+ double Exp(double tau);
+
+ /// generate Gaussian number using Box-Muller method
+ double GausBM( double mean, double sigma);
+
+ /// generate random numbers according to the Accemptance-Complemet-Ratio method
+ double GausACR( double mean, double sigma);
+
+ /// Generate a random number following a Landau distribution
+ /// with location parameter mu and scale parameter sigma:
+ /// Landau( (x-mu)/sigma )
+// double Landau(double mu, double sigma);
+
+ /// Generates a random integer N according to a Poisson law.
+ /// Prob(N) = exp(-mean)*mean^N/Factorial(N)
+// int Poisson(double mean);
+// double PoissonD(double mean);
+
+ /// Generate numbers distributed following a gaussian with mean=0 and sigma=1.
+ /// Using the Box-Muller method
+ void Rannor(double &a, double &b);
+
+ /// Generates random vectors, uniformly distributed over the surface
+ /// of a sphere of given radius.
+ void Sphere(double &x, double &y, double &z, double r);
+
+ /// generate random numbers following a Uniform distribution in the [a,b] interval
+ double Uniform(double a, double b);
+ double Uniform(double a);
+
+ protected:
+ TRandomEngine * fBaseEngine;
+
+ private:
+ // Internal method used by the functions
+ double Rndm() { return fBaseEngine->Rndm(); }
+ // for internal usage
+ double Gaus(double mean, double sigma) { return GausACR(mean,sigma); }
+
+
+ };
+
+
+ template < class Engine, class EngineBaseType>
+ class RandomFunctions { //: public RandomFunctionsImpl<EngineBaseType> {
+
+
+ public:
+
+ //RandomFunctions() {}
+
+ RandomFunctions(Engine & rng) : fEngine(&rng) {
+ fImpl.SetEngine(&rng);
+ }
+
+ /// destructor (no op) we do not mantain the engine)
+ ~RandomFunctions() {}
+
+
+ /// non-virtual method
+ inline double operator() () { return (*fEngine)(); }
+
+
+ ///Generate binomial numbers
+ int Binomial(int ntot, double prob) {
+ return fImpl.Binomial(ntot,prob);
+ }
+
+ /// Return a number distributed following a BreitWigner function with mean and gamma.
+ double BreitWigner(double mean, double gamma) {
+ return fImpl.BreitWigner(mean,gamma);
+ }
+
+ /// Generates random vectors, uniformly distributed over a circle of given radius.
+ /// Input : r = circle radius
+ /// Output: x,y a random 2-d vector of length r
+ void Circle(double &x, double &y, double r) {
+ return fImpl.Circle(x,y,r);
+ }
+
+ /// Returns an exponential deviate.
+ /// exp( -t/tau )
+ double Exp(double tau) {
+ return fImpl.Exp(tau);
+ }
+
+ /// generate Gaussian number using Box-Muller method
+ double GausBM( double mean, double sigma) {
+ return fImpl.GausBM(mean,sigma);
+ }
+
+ /// generate random numbers according to the Accemptance-Complemet-Ratio method
+ double GausACR( double mean, double sigma) {
+ return fImpl.GausACR(mean, sigma);
+ }
+
+ /// Generate a random number following a Landau distribution
+ /// with location parameter mu and scale parameter sigma:
+ /// Landau( (x-mu)/sigma )
+ double Landau(double mu, double sigma) {
+ return fImpl.Landau(mu,sigma);
+ }
+
+ /// Generates a random integer N according to a Poisson law.
+ /// Prob(N) = exp(-mean)*mean^N/Factorial(N)
+ int Poisson(double mean) { return fImpl.Poisson(mean); }
+ double PoissonD(double mean) { return fImpl.PoissonD(mean); }
+
+ /// Generate numbers distributed following a gaussian with mean=0 and sigma=1.
+ /// Using the Box-Muller method
+ void Rannor(double &a, double &b) {
+ return fImpl.Rannor(a,b);
+ }
+
+ /// Generates random vectors, uniformly distributed over the surface
+ /// of a sphere of given radius.
+ void Sphere(double &x, double &y, double &z, double r) {
+ return fImpl.Sphere(x,y,z,r);
+ }
+
+ /// generate random numbers following a Uniform distribution in the [a,b] interval
+ double Uniform(double a, double b) {
+ return (b-a) * Rndm_impl() + a;
+ }
+
+ /// generate random numbers following a Uniform distribution in the [0,a] interval
+ double Uniform(double a) {
+ return a * Rndm_impl() ;
+ }
+
+
+ /// generate Gaussian number using defqault method
+ inline double Gaus( double mean, double sigma) {
+ return fImpl.GausACR(mean,sigma);
+ }
+
+
+ // /// re-implement Gaussian
+ // double GausBM2(double mean, double sigma) {
+ // double y = Rndm_impl();
+ // double z = Rndm_impl();
+ // double x = z * 6.28318530717958623;
+ // double radius = std::sqrt(-2*std::log(y));
+ // double g = radius * std::sin(x);
+ // return mean + g * sigma;
+ // }
+
+
+ /// methods which are only for GSL random generators
+
+
+ /// Gamma functions (not implemented here, requires a GSL random engine)
+ double Gamma( double , double ) {
+ //r.Error("Error: Gamma() requires a GSL Engine type");
+ static_assert(std::is_fundamental<Engine>::value,"Error: Gamma() requires a GSL Engine type");
+ return 0;
+ }
+ double LogNormal(double, double) {
+ static_assert(std::is_fundamental<Engine>::value,"Error: LogNormal() requires a GSL Engine type");
+ return 0;
+ }
+ double ChiSquare(double) {
+ static_assert(std::is_fundamental<Engine>::value,"Error: ChiSquare() requires a GSL Engine type");
+ return 0;
+ }
+ double FDist(double, double) {
+ static_assert(std::is_fundamental<Engine>::value,"Error: FDist() requires a GSL Engine type");
+ return 0;
+ }
+ double tDist(double) {
+ static_assert(std::is_fundamental<Engine>::value,"Error: tDist() requires a GSL Engine type");
+ return 0;
+ }
+ unsigned int NegativeBinomial(double , double ) {
+ static_assert(std::is_fundamental<Engine>::value,"Error: NegativeBinomial() requires a GSL Engine type");
+ return 0;
+ }
+ std::vector<unsigned int> MultiNomial(unsigned int, const std::vector<double> &){
+ static_assert(std::is_fundamental<Engine>::value,"Error: MultiNomial() requires a GSL Engine type");
+ return std::vector<unsigned int>();
+ }
+
+
+ protected:
+
+ Engine & Rng() { assert(fEngine); return *fEngine; }
+
+ /// Internal impelmentation to return random number
+ /// Since this one is not a virtual function is faster than Rndm
+ inline double Rndm_impl() { return (*fEngine)(); }
+
+
+ private:
+
+ Engine * fEngine; //! random number generator engine
+ RandomFunctionsImpl<EngineBaseType> fImpl; //! instance of the class implementing the functions
+
+
+ };
+
+
+
+
+} // namespace Math
+} // namespace ROOT
+
+#endif /* ROOT_Math_RandomFunctions */
diff --git a/ThirdParty/RootMinimizers/inc/Math/TRandomEngine.h b/ThirdParty/RootMinimizers/inc/Math/TRandomEngine.h
new file mode 100644
index 00000000000..8c037d33495
--- /dev/null
+++ b/ThirdParty/RootMinimizers/inc/Math/TRandomEngine.h
@@ -0,0 +1,74 @@
+// @(#)root/mathcore:$Id$
+// Author: L. Moneta Tue Aug 4 2015
+
+/**********************************************************************
+ * *
+ * Copyright (c) 2015 LCG ROOT Math Team, CERN/PH-SFT *
+ * *
+ * *
+ **********************************************************************/
+
+// random engines based on ROOT
+
+#ifndef ROOT_Math_TRandomEngine
+#define ROOT_Math_TRandomEngine
+
+
+
+namespace BA_ROOT {
+
+ namespace Math {
+
+ class RandomBaseEngine {
+ public:
+ virtual double Rndm() = 0;
+ virtual ~RandomBaseEngine() {}
+ };
+
+
+ class TRandomEngine : public RandomBaseEngine {
+ public:
+ virtual ~TRandomEngine() {}
+ };
+
+ class LCGEngine : public TRandomEngine {
+
+
+ public:
+
+ typedef TRandomEngine BaseType;
+
+ LCGEngine() : fSeed(65539) { }
+
+ virtual ~LCGEngine() {}
+
+ void SetSeed(unsigned int seed) { fSeed = seed; }
+
+ virtual double Rndm() {
+ //double Rndm() {
+ return Rndm_impl();
+ }
+ double Rndm_impl() {
+ const double kCONS = 4.6566128730774E-10; // (1/pow(2,31)
+ unsigned int rndm = IntRndm(); // generate integer number
+ if (rndm != 0) return kCONS*rndm;
+ return Rndm_impl();
+ }
+ inline double operator() () { return Rndm_impl(); }
+
+ unsigned int IntRndm() {
+ fSeed = (1103515245 * fSeed + 12345) & 0x7fffffffUL;
+ return fSeed;
+ }
+
+ private:
+ unsigned int fSeed;
+ };
+
+
+ } // end namespace Math
+
+} // end namespace ROOT
+
+
+#endif /* ROOT_Math_TRandomEngine */
diff --git a/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumParameters.h b/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumParameters.h
index 150c11ee768..b9629f272e2 100644
--- a/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumParameters.h
+++ b/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumParameters.h
@@ -25,7 +25,7 @@ class BasicMinimumParameters {
public:
- BasicMinimumParameters(unsigned int n) : fParameters(MnAlgebraicVector(n)), fStepSize(MnAlgebraicVector(n)), fFVal(0.), fValid(false), fHasStep(false) {}
+ BasicMinimumParameters(unsigned int n, double fval) : fParameters(MnAlgebraicVector(n)), fStepSize(MnAlgebraicVector(n)), fFVal(fval), fValid(false), fHasStep(false) {}
BasicMinimumParameters(const MnAlgebraicVector& avec, double fval) :
fParameters(avec), fStepSize(avec.size()), fFVal(fval), fValid(true), fHasStep(false) {}
diff --git a/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumState.h b/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumState.h
index 57f042e7e4c..00a285be220 100644
--- a/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumState.h
+++ b/ThirdParty/RootMinimizers/inc/Minuit2/BasicMinimumState.h
@@ -27,19 +27,23 @@ class BasicMinimumState {
public:
- BasicMinimumState(unsigned int n) :
- fParameters(MinimumParameters(n)), fError(MinimumError(n)),
- fGradient(FunctionGradient(n)), fEDM(0.), fNFcn(0) {}
+ // constructor without parameter values but with function value, edm and nfcn
+ BasicMinimumState(unsigned int n, double fval, double edm, int nfcn) :
+ fParameters(MinimumParameters(n,fval)), fError(MinimumError(n)),
+ fGradient(FunctionGradient(n)), fEDM(edm), fNFcn(nfcn) {}
+
BasicMinimumState(const MinimumParameters& states, const MinimumError& err,
- const FunctionGradient& grad, double edm, int nfcn) :
- fParameters(states), fError(err), fGradient(grad), fEDM(edm), fNFcn(nfcn) {}
+ const FunctionGradient& grad, double edm, int nfcn) :
+ fParameters(states), fError(err), fGradient(grad), fEDM(edm), fNFcn(nfcn) {}
- BasicMinimumState(const MinimumParameters& states, double edm, int nfcn) : fParameters(states), fError(MinimumError(states.Vec().size())), fGradient(FunctionGradient(states.Vec().size())), fEDM(edm), fNFcn(nfcn) {}
+ BasicMinimumState(const MinimumParameters& states, double edm, int nfcn) : fParameters(states), fError(MinimumError(states.Vec().size())),
+ fGradient(FunctionGradient(states.Vec().size())), fEDM(edm), fNFcn(nfcn)
+ {}
~BasicMinimumState() {}
BasicMinimumState(const BasicMinimumState& state) :
- fParameters(state.fParameters), fError(state.fError), fGradient(state.fGradient), fEDM(state.fEDM), fNFcn(state.fNFcn) {}
+ fParameters(state.fParameters), fError(state.fError), fGradient(state.fGradient), fEDM(state.fEDM), fNFcn(state.fNFcn) {}
BasicMinimumState& operator=(const BasicMinimumState& state) {
fParameters = state.fParameters;
@@ -68,6 +72,7 @@ public:
double Edm() const {return fEDM;}
int NFcn() const {return fNFcn;}
+
bool IsValid() const {
if(HasParameters() && HasCovariance())
return Parameters().IsValid() && Error().IsValid();
diff --git a/ThirdParty/RootMinimizers/inc/Minuit2/FCNBase.h b/ThirdParty/RootMinimizers/inc/Minuit2/FCNBase.h
index 61095159f3e..c31c9445b35 100644
--- a/ThirdParty/RootMinimizers/inc/Minuit2/FCNBase.h
+++ b/ThirdParty/RootMinimizers/inc/Minuit2/FCNBase.h
@@ -23,10 +23,12 @@ namespace BA_ROOT {
/**
-@defgroup Minuit Minuit2 Minimization Library
+ \defgroup Minuit Minuit2 Minimization Library
- Object-oriented implementation of the MINUIT minimization package.
- More information is available at the home page of the \ref Minuit2 package.
+ New Object-oriented implementation of the MINUIT minimization package.
+ More information is available at the home page of the \ref Minuit2Page "Minuit2" minimization package".
+
+ \ingroup Math
*/
@@ -38,7 +40,7 @@ Interface (abstract class) defining the function to be minimized, which has to b
@author Fred James and Matthias Winkler; modified by Andras Zsenei and Lorenzo Moneta
-@ingroup Minuit
+\ingroup Minuit
*/
diff --git a/ThirdParty/RootMinimizers/inc/Minuit2/MinimumParameters.h b/ThirdParty/RootMinimizers/inc/Minuit2/MinimumParameters.h
index db500424212..b6688fd9bec 100644
--- a/ThirdParty/RootMinimizers/inc/Minuit2/MinimumParameters.h
+++ b/ThirdParty/RootMinimizers/inc/Minuit2/MinimumParameters.h
@@ -22,15 +22,16 @@ class MinimumParameters {
public:
- MinimumParameters(unsigned int n) :
- fData(MnRefCountedPointer<BasicMinimumParameters>(new BasicMinimumParameters(n))) {}
+ MinimumParameters(unsigned int n, double fval = 0) :
+ fData(MnRefCountedPointer<BasicMinimumParameters>(new BasicMinimumParameters(n,fval))) {}
/** takes the Parameter vector */
MinimumParameters(const MnAlgebraicVector& avec, double fval) :
- fData(MnRefCountedPointer<BasicMinimumParameters>(new BasicMinimumParameters(avec, fval))) {}
+ fData(MnRefCountedPointer<BasicMinimumParameters>(new BasicMinimumParameters(avec, fval))) {}
/** takes the Parameter vector plus step size x1 - x0 = dirin */
- MinimumParameters(const MnAlgebraicVector& avec, const MnAlgebraicVector& dirin, double fval) : fData(MnRefCountedPointer<BasicMinimumParameters>(new BasicMinimumParameters(avec, dirin, fval))) {}
+ MinimumParameters(const MnAlgebraicVector& avec, const MnAlgebraicVector& dirin, double fval) :
+ fData(MnRefCountedPointer<BasicMinimumParameters>(new BasicMinimumParameters(avec, dirin, fval))) {}
~MinimumParameters() {}
diff --git a/ThirdParty/RootMinimizers/inc/Minuit2/MinimumState.h b/ThirdParty/RootMinimizers/inc/Minuit2/MinimumState.h
index bce046d2a75..8168dc34d53 100644
--- a/ThirdParty/RootMinimizers/inc/Minuit2/MinimumState.h
+++ b/ThirdParty/RootMinimizers/inc/Minuit2/MinimumState.h
@@ -31,17 +31,20 @@ class MinimumState {
public:
/** invalid state */
- MinimumState(unsigned int n) :
- fData(MnRefCountedPointer<BasicMinimumState>(new BasicMinimumState(n))) {}
+ MinimumState(unsigned int n) :
+ fData(MnRefCountedPointer<BasicMinimumState>(new BasicMinimumState(n,0.,0.,0.))) {}
+ /** state without parameters and errors (only function value an, edm and nfcn) */
+ MinimumState(double fval, double edm, int nfcn) :
+ fData(MnRefCountedPointer<BasicMinimumState>(new BasicMinimumState(0, fval, edm, nfcn))) {}
/** state with parameters only (from stepping methods like Simplex, Scan) */
- MinimumState(const MinimumParameters& states, double edm, int nfcn) :
- fData(MnRefCountedPointer<BasicMinimumState>(new BasicMinimumState(states, edm, nfcn))) {}
+ MinimumState(const MinimumParameters& states, double edm, int nfcn) :
+ fData(MnRefCountedPointer<BasicMinimumState>(new BasicMinimumState(states, edm, nfcn))) {}
/** state with parameters, Gradient and covariance (from Gradient methods
such as Migrad) */
MinimumState(const MinimumParameters& states, const MinimumError& err,
const FunctionGradient& grad, double edm, int nfcn) :
- fData(MnRefCountedPointer<BasicMinimumState>(new BasicMinimumState(states, err, grad, edm, nfcn))) {}
+ fData(MnRefCountedPointer<BasicMinimumState>(new BasicMinimumState(states, err, grad, edm, nfcn))) {}
~MinimumState() {}
diff --git a/ThirdParty/RootMinimizers/inc/Minuit2/VariableMetricBuilder.h b/ThirdParty/RootMinimizers/inc/Minuit2/VariableMetricBuilder.h
index 5eb6c9afb12..89fa28aca30 100644
--- a/ThirdParty/RootMinimizers/inc/Minuit2/VariableMetricBuilder.h
+++ b/ThirdParty/RootMinimizers/inc/Minuit2/VariableMetricBuilder.h
@@ -23,6 +23,7 @@ namespace BA_ROOT {
/**
Build (find) function minimum using the Variable Metric method (MIGRAD)
+
*/
class VariableMetricBuilder : public MinimumBuilder {
@@ -41,7 +42,7 @@ public:
const VariableMetricEDMEstimator& Estimator() const {return fEstimator;}
const DavidonErrorUpdator& ErrorUpdator() const {return fErrorUpdator;}
- void AddResult(std::vector<MinimumState>& result, const MinimumState & state, bool store = false) const;
+ void AddResult(std::vector<MinimumState>& result, const MinimumState & state) const;
private:
diff --git a/ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h b/ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h
index 232a3c08834..c8194409566 100644
--- a/ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h
+++ b/ThirdParty/RootMinimizers/src/Math/GSLMultiFit.h
@@ -31,6 +31,7 @@
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_multifit_nlin.h"
#include "gsl/gsl_blas.h"
+#include "gsl/gsl_version.h"
#include "GSLMultiFitFunctionWrapper.h"
#include "Math/IFunction.h"
@@ -143,10 +144,8 @@ public:
/// gradient value at the minimum
const double * Gradient() const {
if (fSolver == 0) return 0;
-#ifdef BORNAGAIN_GSL_BIGGEROREQUAL_2
- gsl_matrix * J = gsl_matrix_alloc(fSolver->fdf->n, fSolver->fdf->p);
- gsl_multifit_fdfsolver_jac(fSolver, J);
- gsl_multifit_gradient(J, fSolver->f,fVec);
+#if GSL_MAJOR_VERSION > 1
+ fType->gradient(fSolver->state, fVec);
#else
gsl_multifit_gradient(fSolver->J, fSolver->f,fVec);
#endif
@@ -160,10 +159,11 @@ public:
unsigned int npar = fSolver->fdf->p;
fCov = gsl_matrix_alloc( npar, npar );
static double kEpsrel = 0.0001;
-#ifdef BORNAGAIN_GSL_BIGGEROREQUAL_2
- gsl_matrix * J = gsl_matrix_alloc(fSolver->fdf->n, fSolver->fdf->p);
- gsl_multifit_fdfsolver_jac(fSolver, J);
+#if GSL_MAJOR_VERSION > 1
+ gsl_matrix* J = gsl_matrix_alloc(npar,npar);
+ gsl_multifit_fdfsolver_jac (fSolver, J);
int ret = gsl_multifit_covar(J, kEpsrel, fCov);
+ gsl_matrix_free(J);
#else
int ret = gsl_multifit_covar(fSolver->J, kEpsrel, fCov);
#endif
diff --git a/ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx b/ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx
index 85165c3c5c7..f5214b3b54b 100644
--- a/ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx
+++ b/ThirdParty/RootMinimizers/src/Math/GSLRndmEngines.cxx
@@ -159,14 +159,16 @@ namespace Math {
}
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() );
}
@@ -300,10 +302,12 @@ namespace Math {
// generators
//----------------------------------------------------
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngMT::GSLRngMT() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_mt19937));
+ Initialize();
}
@@ -311,73 +315,92 @@ namespace Math {
GSLRngRanLux::GSLRngRanLux() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_ranlux) );
+ Initialize();
}
// second generation of Ranlux (single precision version - luxury 1)
GSLRngRanLuxS1::GSLRngRanLuxS1() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_ranlxs1) );
+ Initialize();
}
// second generation of Ranlux (single precision version - luxury 2)
GSLRngRanLuxS2::GSLRngRanLuxS2() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_ranlxs2) );
+ Initialize();
}
// double precision version - luxury 1
GSLRngRanLuxD1::GSLRngRanLuxD1() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_ranlxd1) );
+ Initialize();
}
// double precision version - luxury 2
GSLRngRanLuxD2::GSLRngRanLuxD2() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_ranlxd2) );
+ Initialize();
}
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngTaus::GSLRngTaus() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_taus2) );
+ Initialize();
}
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngGFSR4::GSLRngGFSR4() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_gfsr4) );
+ Initialize();
}
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngCMRG::GSLRngCMRG() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_cmrg) );
+ Initialize();
}
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngMRG::GSLRngMRG() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_mrg) );
+ Initialize();
}
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngRand::GSLRngRand() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_rand) );
+ Initialize();
}
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngRanMar::GSLRngRanMar() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_ranmar) );
+ Initialize();
}
- //----------------------------------------------------
+ /////////////////////////////////////////////////////////////////////////////
+
GSLRngMinStd::GSLRngMinStd() : GSLRandomEngine()
{
SetType(new GSLRngWrapper(gsl_rng_minstd) );
+ Initialize();
}
diff --git a/ThirdParty/RootMinimizers/src/Math/RandomFunctions.cxx b/ThirdParty/RootMinimizers/src/Math/RandomFunctions.cxx
new file mode 100644
index 00000000000..bebca8cd5d3
--- /dev/null
+++ b/ThirdParty/RootMinimizers/src/Math/RandomFunctions.cxx
@@ -0,0 +1,365 @@
+// @(#)root/mathcore:$Id$
+// Authors: L. Moneta 8/2015
+
+/**********************************************************************
+ * *
+ * Copyright (c) 2015 , ROOT MathLib Team *
+ * *
+ * *
+ **********************************************************************/
+
+// file for random class
+//
+//
+// Created by: Lorenzo Moneta : Tue 4 Aug 2015
+//
+//
+#include "Math/RandomFunctions.h"
+
+//#include "Math/DistFunc.h"
+
+//#include "TMath.h"
+
+namespace BA_ROOT {
+namespace Math {
+
+
+Int_t RandomFunctionsImpl<TRandomEngine>::Binomial(Int_t ntot, Double_t prob)
+{
+ if (prob < 0 || prob > 1) return 0;
+ Int_t n = 0;
+ for (Int_t i=0;i<ntot;i++) {
+ if (Rndm() > prob) continue;
+ n++;
+ }
+ return n;
+}
+
+ ////////////////////////////////////////////////////////////////////////////////
+/// Return a number distributed following a BreitWigner function with mean and gamma.
+
+Double_t RandomFunctionsImpl<TRandomEngine>::BreitWigner(Double_t mean, Double_t gamma)
+{
+ Double_t rval, displ;
+ rval = 2*Rndm() - 1;
+ displ = 0.5*gamma*std::tan(rval*M_PI/2.);
+
+ return (mean+displ);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Generates random vectors, uniformly distributed over a circle of given radius.
+/// Input : r = circle radius
+/// Output: x,y a random 2-d vector of length r
+
+void RandomFunctionsImpl<TRandomEngine>::Circle(Double_t &x, Double_t &y, Double_t r)
+{
+ Double_t phi = Uniform(0,M_PI*2.0);
+ x = r*std::cos(phi);
+ y = r*std::sin(phi);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Returns an exponential deviate.
+///
+/// exp( -t/tau )
+
+Double_t RandomFunctionsImpl<TRandomEngine>::Exp(Double_t tau)
+{
+ Double_t x = Rndm(); // uniform on ] 0, 1 ]
+ Double_t t = -tau * std::log( x ); // convert to exponential distribution
+ return t;
+}
+
+
+
+
+double RandomFunctionsImpl<TRandomEngine>::GausBM(double mean, double sigma) {
+ double y = Rndm();
+ double z = Rndm();
+ double x = z * 6.28318530717958623;
+ double radius = std::sqrt(-2*std::log(y));
+ double g = radius * std::sin(x);
+ return mean + g * sigma;
+}
+
+
+ // double GausImpl(TRandomEngine * r, double mean, double sigma) {
+ // double y = r->Rndm();
+ // double z = r->Rndm();
+ // double x = z * 6.28318530717958623;
+ // double radius = std::sqrt(-2*std::log(y));
+ // double g = radius * std::sin(x);
+ // return mean + g * sigma;
+ // }
+
+double RandomFunctionsImpl<TRandomEngine>::GausACR(Double_t mean, Double_t sigma)
+{
+ const Double_t kC1 = 1.448242853;
+ const Double_t kC2 = 3.307147487;
+ const Double_t kC3 = 1.46754004;
+ const Double_t kD1 = 1.036467755;
+ const Double_t kD2 = 5.295844968;
+ const Double_t kD3 = 3.631288474;
+ const Double_t kHm = 0.483941449;
+ const Double_t kZm = 0.107981933;
+ const Double_t kHp = 4.132731354;
+ const Double_t kZp = 18.52161694;
+ const Double_t kPhln = 0.4515827053;
+ const Double_t kHm1 = 0.516058551;
+ const Double_t kHp1 = 3.132731354;
+ const Double_t kHzm = 0.375959516;
+ const Double_t kHzmp = 0.591923442;
+ /*zhm 0.967882898*/
+
+ const Double_t kAs = 0.8853395638;
+ const Double_t kBs = 0.2452635696;
+ const Double_t kCs = 0.2770276848;
+ const Double_t kB = 0.5029324303;
+ const Double_t kX0 = 0.4571828819;
+ const Double_t kYm = 0.187308492 ;
+ const Double_t kS = 0.7270572718 ;
+ const Double_t kT = 0.03895759111;
+
+ Double_t result;
+ Double_t rn,x,y,z;
+
+ do {
+ y = Rndm();
+
+ if (y>kHm1) {
+ result = kHp*y-kHp1; break; }
+
+ else if (y<kZm) {
+ rn = kZp*y-1;
+ result = (rn>0) ? (1+rn) : (-1+rn);
+ break;
+ }
+
+ else if (y<kHm) {
+ rn = Rndm();
+ rn = rn-1+rn;
+ z = (rn>0) ? 2-rn : -2-rn;
+ if ((kC1-y)*(kC3+std::abs(z))<kC2) {
+ result = z; break; }
+ else {
+ x = rn*rn;
+ if ((y+kD1)*(kD3+x)<kD2) {
+ result = rn; break; }
+ else if (kHzmp-y<exp(-(z*z+kPhln)/2)) {
+ result = z; break; }
+ else if (y+kHzm<exp(-(x+kPhln)/2)) {
+ result = rn; break; }
+ }
+ }
+
+ while (1) {
+ x = Rndm();
+ y = kYm * Rndm();
+ z = kX0 - kS*x - y;
+ if (z>0)
+ rn = 2+y/x;
+ else {
+ x = 1-x;
+ y = kYm-y;
+ rn = -(2+y/x);
+ }
+ if ((y-kAs+x)*(kCs+x)+kBs<0) {
+ result = rn; break; }
+ else if (y<x+kT)
+ if (rn*rn<4*(kB-log(x))) {
+ result = rn; break; }
+ }
+ } while(0);
+
+ return mean + sigma * result;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Generate a random number following a Landau distribution
+/// with location parameter mu and scale parameter sigma:
+/// Landau( (x-mu)/sigma )
+/// Note that mu is not the mpv(most probable value) of the Landa distribution
+/// and sigma is not the standard deviation of the distribution which is not defined.
+/// For mu =0 and sigma=1, the mpv = -0.22278
+///
+/// The Landau random number generation is implemented using the
+/// function landau_quantile(x,sigma), which provides
+/// the inverse of the landau cumulative distribution.
+/// landau_quantile has been converted from CERNLIB ranlan(G110).
+
+//Double_t RandomFunctionsImpl<TRandomEngine>::Landau(Double_t mu, Double_t sigma)
+//{
+// if (sigma <= 0) return 0;
+// Double_t x = Rndm();
+// Double_t res = mu + BA_ROOT::Math::landau_quantile(x, sigma);
+// return res;
+//}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Generates a random integer N according to a Poisson law.
+/// Prob(N) = exp(-mean)*mean^N/Factorial(N)
+///
+/// Use a different procedure according to the mean value.
+/// The algorithm is the same used by CLHEP.
+/// For lower value (mean < 25) use the rejection method based on
+/// the exponential.
+/// For higher values use a rejection method comparing with a Lorentzian
+/// distribution, as suggested by several authors.
+/// This routine since is returning 32 bits integer will not work for values
+/// larger than 2*10**9.
+/// One should then use the Trandom::PoissonD for such large values.
+
+//Int_t RandomFunctionsImpl<TRandomEngine>::Poisson(Double_t mean)
+//{
+// Int_t n;
+// if (mean <= 0) return 0;
+// if (mean < 25) {
+// Double_t expmean = std::exp(-mean);
+// Double_t pir = 1;
+// n = -1;
+// while(1) {
+// n++;
+// pir *= Rndm();
+// if (pir <= expmean) break;
+// }
+// return n;
+// }
+// // for large value we use inversion method
+// else if (mean < 1E9) {
+// Double_t em, t, y;
+// Double_t sq, alxm, g;
+// Double_t pi = M_PI;
+
+// sq = std::sqrt(2.0*mean);
+// alxm = std::log(mean);
+// g = mean*alxm - TMath::LnGamma(mean + 1.0);
+
+// do {
+// do {
+// y = std::tan(pi*Rndm());
+// em = sq*y + mean;
+// } while( em < 0.0 );
+
+// em = std::floor(em);
+// t = 0.9*(1.0 + y*y)* std::exp(em*alxm - TMath::LnGamma(em + 1.0) - g);
+// } while( Rndm() > t );
+
+// return static_cast<Int_t> (em);
+
+// }
+// else {
+// // use Gaussian approximation vor very large values
+// n = Int_t(Gaus(0,1)*std::sqrt(mean) + mean +0.5);
+// return n;
+// }
+//}
+
+//////////////////////////////////////////////////////////////////////////////////
+///// Generates a random number according to a Poisson law.
+///// Prob(N) = exp(-mean)*mean^N/Factorial(N)
+/////
+///// This function is a variant of RandomFunctionsImpl<TRandomEngine>::Poisson returning a double
+///// instead of an integer.
+
+//Double_t RandomFunctionsImpl<TRandomEngine>::PoissonD(Double_t mean)
+//{
+// Int_t n;
+// if (mean <= 0) return 0;
+// if (mean < 25) {
+// Double_t expmean = std::exp(-mean);
+// Double_t pir = 1;
+// n = -1;
+// while(1) {
+// n++;
+// pir *= Rndm();
+// if (pir <= expmean) break;
+// }
+// return static_cast<Double_t>(n);
+// }
+// // for large value we use inversion method
+// else if (mean < 1E9) {
+// Double_t em, t, y;
+// Double_t sq, alxm, g;
+// Double_t pi = M_PI;
+
+// sq = std::sqrt(2.0*mean);
+// alxm = std::log(mean);
+// g = mean*alxm - TMath::LnGamma(mean + 1.0);
+
+// do {
+// do {
+// y = std::tan(pi*Rndm());
+// em = sq*y + mean;
+// } while( em < 0.0 );
+
+// em = std::floor(em);
+// t = 0.9*(1.0 + y*y)* std::exp(em*alxm - TMath::LnGamma(em + 1.0) - g);
+// } while( Rndm() > t );
+
+// return em;
+
+// } else {
+// // use Gaussian approximation vor very large values
+// return Gaus(0,1)*std::sqrt(mean) + mean +0.5;
+// }
+//}
+
+
+////////////////////////////////////////////////////////////////////////////////
+/// Return 2 numbers distributed following a gaussian with mean=0 and sigma=1.
+
+void RandomFunctionsImpl<TRandomEngine>::Rannor(Double_t &a, Double_t &b)
+{
+ Double_t r, x, y, z;
+
+ y = Rndm();
+ z = Rndm();
+ x = z * 6.28318530717958623;
+ r = std::sqrt(-2*std::log(y));
+ a = r * std::sin(x);
+ b = r * std::cos(x);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Generates random vectors, uniformly distributed over the surface
+/// of a sphere of given radius.
+/// Input : r = sphere radius
+/// Output: x,y,z a random 3-d vector of length r
+/// Method: (based on algorithm suggested by Knuth and attributed to Robert E Knop)
+/// which uses less random numbers than the CERNLIB RN23DIM algorithm
+
+void RandomFunctionsImpl<TRandomEngine>::Sphere(Double_t &x, Double_t &y, Double_t &z, Double_t r)
+{
+ Double_t a=0,b=0,r2=1;
+ while (r2 > 0.25) {
+ a = Rndm() - 0.5;
+ b = Rndm() - 0.5;
+ r2 = a*a + b*b;
+ }
+ z = r* ( -1. + 8.0 * r2 );
+
+ Double_t scale = 8.0 * r * std::sqrt(0.25 - r2);
+ x = a*scale;
+ y = b*scale;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Returns a uniform deviate on the interval (0, x1).
+
+Double_t RandomFunctionsImpl<TRandomEngine>::Uniform(Double_t x1)
+{
+ Double_t ans = Rndm();
+ return x1*ans;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+/// Returns a uniform deviate on the interval (x1, x2).
+
+double RandomFunctionsImpl<TRandomEngine>::Uniform(double a, double b) {
+ return (b-a) * Rndm() + a;
+}
+
+
+ } // namespace Math
+} // namespace ROOT
diff --git a/ThirdParty/RootMinimizers/src/Minuit2/VariableMetricBuilder.cxx b/ThirdParty/RootMinimizers/src/Minuit2/VariableMetricBuilder.cxx
index e79683ff89c..540324e0e6c 100644
--- a/ThirdParty/RootMinimizers/src/Minuit2/VariableMetricBuilder.cxx
+++ b/ThirdParty/RootMinimizers/src/Minuit2/VariableMetricBuilder.cxx
@@ -39,14 +39,14 @@ namespace BA_ROOT {
double inner_product(const LAVector&, const LAVector&);
-void VariableMetricBuilder::AddResult( std::vector<MinimumState>& result, const MinimumState & state, bool store) const {
- if (!store) store = StorageLevel();
- store |= (result.size() == 0);
- if (store)
- result.push_back(state);
- else {
- result.back() = state;
- }
+void VariableMetricBuilder::AddResult( std::vector<MinimumState>& result, const MinimumState & state) const {
+ // // if (!store) store = StorageLevel();
+ // // store |= (result.size() == 0);
+ // if (store)
+ result.push_back(state);
+ // else {
+ // result.back() = state;
+ // }
if (TraceIter() ) TraceIteration(result.size()-1, result.back() );
else {
if (PrintLevel() > 1) {
@@ -68,6 +68,7 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
int printLevel = PrintLevel();
+
#ifdef DEBUG
std::cout<<"VariableMetricBuilder convergence when edm < "<<edmval<<std::endl;
#endif
@@ -103,7 +104,7 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
MnPrint::PrintState(std::cout, seed.State(), "VariableMetric: Initial state ");
}
- AddResult( result, seed.State(), true);
+ AddResult( result, seed.State());
// try first with a maxfxn = 80% of maxfcn
@@ -157,8 +158,7 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
if (printLevel > 1) {
MnPrint::PrintState(std::cout, st, "VariableMetric: After Hessian ");
}
- AddResult( result, st, true);
-
+ AddResult( result, st);
// check new edm
edm = st.Edm();
@@ -263,10 +263,12 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
// keep also prevStep
MnAlgebraicVector prevStep(initialState.Gradient().Vec().size());
+ MinimumState s0 = result.back();
+ assert(s0.IsValid() );
+
do {
- // const MinimumState& s0 = result.back();
- MinimumState s0 = result.back();
+ //MinimumState s0 = result.back();
step = -1.*s0.Error().InvHessian()*s0.Gradient().Vec();
@@ -299,20 +301,22 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
MN_INFO_VAL(gdel);
#endif
MnPosDef psdf;
- s0 = psdf(s0, prec);
- step = -1.*s0.Error().InvHessian()*s0.Gradient().Vec();
+ MinimumState s0new = psdf(s0, prec);
+ step = -1.*s0new.Error().InvHessian()*s0new.Gradient().Vec();
// #ifdef DEBUG
// std::cout << "After MnPosdef - Error " << s0.Error().InvHessian() << " Gradient " << s0.Gradient().Vec() << " step " << step << std::endl;
// #endif
- gdel = inner_product(step, s0.Gradient().Grad());
+ gdel = inner_product(step, s0new.Gradient().Grad());
#ifdef WARNINGMSG
MN_INFO_VAL(gdel);
#endif
if(gdel > 0.) {
- AddResult(result, s0, result.size()<=1);
+ AddResult(result, s0new);
+
return FunctionMinimum(seed, result, fcn.Up());
}
}
+
MnParabolaPoint pp = lsearch(fcn, s0.Parameters(), step, gdel, prec);
// <= needed for case 0 <= 0
@@ -321,10 +325,12 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
MN_INFO_MSG("VariableMetricBuilder: no improvement in line search");
#endif
// no improvement exit (is it really needed LM ? in vers. 1.22 tried alternative )
- // add new state where only fcn changes
- AddResult(result, MinimumState(s0.Parameters(), s0.Error(), s0.Gradient(), s0.Edm(), fcn.NumOfCalls()),
- result.size()<=1);
-
+ // add new state when only fcn changes
+ if (result.size() <= 1 )
+ AddResult(result, MinimumState(s0.Parameters(), s0.Error(), s0.Gradient(), s0.Edm(), fcn.NumOfCalls()));
+ else
+ // no need to re-store the state
+ AddResult(result, MinimumState(pp.Y(), s0.Edm(), fcn.NumOfCalls()));
break;
@@ -352,13 +358,14 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
MN_INFO_MSG("VariableMetricBuilder: matrix not pos.def. : edm is < 0. Make pos def...");
#endif
MnPosDef psdf;
- s0 = psdf(s0, prec);
- edm = Estimator().Estimate(g, s0.Error());
+ MinimumState s0new = psdf(s0, prec);
+ edm = Estimator().Estimate(g, s0new.Error());
if(edm < 0.) {
#ifdef WARNINGMSG
MN_INFO_MSG("VariableMetricBuilder: matrix still not pos.def. : exit iterations ");
#endif
- result.push_back(s0);
+ AddResult(result, s0new);
+
return FunctionMinimum(seed, result, fcn.Up());
}
}
@@ -373,7 +380,13 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
<< " edm = " << edm << std::endl << std::endl;
#endif
- AddResult(result, MinimumState(p, e, g, edm, fcn.NumOfCalls()), result.size() <= 1);
+ // update the state
+ s0 = MinimumState(p, e, g, edm, fcn.NumOfCalls());
+ if (StorageLevel() || result.size() <= 1)
+ AddResult(result, s0);
+ else
+ // use a reduced state for not-final iterations
+ AddResult(result, MinimumState(p.Fval(), edm, fcn.NumOfCalls()));
// correct edm
edm *= (1. + 3.*e.Dcovar());
@@ -386,6 +399,10 @@ FunctionMinimum VariableMetricBuilder::Minimum(const MnFcn& fcn, const GradientC
} while(edm > edmval && fcn.NumOfCalls() < maxfcn); // end of iteration loop
+ // save last result in case of no complete final states
+ if ( ! result.back().IsValid() )
+ result.back() = s0;
+
if(fcn.NumOfCalls() >= maxfcn) {
#ifdef WARNINGMSG
diff --git a/cmake/modules/SearchInstalledSoftware.cmake b/cmake/modules/SearchInstalledSoftware.cmake
index 5d6f9f365c1..f40b7ca453d 100644
--- a/cmake/modules/SearchInstalledSoftware.cmake
+++ b/cmake/modules/SearchInstalledSoftware.cmake
@@ -5,7 +5,10 @@
# --- math packages ---
find_package(Eigen3 REQUIRED)
find_package(FFTW REQUIRED)
-find_package(GSL 1 EXACT REQUIRED) # revert this when issue 1404 is resolved
+#find_package(GSL 1 EXACT REQUIRED) # revert this when issue 1404 is resolved
+find_package(GSL REQUIRED) # revert this when issue 1404 is resolved
+message(INFO "XXX ${GSL_VERSION}")
+
# --- Boost ---
set(Boost_NO_BOOST_CMAKE ON) # prevent shortcut
diff --git a/dev-tools/code-tools/update_minimizers.py b/dev-tools/code-tools/update_minimizers.py
index 0d4571364b7..34692da4cfe 100644
--- a/dev-tools/code-tools/update_minimizers.py
+++ b/dev-tools/code-tools/update_minimizers.py
@@ -5,7 +5,7 @@ from ROOT to our source tree.
import os
import filecmp
-ROOT_SOURCE = "/home/pospelov/software/root/root-6.04.02.source"
+ROOT_SOURCE = "/home/pospelov/software/root/root_v6.06.02.source"
BORNAGAIN_SOURCE = "/home/pospelov/development/BornAgain/BornAgain/ThirdParty/RootMinimizers"
--
GitLab