Robert Kern wrote:
> On Thu, Mar 14, 2013 at 11:00 AM, Neal Becker <[email protected]> wrote:
>> Robert Kern wrote:
>>
>>> On Wed, Mar 13, 2013 at 12:16 AM, Neal Becker <[email protected]> wrote:
>>>> I guess I talked to you about 100 years ago about sharing state between
>>>> numpy
>>>> rng and code I have in c++ that wraps boost::random. So is there a C-api
>>>> for
>>>> this RandomState object I could use to call from c++? Maybe I could do
>>>> something with that.
>>>
>>> There is not one currently. Cython has provisions for sharing such
>>> low-level access to other Cython extensions, but I'm not sure how well
>>> it works for exporting data pointers and function pointers to general
>>> C/++ code. We could probably package the necessities into a struct and
>>> export a pointer to it via a PyCapsule.
>>>
>>
>> One thing this code doesn't do: it requires construction of the wrapper class
>> passing in a RandomState object. It doesn't verify you actually gave it a
>> RandomState object. It's hard to do that. The problem as I see it is to
>> perform this check, I need the RandomStateType object, which unfortunately
>> mtrand.so does not export.
>>
>> The only way to do it is in c++ code:
>>
>> 1. import numpy.random
>> 2. get RandomState class
>> 3. call it to create RandomState instance
>> 4. get the ob_type pointer.
>>
>> Pretty ugly:
>>
>> object mod = object (handle<>
>> (borrowed((PyImport_ImportModule("numpy.random")))));
>> object rs_obj = mod.attr("RandomState");
>> object rs_inst = call<object> (rs_obj.ptr(), 0);
>> RandomStateTypeObj = rs_inst.ptr()->ob_type;
>
> PyObject_IsInstance() should be sufficient.
>
> http://docs.python.org/2/c-api/object.html#PyObject_IsInstance
>
> --
> Robert Kern
Thanks!
For the record, an updated version attached.
#include <boost/python.hpp>
// Include this AFTER Python.h; in projects with multiple source files, you may have to
// #define some things first (see the numpy C-API documentation).
#include <numpy/arrayobject.h>
#include <boost/numpy/dtype.hpp>
#include <boost/python/class.hpp>
#include <boost/python/init.hpp>
#include <boost/python/module.hpp>
#include <boost/python/def.hpp>
#include <boost/python/args.hpp>
#include <boost/python/str.hpp>
#include <boost/python/call.hpp>
#include <boost/numpy.hpp>
#include <ndarray/Array.h>
#include <ndarray/bp/auto.h>
#include <inttypes.h>
#include <stdexcept>
using namespace boost::python;
namespace bp=boost::python;
namespace nd=ndarray;
#define RK_STATE_LEN 624
typedef struct rk_state_
{
unsigned long key[RK_STATE_LEN];
int pos;
int has_gauss; /* !=0: gauss contains a gaussian deviate */
double gauss;
/* The rk_state structure has been extended to store the following
* information for the binomial generator. If the input values of n or p
* are different than nsave and psave, then the other parameters will be
* recomputed. RTK 2005-09-02 */
int has_binomial; /* !=0: following parameters initialized for
binomial */
double psave;
long nsave;
double r;
double q;
double fm;
long m;
double p1;
double xm;
double xl;
double xr;
double c;
double laml;
double lamr;
double p2;
double p3;
double p4;
}
rk_state;
struct RandomState {
PyObject_HEAD
rk_state *internal_state;
};
extern "C" void rk_fill(void *buffer, size_t size, rk_state *state);
extern "C" unsigned long rk_ulong(rk_state *state);
extern "C" unsigned long rk_random(rk_state *state);
inline NPY_TYPES dtype_to_type (bp::object const& dtype, NPY_TYPES _default) {
if (dtype == bp::object())
return _default;
else {
boost::numpy::dtype dt (dtype);
return NPY_TYPES (reinterpret_cast<PyArray_Descr*>(dt.ptr())->type_num);
// PyArray_Descr *tp_descr;
// if (PyArray_DescrConverter(dtype.ptr(), &tp_descr) != NPY_SUCCEED)
// throw bp::error_already_set();
// return NPY_TYPES(tp_descr->type_num);
}
}
template<int> struct type_imap;
template<>
struct type_imap<NPY_INT8> {
typedef npy_int8 type;
};
template<>
struct type_imap<NPY_INT16> {
typedef npy_int16 type;
};
template<>
struct type_imap<NPY_INT32> {
typedef npy_int32 type;
};
template<>
struct type_imap<NPY_INT64> {
typedef npy_int64 type;
};
template<>
struct type_imap<NPY_FLOAT32> {
typedef float type;
};
template<>
struct type_imap<NPY_FLOAT64> {
typedef double type;
};
template<>
struct type_imap<NPY_COMPLEX64> {
typedef std::complex<float> type;
};
template<>
struct type_imap<NPY_COMPLEX128> {
typedef std::complex<double> type;
};
template <class IntType=uint32_t>
class pnseq_generator {
public:
typedef IntType result_type;
typedef IntType input_type;
explicit pnseq_generator(int _width) : width (_width), count (0),
mask (
_width == std::numeric_limits<result_type>::digits ? ~(result_type(0)) : ~(result_type(-1LL) << _width))
{
if (std::numeric_limits<result_type>::digits < width or width <= 0)
throw std::runtime_error ("pnseq_generator::invalid width");
}
explicit pnseq_generator (int _width, int _count, result_type _cache):
width (_width), count (_count), cache (_cache), mask (~(result_type(-1LL) << _width)) {}
template<class Engine>
result_type operator() (Engine& eng) {
// BOOST_STATIC_ASSERT (sizeof (typename Engine::result_type) >= width);
if (count <= width-1) {
cache = eng();
count = std::min (std::numeric_limits<typename Engine::result_type>::digits, std::numeric_limits<result_type>::digits);
}
result_type bits = cache & mask;
cache >>= width;
count -= width;
return bits;
}
// private:
const int width;
int count;
result_type cache;
const result_type mask;
};
static object rs_obj;
template<typename out_type=uint32_t>
struct engine {
typedef out_type result_type;
bp::object rs; // RandomState object
void verify (object const& o) {
if (!PyObject_IsInstance(o.ptr(), rs_obj.ptr()))
throw std::runtime_error ("engine: expected RandomState");
}
engine (bp::object const& _rs) : rs (_rs) { verify (_rs); }
result_type operator() () {
RandomState* r = (RandomState*)(rs.ptr());
// result_type buffer;
// rk_fill ((void*)&buffer, sizeof(buffer), r->internal_state);
if (sizeof(result_type) == sizeof (uint64_t))
return rk_ulong (r->internal_state);
else if (sizeof(result_type) == sizeof (uint32_t))
return rk_random (r->internal_state);
}
};
template<typename result_type=uint32_t>
struct wrapped_pn {
engine<result_type> eng;
pnseq_generator<result_type> pn;
wrapped_pn (bp::object const& _rs, int _width) :
eng (_rs), pn (_width) {}
// for pickle
wrapped_pn (bp::object const& _rs, int width, int count, typename pnseq_generator<result_type>::result_type cache) :
eng (_rs),
pn (width, count, cache) {}
typename pnseq_generator<result_type>::result_type generate1() { return pn (eng); }
template<typename T>
object generate_impl (npy_intp n) {
nd::Array<T,1,1> out (n);
auto i = boost::begin(out);
auto e = boost::end(out);
for (; i < e; ++i)
*i = generate1();
return object (out);
}
object generate (npy_intp size, object const& dtype) {
NPY_TYPES type = dtype_to_type (dtype, NPY_INT);
if (type == NPY_BYTE)
return generate_impl<type_imap<NPY_BYTE>::type> (size);
else if (type == NPY_SHORT)
return generate_impl<type_imap<NPY_SHORT>::type> (size);
else if (type == NPY_INT)
return generate_impl<type_imap<NPY_INT>::type> (size);
else if (type == NPY_LONG)
return generate_impl<type_imap<NPY_LONG>::type> (size);
else
throw std::runtime_error ("pn::generate bad dtype");
}
object get_rs() const { return eng.rs; }
void set_rs(bp::object const& rs) { eng.rs = rs; }
int get_width() const { return pn.width; }
};
struct pn_pickle_suite : bp::pickle_suite {
static bp::tuple getinitargs (wrapped_pn<uint32_t> const& gen) {
return bp::make_tuple (gen.eng.rs, gen.pn.width, gen.pn.count, gen.pn.cache);
}
};
template<typename int_t>
static void export_pn (const char* name) {
class_<wrapped_pn<int_t> > (name,
bp::init<bp::object, int>((bp::arg("rng"),
bp::arg("bits")
),
"__init__(r,bits)\n"
"Construct a pn-generator\n\n"
":Parameters:\n"
" `r` : rng\n"
" Mersenne Twister\n"
" `bits` : int\n"
" #bits\n"
)
)
.def (bp::init<bp::object,int,int,typename pnseq_generator<int_t>::result_type>()) // for pickle
.def ("__call__", &wrapped_pn<int_t>::generate,
bp::default_call_policies(),
(bp::arg ("size"),
bp::arg ("dtype")=object()),
"__call__(cnt[,dtype])\n"
"Without cnt, generate a single r.v.\n"
"With cnt generate a vector of r.v.'s.\n\n"
":Returns: a single r.v. or a vector of r.v.'s\n"
" The returned values will be placed in the lowest bits of an integer.\n"
)
.def ("__call__", &wrapped_pn<int_t>::generate1)
.add_property ("width", &wrapped_pn<int_t>::get_width)
// .def_readwrite ("eng", &wrapped_pn<int_t>::eng)
.add_property ("rs", &wrapped_pn<int_t>::get_rs, &wrapped_pn<int_t>::set_rs)
.def_pickle (pn_pickle_suite())
;
}
BOOST_PYTHON_MODULE (pn) {
boost::numpy::initialize();
import_array();
export_pn<uint64_t> ("pn64");
export_pn<uint32_t> ("pn32");
object mod = object (handle<> ((PyImport_ImportModule("numpy.random"))));
rs_obj = mod.attr("RandomState");
}
_______________________________________________
NumPy-Discussion mailing list
[email protected]
http://mail.scipy.org/mailman/listinfo/numpy-discussion
