jro_soft/compressed_sensing Commit - r20:21

Fixed erros in sfb, idualtree, irls_dn, and irls_dn2

yamaro -

r20:21

parent child

trunk/src/src/RICS_paper_compare_noise.py +18 -18

             #!/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 numpy as np
             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
             from irls_dn2 import *
             ## Calculate Forward Model
             lambda1 = 6.0
             k = 2*math.pi/lambda1
             ## Calculate Magnetic Declination
             #     [~,~,dec] = igrf11magm(350e3, -11-56/60, -76-52/60, 2012);  check this
             # or calculate it with the above function
             dec = -1.24
             # loads rx, ry (Jicamarca antenna positions) #this can be done with numpy.loadtxt()
             rx = np.array( [[127.5000], [91.5000], [127.5000], [19.5000], [91.5000], [-127.5000], [-55.5000], [-220.8240]] )
             ry = np.array( [[127.5000], [91.5000], [91.5000], [55.5000], [-19.5000], [-127.5000], [-127.5000], [-322.2940]] )
             antpos = np.array( [[127.5000, 91.5000, 127.5000, 19.5000, 91.5000, -127.5000, -55.5000, -220.8240],
                             [127.5000, 91.5000, 91.5000, 55.5000, -19.5000, -127.5000, -127.5000, -322.2940]] )
             plt.figure(1)
             plt.plot(rx, ry, 'ro')
             plt.draw()
             # Jicamarca is nominally at a 45 degree angle
             theta = 45  - dec;
             # Rotation matrix from antenna coord to magnetic coord (East North)
             theta_rad = math.radians(theta)        #  trig functions take radians as argument
             val1 = float( math.cos(theta_rad) )
             val2 = float( math.sin(theta_rad) )
             val3 = float( -1*math.sin(theta_rad))
             val4 = float( math.cos(theta_rad) )
             # Rotation matrix from antenna coord to magnetic coord (East North)
             R = np.array( [[val1, val3], [val2, val4]] );
             # Rotate antenna positions to magnetic coord.
             AR = np.dot(R.T, antpos);
             # Only take the East component
             r = AR[0,:]
             r.sort() # ROW VECTOR?
             # Truth model (high and low resolution)
             Nt = (1024.0)*(16.0); # number of pixels in truth image: high resolution
             thbound = 9.0/180*math.pi; # the width of the domain in angle space
             thetat = np.linspace(-thbound, thbound,Nt) # image domain
             thetat = np.transpose(thetat) # transpose  # FUNCIONA??????????????????????????????
             Nr = (256.0); # number of pixels in reconstructed image: low res
             thetar = np.linspace(-thbound, thbound,Nr)  # reconstruction domain
             thetar = np.transpose(thetar) #transpose   # FUNCIONA?????????????????????????????
             # Model for f: Gaussian(s) with amplitudes a, centers mu, widths sig, and
             # background constant b.
             # Triple Gaussian
             # a = np.array([3, 5, 2]);
             # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi, 7.0/180*math.pi]);
             # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi, 0.3/180*math.pi]);
             # b = 0; # background
             # Double Gaussian
             # a = np.array([3, 5]);
             # mu = np.array([-5.0/180*math.pi, 2.0/180*math.pi]);
             # sig = np.array([2.0/180*math.pi, 1.5/180*math.pi]);
             # b = 0; # background
             # Single Gaussian
             a = np.array( [3] );
             mu = np.array( [-3.0/180*math.pi] )
             sig = np.array( [2.0/180*math.pi] )
             b = 0;
             fact = np.zeros(shape=(Nt,1));
             factr = np.zeros(shape=(Nr,1));
             for i in range(0, a.size):
                 temp = (-(thetat-mu[i])**2/(sig[i]**2))
                 tempr = (-(thetar-mu[i])**2/(sig[i]**2))
                 for j in range(0, temp.size):
                     fact[j] = fact[j] + a[i]*math.exp(temp[j]);
                 for m in range(0, tempr.size):
                     factr[m] = factr[m] + a[i]*math.exp(tempr[m]);
             fact = fact + b;
             factr = factr + b;
             # # model for f: Square pulse
             # for j in range(0, fact.size):
             #     if (theta > -5.0/180*math.pi and theta < 2.0/180*math.pi):
             #         fact[j] = 0
             #     else:
             #         fact[j] = 1
             # for k in range(0, factr.size):
             #     if (thetar[k] > -5.0/180*math.pi and thetar[k] < 2/180*math.pi):
             #         fact[k] = 0
             #     else:
             #         fact[k] = 1
             #
             #
             # # model for f: triangle pulse
             # mu = -1.0/180*math.pi;
             # sig = 5.0/180*math.pi;
             # wind1 = theta > mu-sig and theta < mu;
             # wind2 = theta < mu+sig and theta > mu;
             # fact = wind1 * (theta - (mu - sig));
             # factr = wind1 * (thetar - (mu - sig));
             # fact = fact + wind2 * (-(theta-(mu+sig)));
             # factr = factr + wind2 * (-(thetar-(mu+sig)));
             #     fact = fact/(sum(fact)[0]*2*thbound/Nt); # normalize to integral(f)==1
             I = sum(fact)[0];
             fact = fact/I; # normalize to sum(f)==1
             factr = factr/I; # normalize to sum(f)==1
             #plt.figure()
             #plt.plot(thetat,fact,'r');
             #plt.plot(thetar,factr,'k.');
             #xlim([min(thetat) max(thetat)]);
             #x = np.linspace(thetat.min(), thetat.max) ????
             #for i in range(0, thetat.size):
             plt.figure(2)
             plt.plot(thetat, fact, 'r--')
             plt.plot(thetar, factr, 'ro')
             plt.draw()
             # xlim([min(thetat) max(thetat)]); FALTA ARREGLAR ESTO
             ##
             # Control the type and number of inversions with:
             # SNRdBvec: the SNRs that will be used.
             # NN: the number of trials for each SNR
             #SNRdBvec = np.linspace(5,20,10);
             SNRdBvec = np.array([15]);
             NN = 1; # number of trial at each SNR
             # if using vector arguments should be: (4,SNRdBvec.size,NN)
             corr = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
             corrc = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
             rmse = np.zeros(shape=(4,SNRdBvec.size,NN)); # (method, SNR, trial)
             for snri in range(0, SNRdBvec.size): # change 1 for SNRdBvec.size when using SNRdBvec as vector
                 for Ni in range(0, NN):
                     SNRdB = SNRdBvec[snri];
                     SNR = 10**(SNRdB/10.0);
                     # Calculate cross-correlation matrix (Fourier components of image)
                     # This is an inefficient way to do this.
                     R = np.zeros(shape=(r.size, r.size), dtype=object);
                     for i1 in range(0, r.size):
                         for i2 in range(0,r.size):
                             R[i1,i2] = np.dot(fact.T, np.exp(1j*k*np.dot((r[i1]-r[i2]),np.sin(thetat))))
                             R[i1,i2] = sum(R[i1,i2])
                     # Add uncertainty
                     # This is an ad-hoc way of adding "noise".  It models some combination of
                     # receiver noise and finite integration times.  We could use a more
                     # advanced model (like in Yu et al 2000) in the future.
                     # This is a way of adding noise while maintaining the
                     # positive-semi-definiteness of the matrix.
                     U = linalg.cholesky(R.astype(complex), lower=False); # U'*U = R
                     sigma_noise = (np.linalg.norm(U,'fro')/SNR);
                     temp1 = (-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
                     temp2 = 1j*(-1*np.random.rand(U.shape[0], U.shape[1]) + 0.5)
                     temp3 = ((abs(U) > 0).astype(float))   # upper triangle of 1's
                     temp4 = (sigma_noise * (temp1 + temp2))/np.sqrt(2.0)
                     nz = np.multiply(temp4, temp3)
                     #---------------------- Eliminar esto:------------------------------------------
                     #nz = ((abs(np.multiply(temp4, temp3)) > 0).astype(int))
                     #nz = ((abs(np.dot(temp4, temp3)) > 0).astype(int))
                     #nz = np.dot(np.dot(sigma_noise, (temp1 + temp2)/math.sqrt(2), temp3 ));
                     #nz = np.dot(sigma_noise, (np.dot((np.random.rand(8,8) + j*np.random.rand(8,8))/math.sqrt(2.0) , (abs(U) > 0).astype(int))));
                     #--------------------------------------------------------------------------------
                     Unz = U + nz;
                     Rnz = np.dot(Unz.T.conj(),Unz); # the noisy version of R
                     plt.figure(3);
                     plt.pcolor(abs(Rnz));
                     plt.colorbar();
                     # Fourier Inversion ###################
                     f_fourier = np.zeros(shape=(Nr,1), dtype=complex);
                     for i in range(0, thetar.size):
                         th = thetar[i];
                         w = np.exp(1j*k*np.dot(r,np.sin(th)));
                         temp = np.dot(w.T.conj(),U)
                         f_fourier[i] = np.dot(temp, w);
                     f_fourier = f_fourier.real; # get rid of numerical imaginary noise
                     #print f_fourier
                     # Capon Inversion ######################
                     f_capon = np.zeros(shape=(Nr,1));
                     tic_capon = time.time();
                     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 #####################
                     # 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
                     tic_maxent = time.time();
                     lambda0 = 1e-5*np.ones(shape=(M,1)); # initial condition (can be set to anything)
                     toc_maxent = time.time()
                     elapsed_time_maxent = toc_maxent - tic_maxent;
                     # Whitened solution
                     def myfun(lambda1):
                         return y_hysell96(lambda1,gnz,sigma,F,G,Hr);
                     tic_maxEnt = time.time();
                     #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
                     lambda1 = root(myfun,lambda0, method='krylov',  tol=1e-14);
                     #print lambda1
                     #print lambda1.x
                     lambda1 = lambda1.x;
                     toc_maxEnt = time.time();
                     f_maxent = modelf(lambda1, Hr, F);
                     ystar = myfun(lambda1);
                     Lambda = np.sqrt(sum(lambda1**2.*sigma**2)/(4*G));
                     ep = np.multiply(-lambda1,sigma**2)/ (2*Lambda);
                     es = np.dot(Hr, f_maxent) - gnz; # should be same as ep
                     chi2 = np.sum((es/sigma)**2);
-                    # CS inversion using Iteratively Reweighted Least Squares (IRLS)-------------
+                    # --------- CS inversion using Iteratively Reweighted Least Squares (IRLS) -------------
                     # (Use Nr, thetar, gnz, and Hr from MaxEnt above)
-                    Psi = deb4_basis(Nr); ###### REPLACED BY LINEs BELOW (?)
+                    Psi = deb4_basis(Nr);
-                    print 'FINALLY!'
+                    # REMOVE THIS--------------------------------
-                    print Psi.shape
-                    # REMOVE THIS?--------------------------------
                     #wavelet1 = pywt.Wavelet('db4')
                     #Phi, Psi, x = wavelet1.wavefun(level=3)
                     # --------------------------------------------
                     # add "sum to 1" constraint
-                    #         H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
+                    H2 = np.concatenate( (Hr, np.ones(shape=(1,Nr))), axis=0 );
-                    #         N_temp = np.array([[Nr/Nt]]);
+                    N_temp = np.array([[Nr/Nt]]);
-                    #         g2 = np.concatenate( (gnz, N_temp), axis=0 );
+                    g2 = np.concatenate( (gnz, N_temp), axis=0 );
-                    #         H2 = H2.T.conj();
+                    #H2 = H2.T.conj();
-                    #         print 'H2 shape', H2.shape
-                    #         print 'Psi shape', Psi.shape
+                    #Psi = Psi.T.conj(); # to align matrices
-                    #         s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
+                    ####print 'H2 shape', H2.shape
+                    #####print 'Psi shape', Psi.shape
+                    A = np.dot(H2,Psi);
+                    s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
                     #         f_cs = Psi*s;
                     #
                     #         # plot
                     #         plot(thetar,f_cs,'r.-');
                     #         hold on;
                     #         plot(thetat,fact,'k-');
                     #         hold off;
                     # # # Scaling and shifting
                     # # # 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);
                     # To be able to perform dot product (align matrices) done below within the dot calculations
                     #f_fourier = f_fourier.T.conj()
                     #f_capon = f_capon.T.conj()
                     #f_maxent = f_maxent.T.conj()
                     #factr = factr.T.conj()
                     # 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,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.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 - 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()
                     # #     PLOT PARA COMPRESSED SENSING
                     # #
                     # #     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);
                     # # Store Results
                     corr[0, snri, Ni] = corr_fourier;
                     corr[1, snri, Ni] = corr_capon;
                     corr[2, snri, Ni] = corr_maxent;
                     #corr[3, snri, Ni] = corr_cs;
                     rmse[0,snri,Ni] = relrmse_fourier;
                     rmse[1,snri,Ni] = relrmse_capon;
                     rmse[2,snri,Ni] = relrmse_maxent;
                     #rmse[3,snri,Ni] = relrmse_cs;
                     corrc[0,snri,Ni] = corrc_fourier;
                     corrc[1,snri,Ni] = corrc_capon;
                     corrc[2,snri,Ni] = corrc_maxent;
                     #corrc[3,snri,Ni] = corrc_cs;
                     print 'Capon:\t', elapsed_time_capon, 'sec';
                     print 'Maxent:\t',elapsed_time_maxent, 'sec';
                     #print 'CS:\t%3.3f sec\n',elapsed_time_cs;
                     print (NN*(snri+1) + Ni), '/', (SNRdBvec.size*NN);
                     print corr
   No newline at end of file

trunk/src/src/cshift.py 0 -2

             '''
             Created on May 30, 2014
             @author: Yolian Amaro
             '''
             import numpy as np
             def cshift(x, m):
                 # Circular Shift
                 #
                 # USAGE:
                 #    y = cshift(x, m)
                 # INPUT:
                 #    x - N-point vector
                 #    m - amount of shift (pos=left, neg=right)
                 # OUTPUT:
                 #    y - vector x will be shifted by m samples to the left
                 #
                 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
                 # http://taco.poly.edu/WaveletSoftware/
                 N = x.size;
                 n = np.arange(N);
                 n = np.mod(n-m, N);
-                print x.shape
                 y = x[0,n];
                 return y

trunk/src/src/deb4_basis.py +2 -4

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

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

             '''
             Created on Jun 5, 2014
             @author: Yolian Amaro
             '''
             from sfb import *
             def idualtree(w, J, Fsf, sf):
                 # Inverse Dual-tree Complex DWT
                 #
                 # USAGE:
                 #    y = idualtree(w, J, Fsf, sf)
                 # INPUT:
                 #    w - DWT coefficients
                 #    J - number of stages
                 #    Fsf - synthesis filters for the last stage
                 #    sf - synthesis filters for preceeding stages
                 # OUTUT:
                 #    y - output signal
                 # See dualtree
                 #
                 # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY
                 # http://taco.poly.edu/WaveletSoftware/
                 # Tree 1
                 y1 = w[J][0];
                 for j in range (J-1, 0, -1):
                     y1 = sfb(y1, w[j][0], sf[0,0]);
                 y1 = sfb(y1, w[0][0], Fsf[0,0]);
                 # Tree 2
                 y2 = w[J][1];
                 for j in range (J-1, 0, -1):
-                    y2 = sfb(y2, w[j][2], sf[0,1]);
+                    y2 = sfb(y2, w[j][1], sf[0,1]);
                 y2 = sfb(y2, w[0][1], Fsf[0,1]);
                 # normalization
                 y = (y1 + y2)/np.sqrt(2);
                 return y

trunk/src/src/irls_dn.py +16 -13

             '''
             Created on May 27, 2014
             @author: Yolian Amaro
             '''
             #from scipy.sparse import eye
             from scipy import linalg
             import scipy.sparse as sps
             import numpy as np
             def irls_dn(A,b,p,lambda1):
                 # 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)
                 # Note to self:  I found that "warm-starting" didn't really help too much.
                 [M,N] = A.shape;
                 # Initialize and precompute:
                 eps = 1e-2; # damping parameter
-                [Q,R] = linalg.qr(A.T.conj(),0);
-                print A.shape
+                [Q,R] = linalg.qr(A.T.conj(), mode='economic');
-                print R.shape
-                print b.shape
                 c = linalg.solve(R.T.conj(),b); # will be used later also
-                u = Q*c; # minimum 2-norm solution
+                u = np.dot(Q,c); # minimum 2-norm solution
                 I = sps.eye(M);
                 #---------- not needed, defined above--------------
                 # Spacing of floating point numbers
                 #eps = np.spacing(1)
                 #--------------------------------------------------
                 # Loop until damping parameter is small enough
                 while (eps > 1e-7):
                     epschange = 0;
                     # Loop until it's time to change eps
                     while (~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)
                         # Update
                         w = (u**2 + eps)**(1-p/2);
                         # Empty temporary N x N matrix
                         temp = np.zeros(shape=(N,N))
+                        k = 0
                         # Sparse matrix
-                        for i in range (1, N):
+                        for i in range (0, N):
-                            for j in range (1,N):
+                            for j in range (0,N):
                                 if(i==j):
-                                    temp[i,j] = w
+                                    temp[i,j] = w[k]
+                                    k = k+1
                         # Compressed Sparse Matrix
                         W = sps.csr_matrix(temp); #Compressed Sparse Row matrix
                         WAT = W*A.T.conj();
-                        u_new = WAT * ( linalg.solve (A*WAT + lambda1*I), b);
+                        u_new = np.dot(WAT , linalg.solve(np.dot(A,WAT) + np.dot(lambda1,I), b));
                         # See if this subproblem is converging
                         delu = np.linalg.norm(u_new-u)/np.linalg.norm(u);
-                        epschange = delu < (np.sqrt(eps)/100);
+                        epschange = delu < (np.sqrt(eps)/100.0);
                         # Make update
                         u = u_new;
-                    eps = eps/10; # decrease eps
+                    eps = eps/10.0; # decrease eps
                     # Print info
-                    print 'eps =',eps;
+                    #print 'eps =',eps;
                 return u

trunk/src/src/irls_dn2.py +28 -7

             '''
             Created on May 30, 2014
             @author: Yolian Amaro
             '''
             from irls_dn import *
             from scipy.optimize import fsolve
+            import numpy as np
+            from scipy.optimize import root
             def irls_dn2(A,b,p,G):
                 # Minimize ||u||_p subject to ||A*u-b||_2^2 <= G   (0 < p <= 1)
                 # 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
                 # Start with a large lambda1, and do a line search until fidelity <= G.
                 # (Inversions with large lambda1 are really fast anyway).
                 # Then spin up fsolve to localize the root even better
                 # Line Search
-                alpha = 2; # Line search parameter
+                alpha = 2.0; # Line search parameter
                 lambda1 = 1e5; # What's a reasonable but safe initial guess?
                 u = irls_dn(A,b,p,lambda1);
-                fid = np.norm(A*u-b)**2;
+                fid = np.linalg.norm(np.dot(A,u)-b)**2;
                 print '----------------------------------\n';
                 while (fid >= G):
                     lambda1 = lambda1 / alpha; # Balance between speed and accuracy
                     u = irls_dn(A,b,p,lambda1);
-                    fid = np.norm(A*u-b)**2;
+                    fid = np.linalg.norm(np.dot(A,u)-b)**2;
-                    print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid;
+                    print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f\n' % fid;
                 # Refinement using fzero
                 lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
-                f = lambda lambda1: np.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G;
+                def myfun(lambda1):
+                    print "A = ", A.shape
+                    print "b = ", b.shape
+                    lambda1
+                    return np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G;
+                #f = lambda lambda1: np.linalg.norm(A*irls_dn(A,b,p,lambda1) - b)**2 - G; NOOOOOO
             #    opts = optimset('fzero');
             # #     opts.Display = 'iter';
             #     opts.Display = 'none';
             #     opts.TolX = 0.01*lambda1;
-                lambda1 = fsolve(f,lambda0);  # FALTA OPTIMIZE ESTO
+                #sol1 = fsolve(myfun,lambda0.ravel(), args=(), xtol=1e-14, maxfev=100000);
+                print "tolerancia=", 0.01*lambda1
+                #lambda1 = root(myfun,lambda0, method='krylov',  tol=0.01*lambda1);
+                print "lamda1=", lambda1
+                print "lambda0=", lambda0
+                lambda1 = fsolve(myfun,lambda0);  # FALTA OPTIMIZE ESTO
+                print "A = ", A.shape
+                print "b = ", b.shape
+                print "lambda1=", lambda1.shape
                 u = irls_dn(A,b,p,lambda1);
                 return u;

trunk/src/src/sfb.py +4 -12

             '''
             Created on Jun 5, 2014
             @author: Yolian Amaro
             '''
             from multirate import *
             import numpy as np
             from cshift import *
             def sfb(lo, hi, sf):
                 # 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/
                 N = 2*lo.size;
                 L = sf.size/2;
-                #print 'N', N
-                #print 'sf', sf
-                #print 'sf[:,0]', sf[:,0].shape
-                #print 'sf[:,1]', sf[:,1].shape
-                #print 'sbf hi', hi.shape
                 # Need to change format for upfirdn funct:
                 lo = lo.T.conj()
                 lo = lo.reshape(lo.size)
-                print 'sfb hi', hi
+                #print 'sfb hi', hi
                 # Need to change format for upfirdn funct:
                 hi = hi.T.conj()
                 hi = hi.reshape(hi.size)
                 #hi = hi.reshape(1, hi.size)
                 lo = upfirdn(lo, sf[:,0], 2, 1);
                 hi = upfirdn(hi, sf[:,1], 2, 1);
                 y = lo + hi;
                 y[0:L-1] = y[0:L-1] + y[N+ np.arange(0,L-1)]; #CHECK IF ARANGE IS CORRECT
                 y = y[0:N];
-                print 'y en sbf\n', y.shape
+                #print 'y en sbf\n', y.shape
                 y = y.reshape(1, y.size)
-                print 'y en sbf\n', y.shape
+                #print 'y en sbf\n', y.shape
                 y = cshift(y, 1-L/2);
                 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