[Boost-commit] svn:boost r75398 - sandbox/numpy/libs/numpy/example

Subject: [Boost-commit] svn:boost r75398 - sandbox/numpy/libs/numpy/example
From: talljimbo_at_[hidden]
Date: 2011-11-07 22:45:32

Author: jbosch
Date: 2011-11-07 22:45:31 EST (Mon, 07 Nov 2011)
New Revision: 75398
URL: http://svn.boost.org/trac/boost/changeset/75398

Log:
removed ublas dependency from gaussian example
Text files modified: 
   sandbox/numpy/libs/numpy/example/gaussian.cpp |   144 +++++++++++++++++++++++++++++---------- 
   1 files changed, 106 insertions(+), 38 deletions(-)

Modified: sandbox/numpy/libs/numpy/example/gaussian.cpp
==============================================================================
--- sandbox/numpy/libs/numpy/example/gaussian.cpp	(original)
+++ sandbox/numpy/libs/numpy/example/gaussian.cpp	2011-11-07 22:45:31 EST (Mon, 07 Nov 2011)
@@ -2,46 +2,109 @@
 #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'.
  */
 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;
-
-    vector const & get_mu() const { return _mu; }
+    vector2 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 @@
     /**
      *  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 @@
      *  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 @@
         return l;
     }
 
-    vector _mu;
-    matrix _sigma;
-    matrix _cholesky;
+    vector2 _mu;
+    matrix2 _sigma;
+    matrix2 _cholesky;
                         
 };
 
@@ -104,7 +167,7 @@
  *  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_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 @@
 
 /**
  *  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 @@
     }
 }
 
+/**
+ *  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 @@
         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 @@
     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)