initial sandbox import for numpy utilities in boost.python

boost/python/numpy.hpp

@@ -0,0 +1,32 @@
+#ifndef BOOST_PYTHON_NUMPY_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy.hpp
+ *  @brief Main public header file for boost.python.numpy.
+ */
+
+#include <boost/python/numpy/dtype.hpp>
+#include <boost/python/numpy/ndarray.hpp>
+#include <boost/python/numpy/scalars.hpp>
+#include <boost/python/numpy/matrix.hpp>
+#include <boost/python/numpy/ufunc.hpp>
+
+namespace boost { namespace python {
+namespace numpy {
+
+/**
+ *  @brief Initialize the Numpy C-API
+ *
+ *  This must be called before using anything in boost.python.numpy;
+ *  It should probably be the first line inside BOOST_PYTHON_MODULE.
+ *
+ *  @internal This just calls the Numpy C-API functions "import_array()"
+ *            and "import_ufunc()".
+ */
+void initialize();
+
+} // namespace boost::python::numpy
+}} // namespace boost::python
+
+#endif // !BOOST_PYTHON_NUMPY_HPP_INCLUDED

boost/python/numpy/dtype.hpp

@@ -0,0 +1,56 @@
+#ifndef BOOST_PYTHON_NUMPY_DTYPE_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_DTYPE_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy/dtype.hpp
+ *  @brief Object manager for Python's numpy.dtype class.
+ */
+
+#include <boost/python.hpp>
+#include <boost/python/numpy/numpy_object_mgr_traits.hpp>
+
+namespace boost { namespace python {
+namespace numpy {
+
+/**
+ *  @brief A boost.python "object manager" (subclass of object) for numpy.dtype.
+ *
+ *  @todo This could have a lot more interesting accessors.
+ */
+class dtype : public object {
+    static python::detail::new_reference convert(object_cref arg, bool align);
+public:
+
+    /// @brief Convert an arbitrary Python object to a data-type descriptor object.
+    template <typename T>
+    explicit dtype(T arg, bool align=false) : object(convert(arg, align)) {}
+
+    /**
+     *  @brief Get the built-in numpy dtype associated with the given scalar template type.
+     *
+     *  This is perhaps the most useful part of the numpy API: it returns the dtype object
+     *  corresponding to a built-in C++ type.  This should work for any integer or floating point
+     *  type supported by numpy, and will also work for std::complex if 
+     *  sizeof(std::complex<T>) == 2*sizeof(T).
+     *
+     *  It can also be useful for users to add explicit specializations for POD structs
+     *  that return field-based dtypes.
+     */
+    template <typename T> static dtype get_builtin();
+
+    /// @brief Return the size of the data type in bytes.
+    int get_itemsize() const;
+
+    BOOST_PYTHON_FORWARD_OBJECT_CONSTRUCTORS(dtype, object);
+
+};
+
+} // namespace boost::python::numpy
+
+namespace converter {
+NUMPY_OBJECT_MANAGER_TRAITS(python::numpy::dtype);
+} // namespace boost::python::converter
+
+}} // namespace boost::python
+
+#endif // !BOOST_PYTHON_NUMPY_DTYPE_HPP_INCLUDED

boost/python/numpy/internal.hpp

@@ -0,0 +1,29 @@
+#ifndef BOOST_PYTHON_NUMPY_INTERNAL_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_INTERNAL_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy/internal.hpp
+ *  @brief Internal header file to include the Numpy C-API headers.
+ *
+ *  This should only be included by source files in the boost.python.numpy library itself.
+ */
+
+#include <boost/python.hpp>
+#ifdef BOOST_PYTHON_NUMPY_INTERNAL
+#define NO_IMPORT_ARRAY
+#define NO_IMPORT_UFUNC
+#else
+#ifndef BOOST_PYTHON_NUMPY_INTERNAL_MAIN
+ERROR_internal_hpp_is_for_internal_use_only
+#endif
+#endif
+#define PY_ARRAY_UNIQUE_SYMBOL BOOST_NUMPY_ARRAY_API
+#define PY_UFUNC_UNIQUE_SYMBOL BOOST_UFUNC_ARRAY_API
+#include <numpy/arrayobject.h>
+#include <numpy/ufuncobject.h>
+#include <boost/python/numpy.hpp>
+
+#define NUMPY_OBJECT_MANAGER_TRAITS_IMPL(pytype,manager)                \
+    PyTypeObject const * object_manager_traits<manager>::get_pytype() { return &pytype; }
+
+#endif // !BOOST_PYTHON_NUMPY_INTERNAL_HPP_INCLUDED

boost/python/numpy/matrix.hpp

@@ -0,0 +1,62 @@
+#ifndef BOOST_PYTHON_NUMPY_MATRIX_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_MATRIX_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy/matrix.hpp
+ *  @brief Object manager for numpy.matrix.
+ */
+
+#include <boost/python.hpp>
+#include <boost/python/numpy/numpy_object_mgr_traits.hpp>
+#include <boost/python/numpy/ndarray.hpp>
+
+namespace boost { namespace python {
+
+namespace numpy {
+
+/**
+ *  @brief A boost.python "object manager" (subclass of object) for numpy.matrix.
+ *
+ *  @internal numpy.matrix is defined in Python, so object_manager_traits<matrix>::get_pytype()
+ *            is implemented by importing numpy and getting the "matrix" attribute of the module.
+ *            We then just hope that doesn't get destroyed while we need it, because if we put
+ *            a dynamic python object in a static-allocated boost::python::object or handle<>,
+ *            bad things happen when Python shuts down.  I think this solution is safe, but I'd
+ *            love to get that confirmed.
+ */
+class matrix : public ndarray {
+    static object construct(object_cref obj, dtype const & dt, bool copy);
+    static object construct(object_cref obj, bool copy);
+public:
+
+    BOOST_PYTHON_FORWARD_OBJECT_CONSTRUCTORS(matrix, ndarray);
+
+    /// @brief Equivalent to "numpy.matrix(obj,dt,copy)" in Python.
+    explicit matrix(object const & obj, dtype const & dt, bool copy=true) :
+        ndarray(extract<ndarray>(construct(obj, dt, copy))) {}
+
+    /// @brief Equivalent to "numpy.matrix(obj,copy=copy)" in Python.
+    explicit matrix(object const & obj, bool copy=true) :
+        ndarray(extract<ndarray>(construct(obj, copy))) {}
+
+    /// \brief Return a view of the matrix with the given dtype.
+    matrix view(dtype const & dt) const;
+
+    /// \brief Copy the scalar (deep for all non-object fields).
+    matrix copy() const;
+
+    /// \brief Transpose the matrix.
+    matrix transpose() const;
+
+};
+
+} // namespace boost::python::numpy
+
+namespace converter {
+
+NUMPY_OBJECT_MANAGER_TRAITS(python::numpy::matrix);
+
+} // namespace boost::python::converter
+}} // namespace boost::python
+
+#endif // !BOOST_PYTHON_NUMPY_MATRIX_HPP_INCLUDED

boost/python/numpy/ndarray.hpp

@@ -0,0 +1,285 @@
+#ifndef BOOST_PYTHON_NUMPY_NDARRAY_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_NDARRAY_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy/ndarray.hpp
+ *  @brief Object manager and various utilities for numpy.ndarray.
+ */
+
+#include <boost/python.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <boost/type_traits/is_integral.hpp>
+#include <boost/python/numpy/numpy_object_mgr_traits.hpp>
+#include <boost/python/numpy/dtype.hpp>
+
+#include <vector>
+
+namespace boost { namespace python {
+namespace numpy {
+
+/**
+ *  @brief A boost.python "object manager" (subclass of object) for numpy.ndarray.
+ *
+ *  @todo This could have a lot more functionality (like boost::python::numeric::array).
+ *        Right now all that exists is what was needed to move raw data between C++ and Python.
+ */
+class ndarray : public object {
+
+    /**
+     *  @brief An internal struct that's byte-compatible with PyArrayObject.
+     *
+     *  This is just a hack to allow inline access to this stuff while hiding numpy/arrayobject.h
+     *  from the user.
+     */
+    struct array_struct {
+        PyObject_HEAD
+        char * data;
+        int nd;
+        Py_intptr_t * shape;
+        Py_intptr_t * strides;
+        PyObject * base;
+        PyObject * descr;
+        int flags;
+        PyObject * weakreflist;
+    };
+
+    /// @brief Return the held Python object as an array_struct.
+    array_struct * get_struct() const { return reinterpret_cast<array_struct*>(this->ptr()); }
+
+public:
+
+    /**
+     *  @brief Enum to represent (some) of Numpy's internal flags.
+     *
+     *  These don't match the actual Numpy flag values; we can't get those without including 
+     *  numpy/arrayobject.h or copying them directly.  That's very unfortunate.
+     *
+     *  @todo I'm torn about whether this should be an enum.  It's very convenient to not
+     *        make these simple integer values for overloading purposes, but the need to
+     *        define every possible combination and custom bitwise operators is ugly.
+     */
+    enum bitflag {
+        NONE=0x0, C_CONTIGUOUS=0x1, F_CONTIGUOUS=0x2, V_CONTIGUOUS=0x1|0x2, 
+        ALIGNED=0x4, WRITEABLE=0x8, BEHAVED=0x4|0x8,
+        CARRAY_RO=0x1|0x4, CARRAY=0x1|0x4|0x8, CARRAY_MIS=0x1|0x8,
+        FARRAY_RO=0x2|0x4, FARRAY=0x2|0x4|0x8, FARRAY_MIS=0x2|0x8,
+        UPDATE_ALL=0x1|0x2|0x4, VARRAY=0x1|0x2|0x8, ALL=0x1|0x2|0x4|0x8
+    };
+
+    BOOST_PYTHON_FORWARD_OBJECT_CONSTRUCTORS(ndarray, object);
+
+    /// @brief Return a view of the scalar with the given dtype.
+    ndarray view(dtype const & dt) const;
+
+    /// @brief Copy the scalar (deep for all non-object fields).
+    ndarray copy() const;
+
+    /// @brief Return the size of the nth dimension.
+    int const shape(int n) const { return get_shape()[n]; }
+
+    /// @brief Return the stride of the nth dimension.
+    int const strides(int n) const { return get_strides()[n]; }
+    
+    /**
+     *  @brief Return the array's raw data pointer.
+     *
+     *  This returns char so stride math works properly on it.  It's pretty much
+     *  expected that the user will have to reinterpret_cast it.
+     */
+    char * get_data() const { return get_struct()->data; }
+
+    /// @brief Return the array's data-type descriptor object.
+    dtype get_dtype() const;
+
+    /// @brief Return the object that owns the array's data, or None if the array owns its own data.
+    object get_base() const;
+
+    /// @brief Set the object that owns the array's data.  Use with care.
+    void set_base(object const & base);
+
+    /// @brief Return the shape of the array as an array of integers (length == get_nd()).
+    Py_intptr_t const * get_shape() const { return get_struct()->shape; }
+
+    /// @brief Return the stride of the array as an array of integers (length == get_nd()).
+    Py_intptr_t const * get_strides() const { return get_struct()->strides; }
+
+    /// @brief Return the number of array dimensions.
+    int const get_nd() const { return get_struct()->nd; }
+
+    /// @brief Return the array flags.
+    bitflag const get_flags() const;
+
+    /// @brief Reverse the dimensions of the array.
+    ndarray transpose() const;
+
+    /// @brief Eliminate any unit-sized dimensions.
+    ndarray squeeze() const;
+
+    /**
+     *  @brief If the array contains only a single element, return it as an array scalar; otherwise return
+     *         the array.
+     *
+     *  @internal This is simply a call to PyArray_Return();
+     */
+    object scalarize() const;
+};
+
+/**
+ *  @brief Construct a new array with the given shape and data type, with data initialized to zero.
+ */
+ndarray zeros(tuple const & shape, dtype const & dt);
+ndarray zeros(int nd, Py_intptr_t const * shape, dtype const & dt);
+
+/**
+ *  @brief Construct a new array with the given shape and data type, with data left uninitialized.
+ */
+ndarray empty(tuple const & shape, dtype const & dt);
+ndarray empty(int nd, Py_intptr_t const * shape, dtype const & dt);
+
+/**
+ *  @brief Construct a new array from an arbitrary Python sequence.
+ *
+ *  @todo This does't seem to handle ndarray subtypes the same way that "numpy.array" does in Python.
+ */
+ndarray array(object const & obj);
+ndarray array(object const & obj, dtype const & dt);
+
+namespace detail {
+
+ndarray from_data_impl(
+    void * data,
+    dtype const & dt,
+    std::vector<Py_intptr_t> const & shape,
+    std::vector<Py_intptr_t> const & strides,
+    object const & owner,
+    bool writeable
+);
+
+template <typename Container>
+ndarray from_data_impl(
+    void * data,
+    dtype const & dt,
+    Container shape,
+    Container strides,
+    object const & owner,
+    bool writeable,
+    typename boost::enable_if< boost::is_integral<typename Container::value_type> >::type * enabled = NULL
+) {
+    std::vector<Py_intptr_t> shape_(shape.begin(),shape.end());
+    std::vector<Py_intptr_t> strides_(strides.begin(), strides.end());
+    return from_data_impl(data, dt, shape_, strides_, owner, writeable);    
+}
+
+ndarray from_data_impl(
+    void * data,
+    dtype const & dt,
+    object const & shape,
+    object const & strides,
+    object const & owner,
+    bool writeable
+);
+
+} // namespace boost::python::numpy::detail
+
+/**
+ *  @brief Construct a new ndarray object from a raw pointer.
+ *
+ *  @param[in] data    Raw pointer to the first element of the array.
+ *  @param[in] dt      Data type descriptor.  Often retrieved with dtype::get_builtin().
+ *  @param[in] shape   Shape of the array as STL container of integers; must have begin() and end().
+ *  @param[in] strides Shape of the array as STL container of integers; must have begin() and end().
+ *  @param[in] owner   An arbitray Python object that owns that data pointer.  The array object will
+ *                     keep a reference to the object, and decrement it's reference count when the
+ *                     array goes out of scope.  Pass None at your own peril.
+ *
+ *  @todo Should probably take ranges of iterators rather than actual container objects.
+ */
+template <typename Container>
+inline ndarray from_data(
+    void * data,
+    dtype const & dt,
+    Container shape,
+    Container strides,
+    object const & owner
+) {
+    return numpy::detail::from_data_impl(data, dt, shape, strides, owner, true);
+}    
+
+/**
+ *  @brief Construct a new ndarray object from a raw pointer.
+ *
+ *  @param[in] data    Raw pointer to the first element of the array.
+ *  @param[in] dt      Data type descriptor.  Often retrieved with dtype::get_builtin().
+ *  @param[in] shape   Shape of the array as STL container of integers; must have begin() and end().
+ *  @param[in] strides Shape of the array as STL container of integers; must have begin() and end().
+ *  @param[in] owner   An arbitray Python object that owns that data pointer.  The array object will
+ *                     keep a reference to the object, and decrement it's reference count when the
+ *                     array goes out of scope.  Pass None at your own peril.
+ *
+ *  This overload takes a const void pointer and sets the "writeable" flag of the array to false.
+ *
+ *  @todo Should probably take ranges of iterators rather than actual container objects.
+ */
+template <typename Container>
+inline ndarray from_data(
+    void const * data,
+    dtype const & dt,
+    Container shape,
+    Container strides,
+    object const & owner
+) {
+    return numpy::detail::from_data_impl(const_cast<void*>(data), dt, shape, strides, owner, false);
+}    
+
+/**
+ *  @brief Transform an arbitrary object into a numpy array with the given requirements.
+ *
+ *  @param[in] obj     An arbitrary python object to convert.  Arrays that meet the requirements
+ *                     will be passed through directly.
+ *  @param[in] dt      Data type descriptor.  Often retrieved with dtype::get_builtin().
+ *  @param[in] nd_min  Minimum number of dimensions.
+ *  @param[in] nd_max  Maximum number of dimensions.
+ *  @param[in] flags   Bitwise OR of flags specifying additional requirements.
+ */
+ndarray from_object(object const & obj, dtype const & dt, 
+                    int nd_min, int nd_max, ndarray::bitflag flags=ndarray::NONE);
+
+inline ndarray from_object(object const & obj, dtype const & dt, 
+                           int nd, ndarray::bitflag flags=ndarray::NONE) {
+    return from_object(obj, dt, nd, nd, flags);
+}
+
+inline ndarray from_object(object const & obj, dtype const & dt, ndarray::bitflag flags=ndarray::NONE) {
+    return from_object(obj, dt, 0, 0, flags);
+}
+
+ndarray from_object(object const & obj, int nd_min, int nd_max, 
+                    ndarray::bitflag flags=ndarray::NONE);
+
+inline ndarray from_object(object const & obj, int nd, ndarray::bitflag flags=ndarray::NONE) {
+    return from_object(obj, nd, nd, flags);
+}
+
+inline ndarray from_object(object const & obj, ndarray::bitflag flags=ndarray::NONE) {
+    return from_object(obj, 0, 0, flags);
+}
+
+inline ndarray::bitflag operator|(ndarray::bitflag a, ndarray::bitflag b) {
+    return ndarray::bitflag(int(a) | int(b));
+}
+
+inline ndarray::bitflag operator&(ndarray::bitflag a, ndarray::bitflag b) {
+    return ndarray::bitflag(int(a) & int(b));
+}
+
+
+} // namespace boost::python::numpy
+
+namespace converter {
+
+NUMPY_OBJECT_MANAGER_TRAITS(python::numpy::ndarray);
+
+} // namespace boost::python::converter
+}} // namespace boost::python
+
+#endif // !BOOST_PYTHON_NUMPY_NDARRAY_HPP_INCLUDED

boost/python/numpy/numpy_object_mgr_traits.hpp

@@ -0,0 +1,27 @@
+#ifndef BOOST_PYTHON_NUMPY_NUMPY_OBJECT_MGR_TRAITS_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_NUMPY_OBJECT_MGR_TRAITS_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy/numpy_object_mgr_traits.hpp
+ *  @brief Macro that specializes object_manager_traits by requiring a 
+ *         source-file implementation of get_pytype().
+ */
+
+#define NUMPY_OBJECT_MANAGER_TRAITS(manager)                            \
+    template <>                                                         \
+    struct object_manager_traits<manager> {                             \
+        BOOST_STATIC_CONSTANT(bool, is_specialized = true);             \
+        static inline python::detail::new_reference adopt(PyObject* x) { \
+            return python::detail::new_reference(python::pytype_check((PyTypeObject*)get_pytype(), x)); \
+        }                                                               \
+        static bool check(PyObject* x) {                                \
+            return ::PyObject_IsInstance(x, (PyObject*)get_pytype());   \
+        }                                                               \
+        static manager* checked_downcast(PyObject* x) {                 \
+            return python::downcast<manager>((checked_downcast_impl)(x, (PyTypeObject*)get_pytype())); \
+        }                                                               \
+        static PyTypeObject const * get_pytype();                       \
+    }
+
+#endif // !BOOST_PYTHON_NUMPY_NUMPY_OBJECT_MGR_TRAITS_HPP_INCLUDED
+

boost/python/numpy/scalars.hpp

@@ -0,0 +1,55 @@
+#ifndef BOOST_PYTHON_NUMPY_SCALARS_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_SCALARS_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy/scalars.hpp
+ *  @brief Object managers for array scalars (currently only numpy.void is implemented).
+ */
+
+#include <boost/python.hpp>
+#include <boost/python/numpy/numpy_object_mgr_traits.hpp>
+#include <boost/python/numpy/dtype.hpp>
+
+namespace boost { namespace python {
+namespace numpy {
+
+/**
+ *  @brief A boost.python "object manager" (subclass of object) for numpy.void.
+ *
+ *  @todo This could have a lot more functionality.
+ */
+class void_ : public object {
+    static python::detail::new_reference convert(object_cref arg, bool align);
+public:
+
+    /**
+     *  @brief Construct a new array scalar with the given size and void dtype.
+     *
+     *  Data is initialized to zero.  One can create a standalone scalar object
+     *  with a certain dtype "dt" with:
+     *  @code
+     *  void_ scalar = void_(dt.get_itemsize()).view(dt);
+     *  @endcode
+     */
+    explicit void_(Py_ssize_t size);
+
+    BOOST_PYTHON_FORWARD_OBJECT_CONSTRUCTORS(void_, object);
+
+    /// @brief Return a view of the scalar with the given dtype.
+    void_ view(dtype const & dt) const;
+
+    /// @brief Copy the scalar (deep for all non-object fields).
+    void_ copy() const;
+
+};
+
+} // namespace boost::python::numpy
+
+namespace converter {
+
+NUMPY_OBJECT_MANAGER_TRAITS(python::numpy::void_);
+
+} // namespace boost::python::converter
+}} // namespace boost::python
+
+#endif // !BOOST_PYTHON_NUMPY_SCALARS_HPP_INCLUDED

boost/python/numpy/ufunc.hpp

@@ -0,0 +1,193 @@
+#ifndef BOOST_PYTHON_NUMPY_UFUNC_HPP_INCLUDED
+#define BOOST_PYTHON_NUMPY_UFUNC_HPP_INCLUDED
+
+/**
+ *  @file boost/python/numpy/ufunc.hpp
+ *  @brief Utilities to create ufunc-like broadcasting functions out of C++ functors.
+ */
+
+#include <boost/python.hpp>
+#include <boost/python/numpy/numpy_object_mgr_traits.hpp>
+#include <boost/python/numpy/dtype.hpp>
+#include <boost/python/numpy/ndarray.hpp>
+
+namespace boost { namespace python {
+namespace numpy {
+
+/**
+ *  @brief A boost.python "object manager" (subclass of object) for PyArray_MultiIter.
+ *
+ *  multi_iter is a Python object, but a very low-level one.  It should generally only be used
+ *  in loops of the form:
+ *  @code
+ *  while (iter.not_done()) {
+ *      ...
+ *      iter.next();
+ *  }
+ *  @endcode
+ *
+ *  @todo I can't tell if this type is exposed in Python anywhere; if it is, we should use that name.
+ *        It's more dangerous than most object managers, however - maybe it actually belongs in
+ *        a detail namespace?
+ */
+class multi_iter : public object {
+public:
+
+    BOOST_PYTHON_FORWARD_OBJECT_CONSTRUCTORS(multi_iter, object);
+
+    /// @brief Increment the iterator.
+    void next();
+
+    /// @brief Check if the iterator is at its end.
+    bool not_done() const;
+
+    /// @brief Return a pointer to the element of the nth broadcasted array.
+    char * get_data(int n) const;
+
+    /// @brief Return the number of dimensions of the broadcasted array expression.
+    int const get_nd() const;
+    
+    /// @brief Return the shape of the broadcasted array expression as an array of integers.
+    Py_intptr_t const * get_shape() const;
+
+    /// @brief Return the shape of the broadcasted array expression in the nth dimension.
+    Py_intptr_t const shape(int n) const;
+    
+};
+
+/// @brief Construct a multi_iter over a single sequence or scalar object.
+multi_iter make_multi_iter(object const & a1);
+
+/// @brief Construct a multi_iter by broadcasting two objects.
+multi_iter make_multi_iter(object const & a1, object const & a2);
+
+/// @brief Construct a multi_iter by broadcasting three objects.
+multi_iter make_multi_iter(object const & a1, object const & a2, object const & a3);
+
+/**
+ *  @brief Helps wrap a C++ functor taking a single scalar argument as a broadcasting ufunc-like
+ *         Python object.
+ *
+ *  Typical usage looks like this:
+ *  @code
+ *  struct TimesPI {
+ *      typedef double argument_type;
+ *      typedef double result_type;
+ *      double operator()(double input) const { return input * M_PI; }
+ *  };
+ *  
+ *  BOOST_PYTHON_MODULE(example) {
+ *      class_< TimesPI >("TimesPI")
+ *          .def("__call__", unary_ufunc<TimesPI>::make())
+ *          ;
+ *  }
+ *  @endcode
+ *  
+ */
+template <typename TUnaryFunctor, 
+          typename TArgument=typename TUnaryFunctor::argument_type,
+          typename TResult=typename TUnaryFunctor::result_type>
+struct unary_ufunc {
+
+    /**
+     *  @brief A C++ function with object arguments that broadcasts its arguments before
+     *         passing them to the underlying C++ functor.
+     */
+    static object call(TUnaryFunctor & self, object const & input, object const & output) {
+        dtype in_dtype = dtype::get_builtin<TArgument>();
+        dtype out_dtype = dtype::get_builtin<TResult>();
+        ndarray in_array = from_object(input, in_dtype, ndarray::ALIGNED);
+        ndarray out_array = (output != object()) ? 
+            from_object(output, out_dtype, ndarray::ALIGNED | ndarray::WRITEABLE)
+            : zeros(in_array.get_nd(), in_array.get_shape(), out_dtype);
+        multi_iter iter = make_multi_iter(in_array, out_array);
+        while (iter.not_done()) {
+            TArgument * argument = reinterpret_cast<TArgument*>(iter.get_data(0));
+            TResult * result = reinterpret_cast<TResult*>(iter.get_data(1));
+            *result = self(*argument);
+            iter.next();
+        } 
+        return out_array.scalarize();
+    }
+
+    /**
+     *  @brief Construct a boost.python function object from call() with reasonable keyword names.
+     *
+     *  Users will often want to specify their own keyword names with the same signature, but this
+     *  is a convenient shortcut.
+     */
+    static object make() {
+        return make_function(call, default_call_policies(), (arg("input"), arg("output")=object())); 
+    }
+
+};
+
+/**
+ *  @brief Helps wrap a C++ functor taking a pair of scalar arguments as a broadcasting ufunc-like
+ *         Python object.
+ *
+ *  Typical usage looks like this:
+ *  @code
+ *  struct CosSum {
+ *      typedef double first_argument_type;
+ *      typedef double second_argument_type;
+ *      typedef double result_type;
+ *      double operator()(double input1, double input2) const { return std::cos(input1 + input2); }
+ *  };
+ *  
+ *  BOOST_PYTHON_MODULE(example) {
+ *      class_< CosSum >("CosSum")
+ *          .def("__call__", binary_ufunc<CosSum>::make())
+ *          ;
+ *  }
+ *  @endcode
+ *  
+ */
+template <typename TBinaryFunctor, 
+          typename TArgument1=typename TBinaryFunctor::first_argument_type,
+          typename TArgument2=typename TBinaryFunctor::second_argument_type,
+          typename TResult=typename TBinaryFunctor::result_type>
+struct binary_ufunc {
+
+    static object call(TBinaryFunctor & self, object const & input1, object const & input2, 
+                       object const & output) 
+    {
+        dtype in1_dtype = dtype::get_builtin<TArgument1>();
+        dtype in2_dtype = dtype::get_builtin<TArgument2>();
+        dtype out_dtype = dtype::get_builtin<TResult>();
+        ndarray in1_array = from_object(input1, in1_dtype, ndarray::ALIGNED);
+        ndarray in2_array = from_object(input2, in2_dtype, ndarray::ALIGNED);
+        multi_iter iter = make_multi_iter(in1_array, in2_array);
+        ndarray out_array = (output != object()) ? 
+            from_object(output, out_dtype, ndarray::ALIGNED | ndarray::WRITEABLE)
+            : zeros(iter.get_nd(), iter.get_shape(), out_dtype);
+        iter = make_multi_iter(in1_array, in2_array, out_array);
+        while (iter.not_done()) {
+            TArgument1 * argument1 = reinterpret_cast<TArgument1*>(iter.get_data(0));
+            TArgument2 * argument2 = reinterpret_cast<TArgument2*>(iter.get_data(1));
+            TResult * result = reinterpret_cast<TResult*>(iter.get_data(2));
+            *result = self(*argument1, *argument2);
+            iter.next();
+        } 
+        return out_array.scalarize();
+    }
+
+    static object make() {
+        return make_function(
+            call, default_call_policies(), 
+            (arg("input1"), arg("input2"), arg("output")=object())
+        ); 
+    }
+
+};
+
+} // namespace boost::python::numpy
+
+namespace converter {
+
+NUMPY_OBJECT_MANAGER_TRAITS(python::numpy::multi_iter);
+
+} // namespace boost::python::converter
+}} // namespace boost::python
+
+#endif // !BOOST_PYTHON_NUMPY_UFUNC_HPP_INCLUDED

libs/python/numpy/src/dtype.cpp

@@ -0,0 +1,84 @@
+#define BOOST_PYTHON_NUMPY_INTERNAL
+#include <boost/python/numpy/internal.hpp>
+
+#define NUMPY_DTYPE_TRAITS_BUILTIN(ctype,code)                          \
+    template <> struct dtype_traits<ctype> {                            \
+        static dtype get() {                                            \
+            return dtype(                                               \
+                python::detail::new_reference(                          \
+                    reinterpret_cast<PyObject*>(PyArray_DescrFromType(code)) \
+                )                                                       \
+            );                                                          \
+        }                                                               \
+    };                                                                  \
+    template dtype dtype::get_builtin<ctype>()
+
+#define NUMPY_DTYPE_TRAITS_COMPLEX(creal, ctype, code)                  \
+    template <> struct dtype_traits< std::complex<creal> > {            \
+        static dtype get() {                                            \
+            if (sizeof(ctype) != sizeof(std::complex<creal>)) {         \
+                PyErr_SetString(PyExc_TypeError, "Cannot reinterpret std::complex<T> as T[2]"); \
+                throw_error_already_set();                              \
+            }                                                           \
+            return dtype(                                               \
+                python::detail::new_reference(                          \
+                    reinterpret_cast<PyObject*>(PyArray_DescrFromType(code)) \
+                )                                                       \
+            );                                                          \
+        }                                                               \
+    };                                                                  \
+    template dtype dtype::get_builtin< std::complex<creal> >()
+
+namespace boost { namespace python {
+namespace numpy {
+
+template <typename T> struct dtype_traits;
+
+python::detail::new_reference dtype::convert(object const & arg, bool align) {
+    PyArray_Descr* obj=NULL;
+    if (align) {
+        if (PyArray_DescrAlignConverter(arg.ptr(), &obj) < 0)
+            throw_error_already_set();
+    } else {
+        if (PyArray_DescrConverter(arg.ptr(), &obj) < 0)
+            throw_error_already_set();
+    }
+    return python::detail::new_reference(reinterpret_cast<PyObject*>(obj));
+}
+
+int dtype::get_itemsize() const {
+    return reinterpret_cast<PyArray_Descr*>(ptr())->elsize;
+}
+
+template <typename T>
+dtype dtype::get_builtin() { return dtype_traits<T>::get(); }
+
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_ubyte, NPY_UBYTE);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_byte, NPY_BYTE);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_ushort, NPY_USHORT);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_short, NPY_SHORT);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_uint, NPY_UINT);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_int, NPY_INT);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_ulong, NPY_ULONG);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_long, NPY_LONG);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_float, NPY_FLOAT);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_double, NPY_DOUBLE);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_longdouble, NPY_LONGDOUBLE);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_cfloat, NPY_CFLOAT);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_cdouble, NPY_CDOUBLE);
+NUMPY_DTYPE_TRAITS_BUILTIN(npy_clongdouble, NPY_CLONGDOUBLE);
+NUMPY_DTYPE_TRAITS_COMPLEX(float, npy_cfloat, NPY_CFLOAT);
+NUMPY_DTYPE_TRAITS_COMPLEX(double, npy_cdouble, NPY_CDOUBLE);
+NUMPY_DTYPE_TRAITS_COMPLEX(long double, npy_clongdouble, NPY_CLONGDOUBLE);
+
+template <> struct dtype_traits<bool> {
+    static dtype get() {
+	if (sizeof(bool) == sizeof(npy_bool)) return dtype_traits<npy_bool>::get();
+	if (sizeof(bool) == sizeof(npy_ubyte)) return dtype_traits<npy_ubyte>::get();
+	PyErr_SetString(PyExc_TypeError, "Cannot determine numpy dtype corresponding to C++ bool.");
+        throw_error_already_set();
+    }
+};
+template dtype dtype::get_builtin<bool>();
+
+}}} // namespace boost::python::numpy

libs/python/numpy/src/matrix.cpp

@@ -0,0 +1,51 @@
+#define BOOST_PYTHON_NUMPY_INTERNAL
+#include <boost/python/numpy/internal.hpp>
+#include <boost/python/numpy/matrix.hpp>
+
+namespace boost { namespace python {
+namespace numpy { namespace detail {
+inline object get_matrix_type() {
+    object module = import("numpy");
+    return module.attr("matrix");
+}
+}} // namespace numpy::detail
+
+namespace converter {
+
+PyTypeObject const * object_manager_traits<numpy::matrix>::get_pytype() {
+    return reinterpret_cast<PyTypeObject*>(numpy::detail::get_matrix_type().ptr());
+}
+
+} // namespace boost::python::converter
+
+namespace numpy {
+
+object matrix::construct(object const & obj, dtype const & dt, bool copy) {
+    return numpy::detail::get_matrix_type()(obj, dt, copy);
+}
+
+object matrix::construct(object const & obj, bool copy) {
+    return numpy::detail::get_matrix_type()(obj, object(), copy);
+}
+
+matrix matrix::view(dtype const & dt) const {
+    return matrix(
+        python::detail::new_reference(
+            PyObject_CallMethod(this->ptr(), const_cast<char*>("view"), const_cast<char*>("O"), dt.ptr())
+        )
+    );
+}
+
+matrix matrix::copy() const {
+    return matrix(
+        python::detail::new_reference(
+            PyObject_CallMethod(this->ptr(), const_cast<char*>("copy"), const_cast<char*>(""))
+        )
+    );
+}
+
+matrix matrix::transpose() const {
+    return matrix(extract<matrix>(ndarray::transpose()));
+}
+
+}}} // namespace boost::python::numpy

libs/python/numpy/src/ndarray.cpp

@@ -0,0 +1,277 @@
+#define BOOST_PYTHON_NUMPY_INTERNAL
+#include <boost/python/numpy/internal.hpp>
+
+namespace boost { namespace python {
+namespace converter {
+NUMPY_OBJECT_MANAGER_TRAITS_IMPL(PyArray_Type, python::numpy::ndarray)
+} // namespace boost::python::converter
+
+namespace numpy {
+
+namespace detail {
+
+ndarray::bitflag numpy_to_bitflag(int const f) {
+    ndarray::bitflag r = ndarray::NONE;
+    if (f & NPY_C_CONTIGUOUS) r = (r | ndarray::C_CONTIGUOUS);
+    if (f & NPY_F_CONTIGUOUS) r = (r | ndarray::F_CONTIGUOUS);
+    if (f & NPY_ALIGNED) r = (r | ndarray::ALIGNED);
+    if (f & NPY_WRITEABLE) r = (r | ndarray::WRITEABLE);
+    return r;
+}
+
+int const bitflag_to_numpy(ndarray::bitflag f) {
+    int r = 0;
+    if (f & ndarray::C_CONTIGUOUS) r |= NPY_C_CONTIGUOUS;
+    if (f & ndarray::F_CONTIGUOUS) r |= NPY_F_CONTIGUOUS;
+    if (f & ndarray::ALIGNED) r |= NPY_ALIGNED;
+    if (f & ndarray::WRITEABLE) r |= NPY_WRITEABLE;
+    return r;
+}
+
+bool is_c_contiguous(
+    std::vector<Py_intptr_t> const & shape,
+    std::vector<Py_intptr_t> const & strides,
+    int itemsize
+) {
+    std::vector<Py_intptr_t>::const_reverse_iterator j = strides.rbegin();
+    int total = itemsize;
+    for (std::vector<Py_intptr_t>::const_reverse_iterator i = shape.rbegin(); i != shape.rend(); ++i, ++j) {
+        if (total != *j) return false;
+        total *= (*i);
+    }
+    return true;
+}
+
+bool is_f_contiguous(
+    std::vector<Py_intptr_t> const & shape,
+    std::vector<Py_intptr_t> const & strides,
+    int itemsize
+) {
+    std::vector<Py_intptr_t>::const_iterator j = strides.begin();
+    int total = itemsize;
+    for (std::vector<Py_intptr_t>::const_iterator i = shape.begin(); i != shape.end(); ++i, ++j) {
+        if (total != *j) return false;
+        total *= (*i);
+    }
+    return true;
+}
+
+bool is_aligned(
+    std::vector<Py_intptr_t> const & strides,
+    int itemsize
+) {
+    for (std::vector<Py_intptr_t>::const_iterator i = strides.begin(); i != strides.end(); ++i) {
+        if (*i % itemsize) return false;
+    }
+    return true;
+}
+
+inline PyArray_Descr * incref_dtype(dtype const & dt) {
+    Py_INCREF(dt.ptr());
+    return reinterpret_cast<PyArray_Descr*>(dt.ptr());
+}
+
+ndarray from_data_impl(
+    void * data, 
+    dtype const & dt, 
+    object const & shape,
+    object const & strides,
+    object const & owner,
+    bool writeable
+) {
+    std::vector<Py_intptr_t> shape_(len(shape));
+    std::vector<Py_intptr_t> strides_(len(strides));
+    if (shape_.size() != strides_.size()) {
+        PyErr_SetString(PyExc_ValueError, "Length of shape and strides arrays do not match.");
+        throw_error_already_set();
+    }
+    for (std::size_t i = 0; i < shape_.size(); ++i) {
+        shape_[i] = extract<Py_intptr_t>(shape[i]);
+        strides_[i] = extract<Py_intptr_t>(strides[i]);
+    }
+    return from_data_impl(data, dt, shape_, strides_, owner, writeable);
+}
+
+ndarray from_data_impl(
+    void * data, 
+    dtype const & dt, 
+    std::vector<Py_intptr_t> const & shape,
+    std::vector<Py_intptr_t> const & strides,
+    object const & owner,
+    bool writeable
+) {
+    if (shape.size() != strides.size()) {
+        PyErr_SetString(PyExc_ValueError, "Length of shape and strides arrays do not match.");
+        throw_error_already_set();
+    }
+    int itemsize = dt.get_itemsize();
+    int flags = 0;
+    if (writeable) flags |= NPY_WRITEABLE;
+    if (is_c_contiguous(shape, strides, itemsize)) flags |= NPY_C_CONTIGUOUS;
+    if (is_f_contiguous(shape, strides, itemsize)) flags |= NPY_F_CONTIGUOUS;
+    if (is_aligned(strides, itemsize)) flags |= NPY_ALIGNED;
+    ndarray r(
+        python::detail::new_reference(
+            PyArray_NewFromDescr(
+                &PyArray_Type,
+                incref_dtype(dt),
+                shape.size(),
+                const_cast<Py_intptr_t*>(&shape.front()),
+                const_cast<Py_intptr_t*>(&strides.front()),
+                data,
+                flags,
+                NULL
+            )
+        )
+    );
+    r.set_base(owner);
+    return r;
+}
+
+} // namespace detail
+
+ndarray ndarray::view(dtype const & dt) const {
+    return ndarray(
+        python::detail::new_reference(
+            PyObject_CallMethod(this->ptr(), const_cast<char*>("view"), const_cast<char*>("O"), dt.ptr())
+        )
+    );
+}
+
+ndarray ndarray::copy() const {
+    return ndarray(
+        python::detail::new_reference(
+            PyObject_CallMethod(this->ptr(), const_cast<char*>("copy"), const_cast<char*>(""))
+        )
+    );
+}
+
+dtype ndarray::get_dtype() const {
+    return dtype(python::detail::borrowed_reference(get_struct()->descr));
+}
+
+object ndarray::get_base() const {
+    if (get_struct()->base == NULL) return object();
+    return object(python::detail::borrowed_reference(get_struct()->base));
+}
+
+void ndarray::set_base(object const & base) {
+    Py_XDECREF(get_struct()->base);
+    if (base != object()) {
+        Py_INCREF(base.ptr());
+        get_struct()->base = base.ptr();
+    } else {
+        get_struct()->base = NULL;
+    }
+}
+
+ndarray::bitflag const ndarray::get_flags() const {
+    return numpy::detail::numpy_to_bitflag(get_struct()->flags);
+}
+
+ndarray ndarray::transpose() const {
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_Transpose(reinterpret_cast<PyArrayObject*>(this->ptr()), NULL)
+        )
+    );
+}
+
+ndarray ndarray::squeeze() const {
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_Squeeze(reinterpret_cast<PyArrayObject*>(this->ptr()))
+        )
+    );
+}
+
+object ndarray::scalarize() const {
+    Py_INCREF(ptr());
+    return object(python::detail::new_reference(PyArray_Return(reinterpret_cast<PyArrayObject*>(ptr()))));
+}
+
+ndarray zeros(tuple const & shape, dtype const & dt) {
+    int nd = len(shape);
+    Py_intptr_t dims[nd];
+    for (int n=0; n<nd; ++n) dims[n] = extract<Py_intptr_t>(shape[n]);
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_Zeros(nd, dims, detail::incref_dtype(dt), 0)
+        )
+    );
+}
+
+ndarray zeros(int nd, Py_intptr_t const * shape, dtype const & dt) {
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_Zeros(nd, const_cast<Py_intptr_t*>(shape), detail::incref_dtype(dt), 0)
+        )
+    );
+}
+
+ndarray empty(tuple const & shape, dtype const & dt) {
+    int nd = len(shape);
+    Py_intptr_t dims[nd];
+    for (int n=0; n<nd; ++n) dims[n] = extract<Py_intptr_t>(shape[n]);
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_Empty(nd, dims, detail::incref_dtype(dt), 0)
+        )
+    );    
+}
+
+ndarray empty(int nd, Py_intptr_t const * shape, dtype const & dt) {
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_Empty(nd, const_cast<Py_intptr_t*>(shape), detail::incref_dtype(dt), 0)
+        )
+    );
+}
+
+ndarray array(object const & obj) {
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_FromAny(obj.ptr(), NULL, 0, 0, NPY_ENSUREARRAY, NULL)
+        )
+    );
+}
+
+ndarray array(object const & obj, dtype const & dt) {
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_FromAny(obj.ptr(), detail::incref_dtype(dt), 0, 0, NPY_ENSUREARRAY, NULL)
+        )
+    );
+}
+
+ndarray from_object(object const & obj, dtype const & dt, int nd_min, int nd_max, ndarray::bitflag flags) {
+    int requirements = detail::bitflag_to_numpy(flags);
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_FromAny(
+                obj.ptr(),
+                detail::incref_dtype(dt),
+                nd_min, nd_max,
+                requirements,
+                NULL
+            )
+        )
+    );
+}
+
+ndarray from_object(object const & obj, int nd_min, int nd_max, ndarray::bitflag flags) {
+    int requirements = detail::bitflag_to_numpy(flags);
+    return ndarray(
+        python::detail::new_reference(
+            PyArray_FromAny(
+                obj.ptr(),
+                NULL,
+                nd_min, nd_max,
+                requirements,
+                NULL
+            )
+        )
+    );
+}
+
+}}}

libs/python/numpy/src/numpy.cpp

@@ -0,0 +1,13 @@
+#define BOOST_PYTHON_NUMPY_INTERNAL_MAIN
+#include <boost/python/numpy/internal.hpp>
+
+namespace boost { namespace python {
+
+namespace numpy {
+
+void initialize() {
+    import_array();
+    import_ufunc();
+}
+
+}}}

libs/python/numpy/src/scalars.cpp

@@ -0,0 +1,35 @@
+#define BOOST_PYTHON_NUMPY_INTERNAL
+#include <boost/python/numpy/internal.hpp>
+
+namespace boost { namespace python {
+namespace converter {
+NUMPY_OBJECT_MANAGER_TRAITS_IMPL(PyVoidArrType_Type, python::numpy::void_)
+} // namespace boost::python::converter
+
+namespace numpy {
+
+void_::void_(Py_ssize_t size) : 
+    object(
+        python::detail::new_reference(
+            PyObject_CallFunction((PyObject*)&PyVoidArrType_Type, const_cast<char*>("i"), size)
+        )
+    )
+{}
+
+void_ void_::view(dtype const & dt) const {
+    return void_(
+        python::detail::new_reference(
+            PyObject_CallMethod(this->ptr(), const_cast<char*>("view"), const_cast<char*>("O"), dt.ptr())
+        )
+    );
+}
+
+void_ void_::copy() const {
+    return void_(
+        python::detail::new_reference(
+            PyObject_CallMethod(this->ptr(), const_cast<char*>("copy"), const_cast<char*>(""))
+        )
+    );
+}
+
+}}}

libs/python/numpy/src/ufunc.cpp

@@ -0,0 +1,48 @@
+#define BOOST_PYTHON_NUMPY_INTERNAL
+#include <boost/python/numpy/internal.hpp>
+#include <boost/python/numpy/ufunc.hpp>
+
+namespace boost { namespace python {
+namespace converter {
+NUMPY_OBJECT_MANAGER_TRAITS_IMPL(PyArrayMultiIter_Type, python::numpy::multi_iter)
+} // namespace boost::python::converter
+
+namespace numpy {
+
+multi_iter make_multi_iter(object const & a1) {
+    return multi_iter(python::detail::new_reference(PyArray_MultiIterNew(1, a1.ptr())));
+}
+
+multi_iter make_multi_iter(object const & a1, object const & a2) {
+    return multi_iter(python::detail::new_reference(PyArray_MultiIterNew(2, a1.ptr(), a2.ptr())));
+}
+
+multi_iter make_multi_iter(object const & a1, object const & a2, object const & a3) {
+    return multi_iter(python::detail::new_reference(PyArray_MultiIterNew(3, a1.ptr(), a2.ptr(), a3.ptr())));
+}
+
+void multi_iter::next() {
+    PyArray_MultiIter_NEXT(ptr());
+}
+
+bool multi_iter::not_done() const {
+    return PyArray_MultiIter_NOTDONE(ptr());
+}
+
+char * multi_iter::get_data(int i) const {
+    return reinterpret_cast<char*>(PyArray_MultiIter_DATA(ptr(), i));
+}
+
+int const multi_iter::get_nd() const {
+    return reinterpret_cast<PyArrayMultiIterObject*>(ptr())->nd;
+}
+
+Py_intptr_t const * multi_iter::get_shape() const {
+    return reinterpret_cast<PyArrayMultiIterObject*>(ptr())->dimensions;
+}
+
+Py_intptr_t const multi_iter::shape(int n) const {
+    return reinterpret_cast<PyArrayMultiIterObject*>(ptr())->dimensions[n];
+}
+
+}}}