##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
yamaro -
r21:22
parent child
Show More
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (515 lines changed) Show them Hide them
@@ -2,127 +2,133
2 2
3 3 #----------------------------------------------------------
4 4 # Original MATLAB code developed by Brian Harding
5 # Rewritten in python by Yolian Amaro
5 # Rewritten in Python by Yolian Amaro
6 6 # Python version 2.7
7 7 # May 15, 2014
8 8 # Jicamarca Radio Observatory
9 9 #----------------------------------------------------------
10 10
11 import math
12 import numpy as np
11 import time
13 12 import matplotlib.pyplot as plt
14 from scipy import linalg
15 import time
13 from scipy.optimize import root
14
16 15 from y_hysell96 import*
17 16 from deb4_basis import *
18 17 from modelf import *
18 from irls_dn2 import *
19 19 #from scipy.optimize import fsolve
20 from scipy.optimize import root
21 import pywt
22 from irls_dn2 import *
23
20
21
22 #-------------------------------------------------------------------------------------------------
23 # Set parameters
24 #-------------------------------------------------------------------------------------------------
24 25
25 26 ## Calculate Forward Model
26 27 lambda1 = 6.0
27 k = 2*math.pi/lambda1
28
29 ## Calculate Magnetic Declination
30
31 # [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this
32
33 # or calculate it with the above function
28 k = 2*np.pi/lambda1
29
30 ## Magnetic Declination
34 31 dec = -1.24
35
36 # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
32
33 ## Loads Jicamarca antenna positions
34 antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False)
35
36 ## rx and ry -- for plotting purposes
37 37 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
38 38 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
39 39
40 antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240],
41 [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
42
40 ## Plot of antenna positions
43 41 plt.figure(1)
44 42 plt.plot(rx, ry, 'ro')
45 43 plt.draw()
46 44
47 # Jicamarca is nominally at a 45 degree angle
45 ## Jicamarca is nominally at a 45 degree angle
48 46 theta = 45 - dec;
49 47
48 ## Rotation matrix from antenna coord to magnetic coord (East North)
49 theta_rad = np.radians(theta) # trig functions take radians as argument
50 val1 = float( np.cos(theta_rad) )
51 val2 = float( np.sin(theta_rad) )
52 val3 = float( -1*np.sin(theta_rad))
53 val4 = float( np.cos(theta_rad) )
54
50 55 # Rotation matrix from antenna coord to magnetic coord (East North)
51 theta_rad = math.radians(theta) # trig functions take radians as argument
52 val1 = float( math.cos(theta_rad) )
53 val2 = float( math.sin(theta_rad) )
54 val3 = float( -1*math.sin(theta_rad))
55 val4 = float( math.cos(theta_rad) )
56
57 # Rotation matrix from antenna coord to magnetic coord (East North)
58 R = np.array( [[val1, val3], [val2, val4]] );
56 R = np.array( [[val1, val3], [val2, val4]] )
59 57
60 58 # Rotate antenna positions to magnetic coord.
61 AR = np.dot(R.T, antpos);
59 AR = np.dot(R.T, antpos)
62 60
63 61 # Only take the East component
64 62 r = AR[0,:]
65 r.sort() # ROW VECTOR?
63 r.sort()
66 64
67 65 # Truth model (high and low resolution)
68 Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
69 thbound = 9.0/180*math.pi; # the width of the domain in angle space
70 thetat = np.linspace(-thbound, thbound,Nt) # image domain
71 thetat = np.transpose(thetat) # transpose # FUNCIONA??????????????????????????????
72 Nr = (256.0); # number of pixels in reconstructed image: low res
66 Nt = (1024.0)*(16.0) # number of pixels in truth image: high resolution
67 thbound = 9.0/180*np.pi # the width of the domain in angle space
68 thetat = np.linspace(-thbound, thbound,Nt) # image domain
69 thetat = thetat.T # transpose
70 Nr = (256.0) # number of pixels in reconstructed image: low res
73 71 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
74 thetar = np.transpose(thetar) #transpose # FUNCIONA?????????????????????????????
75
76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
77 # background constant b.
72 thetar = thetar.T # transpose
73
74
75 #-------------------------------------------------------------------------------------------------
76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b.
77 #-------------------------------------------------------------------------------------------------
78 78
79 79 # Triple Gaussian
80 80 # a = np.array([3, 5, 2]);
81 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
82 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
83 # b = 0; # background
81 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi, 7.0/180*np.pi]);
82 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi, 0.3/180*np.pi]);
83 # b = 0; # background
84 84
85 85 # Double Gaussian
86 86 # a = np.array([3, 5]);
87 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
88 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
89 # b = 0; # background
87 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi]);
88 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi]);
89 # b = 0; # background
90 90
91 91 # Single Gaussian
92 a = np.array( [3] );
93 mu = np.array( [-3.0/180*math.pi] )
94 sig = np.array( [2.0/180*math.pi] )
95 b = 0;
96
97 fact = np.zeros(shape=(Nt,1));
98 factr = np.zeros(shape=(Nr,1));
99
92 a = np.array( [3] )
93 mu = np.array( [-3.0/180*np.pi] )
94 sig = np.array( [2.0/180*np.pi] )
95 b = 0
96
97 # Empty matrices for factors
98 fact = np.zeros(shape=(Nt,1))
99 factr = np.zeros(shape=(Nr,1))
100
101 # DFT Kernels
100 102 for i in range(0, a.size):
101 103 temp = (-(thetat-mu[i])**2/(sig[i]**2))
102 104 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
103 105 for j in range(0, temp.size):
104 fact[j] = fact[j] + a[i]*math.exp(temp[j]);
106 fact[j] = fact[j] + a[i]*np.exp(temp[j]);
105 107 for m in range(0, tempr.size):
106 factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
108 factr[m] = factr[m] + a[i]*np.exp(tempr[m]);
109
107 110 fact = fact + b;
108 111 factr = factr + b;
109 112
110 # # model for f: Square pulse
113 # #-------------------------------------------------------------------------------------------------
114 # # Model for f: Square pulse
115 # #-------------------------------------------------------------------------------------------------
111 116 # for j in range(0, fact.size):
112 # if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi):
117 # if (theta > -5.0/180*np.pi and theta < 2.0/180*np.pi):
113 118 # fact[j] = 0
114 119 # else:
115 120 # fact[j] = 1
116 121 # for k in range(0, factr.size):
117 # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
118 # fact[k] = 0
122 # if (thetar[k] > -5.0/180*np.pi and thetar[k] < 2/180*np.pi):
123 # factr[k] = 0
119 124 # else:
120 # fact[k] = 1
121 #
122 #
123 # # model for f: triangle pulse
124 # mu = -1.0/180*math.pi;
125 # sig = 5.0/180*math.pi;
125 # factr[k] = 1
126
127 # #-------------------------------------------------------------------------------------------------
128 # # Model for f: Triangle pulse
129 # #-------------------------------------------------------------------------------------------------
130 # mu = -1.0/180*np.pi;
131 # sig = 5.0/180*np.pi;
126 132 # wind1 = theta > mu-sig and theta < mu;
127 133 # wind2 = theta < mu+sig and theta > mu;
128 134 # fact = wind1 * (theta - (mu - sig));
@@ -131,45 +137,43
131 137 # factr = factr + wind2 * (-(thetar-(mu+sig)));
132 138
133 139
134 # fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1
140 # fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1
141
135 142 I = sum(fact)[0];
136 fact = fact/I; # normalize to sum(f)==1
137 factr = factr/I; # normalize to sum(f)==1
138 #plt.figure()
139 #plt.plot(thetat,fact,'r');
140 #plt.plot(thetar,factr,'k.');
141 #xlim([min(thetat) max(thetat)]);
142
143 #x = np.linspace(thetat.min(), thetat.max) ????
144 #for i in range(0, thetat.size):
143 fact = fact/I; # normalize to sum(f)==1
144 factr = factr/I; # normalize to sum(f)==1
145
146 # Plot Gaussian pulse(s)
145 147 plt.figure(2)
146 148 plt.plot(thetat, fact, 'r--')
147 149 plt.plot(thetar, factr, 'ro')
148 150 plt.draw()
149 # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
150
151
152 ##
151
152
153 #-------------------------------------------------------------------------------------------------
153 154 # Control the type and number of inversions with:
154 155 # SNRdBvec: the SNRs that will be used.
155 156 # NN: the number of trials for each SNR
157 #-------------------------------------------------------------------------------------------------
156 158
157 159 #SNRdBvec = np.linspace(5,20,10);
158 SNRdBvec = np.array([15]);
159 NN = 1; # number of trial at each SNR
160
161 # if using vector arguments should be: (4,SNRdBvec.size,NN)
162 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
163 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
164 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
165
166 for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
160 SNRdBvec = np.array([15]); # 15 dB
161 NN = 1; # number of trials at each SNR
162
163 # Statistics simulation (correlation, root mean square)
164 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
165 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
166 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
167
168 # For each SNR and trial
169 for snri in range(0, SNRdBvec.size):
170 SNRdB = SNRdBvec[snri];
171 SNR = 10**(SNRdB/10.0);
172
167 173 for Ni in range(0, NN):
168 SNRdB = SNRdBvec[snri];
169 SNR = 10**(SNRdB/10.0);
170
171 174 # Calculate cross-correlation matrix (Fourier components of image)
172 175 # This is an inefficient way to do this.
176
173 177 R = np.zeros(shape=(r.size, r.size), dtype=object);
174 178
175 179 for i1 in range(0, r.size):
@@ -185,48 +189,42
185 189 # This is a way of adding noise while maintaining the
186 190 # positive-semi-definiteness of the matrix.
187 191
188 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
192 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
189 193
190 194 sigma_noise = (np.linalg.norm(U,'fro')/SNR);
191 195
192 temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
193 temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
194 temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
195 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
196 # temp1 = (-2*np.random.rand(U.shape[0], U.shape[1]) + 2)
197 # temp2 = 1j*(-2*np.random.rand(U.shape[0], U.shape[1]) + 2)
198 # temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
199 # temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
200 #
201 # nz = np.multiply(temp4,temp3)
202
203 nz = np.multiply( sigma_noise * (np.random.randn(U.shape[0]) + 1j*np.random.randn(U.shape[0]))/np.sqrt(2) , (abs(U) > 0).astype(float));
196 204
197 nz = np.multiply(temp4, temp3)
198
199 #---------------------- Eliminar esto:------------------------------------------
200 #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int))
201 #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
202 #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
203 #nz = np.dot(sigma_noise, (np.dot((np.random.rand(8,8) + j*np.random.rand(8,8))/math.sqrt(2.0) , (abs(U) > 0).astype(int))));
204 #--------------------------------------------------------------------------------
205
206 205 Unz = U + nz;
207 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
206 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
208 207 plt.figure(3);
209 208 plt.pcolor(abs(Rnz));
210 209 plt.colorbar();
211
212 # Fourier Inversion ###################
210
211 #-------------------------------------------------------------------------------------------------
212 # Fourier Inversion
213 #-------------------------------------------------------------------------------------------------
213 214 f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
214 215
215 216 for i in range(0, thetar.size):
216 217 th = thetar[i];
217 w = np.exp(1j*k*np.dot(r,np.sin(th)));
218
218 w = np.exp(1j*k*np.dot(r,np.sin(th)));
219 219 temp = np.dot(w.T.conj(),U)
220
221 220 f_fourier[i] = np.dot(temp, w);
222 221
223 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
224
225 #print f_fourier
226
227
228 # Capon Inversion ######################
229
222 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
223
224
225 #-------------------------------------------------------------------------------------------------
226 # Capon Inversion
227 #-------------------------------------------------------------------------------------------------
230 228 f_capon = np.zeros(shape=(Nr,1));
231 229
232 230 tic_capon = time.time();
@@ -236,151 +234,149
236 234 w = np.exp(1j*k*np.dot(r,np.sin(th)));
237 235 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
238 236
239
240 237 toc_capon = time.time()
241 238
242 239 elapsed_time_capon = toc_capon - tic_capon;
243 240
244 f_capon = f_capon.real; # get rid of numerical imaginary noise
245
246 # MaxEnt Inversion #####################
247
248 # create the appropriate sensing matrix (split into real and imaginary # parts)
241 f_capon = f_capon.real; # get rid of numerical imaginary noise
242
243
244 #-------------------------------------------------------------------------------------------------
245 # MaxEnt Inversion
246 #-------------------------------------------------------------------------------------------------
247
248 # Create the appropriate sensing matrix (split into real and imaginary # parts)
249 249 M = (r.size-1)*(r.size);
250 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
251 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
250 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
251 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
252 252
253 # need to re-index our measurements from matrix R into vector g
253 # Need to re-index our measurements from matrix R into vector g
254 254 g = np.zeros(shape=(M,1));
255 gnz = np.zeros(shape=(M,1)); # noisy version of g
255 gnz = np.zeros(shape=(M,1)); # noisy version of g
256 256
257 # triangular indexing to perform this re-indexing
257 # Triangular indexing to perform this re-indexing
258 258 T = np.ones(shape=(r.size,r.size));
259 259 [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
260 260
261 # build H
261 # Build H
262 262 for i1 in range(0, r.size):
263 263 for i2 in range(i1+1, r.size):
264 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward
265 idx1 = 2*idx; # because index starts at 0
264 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward
265 idx1 = 2*idx;
266 266 idx2 = 2*idx+1;
267 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T;
268 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T;
269 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
270 Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
271 g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later
272 g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again
267 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj();
268 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj();
269 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt;
270 Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt;
271 g[idx1] = (R[i1,i2]).real*Nr/Nt;
272 g[idx2] = (R[i1,i2]).imag*Nr/Nt;
273 273 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
274 274 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
275 275
276 # inversion
277 F = Nr/Nt; # normalization
278 sigma = 1; # set to 1 because the difference is accounted for in G
279
280 ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
281 G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2
282
283 tic_maxent = time.time();
284
285 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
276 # Inversion
277 F = Nr/Nt; # normalization
278 sigma = 1; # set to 1 because the difference is accounted for in G
279
280 G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2
281 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
286 282
287 toc_maxent = time.time()
288 elapsed_time_maxent = toc_maxent - tic_maxent;
289 283
290 284 # Whitened solution
291 285 def myfun(lambda1):
292 286 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
293 287
294 tic_maxEnt = time.time();
288
289 tic_maxEnt = time.time(); # start time maxEnt
295 290
296 #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
297 291 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
298 292
299 #print lambda1
300 #print lambda1.x
301
293 toc_maxEnt = time.time()
294 elapsed_time_maxent = toc_maxEnt - tic_maxEnt;
295
296 # Solution
302 297 lambda1 = lambda1.x;
303 298
304 toc_maxEnt = time.time();
305 299 f_maxent = modelf(lambda1, Hr, F);
306 300 ystar = myfun(lambda1);
307 Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G));
301 Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G));
308 302 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
309 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
303 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
310 304 chi2 = np.sum((es/sigma)**2);
311 305
312
313 # --------- CS inversion using Iteratively Reweighted Least Squares (IRLS) -------------
314
306
307 #-------------------------------------------------------------------------------------------------
308 # CS inversion using Iteratively Reweighted Least Squares (IRLS)
309 #-------------------------------------------------------------------------------------------------
315 310 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
316
311
317 312 Psi = deb4_basis(Nr);
318
319 # REMOVE THIS--------------------------------
320 #wavelet1 = pywt.Wavelet('db4')
321 #Phi, Psi, x = wavelet1.wavefun(level=3)
322 # --------------------------------------------
323
313
314 # ------------Remove this-------------------------------------------
315 # wavelet1 = pywt.Wavelet('db4')
316 # Phi, Psi, x = wavelet1.wavefun(level=3)
317 #-------------------------------------------------------------------
318
324 319 # add "sum to 1" constraint
325 320 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
326 N_temp = np.array([[Nr/Nt]]);
327 g2 = np.concatenate( (gnz, N_temp), axis=0 );
328
329 #H2 = H2.T.conj();
330
331 #Psi = Psi.T.conj(); # to align matrices
332
333 ####print 'H2 shape', H2.shape
334 #####print 'Psi shape', Psi.shape
335
336 A = np.dot(H2,Psi);
337
321 g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
322
323 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
331 # Inversion
338 332 s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
339 # f_cs = Psi*s;
340 #
341 # # plot
342 # plot(thetar,f_cs,'r.-');
343 # hold on;
344 # plot(thetat,fact,'k-');
345 # hold off;
346
347
348 # # # Scaling and shifting
349 # # # Only necessary for capon solution
333
334 #print s
335
336 # Brightness function
337 f_cs = np.dot(Psi,s);
338
339 toc_cs = time.time()
340 elapsed_time_cs = toc_cs - tic_cs;
341
342 # Plot
343 plt.figure(4)
344 plt.plot(thetar,f_cs,'r.-');
345 plt.plot(thetat,fact,'k-');
346
347
348 #-------------------------------------------------------------------------------------------------
349 # Scaling and shifting
350 # (Only necessary for Capon solution)
351 #-------------------------------------------------------------------------------------------------
350 352 f_capon = f_capon/np.max(f_capon)*np.max(fact);
351 353
352 354
353 ### analyze stuff ######################
354 # calculate MSE
355 #-------------------------------------------------------------------------------------------------
356 # Analyze stuff
357 #-------------------------------------------------------------------------------------------------
358
359 # Calculate MSE
355 360 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
356 361 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
357 362 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
358 #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
363 rmse_cs = np.sqrt(np.mean((f_cs - factr)**2));
359 364
360 365
361 366 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
362 367 relrmse_capon = rmse_capon / np.linalg.norm(fact);
363 368 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
364 #relrmse_cs = rmse_cs / np.norm(fact);
369 relrmse_cs = rmse_cs / np.linalg.norm(fact);
365 370
366 # To be able to perform dot product (align matrices) done below within the dot calculations
367
368
369 #f_fourier = f_fourier.T.conj()
370 #f_capon = f_capon.T.conj()
371 #f_maxent = f_maxent.T.conj()
372
373 #factr = factr.T.conj()
374
375 # calculate correlation
376
371
372 # Calculate correlation
377 373 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
378 374 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
379 375 corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
380 #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
376 corr_cs = np.dot(f_cs.T.conj(),factr) / (np.linalg.norm(f_cs)*np.linalg.norm(factr));
381 377
382 378
383 # calculate centered correlation
379 # Calculate centered correlation
384 380 f0 = factr - np.mean(factr);
385 381 f1 = f_fourier - np.mean(f_fourier);
386 382
@@ -389,18 +385,19
389 385 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
390 386 f1 = f_maxent - np.mean(f_maxent);
391 387 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
392 #f1 = f_cs - mean(f_cs);
393 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
394
395
396
397 # # # plot stuff #########################
398
388 f1 = f_cs - np.mean(f_cs);
389 corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
390
391
392 #-------------------------------------------------------------------------------------------------
393 # Plot stuff
394 #-------------------------------------------------------------------------------------------------
395
399 396 #---- Capon----
400 plt.figure(4)
401 plt.subplot(2, 1, 1)
402 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
403 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
397 plt.figure(5)
398 plt.subplot(3, 1, 1)
399 plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon')
400 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
404 401 plt.ylabel('Power (arbitrary units)')
405 402 plt.legend(loc='upper right')
406 403
@@ -411,9 +408,9
411 408
412 409
413 410 #---- MaxEnt----
414 plt.subplot(2, 1, 2)
415 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
416 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
411 plt.subplot(3, 1, 2)
412 plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt')
413 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
417 414 plt.ylabel('Power (arbitrary units)')
418 415 plt.legend(loc='upper right')
419 416
@@ -422,7 +419,18
422 419 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
423 420 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
424 421
425 plt.show()
422
423 #---- Compressed Sensing----
424 plt.subplot(3, 1, 3)
425 plt.plot(180/np.pi*thetar, f_cs, 'r', label='CS')
426 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
427 plt.ylabel('Power (arbitrary units)')
428 plt.legend(loc='upper right')
429
430 # formatting y-axis
431 locs,labels = plt.yticks()
432 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)
426 434
427 435
428 436 # # PLOT PARA COMPRESSED SENSING
@@ -448,24 +456,71
448 456 corr[0, snri, Ni] = corr_fourier;
449 457 corr[1, snri, Ni] = corr_capon;
450 458 corr[2, snri, Ni] = corr_maxent;
451 #corr[3, snri, Ni] = corr_cs;
459 corr[3, snri, Ni] = corr_cs;
452 460
453 461 rmse[0,snri,Ni] = relrmse_fourier;
454 462 rmse[1,snri,Ni] = relrmse_capon;
455 463 rmse[2,snri,Ni] = relrmse_maxent;
456 #rmse[3,snri,Ni] = relrmse_cs;
464 rmse[3,snri,Ni] = relrmse_cs;
457 465
458 466 corrc[0,snri,Ni] = corrc_fourier;
459 467 corrc[1,snri,Ni] = corrc_capon;
460 468 corrc[2,snri,Ni] = corrc_maxent;
461 #corrc[3,snri,Ni] = corrc_cs;
469 corrc[3,snri,Ni] = corrc_cs;
462 470
463 471
464 472 print 'Capon:\t', elapsed_time_capon, 'sec';
465 473 print 'Maxent:\t',elapsed_time_maxent, 'sec';
466 #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
474 print 'CS:\t',elapsed_time_cs, 'sec';
467 475
468 476 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
469 477
470 478 print corr
471 No newline at end of file
479
480 print corr.shape
481
482
483 ## Analyze and plot statistics
484
485 metric = corr; # set this to rmse, corr, or corrc
486
487 # Remove outliers (this part was experimental and wasn't used in the paper)
488 # nsig = 3;
489 # for i = 1:4
490 # for snri = 1:length(SNRdBvec)
491 # av = mean(met(i,snri,:));
492 # s = std(met(i,snri,:));
493 # idx = abs(met(i,snri,:) - av) > nsig*s;
494 # met(i,snri,idx) = nan;
495 # if sum(idx)>0
496 # fprintf('i=%i, snr=%i, %i/%i pts removed\n',...
497 # i,round(SNRdBvec(snri)),sum(idx),length(idx));
498 # end
499 # end
500 # end
501
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
514 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
521 plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right')
522 plt.xlabel('SNR')
523 plt.ylabel('Correlation with Truth')
524
525 plt.show()
526
Show file before | Show file after | Hide comments
@@ -50,9 +50,11
50 50 lo = upfirdn(x, af[:,0], 1, 2);
51 51
52 52
53 # VERIFY THIS!!!!!!!!!!!!
54 for i in range(0, L):
55 lo[i] = lo[N/2+i] + lo[i];
53 # # VERIFY THIS!!!!!!!!!!!!
54 # for i in range(0, L):
55 # lo[i] = lo[N/2+i] + lo[i];
56
57 lo[0:L-1] = lo[N/2+np.arange(0,L-1)] + lo[0:L-1]
56 58
57 59 lo = lo[0:N/2];
58 60
@@ -60,8 +62,10
60 62 # highpass filter
61 63 hi = upfirdn(x, af[:,1], 1, 2);
62 64
63 for j in range(0, L):
64 hi[j] = hi[N/2+j] + hi[j];
65 # for j in range(0, L):
66 # hi[j] = hi[N/2+j] + hi[j];
67
68 hi[0:L-1] = hi[N/2+np.arange(0,L-1)] + hi[0:L-1]
65 69
66 70 hi = hi[0:N/2];
67 71
Show file before | Show file after | Hide comments
@@ -9,7 +9,8
9 9 from dualfilt1 import *
10 10 from dualtree import *
11 11 from idualtree import *
12
12
13 # Debauchie 4 Wavelet
13 14 def deb4_basis(N):
14 15
15 16 Psi = np.zeros(shape=(N,2*N+1));
@@ -17,14 +18,15
17 18 J = 4;
18 19 [Faf, Fsf] = FSfarras();
19 20 [af, sf] = dualfilt1();
20
21
21 22 # compute transform of zero vector
22 23 x = np.zeros(shape=(1,N));
23 w = dualtree(x, J, Faf, af);
24 w = dualtree(x, J, Faf, af);
25
24 26
25 27 # Uses both real and imaginary wavelets
26 for i in range (0, J):
27 for j in range (0, 1):
28 for i in range (0, J+1):
29 for j in range (0, 2):
28 30 for k in range (0, (w[i][j]).size):
29 31 w[i][j][0,k] = 1;
30 32 y = idualtree(w, J, Fsf, sf);
Show file before | Show file after | Hide comments
@@ -45,7 +45,7
45 45 #----------------------------------------#
46 46
47 47 # normalization
48 x = x/np.sqrt(2);
48 x = x/np.sqrt(2.0);
49 49
50 50
51 51 w = np.zeros(shape=(J+1), dtype=object)
Show file before | Show file after | Hide comments
@@ -41,6 +41,6
41 41 y2 = sfb(y2, w[0][1], Fsf[0,1]);
42 42
43 43 # normalization
44 y = (y1 + y2)/np.sqrt(2);
44 y = (y1 + y2)/np.sqrt(2.0);
45 45
46 46 return y
Show file before | Show file after | Hide comments
@@ -8,6 +8,7
8 8 from scipy import linalg
9 9 import scipy.sparse as sps
10 10 import numpy as np
11 from numpy.linalg import norm
11 12
12 13 def irls_dn(A,b,p,lambda1):
13 14
@@ -30,6 +31,9
30 31 u = np.dot(Q,c); # minimum 2-norm solution
31 32 I = sps.eye(M);
32 33
34 # Temporary N x N matrix
35 temp = np.zeros(shape=(N,N))
36
33 37 #---------- not needed, defined above--------------
34 38 # Spacing of floating point numbers
35 39 #eps = np.spacing(1)
@@ -39,36 +43,43
39 43 while (eps > 1e-7):
40 44 epschange = 0;
41 45 # Loop until it's time to change eps
42 while (~epschange):
46 while (not(epschange)):
43 47 # main loop
44 48 # u_n = W*A'*(A*W*A'+ lambda*I)^-1 * b
45 49 # where W = diag(1/w)
46 50 # where w = (u.^2 + eps).^(p/2-1)
47 51
48 52 # Update
49 w = (u**2 + eps)**(1-p/2);
53 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 #--------------------------------------------------
50 64
51 # Empty temporary N x N matrix
52 temp = np.zeros(shape=(N,N))
53
54 k = 0
55 # Sparse matrix
56 for i in range (0, N):
57 for j in range (0,N):
58 if(i==j):
59 temp[i,j] = w[k]
60 k = k+1
65 np.fill_diagonal(temp, w)
66 #-----------------------------------------------
61 67
62 68 # Compressed Sparse Matrix
63 69 W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
64 70
65 71
66 72 WAT = W*A.T.conj();
67
68 u_new = np.dot(WAT , linalg.solve(np.dot(A,WAT) + np.dot(lambda1,I), b));
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
79 u_new = np.dot(WAT , linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b));
69 80
70 81 # See if this subproblem is converging
71 delu = np.linalg.norm(u_new-u)/np.linalg.norm(u);
82 delu = norm(u_new-u)/norm(u);
72 83 epschange = delu < (np.sqrt(eps)/100.0);
73 84
74 85 # Make update
@@ -76,8 +87,6
76 87
77 88
78 89 eps = eps/10.0; # decrease eps
79 # Print info
80 #print 'eps =',eps;
81 90
82 91 return u
83 92
Show file before | Show file after | Hide comments
@@ -7,7 +7,11
7 7 from irls_dn import *
8 8 from scipy.optimize import fsolve
9 9 import numpy as np
10 from scipy.optimize import root
10 from scipy.optimize import *
11 from dogleg import *
12 from numpy.linalg import norm
13
14 import matplotlib.pyplot as plt
11 15
12 16 def irls_dn2(A,b,p,G):
13 17
@@ -27,48 +31,42
27 31 alpha = 2.0; # Line search parameter
28 32 lambda1 = 1e5; # What's a reasonable but safe initial guess?
29 33 u = irls_dn(A,b,p,lambda1);
30 fid = np.linalg.norm(np.dot(A,u)-b)**2;
31
34 #print "u\n", u
35
36 fid = norm(np.dot(A,u)-b)**2;
37
32 38 print '----------------------------------\n';
33
39
34 40 while (fid >= G):
35 41 lambda1 = lambda1 / alpha; # Balance between speed and accuracy
36 42 u = irls_dn(A,b,p,lambda1);
37 fid = np.linalg.norm(np.dot(A,u)-b)**2;
43 fid = norm(np.dot(A,u)-b)**2;
38 44 print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
39
40 # Refinement using fzero
45 #print u
46
47 # Refinement using fzero/ brentq
41 48 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
42
49
50
43 51 def myfun(lambda1):
44 print "A = ", A.shape
45 print "b = ", b.shape
46 lambda1
47 return np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G;
48
49 #f = lambda lambda1: np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G; NOOOOOO
50
51
52 # opts = optimset('fzero');
53 # # opts.Display = 'iter';
54 # opts.Display = 'none';
55 # opts.TolX = 0.01*lambda1;
56
57 #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
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 58 print "tolerancia=", 0.01*lambda1
59 59
60 #lambda1 = root(myfun,lambda0, method='krylov', tol=0.01*lambda1);
60 #lambda1 = root(myfun, lambda1, method='krylov', tol=0.01*lambda1);
61 #lambda1 = lambda1.x
61 62
62
63 print "lamda1=", lambda1
64 print "lambda0=", lambda0
65
66 lambda1 = fsolve(myfun,lambda0); # FALTA OPTIMIZE ESTO
67
68 print "A = ", A.shape
69 print "b = ", b.shape
70 print "lambda1=", lambda1.shape
71
63 print "lambda0[0]", lambda0[0]
64 print "lambda0[1]", lambda0[1]
65
66 lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1)
67
68 print "lambda final=", lambda1
69
72 70 u = irls_dn(A,b,p,lambda1);
73 71
74 72
Show file before | Show file after | Hide comments
@@ -45,8 +45,12
45 45
46 46 lo = upfirdn(lo, sf[:,0], 2, 1);
47 47 hi = upfirdn(hi, sf[:,1], 2, 1);
48
49 lo = lo[0:lo.size-1]
50 hi = hi[0:hi.size-1]
51
48 52 y = lo + hi;
49 y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT
53 y[0:L-2] = y[0:L-2] + y[N+ np.arange(0,L-2)]; #CHECK IF ARANGE IS CORRECT
50 54 y = y[0:N];
51 55
52 56 #print 'y en sbf\n', y.shape
General Comments 0
You need to be logged in to leave comments. Login now