Skip to content

Commit b525dc1

Browse files
omazapaomazapa
committed
Math: working in the initial prototype for Scipy optimizers
1 parent 6fc4ddc commit b525dc1

File tree

6 files changed

+379
-0
lines changed

6 files changed

+379
-0
lines changed

‎etc/plugins/ROOT@@Math@@Minimizer/P100_Scipy.C‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
void P100_Scipy()
2+
{
3+
gPluginMgr->AddHandler("ROOT::Math::Minimizer", "Scipy", "ROOT::Math::Experimental::ScipyMinimizer", "Scipy",
4+
"ScipyMinimizer(const char*)");
5+
}

‎math/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,4 +38,8 @@ if(r)
3838
add_subdirectory(rtools)
3939
endif()
4040

41+
if(scipy AND NUMPY_FOUND)
42+
add_subdirectory(scipy)
43+
endif()
44+
4145
add_subdirectory(vecops)

‎math/scipy/CMakeLists.txt‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
############################################################################
2+
# CMakeLists.txt file for building ROOT math/scipy package
3+
############################################################################
4+
# author: [email protected] 2022
5+
6+
include_directories(SYSTEM ${PYTHON_INCLUDE_DIRS_Development_Main} ${NUMPY_INCLUDE_DIRS})
7+
ROOT_STANDARD_LIBRARY_PACKAGE(Scipy
8+
HEADERS
9+
Math/ScipyMinimizer.h
10+
LINKDEF
11+
Math/LinkDef.h
12+
LIBRARIES ${PYTHON_LIBRARIES_Development_Main}
13+
SOURCES
14+
src/ScipyMinimizer.cxx
15+
DEPENDENCIES Core MathCore RIO
16+
)
17+
18+
target_compile_definitions(Scipy PUBLIC USE_ROOT_ERROR ${PYTHON_DEFINITIONS})
19+
target_include_directories(Scipy PUBLIC ${PYTHON_INCLUDE_DIRS})
20+
#ROOT_ADD_TEST_SUBDIRECTORY(test)

‎math/scipy/inc/Math/LinkDef.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
// @(#)root/math/ipopt:$Id$
2+
// Authors: [email protected] 12/2022
3+
4+
#ifdef __CINT__
5+
6+
#pragma link C++ nestedclasses;
7+
#pragma link C++ nestedtypedef;
8+
9+
#pragma link C++ class ROOT::Math::Experimental::ScipyMinimizer;
10+
11+
#endif //__CINT__

‎math/scipy/inc/Math/ScipyMinimizer.h‎

Lines changed: 145 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,145 @@
1+
// @(#)root/math/scipy:$Id$
2+
// Author: [email protected] 2022
3+
4+
/*************************************************************************
5+
* Copyright (C) 1995-2022, Rene Brun and Fons Rademakers. *
6+
* All rights reserved. *
7+
* *
8+
* For the licensing terms see $ROOTSYS/LICENSE. *
9+
* For the list of contributors see $ROOTSYS/README/CREDITS. *
10+
*************************************************************************/
11+
12+
#ifndef ROOT_Math_ScipyMinimizer
13+
#define ROOT_Math_ScipyMinimizer
14+
15+
#include "Math/Minimizer.h"
16+
#include "Math/MinimizerOptions.h"
17+
18+
#include "Math/IFunctionfwd.h"
19+
20+
#include "Math/IParamFunctionfwd.h"
21+
22+
#include "Math/BasicMinimizer.h"
23+
24+
#include "Rtypes.h"
25+
#include "TString.h"
26+
27+
#include <vector>
28+
#include <map>
29+
30+
#ifndef PyObject_HEAD
31+
struct _object;
32+
typedef _object PyObject;
33+
#define Py_single_input 256
34+
#endif
35+
36+
namespace ROOT {
37+
38+
namespace Math {
39+
40+
namespace Experimental {
41+
/**
42+
enumeration specifying the types of Scipy solvers
43+
@ingroup MultiMin
44+
*/
45+
46+
47+
//_____________________________________________________________________________________
48+
/**
49+
* \class ScipyMinimizer
50+
* ScipyMinimizer class.
51+
* Implementation for Scipy ... TODO
52+
53+
See <A HREF="https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html">Scipy doc</A>
54+
from more info on the Scipy minimization algorithms.
55+
56+
@ingroup MultiMin
57+
*/
58+
59+
class ScipyMinimizer : public BasicMinimizer {
60+
private:
61+
PyObject *fMinimize;
62+
PyObject *fTarget;
63+
protected:
64+
PyObject *fGlobalNS;
65+
PyObject *fLocalNS;
66+
public:
67+
/**
68+
Default constructor
69+
*/
70+
ScipyMinimizer();
71+
/**
72+
Constructor with a string giving name of algorithm
73+
*/
74+
ScipyMinimizer(const char *type);
75+
76+
/**
77+
Destructor
78+
*/
79+
virtual ~ScipyMinimizer();
80+
81+
/**
82+
Python eval function
83+
*/
84+
PyObject *Eval(TString code);
85+
86+
/**
87+
Python initialization
88+
*/
89+
void PyInitialize();
90+
91+
/**
92+
Checks if Python was initialized
93+
*/
94+
int PyIsInitialized();
95+
96+
/*
97+
Python finalization
98+
*/
99+
void PyFinalize();
100+
void PyRunString(TString code, TString errorMessage="Failed to run python code", int start=Py_single_input);
101+
void LoadWrappers();
102+
private:
103+
// usually copying is non trivial, so we make this unaccessible
104+
105+
/**
106+
Copy constructor
107+
*/
108+
ScipyMinimizer(const ScipyMinimizer &) : BasicMinimizer() {}
109+
110+
111+
public:
112+
/// set the function to minimize
113+
//virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func);
114+
115+
/// set the function to minimize
116+
//virtual void SetFunction(const ROOT::Math::IMultiGradFunction &func) { BasicMinimizer::SetFunction(func); }
117+
118+
/// method to perform the minimization
119+
virtual bool Minimize();
120+
121+
/// return expected distance reached from the minimum
122+
virtual double Edm() const { return 0; } // not impl. }
123+
124+
/// minimizer provides error and error matrix
125+
virtual bool ProvidesError() const { return false; }
126+
127+
/// return errors at the minimum
128+
virtual const double *Errors() const { return 0; }
129+
130+
/** return covariance matrices elements
131+
if the variable is fixed the matrix is zero
132+
The ordering of the variables is the same as in errors
133+
*/
134+
virtual double CovMatrix(unsigned int, unsigned int) const { return 0; }
135+
protected:
136+
ClassDef(ScipyMinimizer, 0) //
137+
};
138+
139+
} // end namespace Experimental
140+
141+
} // end namespace Math
142+
143+
} // end namespace ROOT
144+
145+
#endif /* ROOT_Math_ScipyMinimizer */

‎math/scipy/src/ScipyMinimizer.cxx‎

Lines changed: 194 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,194 @@
1+
#include <Python.h> // Needs to be included first to avoid redefinition of _POSIX_C_SOURCE
2+
#include <Math/ScipyMinimizer.h>
3+
#include <Fit/ParameterSettings.h>
4+
#include <Math/IFunction.h>
5+
#include <Math/FitMethodFunction.h>
6+
#include <TString.h>
7+
#include <iostream>
8+
9+
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
10+
#include <numpy/arrayobject.h>
11+
12+
using namespace ROOT;
13+
using namespace ROOT::Math;
14+
using namespace ROOT::Math::Experimental;
15+
using namespace ROOT::Fit;
16+
17+
18+
/// function wrapper for the function to be minimized
19+
const ROOT::Math::IMultiGenFunction *gFunction;
20+
/// function wrapper for the gradient of the function to be minimized
21+
const ROOT::Math::IMultiGradFunction *gGradFunction;
22+
23+
PyObject * target_function(PyObject */*self*/, PyObject *args)
24+
{
25+
PyArrayObject *arr = (PyArrayObject *)PyTuple_GetItem(args, 0);
26+
27+
auto params = (double *)PyArray_DATA(arr);
28+
auto r = (*gFunction)(params);
29+
30+
return PyFloat_FromDouble(r);
31+
};
32+
33+
34+
//_______________________________________________________________________
35+
ScipyMinimizer::ScipyMinimizer() : BasicMinimizer()
36+
{
37+
fOptions.SetMinimizerType("Scipy");
38+
fOptions.SetMinimizerAlgorithm("lbfgsb");
39+
if (!PyIsInitialized()) {
40+
PyInitialize();
41+
}
42+
43+
}
44+
45+
46+
47+
//_______________________________________________________________________
48+
ScipyMinimizer::ScipyMinimizer(const char *type)
49+
{
50+
fOptions.SetMinimizerType("Scipy");
51+
fOptions.SetMinimizerAlgorithm(type);
52+
if (!PyIsInitialized()) {
53+
PyInitialize();
54+
}
55+
}
56+
57+
//_______________________________________________________________________
58+
void ScipyMinimizer::PyInitialize()
59+
{
60+
static PyObject *ParamsError;
61+
static PyMethodDef ParamsMethods[] = {
62+
{"target_function", target_function, METH_VARARGS,
63+
"Target function to minimize."},
64+
{NULL, NULL, 0, NULL} /* Sentinel */
65+
};
66+
67+
static struct PyModuleDef paramsmodule = {
68+
PyModuleDef_HEAD_INIT,
69+
"params", /* name of module */
70+
"ROOT Scipy parameters", /* module documentation, may be NULL */
71+
-1, /* size of per-interpreter state of the module,
72+
or -1 if the module keeps state in global variables. */
73+
ParamsMethods
74+
};
75+
76+
auto PyInit_params = [](void)->PyObject *
77+
{
78+
PyObject *m;
79+
80+
m = PyModule_Create(&paramsmodule);
81+
if (m == NULL)
82+
return NULL;
83+
84+
ParamsError = PyErr_NewException("params.error", NULL, NULL);
85+
Py_XINCREF(ParamsError);
86+
if (PyModule_AddObject(m, "error", ParamsError) < 0) {
87+
Py_XDECREF(ParamsError);
88+
Py_CLEAR(ParamsError);
89+
Py_DECREF(m);
90+
return NULL;
91+
}
92+
93+
return m;
94+
};
95+
if (PyImport_AppendInittab("params", PyInit_params) == -1) {
96+
MATH_ERROR_MSG("ScipyMinimizer::Minimize", "could not extend in-built modules table");
97+
exit(1);
98+
}
99+
100+
bool pyIsInitialized = PyIsInitialized();
101+
if (!pyIsInitialized) {
102+
Py_Initialize(); // Python initialization
103+
}
104+
fLocalNS = PyDict_New();
105+
fGlobalNS = PyDict_New();
106+
107+
if (!pyIsInitialized) {
108+
_import_array(); // Numpy initialization
109+
}
110+
//Scipy initialization
111+
PyRunString("from scipy.optimize import minimize");
112+
fMinimize = PyDict_GetItemString(fLocalNS, "minimize");
113+
PyRunString("from params import target_function");
114+
fTarget = PyDict_GetItemString(fLocalNS, "target_function");
115+
}
116+
117+
//_______________________________________________________________________
118+
// Finalize Python interpreter
119+
void ScipyMinimizer::PyFinalize()
120+
{
121+
if (fMinimize) Py_DECREF(fMinimize);
122+
Py_Finalize();
123+
}
124+
125+
//_______________________________________________________________________
126+
int ScipyMinimizer::PyIsInitialized()
127+
{
128+
if (!Py_IsInitialized()) return kFALSE;
129+
return kTRUE;
130+
}
131+
132+
//_______________________________________________________________________
133+
ScipyMinimizer::~ScipyMinimizer()
134+
{
135+
}
136+
137+
//_______________________________________________________________________
138+
bool ScipyMinimizer::Minimize()
139+
{
140+
(gFunction)= ObjFunction();
141+
(gGradFunction) = GradObjFunction();
142+
auto method = fOptions.MinimizerAlgorithm();
143+
144+
std::cout<<"=== Scipy Minimization"<<std::endl;
145+
std::cout<<"=== Method: "<<method<<std::endl;
146+
std::cout<<"=== Initial value: (";
147+
for(uint i=0;i<NDim();i++ )
148+
{
149+
std::cout<<X()[i];
150+
if(i<NDim()-1)
151+
std::cout<<",";
152+
}
153+
std::cout<<")"<<std::endl;
154+
155+
double *values = const_cast<double*>(X());
156+
npy_intp dims[1] = { NDim()};
157+
PyObject *py_array= PyArray_SimpleNewFromData(1, dims, NPY_DOUBLE, values);
158+
159+
PyObject *pargs=PyTuple_New(0);
160+
161+
auto pyvalues = Py_BuildValue("(OOOs)",fTarget,py_array,pargs,method.c_str());
162+
163+
PyObject *result = PyObject_CallObject(fMinimize,pyvalues);
164+
Py_DECREF(pyvalues);
165+
Py_DECREF(py_array);
166+
167+
168+
//if the minimization works
169+
PyObject *pstatus = PyObject_GetAttrString(result,"status");
170+
bool status = PyBool_Check(pstatus );
171+
Py_DECREF(pstatus);
172+
173+
//the x values for the minimum
174+
PyArrayObject *pyx = (PyArrayObject *)PyObject_GetAttrString(result,"x");
175+
const double *x =(const double *)PyArray_DATA(pyx);
176+
Py_DECREF(pyx);
177+
178+
SetFinalValues(x);
179+
auto obj_value = (*gFunction)(x);
180+
SetMinValue(obj_value);
181+
182+
return status;
183+
}
184+
185+
//_______________________________________________________________________
186+
void ScipyMinimizer::PyRunString(TString code, TString errorMessage, int start) {
187+
auto fPyReturn = PyRun_String(code, start, fGlobalNS, fLocalNS);
188+
if (!fPyReturn) {
189+
auto msg = errorMessage + Form(" \ncode = \"%s\"",code.Data());
190+
MATH_ERROR_MSG("ScipyMinimizer::PyRunString", msg.Data() );
191+
PyErr_Print();
192+
exit(1);
193+
}
194+
}

0 commit comments

Comments
 (0)