jro_soft/compressed_sensing Commit - r25:26

Implementation of voltages simulation for estimating correlations (some things need to be fixed)

yamaro -

r25:26

parent child

trunk/src/src/RICS_paper_compare_noise.py +93 -48

              import time
              import matplotlib.pyplot as plt
              from scipy.optimize import root
+             from scipy.optimize import newton_krylov
              from scipy.stats import nanmean
              from y_hysell96 import *
              ## Rotation matrix from antenna coord to magnetic coord (East North) -> 45 degrees
              theta_rad = np.radians(theta)        #  trig functions take radians as argument
              val1 = float( np.cos(theta_rad) )
-             val2 = float( np.sin(theta_rad) )
-             val3 = float( -1*np.sin(theta_rad))
+             val2 = float( -1*np.sin(theta_rad))
+             val3 = float( np.sin(theta_rad) )
              val4 = float( np.cos(theta_rad) )
              # Rotation matrix from antenna coord to magnetic coord (East North)
-             R = np.array( [[val1, val3], [val2, val4]] )
+             R = np.array( [[val1, val2], [val3, val4]] )
              # Rotate antenna positions to magnetic coord.
              AR = np.dot(R.T, antpos)
-             plt.plot(AR[0], 0*AR[1], 'yo')
-             plt.title("Jicamarca Antenna Positions")
+             plt.plot(AR[0], 0*AR[1], 'yo')              # antenna positions only considering magnetic East-West component (x, 0)
+             plt.title("Jicamarca Antenna Positions")
              plt.xlabel("rx")
              plt.ylabel("ry")
              plt.grid(True)
              plt.plot(rx, 0*rx, 'g-')
              plt.draw()
-             print antpos
-             print "\n"
-             print AR
              # Only take the East component
              r = AR[0,:]
              r.sort()
-             # ---NEW---
-             r1 = AR[1,:]    # longitudes / distances
-             d = 40          # separation of antennas in meters
-             theta_w = 15    # angle between antenna and point of observation in degrees
-             # arg = k*(r + d*np.cos(theta_w)) # this is the exponent argument in the corr integral
-             # --------
              # Truth model (high and low resolution)
              fact = np.zeros(shape=(Nt,1))
              factr = np.zeros(shape=(Nr,1))
-             # Constructing Gaussian model
+             # Constructing Gaussian model - brightness function (?)
              for i in range(0, a.size):
                  temp = (-(thetat-mu[i])**2/(sig[i]**2))
                  tempr = (-(thetar-mu[i])**2/(sig[i]**2))
              # Plot Gaussian pulse(s)
              plt.figure(2)
-             plt.plot(thetat, fact, 'r--')
+             plt.plot(thetat, fact, 'bo')
              plt.plot(thetar, factr, 'ro')
-             plt.title("Truth models")
+             plt.title("Image and Reconstruction Domain")
              plt.ylabel("Amplitude")
-             plt.legend(('High res','Low res'), loc='upper right')
+             plt.legend(('High res (img domain)','Low res (reconstruction domain)'), loc='upper right')
              plt.draw()
                  SNR = 10**(SNRdB/10.0);
                  for Ni in range(0, NN):
+                     # -------------------------------------- Voltages estimation----------------------------------------
+                     v = np.zeros(shape=(r.size,1), dtype=object);                           # voltages matrix
+                     R = np.zeros(shape=(r.size, r.size), dtype='complex128');               # corr matrix
+                     mu1, sigma1 = 0, 1                                                      # mean and standard deviation
+                     num = 1
+                     for t in range (0, 10):
+                         num = t+1
+                         print num
+                         n1 =  np.random.normal(mu1, sigma1, [Nt,1])         # random variable
+                         n2 =  np.random.normal(mu1, sigma1, [Nt,1])         # random variable
+                         rho = (n1 + n2*1j)*np.sqrt(fact/2)                  # rho.size is Nt = 16834
+                                                                             # fact is size (Nt, 1) -> column vector
+                         thetat = thetat.reshape(Nt,1)                       # thetat is now (Nt, 1)
+                         for i in range(0, r.size):
+                             v[i] = sum((rho)*np.exp(-1j*k*((r[i])*np.cos(thetat))));     # v is sum ( (Nt, 1)) = 1 value
+                         # Computing correlations
+                         for i1 in range(0, r.size):
+                             for i2 in range(0,r.size):
+                                 R[i1,i2] = R[i1,i2] + (v[i1][0]*v[i2][0].conjugate())/ (np.sqrt(np.abs(v[i1][0])**2 * np.abs(v[i2][0])**2));
+                     R = R / num
+                     print "R from voltages"
+                     print R
+                     plt.figure(3);
+                     plt.pcolor(abs(R));
+                     plt.colorbar();
+                     plt.title("Antenna Correlations")
+                     plt.draw()
+                     plt.show()
+                 #--------------------------------------------------------------------------------------------------
                      # Calculate cross-correlation matrix (Fourier components of image)
                      # This is an inefficient way to do this.
+                     # MATLAB CORR::: R(i1,i2) = fact.' * exp(j*k*(r(i1)-r(i2))*sin(thetat));
                      R = np.zeros(shape=(r.size, r.size), dtype=object);
+                     thetat = thetat.reshape(Nt,)
                      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])
+                             R[i1,i2] = sum(np.dot(fact.T, np.exp(1j*k*((r[i1]-r[i2])*np.sin(thetat)))))
+                     print "R w/o voltages\n", R
+             #         plt.figure(3);
+             #         plt.pcolor(abs(R));
+             #         plt.colorbar();
+             #         plt.title("Antenna Correlations")
+             #         plt.xlabel("Antenna")
+             #         plt.ylabel("Antenna")
+             #         plt.draw()
+             #         plt.show()
                      # 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
                      Unz = U + nz;
                      Rnz = np.dot(Unz.T.conj(),Unz);         # the noisy version of R
-                     plt.figure(3);
-                     plt.pcolor(abs(R));
-                     plt.colorbar();
+             #         plt.figure(3);
+             #         plt.pcolor(abs(R));
+             #         plt.colorbar();
                      #-------------------------------------------------------------------------------------------------
                      # Fourier Inversion
                      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);
+                         w = np.exp(1j*k*np.dot(r,np.sin(th)));
+                         f_fourier[i] = np.dot(np.dot(w.T.conj(),U) , w);
                      f_fourier = f_fourier.real;             # get rid of numerical imaginary noise
                      # 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 1->27
+                             idx = np.where(np.logical_and((i1==i1v), (i2==i2v)))[0];    # kind of awkward (1->27)
                              idx1 = 2*idx;
                              idx2 = 2*idx+1;
                              Hr[idx1,:] = np.cos(k*(r[i1]-r[i2])*np.sin(thetar)).T.conj();
                      tic_maxEnt = time.time();               # start time maxEnt
-                     lambda1 = root(myfun,lambda0, method='krylov',  tol=1e-14);
+                     lambda1 = root(myfun,lambda0, method='krylov',  tol=1e-14); # krylov
+             #         lambda1 = newton_krylov(myfun,lambda0, method='lgmres',  x_tol=1e-14)
                      toc_maxEnt = time.time()
                      elapsed_time_maxent = toc_maxEnt - tic_maxEnt;
                      f_maxent = modelf(lambda1, Hr, F);
                      #------ what is this needed for?-------
-                     ystar = myfun(lambda1);
-                     Lambda = np.sqrt(sum(lambda1**2*sigma**2)/(4*G));
-                     ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
-                     es = np.dot(Hr, f_maxent) - gnz;
-                     chi2 = np.sum((es/sigma)**2);
+             #         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);
                      #--------------------------------------
                      print "Psi ", Psi.shape
                      print "s ", s.shape
                      print "f_CS ", f_cs.shape
-                     # Plot
-                     plt.figure(4)
-                     plt.plot(thetar,f_cs,'r.-');
-                     plt.plot(thetat,fact,'k-');
-                     plt.title('Reconstruction with Compressed Sensing')
+             #         # Plot
+             #         plt.figure(4)
+             #         plt.plot(thetar,f_cs,'r.-');
+             #         plt.plot(thetat,fact,'k-');
+             #         plt.title('Reconstruction with Compressed Sensing')
                      #-------------------------------------------------------------------------------------------------
                      # Calculate centered correlation
                      f0 = factr - np.mean(factr);
                      f1 = f_fourier - np.mean(f_fourier);
                      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));
              ave = np.zeros(shape=(4))
              print metric
-             ave[0] = nanmean(metric[0,:,:]);    # Fourier
-             ave[1] = nanmean(metric[1,:,:]);    # Capon
-             ave[2] = nanmean(metric[2,:,:]);    # MaxEnt
-             ave[3] = nanmean(metric[3,:,:]);    # Compressed Sensing
+             print metric.shape
+             ave[0] = nanmean(metric[0,0,:]);    # Fourier
+             ave[1] = nanmean(metric[1,0,:]);    # Capon
+             ave[2] = nanmean(metric[2,0,:]);    # MaxEnt
+             ave[3] = nanmean(metric[3,0,:]);    # Compressed Sensing
+             print ave
              # Plot based on chosen metric
              plt.figure(6);

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

                  # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
                  # http://taco.poly.edu/WaveletSoftware/
-                 # These cofficients are rounded to 8 decimal places.
+                 # These coefficients are rounded to 8 decimal places.
                  a1 = np.array([
                    [ 0.03516384000000,                  0],

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

                      fid = norm(np.dot(A,u)-b)**2;
                      print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid;
-                 # Refinement using fzero/ brentq
+                 # Refinement using brentq (equiv to fzero in matlab)
                  lambda0 = np.array([lambda1,lambda1*alpha]);    # interval with zero-crossing
                  def myfun(lambda1):

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