Handle also sparse matrices (csr_matrix) in list_factors argument of constructor Faust().

New function FaustCoreCpp::push_back() for pushing sparse matrix as Faust factor.

misc/test/src/Python/test_FaustPy.py

@@ -24,10 +24,13 @@ class TestFaustPy(unittest.TestCase):
         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])
+                        dtype=np.float64)] #.todense() removed
+            #print("factor",i,":", factors[i])
         self.F = Faust(factors)
         self.factors = factors
+        # we keep dense matrices as reference for the tests
+        for i in range(0, num_factors):
+            self.factors[i] =  factors[i].todense()
         print("Tests on random Faust with dims=", self.F.get_nb_rows(),
               self.F.get_nb_cols())
         print("Num. factors:", num_factors)
@@ -72,6 +75,8 @@ class TestFaustPy(unittest.TestCase):
     def testGetFactorAndConstructor(self):
         print("testGetFactorAndConstructor()")
         for ref_fact,i in zip(self.factors,range(0,len(self.factors))):
+            #print("testGetFac() fac ",i, " :", self.F.get_factor(i))
+            #print("testGetFac() reffac",i, ":",ref_fact)
             self.assertTrue((ref_fact == self.F.get_factor(i)).all())
 
     def testGetNumFactors(self):
@@ -253,8 +258,8 @@ class TestFaustPy(unittest.TestCase):
         print("Test Faust.getH()")
         test_Fct = self.F.getH().todense()
         ref_Fct = self.F.todense().conj().T
-        print("test_Fct=", test_Fct)
-        print("ref_Fct=", ref_Fct)
+        #print("test_Fct=", test_Fct)
+        #print("ref_Fct=", ref_Fct)
         ref_Fct[ref_Fct==0] = 1
         test_Fct[test_Fct==0] = 1
         self.assertTrue(((((test_Fct-ref_Fct)/ref_Fct) < 0.01)).all())
@@ -266,14 +271,28 @@ class TestFaustPyCplx(TestFaustPy):
             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 += [np.random.rand(d1, d2).astype(np.complex)]
-                factors[i].imag = [np.random.rand(d1, d2)]
+                d1, d2 = d2, r.randint(TestFaustPy.MIN_DIM_SIZE, TestFaustPy.MAX_DIM_SIZE)
+                if(r.randint(0,1) > 0): # generate a dense complex matrix
+                    factors += [np.random.rand(d1, d2).astype(np.complex)]
+                    factors[i].imag = [np.random.rand(d1, d2)]
+                else:
+                    # generate a sparse matrix
+                    # we can't use scipy.sparse.random directly because only float type is
+                    # supported for random sparse matrix generation
+                    factors += [sparse.random(d1, d2, dtype=np.float64, format='csr',
+                                       density=r.random()).astype(np.complex)]
+                    #print("setUp() i=",i, "d1=",d1, "d2=", d2, "factor.shape[0]=",
+                    #      factors[i].shape[0])
+                    factors[i] *= np.complex(r.random(), r.random())
             self.F = Faust(factors)
             self.factors = factors
-            print("Tests on random Faust with dims=", self.F.get_nb_rows(),
+            # we keep dense matrices as reference for the tests
+            for i in range(0, num_factors):
+                if(not isinstance(factors[i], np.ndarray)):
+                    self.factors[i] =  factors[i].todense()
+            print("Tests on random complex Faust with dims=", self.F.get_nb_rows(),
                   self.F.get_nb_cols())
             print("Num. factors:", num_factors)
             self.r = r

wrapper/python/src/FaustCoreCpp.h

@@ -59,6 +59,7 @@ class FaustCoreCpp
     //FaustCoreCpp(): transform(){}
     void Display() const { transform.display();}
     void push_back(FPP* valueMat,unsigned int nbrow,unsigned int nbcol);
+    void push_back(FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols);
     unsigned int getNbRow() const;
     unsigned int getNbCol() const;
     void multiply(FPP* value_y,int nbrow_y,int nbcol_y,FPP* value_x,int nbrow_x,int nbcol_x/*,bool isTranspose*/)const;

wrapper/python/src/FaustCoreCpp.hpp

@@ -50,6 +50,7 @@
 void FaustCoreCpp<FPP>::push_back(FPP* valueMat, unsigned int nbrow,unsigned int nbcol)
 {
     Faust::MatDense<FPP,Cpu> dense_mat(valueMat,nbrow,nbcol);
+    //TODO: why this auto-conversion in sparse mat ?
     Faust::MatSparse<FPP,Cpu> sparse_mat(dense_mat);
     //sparse_mat.Display();
     this->transform.push_back(&sparse_mat);
@@ -57,6 +58,15 @@ void FaustCoreCpp<FPP>::push_back(FPP* valueMat, unsigned int nbrow,unsigned int
 
 
 
+}
+
+template<typename FPP>
+void FaustCoreCpp<FPP>::push_back(FPP* data, int* row_ptr, int* id_col, int nnz, int nrows, int ncols)
+{
+    //    Faust::MatSparse<FPP,Cpu>::MatSparse(const faust_unsigned_int nnz_, const faust_unsigned_int dim1_, const faust_unsigned_int dim2_, const FPP1* value, const int* row_ptr, const int* id_col) :
+      Faust::MatSparse<FPP, Cpu> sparse_mat(nnz, nrows, ncols, data, row_ptr, id_col);
+//      std::cout << "FaustCoreCpp::push_back()" << nrows << " " << sparse_mat.getNbRow() <<  " " << ncols << " " << sparse_mat.getNbCol() << std::endl;
+      this->transform.push_back(&sparse_mat);
 }
 
 template<typename FPP>

wrapper/python/src/FaustCoreCy.pxd

@@ -44,6 +44,8 @@ cdef extern from "FaustCoreCpp.h" :
         FaustCoreCpp()
         void Display() const
         void push_back(FPP* valueMat,unsigned int nbrow,unsigned int nbcol)
+        void push_back(FPP* data, int* row_ptr, int* id_col, int nnz, int nrows,
+                       int ncols)
         void multiply(FPP* value_y,int nbrow_y,int nbcol_y,FPP* value_x,
                       int nbrow_x, int nbcol_x);#,bool isTranspose*/);
         unsigned int getNbRow() const

wrapper/python/src/FaustCorePy.pyx

@@ -61,24 +61,28 @@ cdef class FaustCore:
     #### CONSTRUCTOR ####
     #def __cinit__(self,np.ndarray[double, mode="fortran", ndim=2] mat):
     def  __cinit__(self,list_factors=None, core=False):
-        #print 'inside cinit'
         cdef double [:,:] data
+        cdef double [:] data1d #only for csr mat factor
+        cdef int [:] indices # only for csr mat
+        cdef int [:] indptr # only for csr mat
         cdef complex [:,:] data_cplx
+        cdef complex [:] data1d_cplx #only for csr mat factor
         cdef unsigned int nbrow
         cdef unsigned int nbcol
         if(list_factors is not None):
-            #print 'longueur list'
-            #print len(list_factors)
-            #print 'avant boucle'
             self._isReal = True
             for i,factor in enumerate(list_factors):
+                # Faust uses row-major order for sparse matrices
+                # and col-major order for dense matrices
+                # but libmatio uses col-major order for sparse matrices
                 if(isinstance(factor, sparse.csc.csc_matrix)):
-                    factor = list_factors[i] = factor.toarray()
+                    factor = list_factors[i] = factor.tocsr()
                     #print("FaustCorePy.pyx __cinit__(),toarray() factor:", factor)
-                if(not isinstance(factor, np.ndarray)):
-                    #print("FaustCorePy.pyx __cinit__(), factor:",factor)
-                    raise ValueError("Faust factors must be a numpy.ndarray or "
-                                     "a scipy.sparse.csr.csr_matrix")
+                if(not isinstance(factor, np.ndarray) and
+                   not isinstance(factor, sparse.csr.csr_matrix)):
+                   #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
@@ -88,27 +92,36 @@ cdef class FaustCore:
             else:
                 self.core_faust_cplx = new FaustCoreCy.FaustCoreCpp[complex]()
             for factor in list_factors:
+                nbrow=factor.shape[0];
+                nbcol=factor.shape[1];
                 if(self._isReal):
-                    data=factor.astype(float,'F')
+                   if(isinstance(factor, np.ndarray)):
+                      data=factor.astype(float,'F')
+                      self.core_faust_dbl.push_back(&data[0,0], nbrow, nbcol)
+                   else: #csr.csr_matrix
+                      data1d=factor.data.astype(float,'F')
+                      indices=factor.indices.astype(np.int32, 'F')
+                      indptr=factor.indptr.astype(np.int32, 'F')
+                      self.core_faust_dbl.push_back(&data1d[0], &indptr[0],
+                                                    &indices[0], factor.nnz, nbrow, nbcol)
                 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()
+                        factor = list_factors[i] = factor.tocsr()
                     #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];
-                #	print nbrow
-                #	print nbcol
-                if(self._isReal):
-                    self.core_faust_dbl.push_back(&data[0,0], nbrow, nbcol)
-                else:
-                    self.core_faust_cplx.push_back(&data_cplx[0,0], nbrow, nbcol)
-                #print(self.__dict__)
-                #print 'apres boucle'
+                    if(isinstance(factor, np.ndarray)):
+                        data_cplx=factor.astype(np.complex128,'F')
+                        self.core_faust_cplx.push_back(&data_cplx[0,0], nbrow, nbcol)
+                    else:
+                        #print("FaustCore, factor dims:", nbrow, nbcol)
+                        data1d_cplx = factor.data.astype(np.complex128, 'F')
+                        indices=factor.indices.astype(np.int32, 'F')
+                        indptr=factor.indptr.astype(np.int32, 'F')
+                        self.core_faust_cplx.push_back(&data1d_cplx[0], &indptr[0],
+                                                    &indices[0], factor.nnz, nbrow, nbcol)
         elif(core): # trick to initialize a new FaustCoreCpp from C++ (see
         # transpose, conj and adjoint)
             pass