From 1e66e33201c3f62a3e129d45884e815429ce2c3d Mon Sep 17 00:00:00 2001 From: Jim Bosch Date: Tue, 8 Nov 2011 03:45:31 +0000 Subject: [PATCH] removed ublas dependency from gaussian example --- libs/numpy/example/gaussian.cpp | 144 +++++++++++++++++++++++--------- 1 file changed, 106 insertions(+), 38 deletions(-) diff --git a/libs/numpy/example/gaussian.cpp b/libs/numpy/example/gaussian.cpp index c20eb30f..c6cf2ebd 100644 --- a/libs/numpy/example/gaussian.cpp +++ b/libs/numpy/example/gaussian.cpp @@ -2,12 +2,85 @@ #include #include -#include -#include namespace bp = boost::python; namespace bn = boost::numpy; +/** + * A 2x2 matrix class, purely for demonstration purposes. + * + * Instead of wrapping this class with Boost.Python, we'll convert it to/from numpy.ndarray. + */ +class matrix2 { +public: + + double & operator()(int i, int j) { + return _data[i*2 + j]; + } + + double const & operator()(int i, int j) const { + return _data[i*2 + j]; + } + + double const * data() const { return _data; } + +private: + double _data[4]; +}; + +/** + * A 2-element vector class, purely for demonstration purposes. + * + * Instead of wrapping this class with Boost.Python, we'll convert it to/from numpy.ndarray. + */ +class vector2 { +public: + + double & operator[](int i) { + return _data[i]; + } + + double const & operator[](int i) const { + return _data[i]; + } + + double const * data() const { return _data; } + + vector2 operator+(vector2 const & other) const { + vector2 r; + r[0] = _data[0] + other[0]; + r[1] = _data[1] + other[1]; + return r; + } + + vector2 operator-(vector2 const & other) const { + vector2 r; + r[0] = _data[0] - other[0]; + r[1] = _data[1] - other[1]; + return r; + } + +private: + double _data[2]; +}; + +/** + * Matrix-vector multiplication. + */ +vector2 operator*(matrix2 const & m, vector2 const & v) { + vector2 r; + r[0] = m(0, 0) * v[0] + m(0, 1) * v[1]; + r[1] = m(1, 0) * v[0] + m(1, 1) * v[1]; + return r; +} + +/** + * Vector inner product. + */ +double dot(vector2 const & v1, vector2 const & v2) { + return v1[0] * v2[0] + v1[1] * v2[1]; +} + /** * This class represents a simple 2-d Gaussian (Normal) distribution, defined by a * mean vector 'mu' and a covariance matrix 'sigma'. @@ -15,33 +88,23 @@ namespace bn = boost::numpy; class bivariate_gaussian { public: - /** - * Boost.NumPy isn't designed to support specific C++ linear algebra or matrix/vector libraries; - * it's intended as a lower-level interface that can be used with any such C++ library. - * - * Here, we'll demonstrate these techniques with boost::ublas, but the same general principles - * should apply to other matrix/vector libraries. - */ - typedef boost::numeric::ublas::c_vector vector; - typedef boost::numeric::ublas::c_matrix matrix; + vector2 const & get_mu() const { return _mu; } - vector const & get_mu() const { return _mu; } - - matrix const & get_sigma() const { return _sigma; } + matrix2 const & get_sigma() const { return _sigma; } /** * Evaluate the density of the distribution at a point defined by a two-element vector. */ - double operator()(vector const & p) const { - vector u = prod(_cholesky, p - _mu); - return 0.5 * _cholesky(0, 0) * _cholesky(1, 1) * std::exp(-0.5 * inner_prod(u, u)) / M_PI; + double operator()(vector2 const & p) const { + vector2 u = _cholesky * (p - _mu); + return 0.5 * _cholesky(0, 0) * _cholesky(1, 1) * std::exp(-0.5 * dot(u, u)) / M_PI; } /** * Evaluate the density of the distribution at an (x, y) point. */ double operator()(double x, double y) const { - vector p; + vector2 p; p[0] = x; p[1] = y; return operator()(p); @@ -50,7 +113,7 @@ public: /** * Construct from a mean vector and covariance matrix. */ - bivariate_gaussian(vector const & mu, matrix const & sigma) + bivariate_gaussian(vector2 const & mu, matrix2 const & sigma) : _mu(mu), _sigma(sigma), _cholesky(compute_inverse_cholesky(sigma)) {} @@ -60,8 +123,8 @@ private: * This evaluates the inverse of the Cholesky factorization of a 2x2 matrix; * it's just a shortcut in evaluating the density. */ - static matrix compute_inverse_cholesky(matrix const & m) { - matrix l; + static matrix2 compute_inverse_cholesky(matrix2 const & m) { + matrix2 l; // First do cholesky factorization: l l^t = m l(0, 0) = std::sqrt(m(0, 0)); l(0, 1) = m(0, 1) / l(0, 0); @@ -73,9 +136,9 @@ private: return l; } - vector _mu; - matrix _sigma; - matrix _cholesky; + vector2 _mu; + matrix2 _sigma; + matrix2 _cholesky; }; @@ -104,7 +167,7 @@ private: * and passing a const pointer to from_data causes NumPy's 'writeable' flag to be set to false. */ static bn::ndarray py_get_mu(bp::object const & self) { - bivariate_gaussian::vector const & mu = bp::extract(self)().get_mu(); + vector2 const & mu = bp::extract(self)().get_mu(); return bn::from_data( mu.data(), bn::dtype::get_builtin(), @@ -114,7 +177,7 @@ static bn::ndarray py_get_mu(bp::object const & self) { ); } static bn::ndarray py_get_sigma(bp::object const & self) { - bivariate_gaussian::matrix const & sigma = bp::extract(self)().get_sigma(); + matrix2 const & sigma = bp::extract(self)().get_sigma(); return bn::from_data( sigma.data(), bn::dtype::get_builtin(), @@ -126,24 +189,25 @@ static bn::ndarray py_get_sigma(bp::object const & self) { /** * To allow the constructor to work, we need to define some from-Python converters from NumPy arrays - * to the ublas types. The rvalue-from-python functionality is not well-documented in Boost.Python + * to the matrix/vector types. The rvalue-from-python functionality is not well-documented in Boost.Python * itself; you can learn more from boost/python/converter/rvalue_from_python_data.hpp. */ /** - * We start with two functions that just copy a NumPy array into ublas objects. These will be used + * We start with two functions that just copy a NumPy array into matrix/vector objects. These will be used * in the templated converted below. The first just uses the operator[] overloads provided by * bp::object. */ -static void copy_ndarray_to_ublas(bn::ndarray const & array, bivariate_gaussian::vector & vec) { +static void copy_ndarray_to_mv2(bn::ndarray const & array, vector2 & vec) { vec[0] = bp::extract(array[0]); vec[1] = bp::extract(array[1]); } + /** * Here, we'll take the alternate approach of using the strides to access the array's memory directly. * This can be much faster for large arrays. */ -static void copy_ndarray_to_ublas(bn::ndarray const & array, bivariate_gaussian::matrix & mat) { +static void copy_ndarray_to_mv2(bn::ndarray const & array, matrix2 & mat) { // Unfortunately, get_strides() can't be inlined, so it's best to call it once up-front. Py_intptr_t const * strides = array.get_strides(); for (int i = 0; i < 2; ++i) { @@ -153,13 +217,17 @@ static void copy_ndarray_to_ublas(bn::ndarray const & array, bivariate_gaussian: } } +/** + * Here's the actual converter. Because we've separated the differences into the above functions, + * we can write a single template class that works for both matrix2 and vector2. + */ template -struct bivariate_gaussian_ublas_from_python { +struct mv2_from_python { /** * Register the converter. */ - bivariate_gaussian_ublas_from_python() { + mv2_from_python() { bp::converter::registry::push_back( &convertible, &construct, @@ -198,9 +266,9 @@ struct bivariate_gaussian_ublas_from_python { typedef bp::converter::rvalue_from_python_storage storage_t; storage_t * storage = reinterpret_cast(data); // Use placement new to initialize the result. - T * ublas_obj = new (storage->storage.bytes) T(); + T * m_or_v = new (storage->storage.bytes) T(); // Fill the result with the values from the NumPy array. - copy_ndarray_to_ublas(*array, *ublas_obj); + copy_ndarray_to_mv2(*array, *m_or_v); // Finish up. data->convertible = storage->storage.bytes; } @@ -212,15 +280,15 @@ BOOST_PYTHON_MODULE(gaussian) { bn::initialize(); // Register the from-python converters - bivariate_gaussian_ublas_from_python< bivariate_gaussian::vector, 1 >(); - bivariate_gaussian_ublas_from_python< bivariate_gaussian::matrix, 2 >(); + mv2_from_python< vector2, 1 >(); + mv2_from_python< matrix2, 2 >(); - typedef double (bivariate_gaussian::*call_vector)(bivariate_gaussian::vector const &) const; + typedef double (bivariate_gaussian::*call_vector)(vector2 const &) const; bp::class_("bivariate_gaussian", bp::init()) // Declare the constructor (wouldn't work without the from-python converters). - .def(bp::init< bivariate_gaussian::vector const &, bivariate_gaussian::matrix const & >()) + .def(bp::init< vector2 const &, matrix2 const & >()) // Use our custom reference-counting getters .add_property("mu", &py_get_mu)