jro_soft/compressed_sensing Commit - r17:18

Partial implementation of wavelet transform. Implemented functions: dualfilt1, FSfarras, irls_dn (partial), irls_dn2 (partial), deb4_basis (partial).

yamaro -

r17:18

parent child

Collapse all files

trunk/src/src/FSfarras.py created 10644 +65 0

@@ -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

trunk/src/src/deb4_basis.py created 10644 +38 0

@@ -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

trunk/src/src/dualfilt1.py created 10644 +65 0

@@ -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

trunk/src/src/irls_dn.py created 10644 +80 0

@@ -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 + \|\|Au-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 = WA'(AWA'+ lambdaI)^-1 b
	43	# where W = diag(1/w)
	44	# where w = (u.^2 + eps).^(p/2-1)
	45
	46	# Update
	47	w = (u2 + 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 (AWAT + lambda1I), 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

trunk/src/src/irls_dn2.py created 10644 +50 0

@@ -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 + \|\|Au-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

trunk/src/src/RICS_paper_compare_noise.py +159 -121

             #!/usr/bin/env python
             #----------------------------------------------------------
             # Original MATLAB code developed by Brian Harding
             # Rewritten in python by Yolian Amaro
             # Python version 2.7
             # May 15, 2014
             # Jicamarca Radio Observatory
             #----------------------------------------------------------
             import math
-            #import cmath
-            #import scipy
-            #import matplotlib
             import numpy as np
-            #from numpy import linalg
             import matplotlib.pyplot as plt
             from scipy import linalg
             import time
             from y_hysell96 import*
             from deb4_basis import *
             from modelf import *
-            from scipy.optimize import fsolve
+            #from scipy.optimize import fsolve
             from scipy.optimize import root
             import pywt
+            from irls_dn2 import *
             ## 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
             dec = -1.24
             # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
             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]] )
             plt.figure(1)
             plt.plot(rx, ry, 'ro')
             plt.draw()
             # Jicamarca is nominally at a 45 degree angle
             theta = 45  - dec;
             # 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]] );
             # Rotate antenna positions to magnetic coord.
             AR = np.dot(R.T, antpos);
             # Only take the East component
             r = AR[0,:]
             r.sort() # ROW VECTOR?
             # 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
             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.
             # 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
             # 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
             # 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));
             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]);
                 for m in range(0, tempr.size):
                     factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
             fact = fact + b;
             factr = factr + b;
             # # 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):
             #         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
             #     else:
             #         fact[k] = 1
             #
             #
             # # model for f: triangle pulse
             # mu = -1.0/180*math.pi;
             # sig = 5.0/180*math.pi;
             # wind1 = theta > mu-sig and theta < mu;
             # wind2 = theta < mu+sig and theta > mu;
             # fact = wind1 * (theta - (mu - sig));
             # 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
             I = sum(fact)[0];
             fact = fact/I; # normalize to sum(f)==1
             factr = factr/I; # normalize to sum(f)==1
-            #plot(thetat,fact,'r'); hold on;
+            #plt.figure()
-            #plot(thetar,factr,'k.'); hold off;
+            #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):
             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
                 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):
                         for i2 in range(0,r.size):
                             R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat))))
                             R[i1,i2] = sum(R[i1,i2])
                     # Add uncertainty
                     # This is an ad-hoc way of adding "noise".  It models some combination of
                     # receiver noise and finite integration times.  We could use a more
                     # advanced model (like in Yu et al 2000) in the future.
                     # 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
                     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)
                     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
                     plt.figure(3);
                     plt.pcolor(abs(Rnz));
                     plt.colorbar();
                     # 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)));
                         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_capon = np.zeros(shape=(Nr,1));
-                    #tic_capon = time.time();
+                    tic_capon = time.time();
                     for i in range(0, thetar.size):
                         th = thetar[i];
                         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()
+                    toc_capon = time.time()
-                    # elapsed_time_capon = toc_capon - tic_capon;
+                    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)
                     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
                     # 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
                     # 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
                     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
                             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
                             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
+                    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)
+                    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();
                     #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
                     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));
                     ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
                     es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
                     chi2 = np.sum((es/sigma)**2);
                     # CS inversion using irls ########################
                     # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
-                    #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
+                    Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
-                    wavelet1 = pywt.Wavelet('db4')
+                    # REMOVE THIS?--------------------------------
-                    Phi, Psi, x = wavelet1.wavefun(level=3)
+                    #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();
-                    #s = irls_dn2(H2*Psi,g2,0.5,G);
+                    print 'H2 shape', H2.shape
+                    print 'Psi shape', Psi.shape
+                    s = irls_dn2(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
+                    # # # Scaling and shifting
-                # # # Only necessary for capon solution
+                    # # # Only necessary for capon solution
-                f_capon = f_capon/np.max(f_capon)*np.max(fact);
-                ### 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));
-                relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
-                relrmse_capon = rmse_capon / np.linalg.norm(fact);
+                    f_capon = f_capon/np.max(f_capon)*np.max(fact);
-                relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
-                #relrmse_cs = rmse_cs / np.norm(fact);
+                    ### analyze stuff ######################
-                factr = factr.T.conj()
+                    # calculate MSE
+                    rmse_fourier = np.sqrt(np.mean((f_fourier - factr)**2));
-                # calculate correlation
+                    rmse_capon = np.sqrt(np.mean((f_capon - factr)**2));
-                corr_fourier = np.dot(f_fourier,factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
+                    rmse_maxent = np.sqrt(np.mean((f_maxent - factr)**2));
-                corr_capon = np.dot(f_capon,factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
+                    #rmse_cs = np.sqrt(np.mean((f_cs - factr).^2));
-                corr_maxent = np.dot(f_maxent,factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
-                #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
+                    relrmse_fourier = rmse_fourier / np.linalg.norm(fact);
-                # calculate centered correlation
+                    relrmse_capon = rmse_capon / np.linalg.norm(fact);
-                f0 = factr - np.mean(factr);
+                    relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
-                f1 = f_fourier - np.mean(f_fourier);
+                    #relrmse_cs = rmse_cs / np.norm(fact);
-                corrc_fourier = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
-                f1 = f_capon - np.mean(f_capon);
+                    # To be able to perform dot product (align matrices) done below within the dot calculations
-                corrc_capon = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
-                f1 = f_maxent - np.mean(f_maxent);
-                corrc_maxent = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
+                    #f_fourier = f_fourier.T.conj()
-                #f1 = f_cs - mean(f_cs);
+                    #f_capon = f_capon.T.conj()
-                #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
+                    #f_maxent = f_maxent.T.conj()
+                    #factr = factr.T.conj()
-                # # # plot stuff #########################
+                    # calculate correlation
-                #---- Capon----
+                    corr_fourier = np.dot(f_fourier.T.conj(),factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
-                plt.figure(4)
+                    corr_capon = np.dot(f_capon.T.conj(),factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
-                plt.subplot(2, 1, 1)
+                    corr_maxent = np.dot(f_maxent.T.conj(),factr) / (np.linalg.norm(f_maxent)*np.linalg.norm(factr));
-                plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
+                    #corr_cs = np.dot(f_cs,factr) / (norm(f_cs)*norm(factr));
-                plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
-                plt.ylabel('Power (arbitrary units)')
-                plt.legend(loc='upper right')
+                    # calculate centered correlation
+                    f0 = factr - np.mean(factr);
-                # formatting y-axis
+                    f1 = f_fourier - np.mean(f_fourier);
-                locs,labels = plt.yticks()
-                plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
+                    corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
-                plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
+                    f1 = f_capon - np.mean(f_capon);
+                    corrc_capon = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
+                    f1 = f_maxent - np.mean(f_maxent);
-                #---- MaxEnt----
+                    corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
-                plt.subplot(2, 1, 2)
+                    #f1 = f_cs - mean(f_cs);
-                plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
+                    #corrc_cs = dot(f0,f1) / (norm(f0)*norm(f1));
-                plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
-                plt.ylabel('Power (arbitrary units)')
-                plt.legend(loc='upper right')
+                    # # # plot stuff #########################
-                # formatting y-axis
-                locs,labels = plt.yticks()
+                    #---- Capon----
-                plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
+                    plt.figure(4)
-                plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
+                    plt.subplot(2, 1, 1)
+                    plt.plot(180/math.pi*thetar, f_capon, 'r', label='Capon')
-                plt.show()
+                    plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
+                    plt.ylabel('Power (arbitrary units)')
+                    plt.legend(loc='upper right')
-            # #     subplot(3,1,2);
-            # #         plot(180/pi*thetar,f_maxent,'r-');
+                    # formatting y-axis
-            # #         hold on;
+                    locs,labels = plt.yticks()
-            # #         plot(180/pi*thetat,fact,'k--');
+                    plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
-            # #         hold off;
+                    plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
-            # #         ylim([min(f_cs) 1.1*max(fact)]);
-            # # #         title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_maxent, corr_maxent, corrc_maxent));
-            # #         ylabel({'Power';'(arbitrary units)'})
+                    #---- MaxEnt----
-            # # #         title 'Maximum Entropy Method'
+                    plt.subplot(2, 1, 2)
-            # #         legend('MaxEnt','Truth');
+                    plt.plot(180/math.pi*thetar, f_maxent, 'r', label='MaxEnt')
-            # #
+                    plt.plot(180/math.pi*thetat,fact, 'k--', label='Truth')
-            # #     subplot(3,1,3);
+                    plt.ylabel('Power (arbitrary units)')
-            # #         plot(180/pi*thetar,f_cs,'r-');
+                    plt.legend(loc='upper right')
-            # #         hold on;
-            # #         plot(180/pi*thetat,fact,'k--');
+                    # formatting y-axis
-            # #         hold off;
+                    locs,labels = plt.yticks()
-            # #         ylim([min(f_cs) 1.1*max(fact)]);
+                    plt.yticks(locs, map(lambda x: "%.1f" % x, locs*1e4))
-            # # #         title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
+                    plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
-            # # #         title 'Compressed Sensing - Debauchies Wavelets'
-            # #         xlabel 'Degrees'
+                    plt.show()
-            # #         ylabel({'Power';'(arbitrary units)'})
-            # #         legend('Comp. Sens.','Truth');
-            # #
+                    # #     PLOT PARA COMPRESSED SENSING
-            # # #     set(gcf,'Position',[749   143   528   881]); # CSL
+                    # #
-            # # #     set(gcf,'Position',[885   -21   528   673]); # macbook
+                    # #     subplot(3,1,3);
-            # #     pause(0.01);            
  No newline at end of file
+                    # #         plot(180/pi*thetar,f_cs,'r-');
+                    # #         hold on;
+                    # #         plot(180/pi*thetat,fact,'k--');
+                    # #         hold off;
+                    # #         ylim([min(f_cs) 1.1*max(fact)]);
+                    # # #         title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
+                    # # #         title 'Compressed Sensing - Debauchies Wavelets'
+                    # #         xlabel 'Degrees'
+                    # #         ylabel({'Power';'(arbitrary units)'})
+                    # #         legend('Comp. Sens.','Truth');
+                    # #
+                    # # #     set(gcf,'Position',[749   143   528   881]); # CSL
+                    # # #     set(gcf,'Position',[885   -21   528   673]); # macbook
+                    # #     pause(0.01);
+                    # # Store Results
+                    corr[0, snri, Ni] = corr_fourier;
+                    corr[1, snri, Ni] = corr_capon;
+                    corr[2, snri, Ni] = corr_maxent;
+                    #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;
+                    corrc[0,snri,Ni] = corrc_fourier;
+                    corrc[1,snri,Ni] = corrc_capon;
+                    corrc[2,snri,Ni] = corrc_maxent;
+                    #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 (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
+                    print corr
+  No newline at end of file

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages