Mentions légales du service

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

Minor changes in pyfaust.tools.omp.

parent 449afc4e
Branches
Tags
No related merge requests found
...@@ -5,12 +5,13 @@ ...@@ -5,12 +5,13 @@
import numpy as np import numpy as np
from scipy.sparse.linalg import spsolve_triangular from scipy.sparse.linalg import spsolve_triangular
from numpy.linalg import solve, lstsq, norm from numpy.linalg import solve, lstsq
from scipy.sparse import hstack, vstack, csr_matrix from scipy.sparse import hstack, vstack, csr_matrix
from numpy import concatenate as cat from numpy import concatenate as cat
from numpy import zeros, argmax, empty, ndarray from numpy import zeros, argmax, empty, ndarray
from pyfaust import Faust from pyfaust import Faust
def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False): def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False):
""" """
Runs the greedy OMP algorithm optimized by Cholesky decomposition. Runs the greedy OMP algorithm optimized by Cholesky decomposition.
...@@ -47,8 +48,8 @@ def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False): ...@@ -47,8 +48,8 @@ def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False):
else: else:
raise Exception("y must be a vector") raise Exception("y must be a vector")
if(Faust.isFaust(D) or isinstance(D, ndarray)): if(Faust.isFaust(D) or isinstance(D, ndarray)):
P = lambda z : D@z P = lambda z: D@z
Pt = lambda z : D.T.conj()@z Pt = lambda z: D.T.conj()@z
else: else:
raise Exception("D must be a Faust or a numpy.ndarray. Here D is " raise Exception("D must be a Faust or a numpy.ndarray. Here D is "
"a:"+str(type(D))) "a:"+str(type(D)))
...@@ -80,15 +81,15 @@ def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False): ...@@ -80,15 +81,15 @@ def omp(y, D, maxiter=None, tol=0, relerr=True, verbose=False):
s = s_initial s = s_initial
R = empty((maxiter+1,maxiter+1)).astype(np.complex) R = empty((maxiter+1,maxiter+1)).astype(np.complex)
oldErr = y.T.conj()@y oldErr = y.T.conj()@y
err_mse = [] # err_mse = []
t=0 t = 0
IN = [] IN = []
DR = Pt(residual) DR = Pt(residual)
done = False done = False
it = 1 it = 1
MAXITER=n MAXITER = n
while not done: while not done:
...@@ -231,3 +232,6 @@ def UpdateCholeskySparse(R,P,Pt,index,m): ...@@ -231,3 +232,6 @@ def UpdateCholeskySparse(R,P,Pt,index,m):
return R return R
greed_omp_chol = omp greed_omp_chol = omp
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment