Add py wrapper for factorization (both palm4msa and hierarchical fact.).

Adding unit tests for these algos also.

misc/test/src/Python/test_FaustPy.py

@@ -5,7 +5,7 @@ import tempfile
 import os
 import sys
 import numpy as np
-from scipy.io import savemat  # , loadmat
+from scipy.io import savemat,loadmat
 from numpy.linalg import norm
 import math
 
@@ -301,6 +301,80 @@ class TestFaustPyCplx(TestFaustPy):
             self.r = r
             self.num_factors = num_factors
 
+class TestFaustFactory(unittest.TestCase):
+
+    def testFactPalm4MSA(self):
+        print("Test FaustFactory.fact_palm4msa()")
+        from pyfaust import FaustFactory, ParamsPalm4MSA, ConstraintScalar,\
+        ConstraintName, StoppingCriterion
+        num_facts = 2
+        is_update_way_R2L = False
+        init_lambda = 1.0
+        init_facts = list()
+        init_facts.append(np.zeros([500,32]))
+        init_facts.append(np.eye(32))
+        #M = np.random.rand(500, 32)
+        M = \
+        loadmat(sys.path[-1]+"/../../../misc/data/mat/config_compared_palm2.mat")['data']
+        # default step_size
+        cons1 = ConstraintScalar(ConstraintName(ConstraintName.SPLIN), 500, 32, 5)
+        cons2 = ConstraintScalar(ConstraintName(ConstraintName.NORMCOL), 32,
+                                 32, 1.0)
+        stop_crit = StoppingCriterion(num_its=200)
+        param = ParamsPalm4MSA(num_facts, is_update_way_R2L, init_lambda,
+                               init_facts, [cons1, cons2], stop_crit,
+                               is_verbose=False, constant_step_size=False)
+        F = FaustFactory.fact_palm4msa(M, param)
+        self.assertEqual(F.shape, M.shape)
+        E = F.todense()-M
+        #print("err.:",norm(F.todense(), "fro"),  norm(E,"fro"), norm (M,"fro"))
+        print("err: ", norm(E,"fro")/norm(M,"fro"))
+        # matrix to factorize and reference relative error come from
+        # misc/test/src/C++/test_palm4MSA.cpp
+        self.assertAlmostEqual(norm(E,"fro")/norm(M,"fro"), 0.98162, places=4)
+
+    def testFactHierarch(self):
+        print("Test FaustFactory.fact_hierarchical()")
+        from pyfaust import FaustFactory, ParamsHierarchicalFact, ConstraintScalar,\
+        ConstraintName, StoppingCriterion
+        num_facts = 4
+        is_update_way_R2L = False
+        init_lambda = 1.0
+        init_facts = list()
+        init_facts.append(np.zeros([500,32]))
+        for i in range(1,num_facts):
+            init_facts.append(np.zeros([32,32]))
+        #M = np.random.rand(500, 32)
+        M = \
+        loadmat(sys.path[-1]+"/../../../misc/data/mat/matrix_hierarchical_fact.mat")['matrix']
+        # default step_size
+        cons1 = ConstraintScalar(ConstraintName(ConstraintName.SPLIN), 500, 32, 5)
+        cons2 = ConstraintScalar(ConstraintName(ConstraintName.SP), 32,
+                                 32, 96)
+        cons3 = ConstraintScalar(ConstraintName(ConstraintName.SP), 32, 32,
+                                 96)
+        cons4 = ConstraintScalar(ConstraintName(ConstraintName.NORMCOL), 32, 32,
+                                 1)
+        cons5 =  ConstraintScalar(ConstraintName(ConstraintName.SP), 32, 32,
+                                 666)
+        cons6 =  ConstraintScalar(ConstraintName(ConstraintName.SP), 32, 32,
+                                 333)
+        stop_crit1 = StoppingCriterion(num_its=200)
+        stop_crit2 = StoppingCriterion(num_its=200)
+        param = ParamsHierarchicalFact(num_facts, is_update_way_R2L, init_lambda,
+                                       init_facts, [cons1, cons2, cons3, cons4,
+                                                   cons5, cons6],
+                                       M.shape[0],M.shape[1],[stop_crit1,
+                                                              stop_crit2],
+                                       is_verbose=False)
+        F = FaustFactory.fact_hierarchical(M, param)
+        self.assertEqual(F.shape, M.shape)
+        E = F.todense()-M
+        #print("err.:",norm(F.todense(), "fro"),  norm(E,"fro"), norm (M,"fro"),
+        print("err: ", norm(E,"fro")/norm(M,"fro"))
+        # matrix to factorize and reference relative error come from
+        # misc/test/src/C++/hierarchicalFactorization.cpp
+        self.assertAlmostEqual(norm(E,"fro")/norm(M,"fro"),0.268513, places=5)
 
 if __name__ == "__main__":
     if(len(sys.argv)> 1):

wrapper/python/pyfaust.py

@@ -116,7 +116,7 @@ class Faust:
             scale: a multiplicative scalar (see examples below).
             core_obj: for internal purpose only. Please don't fill this argument.
 
-        WARNING: filepath and factors arguments are multually exclusive. If you
+        WARNING: filepath and factors arguments are mutually exclusive. If you
         use filepath you must explicitely set argument with the keyword.
 
         Examples:
@@ -195,11 +195,11 @@ class Faust:
             >>> F = Faust.rand(2, 10, .5, is_real=False)
             >>> G = Faust.rand([2, 5], [10, 20], .5, fac_type="dense")
             >>> F
-            Faust size 10x10, density 0.99, nnz_sum 99, 2 factor(s):<br/>
+            Faust size 10x10, density 0.99, nnz_sum 99, 2 factors:<br/>
             FACTOR 0 (complex) SPARSE, size 10x10, density 0.4, nnz 40<br/>
             FACTOR 1 (complex) DENSE, size 10x10, density 0.59, nnz 59<br/>
             >>> G
-            Faust size 19x16, density 1.37171, nnz_sum 417, 4 factor(s):<br/>
+            Faust size 19x16, density 1.37171, nnz_sum 417, 4 factors:<br/>
             FACTOR 0 (real) DENSE, size 19x17, density 0.49226, nnz 159<br/>
             FACTOR 1 (real) DENSE, size 17x10, density 0.517647, nnz 88<br/>
             FACTOR 2 (real) DENSE, size 10x13, density 0.515385, nnz 67<br/>
@@ -457,12 +457,12 @@ class Faust:
             >>> from pyfaust import Faust
             >>> F = Faust.rand([1, 2], [50, 100], .5)
             >>> F.display()
-            Faust size 98x82, density 0.686909, nnz_sum 5520, 2 factor(s):<br/>
+            Faust size 98x82, density 0.686909, nnz_sum 5520, 2 factors:<br/>
             FACTOR 0 (real) SPARSE, size 98x78, density 0.395081, nnz 3020<br/>
             FACTOR 1 (real) SPARSE, size 78x82, density 0.390869, nnz 2500<br/>
 
             >>> F
-            Faust size 98x82, density 0.686909, nnz_sum 5520, 2 factor(s):<br/>
+            Faust size 98x82, density 0.686909, nnz_sum 5520, 2 factors:<br/>
             FACTOR 0 (real) SPARSE, size 98x78, density 0.395081, nnz 3020<br/>
             FACTOR 1 (real) SPARSE, size 78x82, density 0.390869, nnz 2500<br/>
             <!-- >>> -->
@@ -544,7 +544,7 @@ class Faust:
         Examples:
             >>> from pyfaust import Faust
             >>> F = Faust.rand(5, 10**6, .00001, 'sparse')
-               Faust size 1000000x1000000, density 5e-05, nnz_sum 49999995, 5 factor(s):<br/>
+               Faust size 1000000x1000000, density 5e-05, nnz_sum 49999995, 5 factors:<br/>
                 FACTOR 0 (real) SPARSE, size 1000000x1000000, density 1e-05, nnz 9999999<br/>
                 FACTOR 1 (real) SPARSE, size 1000000x1000000, density 1e-05, nnz 9999999<br/>
                 FACTOR 2 (real) SPARSE, size 1000000x1000000, density 1e-05, nnz 9999999<br/>
@@ -582,7 +582,7 @@ class Faust:
             >>> from pyfaust import Faust
             >>> F = Faust.rand(5, 10**6, .00001, 'sparse')
             >>> F
-            Faust size 1000000x1000000, density 5e-05, nnz_sum 49999995, 5 factor(s):<br/>
+            Faust size 1000000x1000000, density 5e-05, nnz_sum 49999995, 5 factors:<br/>
             FACTOR 0 (real) SPARSE, size 1000000x1000000, density 1e-05, nnz 9999999<br/>
             FACTOR 1 (real) SPARSE, size 1000000x1000000, density 1e-05, nnz 9999999<br/>
             FACTOR 2 (real) SPARSE, size 1000000x1000000, density 1e-05, nnz 9999999<br/>
@@ -926,3 +926,233 @@ class Faust:
             return np.dtype(np.float64)
         else:
             return np.dtype(np.complex)
+
+class FaustFactory:
+
+    @staticmethod
+    def fact_palm4msa(M, p):
+        """
+        Factorizes the matrix M using the parameters set in p.
+
+        Args:
+            M: the numpy matrix to factorize.
+            p: the ParamsPalm4MSA instance to define the algorithm parameters.
+
+        Returns:
+            A Faust object as the result of the factorization.
+
+        Examples:
+        >>> from pyfaust import FaustFactory, ParamsPalm4MSA, ConstraintScalar,\
+        >>> ConstraintName, StoppingCriterion
+        >>> import numpy as np
+        >>> num_facts = 2
+        >>> is_update_way_R2L = False
+        >>> init_lambda = 1.0
+        >>> init_facts = list()
+        >>> init_facts.append(np.zeros([500,32]))
+        >>> init_facts.append(np.eye(32))
+        >>> M = np.random.rand(500, 32)
+        >>> cons1 = ConstraintScalar(ConstraintName(ConstraintName.SPLIN), 500, 32, 5)
+        >>> cons2 = ConstraintScalar(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1.0)
+        >>> stop_crit = StoppingCriterion(num_its=200)
+        >>> # default step_size is 1e-16
+        >>> param = ParamsPalm4MSA(num_facts, is_update_way_R2L, init_lambda,
+        >>>                        init_facts, [cons1, cons2], stop_crit,
+        >>>                        is_verbose=False)
+        >>> F = FaustFactory.fact_palm4msa(M, param)
+        >>> F.display()
+        Faust size 500x32, density 0.22025, nnz_sum 3524, 2 factor(s):<br/>
+        FACTOR 0 (real) SPARSE, size 500x32, density 0.15625, nnz 2500<br/>
+        FACTOR 1 (real) SPARSE, size 32x32, density 1, nnz 1024<br/>
+        >>> # verify if NORMCOL constraint is respected by the 2nd factor
+        >>> from numpy.linalg import norm
+        >>> for i in range(0, F.get_factor(1).shape[1]):
+        >>>     if(i<F.get_factor(1).shape[1]-1):
+        >>>             end=' '
+        >>>     else:
+        >>>             end='\n'
+        >>>     print(norm(F.get_factor(1)[:,i]), end=end)
+        1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
+        """
+        if(not isinstance(p, ParamsPalm4MSA)):
+            raise ValueError("p must be a ParamsPalm4MSA object.")
+        return Faust(core_obj=FaustCorePy.FaustFact.fact_palm4MSA(M, p))
+
+    @staticmethod
+    def fact_hierarchical(M, p):
+        """
+        Factorizes the matrix M using the parameters set in p.
+
+        Args:
+            M: the numpy matrix to factorize.
+            p: the ParamsHierarchicalFact instance to define the algorithm parameters.
+
+        Returns:
+            A Faust object as the result of the factorization.
+
+        Examples:
+		>>> from pyfaust import FaustFactory, ParamsHierarchicalFact, ConstraintScalar,\
+		>>> ConstraintName, StoppingCriterion
+		>>> import numpy as np
+		>>> num_facts = 4
+		>>> is_update_way_R2L = False
+		>>> init_lambda = 1.0
+		>>> init_facts = list()
+		>>> init_facts.append(np.zeros([500,32]))
+		>>> for i in range(1,num_facts):
+		>>>     init_facts.append(np.zeros([32,32]))
+		>>> M = np.random.rand(500, 32)
+		>>> cons1 = ConstraintScalar(ConstraintName(ConstraintName.SPLIN), 500, 32, 5)
+		>>> cons2 = ConstraintScalar(ConstraintName(ConstraintName.SP), 32, 32, 96)
+		>>> cons3 = ConstraintScalar(ConstraintName(ConstraintName.SP), 32, 32, 96)
+		>>> cons4 = ConstraintScalar(ConstraintName(ConstraintName.NORMCOL), 32, 32, 1)
+		>>> cons5 =  ConstraintScalar(ConstraintName(ConstraintName.SP), 32, 32, 666)
+		>>> cons6 =  ConstraintScalar(ConstraintName(ConstraintName.SP), 32, 32, 333)
+		>>> stop_crit1 = StoppingCriterion(num_its=200)
+		>>> stop_crit2 = StoppingCriterion(num_its=200)
+		>>> param = ParamsHierarchicalFact(num_facts, is_update_way_R2L, init_lambda,
+		>>>                                init_facts, [cons1, cons2, cons3, cons4,
+		>>>                                            cons5, cons6],
+		>>>                                M.shape[0],M.shape[1],[stop_crit1,
+		>>>                                                       stop_crit2],
+		>>>                                is_verbose=False)
+		>>> F = FaustFactory.fact_hierarchical(M, param)
+		Faust::HierarchicalFact<FPP,DEVICE>::compute_facts : factorisation
+        1/3<br/>
+		Faust::HierarchicalFact<FPP,DEVICE>::compute_facts : factorisation
+        2/3<br/>
+		Faust::HierarchicalFact<FPP,DEVICE>::compute_facts : factorisation
+        3/3<br/>
+		>>> F.display()
+		Faust size 500x32, density 0.189063, nnz_sum 3025, 4 factor(s):
+		- FACTOR 0 (real) SPARSE, size 500x32, density 0.15625, nnz 2500
+		- FACTOR 1 (real) SPARSE, size 32x32, density 0.09375, nnz 96
+		- FACTOR 2 (real) SPARSE, size 32x32, density 0.09375, nnz 96
+		- FACTOR 3 (real) DENSE, size 32x32, density 0.325195, nnz 333
+        """
+        if(not isinstance(p, ParamsHierarchicalFact)):
+            raise ValueError("p must be a ParamsPalm4MSA object.")
+        return Faust(core_obj=FaustCorePy.FaustFact.fact_hierarchical(M, p))
+
+class ParamsFact(object):
+
+    def __init__(self, num_facts, is_update_way_R2L, init_lambda, init_facts,
+                 constraints, step_size, constant_step_size = False, is_verbose = False,
+                 ):
+        self.num_facts = num_facts
+        self.is_update_way_R2L = is_update_way_R2L
+        self.init_lambda = init_lambda
+        self.init_facts = init_facts
+        self.step_size = step_size
+        self.constraints = constraints
+        self.is_verbose = is_verbose
+        self.constant_step_size = constant_step_size
+        #TODO: raise exception if num_facts != len(init_facts)
+        #TODO: likewise for constraints
+
+class ParamsPalm4MSA(ParamsFact):
+
+    def __init__(self, num_facts, is_update_way_R2L, init_lambda, init_facts,
+                 constraints, stop_crit, step_size = 10.0**-16,
+                 constant_step_size = False,
+                 is_verbose = False):
+        super(ParamsPalm4MSA, self).__init__(num_facts, is_update_way_R2L, init_lambda, init_facts,
+                                             constraints, step_size,
+                                             constant_step_size,
+                                             is_verbose)
+        self.stop_crit = stop_crit
+        # TODO: raise exception if stop_crit not a StoppingCriterion object
+
+class ParamsHierarchicalFact(ParamsFact):
+
+    def __init__(self, num_facts, is_update_way_R2L, init_lambda, init_facts,
+                 constraints, data_num_rows, data_num_cols, stop_crits,
+                 step_size = 10.0**-16, constant_step_size = False,
+                 is_verbose = False,
+                 is_fact_side_left = False):
+        super(ParamsHierarchicalFact, self).__init__(num_facts, is_update_way_R2L, init_lambda, init_facts,
+                                                     constraints, step_size,
+                                                     constant_step_size,
+                                                     is_verbose)
+        self.data_num_rows = data_num_rows
+        self.data_num_cols = data_num_cols
+        self.stop_crits = stop_crits
+        self.is_fact_side_left = is_fact_side_left
+        # TODO: raise exception if stop_crits len != 2 or not StoppingCriterion
+        # objects
+
+class StoppingCriterion(object):
+
+    def __init__(self, is_criterion_error = False , error_treshold = 0.3,
+                 num_its = 500,
+                 max_num_its = 1000):
+        self.is_criterion_error = is_criterion_error
+        self.error_treshold = error_treshold
+        self.num_its = num_its
+        self.max_num_its = max_num_its
+        #TODO: check_validity() like C++ code does
+
+class ConstraintName:
+    SP = 0 # Int Constraint
+    SPCOL = 1 # Int Constraint
+    SPLIN=2 # Int Constraint
+    NORMCOL = 3 # Real Constraint
+    SPLINCOL = 4 # Int Constraint
+    CONST = 5 # Mat Constraint
+    SP_POS = 6 # Int Constraint
+    #BLKDIAG = 7 # ?? Constraint #TODO
+    SUPP = 8 # Mat Constraint
+    NORMLIN = 9 # Real Constraint
+
+    def __init__(self, name):
+        if(name < ConstraintName.SP or name > ConstraintName.NORMLIN):
+            raise ValueError("name must be an integer among ConstraintName.SP,"
+                             "ConstraintName.SPCOL, ConstraintName.NORMCOL,"
+                             "ConstraintName.SPLINCOL, ConstraintName.CONST,"
+                             "ConstraintName.SP_POS," # ConstraintName.BLKDIAG,
+                             "ConstraintName.SUPP, ConstraintName.NORMLIN")
+        self.name = name
+
+    def is_int_constraint(self):
+        return self.name in [ ConstraintName.SP, ConstraintName.SPCOL,
+                             ConstraintName.SPLIN, ConstraintName.SPLINCOL,
+                             ConstraintName.SP_POS ]
+
+    def is_real_constraint(self):
+        return self.name in [ ConstraintName.NORMCOL, ConstraintName.NORMLIN ]
+
+    def is_mat_constraint(self):
+        return self.name in [ConstraintName.SUPP, ConstraintName.CONST ]
+
+class ConstraintGeneric(object):
+
+    def __init__(self, name , num_rows, num_cols):
+        self.value = name
+        self.num_rows = num_rows
+        self.num_cols = num_cols
+
+    @property
+    def name(self):
+        return self.value.name
+
+    def is_int_constraint(self):
+        return self.value.is_int_constraint()
+
+    def is_real_constraint(self):
+        return self.value.is_real_constraint()
+
+    def is_mat_constraint(self):
+        return self.value.is_mat_constraint()
+
+class ConstraintScalar(ConstraintGeneric):
+
+    def __init__(self, name, num_rows, num_cols, param):
+        super(ConstraintScalar, self).__init__(name, num_rows, num_cols)
+        self.param = param # raise value error if not np.number (can be complex
+        # or float/double or integer
+
+class ConstraintMat(ConstraintGeneric):
+
+    def __init__(self, name, num_rows, num_cols, param):
+        super(ConstraintMat, self).__init__(name, num_rows, num_cols)
+        self.param = param #TODO: raise ValueError if not np.ndarray

wrapper/python/src/FaustCoreCpp.h

@@ -45,7 +45,6 @@
 #ifndef FAUSTCORECPP_H
 #define FAUSTCORECPP_H
 
-//#include "faust_transform.h"
 #include "faust_MatDense.h"
 #include "faust_TransformHelper.h"
 
@@ -56,7 +55,8 @@ class FaustCoreCpp
     public :
 
 
-    //FaustCoreCpp(): transform(){}
+    FaustCoreCpp(): transform(){}
+    FaustCoreCpp(Faust::TransformHelper<FPP,Cpu> &th);
     void Display() const { transform.display();}
     const char* to_string() const;
     void push_back(FPP* valueMat,unsigned int nbrow,unsigned int nbcol);

wrapper/python/src/FaustCoreCpp.hpp

@@ -46,6 +46,15 @@
 #include <iostream>
 #include <exception>
 
+
+template<typename FPP>
+FaustCoreCpp<FPP>::FaustCoreCpp(Faust::TransformHelper<FPP,Cpu> &th)
+{
+    this->transform = th;
+}
+
+
+
     template<typename FPP>
 void FaustCoreCpp<FPP>::push_back(FPP* valueMat, unsigned int nbrow,unsigned int nbcol)
 {

wrapper/python/src/FaustCoreCy.pxd

@@ -73,3 +73,52 @@ cdef extern from "FaustCoreCpp.h" :
                                              unsigned int min_num_factors, unsigned int max_num_factors,
                                              unsigned int min_dim_size,
                                              unsigned int max_dim_size, float density)
+
+cdef extern from "FaustFact.h":
+    cdef cppclass PyxConstraintGeneric:
+        int name
+        unsigned long num_rows
+        unsigned long num_cols
+
+    cdef cppclass PyxConstraintInt(PyxConstraintGeneric):
+        unsigned long parameter
+
+    cdef cppclass PyxConstraintScalar[FPP](PyxConstraintGeneric):
+        FPP parameter
+
+    cdef cppclass PyxConstraintMat[FPP](PyxConstraintGeneric):
+        FPP* parameter # shape = num_rows, num_cols
+
+    cdef cppclass PyxStoppingCriterion[FPP]:
+        bool is_criterion_error
+        int num_its
+        FPP error_treshold
+        unsigned long max_num_its
+
+    cdef cppclass PyxParamsFact[FPP]:
+        int num_facts
+        bool is_update_way_R2L
+        FPP init_lambda
+        FPP** init_facts # num_facts elts
+        unsigned long* init_fact_sizes
+        FPP step_size
+        PyxConstraintGeneric** constraints # (num_facts-1)*2 elts
+        unsigned int num_constraints
+        bool is_verbose
+        bool constant_step_size
+
+    cdef cppclass PyxParamsFactPalm4MSA[FPP](PyxParamsFact[FPP]):
+        PyxStoppingCriterion stop_crit
+
+    cdef cppclass PyxParamsHierarchicalFact[FPP](PyxParamsFact[FPP]):
+        unsigned int num_rows
+        unsigned int num_cols
+        PyxStoppingCriterion* stop_crits #must be of size 2
+        bool is_fact_side_left
+
+    cdef FaustCoreCpp[FPP]* fact_palm4MSA[FPP](FPP*,unsigned int, unsigned int,
+                                             PyxParamsFactPalm4MSA[FPP]*)
+
+    cdef FaustCoreCpp[FPP]* fact_hierarchical[FPP](FPP*,unsigned int, unsigned int,
+                                             PyxParamsHierarchicalFact[FPP]*)
+

wrapper/python/src/FaustCorePy.pyx

@@ -39,6 +39,8 @@
 
 # importer le module Cython faisant le lien avec la class C++
 cimport FaustCoreCy
+from FaustCoreCy cimport PyxConstraintGeneric, PyxConstraintInt, \
+PyxConstraintMat, PyxConstraintScalar, PyxStoppingCriterion, PyxParamsFactPalm4MSA
 
 import numpy as np
 cimport numpy as np
@@ -47,8 +49,7 @@ import copy
 from libc.stdlib cimport malloc, free;
 from libc.stdio cimport printf
 from libc.string cimport memcpy, strlen;
-from libcpp cimport bool
-from libcpp cimport complex
+from libcpp cimport bool, complex
 from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
 from scipy import sparse
 from re import match
@@ -205,13 +206,14 @@ cdef class FaustCore:
         if not isinstance(M, (np.ndarray) ):
             raise ValueError('input M must a numpy ndarray')
         if(self._isReal):
-            M=M.astype(float,'F')
-            if not M.dtype=='float':
-                raise ValueError('input M must be double array')
+           M=M.astype(float,'F')
+           if not M.dtype=='float':
+               raise ValueError('input M must be double array')
         else:
-            M=M.astype(complex,'F')
-            if(M.dtype not in ['complex', 'complex128', 'complex64'] ): #could fail if complex128 etc.
-                raise ValueError('input M must be complex array')
+           M=M.astype(complex,'F')
+           if(M.dtype not in ['complex', 'complex128', 'complex64'] ): #could fail if complex128 etc.
+               raise ValueError('input M must be complex array')
+        #TODO: raise exception if not real nor complex
         if not M.flags['F_CONTIGUOUS']:
             raise ValueError('input M must be Fortran contiguous (Colmajor)')
 
@@ -219,7 +221,8 @@ cdef class FaustCore:
 
         if (ndim_M > 2) | (ndim_M < 1):
             raise ValueError('input M invalid number of dimensions')
-
+        check_matrix(self._isReal, M)
+        ndim_M=M.ndim
         cdef unsigned int nbrow_x=M.shape[0]
         cdef unsigned int nbcol_x #can't be assigned because we don't know yet if the input vector is 1D or 2D
 
@@ -279,7 +282,7 @@ cdef class FaustCore:
 
         return y
 
-
+    
 
     # print information about the faust (size, number of factor, type of factor (dense/sparse) ...)
     def display(self):
@@ -419,3 +422,319 @@ cdef class FaustCore:
             del self.core_faust_dbl
         else:
             del self.core_faust_cplx
+
+cdef check_matrix(isReal, M):
+        if not isinstance(M, (np.ndarray) ):
+            raise ValueError('input M must a numpy ndarray')
+        if(isReal):
+            M=M.astype(float,'F')
+            if not M.dtype=='float':
+                raise ValueError('input M must be double array')
+        else:
+            M=M.astype(complex,'F')
+            if(M.dtype not in ['complex', 'complex128', 'complex64'] ): #could fail if complex128 etc.
+                raise ValueError('input M must be complex array')
+        #TODO: raise exception if not real nor complex
+        if not M.flags['F_CONTIGUOUS']:
+            raise ValueError('input M must be Fortran contiguous (Colmajor)')
+
+        ndim_M=M.ndim;
+
+        if (ndim_M > 2) | (ndim_M < 1):
+            raise ValueError('input M invalid number of dimensions')
+
+
+cdef class FaustFact:
+
+    @staticmethod
+    def fact_palm4MSA(M, p):
+        isReal = M.dtype in [ 'float', 'float128',
+                             'float16', 'float32',
+                             'float64', 'double']
+        # double == float64
+        #TODO: it not float nor complex, raise exception
+        check_matrix(isReal, M)
+
+        cdef unsigned int M_num_rows=M.shape[0]
+        cdef unsigned int M_num_cols=M.shape[1]
+
+        cdef double[:,:] Mview
+        cdef complex[:,:] Mview_cplx
+
+        cdef double[:,:] tmp_mat
+        cdef complex[:,:] tmp_mat_cplx
+
+        cdef FaustCoreCy.PyxParamsFactPalm4MSA[double] cpp_params
+        cdef FaustCoreCy.PyxParamsFactPalm4MSA[complex] cpp_params_cplx
+        cdef PyxStoppingCriterion[double] cpp_stop_crit
+        # template parameter is always double (never complex) because no need
+        # to have a treshold error of complex type (it wouldn't make sense)
+        cdef PyxConstraintGeneric** cpp_constraints
+
+        cpp_stop_crit.is_criterion_error = p.stop_crit.is_criterion_error
+        cpp_stop_crit.error_treshold = p.stop_crit.error_treshold
+        cpp_stop_crit.num_its = p.stop_crit.num_its
+        cpp_stop_crit.max_num_its = p.stop_crit.max_num_its
+
+        if(isReal):
+            Mview=M
+            cpp_params.num_facts = p.num_facts
+            cpp_params.is_update_way_R2L = p.is_update_way_R2L
+            cpp_params.init_lambda = p.init_lambda
+            cpp_params.step_size = p.step_size
+            cpp_params.stop_crit = cpp_stop_crit
+            cpp_params.init_facts = <double**> \
+                    PyMem_Malloc(sizeof(double*)*p.num_facts)
+            cpp_params.init_fact_sizes = <unsigned long*> \
+            PyMem_Malloc(sizeof(unsigned long)*2*p.num_facts)
+            cpp_params.is_verbose = p.is_verbose
+            cpp_params.constant_step_size = p.constant_step_size
+        else:
+            Mview_cplx=M
+            cpp_params_cplx.num_facts = p.num_facts
+            cpp_params_cplx.num_facts = p.num_facts
+            cpp_params_cplx.is_update_way_R2L = p.is_update_way_R2L
+            cpp_params_cplx.init_lambda = p.init_lambda
+            cpp_params_cplx.step_size = p.step_size
+            cpp_params_cplx.init_facts = <complex**> \
+                    PyMem_Malloc(sizeof(complex*)*p.num_facts)
+            cpp_params_cplx.init_fact_sizes = <unsigned long*> \
+            PyMem_Malloc(sizeof(unsigned long)*2*p.num_facts)
+            cpp_params_cplx.is_verbose = p.is_verbose
+            cpp_params_cplx.constant_step_size = p.constant_step_size
+
+
+        cpp_constraints = \
+        <PyxConstraintGeneric**> \
+        PyMem_Malloc(sizeof(PyxConstraintGeneric*)*len(p.constraints))
+
+        for i in range(0,len(p.constraints)):
+            cons = p.constraints[i]
+            #print("FaustFact.fact_palm4MSA() cons.name =", cons.name)
+            if(cons.is_int_constraint()):
+                #print("FaustFact.fact_palm4MSA() Int Constraint")
+                cpp_constraints[i] = <PyxConstraintInt*> PyMem_Malloc(sizeof(PyxConstraintInt))
+                (<PyxConstraintInt*>cpp_constraints[i]).parameter = cons.param
+            elif(cons.is_real_constraint()):
+                #print("FaustFact.fact_palm4MSA() Real Constraint")
+                cpp_constraints[i] = <PyxConstraintScalar[double]*> \
+                PyMem_Malloc(sizeof(PyxConstraintScalar[double]))
+                (<PyxConstraintScalar[double]*>cpp_constraints[i]).parameter =\
+                        cons.param
+            elif(cons.is_mat_constraint()):
+                #print("FaustFact.fact_palm4MSA() Matrix Constraint")
+                if(isReal):
+                    cpp_constraints[i] = <PyxConstraintMat[double]*> \
+                            PyMem_Malloc(sizeof(PyxConstraintMat[double]))
+                    tmp_mat = cons.param
+                    (<PyxConstraintMat[double]*>cpp_constraints[i]).parameter =\
+                            &tmp_mat[0,0]
+                else:
+                    cpp_constraints[i] = <PyxConstraintMat[complex]*> \
+                            PyMem_Malloc(sizeof(PyxConstraintMat[complex]))
+                    tmp_mat_cplx = cons.param
+                    (<PyxConstraintMat[complex]*>cpp_constraints[i]).parameter =\
+                            &tmp_mat_cplx[0,0]
+
+            else:
+                raise ValueError("Constraint type/name is not recognized.")
+            cpp_constraints[i].name = cons.name
+            cpp_constraints[i].num_rows = cons.num_rows
+            cpp_constraints[i].num_cols = cons.num_cols
+
+
+        if(isReal):
+            cpp_params.constraints = cpp_constraints
+            cpp_params.num_constraints = len(p.constraints)
+        else:
+            cpp_params_cplx.constraints = cpp_constraints
+            cpp_params_cplx.num_constraints = len(p.constraints)
+
+        for i in range(0,p.num_facts):
+            if(isReal):
+                tmp_mat = p.init_facts[i]
+                cpp_params.init_facts[i] = &tmp_mat[0,0]
+                cpp_params.init_fact_sizes[i*2+0] = p.init_facts[i].shape[0]
+                cpp_params.init_fact_sizes[i*2+1] = p.init_facts[i].shape[1]
+            else:
+                tmp_mat_cplx = p.init_facts[i]
+                cpp_params_cplx.init_facts[i] = &tmp_mat_cplx[0,0]
+                cpp_params_cplx.init_fact_sizes[i*2+0] = p.init_facts[i].shape[0]
+                cpp_params_cplx.init_fact_sizes[i*2+1] = p.init_facts[i].shape[1]
+
+        core = FaustCore(core=True)
+        if(isReal):
+            core.core_faust_dbl = FaustCoreCy.fact_palm4MSA(&Mview[0,0], M_num_rows, M_num_cols,
+ #           FaustCoreCy.fact_palm4MSA(&Mview[0,0], M_num_rows, M_num_cols,
+                                      &cpp_params)
+            core._isReal = True
+            #TODO: FPP == complex not yet supported by C++ code
+#        else:
+#            core.core_faust_cplx = FaustCoreCy.fact_palm4MSA(&Mview_cplx[0,0], M_num_rows, M_num_cols,
+#                                     &cpp_params_cplx)
+#            core._isReal = False
+        for i in range(0,len(p.constraints)):
+            PyMem_Free(cpp_constraints[i])
+        PyMem_Free(cpp_constraints)
+        if(isReal):
+            PyMem_Free(cpp_params.init_facts)
+            PyMem_Free(cpp_params.init_fact_sizes)
+        else:
+            PyMem_Free(cpp_params_cplx.init_facts)
+            PyMem_Free(cpp_params_cplx.init_fact_sizes)
+#
+        return core
+
+    @staticmethod
+    def fact_hierarchical(M, p):
+        isReal = M.dtype in [ 'float', 'float128',
+                             'float16', 'float32',
+                             'float64', 'double']
+        # double == float64
+        #TODO: it not float nor complex, raise exception
+        check_matrix(isReal, M)
+
+        cdef unsigned int M_num_rows=M.shape[0]
+        cdef unsigned int M_num_cols=M.shape[1]
+
+        cdef double[:,:] Mview
+        cdef complex[:,:] Mview_cplx
+
+        cdef double[:,:] tmp_mat
+        cdef complex[:,:] tmp_mat_cplx
+
+        cdef FaustCoreCy.PyxParamsHierarchicalFact[double] cpp_params
+        cdef FaustCoreCy.PyxParamsHierarchicalFact[complex] cpp_params_cplx
+        cdef PyxStoppingCriterion[double]* cpp_stop_crits
+        # template parameter is always double (never complex) because no need
+        # to have a treshold error of complex type (it wouldn't make sense)
+        cdef PyxConstraintGeneric** cpp_constraints
+
+        cpp_stop_crits = <PyxStoppingCriterion[double]*>\
+        PyMem_Malloc(sizeof(PyxStoppingCriterion[double])*2)
+
+        cpp_stop_crits[0].is_criterion_error = p.stop_crits[0].is_criterion_error
+        cpp_stop_crits[0].error_treshold = p.stop_crits[0].error_treshold
+        cpp_stop_crits[0].num_its = p.stop_crits[0].num_its
+        cpp_stop_crits[0].max_num_its = p.stop_crits[0].max_num_its
+        cpp_stop_crits[1].is_criterion_error = p.stop_crits[1].is_criterion_error
+        cpp_stop_crits[1].error_treshold = p.stop_crits[1].error_treshold
+        cpp_stop_crits[1].num_its = p.stop_crits[1].num_its
+        cpp_stop_crits[1].max_num_its = p.stop_crits[1].max_num_its
+
+
+
+        if(isReal):
+            Mview=M
+            cpp_params.num_facts = p.num_facts
+            cpp_params.is_update_way_R2L = p.is_update_way_R2L
+            cpp_params.init_lambda = p.init_lambda
+            cpp_params.step_size = p.step_size
+            cpp_params.stop_crits = cpp_stop_crits
+            cpp_params.init_facts = <double**> \
+                    PyMem_Malloc(sizeof(double*)*p.num_facts)
+            cpp_params.init_fact_sizes = <unsigned long*> \
+            PyMem_Malloc(sizeof(unsigned long)*2*p.num_facts)
+            cpp_params.is_verbose = p.is_verbose
+            cpp_params.is_fact_side_left = p.is_fact_side_left
+            cpp_params.constant_step_size = p.constant_step_size
+        else:
+            Mview_cplx=M
+            cpp_params_cplx.num_facts = p.num_facts
+            cpp_params_cplx.num_facts = p.num_facts
+            cpp_params_cplx.is_update_way_R2L = p.is_update_way_R2L
+            cpp_params_cplx.init_lambda = p.init_lambda
+            cpp_params_cplx.step_size = p.step_size
+            cpp_params_cplx.init_facts = <complex**> \
+                    PyMem_Malloc(sizeof(complex*)*p.num_facts)
+            cpp_params_cplx.init_fact_sizes = <unsigned long*> \
+            PyMem_Malloc(sizeof(unsigned long)*2*p.num_facts)
+            cpp_params_cplx.is_verbose = p.is_verbose
+            cpp_params_cplx.is_fact_side_left = p.is_fact_side_left
+            cpp_params_cplx.constant_step_size = p.constant_step_size
+
+        cpp_constraints = \
+        <PyxConstraintGeneric**> \
+        PyMem_Malloc(sizeof(PyxConstraintGeneric*)*len(p.constraints))
+
+        for i in range(0,len(p.constraints)):
+            cons = p.constraints[i]
+            #print("FaustFact.fact_hierarchical() cons.name =", cons.name)
+            if(cons.is_int_constraint()):
+                #print("FaustFact.fact_hierarchical() Int Constraint")
+                cpp_constraints[i] = <PyxConstraintInt*> PyMem_Malloc(sizeof(PyxConstraintInt))
+                (<PyxConstraintInt*>cpp_constraints[i]).parameter = cons.param
+            elif(cons.is_real_constraint()):
+                #print("FaustFact.fact_hierarchical() Real Constraint")
+                cpp_constraints[i] = <PyxConstraintScalar[double]*> \
+                PyMem_Malloc(sizeof(PyxConstraintScalar[double]))
+                (<PyxConstraintScalar[double]*>cpp_constraints[i]).parameter =\
+                        cons.param
+            elif(cons.is_mat_constraint()):
+                #print("FaustFact.fact_hierarchical() Matrix Constraint")
+                if(isReal):
+                    cpp_constraints[i] = <PyxConstraintMat[double]*> \
+                            PyMem_Malloc(sizeof(PyxConstraintMat[double]))
+                    tmp_mat = cons.param
+                    (<PyxConstraintMat[double]*>cpp_constraints[i]).parameter =\
+                            &tmp_mat[0,0]
+                else:
+                    cpp_constraints[i] = <PyxConstraintMat[complex]*> \
+                            PyMem_Malloc(sizeof(PyxConstraintMat[complex]))
+                    tmp_mat_cplx = cons.param
+                    (<PyxConstraintMat[complex]*>cpp_constraints[i]).parameter =\
+                            &tmp_mat_cplx[0,0]
+
+            else:
+                raise ValueError("Constraint type/name is not recognized.")
+            cpp_constraints[i].name = cons.name
+            cpp_constraints[i].num_rows = cons.num_rows
+            cpp_constraints[i].num_cols = cons.num_cols
+
+        if(isReal):
+            cpp_params.constraints = cpp_constraints
+            cpp_params.num_rows = p.data_num_rows
+            cpp_params.num_cols = p.data_num_cols
+            cpp_params.num_constraints = len(p.constraints)
+        else:
+            cpp_params_cplx.constraints = cpp_constraints
+            cpp_params_cplx.num_rows = p.data_num_rows
+            cpp_params_cplx.num_cols = p.data_num_cols
+            cpp_params_cplx.num_constraints = len(p.constraints)
+
+        for i in range(0,p.num_facts):
+            if(isReal):
+                tmp_mat = p.init_facts[i]
+                cpp_params.init_facts[i] = &tmp_mat[0,0]
+                cpp_params.init_fact_sizes[i*2+0] = p.init_facts[i].shape[0]
+                cpp_params.init_fact_sizes[i*2+1] = p.init_facts[i].shape[1]
+            else:
+                tmp_mat_cplx = p.init_facts[i]
+                cpp_params_cplx.init_facts[i] = &tmp_mat_cplx[0,0]
+                cpp_params_cplx.init_fact_sizes[i*2+0] = p.init_facts[i].shape[0]
+                cpp_params_cplx.init_fact_sizes[i*2+1] = p.init_facts[i].shape[1]
+
+        core = FaustCore(core=True)
+        if(isReal):
+            core.core_faust_dbl = FaustCoreCy.fact_hierarchical(&Mview[0,0], M_num_rows, M_num_cols,
+ #           FaustCoreCy.fact_hierarchical(&Mview[0,0], M_num_rows, M_num_cols,
+                                      &cpp_params)
+            core._isReal = True
+            #TODO: FPP == complex not yet supported by C++ code
+#        else:
+#            core.core_faust_cplx = FaustCoreCy.fact_hierarchical(&Mview_cplx[0,0], M_num_rows, M_num_cols,
+#                                     &cpp_params_cplx)
+#            core._isReal = False
+        for i in range(0,len(p.constraints)):
+            PyMem_Free(cpp_constraints[i])
+        PyMem_Free(cpp_constraints)
+        if(isReal):
+            PyMem_Free(cpp_params.init_facts)
+            PyMem_Free(cpp_params.init_fact_sizes)
+        else:
+            PyMem_Free(cpp_params_cplx.init_facts)
+            PyMem_Free(cpp_params_cplx.init_fact_sizes)
+
+        PyMem_Free(cpp_stop_crits)
+
+        return core

wrapper/python/src/FaustFact.h

+#ifndef FAUST_FACT_H
+#define FAUST_FACT_H
+
+
+class PyxConstraintGeneric
+{
+
+    public:
+        int name;
+        unsigned long num_rows;
+        unsigned long num_cols;
+
+        bool is_int_constraint();
+        bool is_real_constraint();
+        bool is_mat_constraint();
+};
+
+class PyxConstraintInt : public PyxConstraintGeneric
+{
+    public:
+        unsigned long parameter;
+};
+
+template<typename FPP>
+class PyxConstraintScalar : public PyxConstraintGeneric
+{
+    public:
+        FPP parameter;
+};
+
+template<typename FPP>
+class PyxConstraintMat : public PyxConstraintGeneric
+{
+    public:
+        FPP* parameter;
+};
+
+template<typename FPP>
+class PyxStoppingCriterion
+{
+    public:
+        bool is_criterion_error;
+        int num_its;
+        FPP error_treshold;
+        unsigned long max_num_its;
+};
+
+template<typename FPP>
+class PyxParamsFact
+{
+    public:
+        int num_facts;
+        bool is_update_way_R2L;
+        FPP init_lambda;
+        FPP** init_facts;// num_facts elts
+        unsigned long* init_fact_sizes;
+        FPP step_size;
+        PyxConstraintGeneric** constraints; // size: (num_facts-1)*2
+        unsigned int num_constraints;
+        bool is_verbose;
+        bool constant_step_size;
+};
+
+template<typename FPP>
+class PyxParamsFactPalm4MSA : public PyxParamsFact<FPP>
+{
+    public:
+        PyxStoppingCriterion<FPP> stop_crit;
+};
+
+template<typename FPP>
+class PyxParamsHierarchicalFact : public PyxParamsFact<FPP>
+{
+    public:
+        unsigned int num_rows;
+        unsigned int num_cols;
+        PyxStoppingCriterion<FPP>* stop_crits; // must be of size 2
+        bool is_fact_side_left;
+};
+
+template<typename FPP>
+FaustCoreCpp<FPP>* fact_palm4MSA(FPP*, unsigned int, unsigned int, PyxParamsFactPalm4MSA<FPP>*, bool);
+
+template<typename FPP>
+FaustCoreCpp<FPP>* fact_hierarchical(FPP*, unsigned int, unsigned int, PyxParamsHierarchicalFact<FPP>*, bool);
+
+#include "FaustFact.hpp"
+
+#endif

wrapper/python/src/FaustFact.hpp

+#include <iostream>
+#include "faust_MatDense.h"
+#include "faust_Params.h"
+#include "faust_ParamsPalm.h"
+#include "faust_StoppingCriterion.h"
+#include "faust_Palm4MSA.h"
+#include "faust_HierarchicalFact.h"
+#include "faust_BlasHandle.h"
+#include "faust_ConstraintGeneric.h"
+
+using namespace std;
+
+bool PyxConstraintGeneric::is_int_constraint()
+{
+    switch(static_cast<faust_constraint_name>(this->name))
+    {
+        case CONSTRAINT_NAME_SP:
+        case CONSTRAINT_NAME_SPCOL:
+        case CONSTRAINT_NAME_SPLIN:
+        case CONSTRAINT_NAME_SP_POS:
+            return true;
+        default:
+            return false;
+    }
+}
+
+bool PyxConstraintGeneric::is_real_constraint()
+{
+    switch(static_cast<faust_constraint_name>(this->name))
+    {
+        case CONSTRAINT_NAME_NORMCOL:
+        case CONSTRAINT_NAME_NORMLIN:
+            return true;
+        default:
+            return false;
+    }
+}
+
+bool PyxConstraintGeneric::is_mat_constraint()
+{
+    switch(static_cast<faust_constraint_name>(this->name))
+    {
+        case CONSTRAINT_NAME_CONST:
+        case CONSTRAINT_NAME_SUPP:
+        case CONSTRAINT_NAME_BLKDIAG:
+            return true;
+        default:
+            return false;
+    }
+}
+
+template<typename FPP>
+void prepare_fact(const FPP* mat, const unsigned int num_rows, const unsigned int num_cols, const PyxParamsFact<FPP>* p,
+        /* out args */vector<const Faust::ConstraintGeneric<FPP,Cpu>*>& cons, vector<Faust::MatDense<FPP,Cpu> > &initFacts)
+{
+    Faust::ConstraintGeneric<FPP, Cpu>* tmp_cons;
+    if(p->is_verbose)
+    {
+        cout << "fact_palm4MSA() mat[0]= " << mat[0] << endl;
+        cout << "p->num_facts: " << p->num_facts << endl;
+        cout << "p->is_update_way_R2L: " << p->is_update_way_R2L << endl;
+        cout << "p->init_lambda: " << p->init_lambda << endl;
+        cout << "p->step_size: " << p->step_size << endl;
+        cout << "p->is_verbose: " << p->is_verbose << endl;
+        cout << "p->constant_step_size: " << p->constant_step_size << endl;
+    }
+    PyxConstraintInt* cons_int;
+    PyxConstraintScalar<FPP>* cons_real;
+    PyxConstraintMat<FPP>* cons_mat;
+    for(int i=0; i < p->num_constraints;i++)
+    {
+        if(p->is_verbose) {
+            cout << "constraint[" << i << "]->name: " << p->constraints[i]->name << endl;
+            cout << "constraint[" << i << "]->num_rows: " << p->constraints[i]->num_rows << endl;
+            cout << "constraint[" << i << "]->num_cols: " << p->constraints[i]->num_cols << endl;
+        }
+        //TODO: make ConstraintGeneric virtual and add a display() function to
+        //avoid this mess and also a function to convert to Faust::Constraint*
+        // corresponding object
+        if(p->constraints[i]->is_int_constraint())
+        {
+            cons_int = static_cast<PyxConstraintInt*>(p->constraints[i]);
+            if(p->is_verbose)
+                cout << "constraint[" << i << "]->parameter: " << cons_int->parameter << endl;
+            tmp_cons = new Faust::ConstraintInt<FPP,Cpu>(static_cast<faust_constraint_name>(p->constraints[i]->name), cons_int->parameter, p->constraints[i]->num_rows, p->constraints[i]->num_cols);
+            cons.push_back(tmp_cons);
+        }
+        else if(p->constraints[i]->is_real_constraint())
+        {
+            cons_real = static_cast<PyxConstraintScalar<FPP>*>(p->constraints[i]);
+            if(p->is_verbose)
+                cout << "constraint[" << i << "]->parameter: " << cons_real->parameter << endl;
+            tmp_cons = new Faust::ConstraintFPP<FPP,Cpu>(static_cast<faust_constraint_name>(p->constraints[i]->name), cons_real->parameter, p->constraints[i]->num_rows, p->constraints[i]->num_cols);
+            cons.push_back(tmp_cons);
+        }
+        else if(p->constraints[i]->is_mat_constraint())
+        {
+            cons_mat = static_cast<PyxConstraintMat<FPP>*>(p->constraints[i]);
+            if(p->is_verbose)
+                cout << "constraint[" << i << "]->parameter: " << cons_mat->parameter[0] << endl;
+            tmp_cons = new Faust::ConstraintMat<FPP,Cpu>(static_cast<faust_constraint_name>(p->constraints[i]->name), Faust::MatDense<FPP, Cpu>(cons_mat->parameter, p->constraints[i]->num_rows, p->constraints[i]->num_cols), p->constraints[i]->num_rows, p->constraints[i]->num_cols);
+            cons.push_back(tmp_cons);
+        }
+        else
+            handleError("FaustFact", "Invalid constraint.");
+    }
+    for(int i=0; i < p->num_facts;i++) {
+        if(p->is_verbose)
+        {
+            cout << "init_facts[" << i << "] ele[0][0]: " << p->init_facts[i][0] << " size: " << p->init_fact_sizes[i*2] << "x"<< p->init_fact_sizes[i*2+1] << endl; 
+        }
+        initFacts.push_back(Faust::MatDense<FPP, Cpu>(p->init_facts[i], p->init_fact_sizes[i*2], p->init_fact_sizes[i*2+1]));
+    }
+
+
+}
+
+template<typename FPP>
+FaustCoreCpp<FPP>* fact_palm4MSA(FPP* mat, unsigned int num_rows, unsigned int num_cols, PyxParamsFactPalm4MSA<FPP>* p)
+{
+    FaustCoreCpp<FPP>* core;
+    Faust::MatDense<FPP,Cpu> inMat(mat, num_rows, num_cols);
+    vector<const Faust::ConstraintGeneric<FPP,Cpu>*> cons;
+    vector<Faust::MatDense<FPP,Cpu> > initFacts;
+    if(p->is_verbose) {
+        cout << "stop_crit.is_criterion_error: " << p->stop_crit.is_criterion_error << endl;
+        cout << "stop_crit.error_treshold: " << p->stop_crit.error_treshold << endl;
+        cout << "stop_crit.num_its: " << p->stop_crit.num_its << endl;
+        cout << "stop_crit.max_num_its: " << p->stop_crit.max_num_its << endl;
+    }
+    prepare_fact(mat, num_rows, num_cols, p, cons, initFacts);
+
+    // set all constructor arguments because they could be at non-default
+    // values
+    Faust::StoppingCriterion<FPP> crit(p->stop_crit.num_its, p->stop_crit.is_criterion_error, p->stop_crit.error_treshold, p->stop_crit.max_num_its);
+
+    Faust::ParamsPalm<FPP,Cpu> params(inMat, p->num_facts, cons, initFacts, crit, p->is_verbose, p->is_update_way_R2L, p->init_lambda, p->constant_step_size, p->step_size);
+
+    if(p->is_verbose) params.Display();
+
+    Faust::BlasHandle<Cpu> blasHandle;
+
+    Faust::Palm4MSA<FPP,Cpu> palm(params,blasHandle,true);
+
+    palm.compute_facts();
+
+
+//    palm.get_facts(faust);
+//
+    const std::vector<Faust::MatDense<FPP,Cpu> >& full_facts = palm.get_facts();
+    Faust::Transform<FPP, Cpu> faust(full_facts, true);
+    if(p->is_verbose) faust.Display();
+
+    Faust::TransformHelper<FPP, Cpu> th(faust);
+    core = new FaustCoreCpp<FPP>(th);
+
+    blasHandle.Destroy();
+
+    for (typename std::vector<const Faust::ConstraintGeneric<FPP,Cpu>*>::iterator it = cons.begin() ; it != cons.end(); ++it)
+        delete *it;
+
+
+    return core;
+
+}
+
+template<typename FPP>
+FaustCoreCpp<FPP>* fact_hierarchical(FPP* mat, unsigned int num_rows, unsigned int num_cols, PyxParamsHierarchicalFact<FPP>* p)
+{
+    FaustCoreCpp<FPP>* core;
+    Faust::MatDense<FPP,Cpu> inMat(mat, num_rows, num_cols);
+    vector<const Faust::ConstraintGeneric<FPP,Cpu>*> cons;
+    vector<std::vector<const Faust::ConstraintGeneric<FPP,Cpu>*>> cons_;
+    vector<const Faust::ConstraintGeneric<FPP,Cpu>*> fact_cons;
+    vector<const Faust::ConstraintGeneric<FPP,Cpu>*> residuum_cons;
+    vector<Faust::MatDense<FPP,Cpu> > initFacts;
+    if(p->is_verbose)
+    {
+        cout << "p->num_rows: " << p->num_rows << endl;
+        cout << "p->num_cols: " << p->num_cols << endl;
+        cout << "p->is_fact_side_left: " << p->is_fact_side_left;
+        cout << "stop_crits[0].is_criterion_error: " << p->stop_crits[0].is_criterion_error << endl;
+        cout << "stop_crits[0].error_treshold: " << p->stop_crits[0].error_treshold << endl;
+        cout << "stop_crits[0].num_its: " << p->stop_crits[0].num_its << endl;
+        cout << "stop_crits[0].max_num_its: " << p->stop_crits[0].max_num_its << endl;
+        cout << "stop_crits[1].is_criterion_error: " << p->stop_crits[1].is_criterion_error << endl;
+        cout << "stop_crits[1].error_treshold: " << p->stop_crits[1].error_treshold << endl;
+        cout << "stop_crits[1].num_its: " << p->stop_crits[1].num_its << endl;
+        cout << "stop_crits[1].max_num_its: " << p->stop_crits[1].max_num_its << endl;
+    }
+    prepare_fact(mat, num_rows, num_cols, p, cons, initFacts);
+
+    // set all constructor arguments because they could be at non-default
+    // values
+    Faust::StoppingCriterion<FPP> crit0(p->stop_crits[0].num_its, p->stop_crits[0].is_criterion_error, p->stop_crits[0].error_treshold, p->stop_crits[0].max_num_its); //2 facts
+    Faust::StoppingCriterion<FPP> crit1(p->stop_crits[1].num_its, p->stop_crits[1].is_criterion_error, p->stop_crits[1].error_treshold, p->stop_crits[1].max_num_its); //global
+
+    for(int i=0;i<p->num_facts-1;i++)
+        fact_cons.push_back(cons[i]);
+
+    for(int i=0;i<p->num_facts-1;i++)
+        residuum_cons.push_back(cons[i+p->num_facts-1]);
+
+    cons_.push_back(fact_cons);
+    cons_.push_back(residuum_cons);
+
+    Faust::Params<FPP,Cpu> params(p->num_rows, p->num_cols, p->num_facts, cons_, initFacts, crit0, crit1, p->is_verbose, p->is_update_way_R2L, p->is_fact_side_left, p->init_lambda, p->constant_step_size, p->step_size);
+    
+    if(p->is_verbose) params.Display();
+
+
+    Faust::BlasHandle<Cpu> blasHandle;
+    Faust::SpBlasHandle<Cpu> spblasHandle;
+
+    Faust::HierarchicalFact<FPP,Cpu> hierFact(inMat, params,blasHandle,spblasHandle);
+
+    hierFact.compute_facts();
+
+    vector<Faust::MatSparse<FPP,Cpu> > facts;
+    hierFact.get_facts(facts);
+    FPP lambda = hierFact.get_lambda();
+    (facts[0]) *= lambda;
+    // transform the sparse matrix into generic one
+    std::vector<Faust::MatGeneric<FPP,Cpu> *> list_fact_generic;
+    list_fact_generic.resize(facts.size());
+    for (int i=0;i<list_fact_generic.size();i++)
+        list_fact_generic[i]=facts[i].Clone();
+
+    Faust::Transform<FPP,Cpu> faust(list_fact_generic);
+
+    for (int i=0;i<list_fact_generic.size();i++)
+        delete list_fact_generic[i];
+
+    if(p->is_verbose) faust.Display();
+
+    Faust::TransformHelper<FPP, Cpu> th(faust);
+    core = new FaustCoreCpp<FPP>(th);
+
+    blasHandle.Destroy();
+    spblasHandle.Destroy();
+
+    for (typename std::vector<const Faust::ConstraintGeneric<FPP,Cpu>*>::iterator it = cons.begin() ; it != cons.end(); ++it)
+        delete *it;
+
+    return core;
+
+}