##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
Cleaned code, fixed bugs/errors, updated plotting part....
yamaro -
r23:24
parent child
Show More
Show file before | Show file after | Hide comments
@@ -4,7 +4,7
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 import pywt
7 #import pywt
8 8 import numpy as np
9 9
10 10 def FSfarras():
Show file before | Show file after | Hide comments
@@ -11,12 +11,12
11 11 import time
12 12 import matplotlib.pyplot as plt
13 13 from scipy.optimize import root
14
15 from y_hysell96 import*
14 from scipy.stats import nanmean
15
16 from y_hysell96 import *
16 17 from deb4_basis import *
17 18 from modelf import *
18 19 from irls_dn2 import *
19 #from scipy.optimize import fsolve
20 20
21 21
22 22 #-------------------------------------------------------------------------------------------------
@@ -136,8 +136,8
136 136 # fact = fact + wind2 * (-(theta-(mu+sig)));
137 137 # factr = factr + wind2 * (-(thetar-(mu+sig)));
138 138
139
140 139 # fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1
140
141 141
142 142 I = sum(fact)[0];
143 143 fact = fact/I; # normalize to sum(f)==1
@@ -232,10 +232,9
232 232 for i in range(0, thetar.size):
233 233 th = thetar[i];
234 234 w = np.exp(1j*k*np.dot(r,np.sin(th)));
235 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
235 f_capon[i] = 1/ ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real
236 236
237 237 toc_capon = time.time()
238
239 238 elapsed_time_capon = toc_capon - tic_capon;
240 239
241 240 f_capon = f_capon.real; # get rid of numerical imaginary noise
@@ -300,7 +299,7
300 299 ystar = myfun(lambda1);
301 300 Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G));
302 301 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
303 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
302 es = np.dot(Hr, f_maxent) - gnz;
304 303 chi2 = np.sum((es/sigma)**2);
305 304
306 305
@@ -310,34 +309,23
310 309 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
311 310
312 311 Psi = deb4_basis(Nr);
313
314 # ------------Remove this-------------------------------------------
315 # wavelet1 = pywt.Wavelet('db4')
316 # Phi, Psi, x = wavelet1.wavefun(level=3)
317 #-------------------------------------------------------------------
318 312
319 313 # add "sum to 1" constraint
320 314 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
321 315 g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
322 316
323 317 tic_cs = time.time();
324
325
326 # plt.figure(4)
327 # plt.imshow(Psi)#, interpolation='nearest')
328 # #plt.xticks([]); plt.yticks([])
329 # plt.show()
330 318
331 319 # Inversion
332 320 s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
333 321
334 #print s
322 toc_cs = time.time()
323 elapsed_time_cs = toc_cs - tic_cs;
335 324
336 325 # Brightness function
337 326 f_cs = np.dot(Psi,s);
338 327
339 toc_cs = time.time()
340 elapsed_time_cs = toc_cs - tic_cs;
328
341 329
342 330 # Plot
343 331 plt.figure(4)
@@ -393,7 +381,7
393 381 # Plot stuff
394 382 #-------------------------------------------------------------------------------------------------
395 383
396 #---- Capon----
384 #---- Capon----#
397 385 plt.figure(5)
398 386 plt.subplot(3, 1, 1)
399 387 plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon')
@@ -407,7 +395,7
407 395 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
408 396
409 397
410 #---- MaxEnt----
398 #---- MaxEnt----#
411 399 plt.subplot(3, 1, 2)
412 400 plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt')
413 401 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
@@ -420,7 +408,7
420 408 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
421 409
422 410
423 #---- Compressed Sensing----
411 #---- Compressed Sensing----#
424 412 plt.subplot(3, 1, 3)
425 413 plt.plot(180/np.pi*thetar, f_cs, 'r', label='CS')
426 414 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
@@ -430,27 +418,7
430 418 # formatting y-axis
431 419 locs,labels = plt.yticks()
432 420 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
433 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
434
435
436 # # PLOT PARA COMPRESSED SENSING
437 # #
438 # # subplot(3,1,3);
439 # # plot(180/pi*thetar,f_cs,'r-');
440 # # hold on;
441 # # plot(180/pi*thetat,fact,'k--');
442 # # hold off;
443 # # ylim([min(f_cs) 1.1*max(fact)]);
444 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
445 # # # title 'Compressed Sensing - Debauchies Wavelets'
446 # # xlabel 'Degrees'
447 # # ylabel({'Power';'(arbitrary units)'})
448 # # legend('Comp. Sens.','Truth');
449 # #
450 # # # set(gcf,'Position',[749 143 528 881]); # CSL
451 # # # set(gcf,'Position',[885 -21 528 673]); # macbook
452 # # pause(0.01);
453
421 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
454 422
455 423 # # Store Results
456 424 corr[0, snri, Ni] = corr_fourier;
@@ -468,20 +436,18
468 436 corrc[2,snri,Ni] = corrc_maxent;
469 437 corrc[3,snri,Ni] = corrc_cs;
470 438
471
439 print "--------Time performace--------"
472 440 print 'Capon:\t', elapsed_time_capon, 'sec';
473 441 print 'Maxent:\t',elapsed_time_maxent, 'sec';
474 print 'CS:\t',elapsed_time_cs, 'sec';
475
476 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
477
478 print corr
479
480 print corr.shape
481
482
483 ## Analyze and plot statistics
484
442 print 'CS:\t',elapsed_time_cs, 'sec\n';
443
444 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN), '\n';
445
446
447 #-------------------------------------------------------------------------------------------------
448 # Analyze and plot statistics
449 #-------------------------------------------------------------------------------------------------
450
485 451 metric = corr; # set this to rmse, corr, or corrc
486 452
487 453 # Remove outliers (this part was experimental and wasn't used in the paper)
@@ -499,28 +465,32
499 465 # end
500 466 # end
501 467
502 # Avg ignoring NaNs
503 def nanmean(data, **args):
504 return numpy.ma.filled(numpy.ma.masked_array(data,numpy.isnan(data)).mean(**args), fill_value=numpy.nan)
505
506 # ave = np.zeros(shape=(4))
507 #
508 # ave[0] = nanmean(metric, axis=0);
509 # ave[1] = nanmean(metric, axis=1);
510 # ave[2] = nanmean(metric, axis=2);
511 # ave[3] = nanmean(metric, axis=3);
512
513 #print ave
468
469 # Avg ignoring NaNs
470 ave = np.zeros(shape=(4))
471
472 ave[0] = nanmean(metric[0,:,:]); # Fourier
473 ave[1] = nanmean(metric[1,:,:]); # Capon
474 ave[2] = nanmean(metric[2,:,:]); # MaxEnt
475 ave[3] = nanmean(metric[3,:,:]); # Compressed Sensing
476
477 # Plot based on chosen metric
514 478 plt.figure(6);
515 f = plt.scatter(SNRdBvec, corr[0], marker='+', color='b', s=60); # Fourier
516 c = plt.scatter(SNRdBvec, corr[1], marker='o', color= 'c', s=60); # Capon
517 me= plt.scatter(SNRdBvec, corr[2], marker='s', color= 'y', s=60); # MaxEnt
518 cs= plt.scatter(SNRdBvec, corr[3], marker='*', color='r', s=60); # Compressed Sensing
519
520
479 f = plt.scatter(SNRdBvec, ave[0], marker='+', color='b', s=60); # Fourier
480 c = plt.scatter(SNRdBvec, ave[1], marker='o', color= 'g', s=60); # Capon
481 me= plt.scatter(SNRdBvec, ave[2], marker='s', color= 'c', s=60); # MaxEnt
482 cs= plt.scatter(SNRdBvec, ave[3], marker='*', color='r', s=60); # Compressed Sensing
483
521 484 plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right')
522 485 plt.xlabel('SNR')
523 486 plt.ylabel('Correlation with Truth')
524 487
488 print "--------Correlations--------"
489 print "Fourier:", ave[0]
490 print "Capon:\t", ave[1]
491 print "MaxEnt:\t", ave[2]
492 print "CS:\t", ave[3]
493
525 494 plt.show()
526 495
496
Show file before | Show file after | Hide comments
@@ -4,7 +4,6
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 import numpy as np
8 7 from FSfarras import *
9 8 from dualfilt1 import *
10 9 from dualtree import *
Show file before | Show file after | Hide comments
@@ -4,7 +4,6
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 #from scipy.sparse import eye
8 7 from scipy import linalg
9 8 import scipy.sparse as sps
10 9 import numpy as np
@@ -27,17 +26,12
27 26 [Q,R] = linalg.qr(A.T.conj(), mode='economic');
28 27
29 28
30 c = linalg.solve(R.T.conj(),b); # will be used later also
31 u = np.dot(Q,c); # minimum 2-norm solution
29 c = linalg.solve(R.T.conj(),b); # will be used later also
30 u = np.dot(Q,c); # minimum 2-norm solution
32 31 I = sps.eye(M);
33 32
34 33 # Temporary N x N matrix
35 temp = np.zeros(shape=(N,N))
36
37 #---------- not needed, defined above--------------
38 # Spacing of floating point numbers
39 #eps = np.spacing(1)
40 #--------------------------------------------------
34 temp = np.zeros(shape=(N,N))
41 35
42 36 # Loop until damping parameter is small enough
43 37 while (eps > 1e-7):
@@ -51,30 +45,12
51 45
52 46 # Update
53 47 w = (u**2 + eps)**(1-p/2.0);
54
55 # #---- Very inefficient- REMOVE THIS PART------
56 # k = 0
57 # # Sparse matrix
58 # for i in range (0, N):
59 # for j in range (0,N):
60 # if(i==j):
61 # temp[i,j] = w[k]
62 # k = k+1
63 #--------------------------------------------------
64
65 48 np.fill_diagonal(temp, w)
66 #-----------------------------------------------
67
49
68 50 # Compressed Sparse Matrix
69 51 W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
70
71
52
72 53 WAT = W*A.T.conj();
73
74 #print "WAT", WAT.shape
75 #print "np.dot(A,WAT)", np.dot(A,WAT).shape
76 #print "np.dot(lambda1,I)", np.dot(lambda1,I).shape
77 #print "linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b)", linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b).shape
78 54
79 55 u_new = np.dot(WAT , linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b));
80 56
Show file before | Show file after | Hide comments
@@ -5,13 +5,7
5 5 '''
6 6
7 7 from irls_dn import *
8 from scipy.optimize import fsolve
9 import numpy as np
10 from scipy.optimize import *
11 from dogleg import *
12 from numpy.linalg import norm
13
14 import matplotlib.pyplot as plt
8 from scipy.optimize import brentq
15 9
16 10 def irls_dn2(A,b,p,G):
17 11
@@ -31,7 +25,7
31 25 alpha = 2.0; # Line search parameter
32 26 lambda1 = 1e5; # What's a reasonable but safe initial guess?
33 27 u = irls_dn(A,b,p,lambda1);
34 #print "u\n", u
28
35 29
36 30 fid = norm(np.dot(A,u)-b)**2;
37 31
@@ -41,34 +35,19
41 35 lambda1 = lambda1 / alpha; # Balance between speed and accuracy
42 36 u = irls_dn(A,b,p,lambda1);
43 37 fid = norm(np.dot(A,u)-b)**2;
44 print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
45 #print u
38 print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid;
46 39
47 40 # Refinement using fzero/ brentq
48 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
49
50
41
42 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
43
51 44 def myfun(lambda1):
52 temp1 = np.dot(A, irls_dn(A,b,p,lambda1))
53 temp2 = norm(temp1-b)
54 temp3 = temp2**2-G
55 #return np.linalg.norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
56 return temp3
57
58 print "tolerancia=", 0.01*lambda1
45 return norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
59 46
60 #lambda1 = root(myfun, lambda1, method='krylov', tol=0.01*lambda1);
61 #lambda1 = lambda1.x
62
63 print "lambda0[0]", lambda0[0]
64 print "lambda0[1]", lambda0[1]
65
47 # Find zero-crossing at given interval (lambda1, lambda1*alpha)
66 48 lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1)
67
68 print "lambda final=", lambda1
69 49
70 50 u = irls_dn(A,b,p,lambda1);
71 51
72
73 52 return u;
74 53
Show file before | Show file after | Hide comments
@@ -4,8 +4,7
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 from multirate import *
8 import numpy as np
7 from multirate import upfirdn
9 8 from cshift import *
10 9
11 10 def sfb(lo, hi, sf):
Show file before | Show file after | Hide comments
@@ -4,7 +4,6
4 4 @author: Yolian Amaro
5 5 '''
6 6
7 import numpy as np
8 7 from modelf import *
9 8
10 9 def y_hysell96(lambda1,g,sigma,F,G,H):
General Comments 0
You need to be logged in to leave comments. Login now