jro_soft/compressed_sensing Commit - r16:17

Plots for Capon and MaxEntropy inversion methods. Compressed Sensing initial/partial implementation.

yamaro -

r16:17

parent child

trunk/src/src/RICS_paper_compare_noise.py +126 -21

              from modelf import *
              from scipy.optimize import fsolve
              from scipy.optimize import root
+             import pywt
              ## Calculate Forward Model
              # 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] );
              # 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
                      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):
                      #print f_fourier
-                     # Capon Inversion %%%%%%%%%%%%%%%%%%%%%%
+                     # Capon Inversion ######################
                      f_capon = np.zeros(shape=(Nr,1));
                      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
-                     # 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

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

              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;

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

              @author: Yolian Amaro
              '''
+             import numpy as np
              from modelf import *
              def y_hysell96(lambda1,g,sigma,F,G,H):
                  # measurement equation
                  y = g + e - np.dot(H, f);
                  return y

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