diff --git a/cmake/modules/SearchInstalledSoftware.cmake b/cmake/modules/SearchInstalledSoftware.cmake index 5a7e7b47f50cc..65d5ea866b6f3 100644 --- a/cmake/modules/SearchInstalledSoftware.cmake +++ b/cmake/modules/SearchInstalledSoftware.cmake @@ -678,8 +678,6 @@ set(Python3_FIND_FRAMEWORK LAST) list(APPEND python_components Interpreter) if(pyroot OR tmva-pymva) list(APPEND python_components Development) -endif() -if(tmva-pymva) list(APPEND python_components NumPy) endif() find_package(Python3 3.8 COMPONENTS ${python_components}) diff --git a/etc/plugins/ROOT@@Math@@Minimizer/P100_Scipy.C b/etc/plugins/ROOT@@Math@@Minimizer/P100_Scipy.C new file mode 100644 index 0000000000000..03b547d306697 --- /dev/null +++ b/etc/plugins/ROOT@@Math@@Minimizer/P100_Scipy.C @@ -0,0 +1,5 @@ +void P100_Scipy() +{ + gPluginMgr->AddHandler("ROOT::Math::Minimizer", "Scipy", "ROOT::Math::Experimental::ScipyMinimizer", "Scipy", + "ScipyMinimizer(const char*)"); +} diff --git a/math/CMakeLists.txt b/math/CMakeLists.txt index 543adbf6b1e8d..becbf16b1e7cd 100644 --- a/math/CMakeLists.txt +++ b/math/CMakeLists.txt @@ -34,4 +34,6 @@ if(r) add_subdirectory(rtools) endif() +add_subdirectory(scipy) + add_subdirectory(vecops) diff --git a/math/scipy/CMakeLists.txt b/math/scipy/CMakeLists.txt new file mode 100644 index 0000000000000..76a57ccffb7e8 --- /dev/null +++ b/math/scipy/CMakeLists.txt @@ -0,0 +1,22 @@ +############################################################################ +# CMakeLists.txt file for building ROOT math/scipy package +############################################################################ +# author: Omar.Zapata@cern.ch 2022 + +include_directories(SYSTEM ${PYTHON_INCLUDE_DIRS_Development_Main} ${NUMPY_INCLUDE_DIRS}) +ROOT_STANDARD_LIBRARY_PACKAGE(Scipy + HEADERS + Math/ScipyMinimizer.h + LINKDEF + Math/LinkDef.h + LIBRARIES + Python3::NumPy + Python3::Python + SOURCES + src/ScipyMinimizer.cxx + DEPENDENCIES Core MathCore RIO + ) + +target_compile_definitions(Scipy PUBLIC USE_ROOT_ERROR ${PYTHON_DEFINITIONS}) +target_include_directories(Scipy PUBLIC ${PYTHON_INCLUDE_DIRS}) +ROOT_ADD_TEST_SUBDIRECTORY(test) diff --git a/math/scipy/inc/Math/LinkDef.h b/math/scipy/inc/Math/LinkDef.h new file mode 100644 index 0000000000000..472370c6d513b --- /dev/null +++ b/math/scipy/inc/Math/LinkDef.h @@ -0,0 +1,11 @@ +// @(#)root/math/ipopt:$Id$ +// Authors: Omar.Zapata@cern.ch 12/2022 + +#ifdef __CINT__ + +#pragma link C++ nestedclasses; +#pragma link C++ nestedtypedef; + +#pragma link C++ class ROOT::Math::Experimental::ScipyMinimizer; + +#endif //__CINT__ diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h new file mode 100644 index 0000000000000..00c85cf01b806 --- /dev/null +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -0,0 +1,173 @@ +// @(#)root/math/scipy:$Id$ +// Author: Omar.Zapata@cern.ch 2023 + +/************************************************************************* + * Copyright (C) 1995-2022, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +#ifndef ROOT_Math_ScipyMinimizer +#define ROOT_Math_ScipyMinimizer + +#include "Math/Minimizer.h" +#include "Math/MinimizerOptions.h" + +#include "Math/IFunctionfwd.h" + +#include "Math/IParamFunctionfwd.h" + +#include "Math/BasicMinimizer.h" + +#include "Rtypes.h" +#include "TString.h" + +#include +#include +#include + +#ifndef PyObject_HEAD +struct _object; +typedef _object PyObject; +#define Py_single_input 256 +#endif + +namespace ROOT { + +namespace Math { + +class GenAlgoOptions; + +namespace Experimental { +/** + enumeration specifying the types of Scipy solvers + @ingroup MultiMin +*/ + +//_____________________________________________________________________________________ +/** + \class ScipyMinimizer + ScipyMinimizer class. + Scipy minimizer implementation using Python C API that supports several methods such as + Nelder-Mead, L-BFGS-B, Powell, CG, BFGS, TNC, COBYLA, SLSQP, trust-constr, + Newton-CG, dogleg, trust-ncg, trust-exact and trust-krylov. + + It supports the Jacobian (Gradients), Hessian and bounds for the variables. + + Support for constraint functions will be implemented in the next releases. + You can find a macro example in the folder $ROOTSYS/tutorial/fit/scipy.C + + To provide extra options to the minimizer, you can use the class GenAlgoOptions + and the method SetExtraOptions(). + Example: + ``` + ROOT::Math::GenAlgoOptions l_bfgs_b_opt; + l_bfgs_b_opt.SetValue("gtol", 1e-3); + l_bfgs_b_opt.SetValue("ftol", 1e7); + minimizer->SetExtraOptions(l_bfgs_b_opt); + ``` + + See Scipy doc + from more info on the Scipy minimization algorithms. + + @ingroup MultiMin +*/ + +class ScipyMinimizer : public BasicMinimizer { +private: + PyObject *fMinimize; + PyObject *fTarget; + PyObject *fJacobian; + PyObject *fHessian; + PyObject *fBoundsMod; + PyObject *fConstraintsList; /// contraints functions + GenAlgoOptions *fExtraOpts; + std::function, double *)> fHessianFunc; + unsigned int fConstN; + unsigned int fCalls; + +protected: + PyObject *fGlobalNS; + PyObject *fLocalNS; + +public: + /** + Default constructor + */ + ScipyMinimizer(); + /** + Constructor with a string giving name of algorithm + */ + ScipyMinimizer(const char *type); + + /** + Destructor + */ + virtual ~ScipyMinimizer(); + + /** + Python initialization + */ + void PyInitialize(); + + /** + Checks if Python was initialized + */ + int PyIsInitialized(); + + /* + Python finalization + */ + void PyFinalize(); + /* + Python code execution + */ + void PyRunString(TString code, TString errorMessage = "Failed to run python code", int start = Py_single_input); + + /* + Number of function calls + */ + virtual unsigned int NCalls() const override; + + /* + Method to add Constraint function, + multiples constraints functions can be added. + type have to be a string "eq" or "ineq" where + eq (means Equal, then fun() = 0) + ineq (means that, it is to be non-negative. fun() >=0) + https://kitchingroup.cheme.cmu.edu/f19-06623/13-constrained-optimization.html + */ + virtual void AddConstraintFunction(std::function &)>, std::string type); + +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + ScipyMinimizer(const ScipyMinimizer &) : BasicMinimizer() {} + + /** + Get extra options from IOptions + */ + void SetExtraOptions(); + +public: + /// method to perform the minimization + virtual bool Minimize() override; + + virtual void SetHessianFunction(std::function, double *)>) override; + +protected: + ClassDef(ScipyMinimizer, 0) // +}; + +} // end namespace Experimental + +} // end namespace Math + +} // end namespace ROOT + +#endif /* ROOT_Math_ScipyMinimizer */ diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx new file mode 100644 index 0000000000000..8e130dfa04120 --- /dev/null +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -0,0 +1,436 @@ +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include // Needs to be included first to avoid redefinition of _POSIX_C_SOURCE +#include +#include +#include +#include +#include +#include +#include +#include + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#include +#include + +using namespace ROOT; +using namespace ROOT::Math; +using namespace ROOT::Math::Experimental; +using namespace ROOT::Fit; + +/// function wrapper for the function to be minimized +const ROOT::Math::IMultiGenFunction *gFunction = nullptr; +/// function wrapper for the gradient of the function to be minimized +const ROOT::Math::IMultiGradFunction *gGradFunction = nullptr; +/// function wrapper for the hessian of the function to be minimized +std::function &, double *)> gfHessianFunction; + +/// simple function for debugging +#define PyPrint(pyo) PyObject_Print(pyo, stdout, Py_PRINT_RAW) + +/// function to wrap into Python the C/C++ target function to be minimized +PyObject *target_function(PyObject * /*self*/, PyObject *args) +{ + PyArrayObject *arr = (PyArrayObject *)PyTuple_GetItem(args, 0); + + auto params = (double *)PyArray_DATA(arr); + auto r = (*gFunction)(params); + + return PyFloat_FromDouble(r); +}; + +/// function to wrap into Python the C/C++ jacobian function +PyObject *jac_function(PyObject * /*self*/, PyObject *args) +{ + PyArrayObject *arr = (PyArrayObject *)PyTuple_GetItem(args, 0); + + uint size = PyArray_SIZE(arr); + auto params = (double *)PyArray_DATA(arr); + double *values = new double[size]; + gGradFunction->Gradient(params, values); + npy_intp dims[1] = {size}; + PyObject *py_array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); + return py_array; +}; + +/// function to wrap into Python the C/C++ hessian function +PyObject *hessian_function(PyObject * /*self*/, PyObject *args) +{ + PyArrayObject *arr = (PyArrayObject *)PyTuple_GetItem(args, 0); + + uint size = PyArray_SIZE(arr); + auto params = (double *)PyArray_DATA(arr); + double *values = new double[size * size]; + std::vector x(params, params + size); + gfHessianFunction(x, values); + npy_intp dims[2] = {size, size}; + PyObject *py_array = PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, values); + return py_array; +}; + +//_______________________________________________________________________ +ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() +{ + fCalls = 0; + fOptions.SetMinimizerType("Scipy"); + fOptions.SetMinimizerAlgorithm("L-BFGS-B"); + PyInitialize(); + fHessianFunc = nullptr; + fConstraintsList = PyList_New(0); + fConstN = 0; + fExtraOpts = nullptr; +} + +//_______________________________________________________________________ +ScipyMinimizer::ScipyMinimizer(const char *type) +{ + fCalls = 0; + fOptions.SetMinimizerType("Scipy"); + fOptions.SetMinimizerAlgorithm(type); + PyInitialize(); + fHessianFunc = nullptr; + fConstraintsList = PyList_New(0); + fConstN = 0; + fExtraOpts = nullptr; +} + +//_______________________________________________________________________ +void ScipyMinimizer::PyInitialize() +{ + static PyObject *ParamsError; + static PyMethodDef ParamsMethods[] = { + {"target_function", target_function, METH_VARARGS, "Target function to minimize."}, + {"jac_function", jac_function, METH_VARARGS, "Jacobian function."}, + {"hessian_function", hessian_function, METH_VARARGS, "Hessianfunction."}, + {NULL, NULL, 0, NULL} /* Sentinel */ + }; + + static struct PyModuleDef paramsmodule = { + PyModuleDef_HEAD_INIT, + "params", /* name of module */ + "ROOT Scipy parameters", /* module documentation, may be NULL */ + -1, /* size of per-interpreter state of the module, + or -1 if the module keeps state in global variables. */ + ParamsMethods, + NULL, /* m_slots */ + NULL, /* m_traverse */ + 0, /* m_clear */ + NULL /* m_free */ + }; + + auto PyInit_params = [](void) -> PyObject * { + PyObject *m; + + m = PyModule_Create(¶msmodule); + if (m == NULL) + return NULL; + + ParamsError = PyErr_NewException("params.error", NULL, NULL); + Py_XINCREF(ParamsError); + if (PyModule_AddObject(m, "error", ParamsError) < 0) { + Py_XDECREF(ParamsError); + Py_CLEAR(ParamsError); + Py_DECREF(m); + return NULL; + } + + return m; + }; + if (PyImport_AppendInittab("params", PyInit_params) == -1) { + MATH_ERROR_MSG("ScipyMinimizer::Minimize", "could not extend in-built modules table"); + exit(1); + } + + bool pyIsInitialized = PyIsInitialized(); + if (!pyIsInitialized) { + Py_Initialize(); // Python initialization + } + fLocalNS = PyDict_New(); + fGlobalNS = PyDict_New(); + + if (!pyIsInitialized) { + _import_array(); // Numpy initialization + } + // Scipy initialization + + PyRunString("from scipy.optimize import minimize, Bounds"); + fMinimize = PyDict_GetItemString(fLocalNS, "minimize"); + fBoundsMod = PyDict_GetItemString(fLocalNS, "Bounds"); + PyRunString("from params import target_function, jac_function, hessian_function"); + fTarget = PyDict_GetItemString(fLocalNS, "target_function"); + fJacobian = PyDict_GetItemString(fLocalNS, "jac_function"); + fHessian = PyDict_GetItemString(fLocalNS, "hessian_function"); +} + +//_______________________________________________________________________ +// Finalize Python interpreter +void ScipyMinimizer::PyFinalize() +{ + Py_Finalize(); +} + +//_______________________________________________________________________ +int ScipyMinimizer::PyIsInitialized() +{ + return Py_IsInitialized(); +} + +//_______________________________________________________________________ +ScipyMinimizer::~ScipyMinimizer() +{ + if (fMinimize) + Py_DECREF(fMinimize); + if (fBoundsMod) + Py_DECREF(fBoundsMod); +} + +//_______________________________________________________________________ +void ScipyMinimizer::SetExtraOptions() +{ + auto constExtraOpts = dynamic_cast(fOptions.ExtraOptions()); + fExtraOpts = const_cast(constExtraOpts); +} + +//_______________________________________________________________________ +bool ScipyMinimizer::Minimize() +{ + (gFunction) = ObjFunction(); + (gGradFunction) = GradObjFunction(); + gfHessianFunction = fHessianFunc; + if (gGradFunction == nullptr) { + fJacobian = Py_None; + } + if (gfHessianFunction == nullptr) { + fHessian = Py_None; + } + auto method = fOptions.MinimizerAlgorithm(); + PyObject *pyoptions = PyDict_New(); + SetExtraOptions(); + if (fExtraOpts) { + for (std::string key : fExtraOpts->GetAllRealKeys()) { + double value = 0; + fExtraOpts->GetRealValue(key.c_str(), value); + PyDict_SetItemString(pyoptions, key.c_str(), PyFloat_FromDouble(value)); + } + for (std::string key : fExtraOpts->GetAllIntKeys()) { + int value = 0; + fExtraOpts->GetIntValue(key.c_str(), value); + PyDict_SetItemString(pyoptions, key.c_str(), PyLong_FromLong(value)); + } + for (std::string key : fExtraOpts->GetAllNamedKeys()) { + std::string value = ""; + fExtraOpts->GetNamedValue(key.c_str(), value); + PyDict_SetItemString(pyoptions, key.c_str(), PyUnicode_FromString(value.c_str())); + } + } + PyDict_SetItemString(pyoptions, "maxiter", PyLong_FromLong(MaxIterations())); + if (PrintLevel() > 0) { + PyDict_SetItemString(pyoptions, "disp", Py_True); + } + std::cout << "=== Scipy Minimization" << std::endl; + std::cout << "=== Method: " << method << std::endl; + std::cout << "=== Initial value: ("; + for (uint i = 0; i < NDim(); i++) { + std::cout << X()[i]; + if (i < NDim() - 1) + std::cout << ","; + } + std::cout << ")" << std::endl; + + double *values = const_cast(X()); + npy_intp dims[1] = {NDim()}; + PyObject *x0 = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); + + PyObject *pybounds_args = PyTuple_New(2); + PyObject *pylimits_lower = PyList_New(NDim()); + PyObject *pylimits_upper = PyList_New(NDim()); + bool foundBounds = false; + for (unsigned int i = 0; i < NDim(); i++) { + ParameterSettings varsettings; + + if (GetVariableSettings(i, varsettings)) { + if (varsettings.HasLowerLimit()) { + foundBounds = true; + PyList_SetItem(pylimits_lower, i, PyFloat_FromDouble(varsettings.LowerLimit())); + } else { + PyList_SetItem(pylimits_lower, i, PyFloat_FromDouble(-NPY_INFINITY)); + } + if (varsettings.HasUpperLimit()) { + foundBounds = true; + PyList_SetItem(pylimits_upper, i, PyFloat_FromDouble(varsettings.UpperLimit())); + } else { + PyList_SetItem(pylimits_upper, i, PyFloat_FromDouble(NPY_INFINITY)); + } + } else { + MATH_ERROR_MSG("ScipyMinimizer::Minimize", Form("Variable index = %d not found", i)); + } + } + PyObject *pybounds = Py_None; + if (foundBounds) { + PyTuple_SetItem(pybounds_args, 0, pylimits_lower); + PyTuple_SetItem(pybounds_args, 1, pylimits_upper); + pybounds = PyObject_CallObject(fBoundsMod, pybounds_args); + } + + // minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, + // callback=None, options=None) + PyObject *args = Py_BuildValue("(OO)", fTarget, x0); + PyObject *kw = + Py_BuildValue("{s:s,s:O,,s:O,s:O,s:O,s:d,s:O}", "method", method.c_str(), "jac", fJacobian, "hess", fHessian, + "bounds", pybounds, "constraints", fConstraintsList, "tol", Tolerance(), "options", pyoptions); + if (PrintLevel() > 0) { + std::cout << "========Minimizer Parameters========\n"; + PyPrint(kw); + std::cout << "====================================\n"; + } + PyObject *result = PyObject_Call(fMinimize, args, kw); + if (result == NULL) { + PyErr_Print(); + return false; + } + if (PrintLevel() > 0) { + std::cout << "========Minimizer Results========\n"; + PyPrint(result); + std::cout << "=================================\n"; + } + Py_DECREF(pybounds); + Py_DECREF(args); + Py_DECREF(kw); + Py_DECREF(x0); + Py_DECREF(fConstraintsList); + fConstraintsList = PyList_New(0); + + // if the minimization works + PyObject *pstatus = PyObject_GetAttrString(result, "status"); + int status = PyLong_AsLong(pstatus); + Py_DECREF(pstatus); + + PyObject *psuccess = PyObject_GetAttrString(result, "success"); + bool success = PyLong_AsLong(psuccess); + Py_DECREF(psuccess); + + // the x values for the minimum + PyArrayObject *pyx = (PyArrayObject *)PyObject_GetAttrString(result, "x"); + const double *x = (const double *)PyArray_DATA(pyx); + Py_DECREF(pyx); + + // number of function evaluations + PyObject *pynfev = PyObject_GetAttrString(result, "nfev"); + long nfev = PyLong_AsLong(pynfev); + Py_DECREF(pynfev); + + PyObject *pymessage = PyObject_GetAttrString(result, "message"); + const char *message = (const char *)PyUnicode_DATA(pymessage); + Py_DECREF(pymessage); + + SetFinalValues(x); + auto obj_value = (*gFunction)(x); + SetMinValue(obj_value); + fCalls = nfev; // number of function evaluations + std::cout << "=== Success: " << success << std::endl; + std::cout << "=== Status: " << status << std::endl; + std::cout << "=== Message: " << message << std::endl; + std::cout << "=== Function calls: " << nfev << std::endl; + if (success) + fStatus = 0; + else + fStatus = status; // suggested by Lorenzo. + + Py_DECREF(result); + return success; +} + +//_______________________________________________________________________ +void ScipyMinimizer::PyRunString(TString code, TString errorMessage, int start) +{ + auto fPyReturn = PyRun_String(code, start, fGlobalNS, fLocalNS); + if (!fPyReturn) { + auto msg = errorMessage + Form(" \ncode = \"%s\"", code.Data()); + MATH_ERROR_MSG("ScipyMinimizer::PyRunString", msg.Data()); + PyErr_Print(); + exit(1); + } +} + +//_______________________________________________________________________ +void ScipyMinimizer::SetHessianFunction(std::function, double *)> func) +{ + fHessianFunc = func; +} + +//_______________________________________________________________________ +unsigned int ScipyMinimizer::NCalls() const +{ + return fCalls; +} + +//_______________________________________________________________________ +void ScipyMinimizer::AddConstraintFunction(std::function &)> func, std::string type) +{ + if (type != "eq" && type != "ineq") { + MATH_ERROR_MSG("ScipyMinimizer::AddConstraintFunction", + Form("Error in constraint type %s, it have to be \"eq\" or \"ineq\"", type.c_str())); + exit(1); + } + static std::function &)> cfunt = func; + auto const_function = [](PyObject * /*self*/, PyObject *args) -> PyObject * { + PyArrayObject *arr = (PyArrayObject *)PyTuple_GetItem(args, 0); + + uint size = PyArray_SIZE(arr); + auto params = (double *)PyArray_DATA(arr); + std::vector x(params, params + size); + auto r = cfunt(x); + return PyFloat_FromDouble(r); + }; + + static const char *name = Form("const_function%d", fConstN); + static const char *name_error = Form("const_function%d.error", fConstN); + static PyObject *ConstError; + static PyMethodDef ConstMethods[] = { + {name, const_function, METH_VARARGS, "Constraint function to minimize."}, {NULL, NULL, 0, NULL} /* Sentinel */ + }; + static struct PyModuleDef constmodule = { + PyModuleDef_HEAD_INIT, + name, /* name of module */ + "ROOT Scipy parameters", /* module documentation, may be NULL */ + -1, /* size of per-interpreter state of the module, + or -1 if the module keeps state in global variables. */ + ConstMethods, + NULL, /* m_slots */ + NULL, /* m_traverse */ + 0, /* m_clear */ + NULL /* m_free */ + }; + + auto PyInit_const = [](void) -> PyObject * { + PyObject *m; + + m = PyModule_Create(&constmodule); + if (m == NULL) + return NULL; + ConstError = PyErr_NewException(name_error, NULL, NULL); + Py_XINCREF(ConstError); + if (PyModule_AddObject(m, "error", ConstError) < 0) { + Py_XDECREF(ConstError); + Py_CLEAR(ConstError); + Py_DECREF(m); + return NULL; + } + return m; + }; + PyImport_AddModule(name); + PyObject *module = PyInit_const(); + PyObject *sys_modules = PyImport_GetModuleDict(); + PyDict_SetItemString(sys_modules, name, module); + + PyRunString(Form("from %s import %s", name, name)); + PyObject *pyconstfun = PyDict_GetItemString(fLocalNS, name); + + PyObject *pyconst = PyDict_New(); + PyDict_SetItemString(pyconst, "type", PyUnicode_FromString(type.c_str())); + PyDict_SetItemString(pyconst, "fun", pyconstfun); + PyList_Append(fConstraintsList, pyconst); + + Py_DECREF(ConstError); + Py_DECREF(module); +} diff --git a/math/scipy/test/CMakeLists.txt b/math/scipy/test/CMakeLists.txt new file mode 100644 index 0000000000000..2cd2c14933e7b --- /dev/null +++ b/math/scipy/test/CMakeLists.txt @@ -0,0 +1,11 @@ +project(scipy-tests) +find_package(ROOT REQUIRED) + +include_directories(${ROOT_INCLUDE_DIRS}) + +set(Libraries Core RIO Net Hist Graf Graf3d Gpad Tree + Rint Postscript Matrix MathCore Thread Scipy) + + +include_directories(${CMAKE_SOURCE_DIR}) +ROOT_ADD_GTEST(scipy-minimizer testScipyMinimizer.cxx LIBRARIES ${Libraries}) diff --git a/math/scipy/test/testScipyMinimizer.cxx b/math/scipy/test/testScipyMinimizer.cxx new file mode 100644 index 0000000000000..153f8d49f9185 --- /dev/null +++ b/math/scipy/test/testScipyMinimizer.cxx @@ -0,0 +1,129 @@ +// @(#)root/math/scipy $Id$ +// Author: Omar Zapata 2023 + +#include + +#include "Math/ScipyMinimizer.h" +#include "Math/Functor.h" +#include "Math/Factory.h" +#include "Math/MultiNumGradFunction.h" +#include "Math/FitMethodFunction.h" + +// target function +double RosenBrock(const double *xx) +{ + const Double_t x = xx[0]; + const Double_t y = xx[1]; + const Double_t tmp1 = y - x * x; + const Double_t tmp2 = 1 - x; + return 100 * tmp1 * tmp1 + tmp2 * tmp2; +} + +// gradient function(jacobian) +double RosenBrockGrad(const double *x, unsigned int ipar) +{ + if (ipar == 0) + return -2 * (1 - x[0]) + 200 * (x[1] - x[0] * x[0]) * (-2 * x[0]); + else + return 200 * (x[1] - x[0] * x[0]); +} +// Hessian function +bool RosenBrockHessian(std::span xx, double *hess) +{ + const double x = xx[0]; + const double y = xx[1]; + + hess[0] = 1200 * x * x - 400 * y + 2; + hess[1] = -400 * x; + hess[2] = -400 * x; + hess[3] = 200; + + return true; +} + +// https://de.mathworks.com/help/optim/ug/example-nonlinear-constrained-minimization.html +// constraint function +double ConstRosenBrock(const std::vector &x) +{ + return x[0] * x[0] - x[1] * x[1] + 1; +} +// NOTE: "dogleg" is problematic, requires better tuning of the parameters +std::string methods[] = {"Nelder-Mead", "L-BFGS-B", "Powell", "CG", "BFGS", + "TNC", "COBYLA", "SLSQP", "trust-constr", "Newton-CG", + "trust-ncg", "trust-exact", "trust-krylov"}; + +// Testing fit using class with gradient +class ScipyFitClass : public ::testing::Test { +public: + ROOT::Math::Experimental::ScipyMinimizer *minimizer = nullptr; + + ScipyFitClass() = default; + void Fit(const char *method, bool useConstraint = false) + { + ROOT::Math::GradFunctor f(&RosenBrock, &RosenBrockGrad, 2); + + minimizer = new ROOT::Math::Experimental::ScipyMinimizer(method); + minimizer->SetMaxFunctionCalls(1000000); + minimizer->SetMaxIterations(100000); + minimizer->SetTolerance(0.001); + + double step[2] = {0.01, 0.01}; + double variable[2] = {0.1, 1.2}; + + minimizer->SetFunction(f); + minimizer->SetHessianFunction(RosenBrockHessian); + + // Set the free variables to be minimized! + minimizer->SetVariable(0, "x", variable[0], step[0]); + minimizer->SetVariable(1, "y", variable[1], step[1]); + if (useConstraint) + minimizer->AddConstraintFunction(ConstRosenBrock, "ineq"); + + minimizer->Minimize(); + } +}; + +TEST_F(ScipyFitClass, Fit) +{ + for (const std::string &text : methods) { + Fit(text.c_str(), false); + EXPECT_EQ(1000000, minimizer->MaxFunctionCalls()); + + EXPECT_EQ(100000, minimizer->MaxIterations()); + + EXPECT_EQ(0.001, minimizer->Tolerance()); + + auto step = minimizer->StepSizes(); + EXPECT_EQ(0.01, step[0]); + EXPECT_EQ(0.01, step[1]); + + auto v1name = minimizer->VariableName(0); + auto v2name = minimizer->VariableName(1); + EXPECT_EQ("x", v1name); + EXPECT_EQ("y", v2name); + + EXPECT_EQ(0, minimizer->VariableIndex("x")); + EXPECT_EQ(1, minimizer->VariableIndex("y")); + + auto opts = minimizer->Options(); + + EXPECT_EQ("Scipy", opts.MinimizerType()); + EXPECT_EQ(text, opts.MinimizerAlgorithm()); + + auto x = minimizer->X(); + ASSERT_NEAR(1, x[0], 0.5); + ASSERT_NEAR(1, x[1], 0.5); + ASSERT_NEAR(0, RosenBrock(x), 0.5); + } +} + +TEST_F(ScipyFitClass, FitContraint) // using constraint function +{ + for (const std::string &text : methods) { + Fit(text.c_str(), true); + auto x = minimizer->X(); + ASSERT_NEAR(1, x[0], 0.5); + ASSERT_NEAR(1, x[1], 0.5); + ASSERT_NEAR(0, RosenBrock(x), 0.5); + } +} diff --git a/tutorials/math/fit/scipy.C b/tutorials/math/fit/scipy.C new file mode 100644 index 0000000000000..d1cbc2a543cb6 --- /dev/null +++ b/tutorials/math/fit/scipy.C @@ -0,0 +1,81 @@ +#include "Math/ScipyMinimizer.h" +#include "Math/MultiNumGradFunction.h" +#include +#include "Math/Functor.h" +#include +#include "Math/MinimizerOptions.h" +#include "TStopwatch.h" + +double RosenBrock(const double *xx) +{ + const Double_t x = xx[0]; + const Double_t y = xx[1]; + const Double_t tmp1 = y - x * x; + const Double_t tmp2 = 1 - x; + return 100 * tmp1 * tmp1 + tmp2 * tmp2; +} + +double RosenBrockGrad(const double *x, unsigned int ipar) +{ + if (ipar == 0) + return -2 * (1 - x[0]) + 200 * (x[1] - x[0] * x[0]) * (-2 * x[0]); + else + return 200 * (x[1] - x[0] * x[0]); +} + +bool RosenBrockHessian(std::span xx, double *hess) +{ + const double x = xx[0]; + const double y = xx[1]; + + hess[0] = 1200 * x * x - 400 * y + 2; + hess[1] = -400 * x; + hess[2] = -400 * x; + hess[3] = 200; + + return true; +} + +// methods that requires hessian to work "dogleg", "trust-ncg","trust-exact","trust-krylov" +using namespace std; +int scipy() +{ + + std::string methods[] = {"Nelder-Mead", "L-BFGS-B", "Powell", "CG", "BFGS", + "TNC", "COBYLA", "SLSQP", "trust-constr", "Newton-CG", + "dogleg", "trust-ncg", "trust-exact", "trust-krylov"}; + TStopwatch t; + for (const std::string &text : methods) { + ROOT::Math::Experimental::ScipyMinimizer minimizer(text.c_str()); + minimizer.SetMaxFunctionCalls(1000000); + minimizer.SetMaxIterations(100000); + minimizer.SetTolerance(1e-3); + ROOT::Math::GradFunctor f(&RosenBrock, &RosenBrockGrad, 2); + double step[2] = {0.01, 0.01}; + double variable[2] = {-1.2, 1.0}; + + minimizer.SetFunction(f); + minimizer.SetHessianFunction(RosenBrockHessian); + + // variables to be minimized! + minimizer.SetVariable(0, "x", variable[0], step[0]); + minimizer.SetVariable(1, "y", variable[1], step[1]); + minimizer.SetVariableLimits(0, -2.0, 2.0); + minimizer.SetVariableLimits(1, -2.0, 2.0); + + t.Reset(); + t.Start(); + minimizer.Minimize(); + t.Stop(); + const double *xs = minimizer.X(); + cout << "Minimum: f(" << xs[0] << "," << xs[1] << "): " << RosenBrock(xs) << endl; + cout << "Cpu Time (sec) = " << t.CpuTime() << endl << "Real Time (sec) = " << t.RealTime() << endl; + cout << endl << "===============" << endl; + } + return 0; +} + +int main() +{ + return scipy(); +}