##// 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 # 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 time
12 import numpy as np
13 import matplotlib.pyplot as plt
12 import matplotlib.pyplot as plt
14 from scipy import linalg
13 from scipy.optimize import root
15 import time
14
16 from y_hysell96 import*
15 from y_hysell96 import*
17 from deb4_basis import *
16 from deb4_basis import *
18 from modelf import *
17 from modelf import *
18 from irls_dn2 import *
19 #from scipy.optimize import fsolve
19 #from scipy.optimize import fsolve
20 from scipy.optimize import root
20
21 import pywt
21
22 from irls_dn2 import *
22 #-------------------------------------------------------------------------------------------------
23
23 # Set parameters
24 #-------------------------------------------------------------------------------------------------
24
25
25 ## Calculate Forward Model
26 ## Calculate Forward Model
26 lambda1 = 6.0
27 lambda1 = 6.0
27 k = 2*math.pi/lambda1
28 k = 2*np.pi/lambda1
28
29
29 ## Calculate Magnetic Declination
30 ## 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
34 dec = -1.24
31 dec = -1.24
35
32
36 # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
33 ## Loads Jicamarca antenna positions
34 antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False)
35
36 ## rx and ry -- for plotting purposes
37 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
37 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
38 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
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],
40 ## Plot of antenna positions
41 [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
42
43 plt.figure(1)
41 plt.figure(1)
44 plt.plot(rx, ry, 'ro')
42 plt.plot(rx, ry, 'ro')
45 plt.draw()
43 plt.draw()
46
44
47 # Jicamarca is nominally at a 45 degree angle
45 ## Jicamarca is nominally at a 45 degree angle
48 theta = 45 - dec;
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 # Rotation matrix from antenna coord to magnetic coord (East North)
55 # Rotation matrix from antenna coord to magnetic coord (East North)
51 theta_rad = math.radians(theta) # trig functions take radians as argument
56 R = np.array( [[val1, val3], [val2, val4]] )
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]] );
59
57
60 # Rotate antenna positions to magnetic coord.
58 # Rotate antenna positions to magnetic coord.
61 AR = np.dot(R.T, antpos);
59 AR = np.dot(R.T, antpos)
62
60
63 # Only take the East component
61 # Only take the East component
64 r = AR[0,:]
62 r = AR[0,:]
65 r.sort() # ROW VECTOR?
63 r.sort()
66
64
67 # Truth model (high and low resolution)
65 # Truth model (high and low resolution)
68 Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
66 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
67 thbound = 9.0/180*np.pi # the width of the domain in angle space
70 thetat = np.linspace(-thbound, thbound,Nt) # image domain
68 thetat = np.linspace(-thbound, thbound,Nt) # image domain
71 thetat = np.transpose(thetat) # transpose # FUNCIONA??????????????????????????????
69 thetat = thetat.T # transpose
72 Nr = (256.0); # number of pixels in reconstructed image: low res
70 Nr = (256.0) # number of pixels in reconstructed image: low res
73 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
71 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
74 thetar = np.transpose(thetar) #transpose # FUNCIONA?????????????????????????????
72 thetar = thetar.T # transpose
75
73
76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
74
77 # background constant b.
75 #-------------------------------------------------------------------------------------------------
76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b.
77 #-------------------------------------------------------------------------------------------------
78
78
79 # Triple Gaussian
79 # Triple Gaussian
80 # a = np.array([3, 5, 2]);
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]);
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*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
82 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi, 0.3/180*np.pi]);
83 # b = 0; # background
83 # b = 0; # background
84
84
85 # Double Gaussian
85 # Double Gaussian
86 # a = np.array([3, 5]);
86 # a = np.array([3, 5]);
87 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
87 # mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi]);
88 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
88 # sig = np.array([2.0/180*np.pi, 1.5/180*np.pi]);
89 # b = 0; # background
89 # b = 0; # background
90
90
91 # Single Gaussian
91 # Single Gaussian
92 a = np.array( [3] );
92 a = np.array( [3] )
93 mu = np.array( [-3.0/180*math.pi] )
93 mu = np.array( [-3.0/180*np.pi] )
94 sig = np.array( [2.0/180*math.pi] )
94 sig = np.array( [2.0/180*np.pi] )
95 b = 0;
95 b = 0
96
96
97 fact = np.zeros(shape=(Nt,1));
97 # Empty matrices for factors
98 factr = np.zeros(shape=(Nr,1));
98 fact = np.zeros(shape=(Nt,1))
99
99 factr = np.zeros(shape=(Nr,1))
100
101 # DFT Kernels
100 for i in range(0, a.size):
102 for i in range(0, a.size):
101 temp = (-(thetat-mu[i])**2/(sig[i]**2))
103 temp = (-(thetat-mu[i])**2/(sig[i]**2))
102 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
104 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
103 for j in range(0, temp.size):
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 for m in range(0, tempr.size):
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 fact = fact + b;
110 fact = fact + b;
108 factr = factr + b;
111 factr = factr + b;
109
112
110 # # model for f: Square pulse
113 # #-------------------------------------------------------------------------------------------------
114 # # Model for f: Square pulse
115 # #-------------------------------------------------------------------------------------------------
111 # for j in range(0, fact.size):
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 # fact[j] = 0
118 # fact[j] = 0
114 # else:
119 # else:
115 # fact[j] = 1
120 # fact[j] = 1
116 # for k in range(0, factr.size):
121 # for k in range(0, factr.size):
117 # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
122 # if (thetar[k] > -5.0/180*np.pi and thetar[k] < 2/180*np.pi):
118 # fact[k] = 0
123 # factr[k] = 0
119 # else:
124 # else:
120 # fact[k] = 1
125 # factr[k] = 1
121 #
126
122 #
127 # #-------------------------------------------------------------------------------------------------
123 # # model for f: triangle pulse
128 # # Model for f: Triangle pulse
124 # mu = -1.0/180*math.pi;
129 # #-------------------------------------------------------------------------------------------------
125 # sig = 5.0/180*math.pi;
130 # mu = -1.0/180*np.pi;
131 # sig = 5.0/180*np.pi;
126 # wind1 = theta > mu-sig and theta < mu;
132 # wind1 = theta > mu-sig and theta < mu;
127 # wind2 = theta < mu+sig and theta > mu;
133 # wind2 = theta < mu+sig and theta > mu;
128 # fact = wind1 * (theta - (mu - sig));
134 # fact = wind1 * (theta - (mu - sig));
@@ -131,45 +137,43
131 # factr = factr + wind2 * (-(thetar-(mu+sig)));
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 I = sum(fact)[0];
142 I = sum(fact)[0];
136 fact = fact/I; # normalize to sum(f)==1
143 fact = fact/I; # normalize to sum(f)==1
137 factr = factr/I; # normalize to sum(f)==1
144 factr = factr/I; # normalize to sum(f)==1
138 #plt.figure()
145
139 #plt.plot(thetat,fact,'r');
146 # Plot Gaussian pulse(s)
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):
145 plt.figure(2)
147 plt.figure(2)
146 plt.plot(thetat, fact, 'r--')
148 plt.plot(thetat, fact, 'r--')
147 plt.plot(thetar, factr, 'ro')
149 plt.plot(thetar, factr, 'ro')
148 plt.draw()
150 plt.draw()
149 # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
151
150
152
151
153 #-------------------------------------------------------------------------------------------------
152 ##
153 # Control the type and number of inversions with:
154 # Control the type and number of inversions with:
154 # SNRdBvec: the SNRs that will be used.
155 # SNRdBvec: the SNRs that will be used.
155 # NN: the number of trials for each SNR
156 # NN: the number of trials for each SNR
157 #-------------------------------------------------------------------------------------------------
156
158
157 #SNRdBvec = np.linspace(5,20,10);
159 #SNRdBvec = np.linspace(5,20,10);
158 SNRdBvec = np.array([15]);
160 SNRdBvec = np.array([15]); # 15 dB
159 NN = 1; # number of trial at each SNR
161 NN = 1; # number of trials at each SNR
160
162
161 # if using vector arguments should be: (4,SNRdBvec.size,NN)
163 # Statistics simulation (correlation, root mean square)
162 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
164 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
163 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
165 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
164 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
166 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
165
167
166 for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
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 for Ni in range(0, NN):
173 for Ni in range(0, NN):
168 SNRdB = SNRdBvec[snri];
169 SNR = 10**(SNRdB/10.0);
170
171 # Calculate cross-correlation matrix (Fourier components of image)
174 # Calculate cross-correlation matrix (Fourier components of image)
172 # This is an inefficient way to do this.
175 # This is an inefficient way to do this.
176
173 R = np.zeros(shape=(r.size, r.size), dtype=object);
177 R = np.zeros(shape=(r.size, r.size), dtype=object);
174
178
175 for i1 in range(0, r.size):
179 for i1 in range(0, r.size):
@@ -185,48 +189,42
185 # This is a way of adding noise while maintaining the
189 # This is a way of adding noise while maintaining the
186 # positive-semi-definiteness of the matrix.
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 sigma_noise = (np.linalg.norm(U,'fro')/SNR);
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)
196 # temp1 = (-2*np.random.rand(U.shape[0], U.shape[1]) + 2)
193 temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
197 # temp2 = 1j*(-2*np.random.rand(U.shape[0], U.shape[1]) + 2)
194 temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
198 # temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
195 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
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 Unz = U + nz;
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 plt.figure(3);
207 plt.figure(3);
209 plt.pcolor(abs(Rnz));
208 plt.pcolor(abs(Rnz));
210 plt.colorbar();
209 plt.colorbar();
211
210
212 # Fourier Inversion ###################
211 #-------------------------------------------------------------------------------------------------
212 # Fourier Inversion
213 #-------------------------------------------------------------------------------------------------
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 temp = np.dot(w.T.conj(),U)
219 temp = np.dot(w.T.conj(),U)
220
221 f_fourier[i] = np.dot(temp, w);
220 f_fourier[i] = np.dot(temp, w);
222
221
223 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
222 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
224
223
225 #print f_fourier
224
226
225 #-------------------------------------------------------------------------------------------------
227
226 # Capon Inversion
228 # Capon Inversion ######################
227 #-------------------------------------------------------------------------------------------------
229
230 f_capon = np.zeros(shape=(Nr,1));
228 f_capon = np.zeros(shape=(Nr,1));
231
229
232 tic_capon = time.time();
230 tic_capon = time.time();
@@ -236,151 +234,149
236 w = np.exp(1j*k*np.dot(r,np.sin(th)));
234 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)
235 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
238
236
239
240 toc_capon = time.time()
237 toc_capon = time.time()
241
238
242 elapsed_time_capon = toc_capon - tic_capon;
239 elapsed_time_capon = toc_capon - tic_capon;
243
240
244 f_capon = f_capon.real; # get rid of numerical imaginary noise
241 f_capon = f_capon.real; # get rid of numerical imaginary noise
245
242
246 # MaxEnt Inversion #####################
243
247
244 #-------------------------------------------------------------------------------------------------
248 # create the appropriate sensing matrix (split into real and imaginary # parts)
245 # MaxEnt Inversion
246 #-------------------------------------------------------------------------------------------------
247
248 # Create the appropriate sensing matrix (split into real and imaginary # parts)
249 M = (r.size-1)*(r.size);
249 M = (r.size-1)*(r.size);
250 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
250 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
251 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
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 g = np.zeros(shape=(M,1));
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 T = np.ones(shape=(r.size,r.size));
258 T = np.ones(shape=(r.size,r.size));
259 [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
260
260
261 # build H
261 # Build H
262 for i1 in range(0, r.size):
262 for i1 in range(0, r.size):
263 for i2 in range(i1+1, r.size):
263 for i2 in range(i1+1, r.size):
264 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
265 idx1 = 2*idx; # because index starts at 0
265 idx1 = 2*idx;
266 idx2 = 2*idx+1;
266 idx2 = 2*idx+1;
267 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.conj();
268 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.conj();
269 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.conj()*Nr/Nt;
270 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.conj()*Nr/Nt;
271 g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later
271 g[idx1] = (R[i1,i2]).real*Nr/Nt;
272 g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again
272 g[idx2] = (R[i1,i2]).imag*Nr/Nt;
273 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
273 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
274 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
274 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
275
275
276 # inversion
276 # Inversion
277 F = Nr/Nt; # normalization
277 F = Nr/Nt; # normalization
278 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
279
279
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 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
282
283 tic_maxent = time.time();
284
285 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 # Whitened solution
284 # Whitened solution
291 def myfun(lambda1):
285 def myfun(lambda1):
292 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
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 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
291 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
298
292
299 #print lambda1
293 toc_maxEnt = time.time()
300 #print lambda1.x
294 elapsed_time_maxent = toc_maxEnt - tic_maxEnt;
301
295
296 # Solution
302 lambda1 = lambda1.x;
297 lambda1 = lambda1.x;
303
298
304 toc_maxEnt = time.time();
305 f_maxent = modelf(lambda1, Hr, F);
299 f_maxent = modelf(lambda1, Hr, F);
306 ystar = myfun(lambda1);
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 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
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 chi2 = np.sum((es/sigma)**2);
304 chi2 = np.sum((es/sigma)**2);
311
305
312
306
313 # --------- CS inversion using Iteratively Reweighted Least Squares (IRLS) -------------
307 #-------------------------------------------------------------------------------------------------
314
308 # CS inversion using Iteratively Reweighted Least Squares (IRLS)
309 #-------------------------------------------------------------------------------------------------
315 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
310 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
316
311
317 Psi = deb4_basis(Nr);
312 Psi = deb4_basis(Nr);
318
313
319 # REMOVE THIS--------------------------------
314 # ------------Remove this-------------------------------------------
320 #wavelet1 = pywt.Wavelet('db4')
315 # wavelet1 = pywt.Wavelet('db4')
321 #Phi, Psi, x = wavelet1.wavefun(level=3)
316 # Phi, Psi, x = wavelet1.wavefun(level=3)
322 # --------------------------------------------
317 #-------------------------------------------------------------------
323
318
324 # add "sum to 1" constraint
319 # add "sum to 1" constraint
325 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
320 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
326 N_temp = np.array([[Nr/Nt]]);
321 g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
327 g2 = np.concatenate( (gnz, N_temp), axis=0 );
322
328
323 tic_cs = time.time();
329 #H2 = H2.T.conj();
324
330
325
331 #Psi = Psi.T.conj(); # to align matrices
326 # plt.figure(4)
332
327 # plt.imshow(Psi)#, interpolation='nearest')
333 ####print 'H2 shape', H2.shape
328 # #plt.xticks([]); plt.yticks([])
334 #####print 'Psi shape', Psi.shape
329 # plt.show()
335
330
336 A = np.dot(H2,Psi);
331 # Inversion
337
338 s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
332 s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
339 # f_cs = Psi*s;
333
340 #
334 #print s
341 # # plot
335
342 # plot(thetar,f_cs,'r.-');
336 # Brightness function
343 # hold on;
337 f_cs = np.dot(Psi,s);
344 # plot(thetat,fact,'k-');
338
345 # hold off;
339 toc_cs = time.time()
346
340 elapsed_time_cs = toc_cs - tic_cs;
347
341
348 # # # Scaling and shifting
342 # Plot
349 # # # Only necessary for capon solution
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 f_capon = f_capon/np.max(f_capon)*np.max(fact);
352 f_capon = f_capon/np.max(f_capon)*np.max(fact);
351
353
352
354
353 ### analyze stuff ######################
355 #-------------------------------------------------------------------------------------------------
354 # calculate MSE
356 # Analyze stuff
357 #-------------------------------------------------------------------------------------------------
358
359 # Calculate MSE
355 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
360 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
356 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
361 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
357 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
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 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
366 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
362 relrmse_capon = rmse_capon / np.linalg.norm(fact);
367 relrmse_capon = rmse_capon / np.linalg.norm(fact);
363 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
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
371
367
372 # Calculate correlation
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
377 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
373 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
378 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
374 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
379 corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
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 f0 = factr - np.mean(factr);
380 f0 = factr - np.mean(factr);
385 f1 = f_fourier - np.mean(f_fourier);
381 f1 = f_fourier - np.mean(f_fourier);
386
382
@@ -389,18 +385,19
389 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
385 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
390 f1 = f_maxent - np.mean(f_maxent);
386 f1 = f_maxent - np.mean(f_maxent);
391 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
387 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
392 #f1 = f_cs - mean(f_cs);
388 f1 = f_cs - np.mean(f_cs);
393 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
389 corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
394
390
395
391
396
392 #-------------------------------------------------------------------------------------------------
397 # # # plot stuff #########################
393 # Plot stuff
398
394 #-------------------------------------------------------------------------------------------------
395
399 #---- Capon----
396 #---- Capon----
400 plt.figure(4)
397 plt.figure(5)
401 plt.subplot(2, 1, 1)
398 plt.subplot(3, 1, 1)
402 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
399 plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon')
403 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
400 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
404 plt.ylabel('Power (arbitrary units)')
401 plt.ylabel('Power (arbitrary units)')
405 plt.legend(loc='upper right')
402 plt.legend(loc='upper right')
406
403
@@ -411,9 +408,9
411
408
412
409
413 #---- MaxEnt----
410 #---- MaxEnt----
414 plt.subplot(2, 1, 2)
411 plt.subplot(3, 1, 2)
415 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
412 plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt')
416 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
413 plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth')
417 plt.ylabel('Power (arbitrary units)')
414 plt.ylabel('Power (arbitrary units)')
418 plt.legend(loc='upper right')
415 plt.legend(loc='upper right')
419
416
@@ -422,7 +419,18
422 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
419 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
423 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
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 # # PLOT PARA COMPRESSED SENSING
436 # # PLOT PARA COMPRESSED SENSING
@@ -448,24 +456,71
448 corr[0, snri, Ni] = corr_fourier;
456 corr[0, snri, Ni] = corr_fourier;
449 corr[1, snri, Ni] = corr_capon;
457 corr[1, snri, Ni] = corr_capon;
450 corr[2, snri, Ni] = corr_maxent;
458 corr[2, snri, Ni] = corr_maxent;
451 #corr[3, snri, Ni] = corr_cs;
459 corr[3, snri, Ni] = corr_cs;
452
460
453 rmse[0,snri,Ni] = relrmse_fourier;
461 rmse[0,snri,Ni] = relrmse_fourier;
454 rmse[1,snri,Ni] = relrmse_capon;
462 rmse[1,snri,Ni] = relrmse_capon;
455 rmse[2,snri,Ni] = relrmse_maxent;
463 rmse[2,snri,Ni] = relrmse_maxent;
456 #rmse[3,snri,Ni] = relrmse_cs;
464 rmse[3,snri,Ni] = relrmse_cs;
457
465
458 corrc[0,snri,Ni] = corrc_fourier;
466 corrc[0,snri,Ni] = corrc_fourier;
459 corrc[1,snri,Ni] = corrc_capon;
467 corrc[1,snri,Ni] = corrc_capon;
460 corrc[2,snri,Ni] = corrc_maxent;
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 print 'Capon:\t', elapsed_time_capon, 'sec';
472 print 'Capon:\t', elapsed_time_capon, 'sec';
465 print 'Maxent:\t',elapsed_time_maxent, 'sec';
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 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
476 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
469
477
470 print corr
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 lo = upfirdn(x, af[:,0], 1, 2);
50 lo = upfirdn(x, af[:,0], 1, 2);
51
51
52
52
53 # VERIFY THIS!!!!!!!!!!!!
53 # # VERIFY THIS!!!!!!!!!!!!
54 for i in range(0, L):
54 # for i in range(0, L):
55 lo[i] = lo[N/2+i] + lo[i];
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 lo = lo[0:N/2];
59 lo = lo[0:N/2];
58
60
@@ -60,8 +62,10
60 # highpass filter
62 # highpass filter
61 hi = upfirdn(x, af[:,1], 1, 2);
63 hi = upfirdn(x, af[:,1], 1, 2);
62
64
63 for j in range(0, L):
65 # for j in range(0, L):
64 hi[j] = hi[N/2+j] + hi[j];
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 hi = hi[0:N/2];
70 hi = hi[0:N/2];
67
71
Show file before | Show file after | Hide comments
@@ -9,7 +9,8
9 from dualfilt1 import *
9 from dualfilt1 import *
10 from dualtree import *
10 from dualtree import *
11 from idualtree import *
11 from idualtree import *
12
12
13 # Debauchie 4 Wavelet
13 def deb4_basis(N):
14 def deb4_basis(N):
14
15
15 Psi = np.zeros(shape=(N,2*N+1));
16 Psi = np.zeros(shape=(N,2*N+1));
@@ -17,14 +18,15
17 J = 4;
18 J = 4;
18 [Faf, Fsf] = FSfarras();
19 [Faf, Fsf] = FSfarras();
19 [af, sf] = dualfilt1();
20 [af, sf] = dualfilt1();
20
21
21 # compute transform of zero vector
22 # compute transform of zero vector
22 x = np.zeros(shape=(1,N));
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 # Uses both real and imaginary wavelets
27 # Uses both real and imaginary wavelets
26 for i in range (0, J):
28 for i in range (0, J+1):
27 for j in range (0, 1):
29 for j in range (0, 2):
28 for k in range (0, (w[i][j]).size):
30 for k in range (0, (w[i][j]).size):
29 w[i][j][0,k] = 1;
31 w[i][j][0,k] = 1;
30 y = idualtree(w, J, Fsf, sf);
32 y = idualtree(w, J, Fsf, sf);
Show file before | Show file after | Hide comments
@@ -45,7 +45,7
45 #----------------------------------------#
45 #----------------------------------------#
46
46
47 # normalization
47 # normalization
48 x = x/np.sqrt(2);
48 x = x/np.sqrt(2.0);
49
49
50
50
51 w = np.zeros(shape=(J+1), dtype=object)
51 w = np.zeros(shape=(J+1), dtype=object)
Show file before | Show file after | Hide comments
@@ -41,6 +41,6
41 y2 = sfb(y2, w[0][1], Fsf[0,1]);
41 y2 = sfb(y2, w[0][1], Fsf[0,1]);
42
42
43 # normalization
43 # normalization
44 y = (y1 + y2)/np.sqrt(2);
44 y = (y1 + y2)/np.sqrt(2.0);
45
45
46 return y
46 return y
Show file before | Show file after | Hide comments
@@ -8,6 +8,7
8 from scipy import linalg
8 from scipy import linalg
9 import scipy.sparse as sps
9 import scipy.sparse as sps
10 import numpy as np
10 import numpy as np
11 from numpy.linalg import norm
11
12
12 def irls_dn(A,b,p,lambda1):
13 def irls_dn(A,b,p,lambda1):
13
14
@@ -30,6 +31,9
30 u = np.dot(Q,c); # minimum 2-norm solution
31 u = np.dot(Q,c); # minimum 2-norm solution
31 I = sps.eye(M);
32 I = sps.eye(M);
32
33
34 # Temporary N x N matrix
35 temp = np.zeros(shape=(N,N))
36
33 #---------- not needed, defined above--------------
37 #---------- not needed, defined above--------------
34 # Spacing of floating point numbers
38 # Spacing of floating point numbers
35 #eps = np.spacing(1)
39 #eps = np.spacing(1)
@@ -39,36 +43,43
39 while (eps > 1e-7):
43 while (eps > 1e-7):
40 epschange = 0;
44 epschange = 0;
41 # Loop until it's time to change eps
45 # Loop until it's time to change eps
42 while (~epschange):
46 while (not(epschange)):
43 # main loop
47 # main loop
44 # u_n = W*A'*(A*W*A'+ lambda*I)^-1 * b
48 # u_n = W*A'*(A*W*A'+ lambda*I)^-1 * b
45 # where W = diag(1/w)
49 # where W = diag(1/w)
46 # where w = (u.^2 + eps).^(p/2-1)
50 # where w = (u.^2 + eps).^(p/2-1)
47
51
48 # Update
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
65 np.fill_diagonal(temp, w)
52 temp = np.zeros(shape=(N,N))
66 #-----------------------------------------------
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
61
67
62 # Compressed Sparse Matrix
68 # Compressed Sparse Matrix
63 W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
69 W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
64
70
65
71
66 WAT = W*A.T.conj();
72 WAT = W*A.T.conj();
67
73
68 u_new = np.dot(WAT , linalg.solve(np.dot(A,WAT) + np.dot(lambda1,I), b));
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 # See if this subproblem is converging
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 epschange = delu < (np.sqrt(eps)/100.0);
83 epschange = delu < (np.sqrt(eps)/100.0);
73
84
74 # Make update
85 # Make update
@@ -76,8 +87,6
76
87
77
88
78 eps = eps/10.0; # decrease eps
89 eps = eps/10.0; # decrease eps
79 # Print info
80 #print 'eps =',eps;
81
90
82 return u
91 return u
83
92
Show file before | Show file after | Hide comments
@@ -7,7 +7,11
7 from irls_dn import *
7 from irls_dn import *
8 from scipy.optimize import fsolve
8 from scipy.optimize import fsolve
9 import numpy as np
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 def irls_dn2(A,b,p,G):
16 def irls_dn2(A,b,p,G):
13
17
@@ -27,48 +31,42
27 alpha = 2.0; # Line search parameter
31 alpha = 2.0; # Line search parameter
28 lambda1 = 1e5; # What's a reasonable but safe initial guess?
32 lambda1 = 1e5; # What's a reasonable but safe initial guess?
29 u = irls_dn(A,b,p,lambda1);
33 u = irls_dn(A,b,p,lambda1);
30 fid = np.linalg.norm(np.dot(A,u)-b)**2;
34 #print "u\n", u
31
35
36 fid = norm(np.dot(A,u)-b)**2;
37
32 print '----------------------------------\n';
38 print '----------------------------------\n';
33
39
34 while (fid >= G):
40 while (fid >= G):
35 lambda1 = lambda1 / alpha; # Balance between speed and accuracy
41 lambda1 = lambda1 / alpha; # Balance between speed and accuracy
36 u = irls_dn(A,b,p,lambda1);
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 print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
44 print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
39
45 #print u
40 # Refinement using fzero
46
47 # Refinement using fzero/ brentq
41 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
48 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
42
49
50
43 def myfun(lambda1):
51 def myfun(lambda1):
44 print "A = ", A.shape
52 temp1 = np.dot(A, irls_dn(A,b,p,lambda1))
45 print "b = ", b.shape
53 temp2 = norm(temp1-b)
46 lambda1
54 temp3 = temp2**2-G
47 return np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G;
55 #return np.linalg.norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
48
56 return temp3
49 #f = lambda lambda1: np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G; NOOOOOO
57
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);
58 print "tolerancia=", 0.01*lambda1
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 "lambda0[0]", lambda0[0]
63 print "lamda1=", lambda1
64 print "lambda0[1]", lambda0[1]
64 print "lambda0=", lambda0
65
65
66 lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1)
66 lambda1 = fsolve(myfun,lambda0); # FALTA OPTIMIZE ESTO
67
67
68 print "lambda final=", lambda1
68 print "A = ", A.shape
69
69 print "b = ", b.shape
70 print "lambda1=", lambda1.shape
71
72 u = irls_dn(A,b,p,lambda1);
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 lo = upfirdn(lo, sf[:,0], 2, 1);
46 lo = upfirdn(lo, sf[:,0], 2, 1);
47 hi = upfirdn(hi, sf[:,1], 2, 1);
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 y = lo + hi;
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 y = y[0:N];
54 y = y[0:N];
51
55
52 #print 'y en sbf\n', y.shape
56 #print 'y en sbf\n', y.shape
General Comments 0
You need to be logged in to leave comments. Login now