##// 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
@@ -1,325 +1,430
1 #!/usr/bin/env python
1 #!/usr/bin/env python
2
2
3 #----------------------------------------------------------
3 #----------------------------------------------------------
4 # Original MATLAB code developed by Brian Harding
4 # Original MATLAB code developed by Brian Harding
5 # Rewritten in python by Yolian Amaro
5 # Rewritten in python by Yolian Amaro
6 # Python version 2.7
6 # Python version 2.7
7 # May 15, 2014
7 # May 15, 2014
8 # Jicamarca Radio Observatory
8 # Jicamarca Radio Observatory
9 #----------------------------------------------------------
9 #----------------------------------------------------------
10
10
11 import math
11 import math
12 #import cmath
12 #import cmath
13 #import scipy
13 #import scipy
14 #import matplotlib
14 #import matplotlib
15 import numpy as np
15 import numpy as np
16 #from numpy import linalg
16 #from numpy import linalg
17 import matplotlib.pyplot as plt
17 import matplotlib.pyplot as plt
18 from scipy import linalg
18 from scipy import linalg
19 import time
19 import time
20 from y_hysell96 import*
20 from y_hysell96 import*
21 from deb4_basis import *
21 from deb4_basis import *
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
28 lambda1 = 6.0
29 lambda1 = 6.0
29 k = 2*math.pi/lambda1
30 k = 2*math.pi/lambda1
30
31
31 ## Calculate Magnetic Declination
32 ## Calculate Magnetic Declination
32
33
33 # [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this
34 # [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this
34
35
35 # or calculate it with the above function
36 # or calculate it with the above function
36 dec = -1.24
37 dec = -1.24
37
38
38 # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
39 # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
39 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
40 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
40 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
41 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
41
42
42 antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240],
43 antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240],
43 [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
44 [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
44
45
45 plt.figure(1)
46 plt.figure(1)
46 plt.plot(rx, ry, 'ro')
47 plt.plot(rx, ry, 'ro')
47 plt.draw()
48 plt.draw()
48
49
49 # Jicamarca is nominally at a 45 degree angle
50 # Jicamarca is nominally at a 45 degree angle
50 theta = 45 - dec;
51 theta = 45 - dec;
51
52
52 # Rotation matrix from antenna coord to magnetic coord (East North)
53 # Rotation matrix from antenna coord to magnetic coord (East North)
53 theta_rad = math.radians(theta) # trig functions take radians as argument
54 theta_rad = math.radians(theta) # trig functions take radians as argument
54 val1 = float( math.cos(theta_rad) )
55 val1 = float( math.cos(theta_rad) )
55 val2 = float( math.sin(theta_rad) )
56 val2 = float( math.sin(theta_rad) )
56 val3 = float( -1*math.sin(theta_rad))
57 val3 = float( -1*math.sin(theta_rad))
57 val4 = float( math.cos(theta_rad) )
58 val4 = float( math.cos(theta_rad) )
58
59
59 # Rotation matrix from antenna coord to magnetic coord (East North)
60 # Rotation matrix from antenna coord to magnetic coord (East North)
60 R = np.array( [[val1, val3], [val2, val4]] );
61 R = np.array( [[val1, val3], [val2, val4]] );
61
62
62 # Rotate antenna positions to magnetic coord.
63 # Rotate antenna positions to magnetic coord.
63 AR = np.dot(R.T, antpos);
64 AR = np.dot(R.T, antpos);
64
65
65 # Only take the East component
66 # Only take the East component
66 r = AR[0,:]
67 r = AR[0,:]
67 r.sort() # ROW VECTOR?
68 r.sort() # ROW VECTOR?
68
69
69 # Truth model (high and low resolution)
70 # Truth model (high and low resolution)
70 Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
71 Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
71 thbound = 9.0/180*math.pi; # the width of the domain in angle space
72 thbound = 9.0/180*math.pi; # the width of the domain in angle space
72 thetat = np.linspace(-thbound, thbound,Nt) # image domain
73 thetat = np.linspace(-thbound, thbound,Nt) # image domain
73 thetat = np.transpose(thetat) # transpose # FUNCIONA??????????????????????????????
74 thetat = np.transpose(thetat) # transpose # FUNCIONA??????????????????????????????
74 Nr = (256.0); # number of pixels in reconstructed image: low res
75 Nr = (256.0); # number of pixels in reconstructed image: low res
75 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
76 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
76 thetar = np.transpose(thetar) #transpose # FUNCIONA?????????????????????????????
77 thetar = np.transpose(thetar) #transpose # FUNCIONA?????????????????????????????
77
78
78 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
79 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
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] );
95 mu = np.array( [-3.0/180*math.pi] )
96 mu = np.array( [-3.0/180*math.pi] )
96 sig = np.array( [2.0/180*math.pi] )
97 sig = np.array( [2.0/180*math.pi] )
97 b = 0;
98 b = 0;
98
99
99 fact = np.zeros(shape=(Nt,1));
100 fact = np.zeros(shape=(Nt,1));
100 factr = np.zeros(shape=(Nr,1));
101 factr = np.zeros(shape=(Nr,1));
101
102
102 for i in range(0, a.size):
103 for i in range(0, a.size):
103 temp = (-(thetat-mu[i])**2/(sig[i]**2))
104 temp = (-(thetat-mu[i])**2/(sig[i]**2))
104 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
105 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
105 for j in range(0, temp.size):
106 for j in range(0, temp.size):
106 fact[j] = fact[j] + a[i]*math.exp(temp[j]);
107 fact[j] = fact[j] + a[i]*math.exp(temp[j]);
107 for m in range(0, tempr.size):
108 for m in range(0, tempr.size):
108 factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
109 factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
109 fact = fact + b;
110 fact = fact + b;
110 factr = factr + b;
111 factr = factr + b;
111
112
112 # # model for f: Square pulse
113 # # model for f: Square pulse
113 # for j in range(0, fact.size):
114 # for j in range(0, fact.size):
114 # if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi):
115 # if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi):
115 # fact[j] = 0
116 # fact[j] = 0
116 # else:
117 # else:
117 # fact[j] = 1
118 # fact[j] = 1
118 # for k in range(0, factr.size):
119 # for k in range(0, factr.size):
119 # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
120 # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
120 # fact[k] = 0
121 # fact[k] = 0
121 # else:
122 # else:
122 # fact[k] = 1
123 # fact[k] = 1
123 #
124 #
124 #
125 #
125 # # model for f: triangle pulse
126 # # model for f: triangle pulse
126 # mu = -1.0/180*math.pi;
127 # mu = -1.0/180*math.pi;
127 # sig = 5.0/180*math.pi;
128 # sig = 5.0/180*math.pi;
128 # wind1 = theta > mu-sig and theta < mu;
129 # wind1 = theta > mu-sig and theta < mu;
129 # wind2 = theta < mu+sig and theta > mu;
130 # wind2 = theta < mu+sig and theta > mu;
130 # fact = wind1 * (theta - (mu - sig));
131 # fact = wind1 * (theta - (mu - sig));
131 # factr = wind1 * (thetar - (mu - sig));
132 # factr = wind1 * (thetar - (mu - sig));
132 # fact = fact + wind2 * (-(theta-(mu+sig)));
133 # fact = fact + wind2 * (-(theta-(mu+sig)));
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
140 #plot(thetat,fact,'r'); hold on;
141 #plot(thetat,fact,'r'); hold on;
141 #plot(thetar,factr,'k.'); hold off;
142 #plot(thetar,factr,'k.'); hold off;
142 #xlim([min(thetat) max(thetat)]);
143 #xlim([min(thetat) max(thetat)]);
143
144
144 #x = np.linspace(thetat.min(), thetat.max) ????
145 #x = np.linspace(thetat.min(), thetat.max) ????
145 #for i in range(0, thetat.size):
146 #for i in range(0, thetat.size):
146 plt.figure(2)
147 plt.figure(2)
147 plt.plot(thetat, fact, 'r--')
148 plt.plot(thetat, fact, 'r--')
148 plt.plot(thetar, factr, 'ro')
149 plt.plot(thetar, factr, 'ro')
149 plt.draw()
150 plt.draw()
150 # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
151 # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
151
152
152
153
153 ##
154 ##
154 # Control the type and number of inversions with:
155 # Control the type and number of inversions with:
155 # SNRdBvec: the SNRs that will be used.
156 # SNRdBvec: the SNRs that will be used.
156 # NN: the number of trials for each SNR
157 # NN: the number of trials for each SNR
157
158
158 #SNRdBvec = np.linspace(5,20,10);
159 #SNRdBvec = np.linspace(5,20,10);
159 SNRdBvec = np.array([15]);
160 SNRdBvec = np.array([15]);
160 NN = 1; # number of trial at each SNR
161 NN = 1; # number of trial at each SNR
161
162
162 # if using vector arguments should be: (4,SNRdBvec.size,NN)
163 # if using vector arguments should be: (4,SNRdBvec.size,NN)
163 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
164 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
164 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
165 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
165 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
166 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
166
167
167 for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
168 for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
168 for Ni in range(0, NN):
169 for Ni in range(0, NN):
169 SNRdB = SNRdBvec[snri];
170 SNRdB = SNRdBvec[snri];
170 SNR = 10**(SNRdB/10.0);
171 SNR = 10**(SNRdB/10.0);
171
172
172 # Calculate cross-correlation matrix (Fourier components of image)
173 # Calculate cross-correlation matrix (Fourier components of image)
173 # This is an inefficient way to do this.
174 # This is an inefficient way to do this.
174 R = np.zeros(shape=(r.size, r.size), dtype=object);
175 R = np.zeros(shape=(r.size, r.size), dtype=object);
175
176
176 for i1 in range(0, r.size):
177 for i1 in range(0, r.size):
177 for i2 in range(0,r.size):
178 for i2 in range(0,r.size):
178 R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat))))
179 R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat))))
179 R[i1,i2] = sum(R[i1,i2])
180 R[i1,i2] = sum(R[i1,i2])
180
181
181 # Add uncertainty
182 # Add uncertainty
182 # This is an ad-hoc way of adding "noise". It models some combination of
183 # This is an ad-hoc way of adding "noise". It models some combination of
183 # receiver noise and finite integration times. We could use a more
184 # receiver noise and finite integration times. We could use a more
184 # advanced model (like in Yu et al 2000) in the future.
185 # advanced model (like in Yu et al 2000) in the future.
185
186
186 # This is a way of adding noise while maintaining the
187 # This is a way of adding noise while maintaining the
187 # positive-semi-definiteness of the matrix.
188 # positive-semi-definiteness of the matrix.
188
189
189 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
190 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
190
191
191 sigma_noise = (np.linalg.norm(U,'fro')/SNR);
192 sigma_noise = (np.linalg.norm(U,'fro')/SNR);
192
193
193 temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
194 temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
194 temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
195 temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
195 temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
196 temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
196 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
197 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
197
198
198 nz = np.multiply(temp4, temp3)
199 nz = np.multiply(temp4, temp3)
199
200
200 #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
201 #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
201
202
202
203
203 #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
204 #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
204 #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))));
205 #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))));
205
206
206 Unz = U + nz;
207 Unz = U + nz;
207 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
208 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
208 plt.figure(3);
209 plt.figure(3);
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):
216 th = thetar[i];
217 th = thetar[i];
217 w = np.exp(1j*k*np.dot(r,np.sin(th)));
218 w = np.exp(1j*k*np.dot(r,np.sin(th)));
218
219
219 temp = np.dot(w.T.conj(),U)
220 temp = np.dot(w.T.conj(),U)
220
221
221 f_fourier[i] = np.dot(temp, w);
222 f_fourier[i] = np.dot(temp, w);
222
223
223 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
224 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
224
225
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
232 #tic_capon = time.time();
233 #tic_capon = time.time();
233
234
234 for i in range(0, thetar.size):
235 for i in range(0, thetar.size):
235 th = thetar[i];
236 th = thetar[i];
236 w = np.exp(1j*k*np.dot(r,np.sin(th)));
237 w = np.exp(1j*k*np.dot(r,np.sin(th)));
237 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
238 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
238
239
239 #toc_capon = time.time()
240 #toc_capon = time.time()
240
241
241 # elapsed_time_capon = toc_capon - tic_capon;
242 # elapsed_time_capon = toc_capon - tic_capon;
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
251
252
252 # 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
253 g = np.zeros(shape=(M,1));
254 g = np.zeros(shape=(M,1));
254 gnz = np.zeros(shape=(M,1)); # noisy version of g
255 gnz = np.zeros(shape=(M,1)); # noisy version of g
255
256
256 # triangular indexing to perform this re-indexing
257 # triangular indexing to perform this re-indexing
257 T = np.ones(shape=(r.size,r.size));
258 T = np.ones(shape=(r.size,r.size));
258 [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
259 [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
259
260
260 # build H
261 # build H
261 for i1 in range(0, r.size):
262 for i1 in range(0, r.size):
262 for i2 in range(i1+1, r.size):
263 for i2 in range(i1+1, r.size):
263 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward
264 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward
264 idx1 = 2*idx; # because index starts at 0
265 idx1 = 2*idx; # because index starts at 0
265 idx2 = 2*idx+1;
266 idx2 = 2*idx+1;
266 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T;
267 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T;
267 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T;
268 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T;
268 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
269 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
269 Ht[idx2,:] = np.sin(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;
270 g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later
271 g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later
271 g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again
272 g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again
272 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
273 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
273 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
274 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
274
275
275 # inversion
276 # inversion
276 F = Nr/Nt; # normalization
277 F = Nr/Nt; # normalization
277 sigma = 1; # set to 1 because the difference is accounted for in G
278 sigma = 1; # set to 1 because the difference is accounted for in G
278
279
279 ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
280 ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below
280 G = np.linalg.norm(g-gnz)**2; # pretend we know in advance the actual value of chi^2
281 G = np.linalg.norm(g-gnz)**2; # pretend we know in advance the actual value of chi^2
281
282
282 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
283 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
283
284
284 # Whitened solution
285 # Whitened solution
285 def myfun(lambda1):
286 def myfun(lambda1):
286 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
287 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
287
288
288 tic_maxEnt = time.time();
289 tic_maxEnt = time.time();
289
290
290 #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
291 #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
291 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
292 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
292
293
293 #print lambda1
294 #print lambda1
294 #print lambda1.x
295 #print lambda1.x
295
296
296 lambda1 = lambda1.x;
297 lambda1 = lambda1.x;
297
298
298 toc_maxEnt = time.time();
299 toc_maxEnt = time.time();
299 f_maxent = modelf(lambda1, Hr, F);
300 f_maxent = modelf(lambda1, Hr, F);
300 ystar = myfun(lambda1);
301 ystar = myfun(lambda1);
301 Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G));
302 Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G));
302 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
303 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
303 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
304 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
304 chi2 = np.sum((es/sigma)**2);
305 chi2 = np.sum((es/sigma)**2);
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
@@ -1,20 +1,20
1 '''
1 '''
2 Created on May 22, 2014
2 Created on May 22, 2014
3
3
4 @author: Yolian Amaro
4 @author: Yolian Amaro
5 '''
5 '''
6
6
7 import numpy as np
7 import numpy as np
8
8
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;
16 E = np.max(exponent);
16 E = np.max(exponent);
17 fnum = np.exp(exponent-E);
17 fnum = np.exp(exponent-E);
18 f = F*fnum/np.sum(fnum);
18 f = F*fnum/np.sum(fnum);
19
19
20 return f No newline at end of file
20 return f
Show file before | Show file after | Hide comments
@@ -1,31 +1,32
1 '''
1 '''
2 Created on May 22, 2014
2 Created on May 22, 2014
3
3
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):
10 # Y_HYSELL96 Implements set of nonlinear equations to solve Hysell96 MaxEnt
11 # Y_HYSELL96 Implements set of nonlinear equations to solve Hysell96 MaxEnt
11 # y(lambda) = 0
12 # y(lambda) = 0
12 # decision variables: lambda
13 # decision variables: lambda
13 # g: data
14 # g: data
14 # sigma: uncertainties (length of g)
15 # sigma: uncertainties (length of g)
15 # F: sum(f)
16 # F: sum(f)
16 # G: desired value for chi^2
17 # G: desired value for chi^2
17 # H: linear operator mapping image (f) to data (g)
18 # H: linear operator mapping image (f) to data (g)
18 # This function is a helper function that returns 0 when a value of lambda
19 # This function is a helper function that returns 0 when a value of lambda
19 # is chosen that satisfies the equations.
20 # is chosen that satisfies the equations.
20
21
21 # model for f
22 # model for f
22 f = modelf(lambda1, H,F);
23 f = modelf(lambda1, H,F);
23
24
24 # solve for Lambda and e
25 # solve for Lambda and e
25 Lambda = np.sqrt(np.sum(np.multiply(lambda1**2,sigma**2))/(4*G)); # positive root (right?)
26 Lambda = np.sqrt(np.sum(np.multiply(lambda1**2,sigma**2))/(4*G)); # positive root (right?)
26 e = np.multiply(-lambda1,sigma**2) / (2*Lambda);
27 e = np.multiply(-lambda1,sigma**2) / (2*Lambda);
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