jro_soft/compressed_sensing Commit - r23:24

Cleaned code, fixed bugs/errors, updated plotting part....

yamaro -

r23:24

parent child

trunk/src/src/FSfarras.py +1 -1

             '''
             Created on May 26, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
             #import pywt
  No newline at end of file
-            import pywt
  No newline at end of file
             import numpy as np
   No newline at end of file
             def FSfarras():
                 #function [af, sf] = FSfarras
   No newline at end of file
                 # Farras filters organized for the dual-tree
                 # complex DWT.
                 #
                 # USAGE:
                 #    [af, sf] = FSfarras
                 # OUTPUT:
                 #    af{i}, i = 1,2 - analysis filters for tree i
                 #    sf{i}, i = 1,2 - synthesis filters for tree i
                 # See farras, dualtree, dualfilt1.
                 #
                 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
                 # http://taco.poly.edu/WaveletSoftware/
                 #
                 # Translated to Python by Yolian Amaro
   No newline at end of file
   No newline at end of file
                 a1 = np.array( [
                       [                0,                  0],
                       [-0.08838834764832,   -0.01122679215254],
                       [ 0.08838834764832,   0.01122679215254],
                       [ 0.69587998903400,    0.08838834764832],
                       [ 0.69587998903400,   0.08838834764832],
                       [ 0.08838834764832,   -0.69587998903400],
                       [-0.08838834764832,   0.69587998903400],
                       [ 0.01122679215254,   -0.08838834764832],
                       [ 0.01122679215254,   -0.08838834764832],
                       [0,                                   0]
                      ] );
   No newline at end of file
                 a2 = np.array([
                       [ 0.01122679215254,                  0],
                       [ 0.01122679215254,                  0],
                       [-0.08838834764832,  -0.08838834764832],
                       [ 0.08838834764832,  -0.08838834764832],
                       [ 0.69587998903400,   0.69587998903400],
                       [ 0.69587998903400,  -0.69587998903400],
                       [ 0.08838834764832,   0.08838834764832],
                       [-0.08838834764832,   0.08838834764832],
                       [                0,   0.01122679215254],
                       [                0,  -0.01122679215254]
                     ]);
   No newline at end of file
   No newline at end of file
                 af = np.array([  [a1,a2]  ], dtype=object)
   No newline at end of file
                 s1 = a1[::-1]
                 s2 = a2[::-1]
   No newline at end of file
                 sf = np.array([  [s1,s2]  ], dtype=object)
   No newline at end of file
                 return af, sf

trunk/src/src/RICS_paper_compare_noise.py +38 -69

             #!/usr/bin/env python
   No newline at end of file
             #----------------------------------------------------------
             # Original MATLAB code developed by Brian Harding
             # Rewritten in Python by Yolian Amaro
             # Python version 2.7
             # May 15, 2014
             # Jicamarca Radio Observatory
             #----------------------------------------------------------
   No newline at end of file
             import time
             import matplotlib.pyplot as plt
             from scipy.optimize import root
             from scipy.stats import nanmean
              No newline at end of file
              No newline at end of file
-            from y_hysell96 import*
  No newline at end of file
+             No newline at end of file
+            from y_hysell96 import *
  No newline at end of file
             from deb4_basis import *
             from modelf import *
             from irls_dn2 import *
              No newline at end of file
-            #from scipy.optimize import fsolve
  No newline at end of file
   No newline at end of file
   No newline at end of file
             #-------------------------------------------------------------------------------------------------
             # Set parameters
             #-------------------------------------------------------------------------------------------------
   No newline at end of file
             ## Calculate Forward Model
             lambda1 = 6.0
             k = 2*np.pi/lambda1
   No newline at end of file
             ## Magnetic Declination
             dec = -1.24
   No newline at end of file
             ## Loads Jicamarca antenna positions
             antpos = np.loadtxt("antpos.txt", comments="#", delimiter=";", unpack=False)
   No newline at end of file
             ## 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]] )
   No newline at end of file
             ## Plot of antenna positions
             plt.figure(1)
             plt.plot(rx, ry, 'ro')
             plt.draw()
   No newline at end of file
             ## Jicamarca is nominally at a 45 degree angle
             theta = 45  - dec;
   No newline at end of file
             ## 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) )
   No newline at end of file
             # Rotation matrix from antenna coord to magnetic coord (East North)
             R = np.array( [[val1, val3], [val2, val4]] )
   No newline at end of file
             # Rotate antenna positions to magnetic coord.
             AR = np.dot(R.T, antpos)
   No newline at end of file
             # Only take the East component
             r = AR[0,:]
             r.sort()
   No newline at end of file
             # Truth model (high and low resolution)
             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 = thetar.T                           # transpose
   No newline at end of file
   No newline at end of file
             #-------------------------------------------------------------------------------------------------
             # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and background constant b.
             #-------------------------------------------------------------------------------------------------
   No newline at end of file
             # Triple Gaussian
             # a = np.array([3, 5, 2]);
             # 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
   No newline at end of file
             # Double Gaussian
             # a = np.array([3, 5]);
             # 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
   No newline at end of file
             # Single Gaussian
             a = np.array( [3] )
             mu = np.array( [-3.0/180*np.pi] )
             sig = np.array( [2.0/180*np.pi] )
             b = 0
   No newline at end of file
             # Empty matrices for factors
             fact = np.zeros(shape=(Nt,1))
             factr = np.zeros(shape=(Nr,1))
   No newline at end of file
             # 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]*np.exp(temp[j]);
                 for m in range(0, tempr.size):
                     factr[m] = factr[m] + a[i]*np.exp(tempr[m]);
   No newline at end of file
             fact = fact + b;
             factr = factr + b;
   No newline at end of file
             # #-------------------------------------------------------------------------------------------------
             # # Model for f: Square pulse
             # #-------------------------------------------------------------------------------------------------
             # for j in range(0, fact.size):
             #     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*np.pi and thetar[k] < 2/180*np.pi):
             #         factr[k] = 0
             #     else:
             #         factr[k] = 1
   No newline at end of file
             # #-------------------------------------------------------------------------------------------------
             # # 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));
             # factr = wind1 * (thetar - (mu - sig));
             # fact = fact + wind2 * (-(theta-(mu+sig)));
             # factr = factr + wind2 * (-(thetar-(mu+sig)));
   No newline at end of file
             # 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
  No newline at end of file
   No newline at end of file
             I = sum(fact)[0];
             fact = fact/I;                              # normalize to sum(f)==1
             factr = factr/I;                            # normalize to sum(f)==1
   No newline at end of file
             # Plot Gaussian pulse(s)
             plt.figure(2)
             plt.plot(thetat, fact, 'r--')
             plt.plot(thetar, factr, 'ro')
             plt.draw()
   No newline at end of file
   No newline at end of file
             #-------------------------------------------------------------------------------------------------
             # Control the type and number of inversions with:
             # SNRdBvec: the SNRs that will be used.
             # NN: the number of trials for each SNR
             #-------------------------------------------------------------------------------------------------
   No newline at end of file
             #SNRdBvec = np.linspace(5,20,10);
             SNRdBvec = np.array([15]);                  # 15 dB
             NN = 1;                                     # number of trials at each SNR
   No newline at end of file
             # 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)
   No newline at end of file
             # For each SNR and trial
             for snri in range(0, SNRdBvec.size):
                 SNRdB = SNRdBvec[snri];
                 SNR = 10**(SNRdB/10.0);
   No newline at end of file
                 for Ni in range(0, NN):
                     # Calculate cross-correlation matrix (Fourier components of image)
                     # This is an inefficient way to do this.
   No newline at end of file
                     R = np.zeros(shape=(r.size, r.size), dtype=object);
   No newline at end of file
                     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])
   No newline at end of file
                     # 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.
   No newline at end of file
                     # This is a way of adding noise while maintaining the
                     # positive-semi-definiteness of the matrix.
   No newline at end of file
                     U = linalg.cholesky(R.astype(complex), lower=False);    # U'*U = R
   No newline at end of file
                     sigma_noise = (np.linalg.norm(U,'fro')/SNR);
   No newline at end of file
             #         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)
   No newline at end of file
                     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));
   No newline at end of file
                     Unz = U + nz;
                     Rnz = np.dot(Unz.T.conj(),Unz);         # the noisy version of R
                     plt.figure(3);
                     plt.pcolor(abs(Rnz));
                     plt.colorbar();
   No newline at end of file
                     #-------------------------------------------------------------------------------------------------
                     # Fourier Inversion
                     #-------------------------------------------------------------------------------------------------
                     f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
   No newline at end of file
                     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);
   No newline at end of file
                     f_fourier = f_fourier.real;             # get rid of numerical imaginary noise
   No newline at end of file
   No newline at end of file
                     #-------------------------------------------------------------------------------------------------
                     # Capon Inversion
                     #-------------------------------------------------------------------------------------------------
                     f_capon = np.zeros(shape=(Nr,1));
   No newline at end of file
                     tic_capon = time.time();
   No newline at end of file
                     for i in range(0, thetar.size):
                         th = thetar[i];
                         w = np.exp(1j*k*np.dot(r,np.sin(th)));
                         f_capon[i] = 1/ ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real
  No newline at end of file
-                        f_capon[i] = np.divide(1, ( np.dot( w.T.conj(), (linalg.solve(Rnz,w)) ) ).real)
  No newline at end of file
   No newline at end of file
                     toc_capon = time.time()
                     elapsed_time_capon = toc_capon - tic_capon;
                     elapsed_time_capon = toc_capon - tic_capon;
  No newline at end of file
                     f_capon = f_capon.real;                 # get rid of numerical imaginary noise
                     f_capon = f_capon.real;                 # get rid of numerical imaginary noise
  No newline at end of file
   No newline at end of file
                     #-------------------------------------------------------------------------------------------------
                     # MaxEnt Inversion
                     #-------------------------------------------------------------------------------------------------
                     #------------------------------------------------------------------------------------------------- 
  No newline at end of file
                     # 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
                     Hr = np.zeros(shape=(M,Nr));            # approximate sensing matrix for reconstruction
  No newline at end of file
                     # 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
  No newline at end of file
                     # 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
                     [i1v,i2v] = np.where(np.triu(T,1) > 0); # converts linear to triangular indexing
  No newline at end of file
                     # 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;
                             idx2 = 2*idx+1;
                             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;
                             gnz[idx2] = (Rnz[i1,i2]).imag*Nr/Nt;
  No newline at end of file
                     # Inversion
                     F = Nr/Nt;                              # normalization
                     sigma = 1;                              # set to 1 because the difference is accounted for in G
                     sigma = 1;                              # set to 1 because the difference is accounted for in G
  No newline at end of file
                     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)
                     lambda0 = 1e-5*np.ones(shape=(M,1));    # initial condition (can be set to anything)
  No newline at end of file
   No newline at end of file
                     # Whitened solution
                     def myfun(lambda1):
                         return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
                         return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
  No newline at end of file
   No newline at end of file
                     tic_maxEnt = time.time();               # start time maxEnt
                     tic_maxEnt = time.time();               # start time maxEnt
  No newline at end of file
                     lambda1 = root(myfun,lambda0, method='krylov',  tol=1e-14);
                     lambda1 = root(myfun,lambda0, method='krylov',  tol=1e-14);
  No newline at end of file
                     toc_maxEnt = time.time()
                     elapsed_time_maxent = toc_maxEnt - tic_maxEnt;
                     elapsed_time_maxent = toc_maxEnt - tic_maxEnt;
  No newline at end of file
                     # Solution
                     lambda1 = lambda1.x;
                     lambda1 = lambda1.x;
  No newline at end of file
                     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;
                     chi2 = np.sum((es/sigma)**2);
  No newline at end of file
-                    es = np.dot(Hr, f_maxent) - gnz;        # should be same as ep
  No newline at end of file
                     chi2 = np.sum((es/sigma)**2);
  No newline at end of file
   No newline at end of file
                     #-------------------------------------------------------------------------------------------------
                     # CS inversion using Iteratively Reweighted Least Squares (IRLS)
                     #-------------------------------------------------------------------------------------------------
                     # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
                     # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
  No newline at end of file
                     Psi = deb4_basis(Nr);
                     Psi = deb4_basis(Nr);
              No newline at end of file
              No newline at end of file
-                    # ------------Remove this-------------------------------------------
              No newline at end of file
-                    # wavelet1 = pywt.Wavelet('db4')
              No newline at end of file
-                    # Phi, Psi, x = wavelet1.wavefun(level=3)
              No newline at end of file
-                    #-------------------------------------------------------------------
  No newline at end of file
                     # add "sum to 1" constraint
                     H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
                     g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
                     g2 = np.concatenate( (gnz, np.array([[Nr/Nt]])), axis=0 );
  No newline at end of file
                     tic_cs = time.time();
                     tic_cs = time.time();
  No newline at end of file
                     # Inversion
              No newline at end of file
              No newline at end of file
-            #         plt.figure(4)
              No newline at end of file
-            #         plt.imshow(Psi)#, interpolation='nearest')
              No newline at end of file
-            #         #plt.xticks([]); plt.yticks([])
              No newline at end of file
-            #         plt.show()
              No newline at end of file
-  No newline at end of file
                     s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
                     s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
  No newline at end of file
                     toc_cs = time.time()
                     elapsed_time_cs = toc_cs - tic_cs;
-                    #print s
  No newline at end of file
+             No newline at end of file
+  No newline at end of file
                     # Brightness function
                     f_cs = np.dot(Psi,s);
                     f_cs = np.dot(Psi,s);
  No newline at end of file
   No newline at end of file
-                    toc_cs = time.time()
              No newline at end of file
-                    elapsed_time_cs = toc_cs - tic_cs;
  No newline at end of file
                     # Plot
                     plt.figure(4)
                     plt.plot(thetar,f_cs,'r.-');
                     plt.plot(thetat,fact,'k-');
                     plt.plot(thetat,fact,'k-');   
  No newline at end of file
   No newline at end of file
                     #-------------------------------------------------------------------------------------------------
                     # Scaling and shifting
                     # (Only necessary for Capon solution)
                     #-------------------------------------------------------------------------------------------------
                     f_capon = f_capon/np.max(f_capon)*np.max(fact);
                     f_capon = f_capon/np.max(f_capon)*np.max(fact);
  No newline at end of file
   No newline at end of file
                     #-------------------------------------------------------------------------------------------------
                     # Analyze stuff
                     #-------------------------------------------------------------------------------------------------
                     #-------------------------------------------------------------------------------------------------
  No newline at end of file
                     # 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));
  No newline at end of file
   No newline at end of file
                     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.linalg.norm(fact);
                     relrmse_cs = rmse_cs / np.linalg.norm(fact);
  No newline at end of file
   No newline at end of file
                     # 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.T.conj(),factr) / (np.linalg.norm(f_cs)*np.linalg.norm(factr));
                     corr_cs = np.dot(f_cs.T.conj(),factr) / (np.linalg.norm(f_cs)*np.linalg.norm(factr));
  No newline at end of file
   No newline at end of file
                     # Calculate centered correlation
                     f0 = factr - np.mean(factr);
                     f1 = f_fourier - np.mean(f_fourier);
                     f1 = f_fourier - np.mean(f_fourier);
  No newline at end of file
                     corrc_fourier = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
                     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);
                     corrc_maxent = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
                     f1 = f_cs - np.mean(f_cs);
                     corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
                     corrc_cs = np.dot(f0.T.conj(),f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));  
  No newline at end of file
   No newline at end of file
                     #-------------------------------------------------------------------------------------------------
                     # Plot stuff
                     #-------------------------------------------------------------------------------------------------
                     #-------------------------------------------------------------------------------------------------
  No newline at end of file
                     #---- Capon----#
                     plt.figure(5)
  No newline at end of file
-                    #---- Capon----
  No newline at end of file
                     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')
                     plt.legend(loc='upper right')
  No newline at end of file
                     # 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)
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
  No newline at end of file
   No newline at end of file
                     #---- MaxEnt----#
                     plt.subplot(3, 1, 2)
  No newline at end of file
-                    #---- MaxEnt----
  No newline at end of file
                     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')
                     plt.legend(loc='upper right')
  No newline at end of file
                     # 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)
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
  No newline at end of file
   No newline at end of file
                     #---- Compressed Sensing----#
                     plt.subplot(3, 1, 3)
  No newline at end of file
-                    #---- Compressed Sensing----
  No newline at end of file
                     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')
                     plt.legend(loc='upper right')
  No newline at end of file
                     # 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)
                     plt.text(0.0, 1.01, '1e-4', fontsize=10, transform = plt.gca().transAxes)
  No newline at end of file
                     # # Store Results
              No newline at end of file
              No newline at end of file
-                    # #     PLOT PARA COMPRESSED SENSING
              No newline at end of file
-                    # #
              No newline at end of file
-                    # #     subplot(3,1,3);
              No newline at end of file
-                    # #         plot(180/pi*thetar,f_cs,'r-');
              No newline at end of file
-                    # #         hold on;
              No newline at end of file
-                    # #         plot(180/pi*thetat,fact,'k--');
              No newline at end of file
-                    # #         hold off;
              No newline at end of file
-                    # #         ylim([min(f_cs) 1.1*max(fact)]);
              No newline at end of file
-                    # # #         title(sprintf('rel. RMSE: #.2e\tCorr: #.3f Corrc: #.3f', relrmse_cs, corr_cs, corrc_cs));
              No newline at end of file
-                    # # #         title 'Compressed Sensing - Debauchies Wavelets'
              No newline at end of file
-                    # #         xlabel 'Degrees'
              No newline at end of file
-                    # #         ylabel({'Power';'(arbitrary units)'})
              No newline at end of file
-                    # #         legend('Comp. Sens.','Truth');
              No newline at end of file
-                    # #
              No newline at end of file
-                    # # #     set(gcf,'Position',[749   143   528   881]); # CSL
              No newline at end of file
-                    # # #     set(gcf,'Position',[885   -21   528   673]); # macbook
              No newline at end of file
-                    # #     pause(0.01);
              No newline at end of file
              No newline at end of file
-  No newline at end of file
                     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;
  No newline at end of file
                     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;     
  No newline at end of file
                     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;
  No newline at end of file
                     print "--------Time performace--------"
                     print 'Capon:\t', elapsed_time_capon, 'sec';
  No newline at end of file
-  No newline at end of file
                     print 'Maxent:\t',elapsed_time_maxent, 'sec';
                     print 'CS:\t',elapsed_time_cs, 'sec\n';
              No newline at end of file
-                    print 'CS:\t',elapsed_time_cs, 'sec';
+             No newline at end of file
                     print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN), '\n';
              No newline at end of file
              No newline at end of file
-                    print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
+             No newline at end of file
              No newline at end of file
              No newline at end of file
                     #-------------------------------------------------------------------------------------------------
-                    print corr
+             No newline at end of file
                     # Analyze and plot statistics
              No newline at end of file
                     #-------------------------------------------------------------------------------------------------
-                    print corr.shape
+             No newline at end of file
   No newline at end of file
              No newline at end of file
              No newline at end of file
-                    ## Analyze and plot statistics
  No newline at end of file
                     metric = corr; # set this to rmse, corr, or corrc
                     metric = corr; # set this to rmse, corr, or corrc
  No newline at end of file
                     # 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
                     # end
  No newline at end of file
   No newline at end of file
             # Avg ignoring NaNs
             ave = np.zeros(shape=(4))
-            def nanmean(data, **args):
+             No newline at end of file
              No newline at end of file
-                return numpy.ma.filled(numpy.ma.masked_array(data,numpy.isnan(data)).mean(**args), fill_value=numpy.nan)
+             No newline at end of file
             ave[0] = nanmean(metric[0,:,:]);    # Fourier
              No newline at end of file
             ave[1] = nanmean(metric[1,:,:]);    # Capon
-            # ave = np.zeros(shape=(4))
+             No newline at end of file
             ave[2] = nanmean(metric[2,:,:]);    # MaxEnt
              No newline at end of file
             ave[3] = nanmean(metric[3,:,:]);    # Compressed Sensing
-            # ave[0] = nanmean(metric, axis=0);
+             No newline at end of file
              No newline at end of file
-            # ave[1] = nanmean(metric, axis=1);
+             No newline at end of file
             # Plot based on chosen metric
  No newline at end of file
-            # ave[2] = nanmean(metric, axis=2);
              No newline at end of file
-            # ave[3] = nanmean(metric, axis=3);
              No newline at end of file
              No newline at end of file
-            #print ave
  No newline at end of file
             plt.figure(6);
             f = plt.scatter(SNRdBvec, ave[0], marker='+', color='b', s=60);     # Fourier
-            f = plt.scatter(SNRdBvec, corr[0], marker='+', color='b', s=60);   # Fourier
+             No newline at end of file
             c = plt.scatter(SNRdBvec, ave[1], marker='o', color= 'g', s=60);    # Capon
-            c = plt.scatter(SNRdBvec, corr[1], marker='o', color= 'c', s=60);     # Capon
+             No newline at end of file
             me= plt.scatter(SNRdBvec, ave[2], marker='s', color= 'c', s=60);    # MaxEnt
-            me= plt.scatter(SNRdBvec, corr[2], marker='s', color= 'y', s=60);    # MaxEnt
+             No newline at end of file
             cs= plt.scatter(SNRdBvec, ave[3], marker='*', color='r', s=60);     # Compressed Sensing
  No newline at end of file
-            cs= plt.scatter(SNRdBvec, corr[3], marker='*', color='r', s=60);        # Compressed Sensing
              No newline at end of file
-  No newline at end of file
   No newline at end of file
             plt.legend((f,c,me,cs),('Fourier','Capon', 'MaxEnt', 'Comp. Sens.'),scatterpoints=1, loc='upper right')
             plt.xlabel('SNR')
             plt.ylabel('Correlation with Truth')
   No newline at end of file
+            print "--------Correlations--------"
              No newline at end of file
+            print "Fourier:", ave[0]
              No newline at end of file
+            print "Capon:\t", ave[1]
              No newline at end of file
+            print "MaxEnt:\t", ave[2]
              No newline at end of file
+            print "CS:\t", ave[3]
              No newline at end of file
+  No newline at end of file
             plt.show()
   No newline at end of file

trunk/src/src/deb4_basis.py 0 -1

             '''
             Created on May 26, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
              No newline at end of file
-            import numpy as np
  No newline at end of file
             from FSfarras import *
             from dualfilt1 import *
             from dualtree import *
             from idualtree import *
   No newline at end of file
             # Debauchie 4 Wavelet
             def deb4_basis(N):
   No newline at end of file
                 Psi = np.zeros(shape=(N,2*N+1));
                 idx = 0;
                 J = 4;
                 [Faf, Fsf] = FSfarras();
                 [af, sf] = dualfilt1();
   No newline at end of file
                 # compute transform of zero vector
                 x = np.zeros(shape=(1,N));
                 w = dualtree(x, J, Faf, af);
   No newline at end of file
   No newline at end of file
                 # Uses both real and imaginary wavelets
                 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);
                             w[i][j][0,k] = 0;
                             # store it
                             Psi[:,idx] = y.T.conj();
                             idx = idx + 1;
   No newline at end of file
                 # Add uniform vector (seems to be useful if there's a background
                 Psi[:,2*N] = 1/np.sqrt(N);
   No newline at end of file
                 return Psi

trunk/src/src/irls_dn.py 0 -23

             '''
             Created on May 27, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
              No newline at end of file
-            #from scipy.sparse import eye
  No newline at end of file
             from scipy import linalg
             import scipy.sparse as sps
             import numpy as np
             from numpy.linalg import norm
   No newline at end of file
             def irls_dn(A,b,p,lambda1):
   No newline at end of file
   No newline at end of file
                 # Minimize lambda*||u||_p + ||A*u-b||_2, 0 < p <= 1
                 # using Iterative Reweighted Least Squares
                 # (see http://math.lanl.gov/Research/Publications/Docs/chartrand-2008-iteratively.pdf
                 # and http://web.eecs.umich.edu/~aey/sparse/sparse11.pdf)
   No newline at end of file
                 # Note to self:  I found that "warm-starting" didn't really help too much.
   No newline at end of file
                 [M,N] = A.shape;
                 # Initialize and precompute:
                 eps = 1e-2; # damping parameter
   No newline at end of file
                 [Q,R] = linalg.qr(A.T.conj(), mode='economic');
   No newline at end of file
   No newline at end of file
                 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);
   No newline at end of file
                 # Temporary N x N matrix
                 temp = np.zeros(shape=(N,N))
              No newline at end of file
              No newline at end of file
-                #---------- not needed, defined above--------------
              No newline at end of file
-                # Spacing of floating point numbers
              No newline at end of file
-                #eps = np.spacing(1)
              No newline at end of file
-                #--------------------------------------------------   
  No newline at end of file
   No newline at end of file
                 # Loop until damping parameter is small enough
                 while (eps > 1e-7):
                     epschange = 0;
                     # Loop until it's time to change eps
                     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)
   No newline at end of file
                         # Update
                         w = (u**2 + eps)**(1-p/2.0);
              No newline at end of file
              No newline at end of file
-            #             #---- Very inefficient- REMOVE THIS PART------
              No newline at end of file
-            #             k = 0
              No newline at end of file
-            #             # Sparse matrix
              No newline at end of file
-            #             for i in range (0, N):
              No newline at end of file
-            #                 for j in range (0,N):
              No newline at end of file
-            #                     if(i==j):
              No newline at end of file
-            #                         temp[i,j] = w[k]
              No newline at end of file
-            #                         k = k+1
              No newline at end of file
-                        #--------------------------------------------------
              No newline at end of file
-  No newline at end of file
                         np.fill_diagonal(temp, w)
              No newline at end of file
-                        #-----------------------------------------------
  No newline at end of file
   No newline at end of file
                         # Compressed Sparse Matrix
                         W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
   No newline at end of file
                         WAT = W*A.T.conj();
                         WAT = W*A.T.conj();
              No newline at end of file
              No newline at end of file
-                        #print "WAT", WAT.shape
              No newline at end of file
-                        #print "np.dot(A,WAT)", np.dot(A,WAT).shape
              No newline at end of file
-                        #print "np.dot(lambda1,I)", np.dot(lambda1,I).shape
              No newline at end of file
-                        #print "linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b)", linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b).shape
  No newline at end of file
                         u_new = np.dot(WAT , linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b));
                         u_new = np.dot(WAT , linalg.solve((np.dot(A,WAT) + np.dot(lambda1,I)), b));
  No newline at end of file
                         # See if this subproblem is converging
                         delu = norm(u_new-u)/norm(u);
                         epschange = delu < (np.sqrt(eps)/100.0);
                         epschange = delu < (np.sqrt(eps)/100.0);
  No newline at end of file
                         # Make update
                         u = u_new;
                         u = u_new;
  No newline at end of file
   No newline at end of file
                     eps = eps/10.0; # decrease eps
                     eps = eps/10.0; # decrease eps
  No newline at end of file
                 return u
                 return u
  No newline at end of file
   No newline at end of file
   No newline at end of file
   No newline at end of file

trunk/src/src/irls_dn2.py +6 -26

             '''
             Created on May 30, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
   No newline at end of file
             from irls_dn import *
             from scipy.optimize import brentq
  No newline at end of file
-            from scipy.optimize import fsolve
              No newline at end of file
-            import numpy as np
              No newline at end of file
-            from scipy.optimize import *
              No newline at end of file
-            from dogleg import *
              No newline at end of file
-            from numpy.linalg import norm
              No newline at end of file
              No newline at end of file
-            import matplotlib.pyplot as plt
  No newline at end of file
   No newline at end of file
             def irls_dn2(A,b,p,G):
   No newline at end of file
                 # Minimize ||u||_p subject to ||A*u-b||_2^2 <= G   (0 < p <= 1)
   No newline at end of file
                 # What this function actually does is finds the lambda1 so that the solution
                 # to the following problem satisfies ||A*u-b||_2^2 <= G:
                 # Minimize lambda1*||u||_p + ||A*u-b||_2
   No newline at end of file
                 # Start with a large lambda1, and do a line search until fidelity <= G.
                 # (Inversions with large lambda1 are really fast anyway).
   No newline at end of file
                 # Then spin up fsolve to localize the root even better
   No newline at end of file
                 # Line Search
   No newline at end of file
                 alpha = 2.0; # Line search parameter
                 lambda1 = 1e5; # What's a reasonable but safe initial guess?
                 u = irls_dn(A,b,p,lambda1);
   No newline at end of file
-                #print "u\n", u
  No newline at end of file
   No newline at end of file
                 fid = norm(np.dot(A,u)-b)**2;
   No newline at end of file
                 print '----------------------------------\n';
   No newline at end of file
                 while (fid >= G):
                     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;
  No newline at end of file
-                    print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
              No newline at end of file
-                    #print u
  No newline at end of file
   No newline at end of file
                 # Refinement using fzero/ brentq
              No newline at end of file
-                lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
  No newline at end of file
   No newline at end of file
+                lambda0 = np.array([lambda1,lambda1*alpha]);    # interval with zero-crossing
  No newline at end of file
   No newline at end of file
                 def myfun(lambda1):
                     return norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
  No newline at end of file
-                    temp1 = np.dot(A, irls_dn(A,b,p,lambda1))
              No newline at end of file
-                    temp2 = norm(temp1-b)
              No newline at end of file
-                    temp3 = temp2**2-G
              No newline at end of file
-                    #return np.linalg.norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
              No newline at end of file
-                    return temp3
  No newline at end of file
                 # Find zero-crossing at given interval (lambda1, lambda1*alpha)
  No newline at end of file
-                print "tolerancia=", 0.01*lambda1
              No newline at end of file
              No newline at end of file
-                #lambda1 = root(myfun, lambda1, method='krylov',  tol=0.01*lambda1);
              No newline at end of file
-                #lambda1 = lambda1.x
              No newline at end of file
              No newline at end of file
-                print "lambda0[0]", lambda0[0]
              No newline at end of file
-                print "lambda0[1]", lambda0[1]
              No newline at end of file
-  No newline at end of file
                 lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1)
              No newline at end of file
              No newline at end of file
-                print "lambda final=", lambda1
  No newline at end of file
   No newline at end of file
                 u = irls_dn(A,b,p,lambda1);
   No newline at end of file
                 return u;
                 return u;
  No newline at end of file
   No newline at end of file

trunk/src/src/sfb.py +1 -2

             '''
             Created on Jun 5, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
             from multirate import upfirdn
  No newline at end of file
-            from multirate import *
              No newline at end of file
-            import numpy as np
  No newline at end of file
             from cshift import *
   No newline at end of file
             def sfb(lo, hi, sf):
   No newline at end of file
                 # Synthesis filter bank
                 #
                 # USAGE:
                 #    y = sfb(lo, hi, sf)
                 # INPUT:
                 #    lo - low frqeuency input
                 #    hi - high frequency input
                 #    sf - synthesis filters
                 #    sf(:, 1) - lowpass filter (even length)
                 #    sf(:, 2) - highpass filter (even length)
                 # OUTPUT:
                 #    y - output signal
                 # See also afb
                 #
                 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
                 # http://taco.poly.edu/WaveletSoftware/
   No newline at end of file
                 N = 2*lo.size;
                 L = sf.size/2;
   No newline at end of file
                 # Need to change format for upfirdn funct:
                 lo = lo.T.conj()
                 lo = lo.reshape(lo.size)
   No newline at end of file
                 #print 'sfb hi', hi
   No newline at end of file
                 # Need to change format for upfirdn funct:
                 hi = hi.T.conj()
                 hi = hi.reshape(hi.size)
   No newline at end of file
                 #hi = hi.reshape(1, hi.size)
   No newline at end of file
   No newline at end of file
                 lo = upfirdn(lo, sf[:,0], 2, 1);
                 hi = upfirdn(hi, sf[:,1], 2, 1);
   No newline at end of file
                 lo = lo[0:lo.size-1]
                 hi = hi[0:hi.size-1]
   No newline at end of file
                 y = lo + hi;
                 y[0:L-2] = y[0:L-2] + y[N+ np.arange(0,L-2)]; #CHECK IF ARANGE IS CORRECT
                 y = y[0:N];
   No newline at end of file
                 #print 'y en sbf\n', y.shape
   No newline at end of file
                 y = y.reshape(1, y.size)
                 #print 'y en sbf\n', y.shape
   No newline at end of file
                 y = cshift(y, 1-L/2);
   No newline at end of file
                 return y;
   No newline at end of file

trunk/src/src/y_hysell96.py 0 -1

             '''
             Created on May 22, 2014
   No newline at end of file
             @author: Yolian Amaro
             '''
              No newline at end of file
-            import numpy as np
  No newline at end of file
             from modelf import *
   No newline at end of file
             def y_hysell96(lambda1,g,sigma,F,G,H):
                 # Y_HYSELL96 Implements set of nonlinear equations to solve Hysell96 MaxEnt
                 # y(lambda) = 0
                 # decision variables: lambda
                 # g: data
                 # sigma: uncertainties (length of g)
                 # F: sum(f)
                 # G: desired value for chi^2
                 # H: linear operator mapping image (f) to data (g)
                 # This function is a helper function that returns 0 when a value of lambda
                 # is chosen that satisfies the equations.
   No newline at end of file
                 # model for f
                 f = modelf(lambda1, H,F);
   No newline at end of file
                 # solve for Lambda and e
                 Lambda = np.sqrt(np.sum(np.multiply(lambda1**2,sigma**2))/(4*G)); # positive root (right?)
                 e = np.multiply(-lambda1,sigma**2) / (2*Lambda);
   No newline at end of file
                 # measurement equation
                 y = g + e - np.dot(H, f);
   No newline at end of file
                 return y

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