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 21) +++ b/trunk/src/src/RICS_paper_compare_noise.py (revision 22) @@ -2,127 +2,133 @@ #---------------------------------------------------------- # Original MATLAB code developed by Brian Harding -# Rewritten in python by Yolian Amaro +# Rewritten in Python by Yolian Amaro # Python version 2.7 # May 15, 2014 # Jicamarca Radio Observatory #---------------------------------------------------------- -import math -import numpy as np +import time import matplotlib.pyplot as plt -from scipy import linalg -import time +from scipy.optimize import root + from y_hysell96 import* from deb4_basis import * from modelf import * +from irls_dn2 import * #from scipy.optimize import fsolve -from scipy.optimize import root -import pywt -from irls_dn2 import * - + + +#------------------------------------------------------------------------------------------------- +# Set parameters +#------------------------------------------------------------------------------------------------- ## Calculate Forward Model lambda1 = 6.0 -k = 2*math.pi/lambda1 - -## Calculate Magnetic Declination - -# [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012); check this - -# or calculate it with the above function +k = 2*np.pi/lambda1 + +## Magnetic Declination dec = -1.24 - -# loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt() + +## Loads Jicamarca antenna positions +antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False) + +## rx and ry -- for plotting purposes rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] ) ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] ) -antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240], - [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] ) - +## Plot of antenna positions plt.figure(1) plt.plot(rx, ry, 'ro') plt.draw() -# Jicamarca is nominally at a 45 degree angle +## Jicamarca is nominally at a 45 degree angle theta = 45 - dec; +## Rotation matrix from antenna coord to magnetic coord (East North) +theta_rad = np.radians(theta) # trig functions take radians as argument +val1 = float( np.cos(theta_rad) ) +val2 = float( np.sin(theta_rad) ) +val3 = float( -1*np.sin(theta_rad)) +val4 = float( np.cos(theta_rad) ) + # Rotation matrix from antenna coord to magnetic coord (East North) -theta_rad = math.radians(theta) # trig functions take radians as argument -val1 = float( math.cos(theta_rad) ) -val2 = float( math.sin(theta_rad) ) -val3 = float( -1*math.sin(theta_rad)) -val4 = float( math.cos(theta_rad) ) - -# Rotation matrix from antenna coord to magnetic coord (East North) -R = np.array( [[val1, val3], [val2, val4]] ); +R = np.array( [[val1, val3], [val2, val4]] ) # Rotate antenna positions to magnetic coord. -AR = np.dot(R.T, antpos); +AR = np.dot(R.T, antpos) # Only take the East component r = AR[0,:] -r.sort() # ROW VECTOR? +r.sort() # Truth model (high and low resolution) -Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution -thbound = 9.0/180*math.pi; # the width of the domain in angle space -thetat = np.linspace(-thbound, thbound,Nt) # image domain -thetat = np.transpose(thetat) # transpose # FUNCIONA?????????????????????????????? -Nr = (256.0); # number of pixels in reconstructed image: low res +Nt = (1024.0)*(16.0) # number of pixels in truth image: high resolution +thbound = 9.0/180*np.pi # the width of the domain in angle space +thetat = np.linspace(-thbound, thbound,Nt) # image domain +thetat = thetat.T # transpose +Nr = (256.0) # number of pixels in reconstructed image: low res thetar = np.linspace(-thbound, thbound,Nr) # reconstruction domain -thetar = np.transpose(thetar) #transpose # FUNCIONA????????????????????????????? - -# Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and -# background constant b. +thetar = thetar.T # transpose + + +#------------------------------------------------------------------------------------------------- +# Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b. +#------------------------------------------------------------------------------------------------- # Triple Gaussian # a = np.array([3, 5, 2]); -# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]); -# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]); -# b = 0; # background +# mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi, 7.0/180*np.pi]); +# sig = np.array([2.0/180*np.pi, 1.5/180*np.pi, 0.3/180*np.pi]); +# b = 0; # background # Double Gaussian # a = np.array([3, 5]); -# mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]); -# sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]); -# b = 0; # background +# mu = np.array([-5.0/180*np.pi, 2.0/180*np.pi]); +# sig = np.array([2.0/180*np.pi, 1.5/180*np.pi]); +# b = 0; # background # Single Gaussian -a = np.array( [3] ); -mu = np.array( [-3.0/180*math.pi] ) -sig = np.array( [2.0/180*math.pi] ) -b = 0; - -fact = np.zeros(shape=(Nt,1)); -factr = np.zeros(shape=(Nr,1)); - +a = np.array( [3] ) +mu = np.array( [-3.0/180*np.pi] ) +sig = np.array( [2.0/180*np.pi] ) +b = 0 + +# Empty matrices for factors +fact = np.zeros(shape=(Nt,1)) +factr = np.zeros(shape=(Nr,1)) + +# DFT Kernels for i in range(0, a.size): temp = (-(thetat-mu[i])**2/(sig[i]**2)) tempr = (-(thetar-mu[i])**2/(sig[i]**2)) for j in range(0, temp.size): - fact[j] = fact[j] + a[i]*math.exp(temp[j]); + fact[j] = fact[j] + a[i]*np.exp(temp[j]); for m in range(0, tempr.size): - factr[m] = factr[m] + a[i]*math.exp(tempr[m]); + factr[m] = factr[m] + a[i]*np.exp(tempr[m]); + fact = fact + b; factr = factr + b; -# # model for f: Square pulse +# #------------------------------------------------------------------------------------------------- +# # Model for f: Square pulse +# #------------------------------------------------------------------------------------------------- # for j in range(0, fact.size): -# if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi): +# if (theta > -5.0/180*np.pi and theta < 2.0/180*np.pi): # fact[j] = 0 # else: # fact[j] = 1 # for k in range(0, factr.size): -# if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi): -# fact[k] = 0 +# if (thetar[k] > -5.0/180*np.pi and thetar[k] < 2/180*np.pi): +# factr[k] = 0 # else: -# fact[k] = 1 -# -# -# # model for f: triangle pulse -# mu = -1.0/180*math.pi; -# sig = 5.0/180*math.pi; +# factr[k] = 1 + +# #------------------------------------------------------------------------------------------------- +# # Model for f: Triangle pulse +# #------------------------------------------------------------------------------------------------- +# mu = -1.0/180*np.pi; +# sig = 5.0/180*np.pi; # wind1 = theta > mu-sig and theta < mu; # wind2 = theta < mu+sig and theta > mu; # fact = wind1 * (theta - (mu - sig)); @@ -131,45 +137,43 @@ # factr = factr + wind2 * (-(thetar-(mu+sig))); -# fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1 +# fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1 + I = sum(fact)[0]; -fact = fact/I; # normalize to sum(f)==1 -factr = factr/I; # normalize to sum(f)==1 -#plt.figure() -#plt.plot(thetat,fact,'r'); -#plt.plot(thetar,factr,'k.'); -#xlim([min(thetat) max(thetat)]); - -#x = np.linspace(thetat.min(), thetat.max) ???? -#for i in range(0, thetat.size): +fact = fact/I; # normalize to sum(f)==1 +factr = factr/I; # normalize to sum(f)==1 + +# Plot Gaussian pulse(s) plt.figure(2) plt.plot(thetat, fact, 'r--') plt.plot(thetar, factr, 'ro') plt.draw() -# xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO - - -## + + +#------------------------------------------------------------------------------------------------- # Control the type and number of inversions with: # SNRdBvec: the SNRs that will be used. # NN: the number of trials for each SNR +#------------------------------------------------------------------------------------------------- #SNRdBvec = np.linspace(5,20,10); -SNRdBvec = np.array([15]); -NN = 1; # number of trial at each SNR - -# if using vector arguments should be: (4,SNRdBvec.size,NN) -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) - -for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector +SNRdBvec = np.array([15]); # 15 dB +NN = 1; # number of trials at each SNR + +# Statistics simulation (correlation, root mean square) +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) + +# For each SNR and trial +for snri in range(0, SNRdBvec.size): + SNRdB = SNRdBvec[snri]; + SNR = 10**(SNRdB/10.0); + for Ni in range(0, NN): - SNRdB = SNRdBvec[snri]; - SNR = 10**(SNRdB/10.0); - # Calculate cross-correlation matrix (Fourier components of image) # This is an inefficient way to do this. + R = np.zeros(shape=(r.size, r.size), dtype=object); for i1 in range(0, r.size): @@ -185,48 +189,42 @@ # This is a way of adding noise while maintaining the # positive-semi-definiteness of the matrix. - U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R + U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R sigma_noise = (np.linalg.norm(U,'fro')/SNR); - temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5) - temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5) - temp3 = ((abs(U) > 0).astype(float)) # upper triangle of 1's - temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0) +# 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(temp4, temp3) - - #---------------------- Eliminar esto:------------------------------------------ - #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int)) - #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int)) - #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 )); - #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)))); - #-------------------------------------------------------------------------------- - Unz = U + nz; - Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R + Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R plt.figure(3); plt.pcolor(abs(Rnz)); plt.colorbar(); - - # Fourier Inversion ################### + + #------------------------------------------------------------------------------------------------- + # Fourier Inversion + #------------------------------------------------------------------------------------------------- f_fourier = np.zeros(shape=(Nr,1), dtype=complex); for i in range(0, thetar.size): th = thetar[i]; - w = np.exp(1j*k*np.dot(r,np.sin(th))); - + w = np.exp(1j*k*np.dot(r,np.sin(th))); temp = np.dot(w.T.conj(),U) - f_fourier[i] = np.dot(temp, w); - f_fourier = f_fourier.real; # get rid of numerical imaginary noise - - #print f_fourier - - - # Capon Inversion ###################### - + f_fourier = f_fourier.real; # get rid of numerical imaginary noise + + + #------------------------------------------------------------------------------------------------- + # Capon Inversion + #------------------------------------------------------------------------------------------------- f_capon = np.zeros(shape=(Nr,1)); tic_capon = time.time(); @@ -236,151 +234,149 @@ w = np.exp(1j*k*np.dot(r,np.sin(th))); f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real) - toc_capon = time.time() elapsed_time_capon = toc_capon - tic_capon; - f_capon = f_capon.real; # get rid of numerical imaginary noise - - # MaxEnt Inversion ##################### - - # create the appropriate sensing matrix (split into real and imaginary # parts) + f_capon = f_capon.real; # get rid of numerical imaginary noise + + + #------------------------------------------------------------------------------------------------- + # MaxEnt Inversion + #------------------------------------------------------------------------------------------------- + + # Create the appropriate sensing matrix (split into real and imaginary # parts) M = (r.size-1)*(r.size); - Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix - Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction + Ht = np.zeros(shape=(M,Nt)); # "true" sensing matrix + Hr = np.zeros(shape=(M,Nr)); # approximate sensing matrix for reconstruction - # need to re-index our measurements from matrix R into vector g + # Need to re-index our measurements from matrix R into vector g g = np.zeros(shape=(M,1)); - gnz = np.zeros(shape=(M,1)); # noisy version of g + gnz = np.zeros(shape=(M,1)); # noisy version of g - # triangular indexing to perform this re-indexing + # Triangular indexing to perform this re-indexing T = np.ones(shape=(r.size,r.size)); [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing - # build H + # 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 - idx1 = 2*idx; # because index starts at 0 + idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0]; # kind of awkward + idx1 = 2*idx; idx2 = 2*idx+1; - Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T; - Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T; - Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt; - Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T*Nr/Nt; - g[idx1] = (R[i1,i2]).real*Nr/Nt; # check this again later - g[idx2] = (R[i1,i2]).imag*Nr/Nt; # check again + Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj(); + Hr[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj(); + Ht[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt; + Ht[idx2,:] = np.sin(k*(r[i1]-r[i2])*np.sin(thetat)).T.conj()*Nr/Nt; + g[idx1] = (R[i1,i2]).real*Nr/Nt; + g[idx2] = (R[i1,i2]).imag*Nr/Nt; gnz[idx1] = (Rnz[i1,i2]).real*Nr/Nt; gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt; - # inversion - F = Nr/Nt; # normalization - sigma = 1; # set to 1 because the difference is accounted for in G - - ##### ADD *10 for consistency with old model, NEED TO VERIFY THIS!!!!? line below - G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2 - - tic_maxent = time.time(); - - lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything) + # Inversion + F = Nr/Nt; # normalization + sigma = 1; # set to 1 because the difference is accounted for in G + + G = np.linalg.norm(g-gnz)**2 ; # pretend we know in advance the actual value of chi^2 + lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything) - toc_maxent = time.time() - elapsed_time_maxent = toc_maxent - tic_maxent; # Whitened solution def myfun(lambda1): return y_hysell96(lambda1,gnz,sigma,F,G,Hr); - tic_maxEnt = time.time(); + + tic_maxEnt = time.time(); # start time maxEnt - #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000); lambda1 = root(myfun,lambda0, method='krylov', tol=1e-14); - #print lambda1 - #print lambda1.x - + toc_maxEnt = time.time() + elapsed_time_maxent = toc_maxEnt - tic_maxEnt; + + # Solution lambda1 = lambda1.x; - toc_maxEnt = time.time(); f_maxent = modelf(lambda1, Hr, F); ystar = myfun(lambda1); - Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G)); + 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; # should be same as ep + es = np.dot(Hr, f_maxent) - gnz; # should be same as ep chi2 = np.sum((es/sigma)**2); - - # --------- CS inversion using Iteratively Reweighted Least Squares (IRLS) ------------- - + + #------------------------------------------------------------------------------------------------- + # CS inversion using Iteratively Reweighted Least Squares (IRLS) + #------------------------------------------------------------------------------------------------- # (Use Nr, thetar, gnz, and Hr from MaxEnt above) - + Psi = deb4_basis(Nr); - - # REMOVE THIS-------------------------------- - #wavelet1 = pywt.Wavelet('db4') - #Phi, Psi, x = wavelet1.wavefun(level=3) - # -------------------------------------------- - + + # ------------Remove this------------------------------------------- + # wavelet1 = pywt.Wavelet('db4') + # Phi, Psi, x = wavelet1.wavefun(level=3) + #------------------------------------------------------------------- + # add "sum to 1" constraint H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 ); - N_temp = np.array([[Nr/Nt]]); - g2 = np.concatenate( (gnz, N_temp), axis=0 ); - - #H2 = H2.T.conj(); - - #Psi = Psi.T.conj(); # to align matrices - - ####print 'H2 shape', H2.shape - #####print 'Psi shape', Psi.shape - - A = np.dot(H2,Psi); - + g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 ); + + tic_cs = time.time(); + + +# plt.figure(4) +# plt.imshow(Psi)#, interpolation='nearest') +# #plt.xticks([]); plt.yticks([]) +# plt.show() + + # Inversion s = irls_dn2(np.dot(H2,Psi),g2,0.5,G); - # f_cs = Psi*s; - # - # # plot - # plot(thetar,f_cs,'r.-'); - # hold on; - # plot(thetat,fact,'k-'); - # hold off; - - - # # # Scaling and shifting - # # # Only necessary for capon solution + + #print s + + # Brightness function + f_cs = np.dot(Psi,s); + + toc_cs = time.time() + elapsed_time_cs = toc_cs - tic_cs; + + # Plot + plt.figure(4) + plt.plot(thetar,f_cs,'r.-'); + plt.plot(thetat,fact,'k-'); + + + #------------------------------------------------------------------------------------------------- + # Scaling and shifting + # (Only necessary for Capon solution) + #------------------------------------------------------------------------------------------------- f_capon = f_capon/np.max(f_capon)*np.max(fact); - ### analyze stuff ###################### - # calculate MSE + #------------------------------------------------------------------------------------------------- + # Analyze stuff + #------------------------------------------------------------------------------------------------- + + # Calculate MSE rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2)); rmse_capon = np.sqrt(np.mean((f_capon - factr)**2)); rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2)); - #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2)); + rmse_cs = np.sqrt(np.mean((f_cs - factr)**2)); relrmse_fourier = rmse_fourier / np.linalg.norm(fact); relrmse_capon = rmse_capon / np.linalg.norm(fact); relrmse_maxent = rmse_maxent / np.linalg.norm(fact); - #relrmse_cs = rmse_cs / np.norm(fact); + relrmse_cs = rmse_cs / np.linalg.norm(fact); - # To be able to perform dot product (align matrices) done below within the dot calculations - - - #f_fourier = f_fourier.T.conj() - #f_capon = f_capon.T.conj() - #f_maxent = f_maxent.T.conj() - - #factr = factr.T.conj() - - # calculate correlation - + + # Calculate correlation corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr)); corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr)); corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr)); - #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr)); + corr_cs = np.dot(f_cs.T.conj(),factr) / (np.linalg.norm(f_cs)*np.linalg.norm(factr)); - # calculate centered correlation + # Calculate centered correlation f0 = factr - np.mean(factr); f1 = f_fourier - np.mean(f_fourier); @@ -389,18 +385,19 @@ corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); f1 = f_maxent - np.mean(f_maxent); corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); - #f1 = f_cs - mean(f_cs); - #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1)); - - - - # # # plot stuff ######################### - + f1 = f_cs - np.mean(f_cs); + corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1)); + + + #------------------------------------------------------------------------------------------------- + # Plot stuff + #------------------------------------------------------------------------------------------------- + #---- Capon---- - plt.figure(4) - plt.subplot(2, 1, 1) - plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon') - plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') + plt.figure(5) + plt.subplot(3, 1, 1) + plt.plot(180/np.pi*thetar, f_capon, 'r', label='Capon') + plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth') plt.ylabel('Power (arbitrary units)') plt.legend(loc='upper right') @@ -411,9 +408,9 @@ #---- MaxEnt---- - plt.subplot(2, 1, 2) - plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt') - plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth') + plt.subplot(3, 1, 2) + plt.plot(180/np.pi*thetar, f_maxent, 'r', label='MaxEnt') + plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth') plt.ylabel('Power (arbitrary units)') plt.legend(loc='upper right') @@ -422,7 +419,18 @@ plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) - plt.show() + + #---- Compressed Sensing---- + plt.subplot(3, 1, 3) + plt.plot(180/np.pi*thetar, f_cs, 'r', label='CS') + plt.plot(180/np.pi*thetat,fact, 'k--', label='Truth') + plt.ylabel('Power (arbitrary units)') + plt.legend(loc='upper right') + + # formatting y-axis + locs,labels = plt.yticks() + plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4)) + plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes) # # PLOT PARA COMPRESSED SENSING @@ -448,24 +456,71 @@ corr[0, snri, Ni] = corr_fourier; corr[1, snri, Ni] = corr_capon; corr[2, snri, Ni] = corr_maxent; - #corr[3, snri, Ni] = corr_cs; + corr[3, snri, Ni] = corr_cs; rmse[0,snri,Ni] = relrmse_fourier; rmse[1,snri,Ni] = relrmse_capon; rmse[2,snri,Ni] = relrmse_maxent; - #rmse[3,snri,Ni] = relrmse_cs; + rmse[3,snri,Ni] = relrmse_cs; corrc[0,snri,Ni] = corrc_fourier; corrc[1,snri,Ni] = corrc_capon; corrc[2,snri,Ni] = corrc_maxent; - #corrc[3,snri,Ni] = corrc_cs; + corrc[3,snri,Ni] = corrc_cs; print 'Capon:\t', elapsed_time_capon, 'sec'; print 'Maxent:\t',elapsed_time_maxent, 'sec'; - #print 'CS:\t%3.3f sec\n',elapsed_time_cs; + print 'CS:\t',elapsed_time_cs, 'sec'; print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN); print corr - \ No newline at end of file + + print corr.shape + + + ## Analyze and plot statistics + + metric = corr; # set this to rmse, corr, or corrc + + # Remove outliers (this part was experimental and wasn't used in the paper) + # nsig = 3; + # for i = 1:4 + # for snri = 1:length(SNRdBvec) + # av = mean(met(i,snri,:)); + # s = std(met(i,snri,:)); + # idx = abs(met(i,snri,:) - av) > nsig*s; + # met(i,snri,idx) = nan; + # if sum(idx)>0 + # fprintf('i=%i, snr=%i, %i/%i pts removed\n',... + # i,round(SNRdBvec(snri)),sum(idx),length(idx)); + # end + # end + # end + + # Avg ignoring NaNs +def nanmean(data, **args): + return numpy.ma.filled(numpy.ma.masked_array(data,numpy.isnan(data)).mean(**args), fill_value=numpy.nan) + +# ave = np.zeros(shape=(4)) +# +# ave[0] = nanmean(metric, axis=0); +# ave[1] = nanmean(metric, axis=1); +# ave[2] = nanmean(metric, axis=2); +# ave[3] = nanmean(metric, axis=3); + +#print ave +plt.figure(6); +f = plt.scatter(SNRdBvec, corr[0], marker='+', color='b', s=60); # Fourier +c = plt.scatter(SNRdBvec, corr[1], marker='o', color= 'c', s=60); # Capon +me= plt.scatter(SNRdBvec, corr[2], marker='s', color= 'y', s=60); # MaxEnt +cs= plt.scatter(SNRdBvec, corr[3], marker='*', color='r', s=60); # Compressed Sensing + + +plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right') +plt.xlabel('SNR') +plt.ylabel('Correlation with Truth') + +plt.show() + Index: trunk/src/src/afb.py =================================================================== diff --git a/trunk/src/src/afb.py b/trunk/src/src/afb.py --- a/trunk/src/src/afb.py (revision 21) +++ b/trunk/src/src/afb.py (revision 22) @@ -50,9 +50,11 @@ lo = upfirdn(x, af[:,0], 1, 2); - # VERIFY THIS!!!!!!!!!!!! - for i in range(0, L): - lo[i] = lo[N/2+i] + lo[i]; +# # VERIFY THIS!!!!!!!!!!!! +# for i in range(0, L): +# lo[i] = lo[N/2+i] + lo[i]; + + lo[0:L-1] = lo[N/2+np.arange(0,L-1)] + lo[0:L-1] lo = lo[0:N/2]; @@ -60,8 +62,10 @@ # highpass filter hi = upfirdn(x, af[:,1], 1, 2); - for j in range(0, L): - hi[j] = hi[N/2+j] + hi[j]; +# for j in range(0, L): +# hi[j] = hi[N/2+j] + hi[j]; + + hi[0:L-1] = hi[N/2+np.arange(0,L-1)] + hi[0:L-1] hi = hi[0:N/2]; Index: trunk/src/src/deb4_basis.py =================================================================== diff --git a/trunk/src/src/deb4_basis.py b/trunk/src/src/deb4_basis.py --- a/trunk/src/src/deb4_basis.py (revision 21) +++ b/trunk/src/src/deb4_basis.py (revision 22) @@ -9,7 +9,8 @@ from dualfilt1 import * from dualtree import * from idualtree import * - + +# Debauchie 4 Wavelet def deb4_basis(N): Psi = np.zeros(shape=(N,2*N+1)); @@ -17,14 +18,15 @@ J = 4; [Faf, Fsf] = FSfarras(); [af, sf] = dualfilt1(); - + # compute transform of zero vector x = np.zeros(shape=(1,N)); - w = dualtree(x, J, Faf, af); + w = dualtree(x, J, Faf, af); + # Uses both real and imaginary wavelets - for i in range (0, J): - for j in range (0, 1): + for i in range (0, J+1): + for j in range (0, 2): for k in range (0, (w[i][j]).size): w[i][j][0,k] = 1; y = idualtree(w, J, Fsf, sf); Index: trunk/src/src/dualtree.py =================================================================== diff --git a/trunk/src/src/dualtree.py b/trunk/src/src/dualtree.py --- a/trunk/src/src/dualtree.py (revision 21) +++ b/trunk/src/src/dualtree.py (revision 22) @@ -45,7 +45,7 @@ #----------------------------------------# # normalization - x = x/np.sqrt(2); + x = x/np.sqrt(2.0); w = np.zeros(shape=(J+1), dtype=object) Index: trunk/src/src/idualtree.py =================================================================== diff --git a/trunk/src/src/idualtree.py b/trunk/src/src/idualtree.py --- a/trunk/src/src/idualtree.py (revision 21) +++ b/trunk/src/src/idualtree.py (revision 22) @@ -41,6 +41,6 @@ y2 = sfb(y2, w[0][1], Fsf[0,1]); # normalization - y = (y1 + y2)/np.sqrt(2); + y = (y1 + y2)/np.sqrt(2.0); return y 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 21) +++ b/trunk/src/src/irls_dn.py (revision 22) @@ -8,6 +8,7 @@ from scipy import linalg import scipy.sparse as sps import numpy as np +from numpy.linalg import norm def irls_dn(A,b,p,lambda1): @@ -30,6 +31,9 @@ u = np.dot(Q,c); # minimum 2-norm solution I = sps.eye(M); + # Temporary N x N matrix + temp = np.zeros(shape=(N,N)) + #---------- not needed, defined above-------------- # Spacing of floating point numbers #eps = np.spacing(1) @@ -39,36 +43,43 @@ while (eps > 1e-7): epschange = 0; # Loop until it's time to change eps - while (~epschange): + while (not(epschange)): # main loop # u_n = W*A'*(A*W*A'+ lambda*I)^-1 * b # where W = diag(1/w) # where w = (u.^2 + eps).^(p/2-1) # Update - w = (u**2 + eps)**(1-p/2); + w = (u**2 + eps)**(1-p/2.0); + +# #---- Very inefficient- REMOVE THIS PART------ +# k = 0 +# # Sparse matrix +# for i in range (0, N): +# for j in range (0,N): +# if(i==j): +# temp[i,j] = w[k] +# k = k+1 + #-------------------------------------------------- - # Empty temporary N x N matrix - temp = np.zeros(shape=(N,N)) - - k = 0 - # Sparse matrix - for i in range (0, N): - for j in range (0,N): - if(i==j): - temp[i,j] = w[k] - k = k+1 + np.fill_diagonal(temp, w) + #----------------------------------------------- # Compressed Sparse Matrix W = sps.csr_matrix(temp); #Compressed Sparse Row matrix WAT = W*A.T.conj(); - - u_new = np.dot(WAT , linalg.solve(np.dot(A,WAT) + np.dot(lambda1,I), b)); + + #print "WAT", WAT.shape + #print "np.dot(A,WAT)", np.dot(A,WAT).shape + #print "np.dot(lambda1,I)", np.dot(lambda1,I).shape + #print "linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b)", linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b).shape + + u_new = np.dot(WAT , linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b)); # See if this subproblem is converging - delu = np.linalg.norm(u_new-u)/np.linalg.norm(u); + delu = norm(u_new-u)/norm(u); epschange = delu < (np.sqrt(eps)/100.0); # Make update @@ -76,8 +87,6 @@ eps = eps/10.0; # decrease eps - # Print info - #print 'eps =',eps; return u 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 21) +++ b/trunk/src/src/irls_dn2.py (revision 22) @@ -7,7 +7,11 @@ from irls_dn import * from scipy.optimize import fsolve import numpy as np -from scipy.optimize import root +from scipy.optimize import * +from dogleg import * +from numpy.linalg import norm + +import matplotlib.pyplot as plt def irls_dn2(A,b,p,G): @@ -27,48 +31,42 @@ alpha = 2.0; # Line search parameter lambda1 = 1e5; # What's a reasonable but safe initial guess? u = irls_dn(A,b,p,lambda1); - fid = np.linalg.norm(np.dot(A,u)-b)**2; - + #print "u\n", u + + fid = norm(np.dot(A,u)-b)**2; + print '----------------------------------\n'; - + while (fid >= G): lambda1 = lambda1 / alpha; # Balance between speed and accuracy u = irls_dn(A,b,p,lambda1); - fid = np.linalg.norm(np.dot(A,u)-b)**2; + fid = norm(np.dot(A,u)-b)**2; print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid; - - # Refinement using fzero + #print u + + # Refinement using fzero/ brentq lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing - + + def myfun(lambda1): - print "A = ", A.shape - print "b = ", b.shape - lambda1 - return np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G; - - #f = lambda lambda1: np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G; NOOOOOO - - -# opts = optimset('fzero'); -# # opts.Display = 'iter'; -# opts.Display = 'none'; -# opts.TolX = 0.01*lambda1; - - #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000); + temp1 = np.dot(A, irls_dn(A,b,p,lambda1)) + temp2 = norm(temp1-b) + temp3 = temp2**2-G + #return np.linalg.norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G; + return temp3 + print "tolerancia=", 0.01*lambda1 - #lambda1 = root(myfun,lambda0, method='krylov', tol=0.01*lambda1); + #lambda1 = root(myfun, lambda1, method='krylov', tol=0.01*lambda1); + #lambda1 = lambda1.x - - print "lamda1=", lambda1 - print "lambda0=", lambda0 - - lambda1 = fsolve(myfun,lambda0); # FALTA OPTIMIZE ESTO - - print "A = ", A.shape - print "b = ", b.shape - print "lambda1=", lambda1.shape - + print "lambda0[0]", lambda0[0] + print "lambda0[1]", lambda0[1] + + lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1) + + print "lambda final=", lambda1 + u = irls_dn(A,b,p,lambda1); Index: trunk/src/src/sfb.py =================================================================== diff --git a/trunk/src/src/sfb.py b/trunk/src/src/sfb.py --- a/trunk/src/src/sfb.py (revision 21) +++ b/trunk/src/src/sfb.py (revision 22) @@ -45,8 +45,12 @@ lo = upfirdn(lo, sf[:,0], 2, 1); hi = upfirdn(hi, sf[:,1], 2, 1); + + lo = lo[0:lo.size-1] + hi = hi[0:hi.size-1] + y = lo + hi; - y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT + y[0:L-2] = y[0:L-2] + y[N+ np.arange(0,L-2)]; #CHECK IF ARANGE IS CORRECT y = y[0:N]; #print 'y en sbf\n', y.shape