Robert Kern wrote:
> On Wed, Mar 13, 2013 at 12:16 AM, Neal Becker <ndbeck...@gmail.com> 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.
>
I did find a way to do this, and the results are good enough. Timing is quite
comparable to my pure c++ implementation.
I used rk_ulong from mtrand.so. I also tried using rk_fill, but it was a bit
slower.
The boost::python c++ code is attached, for posterity.
#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);
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=uint64_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;
};
template<typename out_type=uint64_t>
struct engine {
typedef out_type result_type;
bp::object rs; // RandomState object
engine (bp::object const& _rs) : rs (_rs) {}
result_type operator() () {
RandomState* r = (RandomState*)(rs.ptr());
return rk_ulong (r->internal_state);
// result_type buffer;
// rk_fill ((void*)&buffer, sizeof(buffer), r->internal_state);
// return buffer;
}
};
template<typename result_type=uint64_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 = nd::allocate (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<uint64_t> const& gen) {
return bp::make_tuple (gen.eng.rs, gen.pn.width, gen.pn.count, gen.pn.cache);
}
};
BOOST_PYTHON_MODULE (pn64) {
boost::numpy::initialize();
import_array();
class_<engine<uint64_t> > ("engine", init<bp::object>())
.def ("__call__", &engine<uint64_t>::operator())
;
class_<wrapped_pn<uint64_t> > ("pn",
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,pnseq_generator<uint64_t>::result_type>()) // for pickle
.def ("__call__", &wrapped_pn<uint64_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<uint64_t>::generate1)
.add_property ("width", &wrapped_pn<uint64_t>::get_width)
// .def_readwrite ("eng", &wrapped_pn<uint64_t>::eng)
.add_property ("rs", &wrapped_pn<uint64_t>::get_rs, &wrapped_pn<uint64_t>::set_rs)
.def_pickle (pn_pickle_suite())
;
}
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
