jro_soft/compressed_sensing Commit - r23:24

Cleaned code, fixed bugs/errors, updated plotting part....

yamaro -

r23:24

parent child

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

             @author: Yolian Amaro
             '''
-            import pywt
+            #import pywt
             import numpy as np
             def FSfarras():

trunk/src/src/RICS_paper_compare_noise.py +45 -75

             import time
             import matplotlib.pyplot as plt
             from scipy.optimize import root
+            from scipy.stats import nanmean
-            from y_hysell96 import*
+            from y_hysell96 import *
             from deb4_basis import *
             from modelf import *
             from irls_dn2 import *
-            #from scipy.optimize import fsolve
             #-------------------------------------------------------------------------------------------------
             # fact = fact + wind2 * (-(theta-(mu+sig)));
             # factr = factr + wind2 * (-(thetar-(mu+sig)));
             # fact = fact/(sum(fact)[0]*2*thbound/Nt);  # normalize to integral(f)==1
             I = sum(fact)[0];
             fact = fact/I;                              # normalize to sum(f)==1
                     for i in range(0, thetar.size):
                         th = thetar[i];
                         w = np.exp(1j*k*np.dot(r,np.sin(th)));
-                        f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
+                        f_capon[i] = 1/ ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real
                     toc_capon = time.time()
                     elapsed_time_capon = toc_capon - tic_capon;
                     f_capon = f_capon.real;                 # get rid of numerical imaginary noise
                     ystar = myfun(lambda1);
                     Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G));
                     ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
-                    es = np.dot(Hr, f_maxent) - gnz;        # should be same as ep
+                    es = np.dot(Hr, f_maxent) - gnz;
                     chi2 = np.sum((es/sigma)**2);
                     # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
                     Psi = deb4_basis(Nr);
-                    # ------------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 );
                     g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
                     tic_cs = time.time();
-            #         plt.figure(4)
-            #         plt.imshow(Psi)#, interpolation='nearest')
-            #         #plt.xticks([]); plt.yticks([])
-            #         plt.show()
                     # Inversion
                     s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
-                    #print s
+                    toc_cs = time.time()
+                    elapsed_time_cs = toc_cs - tic_cs;
                     # Brightness function
                     f_cs = np.dot(Psi,s);
-                    toc_cs = time.time()
-                    elapsed_time_cs = toc_cs - tic_cs;
                     # Plot
                     plt.figure(4)
                     # Plot stuff
                     #-------------------------------------------------------------------------------------------------
-                    #---- Capon----
+                    #---- Capon----#
                     plt.figure(5)
                     plt.subplot(3, 1, 1)
                     plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon')
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
-                    #---- MaxEnt----
+                    #---- MaxEnt----#
                     plt.subplot(3, 1, 2)
                     plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt')
                     plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
-                    #---- Compressed Sensing----
+                    #---- Compressed Sensing----#
                     plt.subplot(3, 1, 3)
                     plt.plot(180/np.pi*thetar, f_cs, 'r', label='CS')
                     plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
                     # formatting y-axis
                     locs,labels = plt.yticks()
                     plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
-                    # #     PLOT PARA COMPRESSED SENSING
-                    # #
-                    # #     subplot(3,1,3);
-                    # #         plot(180/pi*thetar,f_cs,'r-');
-                    # #         hold on;
-                    # #         plot(180/pi*thetat,fact,'k--');
-                    # #         hold off;
-                    # #         ylim([min(f_cs) 1.1*max(fact)]);
-                    # # #         title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
-                    # # #         title 'Compressed Sensing - Debauchies Wavelets'
-                    # #         xlabel 'Degrees'
-                    # #         ylabel({'Power';'(arbitrary units)'})
-                    # #         legend('Comp. Sens.','Truth');
-                    # #
-                    # # #     set(gcf,'Position',[749   143   528   881]); # CSL
-                    # # #     set(gcf,'Position',[885   -21   528   673]); # macbook
-                    # #     pause(0.01);
                     # # Store Results
                     corr[0, snri, Ni] = corr_fourier;
                     corrc[2,snri,Ni] = corrc_maxent;
                     corrc[3,snri,Ni] = corrc_cs;
+                    print "--------Time performace--------"
                     print 'Capon:\t', elapsed_time_capon, 'sec';
                     print 'Maxent:\t',elapsed_time_maxent, 'sec';
-                    print 'CS:\t',elapsed_time_cs, 'sec';
+                    print 'CS:\t',elapsed_time_cs, 'sec\n';
-                    print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
+                    print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN), '\n';
-                    print corr
+                    #-------------------------------------------------------------------------------------------------
-                    print corr.shape
+                    # Analyze and plot statistics
+                    #-------------------------------------------------------------------------------------------------
-                    ## Analyze and plot statistics
                     metric = corr; # set this to rmse, corr, or corrc
                     # Remove outliers (this part was experimental and wasn't used in the paper)
                     #     end
                     # end
-                    # Avg ignoring NaNs
-            def nanmean(data, **args):
+            # Avg ignoring NaNs
-                return numpy.ma.filled(numpy.ma.masked_array(data,numpy.isnan(data)).mean(**args), fill_value=numpy.nan)
+            ave = np.zeros(shape=(4))
-            # ave = np.zeros(shape=(4))
+            ave[0] = nanmean(metric[0,:,:]);    # Fourier
+            ave[1] = nanmean(metric[1,:,:]);    # Capon
-            # ave[0] = nanmean(metric, axis=0);
+            ave[2] = nanmean(metric[2,:,:]);    # MaxEnt
-            # ave[1] = nanmean(metric, axis=1);
+            ave[3] = nanmean(metric[3,:,:]);    # Compressed Sensing
-            # ave[2] = nanmean(metric, axis=2);
-            # ave[3] = nanmean(metric, axis=3);
+            # Plot based on chosen metric
-            #print ave
             plt.figure(6);
-            f = plt.scatter(SNRdBvec, corr[0], marker='+', color='b', s=60);   # Fourier
+            f = plt.scatter(SNRdBvec, ave[0], marker='+', color='b', s=60);     # Fourier
-            c = plt.scatter(SNRdBvec, corr[1], marker='o', color= 'c', s=60);     # Capon
+            c = plt.scatter(SNRdBvec, ave[1], marker='o', color= 'g', s=60);    # Capon
-            me= plt.scatter(SNRdBvec, corr[2], marker='s', color= 'y', s=60);    # MaxEnt
+            me= plt.scatter(SNRdBvec, ave[2], marker='s', color= 'c', s=60);    # MaxEnt
-            cs= plt.scatter(SNRdBvec, corr[3], marker='*', color='r', s=60);        # Compressed Sensing
+            cs= plt.scatter(SNRdBvec, ave[3], marker='*', color='r', s=60);     # Compressed Sensing
             plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right')
             plt.xlabel('SNR')
             plt.ylabel('Correlation with Truth')
+            print "--------Correlations--------"
+            print "Fourier:", ave[0]
+            print "Capon:\t", ave[1]
+            print "MaxEnt:\t", ave[2]
+            print "CS:\t", ave[3]
             plt.show()

trunk/src/src/deb4_basis.py 0 -1

             @author: Yolian Amaro
             '''
-            import numpy as np
             from FSfarras import *
             from dualfilt1 import *
             from dualtree import *

trunk/src/src/irls_dn.py +5 -29

             @author: Yolian Amaro
             '''
-            #from scipy.sparse import eye
             from scipy import linalg
             import scipy.sparse as sps
             import numpy as np
                 [Q,R] = linalg.qr(A.T.conj(), mode='economic');
-                c = linalg.solve(R.T.conj(),b); # will be used later also
+                c = linalg.solve(R.T.conj(),b);     # will be used later also
-                u = np.dot(Q,c); # minimum 2-norm solution
+                u = np.dot(Q,c);                    # minimum 2-norm solution
                 I = sps.eye(M);
                 # Temporary N x N matrix
                 temp = np.zeros(shape=(N,N))
-                #---------- not needed, defined above--------------
-                # Spacing of floating point numbers
-                #eps = np.spacing(1)
-                #--------------------------------------------------
                 # Loop until damping parameter is small enough
                 while (eps > 1e-7):
                         # Update
                         w = (u**2 + eps)**(1-p/2.0);
-            #             #---- Very inefficient- REMOVE THIS PART------
-            #             k = 0
-            #             # Sparse matrix
-            #             for i in range (0, N):
-            #                 for j in range (0,N):
-            #                     if(i==j):
-            #                         temp[i,j] = w[k]
-            #                         k = k+1
-                        #--------------------------------------------------
                         np.fill_diagonal(temp, w)
-                        #-----------------------------------------------
                         # Compressed Sparse Matrix
                         W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
                         WAT = W*A.T.conj();
-                        #print "WAT", WAT.shape
-                        #print "np.dot(A,WAT)", np.dot(A,WAT).shape
-                        #print "np.dot(lambda1,I)", np.dot(lambda1,I).shape
-                        #print "linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b)", linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b).shape
                         u_new = np.dot(WAT , linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b));

trunk/src/src/irls_dn2.py +8 -29

             '''
             from irls_dn import *
-            from scipy.optimize import fsolve
+            from scipy.optimize import brentq
-            import numpy as np
-            from scipy.optimize import *
-            from dogleg import *
-            from numpy.linalg import norm
-            import matplotlib.pyplot as plt
             def irls_dn2(A,b,p,G):
                 alpha = 2.0; # Line search parameter
                 lambda1 = 1e5; # What's a reasonable but safe initial guess?
                 u = irls_dn(A,b,p,lambda1);
-                #print "u\n", u
                 fid = norm(np.dot(A,u)-b)**2;
                     lambda1 = lambda1 / alpha; # Balance between speed and accuracy
                     u = irls_dn(A,b,p,lambda1);
                     fid = norm(np.dot(A,u)-b)**2;
-                    print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
+                    print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid;
-                    #print u
                 # Refinement using fzero/ brentq
-                lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
+                lambda0 = np.array([lambda1,lambda1*alpha]);    # interval with zero-crossing
                 def myfun(lambda1):
-                    temp1 = np.dot(A, irls_dn(A,b,p,lambda1))
+                    return norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
-                    temp2 = norm(temp1-b)
-                    temp3 = temp2**2-G
-                    #return np.linalg.norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
-                    return temp3
-                print "tolerancia=", 0.01*lambda1
-                #lambda1 = root(myfun, lambda1, method='krylov',  tol=0.01*lambda1);
+                # Find zero-crossing at given interval (lambda1, lambda1*alpha)
-                #lambda1 = lambda1.x
-                print "lambda0[0]", lambda0[0]
-                print "lambda0[1]", lambda0[1]
                 lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1)
-                print "lambda final=", lambda1
                 u = irls_dn(A,b,p,lambda1);
                 return u;

trunk/src/src/sfb.py +1 -2

             @author: Yolian Amaro
             '''
-            from multirate import *
+            from multirate import upfirdn
-            import numpy as np
             from cshift import *
             def sfb(lo, hi, sf):

trunk/src/src/y_hysell96.py 0 -1

             @author: Yolian Amaro
             '''
-            import numpy as np
             from modelf import *
             def y_hysell96(lambda1,g,sigma,F,G,H):

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