matrix_sparse.hpp

Go to the documentation of this file.

// matrix_sparse.hpp
// 
// Copyright (c) 2009 - Michele De Stefano (micdestefano@users.sourceforge.net)
// 
// Distributed under the MIT License (See accompanying file LICENSE)

#ifndef MDS_UTILS_PYTHON_UBLAS_MATRIX_SPARSE_HPP_INCLUDED
#define MDS_UTILS_PYTHON_UBLAS_MATRIX_SPARSE_HPP_INCLUDED

#include <mds_utils/python/ublas/detail/common.hpp>
#include <mds_utils/ublas/traits.hpp>
#include <mds_utils/python/tuple.hpp>
#include <mds_utils/python/numpy/array_iterator.hpp>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/fusion/container/vector.hpp>
#include <boost/numeric/conversion/cast.hpp>
#include <algorithm>
#include <string>

namespace mds_utils { namespace python { namespace ublas {

    namespace detail {
    
template<class orientation_category>
struct scipy_sparse_mtx;

template<>
struct scipy_sparse_mtx<boost::numeric::ublas::row_major_tag> {
    
    static const std::string    name;
};

const std::string scipy_sparse_mtx<boost::numeric::ublas::row_major_tag>::name("csr_matrix");
    
template<>
struct scipy_sparse_mtx<boost::numeric::ublas::column_major_tag> {
    
    static const std::string    name;
};

const std::string scipy_sparse_mtx<boost::numeric::ublas::column_major_tag>::name("csc_matrix");

    } // detail




template<class M>
inline typename std::enable_if<
    mds_utils::ublas::is_compressed_matrix<M>::value,M>::type
        get(PyObject *po) {

    using namespace mds_utils::python::numpy;
        
    Obj module(PyImport_ImportModule("scipy.sparse"));
    
    module.get_ownership();

    Obj type_obj(module.attr(detail::scipy_sparse_mtx<
            typename M::orientation_category>::name));
    

    if (!PyType_Check(type_obj) ||
        !PyObject_TypeCheck(po,
                reinterpret_cast<PyTypeObject*>(static_cast<PyObject*>(type_obj)))) {
        
        throw std::invalid_argument(
                "Cannot get a uBLAS compressed matrix from the given parameter.\n"
                "Check the data type and the storage layout (i.e. csr or csc).");
    }
    
    Obj     o(po);

    Tuple   shape(o.attr("shape"));

    Obj     data(o.attr("data")),
            indices(o.attr("indices")),
            indptr(o.attr("indptr"));

    size_t
        Nr(get<size_t>(shape[0])),
        Nc(get<size_t>(shape[1])),
        nnz(get<size_t>(o.attr("nnz")));

    if (Nr == 0 || Nc == 0) return M();

    M   mout(Nr,Nc,nnz);

    NDArrayIterator<typename M::value_type,c_storage,NPY_ITER_READONLY>
        it_data(reinterpret_cast<PyArrayObject*>(static_cast<PyObject*>(data))),
        it_data_end(it_data,true);
    NDArrayIterator<typename M::size_type,c_storage,NPY_ITER_READONLY>
        it_inds(reinterpret_cast<PyArrayObject*>(static_cast<PyObject*>(indices))),
        it_inds_end(it_inds,true),
        it_indptr(reinterpret_cast<PyArrayObject*>(static_cast<PyObject*>(indptr))),
        it_indptr_end(it_indptr,true);

    std::copy(it_data,it_data_end,mout.value_data().begin());
    std::copy(it_inds,it_inds_end,mout.index2_data().begin());
    std::copy(it_indptr,it_indptr_end,mout.index1_data().begin());

    mout.set_filled(len(indptr),len(data));
    
    return mout;
}



template<class ublas_mat_T>
inline typename std::enable_if<
        mds_utils::ublas::is_compressed_matrix<ublas_mat_T>::value,
        PyObject*
    >::type to_python(const ublas_mat_T& mat) {
    
    using namespace mds_utils::python::numpy;
    namespace fus = boost::fusion;

    Obj module(PyImport_ImportModule("scipy.sparse"));

    module.get_ownership();

    Obj type_obj(module.attr(detail::scipy_sparse_mtx<
            typename ublas_mat_T::orientation_category>::name));

    if (!PyType_Check(type_obj)) {

        throw std::runtime_error(
                "Cannot create a " + detail::scipy_sparse_mtx<
                    typename ublas_mat_T::orientation_category>::name + " object.");
    }

    int
        flags(detail::numpy_flags<
            typename ublas_mat_T::orientation_category>());
        
    npy_intp
        data_dims[] = {boost::numeric_cast<npy_intp>(mat.filled2())},
        indptr_dims[] = {boost::numeric_cast<npy_intp>(mat.filled1())},
        mat_shape[] = {
            boost::numeric_cast<npy_intp>(mat.size1()),
            boost::numeric_cast<npy_intp>(mat.size2())
        };

    PyObject
        *pdata(PyArray_New(&PyArray_Type,1,data_dims,
            numpy_dtype_traits<
                typename ublas_mat_T::value_type>::typenum,
            NULL,NULL,0,flags,NULL)),
        *pindices(PyArray_New(&PyArray_Type,1,data_dims,
            numpy_dtype_traits<int>::typenum,
            NULL,NULL,0,flags,NULL)),
        *pindptr(PyArray_New(&PyArray_Type,1,indptr_dims,
            numpy_dtype_traits<int>::typenum,
            NULL,NULL,0,flags,NULL));

    if (pdata == NULL) {
        throw std::runtime_error("Could not create the data array.");
    }
    if (pindices == NULL) {
        throw std::runtime_error("Could not create the indices array.");
    }
    if (pindptr == NULL) {
        throw std::runtime_error("Could not create the indptr array.");
    }

    NDArrayIterator<typename ublas_mat_T::value_type,c_storage,NPY_ITER_WRITEONLY>
        it_data(reinterpret_cast<PyArrayObject*>(pdata));

    NDArrayIterator<int,c_storage,NPY_ITER_WRITEONLY>
        it_indices(reinterpret_cast<PyArrayObject*>(pindices)),
        it_indptr(reinterpret_cast<PyArrayObject*>(pindptr));

    size_t  k;
    typename ublas_mat_T::index_array_type::const_iterator
        mat_it_indices(mat.index2_data().begin()),
        mat_it_indptr(mat.index1_data().begin());
    typename ublas_mat_T::value_array_type::const_iterator
        mat_it_data(mat.value_data().begin());

    for (k=0;k < *indptr_dims;++k,++mat_it_indptr,++it_indptr) {
        *it_indptr = *mat_it_indptr;
    }
    for (k=0;k < *data_dims;++k,++mat_it_indices,++mat_it_data,++it_indices,++it_data) {
        *it_indices = *mat_it_indices;
        *it_data = *mat_it_data;
    }


    fus::vector<PyObject*,PyObject*,PyObject*>
        mstruct(pdata,pindices,pindptr);

    Tuple
        tup_mstruct,tup_sizes,tup_args;

    tup_mstruct.set(mstruct);
    tup_sizes.set(mat_shape,mat_shape+2);
    
    fus::vector<Obj,Obj>
        args(tup_mstruct,tup_sizes);
    
    tup_args.set(args);

    Obj out_mtx(type_obj(tup_args));
    
    return out_mtx.transfer();
}

}}}

#endif