diff --git a/cmake/modules/FindLibCmaes.cmake b/cmake/modules/FindLibCmaes.cmake new file mode 100644 index 0000000000000..00757a145447a --- /dev/null +++ b/cmake/modules/FindLibCmaes.cmake @@ -0,0 +1,53 @@ +# - Try to find libcmaes, https://github.com/beniz/libcmaes +# Once done this will define +# +# LIBCMAES_FOUND - system has libcmaes +# LIBCMES_INCLUDE_DIR - the libcmaes include directory +# LIBCMAES_LIBRARIES - Link these to use libcmaes +# LIBCMAES_DEFINITIONS - Compiler switches required for using libcmaes +# Redistribution and use is allowed according to the terms of the BSD license. +# For details see the accompanying COPYING-CMAKE-SCRIPTS file. +# + + +# Copyright (c) 2014, Emmanuel Benazera, +# +# Redistribution and use is allowed according to the terms of the BSD license. +# For details see the accompanying COPYING-CMAKE-SCRIPTS file. + +if ( LIBCMAES_INCLUDE_DIR AND LIBCMAES_LIBRARIES ) + # in cache already + SET(Libcmaes_FIND_QUIETLY TRUE) +endif ( LIBCMAES_INCLUDE_DIR AND LIBCMAES_LIBRARIES ) + +# use pkg-config to get the directories and then use these values +# in the FIND_PATH() and FIND_LIBRARY() calls +if( NOT WIN32 ) + find_package(PkgConfig) + pkg_check_modules(PC_LIBCMAES QUIET libcmaes) + + if( PKG_CONFIG_FOUND ) + set(LIBCMAES_INCLUDE_DIR ${PC_LIBCMAES_INCLUDE_DIRS}) + set(LIBCMAES_LIBRARIES ${PC_LIBCMAES_LIBRARY_DIRS}) + set(LIBCMAES_DEFINITIONS ${PC_LIBCMAES_CFLAGS_OTHER}) + endif( PKG_CONFIG_FOUND ) + +endif( NOT WIN32 ) + +find_path(LIBCMAES_INCLUDE_DIR NAMES cmaes.h + PATHS + ${PC_LIBCMAES_INCLUDEDIR} + ${PC_LIBCMAES_INCLUDE_DIRS} +) + +find_library(LIBCMAES_LIBRARIES NAMES libcmaes + PATHS + ${PC_LIBCMAES_LIBDIR} + ${PC_LIBCMAES_LIBRARY_DIRS} +) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(Libcmaes DEFAULT_MSG LIBCMAES_INCLUDE_DIR LIBCMAES_LIBRARIES ) + +# show the LIBCMAES_INCLUDE_DIR and LIBCMAES_LIBRARIES variables only in the advanced view +mark_as_advanced(LIBCMAES_INCLUDE_DIR LIBCMAES_LIBRARIES ) diff --git a/cmake/modules/RootBuildOptions.cmake b/cmake/modules/RootBuildOptions.cmake index 96ee1aae997c5..948c80f9ad378 100644 --- a/cmake/modules/RootBuildOptions.cmake +++ b/cmake/modules/RootBuildOptions.cmake @@ -132,6 +132,7 @@ ROOT_BUILD_OPTION(http ON "Enable suppport for HTTP server") ROOT_BUILD_OPTION(fcgi OFF "Enable FastCGI suppport in HTTP server") ROOT_BUILD_OPTION(imt ON "Enable support for implicit multi-threading via Intel® Thread Bulding Blocks (TBB)") ROOT_BUILD_OPTION(jemalloc OFF "Use jemalloc memory allocator") +ROOT_BUILD_OPTION(libcmaes OFF "Build the libcmaes minimizer library") ROOT_BUILD_OPTION(libcxx OFF "Build using libc++") ROOT_BUILD_OPTION(macos_native OFF "Disable looking for libraries, includes and binaries in locations other than a native installation (MacOS only)") ROOT_BUILD_OPTION(mathmore ON "Build libMathMore extended math library (requires GSL)") diff --git a/cmake/modules/RootConfiguration.cmake b/cmake/modules/RootConfiguration.cmake index d3e346460d833..d986619cf0c43 100644 --- a/cmake/modules/RootConfiguration.cmake +++ b/cmake/modules/RootConfiguration.cmake @@ -317,6 +317,7 @@ set(gslflags) set(shadowpw ${value${shadowpw}}) set(buildmathmore ${value${mathmore}}) +set(buildcmaes ${value${libcmaes}}) set(buildroofit ${value${roofit}}) set(buildminuit2 ${value${minuit2}}) set(buildunuran ${value${unuran}}) diff --git a/cmake/modules/SearchInstalledSoftware.cmake b/cmake/modules/SearchInstalledSoftware.cmake index d157d1dbdd83d..672d6fd94d28f 100644 --- a/cmake/modules/SearchInstalledSoftware.cmake +++ b/cmake/modules/SearchInstalledSoftware.cmake @@ -729,6 +729,39 @@ if(sqlite) endif() endif() +#---Check for Libcmaes------------------------------------------------------------------- +if(libcmaes) + message(STATUS "Looking for libcmaes") + find_package(LibCmaes) + if(NOT LIBCMAES_FOUND) + message(STATUS "Fetching and compiling libcmaes") + ExternalProject_Add( + eigen3 + PREFIX eigen3 + INSTALL_DIR ${CMAKE_BINARY_DIR} + URL http://bitbucket.org/eigen/eigen/get/3.2.2.tar.gz + CMAKE_ARGS -DCMAKE_INSTALL_PREFIX= + ) + ExternalProject_Add( + cma + DEPENDS eigen3 + PREFIX cma + INSTALL_DIR ${CMAKE_BINARY_DIR} + URL https://github.com/beniz/libcmaes/archive/master.tar.gz + CONFIGURE_COMMAND ./autogen.sh && ./configure --prefix= --with-eigen3-include=${CMAKE_BINARY_DIR}/include/eigen3 --enable-onlylib --disable-surrog + BUILD_IN_SOURCE 1 + ) + set(LIBCMAES_INCLUDE_DIR ${CMAKE_BINARY_DIR}/include/libcmaes ${CMAKE_BINARY_DIR}/include/eigen3) + set(LIBCMAES_LIBRARIES ${CMAKE_BINARY_DIR}/lib/) + #if(fail-on-missing) + # message(FATAL_ERROR "libcmaes not found") + #else() + # message(STATUS "libcmaes not found. Switching off libcmaes option") + # set(libcmaes OFF CACHE BOOL "" FORCE) + #endif() + endif() +endif() + #---Check for Pythia6------------------------------------------------------------------- if(pythia6) message(STATUS "Looking for Pythia6") diff --git a/config/Makefile.in b/config/Makefile.in index ffc31251a2425..52de7df6d978e 100644 --- a/config/Makefile.in +++ b/config/Makefile.in @@ -316,6 +316,11 @@ BUILDROOFIT := @buildroofit@ BUILDMINUIT2 := @buildminuit2@ +BUILDCMAES := @buildcmaes@ +CMAESINCDIR := $(filter-out -I/usr/include, @cmaesincdir@) +CMAESLIBDIR := @cmaeslibdir@ +CMAESLIB := -lcmaes + BUILDUNURAN := @buildunuran@ BUILDVC := @buildvc@ diff --git a/core/base/src/TPluginManager.cxx b/core/base/src/TPluginManager.cxx index 1f654c62d70ea..0cfe7512eab81 100644 --- a/core/base/src/TPluginManager.cxx +++ b/core/base/src/TPluginManager.cxx @@ -375,7 +375,7 @@ void TPluginManager::LoadHandlersFromEnv(TEnv *env) TString ctor = strtok(0, ";\""); if (!ctor.Contains("(")) ctor = strtok(0, ";\""); - AddHandler(s, regexp, clss, plugin, ctor, "TEnv"); + AddHandler(s, regexp, clss, plugin, ctor, "TEnv"); cnt++; } delete [] v; @@ -587,8 +587,8 @@ TPluginHandler *TPluginManager::FindHandler(const char *base, const char *uri) while ((h = (TPluginHandler*) next())) { if (h->CanHandle(base, uri)) { - if (gDebug > 0) - Info("FindHandler", "found plugin for %s", h->GetClass()); + /*if (gDebug > 0) + Info("FindHandler", "found plugin for %s", h->GetClass());*/ return h; } } diff --git a/etc/plugins/ROOT@@Math@@Minimizer/P010_Minuit2Minimizer.C b/etc/plugins/ROOT@@Math@@Minimizer/P010_Minuit2Minimizer.C index 781b5f7e4c03f..2575db17ce019 100644 --- a/etc/plugins/ROOT@@Math@@Minimizer/P010_Minuit2Minimizer.C +++ b/etc/plugins/ROOT@@Math@@Minimizer/P010_Minuit2Minimizer.C @@ -1,5 +1,5 @@ void P010_Minuit2Minimizer() { - gPluginMgr->AddHandler("ROOT::Math::Minimizer", "Minuit2", "ROOT::Minuit2::Minuit2Minimizer", - "Minuit2", "Minuit2Minimizer(const char *)"); + gPluginMgr->AddHandler("ROOT::Math::Minimizer", "Minuit2", "ROOT::Minuit2::Minuit2Minimizer", + "Minuit2", "Minuit2Minimizer(const char *)"); } diff --git a/etc/plugins/ROOT@@Math@@Minimizer/P090_TCMAESMinimizer.C b/etc/plugins/ROOT@@Math@@Minimizer/P090_TCMAESMinimizer.C new file mode 100644 index 0000000000000..161ffdc5a3b8e --- /dev/null +++ b/etc/plugins/ROOT@@Math@@Minimizer/P090_TCMAESMinimizer.C @@ -0,0 +1,5 @@ +void P090_TCMAESMinimizer() +{ + gPluginMgr->AddHandler("ROOT::Math::Minimizer", "cmaes", "ROOT::cmaes::TCMAESMinimizer", + "cmaes_root", "TCMAESMinimizer(const char *)"); +} diff --git a/graf2d/asimage/src/libAfterImage/asfont.c b/graf2d/asimage/src/libAfterImage/asfont.c index 415bf6147686b..c069ddae2ab33 100644 --- a/graf2d/asimage/src/libAfterImage/asfont.c +++ b/graf2d/asimage/src/libAfterImage/asfont.c @@ -63,11 +63,11 @@ # include # include FT_FREETYPE_H # endif -# ifdef HAVE_FREETYPE_FREETYPE +/*# ifdef HAVE_FREETYPE_FREETYPE # include -# else +# else*/ # include -# endif +//# endif # if (FREETYPE_MAJOR == 2) && ((FREETYPE_MINOR == 0) || ((FREETYPE_MINOR == 1) && (FREETYPE_PATCH < 3))) # define FT_KERNING_DEFAULT ft_kerning_default # endif diff --git a/hist/hist/src/TF1.cxx b/hist/hist/src/TF1.cxx index b46690e3eb74b..5d6a44317c0f2 100644 --- a/hist/hist/src/TF1.cxx +++ b/hist/hist/src/TF1.cxx @@ -1766,7 +1766,6 @@ Double_t TF1::GetMinMaxNDim(Double_t *x , bool findmax, Double_t epsilon, Int_t stepSize = (rmax[i] - rmin[i]) / 100; else if (std::abs(x[i]) > 1.) stepSize = 0.1 * x[i]; - // set variable names if (ndim <= 3) { if (i == 0) { diff --git a/math/CMakeLists.txt b/math/CMakeLists.txt index e3c6e62b4ecf3..60e65f776699a 100644 --- a/math/CMakeLists.txt +++ b/math/CMakeLists.txt @@ -13,6 +13,9 @@ add_subdirectory(minuit) if(minuit2) add_subdirectory(minuit2) endif() +if (libcmaes) + add_subdirectory(cmaes) +endif() add_subdirectory(fumili) add_subdirectory(physics) if(mlp) diff --git a/math/cmaes/CMakeLists.txt b/math/cmaes/CMakeLists.txt new file mode 100644 index 0000000000000..d3d3a2331270e --- /dev/null +++ b/math/cmaes/CMakeLists.txt @@ -0,0 +1,13 @@ +############################################################################ +# CMakeLists.txt file for building ROOT math/cmaes package +############################################################################ + +add_definitions(-DWARNINGMSG -DUSE_ROOT_ERROR) + +include_directories(${LIBCMAES_INCLUDE_DIR}) +link_directories(${LIBCMAES_LIBRARIES}) + +ROOT_GENERATE_DICTIONARY(G__cmaes *.h MODULE cmaes LINKDEF LinkDef.h) + +ROOT_LINKER_LIBRARY(cmaes_root *.cxx G__cmaes.cxx DEPENDENCIES MathCore Hist LIBRARIES cmaes) +ROOT_INSTALL_HEADERS() diff --git a/math/cmaes/Module.mk b/math/cmaes/Module.mk new file mode 100644 index 0000000000000..e5c19d24a402d --- /dev/null +++ b/math/cmaes/Module.mk @@ -0,0 +1,140 @@ +# Module.mk for cmaes module. + +MODNAME := cmaes +MODDIR := $(ROOT_SRCDIR)/math/$(MODNAME) +MODDIRS := $(MODDIR)/src +MODDIRI := $(MODDIR)/inc + +CMAESDIR := $(MODDIR) +CMAESDIRS := $(CMAESDIR)/src +CMAESDIRI := $(CMAESDIR)/inc +#CMAESDIRT :=$(call stripsrc,$(CMAESDIR)/test) + +CMAESBASEVERS := cmaes-1_0_0 +CMAESBASESRCS := $(MODDIRS)/$(CMAESBASEVERS).tar.gz +CMAESBASEDIRS := $(MODDIRS)/$(CMAESBASEVERS) +CMAESBASEDIRI := -I$(MODDIRS)/$(CMAESBASEVERS) +CMAESBASEETAG := $(MODDIRS)/headers.d + +##### liblcg_cmaes ##### +ifeq ($(PLATFORM),win32) +CMAESBASELIBA := $(CMAESBASEDIRS)/libcmaesbase.lib +CMAESBASELIB := $(LPATH)/libcmaesbase.lib +ifeq (debug,$(findstring debug,$(ROOTBUILD))) +CMAESBASEBLD = "DEBUG=1" +else +CMAESBASEBLD = "" +endif +else +CMAESBASELIBA := $(CMAESBASEDIRS)/src/.libs/liblcg_cmaes.a +CMAESBASELIB := $(LPATH)/libcmaesbase.a +endif +CMAESBASEDEP := $(CMAESBASELIB) + +##### libcmaes ##### +CMAESL := $(MODDIRI)/LinkDef.h +CMAESDS := $(call stripsrc,$(MODDIRS)/G__cmaes.cxx) +CMAESDO := $(CMAESDS:.cxx=.o) +CMAESDH := $(CMAESDS:.cxx=.h) + +CMAESAH := $(filter-out $(MODDIRI)/LinkDef%,$(wildcard $(MODDIRI)/*.h)) +CMAESBH := $(filter-out $(MODDIRI)/cmaes/LinkDef%,$(wildcard $(MODDIRI)/cmaes/*.h)) +CMAESH := $(CMAESAH) $(CMAESBH) +CMAESS := $(filter-out $(MODDIRS)/G__%,$(wildcard $(MODDIRS)/*.cxx)) +CMAESO := $(call stripsrc,$(CMAESS:.cxx=.o)) + +CMAESDEP := $(CMAESO:.o=.d) $(CMAESDO:.o=.d) + +CMAESLIB := $(LPATH)/libcmaes_root.$(SOEXT) +CMAESMAP := $(CMAESLIB:.$(SOEXT)=.rootmap) + +# use this compiler option if want to optimize object allocation in cmaes +# NOTE: using this option one loses the thread safety. +# It is worth to use only for minimization of cheap (non CPU intensive) functions +#CXXFLAGS += -DMN_USE_STACK_ALLOC + +# used in the main Makefile +ALLHDRS += $(patsubst $(MODDIRI)/%.h,include/%.h,$(CMAESH)) +ALLLIBS += $(CMAESLIB) +ALLMAPS += $(CMAESMAP) + +# include all dependency files +INCLUDEFILES += $(CMAESDEP) + +##### local rules ##### +.PHONY: all-$(MODNAME) clean-$(MODNAME) distclean-$(MODNAME) \ + test-$(MODNAME) + +include/cmaes/%.h: $(CMAESDIRI)/cmaes/%.h + @(if [ ! -d "include/cmaes" ]; then \ + mkdir -p include/cmaes; \ + fi) + cp $< $@ + +include/%.h: $(CMAESDIRI)/%.h + cp $< $@ + +$(CMAESLIB): $(CMAESO) $(CMAESDO) $(ORDER_) $(MAINLIBS) $(CMAESLIBDEP) + @$(MAKELIB) $(PLATFORM) $(LD) "$(LDFLAGS)" \ + "$(SOFLAGS)" libcmaes_root.$(SOEXT) $@ \ + "$(CMAESO) $(CMAESDO)" "-lcmaes -lMathCore" "$(CMAESLIBEXTRA)" + +$(call pcmrule,CMAES) + $(noop) + +$(CMAESDS): $(CMAESH) $(CMAESL) $(ROOTCLINGEXE) $(call pcmdep,CMAES) + $(MAKEDIR) + @echo "Generating dictionary $@..." + $(ROOTCLINGSTAGE2) -f $@ $(call dictModule,CMAES) -I$(CMAESINCDIR) -I/usr/include/eigen3 -c $(CMAESH) $(CMAESL) + +$(CMAESMAP): $(CMAESH) $(CMAESL) $(ROOTCLINGEXE) $(call pcmdep,CMAES) + $(MAKEDIR) + @echo "Generating rootmap $@..." + $(ROOTCLINGSTAGE2) -r $(CMAESDS) $(call dictModule,CMAES) -c $(CMAESH) $(CMAESL) + +all-$(MODNAME): $(CMAESLIB) + +#test-$(MODNAME): all-$(MODNAME) +#ifneq ($(ROOT_OBJDIR),$(ROOT_SRCDIR)) +# @$(INSTALL) $(CMAESDIR)/test $(CMAESDIRT) +#endif + @cd $(CMAESDIRT) && $(MAKE) ROOTCONFIG=../../../bin/root-config + +clean-$(MODNAME): + @rm -f $(CMAESO) $(CMAESDO) + +clean:: clean-$(MODNAME) + +distclean-$(MODNAME): clean-$(MODNAME) + @rm -f $(CMAESDEP) $(CMAESDS) $(CMAESDH) $(CMAESLIB) \ + $(CMAESMAP) + @rm -rf include/cmaes +ifneq ($(ROOT_OBJDIR),$(ROOT_SRCDIR)) + @rm -rf $(CMAESDIRT) +else + @cd $(CMAESDIRT) && $(MAKE) distclean ROOTCONFIG=../../../bin/root-config +endif + +distclean:: distclean-$(MODNAME) + +##### extra rules ###### +$(CMAESO): CXXFLAGS += -DWARNINGMSG -DUSE_ROOT_ERROR +$(CMAESDO): CXXFLAGS += -DWARNINGMSG -DUSE_ROOT_ERROR +#for thread -safet +#$(CMAESO): CXXFLAGS += -DCMAES_THREAD_SAFE +# for openMP +#ifneq ($(USE_PARALLEL_CMAES),) +#ifneq ($(USE_OPENMP),) +#$(CMAESO): CXXFLAGS += -DCMAES_THREAD_SAFE -DCMAES_PARALLEL_OPENMP +#math/cmaes/src/Numerical2PGradientCalculator.o: +$(CMAESO):CXXFLAGS += -D_GLIBCXX_PARALLEL -fopenmp +$(CMAESDO):CXXFLAGS += -D_GLIBCXX_PARALLEL -fopenmp +$(CMAESLIB):LDFLAGS += -fopenmp +#endif + +$(CMAESO):CXXFLAGS += -I$(CMAESINCDIR) -I/usr/include/eigen3 +$(CMAESDO):CXXFLAGS += -I$(CMAESINCDIR) -I/usr/include/eigen3 +$(CMAESDS):CXXFLAGS += -I$(CMAESINCDIR) -I/usr/include/eigen3 +$(CMAESLIB):LDFLAGS += $(CMAESLIBDIR) +# Optimize dictionary with stl containers. +$(CMAESDO): NOOPT = $(OPT) diff --git a/math/cmaes/README.md b/math/cmaes/README.md new file mode 100644 index 0000000000000..d0173c1c4fcf2 --- /dev/null +++ b/math/cmaes/README.md @@ -0,0 +1,128 @@ +### This is support for black-box optimization with CMA-ES from within ROOT (http://root.cern.ch/). + +This work was supported by the ANR-2010-COSI-002 grant of the French NationalA Research Agency. + +** This work is improved when needed, see status here: https://github.com/beniz/libcmaes/issues/13 ** + +** This work may be about to be merged into ROOT repository, follow integration here: +https://github.com/root-mirror/root/pull/40 ** + +libcmaes can be used from CERN's ROOT6 as a seamless replacement or addition to Minuit2 optimizer. It is designed to be used from ROOT6 **exactly** as Minuit2 is used, so code using Minuit2 should be easily run against CMA-ES. + +Below are instructions for testing it out. + +### Building ROOT6 and libcmaes +As for now, the only way to use libcmaes is from ROOT6, using the following special repository, and compiling it from sources (1): https://github.com/beniz/root + +* get ROOT6 from https://github.com/beniz/root/tree/cmaes4root_master, configure & compile it (this will take a while) (2): +```` +git clone https://github.com/beniz/root.git +cd root +``` +It is recommended to build with cmake, see below. + +#### Build with autoconf (configure) #### +For this, you need to have libcmaes installed already, see https://github.com/beniz/libcmaes/wiki +``` +./configure --enable-minuit2 --enable-roofit --enable-python --with-cmaes-incdir=/home/yourusername/include/libcmaes --with-cmaes-libdir=/home/yourusername/lib +make +```` +use make -jx where x is the number of cores on your system in order to minimize the building time. + +#### Build with cmake #### +```` +mkdir mybuild +cd mybuild +cmake ../ -Dall=on -Dtesting=on -Dlibcmaes=on +make +```` +use make -jx where x is the number of cores on your system in order to minimize the building time. + +### Running an example with CMA-ES +To run the basic fitting of a Gaussian, originally taken from Minuit2's tutorial files, do: +```` +root +.L tutorials/fit/cmaesGausFit.C++g +cmaesGausFit() +```` +You should see a plot similar to +![cmaes_gaus_fit_root_errors](https://cloud.githubusercontent.com/assets/3530657/2890890/4d96ae1c-d52d-11e3-9610-f24790b23e98.png) + +To quick test competitiveness against Minuit2: +```` +root +.L tutorials/fit/cmaesFitBench.C + cmaesFitBench() +```` +You should witness a plot similar to +![](http://juban.free.fr/stuff/libcmaes/cmaes_minuit2_competitive.png) + +### Running a benchmark comparison of CMA-ES and Minuit2 + +To run the current benchmark and visualize results, take the following steps: +```` +root +.L tutorials/fit/cmaesFullBench.C +run_experiments(10) +python math/cmaes/test/cmaesFullBench.py +```` + +This should show a series of histograms comparing results from both optimizers on a selection of problems. + +### Options to the CMA-ES minimizers within ROOT +There's built-in control for several hyper-parameters and options of CMA-ES: +* several flavors of the algorithm are available, and can be choosen at creation of the Minimizer object: +```` +TVirtualFitter::SetDefaultFitter(``acmaes''); +```` +or +```` +ROOT::Fit::Fitter fitter; +fitter.Config().SetMinimizer(``cmaes'',''acmaes''); +```` +The available algorithms are: `cmaes, ipop, bipop, acmaes, aipop, abipop, sepcmaes, sepipop, sepbipop`. + +'acmaes' should be the most appropriate in most cases, and 'sepacmaes' when the number of dimensions nears a thousand. + +The options below are not required, but can be used by filling up a MinimizerOptions object beforehand: +```` +const char *fitter = "acmaes" +TVirtualFitter::SetDefaultFitter(fitter); +ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default(fitter); +opts.SetIntValue("lambda",100); +```` +Options below are not activated by default: +* 'sigma': initial step-size +* 'lambda': number of offsprings at each generation +* 'noisy': flag that updates some hyper-parameters if the objective function is noisy +* 'restarts': maximum number of restarts, only applies to ipop, bipop, aipop, abipop, sepipop and sepbipop +* 'ftarget': the objective function target that stops optimization when reached, useful when the final value is known, e.g. 0 +* 'fplot': output file in libcmaes format for later plotting of eigenvalues and state convergence, mostly for debug purposes +* 'lscaling': automatic linear scaling of parameters with auto-selection of step-size sigma, usually recommended if results are not satisfactory +* 'mt_feval': allows for parallel calls of the objective function, faster but the objective function is required to be thread-safe. + +### Using CMA-ES from RooFit +libcmaes support within ROOT extends to RooFit without effort. + +From Roofit, it is enough to set the Minimizer with +```C++ + +``` +and from PyROOT for example +```Python +RooFit.Minimizer("cmaes","acmaes") +``` + +For setting CMA-ES custom options, such as 'sigma', 'lambda' or 'lscaling', it is enough to set the options as explained in the non RooFit case: +```C++ +ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default(fitter); +opts.SetIntValue("lambda",100); +``` +and from PyRoot +```Python +opt = ROOT.Math.MinimizerOptions.Default("cmaes") +opt.SetIntValue("lambda",100) +``` + +(1) more convenient ways will be provided. +(2) we recommend building support for both Minuit2 (i.e. for comparison to CMA-ES) and debug. \ No newline at end of file diff --git a/math/cmaes/inc/CMAESMinimizer.h b/math/cmaes/inc/CMAESMinimizer.h new file mode 100644 index 0000000000000..e0d411aa3bbd6 --- /dev/null +++ b/math/cmaes/inc/CMAESMinimizer.h @@ -0,0 +1,309 @@ +/** + * CMA-ES, Covariance Matrix Evolution Strategy + * Copyright (c) 2014 INRIA + * Author: Emmanuel Benazera + * + * This file is part of libcmaes. + * + * libcmaes is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * libcmaes 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with libcmaes. If not, see . + */ + +#ifndef ROOT_CMAES_CMAESMinimizer +#define ROOT_CMAES_CMAESMinimizer + +#include "Math/Minimizer.h" + +#include "Math/IFunctionfwd.h" + +#include "libcmaes/cmaes.h" + +using namespace libcmaes; + +namespace ROOT +{ + namespace cmaes + { + + class TCMAESMinimizer : public ROOT::Math::Minimizer + { + public: + + /** + * Default constructor + */ + TCMAESMinimizer(); // CMAES_DEFAULT + + /** + * Constructor with a char (used by PM) + */ + TCMAESMinimizer(const char *type); + + /** + * Destructor + */ + virtual ~TCMAESMinimizer(); + + private: + // blocking copy constructor. + TCMAESMinimizer(const TCMAESMinimizer &m); + + /** + * assignment operator. + */ + TCMAESMinimizer& operator = (const TCMAESMinimizer &rhs); + + public: + + // clear resources (parameters) for consecutives minimizations + virtual void Clear(); + + /// 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 variable + virtual bool SetVariableValue(unsigned int ivar, double val); + // set variable values + virtual bool SetVariableValues(const double * val); + /// set the step size of an already existing variable + virtual bool SetVariableStepSize(unsigned int ivar, double step ); + /// set the lower-limit of an already existing variable + virtual bool SetVariableLowerLimit(unsigned int ivar, double lower); + /// set the upper-limit of an already existing variable + virtual bool SetVariableUpperLimit(unsigned int ivar, double upper); + /// set the limits of an already existing variable + virtual bool SetVariableLimits(unsigned int ivar, double lower, double upper); + /// fix an existing variable + virtual bool FixVariable(unsigned int ivar); + /// release an existing variable + //virtual bool ReleaseVariable(unsigned int ivar); + /// query if an existing variable is fixed (i.e. considered constant in the minimization) + /// note that by default all variables are not fixed + virtual bool IsFixedVariable(unsigned int ivar) const; + /// get variable settings in a variable object (like ROOT::Fit::ParamsSettings) + virtual bool GetVariableSettings(unsigned int ivar, ROOT::Fit::ParameterSettings & varObj) const; + /// get name of variables (override if minimizer support storing of variable names) + virtual std::string VariableName(unsigned int ivar) const; + /// get index of variable given a variable given a name + /// return -1 if variable is not found + virtual int VariableIndex(const std::string & name) const; + + // sets libcmaes parameter object based on options. + template + void SetMParameters(CMAParameters &cmaparams, + const int &maxiter, const int &maxfevals, + const int &noisy, const int &nrestarts, + const double &ftarget, + const std::string &fplot, + const bool &withnumgradient, + const bool &mtfeval, + const bool &quiet, + const int &elitist, // 0: off (default), then 1, 2 or 3 + const bool &uh); + + /** + method to perform the minimization. + Return false in case the minimization did not converge. In this case a + status code different than zero is set + (retrieved by the derived method Minimizer::Status() )" + + status = 1 : Covariance was made pos defined -> not yet implemented in libcmaes + status = 2 : Hesse is invalid => Not Applicable, kept for correspondence with Minuit2 + status = 3 : Edm is above max + status = 4 : Reached call limit + status = 5 : Any other failure + */ + virtual bool Minimize(); + + /// return minimum function value + virtual double MinValue() const; + + /// return expected distance reached from the minimum + virtual double Edm() const; + + /// return pointer to X values at the minimum + virtual const double * X() const; + + /// return pointer to gradient values at the minimum + virtual const double * MinGradient() const { return 0; } // not available in CMAES + + /// 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 fDim; } + + /// number of free variables (real dimension of the problem) + /// this is <= Function().NDim() which is the total + virtual unsigned int NFree() const { return fFreeDim; } + + /// minimizer provides error and error matrix + virtual bool ProvidesError() const { return true; } + + /// return errors at the minimum + virtual const double* Errors() const; + + /** + return covariance matrix elements + if the variable is fixed or const the value is zero + The ordering of the variables is the same as in errors and parameter value. + This is different from the direct interface of Minuit2 or TMinuit where the + values were obtained only to variable parameters + */ + virtual double CovMatrix(unsigned int i, unsigned int j) const; + + + /** + Fill the passed array with the covariance matrix elements + if the variable is fixed or const the value is zero. + The array will be filled as cov[i *ndim + j] + The ordering of the variables is the same as in errors and parameter value. + This is different from the direct interface of Minuit2 or TMinuit where the + values were obtained only to variable parameters + */ + virtual bool GetCovMatrix(double * cov) const; + + /** + Fill the passed array with the Hessian matrix elements + The Hessian matrix is the matrix of the second derivatives + and is the inverse of the covariance matrix + If the variable is fixed or const the values for that variables are zero. + The array will be filled as h[i *ndim + j] + */ + //virtual bool GetHessianMatrix(double * h) const; + + + /** + return the status of the covariance matrix + status = -1 : not available (inversion failed or Hesse failed) + status = 0 : available but not positive defined + status = 1 : covariance only approximate + status = 2 : full matrix but forced pos def + status = 3 : full accurate matrix + + */ + virtual int CovMatrixStatus() const { return 3; } + /** + return correlation coefficient between variable i and j. + If the variable is fixed or const the return value is zero + */ + virtual double Correlation(unsigned int i, unsigned int j ) const; + + /** + get global correlation coefficient for the variable i. This is a number between zero and one which gives + the correlation between the i-th variable and that linear combination of all other variables which + is most strongly correlated with i. + If the variable is fixed or const the return value is zero + */ + virtual double GlobalCC(unsigned int i) const; + + /** + get the minos error for parameter i, return false if Minos failed + A minimizaiton must be performed befre, return false if no minimization has been done + In case of Minos failed the status error is updated as following + status += 10 * minosStatus where the minos status is: + status = 1 : maximum number of function calls exceeded when running for lower error + status = 2 : maximum number of function calls exceeded when running for upper error + status = 3 : new minimum found when running for lower error + status = 4 : new minimum found when running for upper error + status = 5 : any other failure + + */ + virtual bool GetMinosError(unsigned int i, double & errLow, double & errUp, int = 0); + + /** + scan a parameter i around the minimum. A minimization must have been done before, + return false if it is not the case + */ + virtual bool Scan(unsigned int i, unsigned int & nstep, double * x, double * y, double xmin = 0, double xmax = 0); + + /** + find the contour points (xi,xj) of the function for parameter i and j around the minimum + The contour will be find for value of the function = Min + ErrorUp(); + */ + virtual bool Contour(unsigned int i, unsigned int j, unsigned int & npoints, double *xi, double *xj); + + + /** + perform a full calculation of the Hessian matrix for error calculation + If a valid minimum exists the calculation is done on the minimum point otherwise is performed + in the current set values of parameters + Status code of minimizer is updated according to the following convention (in case Hesse failed) + status += 100*hesseStatus where hesse status is: + status = 1 : hesse failed + status = 2 : matrix inversion failed + status = 3 : matrix is not pos defined + */ + //virtual bool Hesse(); + + + /// return reference to the objective function + //virtual const ROOT::Math::IGenFunction & Function() const; + + /// print result of minimization + virtual void PrintResults(); + + /// set an object to trace operation for each iteration + /// The object muust implement operator() (unsigned int, MinimumState & state) + //void SetTraceObject(MnTraceObject & obj); + + /// set storage level = 1 : store all iteration states (default) + /// = 0 : store only first and last state to save memory + // N/A + void SetStorageLevel(int level); + + private: + unsigned int fDim = 0; // dimension of the function to be minimized + unsigned int fFreeDim = 0; // Number of free dimensions. + std::string fMinimizer = "cmaes"; // minimizer algo. + const ROOT::Math::IMultiGenFunction *fObjFunc = nullptr; + const ROOT::Math::IMultiGradFunction *fObjFuncGrad = nullptr; + std::vector fLBounds; // Lower bounds of variables + std::vector fUBounds; // Upper bounds of variables + std::vector fVariablesType; // 0 for free variable, 1 for fixed variable, 2 for lower bounded, 3 for upper bounded, 4 for lower and upper bounded. + std::vector fInitialX; + std::vector fNames; // Names of the variables. + std::vector fInitialSigma; // User-set Initial step-size for each variables. + std::map fFixedVariables; // fixed variables and values. + CMASolutions fCMAsols; + CMAParameters> fCMAparams; // params no bounds. + CMAParameters> fCMAparams_b; // params with bounds. + CMAParameters> fCMAparams_l; // params no bounds + linear scaling. + CMAParameters> fCMAparams_lb; // params with bounds + linear scaling. + mutable std::vector fGlobalCC; // vector of global correlation coefficients. + mutable std::vector fValues; // X values. + mutable std::vector fErrors; // X errors. + bool fWithBounds = false; // whether using box-type constraints as required by parameters. + bool fWithGradient = false; // whether to use gradient information when available. + int fWithLinearScaling = 0; // wheter to use linear scaling of objective function parameters. */ + }; + + } // end namespace cmaes +} // end namespace ROOT + +#endif diff --git a/math/cmaes/inc/LinkDef.h b/math/cmaes/inc/LinkDef.h new file mode 100644 index 0000000000000..c5a1c5e2825ba --- /dev/null +++ b/math/cmaes/inc/LinkDef.h @@ -0,0 +1,30 @@ +/** + * CMA-ES, Covariance Matrix Evolution Strategy + * Copyright (c) 2014 INRIA + * Author: Emmanuel Benazera + * + * This file is part of libcmaes. + * + * libcmaes is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * libcmaes 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with libcmaes. If not, see . + */ + +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class TCMAESMinimizer; + +#endif diff --git a/math/cmaes/src/CMAESMinimizer.cxx b/math/cmaes/src/CMAESMinimizer.cxx new file mode 100644 index 0000000000000..0d074f800954b --- /dev/null +++ b/math/cmaes/src/CMAESMinimizer.cxx @@ -0,0 +1,771 @@ +/** + * CMA-ES, Covariance Matrix Evolution Strategy + * Copyright (c) 2014 INRIA + * Author: Emmanuel Benazera + * + * This file is part of libcmaes. + * + * libcmaes is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * libcmaes 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 Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with libcmaes. If not, see . + */ + +#include "CMAESMinimizer.h" +#include "Math/IFunctionfwd.h" // fObjFunc +#include "Math/IOptions.h" +#include "Math/Error.h" +#include "Fit/ParameterSettings.h" + +#include "errstats.h" // libcmaes extras. + +#ifdef USE_ROOT_ERROR +#include "TROOT.h" +#endif + +#include + +using namespace libcmaes; + +namespace ROOT +{ + namespace cmaes + { + + TCMAESMinimizer::TCMAESMinimizer() + :Minimizer(),fDim(0),fFreeDim(0),fWithBounds(false),fWithGradient(false) + { + } + + TCMAESMinimizer::TCMAESMinimizer(const char *type) + :Minimizer(),fDim(0),fFreeDim(0),fWithBounds(false),fWithGradient(false) + { + std::string algoname(type); + // tolower() is not an std function (Windows) + std::transform(algoname.begin(), algoname.end(), algoname.begin(), (int(*)(int)) tolower ); + fMinimizer = algoname; + } + + /*TCMAESMinimizer::TCMAESMinimizer(const TCMAESMinimizer &m) + :Minimizer(),fDim(0),fFreeDim(0),fWithBounds(false),fWithGradient(false) + { + }*/ + + TCMAESMinimizer& TCMAESMinimizer::operator = (const TCMAESMinimizer &rhs) + { + if (this == &rhs) return *this; + return *this; + } + + TCMAESMinimizer::~TCMAESMinimizer() + { + } + + void TCMAESMinimizer::Clear() + { + fCMAsols = CMASolutions(); + fCMAparams = CMAParameters>(); + fCMAparams_b = CMAParameters>(); + fCMAparams_l = CMAParameters>(); + fCMAparams_lb = CMAParameters>(); + fDim = 0; fFreeDim = 0; + fLBounds.clear(); + fUBounds.clear(); + fVariablesType.clear(); + fInitialX.clear(); + fInitialSigma.clear(); + fFixedVariables.clear(); + fNames.clear(); + fGlobalCC.clear(); + fValues.clear(); + fErrors.clear(); + } + + void TCMAESMinimizer::SetFunction(const ROOT::Math::IMultiGenFunction &fun) + { + fObjFunc = &fun; + fDim = fun.NDim(); + } + + void TCMAESMinimizer::SetFunction(const ROOT::Math::IMultiGradFunction &fun) + { + SetFunction(static_cast (fun)); + fObjFuncGrad = &fun; + //fDim = fun.NDim(); + fWithGradient = true; + } + + bool TCMAESMinimizer::SetVariable(unsigned int ivar, const std::string & name, double val, double step) + { + if (ivar > fInitialX.size() ) { + MATH_ERROR_MSG("TCMAESMinimizer::SetVariable","ivar out of range"); + return false; + } + if (ivar == fInitialX.size() ) { + fInitialX.push_back(val); + fNames.push_back(name); + fInitialSigma.push_back(step); + fLBounds.push_back(-std::numeric_limits::max()); + fUBounds.push_back(std::numeric_limits::max()); + if (step==0.){ + fVariablesType.push_back(1); + } + else { + fFreeDim++; + fVariablesType.push_back(0); + } + } + else { + if (step==0.) { + if (fInitialSigma[ivar]!=0.) { //Constraining a free variable. + fFreeDim--; + fVariablesType[ivar] = 1; + } + } + else { + if (fInitialSigma[ivar]==0.) { //Freeing a constrained variable + fFreeDim++; + fVariablesType[ivar] = 0; + } + } + fInitialX[ivar] = val; + fNames[ivar] = name; + fInitialSigma[ivar] = step; + } + return true; + } + + bool TCMAESMinimizer::SetLowerLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower ) + { + if (lower > val) { + MATH_WARN_MSG("TCMAESMinimizer::SetLowerLimitedVariable", "Starting point set into the unfeasible domain"); // fix with val=lower; ? + } + bool r = SetVariable(ivar, name, val, step); + if (!r) return false; + fLBounds[ivar] = lower; + fVariablesType[ivar] = 2; + fWithBounds = true; + return true; + } + + bool TCMAESMinimizer::SetUpperLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double upper ) + { + if (upper > val) { + MATH_WARN_MSG("TCMAESMinimizer::SetUpperLimitedVariable", "Starting point set into the unfeasible domain"); + } + bool r = SetVariable(ivar, name, val, step); + if (!r) return false; + fUBounds[ivar] = upper; + fVariablesType[ivar] = 3; + fWithBounds = true; + return true; + } + + bool TCMAESMinimizer::SetLimitedVariable(unsigned int ivar , const std::string & name , double val , double step , double lower, double upper) + { + if (upper == lower) { + MATH_WARN_MSG("TCMAESMinimizer::SetLimitedVariable","Upper bound equal to lower bound. Variable is constrained to fixed value."); + return SetFixedVariable(ivar, name, val); + } + if (upper < lower) { + MATH_WARN_MSG("TCMAESMinimizer::SetLimitedVariable","Upper bound lesser than lower bound. Bounds exchanged."); + double temp(upper); + upper = lower; + lower = temp; + } + if (val < lower || val > upper) { + MATH_WARN_MSG("TCMAESMinimizer::SetLimitedVariable", "Starting point set into the unfeasible domain"); + } + bool r = SetVariable(ivar, name, val, step); + if (!r) return false; + fLBounds[ivar] = lower; + fUBounds[ivar] = upper; + fVariablesType[ivar] = 4; + fWithBounds = true; + return true; + } + + bool TCMAESMinimizer::SetFixedVariable(unsigned int ivar, const std::string &name, double val) + { + SetVariable(ivar,name,val,0.0); + fFixedVariables.insert(std::pair(ivar,val)); + return true; + } + + bool TCMAESMinimizer::SetVariableValue(unsigned int ivar, double val ) + { + if (ivar >= fInitialX.size() ) { + //TODO string that gives value of ivar and fInitialX.size() + MATH_ERROR_MSG("TCMAESMinimizer::SetVariableValue","ivar out of range"); + return false; + } + if (fVariablesType[ivar] == 2 || fVariablesType[ivar] == 4) { + if (fLBounds[ivar] > val) { + MATH_WARN_MSG("TCMAESMinimizer::SetVariableValue", "Starting point set into the unfeasible domain"); + } + } + if (fVariablesType[ivar] == 3 || fVariablesType[ivar] == 4) { + if (fUBounds[ivar] < val) { + MATH_WARN_MSG("TCMAESMinimizer::SetVariableValue", "Starting point set into the unfeasible domain"); + } + } + fInitialX[ivar] = val; + return true; + } + + bool TCMAESMinimizer::SetVariableValues(const double * x) + { + if (x == NULL) + { + MATH_WARN_MSG("TCMAESMinimizer::SetVariableValues", "No values given, no change to the starting point."); + return false; + } + unsigned int i; + for (i=0; i fInitialX.size()) + return false; + fInitialSigma[ivar] = step; + return true; + } + + bool TCMAESMinimizer::SetVariableLowerLimit(unsigned int ivar, double lower) + { + if (ivar > fLBounds.size()) + return false; + fLBounds[ivar] = lower; + fVariablesType[ivar] = 2; + fWithBounds = true; + return true; + } + + bool TCMAESMinimizer::SetVariableUpperLimit(unsigned int ivar, double upper) + { + if (ivar > fUBounds.size()) + return false; + fUBounds[ivar] = upper; + fVariablesType[ivar] = 3; + fWithBounds = true; + return true; + } + + bool TCMAESMinimizer::SetVariableLimits(unsigned int ivar, double lower, double upper) + { + if (ivar >= fLBounds.size() || ivar >= fUBounds.size()) + return false; + fLBounds[ivar] = lower; + fUBounds[ivar] = upper; + fVariablesType[ivar] = 4; + fWithBounds = true; + return true; + } + + bool TCMAESMinimizer::FixVariable(unsigned int ivar) + { + if (ivar >= fInitialX.size()) + return false; + fFixedVariables.insert(std::pair(ivar,fInitialX.at(ivar))); // XXX: sets initial variable. + return true; + } + + bool TCMAESMinimizer::IsFixedVariable(unsigned int ivar) const + { + std::map::const_iterator mit; + if ((mit=fFixedVariables.find(ivar))!=fFixedVariables.end()) + return true; + return false; + } + + bool TCMAESMinimizer::GetVariableSettings(unsigned int ivar, ROOT::Fit::ParameterSettings &varObj) const + { + if (ivar >= fInitialX.size()) + { + MATH_ERROR_MSG("TCMAESMinimizer::GetVariableSettings","wrong variable index"); + return false; + } + varObj.Set(fNames.at(ivar),fInitialX.at(ivar),false); //XXX: not sure of last param type. + if (fVariablesType.at(ivar) == 4) + varObj.SetLimits(fLBounds.at(ivar),fUBounds.at(ivar)); + else if (fVariablesType.at(ivar) == 3) + varObj.SetUpperLimit(fUBounds.at(ivar)); + else if (fVariablesType.at(ivar) == 2) + varObj.SetLowerLimit(fLBounds.at(ivar)); + return true; + } + + std::string TCMAESMinimizer::VariableName(unsigned int ivar) const + { + if (ivar >= fInitialX.size()) + return std::string(); + return fNames.at(ivar); + } + + int TCMAESMinimizer::VariableIndex(const std::string &name) const + { + for (unsigned int i=0;i + void TCMAESMinimizer::SetMParameters(CMAParameters &cmaparams, + const int &maxiter, const int &maxfevals, + const int &noisy, const int &nrestarts, + const double &ftarget, + const std::string &fplot, + const bool &withnumgradient, + const bool &mtfeval, + const bool &quiet, + const int &elitist, + const bool &uh) + { + cmaparams.set_str_algo(fMinimizer); + if (gDebug > 0 || !quiet) + cmaparams.set_quiet(false); + else cmaparams.set_quiet(true); + for (auto mit=fFixedVariables.begin();mit!=fFixedVariables.end();mit++) + cmaparams.set_fixed_p((*mit).first,(*mit).second); + cmaparams.set_edm(true); // always activate EDM computation. + cmaparams.set_ftolerance(Tolerance()); + cmaparams.set_max_iter(maxiter); + cmaparams.set_max_fevals(maxfevals); + if (noisy > 0) + cmaparams.set_noisy(); + if (nrestarts > 0) + cmaparams.set_restarts(nrestarts); + if (ftarget > 0.0) + cmaparams.set_ftarget(ftarget); + cmaparams.set_fplot(fplot); + cmaparams.set_gradient(withnumgradient); + cmaparams.set_mt_feval(mtfeval); + cmaparams.set_elitism(elitist); + cmaparams.set_uh(uh); + } + + bool TCMAESMinimizer::Minimize() + { + if (!fObjFunc) { + MATH_ERROR_MSG("TCMAESMinimizer::Minimize","Objective function has not been set"); + return false; + } + if (!fDim) { + MATH_ERROR_MSG("TCMAESMinimizer::Minimize","Dimension has not been set"); + return false; + } + if (fDim > fInitialX.size()) { + std::cout << "fDim=" << fDim << " / fInitialX size=" << fInitialX.size() << " / freeDim=" << fFreeDim << std::endl; + MATH_ERROR_MSG("TCMAESMinimizer::Minimize","Dimension larger than initial X size's"); + return false; + } + if (fDim < fInitialX.size()) { + MATH_WARN_MSG("TCMAESMinimizer::Minimize","Dimension smaller than initial X size's"); + } + + ROOT::Math::IOptions *cmaesOpt = ROOT::Math::MinimizerOptions::FindDefault("cmaes"); + //std::cerr << "cmaesOpt ptr: " << cmaesOpt << std::endl; + if (cmaesOpt) + cmaesOpt->Print(std::cout); + + FitFunc ffit = [this](const double *x, const int N) + { + /*std::copy(x,x+N,std::ostream_iterator(std::cout," ")); + std::cout << std::endl;*/ + (void)N; + return (*fObjFunc)(x); + }; + + // gradient function. + //std::cout << "fWithGradient=" << fWithGradient << std::endl; + GradFunc gfit = nullptr; + if (fWithGradient) + { + gfit = [this](const double *x, const int N) + { + dVec grad(N); + fObjFuncGrad->Gradient(x,grad.data()); + return grad; + }; + } + + //debug + /*if (fWithBounds) + { + std::cout << "bounds:\n"; + std::copy(fLBounds.begin(),fLBounds.end(),std::ostream_iterator(std::cout," ")); + std::cout << std::endl; + std::copy(fUBounds.begin(),fUBounds.end(),std::ostream_iterator(std::cout," ")); + std::cout << std::endl; + }*/ + //debug + + double sigma0 = *std::min_element(fInitialSigma.begin(),fInitialSigma.end()); + double sigma0scaled = 1e-1; // default value. + if (!fWithLinearScaling) + sigma0scaled = sigma0; + dVec vscaling = dVec::Constant(fDim,1.0); + for (size_t i=0;iGetValue("lambda",lambda); + cmaesOpt->GetValue("noisy",noisy); + cmaesOpt->GetValue("restarts",nrestarts); + cmaesOpt->GetValue("ftarget",ftarget); + cmaesOpt->GetValue("fplot",fplot); + cmaesOpt->GetValue("lscaling",fWithLinearScaling); + cmaesOpt->GetValue("numgradient",withnumgradient); + cmaesOpt->GetValue("mt_feval",mtfeval); + cmaesOpt->GetValue("quiet",quiet); + cmaesOpt->GetValue("seed",seed); + cmaesOpt->GetValue("elitist",elitist); + cmaesOpt->GetValue("uh",uh); + } + + if (gDebug > 0) + { + std::cout << "Running CMA-ES with dim=" << fDim << " / sigma0=" << sigma0scaled << " / lambda=" << lambda << " / fTol=" << Tolerance() << " / with_bounds=" << fWithBounds << " / with_gradient=" << fWithGradient << " / linear_scaling=" << fWithLinearScaling << " / maxiter=" << maxiter << " / maxfevals=" << maxfevals << " / mtfeval=" << mtfeval << std::endl; + std::cout << "x0="; + std::copy(fInitialX.begin(),fInitialX.end(),std::ostream_iterator(std::cout," ")); + std::cout << std::endl; + } + + if (fWithLinearScaling) + { + if (fWithBounds) + { + Info("CMAESMinimizer","Minimizing with bounds and linear scaling"); + GenoPheno gp(vscaling,vshift,&fLBounds.front(),&fUBounds.front()); + CMAParameters> cmaparams(fDim,&fInitialX.front(),sigma0scaled,lambda,seed,gp); + SetMParameters(cmaparams,maxiter,maxfevals,noisy,nrestarts,ftarget,fplot,withnumgradient,mtfeval,quiet,elitist,uh); + fCMAsols = libcmaes::cmaes>(ffit,cmaparams,CMAStrategy>::_defaultPFunc,fWithGradient?gfit:nullptr); + fCMAparams_lb = cmaparams; + } + else + { + Info("CMAESMinimizer","Minimizing with linear scaling"); + GenoPheno gp(vscaling,vshift); + CMAParameters> cmaparams(fDim,&fInitialX.front(),sigma0scaled,lambda,seed,gp); + SetMParameters(cmaparams,maxiter,maxfevals,noisy,nrestarts,ftarget,fplot,withnumgradient,mtfeval,quiet,elitist,uh); + fCMAsols = libcmaes::cmaes>(ffit,cmaparams,CMAStrategy>::_defaultPFunc,fWithGradient?gfit:nullptr); + fCMAparams_l = cmaparams; + } + } + else + { + if (fWithBounds) + { + Info("CMAESMinimizer","Minimizing with bounds"); + GenoPheno gp(&fLBounds.front(),&fUBounds.front(),fDim); + CMAParameters> cmaparams(fDim,&fInitialX.front(),sigma0scaled,lambda,seed,gp); + SetMParameters(cmaparams,maxiter,maxfevals,noisy,nrestarts,ftarget,fplot,withnumgradient,mtfeval,quiet,elitist,uh); + fCMAsols = libcmaes::cmaes>(ffit,cmaparams,CMAStrategy>::_defaultPFunc,fWithGradient?gfit:nullptr); + fCMAparams_b = cmaparams; + } + else + { + Info("CMAESMinimizer","Minimizing without bounds or linear scaling"); + CMAParameters> cmaparams(fDim,&fInitialX.front(),sigma0scaled,lambda,seed); + SetMParameters(cmaparams,maxiter,maxfevals,noisy,nrestarts,ftarget,fplot,withnumgradient,mtfeval,quiet,elitist,uh); + fCMAsols = libcmaes::cmaes>(ffit,cmaparams,CMAStrategy>::_defaultPFunc,fWithGradient?gfit:nullptr); + fCMAparams = cmaparams; + } + } + Info("CMAESMinimizer","optimization status=%i",fCMAsols.run_status()); + if (fCMAsols.edm() > 10*Tolerance()) // XXX: max edm seems to be left to each minimizer's internal implementation... + fStatus = 3; + else if (fCMAsols.run_status() == 0 || fCMAsols.run_status() == 1) + fStatus = 0; + else if (fCMAsols.run_status() == 7 || fCMAsols.run_status() == 9) + fStatus = 4; // reached budget limit. + else fStatus = 5; + return fCMAsols.run_status() >= 0; // above 0 are partial successes at worst. + } + + double TCMAESMinimizer::MinValue() const + { + return fCMAsols.best_candidate().get_fvalue(); + } + + const double* TCMAESMinimizer::X() const + { + fValues.clear(); + Candidate bc = fCMAsols.best_candidate(); + + dVec x; + if (fWithLinearScaling) + { + if (fWithBounds) + x = bc.get_x_pheno_dvec>(fCMAparams_lb); + else x = bc.get_x_pheno_dvec>(fCMAparams_l); + } + else + { + if (fWithBounds) + x = bc.get_x_pheno_dvec>(fCMAparams_b); + else x = bc.get_x_dvec(); + } + for (int i=0;i<(int)fDim;i++) + fValues.push_back(x(i)); + return &fValues.front(); + } + + double TCMAESMinimizer::Edm() const + { + // XXX: cannot recompute it here as there's no access to the optimizer itself. + // instead this is returning the value computed at the end of last optimization call + // and stored within the solution object. + return fCMAsols.edm(); + } + + const double* TCMAESMinimizer::Errors() const + { + fErrors.clear(); + dVec verrors; + if (fWithLinearScaling) + { + if (fWithBounds) + verrors = fCMAsols.errors(fCMAparams_lb); + else verrors = fCMAsols.errors(fCMAparams_l); + } + else + { + if (fWithBounds) + verrors = fCMAsols.errors(fCMAparams_b); + else verrors = fCMAsols.errors(fCMAparams); + } + for (int i=0;i<(int)fDim;i++) + fErrors.push_back(verrors(i)); + return &fErrors.front(); + } + + unsigned int TCMAESMinimizer::NCalls() const + { + return fCMAsols.nevals(); + } + + double TCMAESMinimizer::CovMatrix(unsigned int i, unsigned int j) const + { + return fCMAsols.cov_ref()(i,j); + } + + bool TCMAESMinimizer::GetCovMatrix(double *cov) const + { + std::copy(fCMAsols.cov_data(),fCMAsols.cov_data()+fCMAsols.cov_ref().size(),cov); + return true; + } + + double TCMAESMinimizer::Correlation(unsigned int i, unsigned int j) const + { + return std::sqrt(std::abs(fCMAsols.cov_ref()(i,i)*fCMAsols.cov_ref()(j,j))); + } + + double TCMAESMinimizer::GlobalCC(unsigned int i) const + { + // original Minuit paper says: + // \rho_k^2 = 1 - [C_{kk}C_{kk}^{-1}]^{-1} + if (fGlobalCC.empty()) // need to pre-compute the vector coefficient + { + dMat covinv = fCMAsols.cov_ref().inverse(); + for (int i=0;i 0.0) + fGlobalCC.push_back(0.0); + else fGlobalCC.push_back(std::sqrt(1.0 - 1.0/denom)); + } + } + return fGlobalCC.at(i); + } + + bool TCMAESMinimizer::GetMinosError(unsigned int i, double &errLow, double &errUp, int j) + { + (void)j; + FitFunc ffit = [this](const double *x, const int N) + { + (void)N; + return (*fObjFunc)(x); + }; + + // runopt is a flag which specifies if only lower or upper error needs to be run. TODO: support for one bound only in libcmaes ? + int samplesize = 10; + if (gDebug > 0) + std::cerr << "Computing 'Minos' confidence interval with profile likelihood on parameter " << i << " / samplesize=" << samplesize << " / with_bounds=" << fWithBounds << std::endl; + pli le; + if (fWithLinearScaling) + { + if (!fWithBounds) + { + le = errstats>::profile_likelihood(ffit,fCMAparams_l,fCMAsols,i,false,samplesize,ErrorDef(),100); + } + else + { + le = errstats>::profile_likelihood(ffit,fCMAparams_lb,fCMAsols,i,false,samplesize); + } + } + else + { + if (!fWithBounds) + { + le = errstats>::profile_likelihood(ffit,fCMAparams,fCMAsols,i,false,samplesize,ErrorDef()); + } + else + { + le = errstats>::profile_likelihood(ffit,fCMAparams_b,fCMAsols,i,false,samplesize); + } + } + errLow = le.get_err_min(); + errUp = le.get_err_max(); + return true; + } + + bool TCMAESMinimizer::Scan(unsigned int i, unsigned int &nstep, double *x, double *y, double xmin, double xmax) + { + std::vector> result; + std::vector params = fValues; + double amin = MinValue(); + result.push_back(std::pair(params[i],amin)); + + double low=xmin, high=xmax; + if (low <= high && nstep-1 >= 2) + { + if (low == 0 && high == 0) + { + low = fValues[i] - 2.0*fErrors.at(i); + high = fValues[i] + 2.0*fErrors.at(i); + } + + if (low == 0 && high == 0 + && (fLBounds[i] > -std::numeric_limits::max() + || fUBounds[i] < std::numeric_limits::max())) + { + if (fLBounds[i] > -std::numeric_limits::max()) + low = fLBounds[i]; + if (fUBounds[i] < std::numeric_limits::max()) + high = fUBounds[i]; + } + + if (fLBounds[i] > -std::numeric_limits::max() + || fUBounds[i] < std::numeric_limits::max()) + { + if (fLBounds[i] > -std::numeric_limits::max()) + low = std::max(low,fLBounds[i]); + if (fUBounds[i] < std::numeric_limits::max()) + high = std::min(high,fUBounds[i]); + } + + double x0 = low; + double stp = (high-low) / static_cast(nstep-2); + for (unsigned int j=0;j(params[i],fval)); + } + } + + for (int s=0;s<(int)nstep;s++) + { + x[s] = result[s].first; + y[s] = result[s].second; + } + return true; + } + + bool TCMAESMinimizer::Contour(unsigned int i, unsigned int j, unsigned int &npoints, double *xi, double *xj) + { + FitFunc ffit = [this](const double *x, const int N) + { + (void)N; + return (*fObjFunc)(x); + }; + + contour ct; + if (fWithLinearScaling) + { + if (!fWithBounds) + { + ct = errstats>::contour_points(ffit,i,j,npoints,ErrorDef(), + fCMAparams_l,fCMAsols,0.1,100); + } + else + { + ct = errstats>::contour_points(ffit,i,j,npoints,ErrorDef(), + fCMAparams_lb,fCMAsols,0.1,100); + } + } + else + { + if (!fWithBounds) + { + ct = errstats>::contour_points(ffit,i,j,npoints,ErrorDef(), + fCMAparams,fCMAsols,0.1,100); + } + else + { + ct = errstats>::contour_points(ffit,i,j,npoints,ErrorDef(), + fCMAparams_b,fCMAsols,0.1,100); + } + } + for (size_t i=0;i 1) + std::cout << "\t(limited)"; + std::cout << std::endl; + } + } + + } +} diff --git a/math/cmaes/test/CMakeLists.txt b/math/cmaes/test/CMakeLists.txt new file mode 100644 index 0000000000000..e41053d5999d0 --- /dev/null +++ b/math/cmaes/test/CMakeLists.txt @@ -0,0 +1,66 @@ +project(cmaes-tests) +find_package(ROOT REQUIRED) + +include(${ROOT_USE_FILE}) +include_directories(${ROOT_INCLUDE_DIRS}) + +set(TestSource + testMinimizer.cxx +) + +set(TestSourceMnTutorial + MnTutorial/Quad1FMain.cxx + MnTutorial/Quad4FMain.cxx + MnTutorial/Quad8FMain.cxx + MnTutorial/Quad12FMain.cxx +) + +set(TestSourceMnSim + MnSim/DemoGaussSim.cxx + MnSim/DemoFumili.cxx + MnSim/PaulTest.cxx + MnSim/PaulTest2.cxx + MnSim/PaulTest3.cxx + MnSim/PaulTest4.cxx + MnSim/ReneTest.cxx + MnSim/ParallelTest.cxx + MnSim/demoMinimizer.cxx +) + + + + +#---For the simple cmaes tests build and defined them--------------- +foreach(file ${TestSourceMnTutorial}) + get_filename_component(testname ${file} NAME_WE) + ROOT_EXECUTABLE(${testname} ${file} LIBRARIES cmaes) + ROOT_ADD_TEST(minuit2-${testname} COMMAND ${testname}) +endforeach() + + +ROOT_LINKER_LIBRARY(Minuit2TestMnSim MnSim/GaussDataGen.cxx MnSim/GaussFcn.cxx MnSim/GaussFcn2.cxx LIBRARIES Minuit2) + +#input text files +configure_file(MnSim/paul.txt paul.txt @COPY_ONLY) +configure_file(MnSim/paul2.txt paul2.txt @COPY_ONLY) +configure_file(MnSim/paul3.txt paul3.txt @COPY_ONLY) +configure_file(MnSim/paul4.txt paul4.txt @COPY_ONLY) + +foreach(file ${TestSourceMnSim}) + get_filename_component(testname ${file} NAME_WE) + ROOT_EXECUTABLE(${testname} ${file} LIBRARIES Minuit2 Minuit2TestMnSim) + ROOT_ADD_TEST(minuit2-${testname} COMMAND ${testname}) +endforeach() + +#for the global tests using ROOT libs (Minuit2 should be taken via the PluginManager) + +set(RootLibraries Core RIO Net Hist Graf Graf3d Gpad Tree + Rint Postscript Matrix Physics MathCore Thread) + +foreach(file ${TestSource}) + get_filename_component(testname ${file} NAME_WE) + ROOT_EXECUTABLE(${testname} ${file} LIBRARIES ${RootLibraries} ) + ROOT_ADD_TEST(minuit2-${testname} COMMAND ${testname}) +endforeach() + + diff --git a/math/cmaes/test/Makefile b/math/cmaes/test/Makefile new file mode 100644 index 0000000000000..077c535c41002 --- /dev/null +++ b/math/cmaes/test/Makefile @@ -0,0 +1,91 @@ +# Makefile for the ROOT test programs. +# This Makefile shows nicely how to compile and link applications +# using the ROOT libraries on all supported platforms. +# +# Copyright (c) 2000 Rene Brun and Fons Rademakers +# +# Author: Fons Rademakers, 29/2/2000 + +ROOTSYS = ../../.. +include $(ROOTSYS)/etc/Makefile.arch + +#------------------------------------------------------------------------------ + +ifeq ($(PLATFORM),win32) +EXTRALIBS = "$(ROOTSYS)/lib/libcmaesroot.lib" +else +EXTRALIBS = -lcmaesroot -lglog -lgflags +CXXFLAGS += -g -L/home/beniz/lib -L/usr/lib -L/usr/lib/x86_64-linux-gnu/ +endif + +# for using with MPI +ifneq ($(USE_MPI),) +CXX=mpic++ +LD=mpic++ +endif +ifneq ($(USE_OPENMP),) +CXXFLAGS += -D_GLIBCXX_PARALLEL -fopenmp +LDFLAGS += -fopenmp +endif + + + +USERFUNCOBJ = testUserFunc.$(ObjSuf) +USERFUNCSRC = testUserFunc.$(SrcSuf) +USERFUNC = testUserFunc$(ExeSuf) + +MINIMIZEROBJ = testMinimizer.$(ObjSuf) +MINIMIZERSRC = testMinimizer.$(SrcSuf) +MINIMIZER = testMinimizer$(ExeSuf) + +NDIMFITOBJ = testNdimFit.$(ObjSuf) +NDIMFITSRC = testNdimFit.$(SrcSuf) +NDIMFIT = testNdimFit$(ExeSuf) + +GAUSFITOBJ = testUnbinGausFit.$(ObjSuf) +GAUSFITSRC = testUnbinGausFit.$(SrcSuf) +GAUSFIT = testUnbinGausFit$(ExeSuf) + + +OBJS = $(USERFUNCOBJ) $(GRAPHOBJ) $(MINIMIZEROBJ) $(NDIMFITOBJ) $(GAUSFITOBJ) + +PROGRAMS = $(USERFUNC) $(GRAPH) $(MINIMIZER) $(NDIMFIT) $(GAUSFIT) + +.SUFFIXES: .$(SrcSuf) .$(ObjSuf) $(ExeSuf) + + +all: $(PROGRAMS) + +$(USERFUNC): $(USERFUNCOBJ) + $(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@ + @echo "$@ done" + +$(MINIMIZER): $(MINIMIZEROBJ) +ifeq ($(PLATFORM),win32) + $(LD) $(LDFLAGS) $^ $(LIBS) "$(ROOTSYS)/lib/libcmaesroot.lib" $(OutPutOpt)$@ +else + $(LD) $(LDFLAGS) $^ $(LIBS) -lMathCore $(OutPutOpt)$@ +endif + @echo "$@ done" + +$(NDIMFIT): $(NDIMFITOBJ) + $(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@ + @echo "$@ done" + +$(GAUSFIT): $(GAUSFITOBJ) + $(LD) $(LDFLAGS) $^ $(LIBS) $(EXTRALIBS) $(OutPutOpt)$@ + @echo "$@ done" + + +clean: + @rm -f $(OBJS) core + +distclean: clean + @rm -f $(PROGRAMS) *Dict.* *.def *.exp \ + *.root *.ps *.so *.lib *.dll *.d .def so_locations + +.SUFFIXES: .$(SrcSuf) + + +.$(SrcSuf).$(ObjSuf): + $(CXX) $(CXXFLAGS) -c $< diff --git a/math/cmaes/test/cmaesFullBench.py b/math/cmaes/test/cmaesFullBench.py new file mode 100644 index 0000000000000..cba3a55f8afbd --- /dev/null +++ b/math/cmaes/test/cmaesFullBench.py @@ -0,0 +1,104 @@ +import sys, csv +import numpy as np +import matplotlib.pyplot as plt +from numpy import * + +datfiles = ["combined.dat","example3D.dat","fit2a.dat","fit2.dat","fit2dhist.dat","gauss2D_fit.dat","gauss_fit.dat","lorentz_fit.dat"] + +fig, axarr = plt.subplots(8,2) + +i = 0 +j = 0 +N = 2 +runs = -1 +for f in datfiles: + dat = loadtxt(f,dtype=float,comments='#') + nalgs = len(dat) + + if runs == -1: + runs = dat[0,3]+dat[0,4] # get the number of runs once. + + hcolor = 'Teal' + aa = 0 + oldrec = None + for a in range(nalgs): + if a == nalgs-1: + hcolor = 'LightGreen' + + aavg1 = [dat[a,2],dat[a,9]] + aastd1 = [0.0,0.0] + aavg2 = [dat[a,5]*1000.0,dat[a,7]] + aastd2 = [dat[a,6]*1000.0,dat[a,8]] + + ## necessary variables + ind = np.arange(N) # the x locations for the groups + width = 0.05 # the width of the bars + + ## the bars + rects1 = axarr[i,j+1].bar(ind+aa*width, aavg1, width, + color=hcolor, + yerr=aastd1, + error_kw=dict(elinewidth=2,ecolor='red')) + if a == 0: + oldrec = rects1 + + #rects2 = ax2.bar(ind+width, minAvg, width, + # color='LightGreen', + # yerr=minStd, + # error_kw=dict(elinewidth=2,ecolor='black')) + + rects12 = axarr[i,j].bar(ind+aa*width, aavg2, width, + color=hcolor, + yerr=aastd2, + error_kw=dict(elinewidth=2,ecolor='red')) + + #rects22 = axarr[i,j].bar(ind+width, minAvg2, width, + # color='LightGreen', + # yerr=minStd2, + # error_kw=dict(elinewidth=2,ecolor='black')) + + #add a legend + if i == 6: + axarr[6,1].legend((oldrec[0], rects12[0]), ('aCMA-ES', 'Minuit2'), fontsize=10 ) + aa = aa + 1 + + + + # axes and labels + axarr[i,j].set_yscale("log",nonposy='clip') + #axarr[i,j].set_xlim(-width,len(ind)+width) + axarr[i,j].set_ylim(0) + #axarr[i,j+1].set_xlim(-width,len(ind)+width) + axarr[i,j+1].set_ylim(top=runs) + if i == 4 or i == 6: + axarr[i,j+1].set_ylim(top=2*runs) + if i == 7: + axarr[i,j+1].set_ylim(top=100) + #axarr[i,j].set_ylabel('Scores') + axarr[i,j].set_title(f + " / " + str(dat[0,0]) + "-D",fontsize=12) + #axarr[i,j].set_xticks(ind+width) + #axarr[i,j+1].set_xticks(ind+width) + plt.setp(axarr[i,j].get_xticklabels(),visible=False) + plt.setp(axarr[i,j+1].get_xticklabels(),visible=False) + if i == 7: + xTickMarks1 = ["CPU avg (ms)","","","","","Budget avg"] + xTickMarks2 = ["Found","","","","","wins"] + xtickNames = axarr[7,j].set_xticklabels(xTickMarks1) + plt.setp(xtickNames, rotation=45, fontsize=10) + xtickNames2 = axarr[7,j+1].set_xticklabels(xTickMarks2) + plt.setp(xtickNames2, rotation=45, fontsize=10) + plt.setp(axarr[i,j].get_xticklabels(),visible=True) + plt.setp(axarr[i,j+1].get_xticklabels(),visible=True) + + #add a legend +# if i == 6: +# axarr[7,1].legend()# (rects1[0], rects12[0]), ('aCMA-ES', 'Minuit2'), fontsize=10 ) + i = i + 1 + j = 0 +# if i == 2: +# i = 0 +# j = j + 1 + +#plt.tight_layout() +plt.suptitle('aCMA-ES / Minuit2 Benchmark Suite / ' + str(len(datfiles)) + ' experiments / ' + str(int(runs)) + ' runs on each\nlambda={auto, 50, 200, auto-aipop-4-restarts, auto-abipop-10-restarts}') #10, 20, 40, 80, 160, 320, 640, 1280}') +plt.show() diff --git a/math/cmaes/test/testMinimizer.cxx b/math/cmaes/test/testMinimizer.cxx new file mode 100644 index 0000000000000..dc42368a9dbff --- /dev/null +++ b/math/cmaes/test/testMinimizer.cxx @@ -0,0 +1,891 @@ +// test of minimization using new minimizer classes + +#include "Math/Minimizer.h" +#include "Math/Factory.h" +#include "Math/Functor.h" + +#include "TVirtualFitter.h" +#include "TSystem.h" + +#include "Math/IFunction.h" +#include "Math/Util.h" +#include +#include + +#include +#include + +#include "TStopwatch.h" +#include "TMatrixD.h" +#include "TVectorD.h" +#include "TRandom3.h" +#include "TMath.h" + +//#define DEBUG + +int gNCall = 0; +int gNCall2 = 0; +int gNmin = 1000; +int gVerbose = 0; +bool useGradient = true; + +bool minos = true; + +double gAbsTolerance = 0.005; + +// Rosenbrok function to be minimize + +typedef void (*FCN)(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag); + + + +// ROSENBROCK function +//______________________________________________________________________________ +void RosenBrock(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t /*iflag*/) +{ + gNCall++; + const Double_t x = par[0]; + const Double_t y = par[1]; + const Double_t tmp1 = y-x*x; + const Double_t tmp2 = 1-x; + f = 100*tmp1*tmp1+tmp2*tmp2; +} + + + +class RosenBrockFunction : public ROOT::Math::IMultiGenFunction { + +public : + + virtual ~RosenBrockFunction() {} + + unsigned int NDim() const { return 2; } + + ROOT::Math::IMultiGenFunction * Clone() const { + return new RosenBrockFunction(); + } + + const double * TrueMinimum() const { + fTrueMin[0] = 1; + fTrueMin[1] = 1; + return fTrueMin; + } + + private: + + inline double DoEval (const double * x) const { +#ifdef USE_FREE_FUNC + double f = 0; + int ierr = 0; + int i = 0; + RosenBrock(i,0,f,const_cast(x),ierr); + return f; +#else + gNCall++; + const Double_t xx = x[0]; + const Double_t yy = x[1]; + const Double_t tmp1 = yy-xx*xx; + const Double_t tmp2 = 1-xx; + return 100*tmp1*tmp1+tmp2*tmp2; +#endif + } + + + mutable double fTrueMin[2]; +}; + + +// TRIGONOMETRIC FLETCHER FUNCTION + +class TrigoFletcherFunction : public ROOT::Math::IMultiGradFunction { + +public : + + + TrigoFletcherFunction(unsigned int dim) : fDim(dim) { + double seed = 3; + A.ResizeTo(dim,dim); + B.ResizeTo(dim,dim); + x0.ResizeTo(dim); + sx0.ResizeTo(dim); + cx0.ResizeTo(dim); + sx.ResizeTo(dim); + cx.ResizeTo(dim); + v0.ResizeTo(dim); + v.ResizeTo(dim); + r.ResizeTo(dim); + A.Randomize(-100.,100,seed); + B.Randomize(-100.,100,seed); + for (unsigned int i = 0; i < dim; i++) { + for (unsigned int j = 0; j < dim; j++) { + A(i,j) = int(A(i,j)); + B(i,j) = int(B(i,j)); + } + } + x0.Randomize(-TMath::Pi(),TMath::Pi(),seed); + // calculate vector Ei + for (unsigned int i = 0; i < fDim ; ++i) { + cx0[i] = std::cos(x0[i]); + sx0[i] = std::sin(x0[i]); + } + v0 = A*sx0+B*cx0; + } + + + unsigned int NDim() const { return fDim; } + + ROOT::Math::IMultiGenFunction * Clone() const { + TrigoFletcherFunction * f = new TrigoFletcherFunction(*this); +// std::cerr <<"cannot clone this function" << std::endl; +// assert(0); + return f; + } + + + void StartPoints(double * x, double * s) { + TRandom3 rndm; + const double stepSize = 0.01; + const double deltaAmp = 0.1; + const double pi = TMath::Pi(); + for (unsigned int i = 0; i < fDim; ++i) { + double delta = rndm.Uniform(-deltaAmp*pi,deltaAmp*pi); + x[i] = x0(i) + 0.1*delta; + if (x[i] <= - pi) x[i] += 2.*pi; + if (x[i] > pi) x[i] -= 2.*pi; + s[i] = stepSize; + } + } + + + const double * TrueMinimum() const { + return x0.GetMatrixArray(); + } + + + void Gradient (const double * x, double * g) const { + gNCall2++; + + for (unsigned int i = 0; i < fDim ; ++i) { + cx [i] = std::cos(x[i]); + sx [i] = std::sin(x[i]); + } + + v = A*sx +B*cx; + r = v0-v; + + + // calculate the grad components + for (unsigned int i = 0; i < fDim ; ++i) { + g[i] = 0; + for (unsigned int k = 0; k < fDim ; ++k) { + g[i] += 2. * r(k) * ( - A(k,i) * cx(i) + B(k,i) * sx(i) ); + } + } + + } + +#ifdef USE_FDF + void FdF (const double * x, double & f, double * g) const { + gNCall++; + + for (unsigned int i = 0; i < fDim ; ++i) { + cx [i] = std::cos(x[i]); + sx [i] = std::sin(x[i]); + } + + v = A*sx +B*cx; + r = v0-v; + + f = r * r; + + + // calculate the grad components + for (unsigned int i = 0; i < fDim ; ++i) { + g[i] = 0; + for (unsigned int k = 0; k < fDim ; ++k) { + g[i] += 2. * r(k) * ( - A(k,i) * cx(i) + B(k,i) * sx(i) ); + } + } + } +#endif + + private: + +// TrigoFletcherFunction(const TrigoFletcherFunction & ) {} +// TrigoFletcherFunction & operator=(const TrigoFletcherFunction &) { return *this; } + + double DoEval (const double * x) const { + gNCall++; + + + for (unsigned int i = 0; i < fDim ; ++i) { + cx [i] = std::cos(x[i]); + sx [i] = std::sin(x[i]); + } + + v = A*sx +B*cx; + r = v0-v; + + return r * r; + } + + + double DoDerivative (const double * x, unsigned int i ) const { + std::vector g(fDim); + Gradient(x,&g[0]); + return g[i]; + } + +private: + + unsigned int fDim; + + TMatrixD A; + TMatrixD B; + TVectorD x0; + mutable TVectorD sx0; + mutable TVectorD cx0; + mutable TVectorD sx; + mutable TVectorD cx; + mutable TVectorD v0; + mutable TVectorD v; + mutable TVectorD r; + + +}; + + +// CHEBYQUAD FUNCTION + +class ChebyQuadFunction : public ROOT::Math::IMultiGradFunction { + +public : + + ChebyQuadFunction(unsigned int n) : + fDim(n), + fvec(std::vector(n) ), + fTrueMin(std::vector(n) ) + { + } + + unsigned int NDim() const { return fDim; } + + ROOT::Math::IMultiGenFunction * Clone() const { + return new ChebyQuadFunction(*this); + } + + const double * TrueMinimum() const { + return &fTrueMin[0]; + } + + // use equally spaced points + void StartPoints(double * x, double * s) { + for (unsigned int i = 0; i < fDim; ++i) { + s[i] = 0.01; + x[i] = double(i)/(double(fDim)+1.0); + } + } + + // compute gradient + + void Gradient(const double * x, double * g) const { + gNCall2++; + unsigned int n = fDim; + // estimate first the fvec + DoCalculatefi(x); + + for (unsigned int j = 0; j < n; ++j) { + g[j] = 0.0; + double t1 = 1.0; + double t2 = 2.0 * x[j] - 1.0; + double t = 2.0 * t2; + double s1 = 0.0; + double s2 = 2.0; + for (unsigned int i = 0; i < n; ++i) { + g[j] += fvec[i] * s2; + double th = 4.0 * t2 + t * s2 - s1; + s1 = s2; + s2 = th; + th = t * t2 - t1; + t1 = t2; + t2 = th; + } + g[j] = 2. * g[j] / double(n); + } + + + } + + private: + + double DoEval (const double * x) const { + + gNCall++; + DoCalculatefi(x); + double f = 0; + for (unsigned int i = 0; i < fDim; ++i) + f += fvec[i] * fvec[i]; + + return f; + + } + + double DoDerivative (const double * x, unsigned int i ) const { + std::vector g(fDim); + Gradient(x,&g[0]); + return g[i]; + } + + void DoCalculatefi(const double * x) const { + // calculate the i- element ; F(X) = Sum {fi] + unsigned int n = fDim; + for (unsigned int i = 0; i < n; ++i) + fvec[i] = 0; + + for (unsigned int j = 0; j < n; ++j) { + double t1 = 1.0; + double t2 = 2.0 * x[j] - 1.0; + double t = 2.0 * t2; + for (unsigned int i = 0; i < n; ++i) { + fvec[i] += t2; + double th = t * t2 - t1; + t1 = t2; + t2 = th; + } + } + + // sum with the integral (integral is zero for odd Cheb polynomial and = 1/(i**2 -1) for the even ones + for (unsigned int i = 1; i <= n; ++i) { + int l = i-1; + fvec[l] /= double ( n ); + if ( ( i % 2 ) == 0 ) { + fvec[l] += 1.0 / ( double ( i*i ) - 1.0 ); + } + } + } + + unsigned int fDim; + mutable std::vector fvec; + mutable std::vector fTrueMin; +}; + +//WOOD function (4 dim function) + +double WoodFunction(const double * par) { + gNCall++; + + const double w = par[0]; + const double x = par[1]; + const double y = par[2]; + const double z = par[3]; + + const double w1 = w-1; + const double x1 = x-1; + const double y1 = y-1; + const double z1 = z-1; + const double tmp1 = x-w*w; + const double tmp2 = z-y*y; + + double f = 100*tmp1*tmp1+w1*w1+90*tmp2*tmp2+y1*y1+10.1*(x1*x1+z1*z1)+19.8*x1*z1; + return f; +} + +//Powell Function (4 dim function) +double PowellFunction(const double * par) { + gNCall++; + + const double w = par[0]; + const double x = par[1]; + const double y = par[2]; + const double z = par[3]; + + const double tmp1 = w+10*x; + const double tmp2 = y-z; + const double tmp3 = x-2*y; + const double tmp4 = w-z; + + double f = tmp1*tmp1+5*tmp2*tmp2+tmp3*tmp3*tmp3*tmp3+10*tmp4*tmp4*tmp4*tmp4; + return f; +} + +double SimpleQuadFunction(const double * par) { + gNCall++; + double x = par[0]; + double y = par[1]; + double f = x * x + 2 * y * y - x*y; + std::cout << "Quadfunc " << gNCall << "\t" << x << " , " << y << " f = " << f << std::endl; + return f; +} + +const double * TrueMinimum(const ROOT::Math::IMultiGenFunction & func) { + + const RosenBrockFunction * fRB = dynamic_cast< const RosenBrockFunction *> (&func); + if (fRB != 0) return fRB->TrueMinimum(); + const TrigoFletcherFunction * fTF = dynamic_cast< const TrigoFletcherFunction *> (&func); + if (fTF != 0) return fTF->TrueMinimum(); +// const ChebyQuadFunction * fCQ = dynamic_cast< const ChebyQuadFunction *> (&func); +// if (fCQ != 0) return fCQ->TrueMinimum(); + return 0; +} + +void printMinimum(const std::vector & x) { + std::cout << "Minimum X values\n"; + std::cout << "\t"; + int pr = std::cout.precision(12); + unsigned int n = x.size(); + for (unsigned int i = 0; i < n; ++i) { + std::cout << x[i]; + if ( i != n-1 ) std::cout << " , "; + if ( i > 0 && i % 6 == 0 ) std::cout << "\n\t"; + } + std::cout << std::endl; + std::cout.precision(pr); +} + +int DoNewMinimization( const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * s0, ROOT::Math::Minimizer * min, double &minval, double &edm, double * minx) { + + int iret = 0; + + min->SetMaxFunctionCalls(1000000); + min->SetMaxIterations(100000); + min->SetTolerance(gAbsTolerance); + if (func.NDim() >= 10) { + min->SetTolerance(0.01); + + } + + min->SetPrintLevel(gVerbose); + // check if func provides gradient + const ROOT::Math::IMultiGradFunction * gfunc = dynamic_cast(&func); + if (gfunc != 0 && useGradient) + min->SetFunction(*gfunc); + else + min->SetFunction(func); + + for (unsigned int i = 0; i < func.NDim(); ++i) { +#ifdef DEBUG + std::cout << "set variable " << i << " to value " << x0[i] << std::endl; +#endif + min->SetVariable(i,"x" + ROOT::Math::Util::ToString(i),x0[i], s0[i]); + } + + bool ret = min->Minimize(); + minval = min->MinValue(); + edm = min->Edm(); + + if (!ret) { + return -1; + } + + const double * xmin = min->X(); + + bool ok = true; + const double * trueMin = TrueMinimum(func); + if (trueMin != 0) { + for (unsigned int i = 0; i < func.NDim(); ++i) + ok &= (std::fabs(xmin[i]-trueMin[i] ) < gAbsTolerance); + } + + if (!ok) iret = -2; + int ncall_migrad = gNCall; + + // test Minos (use the default up of 1) + if (minos) { + + double el,eu; + for (unsigned int i = 0; i < func.NDim(); ++i) { + ret = min->GetMinosError(i,el,eu); + std::cout << "ncalls " << gNCall << "\t"; + if (ret) std::cout << "MINOS error for " << i << " = " << el << " " << eu << std::endl; + else std::cout << "MINOS failed for " << i << std::endl; + } + } + + std::cout << "\nMigrad Ncalls:\t " << ncall_migrad << std::endl; + std::cout << "Minos Ncalls :\t " << gNCall - ncall_migrad << std::endl; + + +// std::cout << "function at the minimum " << func(xmin) << std::endl; + std::copy(xmin,xmin+func.NDim(),minx); + min->Clear(); + + return iret; +} + +int DoOldMinimization( FCN func, TVirtualFitter * min, double &minval, double &edm) { + + int iret = 0; + + assert(min != 0); + min->SetFCN( func ); + + Double_t arglist[100]; + arglist[0] = gVerbose-1; + min->ExecuteCommand("SET PRINT",arglist,1); + + min->SetParameter(0,"x",-1.2,0.01,0,0); + min->SetParameter(1,"y", 1.0,0.01,0,0); + + arglist[0] = 0; + min->ExecuteCommand("SET NOW",arglist,0); + arglist[0] = 1000; // number of function calls + arglist[1] = gAbsTolerance; // tolerance + //min->ExecuteCommand("MIGRAD",arglist,0); + min->ExecuteCommand("MIGRAD",arglist,2); + + if (minos) min->ExecuteCommand("MINOS",arglist,0); + + Double_t parx,pary; + Double_t we,al,bl; + Char_t parName[32]; + min->GetParameter(0,parName,parx,we,al,bl); + min->GetParameter(1,parName,pary,we,al,bl); + + bool ok = ( TMath::Abs(parx-1.) < gAbsTolerance && + TMath::Abs(pary-1.) < gAbsTolerance ); + + + double errdef = 0; + int nvpar, nparx; + min->GetStats(minval,edm,errdef,nvpar,nparx); + if (!ok) iret = -2; + + min->Clear(); // for further use + return iret; + +} + + +int testNewMinimizer( const ROOT::Math::IMultiGenFunction & func, const double * x0, const double * s0, const std::string & minimizer, const std::string & algoType) { + + std::cout << "\n************************************************************\n"; + std::cout << "\tTest new ROOT::Math::Minimizer\n"; + std::cout << "\tMinimizer is " << minimizer << " " << algoType << std::endl; + + int iret = 0; + double minval,edm = 0; + std::vector xmin(func.NDim() ); + + TStopwatch w; + w.Start(); + + ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer(minimizer, algoType); + if (min == 0) { + std::cout << "Error using minimizer " << minimizer << " " << algoType << std::endl; + return -1; + } + + + for (int i = 0; i < gNmin; ++i) { + gNCall = 0; gNCall2 = 0; + iret |= DoNewMinimization(func, x0, s0, min,minval,edm,&xmin[0]); + } + + w.Stop(); + if (iret != 0) std::cout << "\n****** ERROR: Minimization FAILED ! \n"; + int pr = std::cout.precision(18); + std::cout << "\nNCalls: \t" << gNCall << " , " << gNCall2 + << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr); + std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl; + printMinimum(xmin ); + std::cout << "\n************************************************************\n"; + +#ifdef CHECK_WITHMINUIT + // do Minuit after BFGS + if (minimizer == "GSL_BFGS") { + std::cout << "DO cmaes from last point\n"; + gNCall = 0; + iret |= DoNewMinimization(func, &xmin.front(), s0, "cmaes","",minval,edm,&xmin[0]); + int pr = std::cout.precision(18); + std::cout << "\nNCalls: \t" << gNCall << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr); + std::cout << std::endl; + } +#endif + + delete min; + + return iret; +} + + +int testOldMinimizer( FCN func, const std::string & fitter, int n=25) { + + std::cout << "\n************************************************************\n"; + std::cout << "\tTest using TVirtualFitter\n"; + std::cout << "\tFitter is " << fitter << std::endl; + + int iret = 0; + double minval,edm = 0; + + TStopwatch w; + w.Start(); + + TVirtualFitter::SetDefaultFitter(fitter.c_str()); + + TVirtualFitter *min = TVirtualFitter::Fitter(0,n); + + //min->Dump(); + + for (int i = 0; i < gNmin; ++i) { + gNCall = 0; + iret |= DoOldMinimization(func, min,minval,edm); + } + + w.Stop(); + if (iret != 0) std::cout << "\n****** ERROR: Minimization FAILED ! \n"; + int pr = std::cout.precision(18); + std::cout << "\nNCalls: \t" << gNCall << "\tMinValue: \t" << minval << "\tEdm: \t" << edm; std::cout.precision(pr); + std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl; + std::cout << "\n************************************************************\n"; + + return iret; +} + +int testRosenBrock() { + + int iret = 0; + + + std::cout << "\n************************************************************\n"; + std::cout << "\tROSENBROCK function test\n\n"; + + double s0[2] = {0.01,0.01}; + + // minimize using Rosenbrock Function +#ifndef DEBUG + gNmin = 2000; +#endif + + gNmin = 1; + + + RosenBrockFunction fRB; + double xRB[2] = { -1.,1.2}; + + iret |= testNewMinimizer(fRB,xRB,s0,"Minuit",""); + iret |= testNewMinimizer(fRB,xRB,s0,"cmaes",""); + +// iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","ConjugateFR"); +// iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","ConjugatePR"); +// iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","BFGS"); +// iret |= testNewMinimizer(fRB,xRB,s0,"GSLMultiMin","BFGS2"); + +// iret |= testOldMinimizer(RosenBrock,"Minuit",2); +// iret |= testOldMinimizer(RosenBrock,"cmaes",2); + + + return iret; +} + + +int testChebyQuad() { + + int iret = 0; + + // test with ChebyQuad function + + const int n = 8; + ChebyQuadFunction func(n); + +#ifndef DEBUG + gNmin = std::max(1, int(20000/n/n) ); +#endif + gNmin = 1; + + + double s0[n]; + double x0[n]; + func.StartPoints(x0,s0); + + std::cout << "\n************************************************************\n"; + std::cout << "\tCHEBYQUAD function test , n = " << n << "\n\n"; + + +// double x[8] = {0.043153E+00, 0.193091E+00, 0.266329E+00, 0.500000E+00, +// 0.500000E+00, 0.733671E+00, 0.806910E+00, 0.956847E+00 }; +// double x[2] = {0.5, 0.5001}; +// std::cout << "FUNC " << func(x) << std::endl; + double x1[100] = { 0.00712780070646 , 0.0123441993113 , 0.0195428378255 , 0.0283679084192 , 0.0385291131289 , 0.0492202424892 , 0.0591277976178 , + 0.0689433195252 , 0.0791293590525 , 0.088794974369 , 0.0988949579193 , 0.108607151294 , 0.118571075831 , + 0.128605446238 , 0.137918291068 , 0.149177761352 , 0.156665324587 , 0.170851061982 , 0.174688134016 , + 0.192838903364 , 0.193078270803 , 0.209255377225 , 0.217740096779 , 0.225888518345 , 0.241031047421 , + 0.244253844041 , 0.257830449676 , 0.269467652526 , 0.274286498012 , 0.288877029988 , 0.297549406597 , + 0.304950954529 , 0.319230811642 , 0.326387092206 , 0.335229058731 , 0.349178359226 , 0.355905988048 , + 0.365197862755 , 0.379068092603 , 0.385826036925 , 0.394978252826 , 0.408974425717 , 0.415968185065 , + 0.424621041584 , 0.438837361714 , 0.446214149031 , 0.454242324351 , 0.468614308013 , 0.476506553416 , + 0.483916944941 , 0.498229247409 , 0.506794629616 , 0.513736742474 , 0.527712475478 , 0.537073277673 , + 0.543731917673 , 0.557187513963 , 0.567346279639 , 0.57379846397 , 0.586691058785 , 0.597561941009 , + 0.60382873461 , 0.616316037506 , 0.627719652101 , 0.633760038662 , 0.646175283836 , 0.657809344891 , + 0.663569004722 , 0.676314563639 , 0.687674566849 , 0.69332205923 , 0.706839545953 , 0.716907408637 , + 0.723407327715 , 0.738019389561 , 0.744806584048 , 0.754657613362 , 0.769181875619 , 0.772250323489 , + 0.787104833193 , 0.795856360905 , 0.804099304478 , 0.82142178741 , 0.819589601284 , 0.839024540481 , + 0.842457233039 , 0.857393475964 , 0.86408033345 , 0.876329840525 , 0.884867318008 , 0.895744532071 , + 0.905113958163 , 0.915445338697 , 0.925148068352 , 0.935344457785 , 0.945127838313 , 0.955272197168 , + 0.965687518559 , 0.975129521484 , 0.982662007764 }; + + std::cout << "FUNC " << func(x1) << std::endl; + + + iret |= testNewMinimizer(func,x0,s0, "Minuit",""); + iret |= testNewMinimizer(func,x0,s0, "cmaes",""); + +// iret |= testNewMinimizer(func,x0,s0, "GSLMultiMin","ConjugateFR"); +// iret |= testNewMinimizer(func,x0,s0, "GSLMultiMin","ConjugatePR"); +// iret |= testNewMinimizer(func,x0,s0, "GSLMultiMin","BFGS"); + + return iret; +} + + +int testTrigoFletcher() { + + int iret = 0; + + + // test with fletcher trigonometric function +#ifndef DEBUG + gNmin = 2; +#endif + gNmin = 1; + + const int nT = 50; + TrigoFletcherFunction fTrigo(nT); + double sTrigo[nT]; + double xTrigo[nT]; + fTrigo.StartPoints(xTrigo,sTrigo); + + std::cout << "\n************************************************************\n"; + std::cout << "\tTRIGONOMETRIC Fletcher function test , n = " << nT << "\n\n"; + + + iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"cmaes",""); + iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"Minuit",""); + +// iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","ConjugateFR"); +// iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","ConjugatePR"); +// iret |= testNewMinimizer(fTrigo,xTrigo,sTrigo,"GSLMultiMin","BFGS"); + + + return iret; +} + + +int testWood() { + + int iret = 0; + + + // test with Wood function (4d) +// minimum : F(1,1,1,1) = 0. + + +#ifndef DEBUG + gNmin = 1000; +#endif + gNmin = 1; + + ROOT::Math::Functor f(&WoodFunction,4); + + double x0[4] = { -3, -1, -3, -1 }; + double s0[4] = { 0.1, 0.1, 0.1, 0.1}; + + std::cout << "\n************************************************************\n"; + std::cout << "\tWOOD 4 function test \n\n"; + + + iret |= testNewMinimizer(f, x0, s0,"Minuit",""); + iret |= testNewMinimizer(f, x0, s0,"cmaes",""); + + + return iret; +} + +int testPowell() { + + int iret = 0; + + + // test with Powell function (4d) + // minimum is at F(0,0,0,0) = 0 +#ifndef DEBUG + gNmin = 1000; +#endif + gNmin = 1; + + ROOT::Math::Functor f(&PowellFunction,4); + + double x0[4] = { -3, -1, 0, 1 }; + double s0[4] = { 0.1, 0.1, 0.1, 0.1}; + + std::cout << "\n************************************************************\n"; + std::cout << "\tPOWELL function test \n\n"; + + + iret |= testNewMinimizer(f, x0, s0,"Minuit",""); + iret |= testNewMinimizer(f, x0, s0,"cmaes",""); + + + return iret; +} + + +int testQuadFunc() { + + int iret = 0; + + + // test with a simple quadratic function 2d + // minimum is at F(0,0) = 0 +#ifndef DEBUG + gNmin = 1000; +#endif + gNmin = 1; + + ROOT::Math::Functor f(&SimpleQuadFunction,2); + + double x0[4] = { -3, -3 }; + double s0[4] = { 0.1, 0.1}; + + std::cout << "\n************************************************************\n"; + std::cout << "\tSIMPLE QUAD function test \n\n"; + + + iret |= testNewMinimizer(f, x0, s0,"Minuit",""); + iret |= testNewMinimizer(f, x0, s0,"cmaes",""); + + + return iret; +} + + +int main() { + + gSystem->Load("/usr/lib/x86_64-linux-gnu/libglog.so"); + gSystem->Load("/usr/lib/x86_64-linux-gnu/libgflags.so"); + + int iret = 0; + +#ifdef DEBUG + gVerbose = 4; + gNmin = 1; +#endif + + + iret |= testRosenBrock(); + iret |= testChebyQuad(); + iret |= testTrigoFletcher(); + + iret |= testWood(); + iret |= testPowell(); + + iret |= testQuadFunc(); + + + + if (iret != 0) + std::cerr << "testMinim :\t FAILED " << std::endl; + else + std::cerr << "testMinim :\t OK " << std::endl; + return iret; + +} diff --git a/math/cmaes/test/testNdimFit.cxx b/math/cmaes/test/testNdimFit.cxx new file mode 100644 index 0000000000000..b3cecc526f12d --- /dev/null +++ b/math/cmaes/test/testNdimFit.cxx @@ -0,0 +1,285 @@ +#include "TMath.h" +#include "TSystem.h" +#include "TRandom3.h" +#include "TTree.h" +#include "TROOT.h" + +#include "Fit/UnBinData.h" +#include "Fit/Fitter.h" + + +#include "Math/IParamFunction.h" +#include "Math/WrappedTF1.h" +#include "Math/WrappedMultiTF1.h" +#include "Math/WrappedParamFunction.h" +#include "Math/MultiDimParamFunctionAdapter.h" + +#include "TGraphErrors.h" + +#include "TStyle.h" + + +#include "Math/DistFunc.h" + +#include +#include + +#include "TStopwatch.h" + +#define DEBUG + +// test fit with many dimension + +const int N = 10; +const std::string branchType = "x[10]/D"; +const int NPoints = 100000; + +// const int N = 50; +// const std::string branchType = "x[50]/D"; +// const int NPoints = 10000; + + + +double truePar[2*N]; +double iniPar[2*N]; +//const int nfit = 1; +const int strategy = 0; + +double gausnorm(const double *x, const double *p) { + + double invsig = 1./p[1]; + double tmp = (x[0]-p[0]) * invsig; + const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 ); + return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig; +} + +double gausnormN(const double *x, const double *p) { + double f = 1; + for (int i = 0; i < N; ++i) + f *= gausnorm(x+i,p+2*i); + + return f; +} + +struct CMAES { + static std::string name() { return "cmaes"; } + static std::string name2() { return "acmaes"; } +}; + +// fill fit data structure +ROOT::Fit::UnBinData * FillUnBinData(TTree * tree ) { + + // fill the unbin data set from a TTree + ROOT::Fit::UnBinData * d = 0; + + // for the large tree access directly the branches + d = new ROOT::Fit::UnBinData(); + + unsigned int n = tree->GetEntries(); +#ifdef DEBUG + std::cout << "number of unbin data is " << n << " of dim " << N << std::endl; +#endif + d->Initialize(n,N); + TBranch * bx = tree->GetBranch("x"); + double vx[N]; + bx->SetAddress(vx); + std::vector m(N); + for (int unsigned i = 0; i < n; ++i) { + bx->GetEntry(i); + d->Add(vx); + for (int j = 0; j < N; ++j) + m[j] += vx[j]; + } + +#ifdef DEBUG + std::cout << "average values of means :\n"; + for (int j = 0; j < N; ++j) + std::cout << m[j]/n << " "; + std::cout << "\n"; +#endif + + delete tree; + tree = 0; + return d; + +} + + +// unbin fit + +typedef ROOT::Math::IParamMultiFunction Func; +template +int DoUnBinFit(T * tree, Func & func, bool debug = false ) { + + ROOT::Fit::UnBinData * d = FillUnBinData(tree ); + // need to have done Tree->Draw() before fit + //FillUnBinData(d,tree); + + //std::cout << "data size type and size is " << typeid(*d).name() << " " << d->Size() << std::endl; + std::cout << "Fit data size = " << d->Size() << " dimension = " << d->NDim() << std::endl; + + + + //printData(d); + + // create the fitter + //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl; + + ROOT::Fit::Fitter fitter; + fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str()); + + if (debug) + fitter.Config().MinimizerOptions().SetPrintLevel(3); + else + fitter.Config().MinimizerOptions().SetPrintLevel(0); + + // set tolerance 1 for tree to be same as in TTTreePlayer::UnBinFIt + fitter.Config().MinimizerOptions().SetTolerance(1); + + // set strategy (0 to avoid MnHesse + fitter.Config().MinimizerOptions().SetStrategy(strategy); + + + // create the function + + fitter.SetFunction(func); + // need to fix param 0 , normalization in the unbinned fits + //fitter.Config().ParSettings(0).Fix(); + + bool ret = fitter.Fit(*d); + if (!ret) { + std::cout << " Fit Failed " << std::endl; + return -1; + } + if (debug) + fitter.Result().Print(std::cout); + + // check fit result + double chi2 = 0; + for (int i = 0; i < N; ++i) { + double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) ); + chi2 += d*d; + } + double prob = ROOT::Math::chisquared_cdf_c(chi2,N); + int iret = (prob < 1.0E-6) ? -1 : 0; + if (iret != 0) { + std::cout <<"Found difference in fitted values - prob = " << prob << std::endl; + if (!debug) fitter.Result().Print(std::cout); + for (int i = 0; i < N; ++i) { + double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) ); + std::cout << "par_" << i << " = " << fitter.Result().Value(i) << " true = " << truePar[i] << " pull = " << d << std::endl; + } + + } + + delete d; + + return iret; + +} + + +template +int DoFit(TTree * tree, Func & func, bool debug = false ) { + return DoUnBinFit(tree, func, debug); +} +// template +// int DoFit(TH1 * h1, Func & func, bool debug = false, bool copyData = false ) { +// return DoBinFit(h1, func, debug, copyData); +// } +// template +// int DoFit(TGraph * gr, Func & func, bool debug = false, bool copyData = false ) { +// return DoBinFit(gr, func, debug, copyData); +// } + +template +int FitUsingNewFitter(FitObj * fitobj, Func & func ) { + + std::cout << "\n************************************************************\n"; + std::cout << "\tFit using new Fit::Fitter " << typeid(*fitobj).name() << std::endl; + std::cout << "\tMinimizer is " << MinType::name() << " " << MinType::name2() << " func dim = " << func.NDim() << std::endl; + + int iret = 0; + TStopwatch w; w.Start(); + +#ifdef DEBUG + std::cout << "initial Parameters " << iniPar << " " << *iniPar << " " << *(iniPar+1) << std::endl; + func.SetParameters(iniPar); + iret |= DoFit(fitobj,func,true ); + +#else + for (int i = 0; i < nfit; ++i) { + func.SetParameters(iniPar); + iret = DoFit(fitobj,func, false); + if (iret != 0) { + std::cout << "Fit failed " << std::endl; + break; + } + } +#endif + w.Stop(); + std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl; + std::cout << "\n************************************************************\n"; + + return iret; +} + + +int testNdimFit() { + + + std::cout << "\n\n************************************************************\n"; + std::cout << "\t UNBINNED TREE (GAUSSIAN MULTI-DIM) FIT\n"; + std::cout << "************************************************************\n"; + + TTree * t1 = new TTree("t2","a large Tree with gaussian variables"); + double x[N]; + Int_t ev; + t1->Branch("x",x,branchType.c_str()); + t1->Branch("ev",&ev,"ev/I"); + + // generate the true parameters + for (int j = 0; j < N; ++j) { + double mu = double(j)/10.; + double s = 1.0 + double(j)/10.; + truePar[2*j] = mu; + truePar[2*j+1] = s; + } + + + //fill the tree + TRandom3 r; + for (Int_t i=0;iFill(); + + } + //t1.Draw("x"); // to select fit variable + + + for (int i = 0; i f2(&gausnormN,N,2*N,iniPar); + + int iret = 0; + iret |= FitUsingNewFitter(t1,f2); + + return iret; +} + +int main() { + gSystem->Load("/usr/lib/x86_64-linux-gnu/libglog.so"); + gSystem->Load("/usr/lib/x86_64-linux-gnu/libgflags.so"); + return testNdimFit(); +} diff --git a/math/cmaes/test/testUnbinGausFit.cxx b/math/cmaes/test/testUnbinGausFit.cxx new file mode 100644 index 0000000000000..9192ac50dc8fa --- /dev/null +++ b/math/cmaes/test/testUnbinGausFit.cxx @@ -0,0 +1,333 @@ +#include "TMath.h" +#include "TSystem.h" +#include "TRandom3.h" +#include "TTree.h" +#include "TROOT.h" + +#include "Fit/UnBinData.h" +#include "Fit/Fitter.h" + + +#include "Math/IParamFunction.h" +#include "Math/WrappedTF1.h" +#include "Math/WrappedMultiTF1.h" +#include "Math/WrappedParamFunction.h" +#include "Math/MultiDimParamFunctionAdapter.h" + +#include "TGraphErrors.h" + +#include "TStyle.h" + + +#include "Math/DistFunc.h" + +#include +#include + +#include "TStopwatch.h" + +#define DEBUG + +// test fit with many dimension + +const int N = 1; // 1d fit +const int NGaus = 3; +const int NPar = 8; // sum of 3 gaussians +const std::string branchType = "x[1]/D"; + +// for 8 core testing use 1M points +const int NPoints = 100000; +double truePar[NPar]; +double iniPar[NPar]; +//const int nfit = 1; +const int strategy = 0; + +double gausnorm(const double *x, const double *p) { + + double invsig = 1./std::abs(p[1]); + double tmp = (x[0]-p[0]) * invsig; + const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 ); + return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig; +} + +double gausSum(const double *x, const double *p) { + + double f = gausnorm(x,p+2) + + p[0] * gausnorm(x,p+4) + + p[1] * gausnorm(x,p+6); + + double norm = 1. + p[0] + p[1]; + return f/norm; +} + +struct CMAES { + static std::string name() { return "cmaes"; } + static std::string name2() { return "acmaes"; } +}; + +// fill fit data structure +ROOT::Fit::UnBinData * FillUnBinData(TTree * tree ) { + + // fill the unbin data set from a TTree + ROOT::Fit::UnBinData * d = 0; + + // for the large tree access directly the branches + d = new ROOT::Fit::UnBinData(); + + unsigned int n = tree->GetEntries(); +#ifdef DEBUG + std::cout << "number of unbin data is " << n << " of dim " << N << std::endl; +#endif + d->Initialize(n,N); + TBranch * bx = tree->GetBranch("x"); + double vx[N]; + bx->SetAddress(vx); + std::vector m(N); + for (int unsigned i = 0; i < n; ++i) { + bx->GetEntry(i); + d->Add(vx); + for (int j = 0; j < N; ++j) + m[j] += vx[j]; + } + +#ifdef DEBUG + std::cout << "average values of means :\n"; + for (int j = 0; j < N; ++j) + std::cout << m[j]/n << " "; + std::cout << "\n"; +#endif + + return d; +} + + +// unbin fit + +typedef ROOT::Math::IParamMultiFunction Func; +template +int DoUnBinFit(T * tree, Func & func, bool debug = false ) { + + ROOT::Fit::UnBinData * d = FillUnBinData(tree ); + // need to have done Tree->Draw() before fit + //FillUnBinData(d,tree); + + std::cout << "Filled the fit data " << std::endl; + //printData(d); + +#ifdef DEBUG + std::cout << "data size type and size is " << typeid(*d).name() << " " << d->Size() << std::endl; +#endif + + + + // create the fitter + //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl; + + ROOT::Fit::Fitter fitter; + fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str()); + + if (debug) + fitter.Config().MinimizerOptions().SetPrintLevel(3); + else + fitter.Config().MinimizerOptions().SetPrintLevel(1); + + + // set tolerance 1 for tree to be same as in TTTreePlayer::UnBinFIt + fitter.Config().MinimizerOptions().SetTolerance(0.01); + + // set strategy (0 to avoid MnHesse + fitter.Config().MinimizerOptions().SetStrategy(strategy); + + + // create the function + + fitter.SetFunction(func); + // need to set limits to constant term + fitter.Config().ParSettings(0).SetLowerLimit(0.); + fitter.Config().ParSettings(1).SetLowerLimit(0.); + + if (debug) + std::cout << "do fitting... " << std::endl; + + bool ret = fitter.Fit(*d); + if (!ret) { + std::cout << " Fit Failed " << std::endl; + return -1; + } + if (debug) + fitter.Result().Print(std::cout); + + // check fit result + double chi2 = 0; + //if (fitter.Result().Value(0) < 0.5 ) { + for (int i = 0; i < NPar; ++i) { + double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) ); + chi2 += d*d; + } +//} +// else { +// double truePar2[NPar]; +// truePar2[0] = 1.-truePar[0]; +// truePar2[1] = truePar[3]; +// truePar2[2] = truePar[4]; +// truePar2[3] = truePar[1]; +// truePar2[4] = truePar[2]; +// for (int i = 0; i < N; ++i) { +// double d = ( truePar2[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) ); +// chi2 += d*d; +// } +// } + double prob = ROOT::Math::chisquared_cdf_c(chi2,NPar); + int iret = (prob < 1.0E-6) ? -1 : 0; + if (iret != 0) { + std::cout <<"Found difference in fitted values - chi2 = " << chi2 + << " prob = " << prob << std::endl; + fitter.Result().Print(std::cout); + } + + delete d; + + return iret; + +} + + +template +int DoFit(TTree * tree, Func & func, bool debug = false ) { + return DoUnBinFit(tree, func, debug); +} +// template +// int DoFit(TH1 * h1, Func & func, bool debug = false, bool copyData = false ) { +// return DoBinFit(h1, func, debug, copyData); +// } +// template +// int DoFit(TGraph * gr, Func & func, bool debug = false, bool copyData = false ) { +// return DoBinFit(gr, func, debug, copyData); +// } + +template +int FitUsingNewFitter(FitObj * fitobj, Func & func ) { + + std::cout << "\n************************************************************\n"; + std::cout << "\tFit using new Fit::Fitter " << typeid(*fitobj).name() << std::endl; + std::cout << "\tMinimizer is " << MinType::name() << " " << MinType::name2() << " func dim = " << func.NDim() << std::endl; + + int iret = 0; + TStopwatch w; w.Start(); + +#ifdef DEBUG + std::cout << "initial Parameters " << iniPar << " " << *iniPar << " " << *(iniPar+1) << std::endl; + func.SetParameters(iniPar); + iret |= DoFit(fitobj,func,true ); + if (iret != 0) { + std::cout << "Test failed " << std::endl; + } + +#else + for (int i = 0; i < nfit; ++i) { + func.SetParameters(iniPar); + iret = DoFit(fitobj,func, false); + if (iret != 0) { + std::cout << "Test failed " << std::endl; + break; + } + } +#endif + w.Stop(); + std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl; + std::cout << "\n************************************************************\n"; + + return iret; +} + + +int testNdimFit() { + + + std::cout << "\n\n************************************************************\n"; + std::cout << "\t UNBINNED TREE (GAUSSIAN MULTI-DIM) FIT\n"; + std::cout << "************************************************************\n"; + + TTree t1("t2","a large Tree with gaussian variables"); + double x[N]; + Int_t ev; + t1.Branch("x",x,branchType.c_str()); + t1.Branch("ev",&ev,"ev/I"); + + // generate the true parameters +// for (int j = 0; j < NGaus; ++j) { +// double a = j+1; +// double mu = double(j)/NGaus; +// double s = 1.0 + double(j)/NGaus; +// truePar[3*j] = a; +// truePar[3*j+1] = mu; +// truePar[3*j+2] = s; +// tot += a; +// } + truePar[0] = 0.2; // % second gaussian + truePar[1] = 0.05; // % third gaussian ampl + truePar[2] = 0.; // mean first gaussian + truePar[3] = 0.5; // s1 + truePar[4] = 0.; // mean secon gauss + truePar[5] = 1; + truePar[6] = -3; // mean third gaus + truePar[7] = 10; + + + + //fill the tree + TRandom3 r; + double norm = (1+truePar[0] + truePar[1] ); + double a = 1./norm; + double b = truePar[0]/ norm; + double c = truePar[1]/ norm; + assert(a+b+c == 1.); + std::cout << " True amplitude gaussians " << a << " " << b << " " << c << std::endl; + for (Int_t i=0;i f2(&gausSum,1,NPar,iniPar); + + int iret = 0; + iret |= FitUsingNewFitter(&t1,f2); + + return iret; +} + +int main() { + gSystem->Load("/usr/lib/x86_64-linux-gnu/libglog.so"); + gSystem->Load("/usr/lib/x86_64-linux-gnu/libgflags.so"); + return testNdimFit(); +} diff --git a/math/cmaes/test/testUserFunc.cxx b/math/cmaes/test/testUserFunc.cxx new file mode 100644 index 0000000000000..d175326201394 --- /dev/null +++ b/math/cmaes/test/testUserFunc.cxx @@ -0,0 +1,88 @@ +// @(#)root/minuit2:$Id$ +// Author: L. Moneta 10/2005 + +/********************************************************************** + * * + * Copyright (c) 2005 ROOT Foundation, CERN/PH-SFT * + * * + **********************************************************************/ + + +#include "TApplication.h" +#include "TH1.h" +#include "TF1.h" +#include "TRandom3.h" +#include "Math/MinimizerOptions.h" +#include "TMath.h" +#include "TSystem.h" + +#include + +double myfunc( double * x, double * p) { + + return p[0]*TMath::Gaus(x[0],p[1],p[2]); +} + +void testUserFunc(std::string type="cmaes", int n = 1000) { + + + + gRandom = new TRandom3(); + + + ROOT::Math::MinimizerOptions::SetDefaultMinimizer(type.c_str() ); + + + + TH1D * h1 = new TH1D("h1","fit histo 1",100, -5, 5. ); + +// gStyle->SetOptStat(1111111); +// gStyle->SetOptFit(1111111); + + + + + for (int i = 0; i < n; ++i) { + h1->Fill( gRandom->Gaus(0,1) ); + } + + TF1 * f = new TF1("f",myfunc,-10,10,3); + double p[3] = { 100.0, 0.0, 1.0 } ; + f->SetParameters(p); + + h1->Fit(f); + // try fix a parameter + //TVirtualFitter * fitter = TVirtualFitter::GetFitter(); + //std::cout << typeid(*fitter).name() << std::endl; + //fitter->FixParameter(2); + f->FixParameter(2,1.0); + + h1->Fit(f,"V"); + + h1->Draw(); + + + +} + +#ifndef __CINT__ +int main(int argc, char **argv) +{ + gSystem->Load("/usr/lib/x86_64-linux-gnu/libglog.so"); + gSystem->Load("/usr/lib/x86_64-linux-gnu/libgflags.so"); + if (argc > 1) { + TApplication theApp("App", &argc, argv); + testUserFunc( ); + theApp.Run(); + } + else + testUserFunc( ); + return 0; +} +#endif + +//#ifndef __CINT__ +//int main() { +// testUserFunc( ); +//} +//#endif diff --git a/math/mathcore/src/Factory.cxx b/math/mathcore/src/Factory.cxx index f31f4ba8c826c..f7c37141e4ca8 100644 --- a/math/mathcore/src/Factory.cxx +++ b/math/mathcore/src/Factory.cxx @@ -26,6 +26,7 @@ // #define MATH_NO_PLUGIN_MANAGER // #define HAS_MINUIT // #define HAS_MINUIT2 +// #define HAS_CMAES #ifndef MATH_NO_PLUGIN_MANAGER // use ROOT Plug-in manager @@ -41,6 +42,9 @@ #ifdef HAS_MINUIT #include "TMinuitMinimizer.h" #endif +#ifdef HAS_CMAES +#include "TCMAESMinimizer.h" +#endif #ifdef R__HAS_MATHMORE #include "Math/GSLMinimizer.h" #include "Math/GSLNLSMinimizer.h" @@ -52,7 +56,7 @@ #include #include -//#define DEBUG +#define DEBUG #ifdef DEBUG #include #endif @@ -80,6 +84,12 @@ ROOT::Math::Minimizer * ROOT::Math::Factory::CreateMinimizer(const std::string & s1 = "Minuit"; minim = s1.c_str(); } + if (minimizerType.find("cmaes")!=std::string::npos + ||minimizerType.find("ipop")!=std::string::npos) + { + s1 = minimizerType; + minim = s1.c_str(); + } if (minimizerType.empty() ) minim = ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str(); @@ -104,9 +114,9 @@ ROOT::Math::Minimizer * ROOT::Math::Factory::CreateMinimizer(const std::string & else std::cout << "Error creating Minimizer " << minimizerType << " " << algoType << std::endl; #endif - return min; } + std::cout << "returning 0\n"; return 0; } @@ -140,6 +150,12 @@ ROOT::Math::Minimizer * ROOT::Math::Factory::CreateMinimizer(const std::string & min = new TMinuitMinimizer(algoType.c_str()); #endif +#ifdef HAS_CMAES + if (minimizerType.find("cmaes") != std::string::npos + || minimizerType.find("ipop") != std::string::npos) + min = new ROOT::cmaes::TCMAESMinimizer(algoType.c_str()); +#endif + #ifdef R__HAS_MATHMORE // use GSL minimizer if (minimizerType == "GSL") @@ -192,5 +208,3 @@ ROOT::Math::DistSampler * ROOT::Math::Factory::CreateDistSampler(const std::stri return 0; #endif } - - diff --git a/math/mathcore/src/FitConfig.cxx b/math/mathcore/src/FitConfig.cxx index 7938e355a4d5e..1a446fbd28cfb 100644 --- a/math/mathcore/src/FitConfig.cxx +++ b/math/mathcore/src/FitConfig.cxx @@ -157,8 +157,8 @@ void FitConfig::SetParamsSettings(unsigned int npar, const double *params, const if (val == 0) step = 0.3; } else - step = vstep[i]; - + step = vstep[i]; + if (createNew) fSettings.push_back( ParameterSettings("Par_" + ROOT::Math::Util::ToString(i), val, step ) ); else { diff --git a/math/mathcore/src/MinimizerOptions.cxx b/math/mathcore/src/MinimizerOptions.cxx index 67c7a280523d4..40ee5560dc8f1 100644 --- a/math/mathcore/src/MinimizerOptions.cxx +++ b/math/mathcore/src/MinimizerOptions.cxx @@ -178,11 +178,24 @@ void MinimizerOptions::ResetToDefaultOptions() { fMinimType = "Minuit2"; fAlgoType = "Fumili"; } + else if (fMinimType.find("cmaes")!=std::string::npos + || fMinimType.find("ipop")!=std::string::npos) + { + fAlgoType = fMinimType; + fMinimType = "cmaes"; + } else if (fMinimType == "GSLMultiMin" && fAlgoType == "Migrad") fAlgoType = "BFGS2"; - + else if (fMinimType.find("cmaes")!=std::string::npos + || fMinimType.find("ipop")!=std::string::npos) + { + fAlgoType = fMinimType; + fMinimType = "cmaes"; + } + delete fExtraOptions; fExtraOptions = 0; + // check if extra options exists (copy them if needed) if (Minim::gDefaultExtraOptions) fExtraOptions = Minim::gDefaultExtraOptions->Clone(); diff --git a/tutorials/fit/Ifit.C b/tutorials/fit/Ifit.C index 9eae64fa765cf..c93b3bd95f18f 100644 --- a/tutorials/fit/Ifit.C +++ b/tutorials/fit/Ifit.C @@ -36,8 +36,8 @@ Double_t func(float x,float y,Double_t *par) //______________________________________________________________________________ void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { - const Int_t nbins = 5; - Int_t i; + const Int_t nbins = 5; + Int_t i; //calculate chisquare Double_t chisq = 0; diff --git a/tutorials/fit/cmaesFitBench.C b/tutorials/fit/cmaesFitBench.C new file mode 100644 index 0000000000000..fc070af89b756 --- /dev/null +++ b/tutorials/fit/cmaesFitBench.C @@ -0,0 +1,131 @@ +//+ Fitting 1-D histograms with cmaes +// Author: L. Moneta 10/2005 +// Author: E. Benazera 05/2014 + +/********************************************************************** + * * + * Copyright (c) 2014 INRIA * + * Copyright (c) 2005 ROOT Foundation, CERN/PH-SFT * + * * + **********************************************************************/ + +#include "TH1.h" +#include "TF1.h" +#include "TCanvas.h" +#include "TStopwatch.h" +#include "TSystem.h" +#include "TRandom3.h" +#include "TVirtualFitter.h" +#include "TPaveLabel.h" +#include "TStyle.h" +#include "TMath.h" +#include "TROOT.h" +#include "TFrame.h" +#include "TSystem.h" + +//#include "Fit/FitConfig.h" + + +TF1 *fitFcn; +TH1 *histo; + +// Quadratic background function +Double_t background(Double_t *x, Double_t *par) { + return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; +} + +// Lorenzian Peak function +Double_t lorentzianPeak(Double_t *x, Double_t *par) { + return (0.5*par[0]*par[1]/TMath::Pi()) / + TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]); +} + +// Sum of background and peak function +Double_t fitFunction(Double_t *x, Double_t *par) { + return background(x,par) + lorentzianPeak(x,&par[3]); +} + +void DoFit(const char* fitter, TVirtualPad *pad, Int_t npass, Double_t sigma, Int_t lambda) { + printf("\n*********************************************************************************\n"); + printf("\t %s \n",fitter); + printf("*********************************************************************************\n"); + + gRandom = new TRandom3(); + TStopwatch timer; + // timer.Start(); + TVirtualFitter::SetDefaultFitter(fitter); + //ROOT::Fit::FitConfig::SetDefaultMinimizer(fitter); + ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default(fitter); + opts.SetIntValue("lambda",lambda); + opts.SetIntValue("lscaling",1); + opts.SetRealValue("sigma",sigma); + + pad->SetGrid(); + pad->SetLogy(); + fitFcn->SetParameters(1,1,1,6,.03,1); + fitFcn->Update(); + std::string title = std::string(fitter) + " fit bench"; + histo = new TH1D(fitter,title.c_str(),200,0,3); + + timer.Start(); + for (Int_t pass=0;passSetParameters(1,1,1,6,.03,1); + for (Int_t i=0;i<5000;i++) { + histo->Fill(fitFcn->GetRandom()); + } + //histo->Print("all"); + //histo->Fit(fitFcn,"Q0"); // from TH1.cxx: Q: quiet, 0: do not plot + } + + histo->Fit(fitFcn,"V"); // E: use Minos, V: verbose. + timer.Stop(); + + (histo->GetFunction("fitFcn"))->SetLineColor(kRed+3); + gPad->SetFillColor(kYellow-10); + + + Double_t cputime = timer.CpuTime(); + printf("%s, npass=%d : RT=%7.3f s, Cpu=%7.3f s\n",fitter,npass,timer.RealTime(),cputime); + TPaveLabel *p = new TPaveLabel(0.45,0.7,0.88,0.8,Form("%s CPU= %g s",fitter,cputime),"brNDC"); + p->Draw(); + p->SetTextColor(kRed+3); + p->SetFillColor(kYellow-8); + pad->Update(); + + //delete p; + //delete histo; + delete gRandom; +} + +void cmaesFitBench(Int_t npass=20, Double_t sigma=0.1, Int_t lambda=-1) { + TH1::AddDirectory(kFALSE); + TCanvas *c1 = new TCanvas("FitBench","Fitting Demo",10,10,900,900); + c1->Divide(2,2); + c1->SetFillColor(kYellow-9); + // create a TF1 with the range from 0 to 3 and 6 parameters + fitFcn = new TF1("fitFcn",fitFunction,0,3,6); + fitFcn->SetNpx(200); + gStyle->SetOptFit(); + gStyle->SetStatY(0.6); + + //with Minuit2 + c1->cd(1); + DoFit("Minuit2",gPad,npass,sigma,lambda); + + //with Fumili + c1->cd(2); + DoFit("Fumili",gPad,npass,sigma,lambda); + + //with cmaes + c1->cd(3); + DoFit("cmaes",gPad,npass,sigma,lambda); + + //with acmaes + c1->cd(4); + DoFit("acmaes",gPad,npass,sigma,lambda); + + //c1->SaveAs("FitBench.root"); + //delete c1; + delete fitFcn; +} diff --git a/tutorials/fit/cmaesFitBench2D.C b/tutorials/fit/cmaesFitBench2D.C new file mode 100644 index 0000000000000..52b7f2bd4c442 --- /dev/null +++ b/tutorials/fit/cmaesFitBench2D.C @@ -0,0 +1,97 @@ +// Author: L. Moneta 10/2005 +// Author: E. Benazera 05/2014 + +/********************************************************************** + * * + * Copyright (c) 2005 ROOT Foundation, CERN/PH-SFT * + * * + **********************************************************************/ + +#include "TH1.h" +#include "TF1.h" +#include "TH2D.h" +#include "TF2.h" +#include "TCanvas.h" +#include "TStopwatch.h" +#include "TSystem.h" +#include "TRandom3.h" +#include "TVirtualFitter.h" +#include "TPaveLabel.h" +#include "TStyle.h" + + +TF2 *fitFcn; +TH2D *histo; + +// Quadratic background function +Double_t gaus2D(Double_t *x, Double_t *par) { + double t1 = x[0] - par[1]; + double t2 = x[1] - par[2]; + return par[0]* exp( - 0.5 * ( t1*t1/( par[3]*par[3]) + t2*t2 /( par[4]*par[4] ) ) ) ; +} + +// Sum of background and peak function +Double_t fitFunction(Double_t *x, Double_t *par) { + return gaus2D(x,par); +} + +void fillHisto(int n =10000) { + + gRandom = new TRandom3(); + for (int i = 0; i < n; ++i) { + double x = gRandom->Gaus(2,3); + double y = gRandom->Gaus(-1,4); + histo->Fill(x,y,1.); + } +} + +void DoFit(const char* fitter, TVirtualPad *pad, Int_t npass) { + TStopwatch timer; + TVirtualFitter::SetDefaultFitter(fitter); + pad->SetGrid(); + fitFcn->SetParameters(100,0,0,2,7); + fitFcn->Update(); + + timer.Start(); + histo->Fit("fitFcn","0"); + timer.Stop(); + + histo->Draw(); + Double_t cputime = timer.CpuTime(); + printf("%s, npass=%d : RT=%7.3f s, Cpu=%7.3f s\n",fitter,npass,timer.RealTime(),cputime); + TPaveLabel *p = new TPaveLabel(0.5,0.7,0.85,0.8,Form("%s CPU= %g s",fitter,cputime),"brNDC"); + p->Draw(); + pad->Update(); +} + +void cmaesFitBench2D(int n = 100000) { + TH1::AddDirectory(kFALSE); + TCanvas *c1 = new TCanvas("c1","Fitting Demo",10,10,900,900); + c1->Divide(2,2); + // create a TF1 with the range from 0 to 3 and 6 parameters + fitFcn = new TF2("fitFcn",fitFunction,-10,10,-10,10,5); + //fitFcn->SetNpx(200); + gStyle->SetOptFit(); + gStyle->SetStatY(0.6); + + histo = new TH2D("h2","2D Gauss",100,-10,10,100,-10,10); + fillHisto(n); + + int npass=0; + + //with Minuit2 + c1->cd(3); + DoFit("Minuit2",gPad,npass); + + //with Fumili + c1->cd(4); + DoFit("Fumili",gPad,npass); + + //with cmaes + c1->cd(1); + DoFit("cmaes",gPad,npass); + + //with acmaes + c1->cd(2); + DoFit("acmaes",gPad,npass); +} diff --git a/tutorials/fit/cmaesFullBench.C b/tutorials/fit/cmaesFullBench.C new file mode 100644 index 0000000000000..fc04c6913cdf1 --- /dev/null +++ b/tutorials/fit/cmaesFullBench.C @@ -0,0 +1,1037 @@ +// Author: E. Benazera 6/2014 + +#include +#include +#include +#include + +#include "TH1.h" +#include "TF1.h" +#include "TH2D.h" +#include "TF2.h" +#include "TStopwatch.h" +#include "TRandom3.h" +#include "TVirtualFitter.h" +#include "TPaveLabel.h" +#include "TStyle.h" +#include "TSystem.h" +#include "TFile.h" +#include "TList.h" +#include "TMath.h" +#include "TROOT.h" + +class expstats +{ +public: + expstats(const std::string &name, + const int &dim, + const int &lambda=-1): + _name(name),_dim(dim),_lambda(lambda) {} + ~expstats() {} + + void add_exp(const bool &succ, + const double &fmin, + const std::vector &x, + const double &cputime, + const double &budget) + { + if (succ) + ++_succs; + else ++_fails; + _vsuccs.push_back(succ); + _fmin.push_back(fmin); + _x.push_back(x); + _cputime.push_back(cputime); + _cpu_avg = std::accumulate(_cputime.begin(),_cputime.end(),0.0) / static_cast(_cputime.size()); + _cpu_std = stddev(_cputime,_cpu_avg); + _budget.push_back(budget); + _budget_avg = std::accumulate(_budget.begin(),_budget.end(),0.0) / static_cast(_budget.size()); + _budget_std = stddev(_budget,_budget_avg); + } + + void merge(const expstats &stats) + { + for (size_t i=0;i _vsuccs; + std::vector _fmin; + std::vector> _x; + std::vector _cputime; + double _cpu_avg = 0.0; + double _cpu_std = 0.0; + std::vector _budget; + double _budget_avg = 0.0; + double _budget_std = 0.0; + + int _lambda = -1; + + // diff + int _found = 0; /* number of times the correct minima was found. */ + std::vector _fdiff; + std::vector _cputime_diff; + std::vector _cputime_ratio; + std::vector _budget_diff; + std::vector _budget_ratio; + double _cputime_diff_avg = 0.0; + double _cputime_ratio_avg = 1.0; + double _budget_diff_avg = 0.0; + double _budget_ratio_avg = 1.0; + int _isuccs = 0; /* number of times the best f-min among the two algorithms, was found. */ + int _ifails = 0; + int _iequals = 0; /* number of times the same fmin was found. */ +}; + +std::ostream& operator<<(std::ostream &out, const expstats &stats) +{ + return stats.print(out); +} + +typedef std::function ExpFunc; + +class experiment +{ +public: + experiment(const std::string &name) + :_name(name) + { + ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default("cmaes"); + opts.SetRealValue("sigma",0.1); + } + + void set_lambda(const int &lambda) + { + ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default("cmaes"); + opts.SetIntValue("lambda",lambda); + _lambda = lambda; + } + + void set_lscaling(const int &lscaling) + { + ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default("cmaes"); + opts.SetIntValue("lscaling",lscaling); + } + + void set_restarts(const int &nrestarts) + { + ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default("cmaes"); + opts.SetIntValue("restarts",nrestarts); + } + + virtual ~experiment() {} + + virtual void Setup() {} + + virtual void Cleanup() {} + + ExpFunc _ef; + std::string _name; + int _lambda = -1; +}; + +/*- gauss_fit -*/ +class gauss_fit_e : public experiment +{ +public: + gauss_fit_e() + :experiment("gauss_fit") + { + _ef = [this](const std::string &fitter) + { + std::string ename = "gauss_fit"; + + TVirtualFitter::SetDefaultFitter(fitter.c_str() ); + TH1D * h1 = new TH1D("h1f1","Chi2 Fit",100, -5, 5. ); + TH1D * h1bis = new TH1D("h1f2","Likelihood Fit",100, -5, 5. ); + for (int i = 0; i < _n; ++i) { + h1->Fill( _x.at(i) ); + h1bis->Fill( _x.at(i) ); + } + delete ((TF1 *)(gROOT->GetFunction("gaus"))); + TStopwatch timer; + timer.Start(); + TFitResultPtr r1 = h1->Fit("gaus","VS0"); + timer.Stop(); + Double_t cputime1 = timer.CpuTime(); + delete ((TF1 *)(gROOT->GetFunction("gaus"))); + TStopwatch timer2; + timer2.Start(); + TFitResultPtr r2 = h1bis->Fit("gaus","VLS0"); + timer2.Stop(); + Double_t cputime2 = timer2.CpuTime(); + + delete h1; + delete h1bis; + + expstats stats(ename,r1->NTotalParameters(),_lambda); + stats.add_exp(r1->Status()==0,r1->MinFcnValue(),r1->Parameters(),cputime1,r1->NCalls()); + stats.add_exp(r2->Status()==0,r2->MinFcnValue(),r2->Parameters(),cputime2,r2->NCalls()); + std::cout << "gaus_fit stats: " << stats << std::endl; + std::cout << "fmin1=" << r1->MinFcnValue() << std::endl;//" / fmin2=" << r2->MinFcnValue() << std::endl; + return stats; + }; + } + + ~gauss_fit_e() {} + + virtual void Setup() + { + std::cout << "setting up gauss_fit\n"; + for (int i=0;i<_n;i++) + _x.push_back(_grandom.Gaus(0,1)); + } + + int _n = 1000; + TRandom3 _grandom; + std::vector _x; +}; +gauss_fit_e ggauss_fit; + +/*- lorentz_fit -*/ +class lorentz_fit_e : public experiment +{ +public: + lorentz_fit_e() + :experiment("lorentz_fit") + { + _ef = [this](const std::string &fitter) + { + std::string title = "fit bench"; + std::string ename = "lorentz_fit"; + expstats stats(ename,6,_lambda); + TVirtualFitter::SetDefaultFitter(fitter.c_str()); + + TH1D *mhisto =new TH1D("fit",title.c_str(),200,0,3); + for (Int_t pass=0;pass<_npass;pass++) { + TStopwatch timer; + TF1 *fitFcn = new TF1("fitFcn",lorentz_fit_e::fitFunction,0,3,6); + fitFcn->SetParameters(1,1,1,6,.03,1); + TH1D *histo = new TH1D("fit2",title.c_str(),200,0,3); + //histo->Print("all"); + for (Int_t j=0;j<5000;j++) + { + histo->Fill(_xs.at(pass).at(j)); + mhisto->Fill(_xs.at(pass).at(j)); + } + /*timer.Start(); + TFitResultPtr r = histo->Fit(fitFcn,"VS0"); // from TH1.cxx: Q: quiet, 0: do not plot + timer.Stop(); + Double_t cputime = timer.CpuTime(); + stats.add_exp(r->Status()==0,r->MinFcnValue(),r->Parameters(),cputime,r->NCalls());*/ + delete histo; + delete fitFcn; + } + TStopwatch timer; + timer.Start(); + TF1 *fitFcn = new TF1("fitFcn",lorentz_fit_e::fitFunction,0,3,6); + fitFcn->SetParameters(1,1,1,6,.03,1); + TFitResultPtr r = mhisto->Fit(fitFcn,"VS0"); + timer.Stop(); + double cputime = timer.CpuTime(); + stats.add_exp(r->Status()==0,r->MinFcnValue(),r->Parameters(),cputime,r->NCalls()); + printf("%s, npass=%d : RT=%7.3f s, Cpu=%7.3f s\n",fitter.c_str(),_npass,timer.RealTime(),cputime); + //std::cout << "lorentz_fit stats: " << stats << std::endl; + delete mhisto; + delete fitFcn; + return stats; + }; + } + + virtual ~lorentz_fit_e() + { + _xs.clear(); + } + + // Quadratic background function + static Double_t background(Double_t *x, Double_t *par) { + return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; + } + + // Lorenzian Peak function + static Double_t lorentzianPeak(Double_t *x, Double_t *par) { + return (0.5*par[0]*par[1]/TMath::Pi()) / + TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + .25*par[1]*par[1]); + } + + // Sum of background and peak function + static Double_t fitFunction(Double_t *x, Double_t *par) { + return lorentz_fit_e::background(x,par) + lorentz_fit_e::lorentzianPeak(x,&par[3]); + } + + virtual void Setup() + { + TF1 *fitFcn = new TF1("fitFcn",lorentz_fit_e::fitFunction,0,3,6); + fitFcn->SetNpx(200); + fitFcn->SetParameters(1,1,1,6,.03,1); + fitFcn->Update(); + std::vector x; + for (int i=0;i<_npass;i++) + { + for (Int_t j=0;j<5000;j++) { + x.push_back(fitFcn->GetRandom()); + } + _xs.push_back(x); + } + delete fitFcn; + } + + virtual void Cleanup() + { + _xs.clear(); + } + + int _npass = 20; + //TF1 *_fitFcn; + std::vector > _xs; +}; +lorentz_fit_e glorentz_fit; + +/*- gauss2D_fit -*/ +class gauss2D_fit_e : public experiment +{ +public: + gauss2D_fit_e() + :experiment("gauss2D_fit") + { + _ef = [this](const std::string &fitter) + { + std::string ename = "gauss2D_fit"; + int npass = 0; + int n = 100000; + TStopwatch timer; + TVirtualFitter::SetDefaultFitter(fitter.c_str()); + TF2 *fitFcn = new TF2("fitFcn",gauss2D_fit_e::fitFunction2,-10,10,-10,10,5); + fitFcn->SetParameters(100,0,0,2,7); + fitFcn->Update(); + timer.Start(); + TFitResultPtr r = _histo->Fit("fitFcn","VS0"); + timer.Stop(); + Double_t cputime = timer.CpuTime(); + printf("%s : RT=%7.3f s, Cpu=%7.3f s\n",fitter.c_str(),timer.RealTime(),cputime); + delete fitFcn; + expstats stats(ename,r->NTotalParameters(),_lambda); + stats.add_exp(r->Status()==0,r->MinFcnValue(),r->Parameters(),cputime,r->NCalls()); + std::cout << "gauss2D_fit stats: " << stats << std::endl; + return stats; + }; + } + + ~gauss2D_fit_e() + { + if (_histo) + delete _histo; + } + + // Quadratic background function + static Double_t gaus2D(Double_t *x, Double_t *par) { + double t1 = x[0] - par[1]; + double t2 = x[1] - par[2]; + return par[0]* exp( - 0.5 * ( t1*t1/( par[3]*par[3]) + t2*t2 /( par[4]*par[4] ) ) ) ; + } + + // Sum of background and peak function + static Double_t fitFunction2(Double_t *x, Double_t *par) { + return gaus2D(x,par); + } + + static void fillHisto(int n, + TH2D *histo) + { + gRandom = new TRandom3(); + for (int i = 0; i < n; ++i) { + double x = gRandom->Gaus(2,3); + double y = gRandom->Gaus(-1,4); + histo->Fill(x,y,1.); + } + } + + virtual void Setup() + { + _histo = new TH2D("h2","2D Gauss",100,-10,10,100,-10,10); + fillHisto(_n,_histo); + } + + virtual void Cleanup() + { + if (_histo) + delete _histo; + _histo = nullptr; + } + + int _n = 100000; + TH2D *_histo = nullptr; +}; +gauss2D_fit_e ggauss2D_fit; + +/*- fit2a -*/ +class fit2a_e : public experiment +{ +public: + fit2a_e() + :experiment("fit2a") + { + _ef = [this](const std::string &fitter) + { + const Int_t npar = 15; + TVirtualFitter::SetDefaultFitter(fitter.c_str()); + TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar); + TH2F *h2 = new TH2F("h2","From f2",40,-10,10,40,-10,10); + for (int i=0;i<_nentries;i++) + h2->Fill(_x.at(i),_y.at(i)); + + Double_t f2params[npar] = {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; + Float_t ratio = 4*_nentries/100000; + f2params[ 0] *= ratio; + f2params[ 5] *= ratio; + f2params[10] *= ratio; + f2->SetParameters(f2params); + + //Fit h2 with original function f2 + TStopwatch timer; + timer.Start(); + TFitResultPtr r = h2->Fit("f2","SN0"); + timer.Stop(); + Double_t cputime = timer.CpuTime(); + expstats stats("fit2a",r->NTotalParameters(),_lambda); + stats.add_exp(r->Status()==0,r->MinFcnValue(),r->Parameters(),cputime,r->NCalls()); + delete h2; + delete f2; + return stats; + }; + } + + ~fit2a_e() + { + Cleanup(); + } + + static Double_t g2(Double_t *x, Double_t *par) { + Double_t r1 = Double_t((x[0]-par[1])/par[2]); + Double_t r2 = Double_t((x[1]-par[3])/par[4]); + return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2)); + } + + static Double_t fun2(Double_t *x, Double_t *par) { + Double_t *p1 = &par[0]; + Double_t *p2 = &par[5]; + Double_t *p3 = &par[10]; + Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3); + return result; + } + + virtual void Setup() + { + const Int_t npar = 15; + Double_t f2params[npar] = {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; + TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar); + f2->SetParameters(f2params); + _x = std::vector(_nentries); + _y = std::vector(_nentries); + for (int i=0;i<_nentries;i++) + f2->GetRandom2(_x.at(i),_y.at(i)); + delete f2; + } + + virtual void Cleanup() + { + } + + int _nentries = 100000; + std::vector _x,_y; +}; +fit2a_e gfit2a; + +/*- fit2dhist -*/ +class fit2dhist_e : public experiment +{ +public: + fit2dhist_e() + :experiment("fit2dhist") + { + _ef = [this](const std::string &fitter) + { + int nbx1 = 50; + int nby1 = 50; + int nbx2 = 50; + int nby2 = 50; + double xlow1 = 0.; + double ylow1 = 0.; + double xup1 = 10.; + double yup1 = 10.; + double xlow2 = 5.; + double ylow2 = 5.; + double xup2 = 20.; + double yup2 = 20.; + TH2D *h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1); + TH2D *h2 = new TH2D("h2","tails",nbx2,xlow2,xup2,nby2,ylow2,yup2); + double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. }; + TF2 *func1 = new TF2("func",fit2dhist_e::my2Dfunc,xlow2,xup2,ylow2,yup2, 10); + func1->SetParameters(iniParams); + TF2 *func2 = new TF2("func",fit2dhist_e::my2Dfunc,xlow2,xup2,ylow2,yup2, 10); + func2->SetParameters(iniParams); + + // fill up histograms. + int n1 = 1000000; + int n2 = 1000000; + for (int i=0;iFill(_xr1.at(i),_yr1.at(i)); + h2->Fill(_xr2.at(i),_yr2.at(i)); // since n1 == n2. + } + + // scale histograms to same heights (for fitting) + double dx1 = (xup1-xlow1)/double(nbx1); + double dy1 = (yup1-ylow1)/double(nby1); + double dx2 = (xup2-xlow2)/double(nbx2); + double dy2 = (yup2-ylow2)/double(nby2); + // scale histo 2 to scale of 1 + h2->Sumw2(); + h2->Scale( ( double(n1) * dx1 * dy1 ) / ( double(n2) * dx2 * dy2 ) ); + + TStopwatch timer; + timer.Start(); + TVirtualFitter::SetDefaultFitter(fitter.c_str()); + TFitResultPtr r1 = h1->Fit(func1,"VS0"); + timer.Stop(); + Double_t cputime1 = timer.CpuTime(); + TStopwatch timer2; + timer2.Start(); + TFitResultPtr r2 = h2->Fit(func2,"VS0"); + timer2.Stop(); + Double_t cputime2 = timer2.CpuTime(); + expstats stats("fit2dhist",r1->NTotalParameters(),_lambda); + stats.add_exp(r1->Status()==0,r1->MinFcnValue(),r1->Parameters(),cputime1,r1->NCalls()); + stats.add_exp(r2->Status()==0,r2->MinFcnValue(),r2->Parameters(),cputime2,r2->NCalls()); + delete h1; + delete h2; + delete func1; + delete func2; + return stats; + }; + } + + ~fit2dhist_e() + { + Cleanup(); + } + + static double gauss2D(double *x, double *par) { + double z1 = double((x[0]-par[1])/par[2]); + double z2 = double((x[1]-par[3])/par[4]); + return par[0]*exp(-0.5*(z1*z1+z2*z2)); + } + + static double my2Dfunc(double *x, double *par) { + return gauss2D(x,&par[0]) + gauss2D(x,&par[5]); + } + + static void FillHisto2(std::vector &xr, std::vector &yr, int n, double * p) + { + const double mx1 = p[1]; + const double my1 = p[3]; + const double sx1 = p[2]; + const double sy1 = p[4]; + const double mx2 = p[6]; + const double my2 = p[8]; + const double sx2 = p[7]; + const double sy2 = p[9]; + //const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2); + const double w1 = 0.5; + + double x, y; + for (int i = 0; i < n; ++i) { + // generate randoms with larger gaussians + fit2dhist_e::_rndm.Rannor(x,y); + + double r = fit2dhist_e::_rndm.Rndm(1); + if (r < w1) { + x = x*sx1 + mx1; + y = y*sy1 + my1; + } + else { + x = x*sx2 + mx2; + y = y*sy2 + my2; + } + xr.push_back(x); + yr.push_back(y); + } + } + + void Setup() + { + double xlow1 = 0.; + double ylow1 = 0.; + double xup1 = 10.; + double yup1 = 10.; + double xlow2 = 5.; + double ylow2 = 5.; + double xup2 = 20.; + double yup2 = 20.; + double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. }; + int n1 = 1000000; + int n2 = 1000000; + FillHisto2(_xr1,_yr1,n1,iniParams); + FillHisto2(_xr2,_yr2,n2,iniParams); + } + + void Cleanup() + { + } + + static TRandom3 _rndm; + std::vector _xr1,_yr1,_xr2,_yr2; +}; +TRandom3 fit2dhist_e::_rndm = TRandom3(); +fit2dhist_e gfit2dhist; + +/*- combined_fit -*/ +class combined_fit_e : public experiment +{ +public: + struct GlobalChi2 { + GlobalChi2( ROOT::Math::IMultiGenFunction & f1, + ROOT::Math::IMultiGenFunction & f2) : + fChi2_1(&f1), fChi2_2(&f2) {} + + // parameter vector is first background (in common 1 and 2) + // and then is signal (only in 2) + double operator() (const double *par) const { + int iparB[2] = {0,2}; + int iparSB[5] = {1,2,3,4,5}; + double p1[2]; + for (int i = 0; i < 2; ++i) p1[i] = par[iparB[i] ]; + + double p2[5]; + for (int i = 0; i < 5; ++i) p2[i] = par[iparSB[i] ]; + + return (*fChi2_1)(p1) + (*fChi2_2)(p2); + } + + const ROOT::Math::IMultiGenFunction * fChi2_1; + const ROOT::Math::IMultiGenFunction * fChi2_2; + }; + + combined_fit_e() + :experiment("combined_fit") + { + _ef = [this](const std::string &fitter) + { + // perform now global fit + TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100); + + ROOT::Math::WrappedMultiTF1 wfB(*_fB,1); + ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1); + + ROOT::Fit::DataOptions opt; + ROOT::Fit::DataRange rangeB; + // set the data range + rangeB.SetRange(10,90); + ROOT::Fit::BinData dataB(opt,rangeB); + ROOT::Fit::FillData(dataB, _hB); + + ROOT::Fit::DataRange rangeSB; + rangeSB.SetRange(10,50); + ROOT::Fit::BinData dataSB(opt,rangeSB); + ROOT::Fit::FillData(dataSB, _hSB); + + ROOT::Fit::Chi2Function chi2_B(dataB, wfB); + ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB); + + GlobalChi2 globalChi2(chi2_B, chi2_SB); + + ROOT::Fit::Fitter rfitter; + + const int Npar = 6; + double par0[Npar] = { 5,5,-0.1,100, 30,10}; + + // create before the parameter settings in order to fix or set range on them + rfitter.Config().SetParamsSettings(6,par0); + // fix 5-th parameter + //rfitter.Config().ParSettings(4).Fix(); // weird random crash, not yet understood: fitter tries to set variable 4 instead of 5, sometimes... + // set limits on the third and 4-th parameter + rfitter.Config().ParSettings(2).SetLimits(-10,-1.E-4); + rfitter.Config().ParSettings(3).SetLimits(0,10000); + rfitter.Config().ParSettings(3).SetStepSize(5); + + rfitter.Config().MinimizerOptions().SetPrintLevel(0); + if (fitter == "Minuit2") + rfitter.Config().SetMinimizer("Minuit2","Migrad"); + else rfitter.Config().SetMinimizer("cmaes",fitter.c_str()); + + // fit FCN function directly + // (specify optionally data size and flag to indicate that is a chi2 fit) + TStopwatch timer; + timer.Start(); + std::cout << "starting\n"; + rfitter.FitFCN(6,globalChi2,0,dataB.Size()+dataSB.Size(),true); + timer.Stop(); + Double_t cputime = timer.CpuTime(); + ROOT::Fit::FitResult r = rfitter.Result(); + //result.Print(std::cout); + + delete fSB; + expstats stats("combined",r.NTotalParameters(),_lambda); + stats.add_exp(r.Status()==0,r.MinFcnValue(),r.Parameters(),cputime,r.NCalls()); + std::cout << "combined stats: " << stats << std::endl; + std::cout << "fmin=" << r.MinFcnValue() << std::endl; + return stats; + }; + } + + ~combined_fit_e() + { + Cleanup(); + } + + virtual void Setup() + { + _hB = new TH1D("hB","histo B",100,0,100); + _hSB = new TH1D("hSB","histo S+B",100, 0,100); + _fB = new TF1("fB","expo",0,100); + _fB->SetParameters(1,-0.05); + _hB->FillRandom("fB"); + _fS = new TF1("fS","gaus",0,100); + _fS->SetParameters(1,30,5); + _hSB->FillRandom("fB",2000); + _hSB->FillRandom("fS",1000); + } + + virtual void Cleanup() + { + if (_hB) + delete _hB; + if (_hSB) + delete _hSB; + if (_fB) + delete _fB; + if (_fS) + delete _fS; + } + + TH1D *_hB = nullptr; + TH1D *_hSB = nullptr; + TF1 *_fB = nullptr; + TF1 *_fS = nullptr; +}; +combined_fit_e gcombined_fit; + +/*- ex3d -*/ +class ex3d_e : public experiment +{ +public: + ex3d_e() + :experiment("ex3d") + { + _ef = [this](const std::string &fitter) + { + double ev = 0.1; + + // create a 3d binned data structure + ROOT::Fit::BinData data(_n,3); + double xx[3]; + for(int i = 0; i < _n; ++i) { + xx[0] = _x[i]; + xx[1] = _y[i]; + xx[2] = _z[i]; + // add the 3d-data coordinate, the predictor value (v[i]) and its errors + data.Add(xx, _v[i], ev); + } + + TF3 * f3 = new TF3("f3","[0] * sin(x) + [1] * cos(y) + [2] * z",0,10,0,10,0,10); + f3->SetParameters(2,2,2); + ROOT::Fit::Fitter rfitter; + /*if (fitter.find("cmaes")!=std::string::npos) + rfitter.Config().SetMinimizer("cmaes","acmaes");*/ + if (fitter != "Minuit2") + rfitter.Config().SetMinimizer("cmaes",fitter.c_str()); + // wrapped the TF1 in a IParamMultiFunction interface for the Fitter class + ROOT::Math::WrappedMultiTF1 wf(*f3,3); + rfitter.SetFunction(wf); + TStopwatch timer; + timer.Start(); + bool ret = rfitter.Fit(data); + timer.Stop(); + Double_t cputime = timer.CpuTime(); + const ROOT::Fit::FitResult & res = rfitter.Result(); + if (ret) { + // print result (should be around 1) + res.Print(std::cout); + // copy all fit result info (values, chi2, etc..) in TF3 + f3->SetFitResult(res); + // test fit p-value (chi2 probability) + double prob = res.Prob(); + if (prob < 1.E-2) + Error("exampleFit3D","Bad data fit - fit p-value is %f",prob); + else + std::cout << "Good fit : p-value = " << prob << std::endl; + } + else Error("exampleFit3D","3D fit failed"); + expstats stats("example3D",res.NTotalParameters(),_lambda); + stats.add_exp(res.Status()==0,res.MinFcnValue(),res.Parameters(),cputime,res.NCalls()); + delete f3; + return stats; + }; + } + + ~ex3d_e() + { + } + + virtual void Setup() + { + double ev = 0.1; + TRandom2 r; + for (int i = 0; i < _n; ++i) { + _x.push_back(r.Uniform(0,10)); + _y.push_back(r.Uniform(0,10)); + _z.push_back(r.Uniform(0,10)); + _v.push_back(sin(_x[i] ) + cos(_y[i]) + _z[i] + r.Gaus(0,ev)); + } + } + + virtual void Cleanup() + { + _x.clear(); + _y.clear(); + _z.clear(); + _v.clear(); + } + + int _n = 1000; + std::vector _x; + std::vector _y; + std::vector _z; + std::vector _v; +}; +ex3d_e gex3d; + +/*- fit2 -*/ +class fit2_e : public experiment +{ +public: + fit2_e() + :experiment("fit2") + { + _ef = [this](const std::string &fitter) + { + TVirtualFitter::SetDefaultFitter(fitter.c_str()); + TF2 *f2 = new TF2("f2",fit2_e::fun22,-10,10,-10,10, _npar); + + Double_t f2params[15] = + {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; + Float_t ratio = 4*_nentries/100000; + f2params[ 0] *= ratio; + f2params[ 5] *= ratio; + f2params[10] *= ratio; + f2->SetParameters(f2params); + + TH2F *h2 = new TH2F("h2","from f2",40,-10,10,40,-10,10); + for (int i=0;i<_nentries;i++) + h2->Fill(_x.at(i),_y.at(i)); + + //Fit h2 with original function f2 + TStopwatch timer; + timer.Start(); + TFitResultPtr r = h2->Fit("f2","VS0"); + timer.Stop(); + Double_t cputime = timer.CpuTime(); + expstats stats("fit2",r->NTotalParameters(),_lambda); + stats.add_exp(r->Status()==0,r->MinFcnValue(),r->Parameters(),cputime,r->NCalls()); + delete h2; + delete f2; + return stats; + }; + } + + ~fit2_e() + { + } + + virtual void Setup() + { + TF2 *f2 = new TF2("f2",fit2_e::fun22,-10,10,-10,10, _npar); + Double_t f2params[15] = + {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; + f2->SetParameters(f2params); + _x = std::vector(_nentries); + _y = std::vector(_nentries); + for (int i=0;i<_nentries;i++) + f2->GetRandom2(_x.at(i),_y.at(i)); + delete f2; + } + + static Double_t g22(Double_t *x, Double_t *par) { + Double_t r1 = Double_t((x[0]-par[1])/par[2]); + Double_t r2 = Double_t((x[1]-par[3])/par[4]); + return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2)); + } + + static Double_t fun22(Double_t *x, Double_t *par) { + Double_t *p1 = &par[0]; + Double_t *p2 = &par[5]; + Double_t *p3 = &par[10]; + Double_t result = fit2_e::g22(x,p1) + fit2_e::g22(x,p2) + fit2_e::g22(x,p3); + return result; + } + + Int_t _nentries = 100000; + Int_t _npar = 15; + std::vector _x, _y; +}; +fit2_e gfit2; + +void cmaesFullBench(const int &n=100, + const int &lscaling=1) +{ + std::cout << "Proceeding with " << n << " runs on every problems\n"; + if (lscaling > 0) + std::cout << "Linear scaling of parameters in ON\n"; + //std::vector lambdas = {-1, 5, 10, 20, 40, 80, 160, 320, 640, 1280}; + std::vector lambdas = {-1, 50, 200, -2, -3}; + std::vector acmaes_stats; + std::vector minuit2_stats; + std::map mexperiments; + mexperiments.insert(std::pair(ggauss_fit._name,&ggauss_fit)); + mexperiments.insert(std::pair(glorentz_fit._name,&glorentz_fit)); + mexperiments.insert(std::pair(gfit2._name,&gfit2)); + mexperiments.insert(std::pair(ggauss2D_fit._name,&ggauss2D_fit)); + mexperiments.insert(std::pair(gfit2a._name,&gfit2a)); + mexperiments.insert(std::pair(gfit2dhist._name,&gfit2dhist)); + mexperiments.insert(std::pair(gcombined_fit._name,&gcombined_fit)); + mexperiments.insert(std::pair(gex3d._name,&gex3d)); + int nexp = mexperiments.size(); + int cn = 0; + std::map::iterator mit = mexperiments.begin(); + while(mit!=mexperiments.end()) + { + std::cout << "Running " << (*mit).first << std::endl; + for (int i=0;iSetup(); + for (int j=0;j<(int)lambdas.size();j++) + { + std::string fitter_name = "acmaes"; + if (lambdas.at(j) >= -1) + (*mit).second->set_lambda(lambdas.at(j)); + else if (lambdas.at(j) == -2) + { + (*mit).second->set_lambda(-1); + (*mit).second->set_restarts(4); + fitter_name = "aipop"; + } + else if (lambdas.at(j) == -3) + { + (*mit).second->set_lambda(-1); + (*mit).second->set_restarts(10); + fitter_name = "abipop"; + } + (*mit).second->set_lscaling(lscaling); + + if (i == 0) + { + acmaes_stats.push_back((*mit).second->_ef(fitter_name)); + } + else + { + acmaes_stats.at(cn*(lambdas.size())+j).merge((*mit).second->_ef(fitter_name)); + } + } + (*mit).second->set_lambda(-1); // N/A to Minuit2 + if (i == 0) + minuit2_stats.push_back((*mit).second->_ef("Minuit2")); + else minuit2_stats.back().merge((*mit).second->_ef("Minuit2")); + (*mit).second->Cleanup(); + } + ++cn; + ++mit; + } + + std::cout << "nexp=" << nexp << " / stats size=" << acmaes_stats.size() << std::endl; + + for (size_t i=0;i +#include + +bool libloaded = false; + +void testGausFit( std::string type = "cmaes", int n = 1000) { + + gRandom = new TRandom3(); + + TVirtualFitter::SetDefaultFitter(type.c_str() ); + + std::string name; + name = "h1_" + type; + TH1D * h1 = new TH1D(name.c_str(),"Chi2 Fit",100, -5, 5. ); + name = "h1bis_" + type; + TH1D * h1bis = new TH1D(name.c_str(),"Likelihood Fit",100, -5, 5. ); + /*name = "h2_" + type; + TH1D * h2 = new TH1D(name.c_str(),"Chi2 Fit with Minos Error",100, -5, 5. ); + name = "h3_" + type; + TH1D * h3 = new TH1D(name.c_str(),"Chi2 Fit with Integral and Minos",100, -5, 5. ); + name = "h4_" + type; + TH1D * h4 = new TH1D(name.c_str(),"Likelihood Fit with Minos Error",100, -5, 5. );*/ + + gStyle->SetOptStat(1111111); + gStyle->SetOptFit(1111111); + + for (int i = 0; i < n; ++i) { + double x = gRandom->Gaus(0,1); + h1->Fill( x ); + h1bis->Fill( x ); + /* h2->Fill( x ); + h3->Fill( x ); + h4->Fill( x ); */ + } + + std::string cname = type + "Canvas" ; + std::string ctitle = type + " Gaussian Fit" ; + TCanvas *c1 = new TCanvas(cname.c_str(),cname.c_str(),10,10,900,900); + c1->Divide(2,1); + + c1->cd(1); + cout << "\nDo Fit 1\n"; + h1->Fit("gaus","Q"); + h1->Draw(); + c1->cd(2); + cout << "\nDo Fit 1bis\n"; + h1bis->Fit("gaus","VLE"); + h1bis->Draw(); + /*c1->cd(2); + cout << "\nDo Fit 2\n"; + h2->Fit("gaus","VE"); + h2->Draw(); + c1->cd(3); + cout << "\nDo Fit 3\n"; + h3->Fit("gaus","IE"); + h3->Draw(); + c1->cd(4); + cout << "\nDo Fit 4\n"; + h4->Fit("gaus","VLE"); + h4->Draw();*/ + +} + +void cmaesGausFit() { + if (!libloaded) + { + libloaded = true; + } + int n = 1000; + testGausFit("cmaes",n); + //testGausFit("acmaes",n); +} + + + diff --git a/tutorials/fit/cmaes_fit2a.C b/tutorials/fit/cmaes_fit2a.C new file mode 100644 index 0000000000000..3911b4316e057 --- /dev/null +++ b/tutorials/fit/cmaes_fit2a.C @@ -0,0 +1,82 @@ +#include "TF2.h" +#include "TH2.h" +#include "TCutG.h" +#include "TMath.h" +#include "TCanvas.h" +#include "TStyle.h" +#include "TVirtualFitter.h" + + +//+ Fitting a 2-D histogram (a variant) +// This tutorial illustrates : +// - how to create a 2-d function +// - fill a 2-d histogram randomly from this function +// - fit the histogram +// - display the fitted function on top of the histogram (lego-plot) +// using a surface plot in a sub-range of the histogram. +// +// This example can be executed via the interpreter or/and the compiler +// root > .x fit2a.C +// root > .x fit2a.C++ +//Author: Rene Brun +//Author: Emmanuel Benazera + +bool libloaded = false; + +Double_t g2(Double_t *x, Double_t *par) { + Double_t r1 = Double_t((x[0]-par[1])/par[2]); + Double_t r2 = Double_t((x[1]-par[3])/par[4]); + return par[0]*TMath::Exp(-0.5*(r1*r1+r2*r2)); +} +Double_t fun2(Double_t *x, Double_t *par) { + Double_t *p1 = &par[0]; + Double_t *p2 = &par[5]; + Double_t *p3 = &par[10]; + Double_t result = g2(x,p1) + g2(x,p2) + g2(x,p3); + return result; +} + +TCanvas *fit2a(const char *fitter="cmaes") { + if (!libloaded) + { + libloaded = true; + } + + TVirtualFitter::SetDefaultFitter(fitter); + + ROOT::Math::IOptions &opts = ROOT::Math::MinimizerOptions::Default(fitter); + //opts.SetIntValue("lambda",100); + opts.SetNamedValue("fplot","fit2a.dat"); + + TCanvas *c = new TCanvas(); + gStyle->SetOptStat(kTRUE); + gStyle->SetPalette(1); + const Int_t npar = 15; + Double_t f2params[npar] = {100,-3,3,-3,3,160,0,0.8,0,0.9,40,4,0.7,4,0.7}; + TF2 *f2 = new TF2("f2",fun2,-10,10,-10,10, npar); + f2->SetParameters(f2params); + + //Create an histogram and fill it randomly with f2 + TH2F *h2 = new TH2F("h2","From f2",40,-10,10,40,-10,10); + Int_t nentries = 100000; + h2->FillRandom("f2",nentries); + //Fit h2 with original function f2 + Float_t ratio = 4*nentries/100000; + f2params[ 0] *= ratio; + f2params[ 5] *= ratio; + f2params[10] *= ratio; + f2->SetParameters(f2params); + h2->Fit("f2","N"); + TCutG *cutg = new TCutG("cutg",5); + cutg->SetPoint(0,-7,-7); + cutg->SetPoint(1, 2,-7); + cutg->SetPoint(2, 2, 2); + cutg->SetPoint(3,-7, 2); + cutg->SetPoint(4,-7,-7); + h2->Draw("lego2 0"); + h2->SetFillColor(38); + f2->SetNpx(80); + f2->SetNpy(80); + f2->Draw("surf1 same bb [cutg]"); + return c; +} diff --git a/tutorials/fit/cmaes_fit2dHist.C b/tutorials/fit/cmaes_fit2dHist.C new file mode 100644 index 0000000000000..607aaff2f0e02 --- /dev/null +++ b/tutorials/fit/cmaes_fit2dHist.C @@ -0,0 +1,244 @@ +// +//+ Example to fit two histograms at the same time via TVirtualFitter +// +// To execute this tutorial, you can do: +// +// root > .x fit2dHist.C (executing via CINT, slow) +// or +// root > .x fit2dHist.C+ (executing via ACLIC , fast, with Minuit) +// root > .x fit2dHist.C+(2) (executing via ACLIC , fast, with Minuit2) +// or using the option to fit independently the 2 histos +// root > .x fit2dHist.C+(10) (via ACLIC, fast, independent fits with Minuit) +// root > .x fit2dHist.C+(12) (via ACLIC, fast, independent fits with Minuit2) +// +// Note that you can also execute this script in batch with eg, +// root -b -q "fit2dHist.C+(12)" +// +// or execute interactively from the shell +// root fit2dHist.C+ +// root "fit2dHist.C+(12)" +// +// Authors: Lorenzo Moneta, Rene Brun 18/01/2006 + +#include "TH2D.h" +#include "TF2.h" +#include "TCanvas.h" +#include "TStyle.h" +#include "TRandom3.h" +#include "TVirtualFitter.h" +#include "TList.h" +#include "TSystem.h" + +#include + +double gauss2D(double *x, double *par) { + double z1 = double((x[0]-par[1])/par[2]); + double z2 = double((x[1]-par[3])/par[4]); + return par[0]*exp(-0.5*(z1*z1+z2*z2)); +} +double my2Dfunc(double *x, double *par) { + return gauss2D(x,&par[0]) + gauss2D(x,&par[5]); +} + + +// data need to be globals to be visible by fcn +TRandom3 rndm; +TH2D *h1, *h2; +Int_t npfits; + +void myFcn(Int_t & /*nPar*/, Double_t * /*grad*/ , Double_t &fval, Double_t *p, Int_t /*iflag */ ) +{ + TAxis *xaxis1 = h1->GetXaxis(); + TAxis *yaxis1 = h1->GetYaxis(); + TAxis *xaxis2 = h2->GetXaxis(); + TAxis *yaxis2 = h2->GetYaxis(); + + int nbinX1 = h1->GetNbinsX(); + int nbinY1 = h1->GetNbinsY(); + int nbinX2 = h2->GetNbinsX(); + int nbinY2 = h2->GetNbinsY(); + + double chi2 = 0; + double x[2]; + double tmp; + npfits = 0; + for (int ix = 1; ix <= nbinX1; ++ix) { + x[0] = xaxis1->GetBinCenter(ix); + for (int iy = 1; iy <= nbinY1; ++iy) { + if ( h1->GetBinError(ix,iy) > 0 ) { + x[1] = yaxis1->GetBinCenter(iy); + tmp = (h1->GetBinContent(ix,iy) - my2Dfunc(x,p))/h1->GetBinError(ix,iy); + chi2 += tmp*tmp; + npfits++; + } + } + } + for (int ix = 1; ix <= nbinX2; ++ix) { + x[0] = xaxis2->GetBinCenter(ix); + for (int iy = 1; iy <= nbinY2; ++iy) { + if ( h2->GetBinError(ix,iy) > 0 ) { + x[1] = yaxis2->GetBinCenter(iy); + tmp = (h2->GetBinContent(ix,iy) - my2Dfunc(x,p))/h2->GetBinError(ix,iy); + chi2 += tmp*tmp; + npfits++; + } + } + } + fval = chi2; +} + +void FillHisto(TH2D * h, int n, double * p) { + + + const double mx1 = p[1]; + const double my1 = p[3]; + const double sx1 = p[2]; + const double sy1 = p[4]; + const double mx2 = p[6]; + const double my2 = p[8]; + const double sx2 = p[7]; + const double sy2 = p[9]; + //const double w1 = p[0]*sx1*sy1/(p[5]*sx2*sy2); + const double w1 = 0.5; + + double x, y; + for (int i = 0; i < n; ++i) { + // generate randoms with larger gaussians + rndm.Rannor(x,y); + + double r = rndm.Rndm(1); + if (r < w1) { + x = x*sx1 + mx1; + y = y*sy1 + my1; + } + else { + x = x*sx2 + mx2; + y = y*sy2 + my2; + } + h->Fill(x,y); + + } +} + + + + +int fit2dHist(int option=1) { + + // create two histograms + + int nbx1 = 50; + int nby1 = 50; + int nbx2 = 50; + int nby2 = 50; + double xlow1 = 0.; + double ylow1 = 0.; + double xup1 = 10.; + double yup1 = 10.; + double xlow2 = 5.; + double ylow2 = 5.; + double xup2 = 20.; + double yup2 = 20.; + + h1 = new TH2D("h1","core",nbx1,xlow1,xup1,nby1,ylow1,yup1); + h2 = new TH2D("h2","tails",nbx2,xlow2,xup2,nby2,ylow2,yup2); + + double iniParams[10] = { 100, 6., 2., 7., 3, 100, 12., 3., 11., 2. }; + // create fit function + TF2 * func = new TF2("func",my2Dfunc,xlow2,xup2,ylow2,yup2, 10); + func->SetParameters(iniParams); + + // fill Histos + int n1 = 1000000; + int n2 = 1000000; + FillHisto(h1,n1,iniParams); + FillHisto(h2,n2,iniParams); + + // scale histograms to same heights (for fitting) + double dx1 = (xup1-xlow1)/double(nbx1); + double dy1 = (yup1-ylow1)/double(nby1); + double dx2 = (xup2-xlow2)/double(nbx2); + double dy2 = (yup2-ylow2)/double(nby2); + // scale histo 2 to scale of 1 + h2->Sumw2(); + h2->Scale( ( double(n1) * dx1 * dy1 ) / ( double(n2) * dx2 * dy2 ) ); + + bool global = false; + if (option > 10) global = true; + if (global) { + // fill data structure for fit (coordinates + values + errors) + std::cout << "Do global fit" << std::endl; + // fit now all the function together + + //The default minimizer is Minuit, you can also try Minuit2 + if (option%10 == 2) + TVirtualFitter::SetDefaultFitter("Minuit2"); + else TVirtualFitter::SetDefaultFitter("Minuit"); + TVirtualFitter * minuit = TVirtualFitter::Fitter(0,10); + for (int i = 0; i < 10; ++i) { + minuit->SetParameter(i, func->GetParName(i), func->GetParameter(i), 0.01, 0,0); + } + minuit->SetFCN(myFcn); + + double arglist[100]; + arglist[0] = 0; + // set print level + minuit->ExecuteCommand("SET PRINT",arglist,2); + + // minimize + arglist[0] = 5000; // number of function calls + arglist[1] = 0.01; // tolerance + minuit->ExecuteCommand("MIGRAD",arglist,2); + + //get result + double minParams[10]; + double parErrors[10]; + for (int i = 0; i < 10; ++i) { + minParams[i] = minuit->GetParameter(i); + parErrors[i] = minuit->GetParError(i); + } + double chi2, edm, errdef; + int nvpar, nparx; + minuit->GetStats(chi2,edm,errdef,nvpar,nparx); + + func->SetParameters(minParams); + func->SetParErrors(parErrors); + func->SetChisquare(chi2); + int ndf = npfits-nvpar; + func->SetNDF(ndf); + + // add to list of functions + h1->GetListOfFunctions()->Add(func); + h2->GetListOfFunctions()->Add(func); + } + else { + // fit independently + TVirtualFitter::SetDefaultFitter("cmaes"); + h1->Fit(func); + h2->Fit(func); + } + + // Create a new canvas. + TCanvas * c1 = new TCanvas("c1","Two HIstogram Fit example",100,10,900,800); + c1->Divide(2,2); + gStyle->SetOptFit(); + gStyle->SetStatY(0.6); + + c1->cd(1); + h1->Draw(); + func->SetRange(xlow1,ylow1,xup1,yup1); + func->DrawCopy("cont1 same"); + c1->cd(2); + h1->Draw("lego"); + func->DrawCopy("surf1 same"); + c1->cd(3); + func->SetRange(xlow2,ylow2,xup2,yup2); + h2->Draw(); + func->DrawCopy("cont1 same"); + c1->cd(4); + h2->Draw("lego"); + gPad->SetLogz(); + func->Draw("surf1 same"); + + return 0; +} diff --git a/tutorials/fit/minuit2GausFit.C b/tutorials/fit/minuit2GausFit.C index e760c6843c210..8f40ba7afdbea 100644 --- a/tutorials/fit/minuit2GausFit.C +++ b/tutorials/fit/minuit2GausFit.C @@ -77,7 +77,6 @@ void minuit2GausFit() { int n = 1000; testGausFit("Minuit2",n); testGausFit("Fumili2",n); - }