Index: trunk/src/src/RICS_paper_compare_noise.py =================================================================== diff --git a/trunk/src/src/RICS_paper_compare_noise.py b/trunk/src/src/RICS_paper_compare_noise.py --- a/trunk/src/src/RICS_paper_compare_noise.py (revision 24) +++ b/trunk/src/src/RICS_paper_compare_noise.py (revision 25) @@ -40,12 +40,17 @@ ## Plot of antenna positions plt.figure(1) plt.plot(rx, ry, 'ro') +plt.title("Jicamarca Antenna Positions") +plt.xlabel("rx") +plt.ylabel("ry") +plt.grid(True) +plt.plot(rx, rx, 'b-') plt.draw() ## Jicamarca is nominally at a 45 degree angle theta = 45 - dec; -## Rotation matrix from antenna coord to magnetic coord (East North) +## Rotation matrix from antenna coord to magnetic coord (East North) -> 45 degrees theta_rad = np.radians(theta) # trig functions take radians as argument val1 = float( np.cos(theta_rad) ) val2 = float( np.sin(theta_rad) ) @@ -58,9 +63,29 @@ # Rotate antenna positions to magnetic coord. AR = np.dot(R.T, antpos) +plt.plot(AR[0], 0*AR[1], 'yo') +plt.title("Jicamarca Antenna Positions") +plt.xlabel("rx") +plt.ylabel("ry") +plt.grid(True) +plt.plot(rx, 0*rx, 'g-') +plt.draw() + +print antpos +print "\n" +print AR + # Only take the East component r = AR[0,:] r.sort() + +# ---NEW--- +r1 = AR[1,:] # longitudes / distances +d = 40 # separation of antennas in meters +theta_w = 15 # angle between antenna and point of observation in degrees +# arg = k*(r + d*np.cos(theta_w)) # this is the exponent argument in the corr integral +# -------- + # Truth model (high and low resolution) Nt = (1024.0)*(16.0) # number of pixels in truth image: high resolution @@ -89,16 +114,16 @@ # b = 0; # background # Single Gaussian -a = np.array( [3] ) -mu = np.array( [-3.0/180*np.pi] ) -sig = np.array( [2.0/180*np.pi] ) -b = 0 +a = np.array( [3] ) # amplitude(s) +mu = np.array( [-3.0/180*np.pi] ) # center +sig = np.array( [2.0/180*np.pi] ) # width +b = 0 # background constant # Empty matrices for factors fact = np.zeros(shape=(Nt,1)) factr = np.zeros(shape=(Nr,1)) -# DFT Kernels +# Constructing Gaussian model for i in range(0, a.size): temp = (-(thetat-mu[i])**2/(sig[i]**2)) tempr = (-(thetar-mu[i])**2/(sig[i]**2)) @@ -106,9 +131,9 @@ fact[j] = fact[j] + a[i]*np.exp(temp[j]); for m in range(0, tempr.size): factr[m] = factr[m] + a[i]*np.exp(tempr[m]); - fact = fact + b; factr = factr + b; + # #------------------------------------------------------------------------------------------------- # # Model for f: Square pulse @@ -135,8 +160,8 @@ # factr = wind1 * (thetar - (mu - sig)); # fact = fact + wind2 * (-(theta-(mu+sig))); # factr = factr + wind2 * (-(thetar-(mu+sig))); - -# fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1 +# +# fact = np.divide(fact, (np.sum(fact)*2*thbound/Nt)); # normalize to integral(f)==1 I = sum(fact)[0]; @@ -147,6 +172,9 @@ plt.figure(2) plt.plot(thetat, fact, 'r--') plt.plot(thetar, factr, 'ro') +plt.title("Truth models") +plt.ylabel("Amplitude") +plt.legend(('High res','Low res'), loc='upper right') plt.draw() @@ -160,7 +188,7 @@ SNRdBvec = np.array([15]); # 15 dB NN = 1; # number of trials at each SNR -# Statistics simulation (correlation, root mean square) +# Statistics simulation (correlation, root mean square error) corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial) @@ -190,22 +218,15 @@ # positive-semi-definiteness of the matrix. U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R - + sigma_noise = (np.linalg.norm(U,'fro')/SNR); - -# temp1 = (-2*np.random.rand(U.shape[0], U.shape[1]) + 2) -# temp2 = 1j*(-2*np.random.rand(U.shape[0], U.shape[1]) + 2) -# temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's -# temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0) -# -# nz = np.multiply(temp4,temp3) - - 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)); + + 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)); Unz = U + nz; Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R plt.figure(3); - plt.pcolor(abs(Rnz)); + plt.pcolor(abs(R)); plt.colorbar(); #------------------------------------------------------------------------------------------------- @@ -260,7 +281,7 @@ # Build H for i1 in range(0, r.size): for i2 in range(i1+1, r.size): - idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward + idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward 1->27 idx1 = 2*idx; idx2 = 2*idx+1; Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj(); @@ -296,11 +317,14 @@ lambda1 = lambda1.x; f_maxent = modelf(lambda1, Hr, F); + + #------ what is this needed for?------- ystar = myfun(lambda1); Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G)); ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda); es = np.dot(Hr, f_maxent) - gnz; chi2 = np.sum((es/sigma)**2); + #-------------------------------------- #------------------------------------------------------------------------------------------------- @@ -326,11 +350,16 @@ f_cs = np.dot(Psi,s); + print "shapes:" + print "Psi ", Psi.shape + print "s ", s.shape + print "f_CS ", f_cs.shape # Plot plt.figure(4) plt.plot(thetar,f_cs,'r.-'); - plt.plot(thetat,fact,'k-'); + plt.plot(thetat,fact,'k-'); + plt.title('Reconstruction with Compressed Sensing') #------------------------------------------------------------------------------------------------- @@ -441,7 +470,7 @@ print 'Maxent:\t',elapsed_time_maxent, 'sec'; print 'CS:\t',elapsed_time_cs, 'sec\n'; - print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN), '\n'; + print (NN*(snri)+Ni+1), '/', (SNRdBvec.size*NN), '\n'; #------------------------------------------------------------------------------------------------- @@ -469,6 +498,8 @@ # Avg ignoring NaNs ave = np.zeros(shape=(4)) +print metric + ave[0] = nanmean(metric[0,:,:]); # Fourier ave[1] = nanmean(metric[1,:,:]); # Capon ave[2] = nanmean(metric[2,:,:]); # MaxEnt Index: trunk/src/src/irls_dn.py =================================================================== diff --git a/trunk/src/src/irls_dn.py b/trunk/src/src/irls_dn.py --- a/trunk/src/src/irls_dn.py (revision 24) +++ b/trunk/src/src/irls_dn.py (revision 25) @@ -23,12 +23,14 @@ # Initialize and precompute: eps = 1e-2; # damping parameter + # Compute the qr factorization of a matrix. + # Factor the matrix A as qr, where Q is orthonormal and R is upper-triangular. [Q,R] = linalg.qr(A.T.conj(), mode='economic'); c = linalg.solve(R.T.conj(),b); # will be used later also u = np.dot(Q,c); # minimum 2-norm solution - I = sps.eye(M); + I = sps.eye(M); # M as sparse matrix # Temporary N x N matrix temp = np.zeros(shape=(N,N)) Index: trunk/src/src/irls_dn2.py =================================================================== diff --git a/trunk/src/src/irls_dn2.py b/trunk/src/src/irls_dn2.py --- a/trunk/src/src/irls_dn2.py (revision 24) +++ b/trunk/src/src/irls_dn2.py (revision 25) @@ -22,8 +22,8 @@ # Line Search - alpha = 2.0; # Line search parameter - lambda1 = 1e5; # What's a reasonable but safe initial guess? + alpha = 2.0; # Line search parameter + lambda1 = 1e5; # What's a reasonable but safe initial guess? u = irls_dn(A,b,p,lambda1); @@ -32,13 +32,12 @@ print '----------------------------------\n'; while (fid >= G): - lambda1 = lambda1 / alpha; # Balance between speed and accuracy + lambda1 = lambda1 / alpha; # Balance between speed and accuracy u = irls_dn(A,b,p,lambda1); fid = norm(np.dot(A,u)-b)**2; print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid; # Refinement using fzero/ brentq - lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing def myfun(lambda1):