-
Notifications
You must be signed in to change notification settings - Fork 0
/
func.hpp
68 lines (56 loc) · 1.97 KB
/
func.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
namespace callpy {
#define PY_SSIZE_T_CLEAN
#include "Python.h"
#include "object.h"
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/arrayobject.h>
#include <math.h>
void foo(int n, const double *x, double *c_lp, double *c_glp) {
PyGILState_STATE gstate;
gstate = PyGILState_Ensure();
_import_array();
PyObject *main = PyImport_ImportModule("__main__");
npy_intp shape[1] = {n};
PyObject *np_x = PyArray_SimpleNewFromData(1, shape, NPY_DOUBLE, (void*) x);
if (c_lp != NULL)
{
PyObject *lp = PyObject_GetAttrString(main, "lp");
PyArrayObject *lp_val = (PyArrayObject*) PyObject_CallFunctionObjArgs(lp, np_x, NULL);
PyArray_ScalarAsCtype((PyObject*) lp_val, (void*) &c_lp);
Py_DECREF(lp);
Py_DECREF(lp_val);
}
if (c_glp != NULL)
{
PyObject *glp = PyObject_GetAttrString(main, "glp");
PyArrayObject *glp_val = (PyArrayObject*) PyObject_CallFunctionObjArgs(glp, np_x, NULL);
double *c_glp_val = (double*) PyArray_DATA(glp_val);
for (int i=0; i<n; i )
c_glp[i] = c_glp_val[i];
Py_DECREF(glp);
Py_DECREF(glp_val);
}
PyGILState_Release(gstate);
}
}
template <typename T0__>
typename boost::math::tools::promote_args<T0__>::type
userfunc2(const Eigen::Matrix<T0__, Eigen::Dynamic, 1>& theta, std::ostream* pstream__) {
throw std::logic_error("not implemented");
}
template <> double userfunc2(const Eigen::Matrix<double, -1, 1>& theta, std::ostream* pstream__) {
double out;
callpy::foo(theta.rows(), theta.data(), &out, NULL);
return out;
}
template <> var userfunc2(const Eigen::Matrix<var, -1, 1>& theta, std::ostream* pstream__) {
const Eigen::Matrix<double, -1, 1> theta_val = value_of(theta);
double fa;
double *ga = new double[theta.rows()];
callpy::foo(theta.rows(), theta_val.data(), &fa, ga);
std::vector<double> std_ga(theta.rows());
std_ga.assign(ga, ga theta.rows());
delete ga;
vector<var> theta_std(theta.data(), theta.data() theta.rows());
return precomputed_gradients(fa, theta_std, std_ga);
}