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 16) +++ b/trunk/src/src/RICS_paper_compare_noise.py (revision 17) @@ -22,6 +22,7 @@ from modelf import * from scipy.optimize import fsolve from scipy.optimize import root +import pywt ## Calculate Forward Model @@ -79,16 +80,16 @@ # background constant b. # Triple Gaussian -# a = np.array([3, 5, 2]); -# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]); -# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]); -# b = 0; # background +# a = np.array([3, 5, 2]); +# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]); +# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]); +# b = 0; # background # Double Gaussian -# a = np.array([3, 5]); -# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]); -# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]); -# b = 0; # background +# a = np.array([3, 5]); +# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]); +# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]); +# b = 0; # background # Single Gaussian a = np.array( [3] ); @@ -133,7 +134,7 @@ # factr = factr + wind2 * (-(thetar-(mu+sig))); -# fact = fact/(sum(fact)[0]*2*thbound/Nt); % normalize to integral(f)==1 +# 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 factr = factr/I; # normalize to sum(f)==1 @@ -209,7 +210,7 @@ plt.pcolor(abs(Rnz)); plt.colorbar(); - # Fourier Inversion %%%%%%%%%%%%%%%%%%% + # Fourier Inversion ################### f_fourier = np.zeros(shape=(Nr,1), dtype=complex); for i in range(0, thetar.size): @@ -225,7 +226,7 @@ #print f_fourier - # Capon Inversion %%%%%%%%%%%%%%%%%%%%%% + # Capon Inversion ###################### f_capon = np.zeros(shape=(Nr,1)); @@ -242,9 +243,9 @@ f_capon = f_capon.real; # get rid of numerical imaginary noise - # MaxEnt Inversion %%%%%%%%%%%%%%%%%%%%% + # MaxEnt Inversion ##################### - # create the appropriate sensing matrix (split into real and imaginary % parts) + # create the appropriate sensing matrix (split into real and imaginary # parts) M = (r.size-1)*(r.size); Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction @@ -305,21 +306,125 @@ - # CS inversion using irls %%%%%%%%%%%%%%%%%%%%%%%% + # CS inversion using irls ######################## # (Use Nr, thetar, gnz, and Hr from MaxEnt above) - Psi = deb4_basis(Nr); + #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?) + + wavelet1 = pywt.Wavelet('db4') + Phi, Psi, x = wavelet1.wavefun(level=3) - # # add "sum to 1" constraint - # H2 = [Hr; ones(1,Nr)]; - # g2 = [gnz; Nr/Nt]; - # - # s = irls_dn2(H2*Psi,g2,0.5,G); + # add "sum to 1" constraint + H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); + N_temp = np.array([[Nr/Nt]]); + g2 = np.concatenate( (gnz, N_temp), axis=0 ); + + + #s = irls_dn2(H2*Psi,g2,0.5,G); # f_cs = Psi*s; # # # plot # plot(thetar,f_cs,'r.-'); # hold on; # plot(thetat,fact,'k-'); - # hold off; \ No newline at end of file + # hold off; + + + # # # Scaling and shifting + # # # Only necessary for capon solution + + + f_capon = f_capon/np.max(f_capon)*np.max(fact); + + + ### analyze stuff ###################### + # calculate MSE + rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2)); + rmse_capon = np.sqrt(np.mean((f_capon - factr)**2)); + rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2)); + #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2)); + + relrmse_fourier = rmse_fourier / np.linalg.norm(fact); + relrmse_capon = rmse_capon / np.linalg.norm(fact); + relrmse_maxent = rmse_maxent / np.linalg.norm(fact); + #relrmse_cs = rmse_cs / np.norm(fact); + + factr = factr.T.conj() + + # calculate correlation + corr_fourier = np.dot(f_fourier,factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr)); + corr_capon = np.dot(f_capon,factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr)); + corr_maxent = np.dot(f_maxent,factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr)); + #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr)); + + # calculate centered correlation + f0 = factr - np.mean(factr); + f1 = f_fourier - np.mean(f_fourier); + corrc_fourier = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); + f1 = f_capon - np.mean(f_capon); + corrc_capon = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); + f1 = f_maxent - np.mean(f_maxent); + corrc_maxent = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); + #f1 = f_cs - mean(f_cs); + #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1)); + + + + # # # plot stuff ######################### + + #---- Capon---- + plt.figure(4) + plt.subplot(2, 1, 1) + plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon') + plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') + plt.ylabel('Power (arbitrary units)') + plt.legend(loc='upper right') + + # 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) + + + #---- MaxEnt---- + plt.subplot(2, 1, 2) + plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt') + plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') + plt.ylabel('Power (arbitrary units)') + plt.legend(loc='upper right') + + # 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) + + plt.show() + + +# # subplot(3,1,2); +# # plot(180/pi*thetar,f_maxent,'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_maxent, corr_maxent, corrc_maxent)); +# # ylabel({'Power';'(arbitrary units)'}) +# # # title 'Maximum Entropy Method' +# # legend('MaxEnt','Truth'); +# # +# # 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); \ No newline at end of file Index: trunk/src/src/modelf.py =================================================================== diff --git a/trunk/src/src/modelf.py b/trunk/src/src/modelf.py --- a/trunk/src/src/modelf.py (revision 16) +++ b/trunk/src/src/modelf.py (revision 17) @@ -9,7 +9,7 @@ def modelf(lambda1, H, F): # THRESH = 700; - exponent = H.T.conj()*lambda1; + exponent = np.dot(H.T.conj(),lambda1); # exponent(exponent>THRESH) = THRESH; % numerical stabilizer # exponent(exponent<-THRESH) = -THRESH; 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 16) +++ b/trunk/src/src/y_hysell96.py (revision 17) @@ -4,6 +4,7 @@ @author: Yolian Amaro ''' +import numpy as np from modelf import * def y_hysell96(lambda1,g,sigma,F,G,H): @@ -27,5 +28,5 @@ # measurement equation y = g + e - np.dot(H, f); - + return y