##// END OF EJS Templates
jro_soft » compressed_sensing
RSS
Partial implementation of wavelet transform. Implemented functions: dualfilt1, FSfarras, irls_dn (partial), irls_dn2 (partial), deb4_basis (partial).
yamaro -
r17:18
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,65
1 '''
2 Created on May 26, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import pywt
8 import numpy as np
9
10 def FSfarras():
11 #function [af, sf] = FSfarras
12
13 # Farras filters organized for the dual-tree
14 # complex DWT.
15 #
16 # USAGE:
17 # [af, sf] = FSfarras
18 # OUTPUT:
19 # af{i}, i = 1,2 - analysis filters for tree i
20 # sf{i}, i = 1,2 - synthesis filters for tree i
21 # See farras, dualtree, dualfilt1.
22 #
23 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
24 # http://taco.poly.edu/WaveletSoftware/
25 #
26 # Translated to Python by Yolian Amaro
27
28
29
30 a1 = np.array( [
31 [ 0, 0],
32 [-0.08838834764832, -0.01122679215254],
33 [ 0.08838834764832, 0.01122679215254],
34 [ 0.69587998903400, 0.08838834764832],
35 [ 0.69587998903400, 0.08838834764832],
36 [ 0.08838834764832, -0.69587998903400],
37 [-0.08838834764832, 0.69587998903400],
38 [ 0.01122679215254, -0.08838834764832],
39 [ 0.01122679215254, -0.08838834764832],
40 [0, 0]
41 ] );
42
43 a2 = np.array([
44 [ 0.01122679215254, 0],
45 [ 0.01122679215254, 0],
46 [-0.08838834764832, -0.08838834764832],
47 [ 0.08838834764832, -0.08838834764832],
48 [ 0.69587998903400, 0.69587998903400],
49 [ 0.69587998903400, -0.69587998903400],
50 [ 0.08838834764832, 0.08838834764832],
51 [-0.08838834764832, 0.08838834764832],
52 [ 0, 0.01122679215254],
53 [ 0, -0.01122679215254]
54 ]);
55
56 #print a2.shape
57
58 af = np.array([ [a1,a2] ], dtype=object)
59
60 s1 = a1[::-1]
61 s2 = a2[::-1]
62
63 sf = np.array([ [s1,s2] ], dtype=object)
64
65 return af, sf
Show file before | Show file after | Hide comments
@@ -0,0 +1,38
1 '''
2 Created on May 26, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import numpy as np
8 import FSfarras
9 import dualfilt1
10
11 def deb4_basis(N):
12
13 Psi = np.zeros(shape=(N,2*N+1));
14 idx = 1;
15
16 J = 4;
17 [Faf, Fsf] = FSfarras;
18 [af, sf] = dualfilt1;
19 # #
20 # # # compute transform of zero vector
21 # # x = zeros(1,N);
22 # # w = dualtree(x, J, Faf, af);
23 # #
24 # # # Uses both real and imaginary wavelets
25 # # for i in range (1, J+1):
26 # # for j in range (1, 2):
27 # # for k in range (1, (w[i][j]).size):
28 # # w[i][j](k) = 1;
29 # # y = idualtree(w, J, Fsf, sf);
30 # # w[i][j](k) = 0;
31 # # # store it
32 # # Psi(:,idx) = y.T.conj();
33 # # idx = idx + 1;
34 # #
35 # # # Add uniform vector (seems to be useful if there's a background
36 # # Psi(:,2*N+1) = 1/np.sqrt(N);
37 #
38 # return Psi No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,65
1 '''
2 Created on May 29, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 import numpy as np
8
9 def dualfilt1():
10
11 # Kingsbury Q-filters for the dual-tree complex DWT
12 #
13 # USAGE:
14 # [af, sf] = dualfilt1
15 # OUTPUT:
16 # af{i}, i = 1,2 - analysis filters for tree i
17 # sf{i}, i = 1,2 - synthesis filters for tree i
18 # note: af{2} is the reverse of af{1}
19 # REFERENCE:
20 # N. G. Kingsbury, "A dual-tree complex wavelet
21 # transform with improved orthogonality and symmetry
22 # properties", Proceedings of the IEEE Int. Conf. on
23 # Image Proc. (ICIP), 2000
24 # See dualtree
25 #
26 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
27 # http://taco.poly.edu/WaveletSoftware/
28
29 # These cofficients are rounded to 8 decimal places.
30
31 a1 = np.array([
32 [ 0.03516384000000, 0],
33 [ 0, 0],
34 [-0.08832942000000, -0.11430184000000],
35 [ 0.23389032000000, 0],
36 [ 0.76027237000000, 0.58751830000000],
37 [ 0.58751830000000, -0.76027237000000],
38 [ 0, 0.23389032000000],
39 [-0.11430184000000, 0.08832942000000],
40 [ 0, 0],
41 [ 0, -0.03516384000000]
42 ]);
43
44 a2 = np.array([
45 [ 0, -0.03516384000000],
46 [ 0, 0],
47 [-0.11430184000000, 0.08832942000000],
48 [ 0, 0.23389032000000],
49 [ 0.58751830000000, -0.76027237000000],
50 [ 0.76027237000000, 0.58751830000000],
51 [ 0.23389032000000, 0],
52 [ -0.08832942000000, -0.11430184000000],
53 [ 0, 0],
54 [ 0.03516384000000, 0]
55 ]);
56
57 af = np.array([ [a1,a2] ], dtype=object)
58
59 s1 = a1[::-1]
60 s2 = a2[::-1]
61
62 sf = np.array([ [s1,s2] ], dtype=object)
63
64
65 return af, sf
Show file before | Show file after | Hide comments
@@ -0,0 +1,80
1 '''
2 Created on May 27, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 #from scipy.sparse import eye
8 from scipy import linalg
9 import scipy.sparse as sps
10 import numpy as np
11
12 def irls_dn(A,b,p,lambda1):
13
14
15 # Minimize lambda*||u||_p + ||A*u-b||_2, 0 < p <= 1
16 # using Iterative Reweighted Least Squares
17 # (see http://math.lanl.gov/Research/Publications/Docs/chartrand-2008-iteratively.pdf
18 # and http://web.eecs.umich.edu/~aey/sparse/sparse11.pdf)
19
20 # Note to self: I found that "warm-starting" didn't really help too much.
21
22 [M,N] = A.shape;
23 # Initialize and precompute:
24 eps = 1e-2; # damping parameter
25 [Q,R] = linalg.qr(A.T.conj(),0);
26 print A.shape
27 print R.shape
28 print b.shape
29 c = linalg.solve(R.T.conj(),b); # will be used later also
30 u = Q*c; # minimum 2-norm solution
31 I = sps.eye(M);
32
33 # Spacing of floating point numbers
34 eps = np.spacing(1)
35
36 # Loop until damping parameter is small enough
37 while (eps > 1e-7):
38 epschange = 0;
39 # Loop until it's time to change eps
40 while (~epschange):
41 # main loop
42 # u_n = W*A'*(A*W*A'+ lambda*I)^-1 * b
43 # where W = diag(1/w)
44 # where w = (u.^2 + eps).^(p/2-1)
45
46 # Update
47 w = (u**2 + eps)**(1-p/2);
48
49 # Empty temporary N x N matrix
50 temp = np.zeros(shape=(N,N))
51
52 # Sparse matrix
53 for i in range (1, N):
54 for j in range (1,N):
55 if(i==j):
56 temp[i,j] = w
57
58 # Compressed Sparse Matrix
59 W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
60
61
62 WAT = W*A.T.conj();
63 u_new = WAT * ( linalg.solve (A*WAT + lambda1*I), b);
64
65 # See if this subproblem is converging
66 delu = np.linalg.norm(u_new-u)/np.linalg.norm(u);
67 epschange = delu < (np.sqrt(eps)/100);
68
69 # Make update
70 u = u_new;
71
72
73 eps = eps/10; # decrease eps
74 # Print info
75 print 'eps =',eps;
76
77 return u
78
79
80
Show file before | Show file after | Hide comments
@@ -0,0 +1,50
1 '''
2 Created on May 27, 2014
3
4 @author: Yolian Amaro
5 '''
6
7 from irls_dn import *
8
9 def irls_dn2(A,b,p,G):
10
11 # Minimize ||u||_p subject to ||A*u-b||_2^2 <= G (0 < p <= 1)
12
13 # What this function actually does is finds the lambda1 so that the solution
14 # to the following problem satisfies ||A*u-b||_2^2 <= G:
15 # Minimize lambda1*||u||_p + ||A*u-b||_2
16
17 # Start with a large lambda1, and do a line search until fidelity <= G.
18 # (Inversions with large lambda1 are really fast anyway).
19
20 # Then spin up fzero to localize the root even better
21
22 # Line Search
23
24 alpha = 2; # Line search parameter
25 lambda1 = 1e5; # What's a reasonable but safe initial guess?
26 u = irls_dn(A,b,p,lambda1);
27 # fid = np.norm(A*u-b)^2;
28 #
29 # print '----------------------------------\n';
30 #
31 # while (fid >= G)
32 # lambda1 = lambda1 / alpha; # Balance between speed and accuracy
33 # u = irls_dn(A,b,p,lambda1);
34 # fid = np.norm(A*u-b)^2;
35 # print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid);
36 #
37 # # Refinement using fzero
38 # lambda10 = [lambda1 lambda1*alpha]; # interval with zero-crossing
39 # f = @(lambda1) np.norm(A*irls_dn(A,b,p,lambda1) - b)^2 - G;
40 # opts = optimset('fzero');
41 # # opts.Display = 'iter';
42 # opts.Display = 'none';
43 # opts.TolX = 0.01*lambda1;
44 # lambda1 = fzero(f,lambda10,opts);
45 #
46 # u = irls_dn(A,b,p,lambda1);
47 #
48 #
49 # return u;
50
Show file before | Show file after | Hide comments
@@ -1,430 +1,468
1 1 #!/usr/bin/env python
2 2
3 3 #----------------------------------------------------------
4 4 # Original MATLAB code developed by Brian Harding
5 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 11 import math
12 #import cmath
13 #import scipy
14 #import matplotlib
15 12 import numpy as np
16 #from numpy import linalg
17 13 import matplotlib.pyplot as plt
18 14 from scipy import linalg
19 15 import time
20 16 from y_hysell96 import*
21 17 from deb4_basis import *
22 18 from modelf import *
23 from scipy.optimize import fsolve
19 #from scipy.optimize import fsolve
24 20 from scipy.optimize import root
25 21 import pywt
22 from irls_dn2 import *
26 23
27 24
28 25 ## Calculate Forward Model
29 26 lambda1 = 6.0
30 27 k = 2*math.pi/lambda1
31 28
32 29 ## Calculate Magnetic Declination
33 30
34 31 # [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this
35 32
36 33 # or calculate it with the above function
37 34 dec = -1.24
38 35
39 36 # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
40 37 rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
41 38 ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
42 39
43 40 antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240],
44 41 [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
45 42
46 43 plt.figure(1)
47 44 plt.plot(rx, ry, 'ro')
48 45 plt.draw()
49 46
50 47 # Jicamarca is nominally at a 45 degree angle
51 48 theta = 45 - dec;
52 49
53 50 # Rotation matrix from antenna coord to magnetic coord (East North)
54 51 theta_rad = math.radians(theta) # trig functions take radians as argument
55 52 val1 = float( math.cos(theta_rad) )
56 53 val2 = float( math.sin(theta_rad) )
57 54 val3 = float( -1*math.sin(theta_rad))
58 55 val4 = float( math.cos(theta_rad) )
59 56
60 57 # Rotation matrix from antenna coord to magnetic coord (East North)
61 58 R = np.array( [[val1, val3], [val2, val4]] );
62 59
63 60 # Rotate antenna positions to magnetic coord.
64 61 AR = np.dot(R.T, antpos);
65 62
66 63 # Only take the East component
67 64 r = AR[0,:]
68 65 r.sort() # ROW VECTOR?
69 66
70 67 # Truth model (high and low resolution)
71 68 Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
72 69 thbound = 9.0/180*math.pi; # the width of the domain in angle space
73 70 thetat = np.linspace(-thbound, thbound,Nt) # image domain
74 71 thetat = np.transpose(thetat) # transpose # FUNCIONA??????????????????????????????
75 72 Nr = (256.0); # number of pixels in reconstructed image: low res
76 73 thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain
77 74 thetar = np.transpose(thetar) #transpose # FUNCIONA?????????????????????????????
78 75
79 76 # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
80 77 # background constant b.
81 78
82 79 # Triple Gaussian
83 80 # a = np.array([3, 5, 2]);
84 81 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
85 82 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
86 83 # b = 0; # background
87 84
88 85 # Double Gaussian
89 86 # a = np.array([3, 5]);
90 87 # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
91 88 # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
92 89 # b = 0; # background
93 90
94 91 # Single Gaussian
95 92 a = np.array( [3] );
96 93 mu = np.array( [-3.0/180*math.pi] )
97 94 sig = np.array( [2.0/180*math.pi] )
98 95 b = 0;
99 96
100 97 fact = np.zeros(shape=(Nt,1));
101 98 factr = np.zeros(shape=(Nr,1));
102 99
103 100 for i in range(0, a.size):
104 101 temp = (-(thetat-mu[i])**2/(sig[i]**2))
105 102 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
106 103 for j in range(0, temp.size):
107 104 fact[j] = fact[j] + a[i]*math.exp(temp[j]);
108 105 for m in range(0, tempr.size):
109 106 factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
110 107 fact = fact + b;
111 108 factr = factr + b;
112 109
113 110 # # model for f: Square pulse
114 111 # for j in range(0, fact.size):
115 112 # if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi):
116 113 # fact[j] = 0
117 114 # else:
118 115 # fact[j] = 1
119 116 # for k in range(0, factr.size):
120 117 # if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
121 118 # fact[k] = 0
122 119 # else:
123 120 # fact[k] = 1
124 121 #
125 122 #
126 123 # # model for f: triangle pulse
127 124 # mu = -1.0/180*math.pi;
128 125 # sig = 5.0/180*math.pi;
129 126 # wind1 = theta > mu-sig and theta < mu;
130 127 # wind2 = theta < mu+sig and theta > mu;
131 128 # fact = wind1 * (theta - (mu - sig));
132 129 # factr = wind1 * (thetar - (mu - sig));
133 130 # fact = fact + wind2 * (-(theta-(mu+sig)));
134 131 # factr = factr + wind2 * (-(thetar-(mu+sig)));
135 132
136 133
137 134 # fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1
138 135 I = sum(fact)[0];
139 136 fact = fact/I; # normalize to sum(f)==1
140 137 factr = factr/I; # normalize to sum(f)==1
141 #plot(thetat,fact,'r'); hold on;
142 #plot(thetar,factr,'k.'); hold off;
138 #plt.figure()
139 #plt.plot(thetat,fact,'r');
140 #plt.plot(thetar,factr,'k.');
143 141 #xlim([min(thetat) max(thetat)]);
144 142
145 143 #x = np.linspace(thetat.min(), thetat.max) ????
146 144 #for i in range(0, thetat.size):
147 145 plt.figure(2)
148 146 plt.plot(thetat, fact, 'r--')
149 147 plt.plot(thetar, factr, 'ro')
150 148 plt.draw()
151 149 # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
152 150
153 151
154 152 ##
155 153 # Control the type and number of inversions with:
156 154 # SNRdBvec: the SNRs that will be used.
157 155 # NN: the number of trials for each SNR
158 156
159 157 #SNRdBvec = np.linspace(5,20,10);
160 158 SNRdBvec = np.array([15]);
161 159 NN = 1; # number of trial at each SNR
162 160
163 161 # if using vector arguments should be: (4,SNRdBvec.size,NN)
164 162 corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
165 163 corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
166 164 rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
167 165
168 166 for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
169 167 for Ni in range(0, NN):
170 168 SNRdB = SNRdBvec[snri];
171 169 SNR = 10**(SNRdB/10.0);
172 170
173 171 # Calculate cross-correlation matrix (Fourier components of image)
174 172 # This is an inefficient way to do this.
175 173 R = np.zeros(shape=(r.size, r.size), dtype=object);
176 174
177 175 for i1 in range(0, r.size):
178 176 for i2 in range(0,r.size):
179 177 R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat))))
180 178 R[i1,i2] = sum(R[i1,i2])
181 179
182 180 # Add uncertainty
183 181 # This is an ad-hoc way of adding "noise". It models some combination of
184 182 # receiver noise and finite integration times. We could use a more
185 183 # advanced model (like in Yu et al 2000) in the future.
186 184
187 185 # This is a way of adding noise while maintaining the
188 186 # positive-semi-definiteness of the matrix.
189 187
190 188 U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
191 189
192 190 sigma_noise = (np.linalg.norm(U,'fro')/SNR);
193 191
194 192 temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
195 193 temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
196 194 temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's
197 195 temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
198 196
199 197 nz = np.multiply(temp4, temp3)
200
198
199 #---------------------- Eliminar esto:------------------------------------------
200 #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int))
201 201 #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
202
203
204 202 #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
205 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))));
206
204 #--------------------------------------------------------------------------------
205
207 206 Unz = U + nz;
208 207 Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
209 208 plt.figure(3);
210 209 plt.pcolor(abs(Rnz));
211 210 plt.colorbar();
212 211
213 212 # Fourier Inversion ###################
214 213 f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
215 214
216 215 for i in range(0, thetar.size):
217 216 th = thetar[i];
218 217 w = np.exp(1j*k*np.dot(r,np.sin(th)));
219 218
220 219 temp = np.dot(w.T.conj(),U)
221 220
222 221 f_fourier[i] = np.dot(temp, w);
223 222
224 223 f_fourier = f_fourier.real; # get rid of numerical imaginary noise
225 224
226 225 #print f_fourier
227 226
228 227
229 228 # Capon Inversion ######################
230 229
231 230 f_capon = np.zeros(shape=(Nr,1));
232 231
233 #tic_capon = time.time();
232 tic_capon = time.time();
234 233
235 234 for i in range(0, thetar.size):
236 235 th = thetar[i];
237 236 w = np.exp(1j*k*np.dot(r,np.sin(th)));
238 237 f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
239 238
240 #toc_capon = time.time()
241
242 # elapsed_time_capon = toc_capon - tic_capon;
239
240 toc_capon = time.time()
241
242 elapsed_time_capon = toc_capon - tic_capon;
243 243
244 244 f_capon = f_capon.real; # get rid of numerical imaginary noise
245 245
246 246 # MaxEnt Inversion #####################
247 247
248 248 # create the appropriate sensing matrix (split into real and imaginary # parts)
249 249 M = (r.size-1)*(r.size);
250 250 Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix
251 251 Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction
252 252
253 253 # need to re-index our measurements from matrix R into vector g
254 254 g = np.zeros(shape=(M,1));
255 255 gnz = np.zeros(shape=(M,1)); # noisy version of g
256 256
257 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 261 # build H
262 262 for i1 in range(0, r.size):
263 263 for i2 in range(i1+1, r.size):
264 264 idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward
265 265 idx1 = 2*idx; # because index starts at 0
266 266 idx2 = 2*idx+1;
267 267 Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T;
268 268 Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T;
269 269 Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
270 270 Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt;
271 271 g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later
272 272 g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again
273 273 gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt;
274 274 gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
275 275
276 276 # inversion
277 277 F = Nr/Nt; # normalization
278 278 sigma = 1; # set to 1 because the difference is accounted for in G
279 279
280 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
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();
282 284
283 285 lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
284 286
287 toc_maxent = time.time()
288 elapsed_time_maxent = toc_maxent - tic_maxent;
289
285 290 # Whitened solution
286 291 def myfun(lambda1):
287 292 return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
288 293
289 294 tic_maxEnt = time.time();
290 295
291 296 #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
292 297 lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14);
293 298
294 299 #print lambda1
295 300 #print lambda1.x
296 301
297 302 lambda1 = lambda1.x;
298 303
299 304 toc_maxEnt = time.time();
300 305 f_maxent = modelf(lambda1, Hr, F);
301 306 ystar = myfun(lambda1);
302 307 Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G));
303 308 ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
304 309 es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
305 310 chi2 = np.sum((es/sigma)**2);
306 311
307
308
312
309 313 # CS inversion using irls ########################
310 314
311 315 # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
312 316
313 #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
314
315 wavelet1 = pywt.Wavelet('db4')
316 Phi, Psi, x = wavelet1.wavefun(level=3)
317 Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
318
319 # REMOVE THIS?--------------------------------
320 #wavelet1 = pywt.Wavelet('db4')
321 #Phi, Psi, x = wavelet1.wavefun(level=3)
322 # --------------------------------------------
317 323
318 324 # add "sum to 1" constraint
319 325 H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
320 326 N_temp = np.array([[Nr/Nt]]);
321 327 g2 = np.concatenate( (gnz, N_temp), axis=0 );
322
323
324 #s = irls_dn2(H2*Psi,g2,0.5,G);
328 H2 = H2.T.conj();
329
330 print 'H2 shape', H2.shape
331 print 'Psi shape', Psi.shape
332
333 s = irls_dn2(H2*Psi,g2,0.5,G);
325 334 # f_cs = Psi*s;
326 335 #
327 336 # # plot
328 337 # plot(thetar,f_cs,'r.-');
329 338 # hold on;
330 339 # plot(thetat,fact,'k-');
331 340 # hold off;
332 341
333 342
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));
343 # # # Scaling and shifting
344 # # # Only necessary for capon solution
347 345
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
346
347 f_capon = f_capon/np.max(f_capon)*np.max(fact);
348
349
350 ### analyze stuff ######################
351 # calculate MSE
352 rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
353 rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
354 rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
355 #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
356
357
358 relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
359 relrmse_capon = rmse_capon / np.linalg.norm(fact);
360 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
361 #relrmse_cs = rmse_cs / np.norm(fact);
362
363 # To be able to perform dot product (align matrices) done below within the dot calculations
364
365
366 #f_fourier = f_fourier.T.conj()
367 #f_capon = f_capon.T.conj()
368 #f_maxent = f_maxent.T.conj()
369
370 #factr = factr.T.conj()
371
372 # calculate correlation
373
374 corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
375 corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
376 corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
377 #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
378
379
380 # calculate centered correlation
381 f0 = factr - np.mean(factr);
382 f1 = f_fourier - np.mean(f_fourier);
383
384 corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
385 f1 = f_capon - np.mean(f_capon);
386 corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
387 f1 = f_maxent - np.mean(f_maxent);
388 corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
389 #f1 = f_cs - mean(f_cs);
390 #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
391
392
393
394 # # # plot stuff #########################
395
396 #---- Capon----
397 plt.figure(4)
398 plt.subplot(2, 1, 1)
399 plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
400 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
401 plt.ylabel('Power (arbitrary units)')
402 plt.legend(loc='upper right')
403
404 # formatting y-axis
405 locs,labels = plt.yticks()
406 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
407 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
408
409
410 #---- MaxEnt----
411 plt.subplot(2, 1, 2)
412 plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
413 plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
414 plt.ylabel('Power (arbitrary units)')
415 plt.legend(loc='upper right')
416
417 # formatting y-axis
418 locs,labels = plt.yticks()
419 plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
420 plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
421
422 plt.show()
423
424
425 # # PLOT PARA COMPRESSED SENSING
426 # #
427 # # subplot(3,1,3);
428 # # plot(180/pi*thetar,f_cs,'r-');
429 # # hold on;
430 # # plot(180/pi*thetat,fact,'k--');
431 # # hold off;
432 # # ylim([min(f_cs) 1.1*max(fact)]);
433 # # # title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
434 # # # title 'Compressed Sensing - Debauchies Wavelets'
435 # # xlabel 'Degrees'
436 # # ylabel({'Power';'(arbitrary units)'})
437 # # legend('Comp. Sens.','Truth');
438 # #
439 # # # set(gcf,'Position',[749 143 528 881]); # CSL
440 # # # set(gcf,'Position',[885 -21 528 673]); # macbook
441 # # pause(0.01);
442
443
444 # # Store Results
445 corr[0, snri, Ni] = corr_fourier;
446 corr[1, snri, Ni] = corr_capon;
447 corr[2, snri, Ni] = corr_maxent;
448 #corr[3, snri, Ni] = corr_cs;
449
450 rmse[0,snri,Ni] = relrmse_fourier;
451 rmse[1,snri,Ni] = relrmse_capon;
452 rmse[2,snri,Ni] = relrmse_maxent;
453 #rmse[3,snri,Ni] = relrmse_cs;
454
455 corrc[0,snri,Ni] = corrc_fourier;
456 corrc[1,snri,Ni] = corrc_capon;
457 corrc[2,snri,Ni] = corrc_maxent;
458 #corrc[3,snri,Ni] = corrc_cs;
459
460
461 print 'Capon:\t', elapsed_time_capon, 'sec';
462 print 'Maxent:\t',elapsed_time_maxent, 'sec';
463 #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
464
465 print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
466
467 print corr
468 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now