jro_soft/compressed_sensing Commit - r16:17

Plots for Capon and MaxEntropy inversion methods. Compressed Sensing initial/partial implementation.

yamaro -

r16:17

parent child

trunk/src/src/RICS_paper_compare_noise.py +126 -21

              #!/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 root
+             import pywt
              ## 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
+             # 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
+             # 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
+             #     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;
              #plot(thetar,factr,'k.'); hold off;
              #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)
                      #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 %%%%%%%%%%%%%%%%%%%
+                     # 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 %%%%%%%%%%%%%%%%%%%%%%
+                     # Capon Inversion ######################
                      f_capon = np.zeros(shape=(Nr,1));
                      #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()
                      # elapsed_time_capon = toc_capon - tic_capon;
                      f_capon = f_capon.real; # get rid of numerical imaginary noise
-                     # MaxEnt Inversion %%%%%%%%%%%%%%%%%%%%%
+                     # MaxEnt Inversion #####################
-                     # create the appropriate sensing matrix (split into real and imaginary % parts)
+                     # 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
                      lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
                      # 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 %%%%%%%%%%%%%%%%%%%%%%%%
+                     # CS inversion using irls ########################
                      # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
-                     Psi = deb4_basis(Nr);
+                     #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?)
+                     wavelet1 = pywt.Wavelet('db4')
+                     Phi, Psi, x = wavelet1.wavefun(level=3)
-                     #         # add "sum to 1" constraint
-                     #         H2 = [Hr; ones(1,Nr)];
-                     #         g2 = [gnz; Nr/Nt];
+                     #
-                     #         s = irls_dn2(H2*Psi,g2,0.5,G);
+                     # 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 );
+                     #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;                   
  No newline at end of file
+                     #         hold off;
+                 # # # Scaling and shifting
+                 # # # 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);
+                 relrmse_maxent = rmse_maxent / np.linalg.norm(fact);
+                 #relrmse_cs = rmse_cs / np.norm(fact);
+                 factr = factr.T.conj()
+                 # calculate correlation
+                 corr_fourier = np.dot(f_fourier,factr) / (np.linalg.norm(f_fourier)*np.linalg.norm(factr));
+                 corr_capon = np.dot(f_capon,factr) / (np.linalg.norm(f_capon)*np.linalg.norm(factr));
+                 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));
+                 # calculate centered correlation
+                 f0 = factr - np.mean(factr);
+                 f1 = f_fourier - np.mean(f_fourier);
+                 corrc_fourier = np.dot(f0,f1) / (np.linalg.norm(f0)*np.linalg.norm(f1));
+                 f1 = f_capon - np.mean(f_capon);
+                 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));
+                 #f1 = f_cs - mean(f_cs);
+                 #corrc_cs = dot(f0,f1) / (norm(f0)*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.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)
+                 #---- 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.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)
+                 plt.show()
+             # #     subplot(3,1,2);
+             # #         plot(180/pi*thetar,f_maxent,'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_maxent, corr_maxent, corrc_maxent));
+             # #         ylabel({'Power';'(arbitrary units)'})
+             # # #         title 'Maximum Entropy Method'
+             # #         legend('MaxEnt','Truth');
+             # #
+             # #     subplot(3,1,3);
+             # #         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);            
  No newline at end of file

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

              '''
              Created on May 22, 2014
              @author: Yolian Amaro
              '''
              import numpy as np
              def modelf(lambda1, H, F):
                  # THRESH = 700;
-                 exponent = H.T.conj()*lambda1;
+                 exponent = np.dot(H.T.conj(),lambda1);
                  # exponent(exponent>THRESH) = THRESH; % numerical stabilizer
                  # exponent(exponent<-THRESH) = -THRESH;
                  E = np.max(exponent);
                  fnum = np.exp(exponent-E);
                  f = F*fnum/np.sum(fnum);
                  return f
  No newline at end of file

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

              '''
              Created on May 22, 2014
              @author: Yolian Amaro
              '''
+             import numpy as np
              from modelf import *
              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.
                  # model for f
                  f = modelf(lambda1, H,F);
                  # 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);
                  # measurement equation
                  y = g + e - np.dot(H, f);
                  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