boost/python
2
0
Fork 0
mirror of https://github.com/boostorg/python.git synced 2026-01-26 18:52:26 +00:00
Code Issues Projects Releases Wiki Activity

removed ublas dependency from gaussian example

Browse Source
This commit is contained in:
Jim Bosch
2011-11-08 03:45:31 +00:00
parent 8f909d55ac
commit 1e66e33201
1 changed files with 106 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

144
libs/numpy/example/gaussian.cpp
View File

@@ -2,12 +2,85 @@
#include <memory>
#include <boost/numpy.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
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<double,2> vector;
typedef boost::numeric::ublas::c_matrix<double,2,2> 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<bivariate_gaussian const &>(self)().get_mu();
vector2 const & mu = bp::extract<bivariate_gaussian const &>(self)().get_mu();
return bn::from_data(
mu.data(),
bn::dtype::get_builtin<double>(),
@@ -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<bivariate_gaussian const &>(self)().get_sigma();
matrix2 const & sigma = bp::extract<bivariate_gaussian const &>(self)().get_sigma();
return bn::from_data(
sigma.data(),
bn::dtype::get_builtin<double>(),
@@ -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<double>(array[0]);
vec[1] = bp::extract<double>(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 <typename T, int N>
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<T> storage_t;
storage_t * storage = reinterpret_cast<storage_t*>(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>("bivariate_gaussian", bp::init<bivariate_gaussian const &>())
// 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)