From 8fe3e78f3f0460e2771e76a5d16c9f0f7546c724 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Fri, 16 Dec 2022 03:32:20 -0500 Subject: [PATCH 01/25] Math: working in the initial prototype for Scipy optimizers --- .../ROOT@@Math@@Minimizer/P100_Scipy.C | 5 + math/CMakeLists.txt | 4 + math/scipy/CMakeLists.txt | 20 ++ math/scipy/inc/Math/LinkDef.h | 11 + math/scipy/inc/Math/ScipyMinimizer.h | 145 +++++++++++++ math/scipy/src/ScipyMinimizer.cxx | 194 ++++++++++++++++++ 6 files changed, 379 insertions(+) create mode 100644 etc/plugins/ROOT@@Math@@Minimizer/P100_Scipy.C create mode 100644 math/scipy/CMakeLists.txt create mode 100644 math/scipy/inc/Math/LinkDef.h create mode 100644 math/scipy/inc/Math/ScipyMinimizer.h create mode 100644 math/scipy/src/ScipyMinimizer.cxx 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..2a72ef78b012c 100644 --- a/math/CMakeLists.txt +++ b/math/CMakeLists.txt @@ -34,4 +34,8 @@ if(r) add_subdirectory(rtools) endif() +if(scipy AND NUMPY_FOUND) + add_subdirectory(scipy) +endif() + add_subdirectory(vecops) diff --git a/math/scipy/CMakeLists.txt b/math/scipy/CMakeLists.txt new file mode 100644 index 0000000000000..69ab2a5e88efa --- /dev/null +++ b/math/scipy/CMakeLists.txt @@ -0,0 +1,20 @@ +############################################################################ +# 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 ${PYTHON_LIBRARIES_Development_Main} + 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..af9a5af4df464 --- /dev/null +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -0,0 +1,145 @@ +// @(#)root/math/scipy:$Id$ +// Author: Omar.Zapata@cern.ch 2022 + +/************************************************************************* + * 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 + +#ifndef PyObject_HEAD +struct _object; +typedef _object PyObject; +#define Py_single_input 256 +#endif + +namespace ROOT { + +namespace Math { + +namespace Experimental { +/** + enumeration specifying the types of Scipy solvers + @ingroup MultiMin +*/ + + +//_____________________________________________________________________________________ +/** + * \class ScipyMinimizer + * ScipyMinimizer class. + * Implementation for Scipy ... TODO + + See Scipy doc + from more info on the Scipy minimization algorithms. + + @ingroup MultiMin +*/ + +class ScipyMinimizer : public BasicMinimizer { +private: + PyObject *fMinimize; + PyObject *fTarget; +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 eval function + */ + PyObject *Eval(TString code); + + /** + Python initialization + */ + void PyInitialize(); + + /** + Checks if Python was initialized + */ + int PyIsInitialized(); + + /* + Python finalization + */ + void PyFinalize(); + void PyRunString(TString code, TString errorMessage="Failed to run python code", int start=Py_single_input); + void LoadWrappers(); +private: + // usually copying is non trivial, so we make this unaccessible + + /** + Copy constructor + */ + ScipyMinimizer(const ScipyMinimizer &) : BasicMinimizer() {} + + +public: + /// set the function to minimize + //virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func); + + /// set the function to minimize + //virtual void SetFunction(const ROOT::Math::IMultiGradFunction &func) { BasicMinimizer::SetFunction(func); } + + /// method to perform the minimization + virtual bool Minimize(); + + /// return expected distance reached from the minimum + virtual double Edm() const { return 0; } // not impl. } + + /// minimizer provides error and error matrix + virtual bool ProvidesError() const { return false; } + + /// return errors at the minimum + virtual const double *Errors() const { return 0; } + + /** return covariance matrices elements + if the variable is fixed the matrix is zero + The ordering of the variables is the same as in errors + */ + virtual double CovMatrix(unsigned int, unsigned int) const { return 0; } +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..8389e4e4345e8 --- /dev/null +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -0,0 +1,194 @@ +#include // Needs to be included first to avoid redefinition of _POSIX_C_SOURCE +#include +#include +#include +#include +#include +#include + +#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION +#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; +/// function wrapper for the gradient of the function to be minimized +const ROOT::Math::IMultiGradFunction *gGradFunction; + +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); +}; + + +//_______________________________________________________________________ +ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() +{ + fOptions.SetMinimizerType("Scipy"); + fOptions.SetMinimizerAlgorithm("lbfgsb"); + if (!PyIsInitialized()) { + PyInitialize(); + } + +} + + + +//_______________________________________________________________________ +ScipyMinimizer::ScipyMinimizer(const char *type) +{ + fOptions.SetMinimizerType("Scipy"); + fOptions.SetMinimizerAlgorithm(type); + if (!PyIsInitialized()) { + PyInitialize(); + } +} + +//_______________________________________________________________________ +void ScipyMinimizer::PyInitialize() +{ + static PyObject *ParamsError; + static PyMethodDef ParamsMethods[] = { + {"target_function", target_function, METH_VARARGS, + "Target function to minimize."}, + {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 + }; + + 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"); + fMinimize = PyDict_GetItemString(fLocalNS, "minimize"); + PyRunString("from params import target_function"); + fTarget = PyDict_GetItemString(fLocalNS, "target_function"); +} + +//_______________________________________________________________________ +// Finalize Python interpreter +void ScipyMinimizer::PyFinalize() +{ + if (fMinimize) Py_DECREF(fMinimize); + Py_Finalize(); +} + +//_______________________________________________________________________ +int ScipyMinimizer::PyIsInitialized() +{ + if (!Py_IsInitialized()) return kFALSE; + return kTRUE; +} + +//_______________________________________________________________________ +ScipyMinimizer::~ScipyMinimizer() +{ +} + +//_______________________________________________________________________ +bool ScipyMinimizer::Minimize() +{ + (gFunction)= ObjFunction(); + (gGradFunction) = GradObjFunction(); + auto method = fOptions.MinimizerAlgorithm(); + + std::cout<<"=== Scipy Minimization"<(X()); + npy_intp dims[1] = { NDim()}; + PyObject *py_array= PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); + + PyObject *pargs=PyTuple_New(0); + + auto pyvalues = Py_BuildValue("(OOOs)",fTarget,py_array,pargs,method.c_str()); + + PyObject *result = PyObject_CallObject(fMinimize,pyvalues); + Py_DECREF(pyvalues); + Py_DECREF(py_array); + + + //if the minimization works + PyObject *pstatus = PyObject_GetAttrString(result,"status"); + bool status = PyBool_Check(pstatus ); + Py_DECREF(pstatus); + + //the x values for the minimum + PyArrayObject *pyx = (PyArrayObject *)PyObject_GetAttrString(result,"x"); + const double *x =(const double *)PyArray_DATA(pyx); + Py_DECREF(pyx); + + SetFinalValues(x); + auto obj_value = (*gFunction)(x); + SetMinValue(obj_value); + + return status; +} + +//_______________________________________________________________________ +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); + } +} From 12476281078bc2e0b082cea79c0dbf85c07003e9 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Fri, 16 Dec 2022 03:56:41 -0500 Subject: [PATCH 02/25] Math Scipy: extracted informatio for the results --- math/scipy/inc/Math/ScipyMinimizer.h | 22 +++--- math/scipy/src/ScipyMinimizer.cxx | 105 ++++++++++++++------------- 2 files changed, 65 insertions(+), 62 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index af9a5af4df464..11162cb42f94c 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -43,7 +43,6 @@ namespace Experimental { @ingroup MultiMin */ - //_____________________________________________________________________________________ /** * \class ScipyMinimizer @@ -58,11 +57,13 @@ namespace Experimental { class ScipyMinimizer : public BasicMinimizer { private: - PyObject *fMinimize; - PyObject *fTarget; + PyObject *fMinimize; + PyObject *fTarget; + protected: - PyObject *fGlobalNS; - PyObject *fLocalNS; + PyObject *fGlobalNS; + PyObject *fLocalNS; + public: /** Default constructor @@ -91,14 +92,15 @@ class ScipyMinimizer : public BasicMinimizer { /** Checks if Python was initialized */ - int PyIsInitialized(); + int PyIsInitialized(); /* Python finalization */ void PyFinalize(); - void PyRunString(TString code, TString errorMessage="Failed to run python code", int start=Py_single_input); + void PyRunString(TString code, TString errorMessage = "Failed to run python code", int start = Py_single_input); void LoadWrappers(); + private: // usually copying is non trivial, so we make this unaccessible @@ -107,13 +109,12 @@ class ScipyMinimizer : public BasicMinimizer { */ ScipyMinimizer(const ScipyMinimizer &) : BasicMinimizer() {} - public: /// set the function to minimize - //virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func); + // virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func); /// set the function to minimize - //virtual void SetFunction(const ROOT::Math::IMultiGradFunction &func) { BasicMinimizer::SetFunction(func); } + // virtual void SetFunction(const ROOT::Math::IMultiGradFunction &func) { BasicMinimizer::SetFunction(func); } /// method to perform the minimization virtual bool Minimize(); @@ -132,6 +133,7 @@ class ScipyMinimizer : public BasicMinimizer { The ordering of the variables is the same as in errors */ virtual double CovMatrix(unsigned int, unsigned int) const { return 0; } + protected: ClassDef(ScipyMinimizer, 0) // }; diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 8389e4e4345e8..a592611618d44 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -14,13 +14,12 @@ 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; /// function wrapper for the gradient of the function to be minimized const ROOT::Math::IMultiGradFunction *gGradFunction; -PyObject * target_function(PyObject */*self*/, PyObject *args) +PyObject *target_function(PyObject * /*self*/, PyObject *args) { PyArrayObject *arr = (PyArrayObject *)PyTuple_GetItem(args, 0); @@ -30,7 +29,6 @@ PyObject * target_function(PyObject */*self*/, PyObject *args) return PyFloat_FromDouble(r); }; - //_______________________________________________________________________ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() { @@ -39,11 +37,8 @@ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() if (!PyIsInitialized()) { PyInitialize(); } - } - - //_______________________________________________________________________ ScipyMinimizer::ScipyMinimizer(const char *type) { @@ -59,22 +54,17 @@ void ScipyMinimizer::PyInitialize() { static PyObject *ParamsError; static PyMethodDef ParamsMethods[] = { - {"target_function", target_function, METH_VARARGS, - "Target function to minimize."}, - {NULL, NULL, 0, NULL} /* Sentinel */ + {"target_function", target_function, METH_VARARGS, "Target function to minimize."}, + {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 - }; + 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}; - auto PyInit_params = [](void)->PyObject * - { + auto PyInit_params = [](void) -> PyObject * { PyObject *m; m = PyModule_Create(¶msmodule); @@ -99,7 +89,7 @@ void ScipyMinimizer::PyInitialize() bool pyIsInitialized = PyIsInitialized(); if (!pyIsInitialized) { - Py_Initialize(); // Python initialization + Py_Initialize(); // Python initialization } fLocalNS = PyDict_New(); fGlobalNS = PyDict_New(); @@ -107,7 +97,7 @@ void ScipyMinimizer::PyInitialize() if (!pyIsInitialized) { _import_array(); // Numpy initialization } - //Scipy initialization + // Scipy initialization PyRunString("from scipy.optimize import minimize"); fMinimize = PyDict_GetItemString(fLocalNS, "minimize"); PyRunString("from params import target_function"); @@ -118,76 +108,87 @@ void ScipyMinimizer::PyInitialize() // Finalize Python interpreter void ScipyMinimizer::PyFinalize() { - if (fMinimize) Py_DECREF(fMinimize); + if (fMinimize) + Py_DECREF(fMinimize); Py_Finalize(); } //_______________________________________________________________________ int ScipyMinimizer::PyIsInitialized() { - if (!Py_IsInitialized()) return kFALSE; + if (!Py_IsInitialized()) + return kFALSE; return kTRUE; } //_______________________________________________________________________ -ScipyMinimizer::~ScipyMinimizer() -{ -} +ScipyMinimizer::~ScipyMinimizer() {} //_______________________________________________________________________ bool ScipyMinimizer::Minimize() { - (gFunction)= ObjFunction(); + (gFunction) = ObjFunction(); (gGradFunction) = GradObjFunction(); auto method = fOptions.MinimizerAlgorithm(); - std::cout<<"=== Scipy Minimization"<(X()); - npy_intp dims[1] = { NDim()}; - PyObject *py_array= PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); + double *values = const_cast(X()); + npy_intp dims[1] = {NDim()}; + PyObject *py_array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); - PyObject *pargs=PyTuple_New(0); + PyObject *pargs = PyTuple_New(0); - auto pyvalues = Py_BuildValue("(OOOs)",fTarget,py_array,pargs,method.c_str()); + auto pyvalues = Py_BuildValue("(OOOs)", fTarget, py_array, pargs, method.c_str()); - PyObject *result = PyObject_CallObject(fMinimize,pyvalues); + PyObject *result = PyObject_CallObject(fMinimize, pyvalues); Py_DECREF(pyvalues); Py_DECREF(py_array); - - //if the minimization works - PyObject *pstatus = PyObject_GetAttrString(result,"status"); - bool status = PyBool_Check(pstatus ); + // if the minimization works + PyObject *pstatus = PyObject_GetAttrString(result, "status"); + bool status = PyBool_Check(pstatus); Py_DECREF(pstatus); - //the x values for the minimum - PyArrayObject *pyx = (PyArrayObject *)PyObject_GetAttrString(result,"x"); - const double *x =(const double *)PyArray_DATA(pyx); + // 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); + std::cout << "=== Status: " << status << std::endl; + std::cout << "=== Message: " << message << std::endl; + std::cout << "=== Function calls: " << nfev << std::endl; return status; } //_______________________________________________________________________ -void ScipyMinimizer::PyRunString(TString code, TString errorMessage, int start) { +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() ); + auto msg = errorMessage + Form(" \ncode = \"%s\"", code.Data()); + MATH_ERROR_MSG("ScipyMinimizer::PyRunString", msg.Data()); PyErr_Print(); exit(1); } From 19cd43a7a82385706494382388455195d38745b4 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Sat, 17 Dec 2022 05:13:34 -0500 Subject: [PATCH 03/25] implemented support for jacobian --- math/scipy/inc/Math/ScipyMinimizer.h | 5 ++++- math/scipy/src/ScipyMinimizer.cxx | 20 +++++++++++++++++--- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index 11162cb42f94c..38fc5c12a7a50 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -59,6 +59,7 @@ class ScipyMinimizer : public BasicMinimizer { private: PyObject *fMinimize; PyObject *fTarget; + PyObject *fJacobian; protected: PyObject *fGlobalNS; @@ -98,8 +99,10 @@ class ScipyMinimizer : public BasicMinimizer { Python finalization */ void PyFinalize(); + /* + Python code execution + */ void PyRunString(TString code, TString errorMessage = "Failed to run python code", int start = Py_single_input); - void LoadWrappers(); private: // usually copying is non trivial, so we make this unaccessible diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index a592611618d44..2694d58a0d303 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -29,6 +29,19 @@ PyObject *target_function(PyObject * /*self*/, PyObject *args) return PyFloat_FromDouble(r); }; +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[size]; + gGradFunction->Gradient(params, values); + npy_intp dims[1] = {size}; + PyObject *py_array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); + return py_array; +}; + //_______________________________________________________________________ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() { @@ -55,6 +68,7 @@ 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."}, {NULL, NULL, 0, NULL} /* Sentinel */ }; @@ -100,8 +114,9 @@ void ScipyMinimizer::PyInitialize() // Scipy initialization PyRunString("from scipy.optimize import minimize"); fMinimize = PyDict_GetItemString(fLocalNS, "minimize"); - PyRunString("from params import target_function"); + PyRunString("from params import target_function, jac_function"); fTarget = PyDict_GetItemString(fLocalNS, "target_function"); + fJacobian = PyDict_GetItemString(fLocalNS, "jac_function"); } //_______________________________________________________________________ @@ -130,7 +145,6 @@ bool ScipyMinimizer::Minimize() (gFunction) = ObjFunction(); (gGradFunction) = GradObjFunction(); auto method = fOptions.MinimizerAlgorithm(); - std::cout << "=== Scipy Minimization" << std::endl; std::cout << "=== Method: " << method << std::endl; std::cout << "=== Initial value: ("; @@ -147,7 +161,7 @@ bool ScipyMinimizer::Minimize() PyObject *pargs = PyTuple_New(0); - auto pyvalues = Py_BuildValue("(OOOs)", fTarget, py_array, pargs, method.c_str()); + auto pyvalues = Py_BuildValue("(OOOsO)", fTarget, py_array, pargs, method.c_str(), fJacobian); PyObject *result = PyObject_CallObject(fMinimize, pyvalues); Py_DECREF(pyvalues); From 58ea062d112ae718a03e287942a9772dda7d48cf Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Wed, 21 Dec 2022 02:25:07 -0500 Subject: [PATCH 04/25] Scipy: added support for extra options --- math/scipy/inc/Math/ScipyMinimizer.h | 27 ++++++------- math/scipy/src/ScipyMinimizer.cxx | 59 ++++++++++++++++++++++------ 2 files changed, 59 insertions(+), 27 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index 38fc5c12a7a50..e823591578f4b 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -21,6 +21,8 @@ #include "Math/BasicMinimizer.h" +#include "Math/GenAlgoOptions.h" + #include "Rtypes.h" #include "TString.h" @@ -37,6 +39,8 @@ namespace ROOT { namespace Math { +class GenAlgoOptions; + namespace Experimental { /** enumeration specifying the types of Scipy solvers @@ -60,6 +64,7 @@ class ScipyMinimizer : public BasicMinimizer { PyObject *fMinimize; PyObject *fTarget; PyObject *fJacobian; + GenAlgoOptions fExtraOpts; protected: PyObject *fGlobalNS; @@ -111,6 +116,7 @@ class ScipyMinimizer : public BasicMinimizer { Copy constructor */ ScipyMinimizer(const ScipyMinimizer &) : BasicMinimizer() {} + void SetAlgoExtraOptions(); public: /// set the function to minimize @@ -120,22 +126,13 @@ class ScipyMinimizer : public BasicMinimizer { // virtual void SetFunction(const ROOT::Math::IMultiGradFunction &func) { BasicMinimizer::SetFunction(func); } /// method to perform the minimization - virtual bool Minimize(); - - /// return expected distance reached from the minimum - virtual double Edm() const { return 0; } // not impl. } + virtual bool Minimize() override; - /// minimizer provides error and error matrix - virtual bool ProvidesError() const { return false; } - - /// return errors at the minimum - virtual const double *Errors() const { return 0; } - - /** return covariance matrices elements - if the variable is fixed the matrix is zero - The ordering of the variables is the same as in errors - */ - virtual double CovMatrix(unsigned int, unsigned int) const { return 0; } + template + void SetExtraOption(const char *key, T value) + { + fExtraOpts.SetValue(key, value); + } protected: ClassDef(ScipyMinimizer, 0) // diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 2694d58a0d303..c818521e1b4e8 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -3,6 +3,7 @@ #include #include #include +#include "Math/GenAlgoOptions.h" #include #include @@ -15,9 +16,11 @@ using namespace ROOT::Math::Experimental; using namespace ROOT::Fit; /// function wrapper for the function to be minimized -const ROOT::Math::IMultiGenFunction *gFunction; +const ROOT::Math::IMultiGenFunction *gFunction = nullptr; /// function wrapper for the gradient of the function to be minimized -const ROOT::Math::IMultiGradFunction *gGradFunction; +const ROOT::Math::IMultiGradFunction *gGradFunction = nullptr; + +#define PyPrint(pyo) PyObject_Print(pyo, stdout, Py_PRINT_RAW) PyObject *target_function(PyObject * /*self*/, PyObject *args) { @@ -46,10 +49,12 @@ PyObject *jac_function(PyObject * /*self*/, PyObject *args) ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() { fOptions.SetMinimizerType("Scipy"); - fOptions.SetMinimizerAlgorithm("lbfgsb"); + fOptions.SetMinimizerAlgorithm("L-BFGS-B"); if (!PyIsInitialized()) { PyInitialize(); } + // set extra options + SetAlgoExtraOptions(); } //_______________________________________________________________________ @@ -60,6 +65,19 @@ ScipyMinimizer::ScipyMinimizer(const char *type) if (!PyIsInitialized()) { PyInitialize(); } + // set extra options + SetAlgoExtraOptions(); +} + +//_______________________________________________________________________ +void ScipyMinimizer::SetAlgoExtraOptions() +{ + std::string type = fOptions.MinimizerAlgorithm(); + if (type == "L-BFGS-B") { + fExtraOpts.SetValue("gtol", 1e-10); + fExtraOpts.SetValue("eps", 1.0); + } + SetExtraOptions(fExtraOpts); } //_______________________________________________________________________ @@ -144,7 +162,19 @@ bool ScipyMinimizer::Minimize() { (gFunction) = ObjFunction(); (gGradFunction) = GradObjFunction(); + if (gGradFunction == nullptr) { + fJacobian = Py_None; + } auto method = fOptions.MinimizerAlgorithm(); + PyObject *pyoptions = PyDict_New(); + if (method == "L-BFGS-B") { + for (std::string key : fExtraOpts.GetAllRealKeys()) { + double value = 0; + fExtraOpts.GetRealValue(key.c_str(), value); + PyDict_SetItemString(pyoptions, key.c_str(), PyFloat_FromDouble(value)); + } + } + std::cout << "=== Scipy Minimization" << std::endl; std::cout << "=== Method: " << method << std::endl; std::cout << "=== Initial value: ("; @@ -157,15 +187,20 @@ bool ScipyMinimizer::Minimize() double *values = const_cast(X()); npy_intp dims[1] = {NDim()}; - PyObject *py_array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); - - PyObject *pargs = PyTuple_New(0); - - auto pyvalues = Py_BuildValue("(OOOsO)", fTarget, py_array, pargs, method.c_str(), fJacobian); - - PyObject *result = PyObject_CallObject(fMinimize, pyvalues); - Py_DECREF(pyvalues); - Py_DECREF(py_array); + PyObject *x0 = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); + + // minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, + // callback=None, options=None) + auto args = Py_BuildValue("(OO)", fTarget, x0); + auto kw = Py_BuildValue("{s:s,s:O,s:d,s:O}", "method", method.c_str(), "jac", fJacobian, "tol", Tolerance(), + "options", pyoptions); + + //PyPrint(kw); + PyObject *result = PyObject_Call(fMinimize, args, kw); + //PyPrint(result); + Py_DECREF(args); + Py_DECREF(kw); + Py_DECREF(x0); // if the minimization works PyObject *pstatus = PyObject_GetAttrString(result, "status"); From b8a1dfa60e2e487b2d707396c89fd0bfd459b3fc Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Wed, 21 Dec 2022 03:21:37 -0500 Subject: [PATCH 05/25] Scipy: pointer to values in the jacobian can not be free, It belongs to numpy --- math/scipy/src/ScipyMinimizer.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index c818521e1b4e8..489b6e3b55e5e 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -38,7 +38,7 @@ PyObject *jac_function(PyObject * /*self*/, PyObject *args) uint size = PyArray_SIZE(arr); auto params = (double *)PyArray_DATA(arr); - double values[size]; + double *values=new double[size]; gGradFunction->Gradient(params, values); npy_intp dims[1] = {size}; PyObject *py_array = PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values); @@ -74,8 +74,8 @@ void ScipyMinimizer::SetAlgoExtraOptions() { std::string type = fOptions.MinimizerAlgorithm(); if (type == "L-BFGS-B") { - fExtraOpts.SetValue("gtol", 1e-10); - fExtraOpts.SetValue("eps", 1.0); + fExtraOpts.SetValue("gtol", 1e-8); + fExtraOpts.SetValue("eps", 0.01); } SetExtraOptions(fExtraOpts); } From e6098022afa4ebc19f2e955692e36de2f66b2ab8 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Wed, 21 Dec 2022 13:29:04 -0500 Subject: [PATCH 06/25] Scipy: printing python error when minimization fails --- math/scipy/inc/Math/ScipyMinimizer.h | 5 ----- math/scipy/src/ScipyMinimizer.cxx | 9 +++++---- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index e823591578f4b..ef5061534fb2d 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -119,11 +119,6 @@ class ScipyMinimizer : public BasicMinimizer { void SetAlgoExtraOptions(); public: - /// set the function to minimize - // virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func); - - /// set the function to minimize - // virtual void SetFunction(const ROOT::Math::IMultiGradFunction &func) { BasicMinimizer::SetFunction(func); } /// method to perform the minimization virtual bool Minimize() override; diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 489b6e3b55e5e..7899f129aa98a 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -73,10 +73,6 @@ ScipyMinimizer::ScipyMinimizer(const char *type) void ScipyMinimizer::SetAlgoExtraOptions() { std::string type = fOptions.MinimizerAlgorithm(); - if (type == "L-BFGS-B") { - fExtraOpts.SetValue("gtol", 1e-8); - fExtraOpts.SetValue("eps", 0.01); - } SetExtraOptions(fExtraOpts); } @@ -197,6 +193,11 @@ bool ScipyMinimizer::Minimize() //PyPrint(kw); PyObject *result = PyObject_Call(fMinimize, args, kw); + if(result == NULL) + { + PyErr_Print(); + return false; + } //PyPrint(result); Py_DECREF(args); Py_DECREF(kw); From 07ae9b386d00baf120cc05ea1a3d5efc1dc6d600 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Thu, 22 Dec 2022 01:37:43 -0500 Subject: [PATCH 07/25] Scipy: working in the hessian implementation --- math/scipy/inc/Math/ScipyMinimizer.h | 11 +++---- math/scipy/src/ScipyMinimizer.cxx | 48 +++++++++++++++++++++++----- 2 files changed, 45 insertions(+), 14 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index ef5061534fb2d..2c6a938212170 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -26,6 +26,7 @@ #include "Rtypes.h" #include "TString.h" +#include #include #include @@ -64,7 +65,10 @@ class ScipyMinimizer : public BasicMinimizer { PyObject *fMinimize; PyObject *fTarget; PyObject *fJacobian; + PyObject *fHessian; + GenAlgoOptions fExtraOpts; + std::function &, double *)> fHessianFunc; protected: PyObject *fGlobalNS; @@ -85,11 +89,6 @@ class ScipyMinimizer : public BasicMinimizer { */ virtual ~ScipyMinimizer(); - /** - Python eval function - */ - PyObject *Eval(TString code); - /** Python initialization */ @@ -119,10 +118,10 @@ class ScipyMinimizer : public BasicMinimizer { void SetAlgoExtraOptions(); public: - /// method to perform the minimization virtual bool Minimize() override; + virtual void SetHessianFunction(std::function &, double *)>) override; template void SetExtraOption(const char *key, T value) { diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 7899f129aa98a..f7fb2192a7fb1 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -20,6 +20,8 @@ const ROOT::Math::IMultiGenFunction *gFunction = nullptr; /// function wrapper for the gradient of the function to be minimized const ROOT::Math::IMultiGradFunction *gGradFunction = nullptr; +std::function &, double *)> gfHessianFunction; + #define PyPrint(pyo) PyObject_Print(pyo, stdout, Py_PRINT_RAW) PyObject *target_function(PyObject * /*self*/, PyObject *args) @@ -38,13 +40,30 @@ PyObject *jac_function(PyObject * /*self*/, PyObject *args) uint size = PyArray_SIZE(arr); auto params = (double *)PyArray_DATA(arr); - double *values=new double[size]; + 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; }; +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); + // std::cout<<"---------------"< &, double *) -> bool { return false; }; // set extra options SetAlgoExtraOptions(); } @@ -65,6 +85,7 @@ ScipyMinimizer::ScipyMinimizer(const char *type) if (!PyIsInitialized()) { PyInitialize(); } + fHessianFunc = [](const std::vector &, double *) -> bool { return false; }; // set extra options SetAlgoExtraOptions(); } @@ -83,6 +104,7 @@ void ScipyMinimizer::PyInitialize() 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 */ }; @@ -128,9 +150,10 @@ void ScipyMinimizer::PyInitialize() // Scipy initialization PyRunString("from scipy.optimize import minimize"); fMinimize = PyDict_GetItemString(fLocalNS, "minimize"); - PyRunString("from params import target_function, jac_function"); + 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"); } //_______________________________________________________________________ @@ -158,9 +181,13 @@ bool ScipyMinimizer::Minimize() { (gFunction) = ObjFunction(); (gGradFunction) = GradObjFunction(); + gfHessianFunction = fHessianFunc; if (gGradFunction == nullptr) { fJacobian = Py_None; } + if (!gfHessianFunction) { + fHessian = Py_None; + } auto method = fOptions.MinimizerAlgorithm(); PyObject *pyoptions = PyDict_New(); if (method == "L-BFGS-B") { @@ -188,17 +215,16 @@ bool ScipyMinimizer::Minimize() // minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, // callback=None, options=None) auto args = Py_BuildValue("(OO)", fTarget, x0); - auto kw = Py_BuildValue("{s:s,s:O,s:d,s:O}", "method", method.c_str(), "jac", fJacobian, "tol", Tolerance(), - "options", pyoptions); + auto kw = Py_BuildValue("{s:s,s:O,,s:O,s:d,s:O}", "method", method.c_str(), "jac", fJacobian, "hess", fHessian, + "tol", Tolerance(), "options", pyoptions); - //PyPrint(kw); + // PyPrint(kw); PyObject *result = PyObject_Call(fMinimize, args, kw); - if(result == NULL) - { + if (result == NULL) { PyErr_Print(); return false; } - //PyPrint(result); + // PyPrint(result); Py_DECREF(args); Py_DECREF(kw); Py_DECREF(x0); @@ -243,3 +269,9 @@ void ScipyMinimizer::PyRunString(TString code, TString errorMessage, int start) exit(1); } } + +//_______________________________________________________________________ +void ScipyMinimizer::SetHessianFunction(std::function &, double *)> func) +{ + fHessianFunc = func; +} From d1a201c0ed1224a977dc53304d7db93a23cdc97c Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Wed, 12 Apr 2023 22:09:29 -0500 Subject: [PATCH 08/25] Scipy: added support for bounds --- math/scipy/inc/Math/ScipyMinimizer.h | 4 +-- math/scipy/src/ScipyMinimizer.cxx | 46 ++++++++++++++++++++++++---- 2 files changed, 42 insertions(+), 8 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index 2c6a938212170..4555e5c74cfd9 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -1,5 +1,5 @@ // @(#)root/math/scipy:$Id$ -// Author: Omar.Zapata@cern.ch 2022 +// Author: Omar.Zapata@cern.ch 2023 /************************************************************************* * Copyright (C) 1995-2022, Rene Brun and Fons Rademakers. * @@ -66,7 +66,7 @@ class ScipyMinimizer : public BasicMinimizer { PyObject *fTarget; PyObject *fJacobian; PyObject *fHessian; - + PyObject *fBoundsMod; GenAlgoOptions fExtraOpts; std::function &, double *)> fHessianFunc; diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index f7fb2192a7fb1..9362c829ffda6 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -1,14 +1,16 @@ +#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 "Math/GenAlgoOptions.h" +#include #include #include #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include +#include using namespace ROOT; using namespace ROOT::Math; @@ -148,8 +150,10 @@ void ScipyMinimizer::PyInitialize() _import_array(); // Numpy initialization } // Scipy initialization - PyRunString("from scipy.optimize import minimize"); + + 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"); @@ -162,6 +166,8 @@ void ScipyMinimizer::PyFinalize() { if (fMinimize) Py_DECREF(fMinimize); + if (fBoundsMod) + Py_DECREF(fBoundsMod); Py_Finalize(); } @@ -212,19 +218,47 @@ bool ScipyMinimizer::Minimize() 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()); + for (unsigned int i = 0; i < NDim(); i++) { + ParameterSettings varsettings; + + if (GetVariableSettings(i, varsettings)) { + if (varsettings.HasLowerLimit()) { + PyList_SetItem(pylimits_lower, i, PyFloat_FromDouble(varsettings.LowerLimit())); + } else { + PyList_SetItem(pylimits_lower, i, PyFloat_FromDouble(-NPY_INFINITY)); + } + if (varsettings.HasUpperLimit()) { + 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)); + } + } + PyTuple_SetItem(pybounds_args, 0, pylimits_lower); + PyTuple_SetItem(pybounds_args, 1, pylimits_upper); + + PyObject *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) - auto args = Py_BuildValue("(OO)", fTarget, x0); - auto kw = Py_BuildValue("{s:s,s:O,,s:O,s:d,s:O}", "method", method.c_str(), "jac", fJacobian, "hess", fHessian, - "tol", Tolerance(), "options", pyoptions); + PyObject *args = Py_BuildValue("(OO)", fTarget, x0); + PyObject *kw = Py_BuildValue("{s:s,s:O,,s:O,s:O,s:d,s:O}", "method", method.c_str(), "jac", fJacobian, "hess", + fHessian, "bounds", pybounds, "tol", Tolerance(), "options", pyoptions); - // PyPrint(kw); PyObject *result = PyObject_Call(fMinimize, args, kw); if (result == NULL) { PyErr_Print(); return false; } // PyPrint(result); + Py_DECREF(pylimits_lower); + Py_DECREF(pylimits_upper); + Py_DECREF(pybounds); Py_DECREF(args); Py_DECREF(kw); Py_DECREF(x0); From 7cd71e60e7ab7fa6be9675adf786ee8cc5908ee5 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Thu, 13 Apr 2023 00:06:37 -0500 Subject: [PATCH 09/25] Scipy: removed compilation warnings --- math/scipy/src/ScipyMinimizer.cxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 9362c829ffda6..eca3e404a1d16 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -114,7 +114,12 @@ void ScipyMinimizer::PyInitialize() "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}; + ParamsMethods, + NULL, /* m_slots */ + NULL, /* m_traverse */ + 0, /* m_clear */ + NULL /* m_free */ + }; auto PyInit_params = [](void) -> PyObject * { PyObject *m; From 8dd259c86ae209ab910beb6b0b5ac8af58663432 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Thu, 13 Apr 2023 21:52:51 -0500 Subject: [PATCH 10/25] added some documentations, removed some comments, added some deletes in the destructor --- math/scipy/inc/Math/ScipyMinimizer.h | 17 ++++++++++---- math/scipy/src/ScipyMinimizer.cxx | 34 +++++++++++++--------------- 2 files changed, 28 insertions(+), 23 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index 4555e5c74cfd9..c700cabe077d2 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -50,12 +50,19 @@ namespace Experimental { //_____________________________________________________________________________________ /** - * \class ScipyMinimizer - * ScipyMinimizer class. - * Implementation for Scipy ... TODO + \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. - See Scipy doc - from more info on the Scipy minimization algorithms. + 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 + + See Scipy doc + from more info on the Scipy minimization algorithms. @ingroup MultiMin */ diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index eca3e404a1d16..8a4c30fac3f60 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -21,11 +21,14 @@ using namespace ROOT::Fit; 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); @@ -36,6 +39,7 @@ PyObject *target_function(PyObject * /*self*/, PyObject *args) 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); @@ -49,6 +53,7 @@ PyObject *jac_function(PyObject * /*self*/, PyObject *args) 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); @@ -60,9 +65,6 @@ PyObject *hessian_function(PyObject * /*self*/, PyObject *args) gfHessianFunction(x, values); npy_intp dims[2] = {size, size}; PyObject *py_array = PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, values); - // std::cout<<"---------------"< &, double *) -> bool { return false; }; // set extra options SetAlgoExtraOptions(); @@ -84,9 +84,7 @@ ScipyMinimizer::ScipyMinimizer(const char *type) { fOptions.SetMinimizerType("Scipy"); fOptions.SetMinimizerAlgorithm(type); - if (!PyIsInitialized()) { - PyInitialize(); - } + PyInitialize(); fHessianFunc = [](const std::vector &, double *) -> bool { return false; }; // set extra options SetAlgoExtraOptions(); @@ -169,23 +167,23 @@ void ScipyMinimizer::PyInitialize() // Finalize Python interpreter void ScipyMinimizer::PyFinalize() { - if (fMinimize) - Py_DECREF(fMinimize); - if (fBoundsMod) - Py_DECREF(fBoundsMod); Py_Finalize(); } //_______________________________________________________________________ int ScipyMinimizer::PyIsInitialized() { - if (!Py_IsInitialized()) - return kFALSE; - return kTRUE; + return Py_IsInitialized(); } //_______________________________________________________________________ -ScipyMinimizer::~ScipyMinimizer() {} +ScipyMinimizer::~ScipyMinimizer() +{ + if (fMinimize) + Py_DECREF(fMinimize); + if (fBoundsMod) + Py_DECREF(fBoundsMod); +} //_______________________________________________________________________ bool ScipyMinimizer::Minimize() From 8f6ee823068d2ca722a79f938a9ca85f124ea0a9 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Thu, 13 Apr 2023 21:57:24 -0500 Subject: [PATCH 11/25] added scipy.C macro in fit tutorial --- tutorials/math/fit/scipy.C | 83 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 tutorials/math/fit/scipy.C diff --git a/tutorials/math/fit/scipy.C b/tutorials/math/fit/scipy.C new file mode 100644 index 0000000000000..e97bd1520f18d --- /dev/null +++ b/tutorials/math/fit/scipy.C @@ -0,0 +1,83 @@ +#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(const std::vector &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); + minimizer.SetExtraOption("gtol",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() < Date: Sat, 15 Apr 2023 13:38:03 -0500 Subject: [PATCH 12/25] Scipy: implemented method NCalls --- math/scipy/inc/Math/ScipyMinimizer.h | 6 +++++- math/scipy/src/ScipyMinimizer.cxx | 9 +++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index c700cabe077d2..ef3c004ce97bb 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -76,7 +76,7 @@ class ScipyMinimizer : public BasicMinimizer { PyObject *fBoundsMod; GenAlgoOptions fExtraOpts; std::function &, double *)> fHessianFunc; - + unsigned int fCalls; protected: PyObject *fGlobalNS; PyObject *fLocalNS; @@ -115,6 +115,10 @@ class ScipyMinimizer : public BasicMinimizer { */ 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; private: // usually copying is non trivial, so we make this unaccessible diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 8a4c30fac3f60..1661199192a48 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -71,6 +71,7 @@ PyObject *hessian_function(PyObject * /*self*/, PyObject *args) //_______________________________________________________________________ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() { + fCalls = 0; fOptions.SetMinimizerType("Scipy"); fOptions.SetMinimizerAlgorithm("L-BFGS-B"); PyInitialize(); @@ -82,6 +83,7 @@ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() //_______________________________________________________________________ ScipyMinimizer::ScipyMinimizer(const char *type) { + fCalls = 0; fOptions.SetMinimizerType("Scipy"); fOptions.SetMinimizerAlgorithm(type); PyInitialize(); @@ -288,6 +290,7 @@ bool ScipyMinimizer::Minimize() SetFinalValues(x); auto obj_value = (*gFunction)(x); SetMinValue(obj_value); + fCalls = nfev; //number of function evaluations std::cout << "=== Status: " << status << std::endl; std::cout << "=== Message: " << message << std::endl; @@ -312,3 +315,9 @@ void ScipyMinimizer::SetHessianFunction(std::function Date: Mon, 17 Apr 2023 21:51:26 -0500 Subject: [PATCH 13/25] Scipy: implemented maxiter and verbose option, returned the proper values in the minimizer as well. --- math/scipy/src/ScipyMinimizer.cxx | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 1661199192a48..2ffc804bffd03 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -208,7 +208,11 @@ bool ScipyMinimizer::Minimize() PyDict_SetItemString(pyoptions, key.c_str(), PyFloat_FromDouble(value)); } } - + 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: ("; @@ -270,9 +274,14 @@ bool ScipyMinimizer::Minimize() // if the minimization works PyObject *pstatus = PyObject_GetAttrString(result, "status"); - bool status = PyBool_Check(pstatus); + 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); @@ -292,10 +301,11 @@ bool ScipyMinimizer::Minimize() 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; - return status; + return success; } //_______________________________________________________________________ From 2c83094bee5e4deaa49b9deb74504366f80da9f8 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Thu, 20 Apr 2023 00:13:42 -0500 Subject: [PATCH 14/25] Scipy: working in Constraint functions, it was possible to load the function before call Py_Initilaize() --- math/scipy/inc/Math/ScipyMinimizer.h | 16 ++++- math/scipy/src/ScipyMinimizer.cxx | 102 +++++++++++++++++++++++---- 2 files changed, 102 insertions(+), 16 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index ef3c004ce97bb..c5dccba33d487 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -74,9 +74,12 @@ class ScipyMinimizer : public BasicMinimizer { 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; @@ -118,7 +121,18 @@ class ScipyMinimizer : public BasicMinimizer { /* Number of function calls */ - virtual unsigned int NCalls() const override; + 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 diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 2ffc804bffd03..67a1b9e8fe2de 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -27,7 +27,6 @@ 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) { @@ -78,6 +77,8 @@ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() fHessianFunc = [](const std::vector &, double *) -> bool { return false; }; // set extra options SetAlgoExtraOptions(); + fConstraintsList = PyList_New(0); + fConstN = 0; } //_______________________________________________________________________ @@ -90,6 +91,8 @@ ScipyMinimizer::ScipyMinimizer(const char *type) fHessianFunc = [](const std::vector &, double *) -> bool { return false; }; // set extra options SetAlgoExtraOptions(); + fConstraintsList = PyList_New(0); + fConstN = 0; } //_______________________________________________________________________ @@ -110,16 +113,18 @@ void ScipyMinimizer::PyInitialize() {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 */ - }; + 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; @@ -209,8 +214,7 @@ bool ScipyMinimizer::Minimize() } } PyDict_SetItemString(pyoptions, "maxiter", PyLong_FromLong(MaxIterations())); - if(PrintLevel()>0) - { + if (PrintLevel() > 0) { PyDict_SetItemString(pyoptions, "disp", Py_True); } std::cout << "=== Scipy Minimization" << std::endl; @@ -281,7 +285,6 @@ bool ScipyMinimizer::Minimize() 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); @@ -299,7 +302,7 @@ bool ScipyMinimizer::Minimize() SetFinalValues(x); auto obj_value = (*gFunction)(x); SetMinValue(obj_value); - fCalls = nfev; //number of function evaluations + fCalls = nfev; // number of function evaluations std::cout << "=== Success: " << success << std::endl; std::cout << "=== Status: " << status << std::endl; @@ -331,3 +334,72 @@ 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); + PyPrint(fConstraintsList); +} From 2405528a3ea6cb57f439ed913c8f60ad01083781 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Thu, 20 Apr 2023 18:01:23 -0500 Subject: [PATCH 15/25] Scipy: passed constraints list to the minimizer function --- math/scipy/src/ScipyMinimizer.cxx | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 67a1b9e8fe2de..209689023efca 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -260,8 +260,9 @@ bool ScipyMinimizer::Minimize() // 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:d,s:O}", "method", method.c_str(), "jac", fJacobian, "hess", - fHessian, "bounds", pybounds, "tol", Tolerance(), "options", pyoptions); + 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); PyObject *result = PyObject_Call(fMinimize, args, kw); if (result == NULL) { @@ -275,6 +276,8 @@ bool ScipyMinimizer::Minimize() 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"); @@ -401,5 +404,8 @@ void ScipyMinimizer::AddConstraintFunction(std::function Date: Sat, 29 Apr 2023 21:24:37 -0500 Subject: [PATCH 16/25] Scipy: removed Py_DECREF for lower and upper bounds, calling it to pybounds is enogh --- math/scipy/src/ScipyMinimizer.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 209689023efca..bb6c2e557d549 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -270,8 +270,6 @@ bool ScipyMinimizer::Minimize() return false; } // PyPrint(result); - Py_DECREF(pylimits_lower); - Py_DECREF(pylimits_upper); Py_DECREF(pybounds); Py_DECREF(args); Py_DECREF(kw); From 9ef0d94eb46fcad3f727daead15ca9b2c5936bde Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Mon, 1 May 2023 22:34:12 -0500 Subject: [PATCH 17/25] Scipy: assigned fStatus --- math/scipy/src/ScipyMinimizer.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index bb6c2e557d549..eb461e212d634 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -309,6 +309,7 @@ bool ScipyMinimizer::Minimize() std::cout << "=== Status: " << status << std::endl; std::cout << "=== Message: " << message << std::endl; std::cout << "=== Function calls: " << nfev << std::endl; + if(success) fStatus=0; return success; } From 93d4c2becae6f169bb547b476e4b30aa687cd574 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Sun, 7 May 2023 12:39:19 -0500 Subject: [PATCH 18/25] Scipy: fixed bounds, if there is not bounds It will pass Py_None added verbose output, to check parameters passed to the minimizer and output returned --- math/scipy/src/ScipyMinimizer.cxx | 41 ++++++++++++++++++++++--------- 1 file changed, 29 insertions(+), 12 deletions(-) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index eb461e212d634..d93c8f9007346 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -74,7 +74,7 @@ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() fOptions.SetMinimizerType("Scipy"); fOptions.SetMinimizerAlgorithm("L-BFGS-B"); PyInitialize(); - fHessianFunc = [](const std::vector &, double *) -> bool { return false; }; + fHessianFunc = nullptr; // set extra options SetAlgoExtraOptions(); fConstraintsList = PyList_New(0); @@ -88,7 +88,7 @@ ScipyMinimizer::ScipyMinimizer(const char *type) fOptions.SetMinimizerType("Scipy"); fOptions.SetMinimizerAlgorithm(type); PyInitialize(); - fHessianFunc = [](const std::vector &, double *) -> bool { return false; }; + fHessianFunc = nullptr; // set extra options SetAlgoExtraOptions(); fConstraintsList = PyList_New(0); @@ -98,7 +98,6 @@ ScipyMinimizer::ScipyMinimizer(const char *type) //_______________________________________________________________________ void ScipyMinimizer::SetAlgoExtraOptions() { - std::string type = fOptions.MinimizerAlgorithm(); SetExtraOptions(fExtraOpts); } @@ -201,7 +200,7 @@ bool ScipyMinimizer::Minimize() if (gGradFunction == nullptr) { fJacobian = Py_None; } - if (!gfHessianFunction) { + if (gfHessianFunction == nullptr) { fHessian = Py_None; } auto method = fOptions.MinimizerAlgorithm(); @@ -234,16 +233,19 @@ bool ScipyMinimizer::Minimize() 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)); @@ -252,24 +254,37 @@ bool ScipyMinimizer::Minimize() MATH_ERROR_MSG("ScipyMinimizer::Minimize", Form("Variable index = %d not found", i)); } } - PyTuple_SetItem(pybounds_args, 0, pylimits_lower); - PyTuple_SetItem(pybounds_args, 1, pylimits_upper); - - PyObject *pybounds = PyObject_CallObject(fBoundsMod, pybounds_args); - + 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; } - // PyPrint(result); + if(PrintLevel()>0) + { + std::cout<<"========Minimizer Results========\n"; + PyPrint(result); + std::cout<<"=================================\n"; + } Py_DECREF(pybounds); Py_DECREF(args); Py_DECREF(kw); @@ -304,12 +319,14 @@ bool ScipyMinimizer::Minimize() 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; } From 96412ebe658b94f9502f327fd391850ba7215c38 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Mon, 5 Jun 2023 00:13:27 -0500 Subject: [PATCH 19/25] Scipy: added tests --- math/scipy/CMakeLists.txt | 2 +- math/scipy/src/ScipyMinimizer.cxx | 1 - math/scipy/test/CMakeLists.txt | 11 +++ math/scipy/test/testScipyMinimizer.cxx | 130 +++++++++++++++++++++++++ 4 files changed, 142 insertions(+), 2 deletions(-) create mode 100644 math/scipy/test/CMakeLists.txt create mode 100644 math/scipy/test/testScipyMinimizer.cxx diff --git a/math/scipy/CMakeLists.txt b/math/scipy/CMakeLists.txt index 69ab2a5e88efa..0e25d4c232fe8 100644 --- a/math/scipy/CMakeLists.txt +++ b/math/scipy/CMakeLists.txt @@ -17,4 +17,4 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Scipy target_compile_definitions(Scipy PUBLIC USE_ROOT_ERROR ${PYTHON_DEFINITIONS}) target_include_directories(Scipy PUBLIC ${PYTHON_INCLUDE_DIRS}) -#ROOT_ADD_TEST_SUBDIRECTORY(test) +ROOT_ADD_TEST_SUBDIRECTORY(test) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index d93c8f9007346..4af02f9e47ad0 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -423,5 +423,4 @@ void ScipyMinimizer::AddConstraintFunction(std::function + +#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(const std::vector &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); + } +} \ No newline at end of file From 32e3d5847d19005fd74ab0df085607fb4765e446 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Mon, 12 Jun 2023 11:47:34 -0500 Subject: [PATCH 20/25] Scipy: remove method to pass extra options, this is not defined yet --- math/scipy/inc/Math/ScipyMinimizer.h | 5 ----- tutorials/math/fit/scipy.C | 1 - 2 files changed, 6 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index c5dccba33d487..4ee3ba2d2ed03 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -147,11 +147,6 @@ class ScipyMinimizer : public BasicMinimizer { virtual bool Minimize() override; virtual void SetHessianFunction(std::function &, double *)>) override; - template - void SetExtraOption(const char *key, T value) - { - fExtraOpts.SetValue(key, value); - } protected: ClassDef(ScipyMinimizer, 0) // diff --git a/tutorials/math/fit/scipy.C b/tutorials/math/fit/scipy.C index e97bd1520f18d..f99c2cac972f1 100644 --- a/tutorials/math/fit/scipy.C +++ b/tutorials/math/fit/scipy.C @@ -50,7 +50,6 @@ int scipy() minimizer.SetMaxFunctionCalls(1000000); minimizer.SetMaxIterations(100000); minimizer.SetTolerance(1e-3); - minimizer.SetExtraOption("gtol",1e-3); ROOT::Math::GradFunctor f(&RosenBrock,&RosenBrockGrad,2); double step[2] = {0.01,0.01}; double variable[2] = { -1.2,1.0}; From 8f9b91aee0f79198c7561b7cf6e9bb677d47fa69 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Mon, 12 Jun 2023 12:33:52 -0500 Subject: [PATCH 21/25] Scipy: fixed style with clang-format --- math/scipy/src/ScipyMinimizer.cxx | 37 ++++++++--------- math/scipy/test/testScipyMinimizer.cxx | 13 +++--- tutorials/math/fit/scipy.C | 57 +++++++++++++------------- 3 files changed, 52 insertions(+), 55 deletions(-) diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index 4af02f9e47ad0..abcd706be0c3b 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -74,7 +74,7 @@ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() fOptions.SetMinimizerType("Scipy"); fOptions.SetMinimizerAlgorithm("L-BFGS-B"); PyInitialize(); - fHessianFunc = nullptr; + fHessianFunc = nullptr; // set extra options SetAlgoExtraOptions(); fConstraintsList = PyList_New(0); @@ -239,13 +239,13 @@ bool ScipyMinimizer::Minimize() if (GetVariableSettings(i, varsettings)) { if (varsettings.HasLowerLimit()) { - foundBounds=true; + 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; + foundBounds = true; PyList_SetItem(pylimits_upper, i, PyFloat_FromDouble(varsettings.UpperLimit())); } else { PyList_SetItem(pylimits_upper, i, PyFloat_FromDouble(NPY_INFINITY)); @@ -255,35 +255,32 @@ bool ScipyMinimizer::Minimize() } } PyObject *pybounds = Py_None; - if(foundBounds) - { + 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"; + 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"; + 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"; + if (PrintLevel() > 0) { + std::cout << "========Minimizer Results========\n"; PyPrint(result); - std::cout<<"=================================\n"; + std::cout << "=================================\n"; } Py_DECREF(pybounds); Py_DECREF(args); @@ -323,8 +320,10 @@ bool ScipyMinimizer::Minimize() 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. + if (success) + fStatus = 0; + else + fStatus = status; // suggested by Lorenzo. Py_DECREF(result); return success; diff --git a/math/scipy/test/testScipyMinimizer.cxx b/math/scipy/test/testScipyMinimizer.cxx index 50132bf2a066b..4cb5d9414801c 100644 --- a/math/scipy/test/testScipyMinimizer.cxx +++ b/math/scipy/test/testScipyMinimizer.cxx @@ -48,9 +48,9 @@ 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"}; +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 { @@ -78,15 +78,15 @@ class ScipyFitClass : public ::testing::Test { minimizer->SetVariable(1, "y", variable[1], step[1]); if (useConstraint) minimizer->AddConstraintFunction(ConstRosenBrock, "ineq"); - + minimizer->Minimize(); } }; -TEST_F(ScipyFitClass, Fit) +TEST_F(ScipyFitClass, Fit) { for (const std::string &text : methods) { - Fit(text.c_str(),false); + Fit(text.c_str(), false); EXPECT_EQ(1000000, minimizer->MaxFunctionCalls()); EXPECT_EQ(100000, minimizer->MaxIterations()); @@ -117,7 +117,6 @@ TEST_F(ScipyFitClass, Fit) } } - TEST_F(ScipyFitClass, FitContraint) // using constraint function { for (const std::string &text : methods) { diff --git a/tutorials/math/fit/scipy.C b/tutorials/math/fit/scipy.C index f99c2cac972f1..2ba25cb3413d0 100644 --- a/tutorials/math/fit/scipy.C +++ b/tutorials/math/fit/scipy.C @@ -6,14 +6,13 @@ #include "Math/MinimizerOptions.h" #include "TStopwatch.h" - -double RosenBrock(const double *xx ) +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; + 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) @@ -28,49 +27,49 @@ bool RosenBrockHessian(const std::vector &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[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"}; +{ + + 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) - { + 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}; - + 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.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(); + 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() < Date: Mon, 3 Jul 2023 18:09:38 -0500 Subject: [PATCH 22/25] Scipy: Added support for extra options using GenAlgoOptions --- math/scipy/inc/Math/ScipyMinimizer.h | 10 ++++---- math/scipy/src/ScipyMinimizer.cxx | 36 ++++++++++++++++++---------- 2 files changed, 29 insertions(+), 17 deletions(-) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index 4ee3ba2d2ed03..37b7c12b449ef 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -21,8 +21,6 @@ #include "Math/BasicMinimizer.h" -#include "Math/GenAlgoOptions.h" - #include "Rtypes.h" #include "TString.h" @@ -75,7 +73,7 @@ class ScipyMinimizer : public BasicMinimizer { PyObject *fHessian; PyObject *fBoundsMod; PyObject *fConstraintsList; /// contraints functions - GenAlgoOptions fExtraOpts; + GenAlgoOptions *fExtraOpts; std::function &, double *)> fHessianFunc; unsigned int fConstN; unsigned int fCalls; @@ -140,7 +138,11 @@ class ScipyMinimizer : public BasicMinimizer { Copy constructor */ ScipyMinimizer(const ScipyMinimizer &) : BasicMinimizer() {} - void SetAlgoExtraOptions(); + + /** + Get extra options from IOptions + */ + void SetExtraOptions(); public: /// method to perform the minimization diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index abcd706be0c3b..bd3eea3e5882c 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -75,10 +75,9 @@ ScipyMinimizer::ScipyMinimizer() : BasicMinimizer() fOptions.SetMinimizerAlgorithm("L-BFGS-B"); PyInitialize(); fHessianFunc = nullptr; - // set extra options - SetAlgoExtraOptions(); fConstraintsList = PyList_New(0); fConstN = 0; + fExtraOpts = nullptr; } //_______________________________________________________________________ @@ -89,16 +88,9 @@ ScipyMinimizer::ScipyMinimizer(const char *type) fOptions.SetMinimizerAlgorithm(type); PyInitialize(); fHessianFunc = nullptr; - // set extra options - SetAlgoExtraOptions(); fConstraintsList = PyList_New(0); fConstN = 0; -} - -//_______________________________________________________________________ -void ScipyMinimizer::SetAlgoExtraOptions() -{ - SetExtraOptions(fExtraOpts); + fExtraOpts = nullptr; } //_______________________________________________________________________ @@ -191,6 +183,13 @@ ScipyMinimizer::~ScipyMinimizer() Py_DECREF(fBoundsMod); } +//_______________________________________________________________________ +void ScipyMinimizer::SetExtraOptions() +{ + auto constExtraOpts = dynamic_cast(fOptions.ExtraOptions()); + fExtraOpts = const_cast(constExtraOpts); +} + //_______________________________________________________________________ bool ScipyMinimizer::Minimize() { @@ -205,12 +204,23 @@ bool ScipyMinimizer::Minimize() } auto method = fOptions.MinimizerAlgorithm(); PyObject *pyoptions = PyDict_New(); - if (method == "L-BFGS-B") { - for (std::string key : fExtraOpts.GetAllRealKeys()) { + SetExtraOptions(); + if (fExtraOpts) { + for (std::string key : fExtraOpts->GetAllRealKeys()) { double value = 0; - fExtraOpts.GetRealValue(key.c_str(), value); + 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) { From 706ae0e431bc2171905567f22ed4ef1286730b52 Mon Sep 17 00:00:00 2001 From: Omar Zapata Date: Tue, 4 Jul 2023 01:28:44 -0500 Subject: [PATCH 23/25] Scipy: added documentation for extra options --- math/scipy/inc/Math/ScipyMinimizer.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index 37b7c12b449ef..e01c3a0dac817 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -59,6 +59,16 @@ namespace Experimental { 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. From 2cee20e877a1ed2877790ed2db33d9fe6d04f03a Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 29 Jul 2025 13:24:08 +0200 Subject: [PATCH 24/25] [math] Make ScipyMinimizer compile again --- math/scipy/CMakeLists.txt | 4 +++- math/scipy/inc/Math/ScipyMinimizer.h | 4 ++-- math/scipy/src/ScipyMinimizer.cxx | 3 ++- math/scipy/test/testScipyMinimizer.cxx | 4 ++-- tutorials/math/fit/scipy.C | 4 ++-- 5 files changed, 11 insertions(+), 8 deletions(-) diff --git a/math/scipy/CMakeLists.txt b/math/scipy/CMakeLists.txt index 0e25d4c232fe8..76a57ccffb7e8 100644 --- a/math/scipy/CMakeLists.txt +++ b/math/scipy/CMakeLists.txt @@ -9,7 +9,9 @@ ROOT_STANDARD_LIBRARY_PACKAGE(Scipy Math/ScipyMinimizer.h LINKDEF Math/LinkDef.h - LIBRARIES ${PYTHON_LIBRARIES_Development_Main} + LIBRARIES + Python3::NumPy + Python3::Python SOURCES src/ScipyMinimizer.cxx DEPENDENCIES Core MathCore RIO diff --git a/math/scipy/inc/Math/ScipyMinimizer.h b/math/scipy/inc/Math/ScipyMinimizer.h index e01c3a0dac817..00c85cf01b806 100644 --- a/math/scipy/inc/Math/ScipyMinimizer.h +++ b/math/scipy/inc/Math/ScipyMinimizer.h @@ -84,7 +84,7 @@ class ScipyMinimizer : public BasicMinimizer { PyObject *fBoundsMod; PyObject *fConstraintsList; /// contraints functions GenAlgoOptions *fExtraOpts; - std::function &, double *)> fHessianFunc; + std::function, double *)> fHessianFunc; unsigned int fConstN; unsigned int fCalls; @@ -158,7 +158,7 @@ class ScipyMinimizer : public BasicMinimizer { /// method to perform the minimization virtual bool Minimize() override; - virtual void SetHessianFunction(std::function &, double *)>) override; + virtual void SetHessianFunction(std::function, double *)>) override; protected: ClassDef(ScipyMinimizer, 0) // diff --git a/math/scipy/src/ScipyMinimizer.cxx b/math/scipy/src/ScipyMinimizer.cxx index bd3eea3e5882c..8e130dfa04120 100644 --- a/math/scipy/src/ScipyMinimizer.cxx +++ b/math/scipy/src/ScipyMinimizer.cxx @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -352,7 +353,7 @@ void ScipyMinimizer::PyRunString(TString code, TString errorMessage, int start) } //_______________________________________________________________________ -void ScipyMinimizer::SetHessianFunction(std::function &, double *)> func) +void ScipyMinimizer::SetHessianFunction(std::function, double *)> func) { fHessianFunc = func; } diff --git a/math/scipy/test/testScipyMinimizer.cxx b/math/scipy/test/testScipyMinimizer.cxx index 4cb5d9414801c..153f8d49f9185 100644 --- a/math/scipy/test/testScipyMinimizer.cxx +++ b/math/scipy/test/testScipyMinimizer.cxx @@ -28,7 +28,7 @@ double RosenBrockGrad(const double *x, unsigned int ipar) return 200 * (x[1] - x[0] * x[0]); } // Hessian function -bool RosenBrockHessian(const std::vector &xx, double *hess) +bool RosenBrockHessian(std::span xx, double *hess) { const double x = xx[0]; const double y = xx[1]; @@ -126,4 +126,4 @@ TEST_F(ScipyFitClass, FitContraint) // using constraint function ASSERT_NEAR(1, x[1], 0.5); ASSERT_NEAR(0, RosenBrock(x), 0.5); } -} \ No newline at end of file +} diff --git a/tutorials/math/fit/scipy.C b/tutorials/math/fit/scipy.C index 2ba25cb3413d0..d1cbc2a543cb6 100644 --- a/tutorials/math/fit/scipy.C +++ b/tutorials/math/fit/scipy.C @@ -23,7 +23,7 @@ double RosenBrockGrad(const double *x, unsigned int ipar) return 200 * (x[1] - x[0] * x[0]); } -bool RosenBrockHessian(const std::vector &xx, double *hess) +bool RosenBrockHessian(std::span xx, double *hess) { const double x = xx[0]; const double y = xx[1]; @@ -78,4 +78,4 @@ int scipy() int main() { return scipy(); -} \ No newline at end of file +} From 2e498fe6f0c412e179508ac455ef749b2784a41b Mon Sep 17 00:00:00 2001 From: Jonas Rembser Date: Tue, 29 Jul 2025 10:42:06 +0200 Subject: [PATCH 25/25] [ci] Build with `scipy=ON` --- cmake/modules/SearchInstalledSoftware.cmake | 2 -- math/CMakeLists.txt | 4 +--- 2 files changed, 1 insertion(+), 5 deletions(-) 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/math/CMakeLists.txt b/math/CMakeLists.txt index 2a72ef78b012c..becbf16b1e7cd 100644 --- a/math/CMakeLists.txt +++ b/math/CMakeLists.txt @@ -34,8 +34,6 @@ if(r) add_subdirectory(rtools) endif() -if(scipy AND NUMPY_FOUND) - add_subdirectory(scipy) -endif() +add_subdirectory(scipy) add_subdirectory(vecops)