Ravi wrote:
> On Wednesday 21 January 2009 10:22:36 Neal Becker wrote:
>> > [http://mail.python.org/pipermail/cplusplus-sig/2008-
October/013825.html
>>
>> Thanks for reminding me about this!
>>
>> Do you have a current version of the code? I grabbed the files from the
>> above message, but I see some additional subsequent messages with more
>> patches.
>
> That is the latest publicly posted code. Since then, there is just one
> minor patch (attached) which enables use of row-major (c-contiguous)
> arrays.
>
> This does *not* work with strided arrays which would be a fair bit of
> effort to support. Further, you will have to work with the numpy iterator
> interface, which, while well-designed, is a great illustration of the
> effort required to support OO programming in an non-OO language, and is
> pretty tedious to map to the ublas storage iterator interface. If you do
> implement it, I would very much like to take a look at it.
>
> Regards,
> Ravi
I'm only interested in simple strided 1-d vectors. In that case, I think
your code already works. If you have c++ code using the iterator interface,
the iterators dereference will use (*array )[index]. This will use
operator[], which will call PyArray_GETPTR. So I think this will obey
strides.
Unfortunately, it will also be slow. I suggest something like the enclosed.
I have done some simple tests, and it seems to work.
// Copyright Ravikiran Rajagopal 2003.
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef NUMPY_RR2003_H
#define NUMPY_RR2003_H
// Meta-programming magic to interface ublas and numpy
#include <boost/iterator/iterator_facade.hpp>
#include <boost/mpl/at.hpp>
#include <boost/mpl/if.hpp>
#include <boost/mpl/int.hpp>
#include <boost/mpl/map.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/python/refcount.hpp>
#include <boost/type_traits/add_const.hpp>
#include <boost/type_traits/add_reference.hpp>
#include <boost/type_traits/is_const.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/utility/enable_if.hpp>
#include <algorithm>
#include <complex>
#include <exception>
#include <memory>
#ifndef IN_FILE_IMPORTING_ARRAY
#define NO_IMPORT_ARRAY
#endif
#define PY_ARRAY_UNIQUE_SYMBOL Wicket_Numpy
#include <numpy/arrayobject.h>
namespace numpy {
// Create a type map for numpy supported types
typedef boost::mpl::map<
boost::mpl::pair< short, boost::mpl::int_<PyArray_SHORT> >
, boost::mpl::pair< int, boost::mpl::int_<PyArray_INT> >
, boost::mpl::pair< long, boost::mpl::int_<PyArray_LONG> >
, boost::mpl::pair< long long, boost::mpl::int_<PyArray_LONGLONG> >
, boost::mpl::pair< unsigned short, boost::mpl::int_<PyArray_USHORT> >
, boost::mpl::pair< unsigned int, boost::mpl::int_<PyArray_UINT> >
, boost::mpl::pair< unsigned long, boost::mpl::int_<PyArray_ULONG> >
, boost::mpl::pair< unsigned long long, boost::mpl::int_<PyArray_ULONGLONG> >
, boost::mpl::pair< float, boost::mpl::int_<PyArray_FLOAT> >
, boost::mpl::pair< double, boost::mpl::int_<PyArray_DOUBLE> >
, boost::mpl::pair< long double, boost::mpl::int_<PyArray_LONGDOUBLE> >
, boost::mpl::pair< std::complex<float>, boost::mpl::int_<PyArray_CFLOAT> >
, boost::mpl::pair< std::complex<double>, boost::mpl::int_<PyArray_CDOUBLE> >
, boost::mpl::pair< std::complex<long double>, boost::mpl::int_<PyArray_CLONGDOUBLE> >
> numpy_typemap;
namespace detail {
// Storage shared with a numpy array
template<class T>
class numpy_storage_array:
public boost::numeric::ublas::storage_array<numpy_storage_array<T> >
{
typedef numpy_storage_array<T> self_type;
template < class Value >
class simple_iter
: public boost::iterator_facade< simple_iter<Value>,
Value,
boost::random_access_traversal_tag,
typename boost::add_reference<Value>::type,
npy_intp >
{
private:
struct enabler {};
typedef typename boost::mpl::if_<typename boost::is_const<Value>::type,
const numpy_storage_array<T>,
numpy_storage_array<T>
>::type array_type;
public:
typedef boost::iterator_facade< simple_iter<Value>,
Value,
boost::random_access_traversal_tag,
typename boost::add_reference<Value>::type,
npy_intp > self_t;
typedef typename self_t::pointer pointer_t;
explicit simple_iter( array_type *a )
: pt (m_stride >= 0 ? (pointer_t)(PyArray_DATA (a->array_)) : ((pointer_t)(PyArray_DATA (a->array_ + a->size())))),
m_stride (a->stride()) {}
// template <class OtherValue>
// simple_iter( simple_iter<OtherValue> const &other,
// typename boost::enable_if<
// boost::is_convertible< OtherValue*, Value* >, enabler
// >::type = enabler() )
// : pt (&array( other.array ), m_stride (other.array->stride()), index( other.index ) {}
private:
friend class boost::iterator_core_access;
template <class> friend class simple_iter;
void increment() { pt += m_stride; }
void decrement() { pt -= m_stride; }
template <class OtherValue>
bool equal( simple_iter<OtherValue> const &other ) const
{
return pt == other.pt;
}
// Value &dereference() const { return ( *array )[index]; }
Value & dereference() const { return *pt; }
void advance( npy_intp n ) { pt += n*m_stride; }
template <class OtherValue>
npy_intp distance_to( simple_iter<OtherValue> const &other ) const
{
return (other.pt - pt)/m_stride;
}
private:
pointer_t pt;
// npy_intp index;
npy_intp m_stride;
};
public:
typedef typename std::allocator<T>::size_type size_type;
typedef typename std::allocator<T>::difference_type difference_type;
typedef T value_type;
typedef const T &const_reference;
typedef T &reference;
typedef const T *const_pointer;
typedef T *pointer;
typedef simple_iter<typename boost::add_const<T>::type> const_iterator;
typedef simple_iter<T> iterator;
private:
// No default constructor; must initialize object using a numpy array!
explicit BOOST_UBLAS_INLINE numpy_storage_array();
public:
// Construction and destruction
explicit BOOST_UBLAS_INLINE numpy_storage_array( PyObject *arr )
: boost::numeric::ublas::storage_array<numpy_storage_array<T> >()
{
if ( !PyArray_Check( arr )
|| ( PyArray_TYPE( arr ) != boost::mpl::at<numpy_typemap, T>::type::value )
|| ( PyArray_NDIM( arr ) != 1 )
)
throw std::bad_alloc();
array_ = ( PyArrayObject * )arr;
// Get a reference to the right sort of array
#ifdef NUMPY_STORAGE_ARRAY_MUST_BE_CONTIGUOUS
array_ = PyArray_GETCONTIGUOUS( array_ );
#else
Py_INCREF( array_ );
#endif
}
BOOST_UBLAS_INLINE numpy_storage_array(const numpy_storage_array &c)
: boost::numeric::ublas::storage_array<numpy_storage_array<T> >(),
array_( c.array_ )
{
Py_INCREF( array_ );
}
BOOST_UBLAS_INLINE ~numpy_storage_array()
{
Py_DECREF( array_ );
}
BOOST_UBLAS_INLINE void resize( size_type size_a )
{
if ( size_a != size() )
{
npy_intp newsize = size_a;
PyArray_Dims newshape = { 1, &newsize };
PyObject *dummy = PyArray_Resize( array_, &newshape, 1, NPY_ANYORDER );
if ( dummy == 0 )
throw std::bad_alloc();
Py_DECREF( dummy );
// Note that all iterators are invalid at this point, i.e., iterators of other
// arrays which refer to this underlying storage could also be invalid!
}
}
// Random Access Container
BOOST_UBLAS_INLINE size_type size() const { return PyArray_SIZE( array_ ); }
BOOST_UBLAS_INLINE size_type max_size() const { return size(); }
BOOST_UBLAS_INLINE bool empty() const { return size() == 0; }
// Element access
BOOST_UBLAS_INLINE const_reference operator[](size_type i) const
{
BOOST_UBLAS_CHECK( i < size(), boost::numeric::ublas::bad_index() );
return *( static_cast<T*>( PyArray_GETPTR1( array_, i ) ) );
}
BOOST_UBLAS_INLINE reference operator[](size_type i)
{
BOOST_UBLAS_CHECK( i < size(), boost::numeric::ublas::bad_index() );
return *( static_cast<T*>( PyArray_GETPTR1( array_, i ) ) );
}
// Assignment
BOOST_UBLAS_INLINE numpy_storage_array &operator=( const numpy_storage_array &a )
{
if ( this != &a )
{
resize( a.size() );
std::copy( a.begin(), a.end(), begin() );
}
return *this;
}
BOOST_UBLAS_INLINE numpy_storage_array &assign_temporary( numpy_storage_array &a )
{
swap( a );
return *this;
}
// Swapping
BOOST_UBLAS_INLINE void swap( numpy_storage_array &a )
{
if ( this != &a )
{
std::swap( array_, a.array_ );
}
}
BOOST_UBLAS_INLINE friend void swap( numpy_storage_array &a1, numpy_storage_array &a2 )
{
a1.swap( a2 );
}
BOOST_UBLAS_INLINE const_iterator begin() const
{
return const_iterator(this);
}
BOOST_UBLAS_INLINE const_iterator end() const
{
const_iterator it = begin();
it.advance (size());
return it;
}
BOOST_UBLAS_INLINE iterator begin()
{
return iterator( this );
}
BOOST_UBLAS_INLINE iterator end()
{
iterator it = begin();
it.advance (size());
return it;
}
// Reverse iterators
typedef std::reverse_iterator<const_iterator> const_reverse_iterator;
typedef std::reverse_iterator<iterator> reverse_iterator;
BOOST_UBLAS_INLINE const_reverse_iterator rbegin() const
{
return const_reverse_iterator( end() );
}
BOOST_UBLAS_INLINE const_reverse_iterator rend() const
{
return const_reverse_iterator( begin() );
}
BOOST_UBLAS_INLINE reverse_iterator rbegin()
{
return reverse_iterator( end() );
}
BOOST_UBLAS_INLINE reverse_iterator rend()
{
return reverse_iterator( begin() );
}
// Underlying array object
PyArrayObject *get_array() const { return array_; }
npy_intp stride () const { return PyArray_STRIDE (array_, 0)/npy_intp(sizeof(T)); }
#if 0
private:
friend class boost::serialization::access;
// Serialization
template<class Archive>
void serialize(Archive & ar, const unsigned int version)
{
serialization::collection_size_type s(size_);
ar & serialization::make_nvp("size",s);
if ( Archive::is_loading::value ) {
resize(s);
}
ar & serialization::make_array(data_, s);
}
#endif
private:
PyArrayObject *array_;
}; // numpy_storage_array
} // namespace detail
template <typename data_t>
struct array_from_py
{
typedef typename
boost::numeric::ublas::vector< data_t, detail::numpy_storage_array<data_t> > type;
};
template <typename data_t>
struct const_array_from_py
{
private:
typedef typename array_from_py<data_t>::type o_type;
public:
typedef typename boost::add_const<o_type>::type type;
};
} // namespace numpy
#endif
#define IN_FILE_IMPORTING_ARRAY
#include "numpy.hpp"
#include "numpyregister.hpp"
#include <boost/python/def.hpp>
#include <boost/python/extract.hpp>
#include <boost/python/module.hpp>
#include <boost/python/scope.hpp>
//#include <boost/python/return_value_policy.hpp>
//#include <boost/python/return_by_value.hpp>
#include "reference_existing_object_1.hpp"
//#include <iostream>
//#include <boost/shared_ptr.hpp>
#include <boost/python.hpp>
#include <complex>
#include <boost/range.hpp>
namespace ublas = boost::numeric::ublas;
namespace bp = boost::python;
typedef std::complex<double> complex_t;
template<typename out_t, typename in_t>
out_t sum (in_t const& in) {
out_t out (boost::size (in));
typedef typename boost::range_value<out_t>::type value_t;
typename boost::range_const_iterator<in_t>::type i = boost::begin (in);
typename boost::range_const_iterator<in_t>::type e = boost::end (in);
typename boost::range_iterator<out_t>::type o = boost::begin (out);
value_t s = 0;
for (; i != e; ++i, ++o) {
s += *i;
*o = s;
}
return out;
}
BOOST_PYTHON_MODULE (test1) {
import_array();
import_smart_ptr_deleter_type();
numpy::register_default_ublas_to_python( true );
numpy::register_ublas_from_python_converters();
bp::def ("sum", &sum<ublas::vector<complex_t>,numpy::array_from_py<complex_t>::type>);
}
import numpy
import test1
import sys
u = numpy.array (xrange (4), dtype=complex)
print sys.getrefcount (u)
v = test1.sum (u)
print sys.getrefcount (u), sys.getrefcount (v)
w = test1.sum (u[::2])
x = test1.sum (u[::-1])
y = test1.sum (u[::-2])
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
