Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 3ddb1815 authored by hhakim's avatar hhakim
Browse files

Implement sparse matio/matlab file output for Faust factors with support of...

Implement sparse matio/matlab file output for Faust factors with support of real and complex scalars, Faust's conjugate and transpose.

Main addon: Faust::Transform::MatSparse::toMatIOVar() (writing the sparse mat to a matio var of type mat_sparse_t). The old toMatIOVar() has been renamed toMatIOVarDense().
parent 5c9d2cf0
No related branches found
No related tags found
No related merge requests found
......@@ -11,19 +11,21 @@ import math
class TestFaustPy(unittest.TestCase):
MAX_NUM_FACTORS = 32 # for the tested Faust
MAX_NUM_FACTORS = 8 # for the tested Faust
MAX_DIM_SIZE = 512
MIN_DIM_SIZE = 3
def setUp(self):
""" Initializes the tests objects """
r = random.Random() # initialized from time or system
num_factors = r.randint(1, TestFaustPy.MAX_NUM_FACTORS)
factors = []
d2 = r.randint(1, TestFaustPy.MAX_DIM_SIZE)
d2 = r.randint(TestFaustPy.MIN_DIM_SIZE, TestFaustPy.MAX_DIM_SIZE)
for i in range(0, num_factors):
d1, d2 = d2, r.randint(1, TestFaustPy.MAX_DIM_SIZE)
factors += [sparse.random(d1, d2, density=0.1, format='csr',
dtype=np.float64).todense()]
print("factor",i,":", factors[i])
self.F = Faust(factors)
self.factors = factors
print("Tests on random Faust with dims=", self.F.get_nb_rows(),
......
......@@ -718,7 +718,7 @@ t_print_file.stop();
template<typename FPP>
matvar_t* Faust::MatDense<FPP, Cpu>::toMatIOVar(bool transpose, bool conjugate) const
{
matvar_t *var = NULL; //TODO: should be nullptr in C++11
matvar_t *var = NULL;
size_t dims[2];
int opt = typeid(mat(0,0))==typeid(complex<double>(1.0,1.0))?MAT_F_COMPLEX:0;
mat_complex_split_t z = {NULL,NULL};
......@@ -751,7 +751,8 @@ matvar_t* Faust::MatDense<FPP, Cpu>::toMatIOVar(bool transpose, bool conjugate)
{
dims[0] = this->getNbRow();
dims[1] = this->getNbCol();
if(!opt) //we don't copy the data again, we use it directly (col-major order organized)
if(!opt) // we use directly the data pointer (col-major order organized)
// but without the MAT_F_DONT_COPY_DATA flag, MatIO copies the data internally
var = Mat_VarCreate(NULL, MAT_C_DOUBLE, MAT_T_DOUBLE, 2, dims, (FPP*) mat.data(), opt);
else {
Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> dst_re(mat.rows(), mat.cols());
......
......@@ -252,6 +252,13 @@ namespace Faust
// M = (*this)' * M if opThis='T'
void multiply(Faust::MatDense<FPP,Cpu> & M, char opThis) const;
matvar_t* toMatIOVar(bool transpose, bool conjugate) const;
//! \brief Converts the Matrix to a matio variable, especially useful for writing into a file with libmatio. The variable is encoded in full matrix format (on the contrary to toMatVarIOVar() which keeps the sparse representation).
// \param transpose: set to true to obtain the matio variable for the transpose Matrix.
// \param conjugate: set it to true to obtain the matio variable for the conjugate Matrix.
// \return The matio variable matvar_t if it succeeded or NULL otherwise.
// \see Faust::Transform::save_mat_file()
// \see toMatIOVar()
matvar_t* toMatIOVarDense(bool transpose, bool conjugate) const;
FPP normL1() const;
FPP normL1(faust_unsigned_int&) const;
Faust::Vect<FPP,Cpu> get_col(faust_unsigned_int id) const;
......
......@@ -43,7 +43,7 @@
#include <iostream>
#include <fstream>
#include <iomanip>
#include <complex>
using namespace std;
template<typename FPP>
......@@ -569,16 +569,112 @@ void Faust::MatSparse<FPP,Cpu>::init_from_file(const char* filename)
}
template<typename FPP>
matvar_t* Faust::MatSparse<FPP, Cpu>::toMatIOVar(bool transpose, bool conjugate) const
matvar_t* Faust::MatSparse<FPP, Cpu>::toMatIOVarDense(bool transpose, bool conjugate) const
{
matvar_t *var = NULL; //TODO: should be nullptr in C++11
matvar_t *var = NULL;
Faust::MatDense<FPP,Cpu> dense_factor;
//TODO: shouldn't be processed as dense factor, we lose compression of sparse matrix here
// shouldn't be processed as dense factor, we lose compression of sparse matrix here
// (see toMatIOVar())
dense_factor = (*this); //data is copied with operator = redef.
var = dense_factor.toMatIOVar(transpose, conjugate);
return var;
}
template<typename FPP>
matvar_t* Faust::MatSparse<FPP, Cpu>::toMatIOVar(bool transpose, bool conjugate) const {
//TODO: refactor this function because it is a bit too long
matvar_t* var = NULL;
Eigen::SparseMatrix<FPP,Eigen::RowMajor> mat_;
size_t dims[2];
mat_sparse_t sparse = {0,};
// sparse.njc = (int) this->getNbCol()+1;
sparse.nzmax = (int) this->nnz;
sparse.ndata = (int) this->nnz;
int* jc;
int* ir = new int[sparse.nzmax];
double* data;
mat_complex_split_t z = {0,};
int nir = 0; //incremented later row by row
int i = 0;
int opt = typeid(getValuePtr()[0])==typeid(complex<double>(1.0,1.0))?MAT_F_COMPLEX:0;
if(opt) {
z.Re = new double[sparse.nzmax];
z.Im = new double[sparse.nzmax];
}
else
data = new double[sparse.nzmax];
if(transpose)
{
mat_ = mat.transpose();
dims[0]=this->getNbCol();
dims[1]= this->getNbRow();
}
else
{
mat_ = mat;
dims[0]=this->getNbRow();
dims[1]= this->getNbCol();
}
sparse.njc = dims[1]+1;
jc = new int[sparse.njc];
jc[sparse.njc-1] = this->nnz;
for(int j=0;j<sparse.njc-1;j++) jc[j] = -1;
// we use the transpose matrix because we are in row-major order but MatIO wants col-major order
// and the sparse matrix iterator respects the row-major order
Eigen::SparseMatrix<FPP,Eigen::RowMajor> st;
st = mat_.transpose();
for (int k=0; k<st.outerSize(); ++k)
for (typename Eigen::SparseMatrix<FPP,Eigen::RowMajor>::InnerIterator it(st,k); it; ++it)
{
// std::cout << "row:" << it.row() << " col:" << it.col() << " val:" << it.value() << std::endl;
if(it.row() > 0){
i=1;
while(it.row()>=i && jc[it.row()-i] < 0) {
jc[it.row()-i] = nir;
i++;
}
}
if(jc[it.row()] < 0) jc[it.row()] = nir;
ir[nir] = it.col();
if(opt) {
((double*)z.Re)[nir] = std::real((complex<double>)it.value());
if(conjugate)
((double*)z.Im)[nir] = -std::imag((complex<double>)it.value());
else
((double*)z.Im)[nir] = std::imag((complex<double>)it.value());
}
else
data[nir] = std::real(complex<double>(it.value()));
nir++;
}
i=1;
while(i<=st.rows() && jc[st.rows()-i] < 0) {
jc[st.rows()-i] = nir;
i++;
}
sparse.ir = ir;
sparse.jc = jc;
sparse.nir = nir;
if(opt) {
sparse.data = &z;
}
else
sparse.data = data;
var = Mat_VarCreate(NULL /* no-name */, MAT_C_SPARSE, MAT_T_DOUBLE, 2, dims, &sparse, opt);
// if(var != NULL)
// Mat_VarPrint(var,1);
delete[] jc;
delete[] ir;
delete[] data;
return var;
}
template<typename FPP>
FPP Faust::MatSparse<FPP, Cpu>::normL1(faust_unsigned_int& col_id) const
{
......
......@@ -144,7 +144,7 @@ namespace Faust
//! \brief Converts the Matrix to a matio variable, especially useful for writing into a file with libmatio.
// \param transpose: set to true to obtain the matio variable for the transpose Matrix.
// \param conjugate: set it to true to obtain the matio variable for the conjugate Matrix.
// \return The matio variable matvar_t if it succeeded or nullptr otherwise.
// \return The matio variable matvar_t if it succeeded or NULL otherwise.
// \see Faust::Transform::save_mat_file()
virtual matvar_t* toMatIOVar(bool transpose, bool conjugate) const=0;
......
......@@ -63,16 +63,17 @@ void Faust::MatGeneric<FPP,DEVICE>::setOp(const char op, faust_unsigned_int& nbR
//template <typename FPP, Device DEVICE>
template<typename FPP>
template<typename FPP>
Faust::MatGeneric<FPP,Cpu>* Faust::optimize(Faust::MatDense<FPP,Cpu> const & M,Faust::MatSparse<FPP,Cpu> const & S)
{
//std::cout<<"DEBUT OPTIMIZE "<<std::endl;
if ( (M.getNbCol() != S.getNbCol()) | (M.getNbRow() != S.getNbRow()) )
handleError("Faust::MatGeneric::", " Faust::optimize : matrix M and S have not the same size");
Faust::Vect<FPP,Cpu> x_dense(M.getNbCol());
for (int i=0;i<M.getNbCol();i++)
{
x_dense[i]=i*0.005;
......@@ -84,22 +85,19 @@ Faust::MatGeneric<FPP,Cpu>* Faust::optimize(Faust::MatDense<FPP,Cpu> const & M,F
int nb_mult=10;
Faust::Timer t_dense,t_sparse;
for (int i=0;i<nb_mult;i++)
{
{
x_sparse=x;
x_dense=x;
x_dense=x;
t_sparse.start();
S.multiply(x_sparse,'N');
S.multiply(x_sparse,'N');
t_sparse.stop();
t_dense.start();
M.multiply(x_dense,'N');
M.multiply(x_dense,'N');
t_dense.stop();
}
//float density = ((float)S.getNonZeros())/((float)(S.getNbRow()*S.getNbCol()));
float density = S.density();
//float density = ((float)S.getNonZeros())/((float)(S.getNbRow()*S.getNbCol()));
float density = S.density();
//std::cout<<" density "<<density<<std::endl;
//std::cout<<" tps sparse "<<t_sparse.get_time()<<std::endl;
//std::cout<<" tps dense "<<t_dense.get_time()<<std::endl;
......@@ -111,10 +109,10 @@ Faust::MatGeneric<FPP,Cpu>* Faust::optimize(Faust::MatDense<FPP,Cpu> const & M,F
return new MatSparse<FPP,Cpu>(S);
}else
{
//std::cout<<" CHOICE DENSE "<<t_dense.get_time()<<std::endl;
//std::cout<<" CHOICE DENSE "<<t_dense.get_time()<<std::endl;
return new MatDense<FPP,Cpu>(M);
}
//std::cout<<"FIN OPTIMIZE "<<t_sparse.get_time()<<std::endl;
......@@ -123,7 +121,6 @@ Faust::MatGeneric<FPP,Cpu>* Faust::optimize(Faust::MatDense<FPP,Cpu> const & M,F
template<typename FPP,Device DEVICE>
Faust::MatGeneric<FPP,DEVICE>::~MatGeneric()
{
}
......
......@@ -49,6 +49,7 @@ from libc.string cimport memcpy;
from libcpp cimport bool
from libcpp cimport complex
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
from scipy import sparse
cdef class FaustCore:
......@@ -70,9 +71,14 @@ cdef class FaustCore:
#print len(list_factors)
#print 'avant boucle'
self._isReal = True
for factor in list_factors:
for i,factor in enumerate(list_factors):
if(isinstance(factor, sparse.csc.csc_matrix)):
factor = list_factors[i] = factor.toarray()
#print("FaustCorePy.pyx __cinit__(),toarray() factor:", factor)
if(not isinstance(factor, np.ndarray)):
raise ValueError('Faust factors must be numpy.ndarray')
#print("FaustCorePy.pyx __cinit__(), factor:",factor)
raise ValueError("Faust factors must be a numpy.ndarray or "
"a scipy.sparse.csr.csr_matrix")
if(isinstance(factor[0,0], np.complex)):
# str(factor.dtype) == complex128/64/complex_
self._isReal = False
......@@ -85,6 +91,13 @@ cdef class FaustCore:
if(self._isReal):
data=factor.astype(float,'F')
else:
if(isinstance(factor, sparse.csc.csc_matrix)):
#TODO: understand how is it possible to have a sparse
# mat here and fix it (because it should have been
# converted above already)
factor = list_factors[i] = factor.toarray()
#print('FaustCorePy.pyx type factor=',type(factor))
#print("FaustCorePy.pyx factor=",factor)
data_cplx=factor.astype(np.complex128,'F')
nbrow=factor.shape[0];
nbcol=factor.shape[1];
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment