jro_soft/compressed_sensing Commit - r24:25

Update

yamaro -

r24:25

parent child

trunk/src/src/RICS_paper_compare_noise.py +55 -24

             ## Plot of antenna positions
             plt.figure(1)
             plt.plot(rx, ry, 'ro')
+            plt.title("Jicamarca Antenna Positions")
+            plt.xlabel("rx")
+            plt.ylabel("ry")
+            plt.grid(True)
+            plt.plot(rx, rx, 'b-')
             plt.draw()
             ## Jicamarca is nominally at a 45 degree angle
             theta = 45  - dec;
-            ## Rotation matrix from antenna coord to magnetic coord (East North)
+            ## 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) )
             # 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.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)
             Nt = (1024.0)*(16.0)                        # number of pixels in truth image: high resolution
             # b = 0;                                    # background
             # Single Gaussian
-            a = np.array( [3] )
+            a = np.array( [3] )                 # amplitude(s)
-            mu = np.array( [-3.0/180*np.pi] )
+            mu = np.array( [-3.0/180*np.pi] )   # center
-            sig = np.array( [2.0/180*np.pi] )
+            sig = np.array( [2.0/180*np.pi] )   # width
-            b = 0
+            b = 0                               # background constant
             # Empty matrices for factors
             fact = np.zeros(shape=(Nt,1))
             factr = np.zeros(shape=(Nr,1))
-            # DFT Kernels
+            # Constructing Gaussian model
             for i in range(0, a.size):
                 temp = (-(thetat-mu[i])**2/(sig[i]**2))
                 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
                     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]);
             fact = fact + b;
             factr = factr + b;
             # #-------------------------------------------------------------------------------------------------
             # # Model for f: Square pulse
             # 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 = np.divide(fact, (np.sum(fact)*2*thbound/Nt));  # normalize to integral(f)==1
             I = sum(fact)[0];
             plt.figure(2)
             plt.plot(thetat, fact, 'r--')
             plt.plot(thetar, factr, 'ro')
+            plt.title("Truth models")
+            plt.ylabel("Amplitude")
+            plt.legend(('High res','Low res'), loc='upper right')
             plt.draw()
             SNRdBvec = np.array([15]);                  # 15 dB
             NN = 1;                                     # number of trials at each SNR
-            # Statistics simulation (correlation, root mean square)
+            # Statistics simulation (correlation, root mean square error)
             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)
                     # 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 = (-2*np.random.rand(U.shape[0], U.shape[1]) + 2)
+                    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));
-            #         temp2 = 1j*(-2*np.random.rand(U.shape[0], U.shape[1]) + 2)
-            #         temp3 = ((abs(U) > 0).astype(float))                    # upper triangle of 1's
-            #         temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
-            #         nz = np.multiply(temp4,temp3)
-                    nz = np.multiply( sigma_noise * (np.random.randn(U.shape[0]) + 1j*np.random.randn(U.shape[0]))/np.sqrt(2) , (abs(U) > 0).astype(float));
                     Unz = U + nz;
                     Rnz = np.dot(Unz.T.conj(),Unz);         # the noisy version of R
                     plt.figure(3);
-                    plt.pcolor(abs(Rnz));
+                    plt.pcolor(abs(R));
                     plt.colorbar();
                     #-------------------------------------------------------------------------------------------------
                     # 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
+                            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();
                     lambda1 = lambda1.x;
                     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);
+                    #--------------------------------------
                     #-------------------------------------------------------------------------------------------------
                     f_cs = np.dot(Psi,s);
+                    print "shapes:"
+                    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')
                     #-------------------------------------------------------------------------------------------------
                     print 'Maxent:\t',elapsed_time_maxent, 'sec';
                     print 'CS:\t',elapsed_time_cs, 'sec\n';
-                    print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN), '\n';
+                    print (NN*(snri)+Ni+1), '/', (SNRdBvec.size*NN), '\n';
                     #-------------------------------------------------------------------------------------------------
             # Avg ignoring NaNs
             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

trunk/src/src/irls_dn.py +3 -1

                 # Initialize and precompute:
                 eps = 1e-2; # damping parameter
+                # Compute the qr factorization of a matrix.
+                # Factor the matrix A as qr, where Q is orthonormal and R is upper-triangular.
                 [Q,R] = linalg.qr(A.T.conj(), mode='economic');
                 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);
+                I = sps.eye(M);                     # M as sparse matrix
                 # Temporary N x N matrix
                 temp = np.zeros(shape=(N,N))

trunk/src/src/irls_dn2.py +3 -4

                 # Line Search
-                alpha = 2.0; # Line search parameter
+                alpha = 2.0;                                    # Line search parameter
-                lambda1 = 1e5; # What's a reasonable but safe initial guess?
+                lambda1 = 1e5;                                  # What's a reasonable but safe initial guess?
                 u = irls_dn(A,b,p,lambda1);
                 print '----------------------------------\n';
                 while (fid >= G):
-                    lambda1 = lambda1 / alpha; # Balance between speed and accuracy
+                    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;
                 # Refinement using fzero/ brentq
                 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