Mentions légales du service

Skip to content
Snippets Groups Projects
Commit 5ec8d880 authored by hhakim's avatar hhakim
Browse files

Fix correctness of poly module functions and display of Fausts they produce. 

Major bug is a regression introduced by e4cbb0c3.
parent f8a6c144
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/faust_linear_operator/CPU/faust_TransformHelperPoly.h
+ 11
− 10
View file @ 5ec8d880
......@@ -67,6 +67,7 @@ namespace Faust
bool on_gpu=false);
TransformHelperPoly(uint K, const TransformHelperPoly<FPP>& src);
TransformHelperPoly() : TransformHelper<FPP, Cpu>(), mul_and_combi_lin_on_gpu(false), T0_is_arbitrary(false), L(nullptr), rR(nullptr) {}
public:
faust_unsigned_int getNbRow() const;
......@@ -98,7 +99,7 @@ namespace Faust
TransformHelper<FPP,Cpu>* clone();
void update_total_nnz();
MatDense<FPP, Cpu> multiply(const MatSparse<FPP,Cpu> &A, const bool transpose=false, const bool conjugate=false);
MatDense<FPP, Cpu> multiply(const MatSparse<FPP,Cpu> &A);
TransformHelper<FPP, Cpu>* multiply(const TransformHelper<FPP, Cpu>*) const;
TransformHelper<FPP, Cpu>* multiply(FPP& scalar);
faust_unsigned_int getNBytes() const;
......@@ -125,15 +126,15 @@ namespace Faust
TransformHelper<FPP,Cpu>* swap_rows(const faust_unsigned_int id1, const faust_unsigned_int id2, const bool permutation=false, const bool inplace=false, const bool check_transpose=true);
TransformHelper<FPP, Cpu>* slice(faust_unsigned_int start_row_id, faust_unsigned_int end_row_id,
faust_unsigned_int start_col_id, faust_unsigned_int end_col_id);
Vect<FPP,Cpu> multiply(const Vect<FPP,Cpu> &x, const bool transpose=false, const bool conjugate=false);
Vect<FPP,Cpu> multiply(const FPP* x, const bool transpose=false, const bool conjugate=false);
void multiply(const FPP* x, FPP* y, const bool transpose=false, const bool conjugate=false);
void multiply_cpu(const FPP* x, FPP* y, const bool transpose=false, const bool conjugate=false);
void multiply_gpu(const FPP* x, FPP* y, const bool transpose=false, const bool conjugate=false);
MatDense<FPP, Cpu> multiply(const MatDense<FPP,Cpu> &X, const bool transpose=false, const bool conjugate=false);
void multiply(const FPP* X, int n, FPP* out, const bool transpose=false, const bool conjugate=false);
void multiply_cpu(const FPP* X, int n, FPP* out, const bool transpose=false, const bool conjugate=false);
void multiply_gpu(const FPP* X, int n, FPP* out, const bool transpose=false, const bool conjugate=false);
Vect<FPP,Cpu> multiply(const Vect<FPP,Cpu> &x);
Vect<FPP,Cpu> multiply(const FPP* x);
void multiply(const FPP* x, FPP* y);
void multiply_cpu(const FPP* x, FPP* y);
void multiply_gpu(const FPP* x, FPP* y);
MatDense<FPP, Cpu> multiply(const MatDense<FPP,Cpu> &X);
void multiply(const FPP* X, int n, FPP* out);
void multiply_cpu(const FPP* X, int n, FPP* out);
void multiply_gpu(const FPP* X, int n, FPP* out);
TransformHelper<FPP, Cpu>* next(uint K);
TransformHelper<FPP, Cpu>* next();
Vect<FPP, Cpu> poly(MatDense<FPP,Cpu> & basisX, Vect<FPP, Cpu> coeffs);
......
src/faust_linear_operator/CPU/faust_TransformHelperPoly.hpp
+ 26
− 26
View file @ 5ec8d880
......@@ -11,7 +11,7 @@ namespace Faust
MatSparse<FPP, Cpu>* rR /* = nullptr*/,
MatSparse<FPP,Cpu> *T0 /*= nullptr*/,
BasisLaziness laziness /*= INSTANTIATE_COMPUTE_AND_FREE */,
bool on_gpu /*= false*/) : TransformHelper<FPP,Cpu>(), T0_is_arbitrary(false)
bool on_gpu /*= false*/) : TransformHelperPoly()
{
// assuming L is symmetric
this->L = L;
......@@ -33,7 +33,6 @@ namespace Faust
if(T0 != nullptr)
{
T0->Display();
// T0 is arbitrary and need to be stored
// best choice is to instantiate the factor directly
this->basisChebyshevT0(T0);
......@@ -48,7 +47,7 @@ namespace Faust
}
template<typename FPP>
TransformHelperPoly<FPP>::TransformHelperPoly(uint K, const TransformHelperPoly<FPP>& src) : TransformHelper<FPP, Cpu>()
TransformHelperPoly<FPP>::TransformHelperPoly(uint K, const TransformHelperPoly<FPP>& src) : TransformHelperPoly()
{
if(K+1 < src.size())
throw std::runtime_error("The src TransformHelperPoly size can't be greater than K+1.");
......@@ -113,23 +112,23 @@ namespace Faust
}
template<typename FPP>
Vect<FPP,Cpu> TransformHelperPoly<FPP>::multiply(const Vect<FPP,Cpu> &x, const bool transpose/*=false*/, const bool conjugate/*=false*/)
Vect<FPP,Cpu> TransformHelperPoly<FPP>::multiply(const Vect<FPP,Cpu> &x)
{
return std::move(this->multiply(x.getData(), transpose, conjugate));
return std::move(this->multiply(x.getData()));
}
template<typename FPP>
Vect<FPP,Cpu> TransformHelperPoly<FPP>::multiply(const FPP* x, const bool transpose/*=false*/, const bool conjugate/*=false*/)
Vect<FPP,Cpu> TransformHelperPoly<FPP>::multiply(const FPP* x)
{
int d = L->getNbRow();
uint K = this->size()-1;
Vect<FPP, Cpu> v0(d*(K+1));
multiply(x, v0.getData(), transpose, conjugate);
multiply(x, v0.getData());
return std::move(v0);
}
template<typename FPP>
void TransformHelperPoly<FPP>::multiply_cpu(const FPP* x, FPP* y, const bool transpose/*=false*/, const bool conjugate/*=false*/)
void TransformHelperPoly<FPP>::multiply_cpu(const FPP* x, FPP* y)
{
/**
* Recurrence relation (k=1 to K):
......@@ -166,7 +165,7 @@ namespace Faust
}
template<typename FPP>
void TransformHelperPoly<FPP>::multiply_gpu(const FPP* x, FPP* y, const bool transpose/*=false*/, const bool conjugate/*=false*/)
void TransformHelperPoly<FPP>::multiply_gpu(const FPP* x, FPP* y)
{
#ifdef USE_GPU_MOD
int d = L->getNbRow();
......@@ -202,15 +201,15 @@ namespace Faust
}
template<typename FPP>
void TransformHelperPoly<FPP>::multiply(const FPP* x, FPP* y, const bool transpose/*=false*/, const bool conjugate/*=false*/)
void TransformHelperPoly<FPP>::multiply(const FPP* x, FPP* y)
{
if(this->mul_and_combi_lin_on_gpu)
{
multiply_gpu(x, y, transpose, conjugate);
multiply_gpu(x, y);
}
else
{
multiply_cpu(x, y, transpose, conjugate);
multiply_cpu(x, y);
}
}
......@@ -352,32 +351,32 @@ namespace Faust
// }
template<typename FPP>
MatDense<FPP, Cpu> TransformHelperPoly<FPP>::multiply(const MatDense<FPP,Cpu> &X, const bool transpose/*=false*/, const bool conjugate/*=false*/)
MatDense<FPP, Cpu> TransformHelperPoly<FPP>::multiply(const MatDense<FPP,Cpu> &X)
{
// std::cout << "TransformHelperPoly<FPP>::multiply(MatDense)" << std::endl;
int d = L->getNbRow();
uint K = this->size()-1;
int n = X.getNbCol();
MatDense<FPP,Cpu> Y(d*(K+1), n);
multiply(X.getData(), n, Y.getData(), transpose, conjugate);
multiply(X.getData(), n, Y.getData());
return Y;
}
template<typename FPP>
void TransformHelperPoly<FPP>::multiply(const FPP* X, int n, FPP* Y, const bool transpose/*=false*/, const bool conjugate/*=false*/)
void TransformHelperPoly<FPP>::multiply(const FPP* X, int n, FPP* Y)
{
if(this->mul_and_combi_lin_on_gpu)
{
multiply_gpu(X, n, Y, transpose, conjugate);
multiply_gpu(X, n, Y);
}
else
{
multiply_cpu(X, n, Y, transpose, conjugate);
multiply_cpu(X, n, Y);
}
}
template<typename FPP>
void TransformHelperPoly<FPP>::multiply_cpu(const FPP* X, int n, FPP* Y, const bool transpose/*=false*/, const bool conjugate/*=false*/)
void TransformHelperPoly<FPP>::multiply_cpu(const FPP* X, int n, FPP* Y)
{
int d = L->getNbRow();
uint K = this->size()-1;
......@@ -390,12 +389,12 @@ namespace Faust
// memcpy(V0_ord.getData()+scale*i, y.getData(), sizeof(FPP)*scale);
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>> x_vec(const_cast<FPP*>(X+i*d), d);
Eigen::Map<Eigen::Matrix<FPP, Eigen::Dynamic, 1>> y_vec(const_cast<FPP*>(Y+i*scale), scale);
multiply(x_vec.data(), y_vec.data(), transpose, conjugate);
multiply(x_vec.data(), y_vec.data());
}
}
template<typename FPP>
void TransformHelperPoly<FPP>::multiply_gpu(const FPP* X, int n, FPP* Y, const bool transpose/*=false*/, const bool conjugate/*=false*/)
void TransformHelperPoly<FPP>::multiply_gpu(const FPP* X, int n, FPP* Y)
{
#ifdef USE_GPU_MOD
int d = L->getNbRow();
......@@ -441,7 +440,7 @@ namespace Faust
}
template<typename FPP>
MatDense<FPP, Cpu> TransformHelperPoly<FPP>::multiply(const MatSparse<FPP,Cpu> &A, const bool transpose/*=false*/, const bool conjugate/*=false*/)
MatDense<FPP, Cpu> TransformHelperPoly<FPP>::multiply(const MatSparse<FPP,Cpu> &A)
{
MatDense<FPP, Cpu> M(this->getNbRow(), A.getNbCol());
M.setZeros();
......@@ -459,7 +458,7 @@ namespace Faust
{
auto id = *i;
auto vect = A.get_col(id);
this->multiply(vect.getData(), M.getData()+M.getNbRow()*id, transpose, conjugate);
this->multiply(vect.getData(), M.getData()+M.getNbRow()*id);
}
return M;
}
......@@ -1080,7 +1079,7 @@ namespace Faust
str<<"Faust size ";
str << this->get_fact_nb_rows(0) << "x" << this->get_fact_nb_cols(this->size()-1);
auto density = (double)this->get_total_nnz()/this->getNbRow()/this->getNbCol();
str <<", density "<< density << ", nnz_sum "<<this->get_total_nnz() << ", " << this->size() << " factor(s): "<< std::endl;
str <<", density "<< density << ", nnz_sum "<<this->get_total_nnz() << ", " << this->size() << " factor(s):"<< std::endl;
for (int i=0 ; i<this->size() ; i++)
{
str << "- ";
......@@ -1099,9 +1098,10 @@ namespace Faust
else
str << nrows << "x" << ncols;
str << ", density "<< density <<", nnz "<< this->get_fact_nnz(i);
str <<std::endl;
if (! T0_is_arbitrary || is_fact_created[this->size()-1] && this->get_gen_fact_nonconst(this->size()-1)->is_id())
str <<" identity matrix flag" << std::endl;
if(i < this->size() - 1)
str <<std::endl;
else if (! T0_is_arbitrary || is_fact_created[this->size()-1] && this->get_gen_fact_nonconst(this->size()-1)->is_id())
str << std::endl << " identity matrix flag" /*<< std::endl*/;
}
return str.str();
}
......
wrapper/python/pyfaust/poly.py
+ 44
− 41
View file @ 5ec8d880
......@@ -82,8 +82,7 @@ def Chebyshev(L, K, dev='cpu', T0=None, impl='native'):
def basis(L, K, basis_name, dev='cpu', T0=None, **kwargs):
"""
Builds the Faust of the polynomial basis defined on the sparse matrix L.
"""Builds the Faust of the polynomial basis defined on the sparse matrix L.
Args:
L: the sparse scipy square matrix in CSR format (scipy.sparse.csr_matrix).
......@@ -100,28 +99,33 @@ def basis(L, K, basis_name, dev='cpu', T0=None, **kwargs):
Example:
>>> from pyfaust.poly import basis
>>> from scipy.sparse import random
>>> from numpy.random import seed
>>> seed(42) # for reproducibility
>>> L = random(50, 50, .02, format='csr')
>>> L = L@L.T
>>> K = 3
>>> F = basis(L, K, 'chebyshev')
>>> F
Faust size 200x50, density 0.0687, nnz_sum 687, 4 factor(s):
- FACTOR 0 (real) SPARSE, size 200x150, density 0.0093, nnz 279
- FACTOR 1 (real) SPARSE, size 150x100, density 0.0152667, nnz 229
- FACTOR 2 (real) SPARSE, size 100x50, density 0.0258, nnz 129
- FACTOR 3 (real) SPARSE, size 50x50, density 0.02, nnz 50
Faust size 200x50, density 0.0699, nnz_sum 699, 4 factor(s):
- FACTOR 0 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 1 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 2 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 3 (double) SPARSE, size 50x50, density 0.02, nnz 50
identity matrix flag
Generate the next basis (the one with one additional dimension,
whose the polynomial greatest degree is K+1 = 4)
>>> G = next(F)
>>> G
Faust size 250x50, density 0.08128, nnz_sum 1016, 5 factor(s):
- FACTOR 0 (real) SPARSE, size 250x200, density 0.00658, nnz 329
- FACTOR 1 (real) SPARSE, size 200x150, density 0.0093, nnz 279
- FACTOR 2 (real) SPARSE, size 150x100, density 0.0152667, nnz 229
- FACTOR 3 (real) SPARSE, size 100x50, density 0.0258, nnz 129
- FACTOR 4 (real) SPARSE, size 50x50, density 0.02, nnz 50
Faust size 250x50, density 0.08256, nnz_sum 1032, 5 factor(s):
- FACTOR 0 (double) SPARSE, size 250x200, density 0.00666, nnz 333
- FACTOR 1 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 2 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 3 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 4 (double) SPARSE, size 50x50, density 0.02, nnz 50
identity matrix flag
The factors 0 to 3 of G are views of the same factors of F.
......@@ -132,13 +136,11 @@ def basis(L, K, basis_name, dev='cpu', T0=None, **kwargs):
>>> F2 = basis(L, K, 'chebyshev', T0=random(50,2, .3, format='csr'))
>>> F2
Faust size 200x2, density 1.7125, nnz_sum 685, 4 factor(s):
- FACTOR 0 (real) SPARSE, size 200x150, density 0.0095, nnz 285
- FACTOR 1 (real) SPARSE, size 150x100, density 0.0156667, nnz 235
- FACTOR 2 (real) SPARSE, size 100x50, density 0.027, nnz 135
- FACTOR 3 (real) SPARSE, size 50x2, density 0.3, nnz 30
Faust size 200x2, density 1.7475, nnz_sum 699, 4 factor(s):
- FACTOR 0 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 1 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 2 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 3 (double) SPARSE, size 50x2, density 0.5, nnz 50
"""
# impl (optional): 'native' (by default) for the C++ impl., "py" for the Python impl.
......@@ -188,17 +190,19 @@ def poly(coeffs, basis='chebyshev', L=None, X=None, dev='cpu', out=None,
>>> import numpy as np
>>> from pyfaust.poly import basis, poly
>>> from scipy.sparse import random
>>> np.random.seed(42) # for reproducibility
>>> L = random(50, 50, .02, format='csr')
>>> L = L@L.T
>>> coeffs = np.array([.5, 1, 2, 3])
>>> G = poly(coeffs, 'chebyshev', L)
>>> G
Faust size 50x50, density 0.3608, nnz_sum 902, 5 factor(s):
- FACTOR 0 (real) SPARSE, size 50x200, density 0.02, nnz 200
- FACTOR 1 (real) SPARSE, size 200x150, density 0.00946667, nnz 284
- FACTOR 2 (real) SPARSE, size 150x100, density 0.0156, nnz 234
- FACTOR 3 (real) SPARSE, size 100x50, density 0.0268, nnz 134
- FACTOR 4 (real) SPARSE, size 50x50, density 0.02, nnz 50
Faust size 50x50, density 0.3596, nnz_sum 899, 5 factor(s):
- FACTOR 0 (double) SPARSE, size 50x200, density 0.02, nnz 200
- FACTOR 1 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 2 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 3 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 4 (double) SPARSE, size 50x50, density 0.02, nnz 50
identity matrix flag
Which is equivalent to do as below (in two times):
......@@ -207,18 +211,19 @@ def poly(coeffs, basis='chebyshev', L=None, X=None, dev='cpu', out=None,
>>> coeffs = np.array([.5, 1, 2, 3])
>>> G = poly(coeffs, F)
>>> G
Faust size 50x50, density 0.3608, nnz_sum 902, 5 factor(s):
- FACTOR 0 (real) SPARSE, size 50x200, density 0.02, nnz 200
- FACTOR 1 (real) SPARSE, size 200x150, density 0.00946667, nnz 284
- FACTOR 2 (real) SPARSE, size 150x100, density 0.0156, nnz 234
- FACTOR 3 (real) SPARSE, size 100x50, density 0.0268, nnz 134
- FACTOR 4 (real) SPARSE, size 50x50, density 0.02, nnz 50
Faust size 50x50, density 0.3596, nnz_sum 899, 5 factor(s):
- FACTOR 0 (double) SPARSE, size 50x200, density 0.02, nnz 200
- FACTOR 1 (double) SPARSE, size 200x150, density 0.00943333, nnz 283
- FACTOR 2 (double) SPARSE, size 150x100, density 0.0155333, nnz 233
- FACTOR 3 (double) SPARSE, size 100x50, density 0.0266, nnz 133
- FACTOR 4 (double) SPARSE, size 50x50, density 0.02, nnz 50
identity matrix flag
Above G is a Faust because F is too.
Below the full array of the Faust F is passed, so an array is returned into GA.
>>> GA = poly(coeffs, F.toarray())
>>> type(GA)
numpy.ndarray
<class 'numpy.ndarray'>
But of course they are equal:
......@@ -561,19 +566,16 @@ def expm_multiply(A, B, t, K=10, tradeoff='time', dev='cpu', **kwargs):
>>> import numpy as np
>>> from scipy.sparse import random
>>> from pyfaust.poly import expm_multiply as fexpm_multiply
>>> np.random.seed(42) # for reproducibility
>>> L = random(5, 5, .2, format='csr')
>>> L = L@L.T
>>> x = np.random.rand(L.shape[1])
>>> t = np.linspace(start=-.5, stop=-0.1, num=3, endpoint=True)
>>> y = fexpm_multiply(L, x, t)
>>> y
array([[ 0.20063382, 0.39176039, 0.62490929, 0.60165209,
-0.00082166],
[ 0.20063382, 0.44945087, 0.62490929, 0.6401542 ,
0.02325689],
[ 0.20063382, 0.51456348, 0.62490929, 0.68279266,
0.05458717]])
array([[0.12706927, 0.21111065, 0.36636184, 0.41736344, 0.78517596],
[0.13190047, 0.24040598, 0.36636184, 0.4324354 , 0.78517596],
[0.13691536, 0.27376655, 0.36636184, 0.44805164, 0.78517596]])
"""
if not isinstance(A, csr_matrix):
......@@ -679,13 +681,14 @@ def invm_multiply(A, B, rel_err=1e-6, tradeoff='time', max_K=np.inf, dev='cpu',
>>> from scipy.sparse import random
>>> from pyfaust.poly import invm_multiply
>>> from numpy.linalg import norm, inv
>>> np.random.seed(42) # for reproducibility
>>> A = random(64, 64, .1, format='csr')
>>> A = A@A.T
>>> B = np.random.rand(A.shape[1],2)
>>> A_inv_B = invm_multiply(A, B, rel_err=1e-2, max_K=2048)
>>> A_inv_B_ref = inv(A.toarray())@B
>>> print("err:", norm(A_inv_B-A_inv_B_ref)/norm(A_inv_B))
err: 0.011045280586803805
err: 0.027283169722939017
"""
eps = rel_err
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment