numpy_array.hpp

Go to the documentation of this file.

// mds_utils/python/ublas/numpy_array.hpp
//
// Copyright (c) 2014 - Michele De Stefano (micdestefano@users.sourceforge.net)
//
// Distributed under the MIT License (See accompanying file COPYING)

#ifndef MDS_UTILS_PYTHON_UBLAS_NUMPY_ARRAY_HPP_INCLUDED
#define MDS_UTILS_PYTHON_UBLAS_NUMPY_ARRAY_HPP_INCLUDED

#include <numpy/arrayobject.h>
#include <mds_utils/python/obj.hpp>
#include <mds_utils/python/numpy/array_iterator.hpp>
#include <boost/numeric/ublas/traits.hpp>
#include <boost/numeric/ublas/detail/vector_assign.hpp>


namespace boost { namespace numeric { namespace ublas {

template<class T>
class NumPy1DArray : public mds_utils::python::Obj, public vector_container< NumPy1DArray<T> > {

    typedef NumPy1DArray    self_type;

public:

    typedef npy_intp                                                    size_type;
    typedef npy_intp                                                    difference_type;
    typedef T                                                           value_type;
    typedef const T&                                                    const_reference;
    typedef T&                                                          reference;
    typedef T*                                                          pointer;
    typedef const T*                                                    const_pointer;
    typedef pointer                                                     array_type;
    typedef const vector_reference<const self_type>                     const_closure_type;
    typedef vector_reference<self_type>                                 closure_type;
    typedef self_type                                                   vector_temporary_type;
    typedef dense_tag                                                   storage_category;



    PyArrayObject* getPyArrayObject() const {
        return reinterpret_cast<PyArrayObject*>(getPyObject());
    }

    operator PyArrayObject*() const {
        return getPyArrayObject();
    }


    void incref() {
        PyArray_INCREF(reinterpret_cast<PyArrayObject*>(getPyObject()));
        ++decref_on_destroy;
    }

    void decref() {
        if (decref_on_destroy > 0) {
            PyArray_XDECREF(reinterpret_cast<PyArrayObject*>(getPyObject()));
            --decref_on_destroy;
        }
    }

    PyObject* transfer() {
        PyObject *pret(getPyObject());
        PyArray_INCREF(reinterpret_cast<PyArrayObject*>(pret));
        reset();
        return pret;
    }


    // Construction and destruction

    BOOST_UBLAS_INLINE
    NumPy1DArray() {}


    explicit BOOST_UBLAS_INLINE
    NumPy1DArray(size_type n) :
        mds_utils::python::Obj(PyArray_SimpleNew(1,&n,
                mds_utils::python::numpy::numpy_dtype_traits<T>::typenum)) {

        get_ownership();
    }


    template<class FwIt>
    BOOST_UBLAS_INLINE
    NumPy1DArray(FwIt b,FwIt e) : NumPy1DArray(std::distance(b,e)) {

        using namespace mds_utils::python::numpy;

        if (b != e) std::copy(b,e,NDArrayIterator<T,c_storage,
                NPY_ITER_WRITEONLY>(getPyArrayObject()));
    }


    explicit BOOST_UBLAS_INLINE
    NumPy1DArray(size_type n,const_reference init) : NumPy1DArray(n) {

        using namespace mds_utils::python::numpy;

        NDArrayIterator<T,c_storage,NPY_ITER_WRITEONLY>
            b(getPyArrayObject()),
            e(b,true);

        if (b != e) std::fill(b,e,init);
    }


    BOOST_UBLAS_INLINE
    size_type size() const {
        PyArrayObject   *po(getPyArrayObject());
        if (po == NULL) return 0;
        return PyArray_SIZE(po);
    }

    BOOST_UBLAS_INLINE
    NumPy1DArray(const NumPy1DArray &rhs) : NumPy1DArray(rhs.size()) {

        if (PyArray_CopyInto(getPyArrayObject(),rhs.getPyArrayObject())) {
            throw std::runtime_error("Cannot copy the NumPy array.");
        }
    }

    BOOST_UBLAS_INLINE
    NumPy1DArray(NumPy1DArray&& rhs) : mds_utils::python::Obj(std::move(rhs)) {}

    template<class AE>
    BOOST_UBLAS_INLINE
    NumPy1DArray(const vector_expression<AE> &ae) : NumPy1DArray(ae().size()) {

         vector_assign<scalar_assign>(*this,ae);
    }


    explicit NumPy1DArray(PyObject *rhs) : mds_utils::python::Obj(rhs) {

        if (!PyArray_Check(rhs)) {
            throw std::invalid_argument("Cannot convert to "
                "boost::numeric::ublas::NumPy1DArrayr: rhs is not a NumPy array.");
        }

        incref();
    }


private:

    void check_index(size_type i) const {
        assert(getPyObject() != NULL && i >= 0 && i < size());
    }

public:

    // -----------------
    // Storage accessors
    // -----------------


    BOOST_UBLAS_INLINE
    pointer find_element(size_type i) {
        return const_cast<pointer>(const_cast<const self_type&>(*this).find_element(i));
    }


    BOOST_UBLAS_INLINE
    const_pointer find_element(size_type i) const {
        check_index(i);
        return const_cast<const pointer>(reinterpret_cast<pointer>(PyArray_GETPTR1(getPyArrayObject(),i)));
    }


    // --------------
        // Element access
    // --------------

    BOOST_UBLAS_INLINE
    const_reference operator () (size_type i) const {
        check_index(i);
        return *reinterpret_cast<pointer>(PyArray_GETPTR1(getPyArrayObject(),i));
    }


    BOOST_UBLAS_INLINE
    reference operator () (size_type i) {
        check_index(i);
        return *reinterpret_cast<pointer>(PyArray_GETPTR1(getPyArrayObject(),i));
    }


    BOOST_UBLAS_INLINE
    const_reference operator [] (size_type i) const {
        return (*this)(i);
    }


    BOOST_UBLAS_INLINE
    reference operator [] (size_type i) {
        return (*this) (i);
    }

    // -------
    // Zeroing
    // -------

    BOOST_UBLAS_INLINE
    void clear() {

        using namespace mds_utils::python::numpy;

        NDArrayIterator<T,c_storage,NPY_ITER_WRITEONLY>
            b(getPyArrayObject()),
            e(b,true);

        if (b != e) std::fill(b,e,static_cast<T>(0));
    }

    // Assignment

    template<class AE>
    BOOST_UBLAS_INLINE
    NumPy1DArray& assign(const vector_expression<AE> &ae) {
        vector_assign<scalar_assign>(*this, ae);
        return *this;
    }

    BOOST_UBLAS_INLINE
    NumPy1DArray& assign_temporary(NumPy1DArray& rhs) {
        *this = rhs.getPyObject();
        incref();
        return *this;
    }


    BOOST_UBLAS_INLINE
    NumPy1DArray& operator =(const NumPy1DArray& rhs) {

        size_type   n(rhs.size());
        
        if (size() != n) {
            reset();

            m_po = PyArray_SimpleNew(1,&n,
                    mds_utils::python::numpy::numpy_dtype_traits<T>::typenum);

            get_ownership();
        }

        if (PyArray_CopyInto(getPyArrayObject(),rhs.getPyArrayObject())) {
            throw std::runtime_error("Cannot copy the NumPy array.");
        }

        return *this;
    }


    BOOST_UBLAS_INLINE
    NumPy1DArray& operator =(NumPy1DArray&& rhs) {

        mds_utils::python::Obj::operator =(std::move(rhs));
        
        return *this;
    }


    NumPy1DArray& operator =(PyObject *rhs) {

        if (!PyArray_Check(rhs) && getPyObject() != rhs) {
            throw std::invalid_argument("Cannot assign: rhs is not a NumPy array.");
        }

        mds_utils::python::Obj::operator =(rhs);

        return *this;
    }


    template<class val_T>
    NumPy1DArray& operator =(const val_T& val) {

        using namespace std;

        if (size() == 0) return *this;

        for_each(begin(),end(),[&val](T& x){ x = val; });

        return *this;
    }


    template<class C>
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator =(const vector_container<C>& rhs) {
        assert(size() == rhs.size());

        assign(rhs);
        return  *this;
    }


    template<class AE>
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator =(const vector_expression<AE> &ae) {
        self_type   temporary(ae);
        return assign_temporary(temporary);
    }


    // -------------------
    // Computed assignment
    // -------------------

    template<class AE>
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator += (const vector_expression<AE> &ae) {
        self_type temporary(*this + ae);
        return assign_temporary (temporary);
    }

    template<class C>          // Container assignment without temporary
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator += (const vector_container<C>& v) {
        plus_assign (v);
        return *this;
    }


    template<class AE>
    BOOST_UBLAS_INLINE
    NumPy1DArray& plus_assign(const vector_expression<AE>& ae) {
        vector_assign<scalar_plus_assign> (*this, ae);
        return *this;
    }


    template<class AE>
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator -= (const vector_expression<AE>& ae) {
        self_type temporary (*this - ae);
        return assign_temporary (temporary);
    }


    template<class C>          // Container assignment without temporary
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator -= (const vector_container<C>& v) {
        minus_assign (v);
        return *this;
    }


    template<class AE>
    BOOST_UBLAS_INLINE
    NumPy1DArray& minus_assign (const vector_expression<AE>& ae) {
        vector_assign<scalar_minus_assign> (*this, ae);
        return *this;
    }


    template<class AT>
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator *= (const AT &at) {
        vector_assign_scalar<scalar_multiplies_assign> (*this, at);
        return *this;
    }


    template<class AT>
    BOOST_UBLAS_INLINE
    NumPy1DArray& operator /= (const AT &at) {
        vector_assign_scalar<scalar_divides_assign> (*this, at);
        return *this;
    }


    // --------
    // Swapping
    // --------

    BOOST_UBLAS_INLINE
    void swap(NumPy1DArray& rhs) {
        using namespace mds_utils::python::numpy;

        if (this != &rhs) {
            NDArrayIterator<T,c_storage,NPY_ITER_READWRITE>
                it(getPyArrayObject()),
                e(it,true),
                it_rhs(rhs.getPyArrayObject());

            T   aux;

            for (;it != e;++it,++it_rhs) {
                aux = *it;
                *it = *it_rhs;
                *it_rhs = aux;
            }
        }
    }

    BOOST_UBLAS_INLINE
    friend void swap(NumPy1DArray& v1,NumPy1DArray& v2) {
        v1.swap(v2);
    }

private:
    // Use the storage array iterator
    typedef mds_utils::python::numpy::NDArrayIterator<T,
            mds_utils::python::numpy::c_storage,NPY_ITER_READONLY>
                const_subiterator_type;

    typedef mds_utils::python::numpy::NDArrayIterator<T,
            mds_utils::python::numpy::c_storage,NPY_ITER_READWRITE>
                subiterator_type;

public:


    typedef indexed_iterator<self_type,dense_random_access_iterator_tag>            iterator;
    typedef indexed_const_iterator<self_type, dense_random_access_iterator_tag>     const_iterator;

    // --------------
    // Element lookup
    // --------------

    BOOST_UBLAS_INLINE
    const_iterator find(size_type i) const {
        return const_iterator(*this,i);
    }

    BOOST_UBLAS_INLINE
    iterator find(size_type i) {
        return iterator(*this,i);
    }

    BOOST_UBLAS_INLINE
    const_iterator begin() const {
        return find(0);
    }

    BOOST_UBLAS_INLINE
    const_iterator end() const {
        return find(size());
    }

    BOOST_UBLAS_INLINE
    iterator begin () {
        return find (0);
    }

    BOOST_UBLAS_INLINE
    iterator end () {
        return find (size());
    }

    // Reverse iterator
    typedef reverse_iterator_base<const_iterator>   const_reverse_iterator;
    typedef reverse_iterator_base<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 ());
    }

};


}}}

#endif /* MDS_UTILS_PYTHON_UBLAS_NUMPY_ARRAY_HPP_INCLUDED */