Index: trunk/src/src/FSfarras.py =================================================================== diff --git a/trunk/src/src/FSfarras.py b/trunk/src/src/FSfarras.py --- a/trunk/src/src/FSfarras.py (revision 23) +++ b/trunk/src/src/FSfarras.py (revision 24) @@ -4,7 +4,7 @@ @author: Yolian Amaro ''' -import pywt +#import pywt import numpy as np def FSfarras(): Index: trunk/src/src/RICS_paper_compare_noise.py =================================================================== diff --git a/trunk/src/src/RICS_paper_compare_noise.py b/trunk/src/src/RICS_paper_compare_noise.py --- a/trunk/src/src/RICS_paper_compare_noise.py (revision 23) +++ b/trunk/src/src/RICS_paper_compare_noise.py (revision 24) @@ -11,12 +11,12 @@ import time import matplotlib.pyplot as plt from scipy.optimize import root - -from y_hysell96 import* +from scipy.stats import nanmean + +from y_hysell96 import * from deb4_basis import * from modelf import * from irls_dn2 import * -#from scipy.optimize import fsolve #------------------------------------------------------------------------------------------------- @@ -136,8 +136,8 @@ # 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 @@ -232,10 +232,9 @@ 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 @@ -300,7 +299,7 @@ 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); @@ -310,34 +309,23 @@ # (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) @@ -393,7 +381,7 @@ # Plot stuff #------------------------------------------------------------------------------------------------- - #---- Capon---- + #---- Capon----# plt.figure(5) plt.subplot(3, 1, 1) plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon') @@ -407,7 +395,7 @@ 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') @@ -420,7 +408,7 @@ 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') @@ -430,27 +418,7 @@ # 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); - + plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) # # Store Results corr[0, snri, Ni] = corr_fourier; @@ -468,20 +436,18 @@ 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 (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN); - - print corr - - print corr.shape - - - ## Analyze and plot statistics - + print 'CS:\t',elapsed_time_cs, 'sec\n'; + + print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN), '\n'; + + + #------------------------------------------------------------------------------------------------- + # 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) @@ -499,28 +465,32 @@ # end # end - # Avg ignoring NaNs -def nanmean(data, **args): - return numpy.ma.filled(numpy.ma.masked_array(data,numpy.isnan(data)).mean(**args), fill_value=numpy.nan) - -# ave = np.zeros(shape=(4)) -# -# ave[0] = nanmean(metric, axis=0); -# ave[1] = nanmean(metric, axis=1); -# ave[2] = nanmean(metric, axis=2); -# ave[3] = nanmean(metric, axis=3); - -#print ave + +# Avg ignoring NaNs +ave = np.zeros(shape=(4)) + +ave[0] = nanmean(metric[0,:,:]); # Fourier +ave[1] = nanmean(metric[1,:,:]); # Capon +ave[2] = nanmean(metric[2,:,:]); # MaxEnt +ave[3] = nanmean(metric[3,:,:]); # Compressed Sensing + +# Plot based on chosen metric plt.figure(6); -f = plt.scatter(SNRdBvec, corr[0], marker='+', color='b', s=60); # Fourier -c = plt.scatter(SNRdBvec, corr[1], marker='o', color= 'c', s=60); # Capon -me= plt.scatter(SNRdBvec, corr[2], marker='s', color= 'y', s=60); # MaxEnt -cs= plt.scatter(SNRdBvec, corr[3], marker='*', color='r', s=60); # Compressed Sensing - - +f = plt.scatter(SNRdBvec, ave[0], marker='+', color='b', s=60); # Fourier +c = plt.scatter(SNRdBvec, ave[1], marker='o', color= 'g', s=60); # Capon +me= plt.scatter(SNRdBvec, ave[2], marker='s', color= 'c', s=60); # MaxEnt +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() + Index: trunk/src/src/deb4_basis.py =================================================================== diff --git a/trunk/src/src/deb4_basis.py b/trunk/src/src/deb4_basis.py --- a/trunk/src/src/deb4_basis.py (revision 23) +++ b/trunk/src/src/deb4_basis.py (revision 24) @@ -4,7 +4,6 @@ @author: Yolian Amaro ''' -import numpy as np from FSfarras import * from dualfilt1 import * from dualtree import * Index: trunk/src/src/irls_dn.py =================================================================== diff --git a/trunk/src/src/irls_dn.py b/trunk/src/src/irls_dn.py --- a/trunk/src/src/irls_dn.py (revision 23) +++ b/trunk/src/src/irls_dn.py (revision 24) @@ -4,7 +4,6 @@ @author: Yolian Amaro ''' -#from scipy.sparse import eye from scipy import linalg import scipy.sparse as sps import numpy as np @@ -27,17 +26,12 @@ [Q,R] = linalg.qr(A.T.conj(), mode='economic'); - c = linalg.solve(R.T.conj(),b); # will be used later also - u = np.dot(Q,c); # minimum 2-norm solution + c = linalg.solve(R.T.conj(),b); # will be used later also + 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) - #-------------------------------------------------- + temp = np.zeros(shape=(N,N)) # Loop until damping parameter is small enough while (eps > 1e-7): @@ -51,30 +45,12 @@ # 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)); Index: trunk/src/src/irls_dn2.py =================================================================== diff --git a/trunk/src/src/irls_dn2.py b/trunk/src/src/irls_dn2.py --- a/trunk/src/src/irls_dn2.py (revision 23) +++ b/trunk/src/src/irls_dn2.py (revision 24) @@ -5,13 +5,7 @@ ''' from irls_dn import * -from scipy.optimize import fsolve -import numpy as np -from scipy.optimize import * -from dogleg import * -from numpy.linalg import norm - -import matplotlib.pyplot as plt +from scipy.optimize import brentq def irls_dn2(A,b,p,G): @@ -31,7 +25,7 @@ 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; @@ -41,34 +35,19 @@ 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 u + print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid; # 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)) - 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 + return norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G; - #lambda1 = root(myfun, lambda1, method='krylov', tol=0.01*lambda1); - #lambda1 = lambda1.x - - print "lambda0[0]", lambda0[0] - print "lambda0[1]", lambda0[1] - + # Find zero-crossing at given interval (lambda1, lambda1*alpha) lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1) - - print "lambda final=", lambda1 u = irls_dn(A,b,p,lambda1); - return u; Index: trunk/src/src/sfb.py =================================================================== diff --git a/trunk/src/src/sfb.py b/trunk/src/src/sfb.py --- a/trunk/src/src/sfb.py (revision 23) +++ b/trunk/src/src/sfb.py (revision 24) @@ -4,8 +4,7 @@ @author: Yolian Amaro ''' -from multirate import * -import numpy as np +from multirate import upfirdn from cshift import * def sfb(lo, hi, sf): Index: trunk/src/src/y_hysell96.py =================================================================== diff --git a/trunk/src/src/y_hysell96.py b/trunk/src/src/y_hysell96.py --- a/trunk/src/src/y_hysell96.py (revision 23) +++ b/trunk/src/src/y_hysell96.py (revision 24) @@ -4,7 +4,6 @@ @author: Yolian Amaro ''' -import numpy as np from modelf import * def y_hysell96(lambda1,g,sigma,F,G,H):