jro_soft/compressed_sensing Commit - r20:21

Fixed erros in sfb, idualtree, irls_dn, and irls_dn2

yamaro -

r20:21

parent child

trunk/src/src/RICS_paper_compare_noise.py +18 -18

                     chi2 = np.sum((es/sigma)**2);
-                    # CS inversion using Iteratively Reweighted Least Squares (IRLS)-------------
+                    # --------- CS inversion using Iteratively Reweighted Least Squares (IRLS) -------------
                     # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
-                    Psi = deb4_basis(Nr); ###### REPLACED BY LINEs BELOW (?)
+                    Psi = deb4_basis(Nr);
-                    print 'FINALLY!'
+                    # REMOVE THIS--------------------------------
-                    print Psi.shape
-                    # REMOVE THIS?--------------------------------
                     #wavelet1 = pywt.Wavelet('db4')
                     #Phi, Psi, x = wavelet1.wavefun(level=3)
                     # --------------------------------------------
                     # add "sum to 1" constraint
-                    #         H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
+                    H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
-                    #         N_temp = np.array([[Nr/Nt]]);
+                    N_temp = np.array([[Nr/Nt]]);
-                    #         g2 = np.concatenate( (gnz, N_temp), axis=0 );
+                    g2 = np.concatenate( (gnz, N_temp), axis=0 );
-                    #         H2 = H2.T.conj();
+                    #H2 = H2.T.conj();
-                    #         print 'H2 shape', H2.shape
-                    #         print 'Psi shape', Psi.shape
+                    #Psi = Psi.T.conj(); # to align matrices
-                    #         s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
+                    ####print 'H2 shape', H2.shape
+                    #####print 'Psi shape', Psi.shape
+                    A = np.dot(H2,Psi);
+                    s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
                     #         f_cs = Psi*s;
                     #
                     #         # plot
                     # # # Scaling and shifting
                     # # # Only necessary for capon solution
                     f_capon = f_capon/np.max(f_capon)*np.max(fact);

trunk/src/src/cshift.py 0 -2

                 N = x.size;
                 n = np.arange(N);
                 n = np.mod(n-m, N);
-                print x.shape
                 y = x[0,n];

trunk/src/src/deb4_basis.py +2 -4

             def deb4_basis(N):
                 Psi = np.zeros(shape=(N,2*N+1));
-                idx = 1;
+                idx = 0;
                 J = 4;
                 [Faf, Fsf] = FSfarras();
                 [af, sf] = dualfilt1();
                 # compute transform of zero vector
                 x = np.zeros(shape=(1,N));
                 w = dualtree(x, J, Faf, af);
                 # Uses both real and imaginary wavelets
                 for i in range (0, J):
                             idx = idx + 1;
                 # Add uniform vector (seems to be useful if there's a background
-                Psi[:,2*N+1] = 1/np.sqrt(N);
+                Psi[:,2*N] = 1/np.sqrt(N);
                 return Psi

trunk/src/src/idualtree.py +1 -1

                 y2 = w[J][1];
                 for j in range (J-1, 0, -1):
-                    y2 = sfb(y2, w[j][2], sf[0,1]);
+                    y2 = sfb(y2, w[j][1], sf[0,1]);
                 y2 = sfb(y2, w[0][1], Fsf[0,1]);

trunk/src/src/irls_dn.py +16 -13

                 # and http://web.eecs.umich.edu/~aey/sparse/sparse11.pdf)
                 # Note to self:  I found that "warm-starting" didn't really help too much.
                 [M,N] = A.shape;
                 # Initialize and precompute:
                 eps = 1e-2; # damping parameter
-                [Q,R] = linalg.qr(A.T.conj(),0);
-                print A.shape
+                [Q,R] = linalg.qr(A.T.conj(), mode='economic');
-                print R.shape
-                print b.shape
                 c = linalg.solve(R.T.conj(),b); # will be used later also
-                u = Q*c; # minimum 2-norm solution
+                u = np.dot(Q,c); # minimum 2-norm solution
                 I = sps.eye(M);
                 #---------- not needed, defined above--------------
                         # Empty temporary N x N matrix
                         temp = np.zeros(shape=(N,N))
+                        k = 0
                         # Sparse matrix
-                        for i in range (1, N):
+                        for i in range (0, N):
-                            for j in range (1,N):
+                            for j in range (0,N):
                                 if(i==j):
-                                    temp[i,j] = w
+                                    temp[i,j] = w[k]
+                                    k = k+1
                         # Compressed Sparse Matrix
                         W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
                         WAT = W*A.T.conj();
-                        u_new = WAT * ( linalg.solve (A*WAT + lambda1*I), b);
+                        u_new = np.dot(WAT , linalg.solve(np.dot(A,WAT) + np.dot(lambda1,I), b));
                         # See if this subproblem is converging
                         delu = np.linalg.norm(u_new-u)/np.linalg.norm(u);
-                        epschange = delu < (np.sqrt(eps)/100);
+                        epschange = delu < (np.sqrt(eps)/100.0);
                         # Make update
                         u = u_new;
-                    eps = eps/10; # decrease eps
+                    eps = eps/10.0; # decrease eps
                     # Print info
-                    print 'eps =',eps;
+                    #print 'eps =',eps;
                 return u

trunk/src/src/irls_dn2.py +28 -7

             from irls_dn import *
             from scipy.optimize import fsolve
+            import numpy as np
+            from scipy.optimize import root
             def irls_dn2(A,b,p,G):
                 # Line Search
-                alpha = 2; # Line search parameter
+                alpha = 2.0; # Line search parameter
                 lambda1 = 1e5; # What's a reasonable but safe initial guess?
                 u = irls_dn(A,b,p,lambda1);
-                fid = np.norm(A*u-b)**2;
+                fid = np.linalg.norm(np.dot(A,u)-b)**2;
                 print '----------------------------------\n';
                 while (fid >= G):
                     lambda1 = lambda1 / alpha; # Balance between speed and accuracy
                     u = irls_dn(A,b,p,lambda1);
-                    fid = np.norm(A*u-b)**2;
+                    fid = np.linalg.norm(np.dot(A,u)-b)**2;
-                    print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid;
+                    print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
                 # Refinement using fzero
                 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
-                f = lambda lambda1: np.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G;
+                def myfun(lambda1):
+                    print "A = ", A.shape
+                    print "b = ", b.shape
+                    lambda1
+                    return np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G;
+                #f = lambda lambda1: np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G; NOOOOOO
             #    opts = optimset('fzero');
             #     opts.Display = 'none';
             #     opts.TolX = 0.01*lambda1;
-                lambda1 = fsolve(f,lambda0);  # FALTA OPTIMIZE ESTO
+                #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
+                print "tolerancia=", 0.01*lambda1
+                #lambda1 = root(myfun,lambda0, method='krylov',  tol=0.01*lambda1);
+                print "lamda1=", lambda1
+                print "lambda0=", lambda0
+                lambda1 = fsolve(myfun,lambda0);  # FALTA OPTIMIZE ESTO
+                print "A = ", A.shape
+                print "b = ", b.shape
+                print "lambda1=", lambda1.shape
                 u = irls_dn(A,b,p,lambda1);

trunk/src/src/sfb.py +4 -12

                 N = 2*lo.size;
                 L = sf.size/2;
-                #print 'N', N
-                #print 'sf', sf
-                #print 'sf[:,0]', sf[:,0].shape
-                #print 'sf[:,1]', sf[:,1].shape
-                #print 'sbf hi', hi.shape
                 # Need to change format for upfirdn funct:
                 lo = lo.T.conj()
                 lo = lo.reshape(lo.size)
-                print 'sfb hi', hi
+                #print 'sfb hi', hi
                 # Need to change format for upfirdn funct:
                 hi = hi.T.conj()
                 hi = hi.reshape(hi.size)
                 #hi = hi.reshape(1, hi.size)
                 lo = upfirdn(lo, sf[:,0], 2, 1);
                 hi = upfirdn(hi, sf[:,1], 2, 1);
                 y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT
                 y = y[0:N];
-                print 'y en sbf\n', y.shape
+                #print 'y en sbf\n', y.shape
                 y = y.reshape(1, y.size)
-                print 'y en sbf\n', y.shape
+                #print 'y en sbf\n', y.shape
                 y = cshift(y, 1-L/2);

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages