##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
Plots for Capon and MaxEntropy inversion methods. Compressed Sensing initial/partial implementation.
yamaro -
r16:17
parent child
Show More
Show file before | Show file after | Hide comments
@@ -22,6 +22,7
22 from modelf import *
22 from modelf import *
23 from scipy.optimize import fsolve
23 from scipy.optimize import fsolve
24 from scipy.optimize import root
24 from scipy.optimize import root
25 import pywt
25
26
26
27
27 ## Calculate Forward Model
28 ## Calculate Forward Model
@@ -79,16 +80,16
79 # background constant b.
80 # background constant b.
80
81
81 # Triple Gaussian
82 # Triple Gaussian
82 # a = np.array([3, 5, 2]);
83 # a = np.array([3, 5, 2]);
83 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
84 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
84 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
85 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
85 # b = 0; # background
86 # b = 0; # background
86
87
87 # Double Gaussian
88 # Double Gaussian
88 # a = np.array([3, 5]);
89 # a = np.array([3, 5]);
89 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
90 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
90 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
91 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
91 # b = 0; # background
92 # b = 0; # background
92
93
93 # Single Gaussian
94 # Single Gaussian
94 a = np.array( [3] );
95 a = np.array( [3] );
@@ -133,7 +134,7
133 # factr = factr + wind2 * (-(thetar-(mu+sig)));
134 # factr = factr + wind2 * (-(thetar-(mu+sig)));
134
135
135
136
136 # fact = fact/(sum(fact)[0]*2*thbound/Nt); % normalize to integral(f)==1
137 # fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1
137 I = sum(fact)[0];
138 I = sum(fact)[0];
138 fact = fact/I; # normalize to sum(f)==1
139 fact = fact/I; # normalize to sum(f)==1
139 factr = factr/I; # normalize to sum(f)==1
140 factr = factr/I; # normalize to sum(f)==1
@@ -209,7 +210,7
209 plt.pcolor(abs(Rnz));
210 plt.pcolor(abs(Rnz));
210 plt.colorbar();
211 plt.colorbar();
211
212
212 # Fourier Inversion %%%%%%%%%%%%%%%%%%%
213 # Fourier Inversion ###################
213 f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
214 f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
214
215
215 for i in range(0, thetar.size):
216 for i in range(0, thetar.size):
@@ -225,7 +226,7
225 #print f_fourier
226 #print f_fourier
226
227
227
228
228 # Capon Inversion %%%%%%%%%%%%%%%%%%%%%%
229 # Capon Inversion ######################
229
230
230 f_capon = np.zeros(shape=(Nr,1));
231 f_capon = np.zeros(shape=(Nr,1));
231
232
@@ -242,9 +243,9
242
243
243 f_capon = f_capon.real; # get rid of numerical imaginary noise
244 f_capon = f_capon.real; # get rid of numerical imaginary noise
244
245
245 # MaxEnt Inversion %%%%%%%%%%%%%%%%%%%%%
246 # MaxEnt Inversion #####################
246
247
247 # create the appropriate sensing matrix (split into real and imaginary % parts)
248 # create the appropriate sensing matrix (split into real and imaginary # parts)
248 M = (r.size-1)*(r.size);
249 M = (r.size-1)*(r.size);
249 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
250 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
250 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
251 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
@@ -305,21 +306,125
305
306
306
307
307
308
308 # CS inversion using irls %%%%%%%%%%%%%%%%%%%%%%%%
309 # CS inversion using irls ########################
309
310
310 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
311 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
311
312
312 Psi = deb4_basis(Nr);
313 #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
314
315 wavelet1 = pywt.Wavelet('db4')
316 Phi, Psi, x = wavelet1.wavefun(level=3)
313
317
314 # # add "sum to 1" constraint
318 # add "sum to 1" constraint
315 # H2 = [Hr; ones(1,Nr)];
319 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
316 # g2 = [gnz; Nr/Nt];
320 N_temp = np.array([[Nr/Nt]]);
317 #
321 g2 = np.concatenate( (gnz, N_temp), axis=0 );
318 # s = irls_dn2(H2*Psi,g2,0.5,G);
322
323
324 #s = irls_dn2(H2*Psi,g2,0.5,G);
319 # f_cs = Psi*s;
325 # f_cs = Psi*s;
320 #
326 #
321 # # plot
327 # # plot
322 # plot(thetar,f_cs,'r.-');
328 # plot(thetar,f_cs,'r.-');
323 # hold on;
329 # hold on;
324 # plot(thetat,fact,'k-');
330 # plot(thetat,fact,'k-');
325 # hold off; No newline at end of file
331 # hold off;
332
333
334 # # # Scaling and shifting
335 # # # Only necessary for capon solution
336
337
338 f_capon = f_capon/np.max(f_capon)*np.max(fact);
339
340
341 ### analyze stuff ######################
342 # calculate MSE
343 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
344 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
345 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
346 #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
347
348 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
349 relrmse_capon = rmse_capon / np.linalg.norm(fact);
350 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
351 #relrmse_cs = rmse_cs / np.norm(fact);
352
353 factr = factr.T.conj()
354
355 # calculate correlation
356 corr_fourier = np.dot(f_fourier,factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
357 corr_capon = np.dot(f_capon,factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
358 corr_maxent = np.dot(f_maxent,factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
359 #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
360
361 # calculate centered correlation
362 f0 = factr - np.mean(factr);
363 f1 = f_fourier - np.mean(f_fourier);
364 corrc_fourier = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
365 f1 = f_capon - np.mean(f_capon);
366 corrc_capon = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
367 f1 = f_maxent - np.mean(f_maxent);
368 corrc_maxent = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
369 #f1 = f_cs - mean(f_cs);
370 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
371
372
373
374 # # # plot stuff #########################
375
376 #---- Capon----
377 plt.figure(4)
378 plt.subplot(2, 1, 1)
379 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
380 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
381 plt.ylabel('Power (arbitrary units)')
382 plt.legend(loc='upper right')
383
384 # formatting y-axis
385 locs,labels = plt.yticks()
386 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
387 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
388
389
390 #---- MaxEnt----
391 plt.subplot(2, 1, 2)
392 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
393 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
394 plt.ylabel('Power (arbitrary units)')
395 plt.legend(loc='upper right')
396
397 # formatting y-axis
398 locs,labels = plt.yticks()
399 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
400 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
401
402 plt.show()
403
404
405 # # subplot(3,1,2);
406 # # plot(180/pi*thetar,f_maxent,'r-');
407 # # hold on;
408 # # plot(180/pi*thetat,fact,'k--');
409 # # hold off;
410 # # ylim([min(f_cs) 1.1*max(fact)]);
411 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_maxent, corr_maxent, corrc_maxent));
412 # # ylabel({'Power';'(arbitrary units)'})
413 # # # title 'Maximum Entropy Method'
414 # # legend('MaxEnt','Truth');
415 # #
416 # # subplot(3,1,3);
417 # # plot(180/pi*thetar,f_cs,'r-');
418 # # hold on;
419 # # plot(180/pi*thetat,fact,'k--');
420 # # hold off;
421 # # ylim([min(f_cs) 1.1*max(fact)]);
422 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
423 # # # title 'Compressed Sensing - Debauchies Wavelets'
424 # # xlabel 'Degrees'
425 # # ylabel({'Power';'(arbitrary units)'})
426 # # legend('Comp. Sens.','Truth');
427 # #
428 # # # set(gcf,'Position',[749 143 528 881]); # CSL
429 # # # set(gcf,'Position',[885 -21 528 673]); # macbook
430 # # pause(0.01); No newline at end of file
Show file before | Show file after | Hide comments
@@ -9,7 +9,7
9 def modelf(lambda1, H, F):
9 def modelf(lambda1, H, F):
10
10
11 # THRESH = 700;
11 # THRESH = 700;
12 exponent = H.T.conj()*lambda1;
12 exponent = np.dot(H.T.conj(),lambda1);
13
13
14 # exponent(exponent>THRESH) = THRESH; % numerical stabilizer
14 # exponent(exponent>THRESH) = THRESH; % numerical stabilizer
15 # exponent(exponent<-THRESH) = -THRESH;
15 # exponent(exponent<-THRESH) = -THRESH;
Show file before | Show file after | Hide comments
@@ -4,6 +4,7
4 @author: Yolian Amaro
4 @author: Yolian Amaro
5 '''
5 '''
6
6
7 import numpy as np
7 from modelf import *
8 from modelf import *
8
9
9 def y_hysell96(lambda1,g,sigma,F,G,H):
10 def y_hysell96(lambda1,g,sigma,F,G,H):
@@ -27,5 +28,5
27
28
28 # measurement equation
29 # measurement equation
29 y = g + e - np.dot(H, f);
30 y = g + e - np.dot(H, f);
30
31
31 return y
32 return y
General Comments 0
You need to be logged in to leave comments. Login now