Index: trunk/src/src/FSfarras.py =================================================================== diff --git a/trunk/src/src/FSfarras.py b/trunk/src/src/FSfarras.py new file mode 10644 --- /dev/null (revision 0) +++ b/trunk/src/src/FSfarras.py (revision 18) @@ -0,0 +1,65 @@ +''' +Created on May 26, 2014 + +@author: Yolian Amaro +''' + +import pywt +import numpy as np + +def FSfarras(): + #function [af, sf] = FSfarras + + # 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 + + + + 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] + ] ); + + 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] + ]); + + #print a2.shape + + af = np.array([ [a1,a2] ], dtype=object) + + s1 = a1[::-1] + s2 = a2[::-1] + + sf = np.array([ [s1,s2] ], dtype=object) + + return af, sf Index: trunk/src/src/RICS_paper_compare_noise.py =================================================================== diff --git a/trunk/src/src/RICS_paper_compare_noise.py b/trunk/src/src/RICS_paper_compare_noise.py --- a/trunk/src/src/RICS_paper_compare_noise.py (revision 17) +++ b/trunk/src/src/RICS_paper_compare_noise.py (revision 18) @@ -9,20 +9,17 @@ #---------------------------------------------------------- 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 fsolve from scipy.optimize import root import pywt +from irls_dn2 import * ## Calculate Forward Model @@ -138,8 +135,9 @@ 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; +#plt.figure() +#plt.plot(thetat,fact,'r'); +#plt.plot(thetar,factr,'k.'); #xlim([min(thetat) max(thetat)]); #x = np.linspace(thetat.min(), thetat.max) ???? @@ -197,13 +195,14 @@ 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); @@ -230,16 +229,17 @@ f_capon = np.zeros(shape=(Nr,1)); - #tic_capon = time.time(); + 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; + + toc_capon = time.time() + + elapsed_time_capon = toc_capon - tic_capon; f_capon = f_capon.real; # get rid of numerical imaginary noise @@ -278,10 +278,15 @@ 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 + 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); @@ -304,24 +309,28 @@ es = np.dot(Hr, f_maxent) - gnz; # should be same as ep chi2 = np.sum((es/sigma)**2); - - + # CS inversion using irls ######################## # (Use Nr, thetar, gnz, and Hr from MaxEnt above) - #Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?) - - wavelet1 = pywt.Wavelet('db4') - Phi, Psi, x = wavelet1.wavefun(level=3) + Psi = deb4_basis(Nr); ###### REPLACED BY LINE BELOW (?) + + # 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 ); N_temp = np.array([[Nr/Nt]]); g2 = np.concatenate( (gnz, N_temp), axis=0 ); - - - #s = irls_dn2(H2*Psi,g2,0.5,G); + H2 = H2.T.conj(); + + print 'H2 shape', H2.shape + print 'Psi shape', Psi.shape + + s = irls_dn2(H2*Psi,g2,0.5,G); # f_cs = Psi*s; # # # plot @@ -331,100 +340,129 @@ # 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)); + # # # Scaling and shifting + # # # Only necessary for capon solution - 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 + + 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 Index: trunk/src/src/deb4_basis.py =================================================================== diff --git a/trunk/src/src/deb4_basis.py b/trunk/src/src/deb4_basis.py new file mode 10644 --- /dev/null (revision 0) +++ b/trunk/src/src/deb4_basis.py (revision 18) @@ -0,0 +1,38 @@ +''' +Created on May 26, 2014 + +@author: Yolian Amaro +''' + +import numpy as np +import FSfarras +import dualfilt1 + +def deb4_basis(N): + + Psi = np.zeros(shape=(N,2*N+1)); + idx = 1; + + J = 4; + [Faf, Fsf] = FSfarras; + [af, sf] = dualfilt1; +# # +# # # compute transform of zero vector +# # x = zeros(1,N); +# # w = dualtree(x, J, Faf, af); +# # +# # # Uses both real and imaginary wavelets +# # for i in range (1, J+1): +# # for j in range (1, 2): +# # for k in range (1, (w[i][j]).size): +# # w[i][j](k) = 1; +# # y = idualtree(w, J, Fsf, sf); +# # w[i][j](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); +# +# return Psi \ No newline at end of file Index: trunk/src/src/dualfilt1.py =================================================================== diff --git a/trunk/src/src/dualfilt1.py b/trunk/src/src/dualfilt1.py new file mode 10644 --- /dev/null (revision 0) +++ b/trunk/src/src/dualfilt1.py (revision 18) @@ -0,0 +1,65 @@ +''' +Created on May 29, 2014 + +@author: Yolian Amaro +''' + +import numpy as np + +def dualfilt1(): + + # Kingsbury Q-filters for the dual-tree complex DWT + # + # USAGE: + # [af, sf] = dualfilt1 + # OUTPUT: + # af{i}, i = 1,2 - analysis filters for tree i + # sf{i}, i = 1,2 - synthesis filters for tree i + # note: af{2} is the reverse of af{1} + # REFERENCE: + # N. G. Kingsbury, "A dual-tree complex wavelet + # transform with improved orthogonality and symmetry + # properties", Proceedings of the IEEE Int. Conf. on + # Image Proc. (ICIP), 2000 + # See dualtree + # + # WAVELET SOFTWARE AT POLYTECHNIC UNIVERSITY, BROOKLYN, NY + # http://taco.poly.edu/WaveletSoftware/ + + # These cofficients are rounded to 8 decimal places. + + a1 = np.array([ + [ 0.03516384000000, 0], + [ 0, 0], + [-0.08832942000000, -0.11430184000000], + [ 0.23389032000000, 0], + [ 0.76027237000000, 0.58751830000000], + [ 0.58751830000000, -0.76027237000000], + [ 0, 0.23389032000000], + [-0.11430184000000, 0.08832942000000], + [ 0, 0], + [ 0, -0.03516384000000] + ]); + + a2 = np.array([ + [ 0, -0.03516384000000], + [ 0, 0], + [-0.11430184000000, 0.08832942000000], + [ 0, 0.23389032000000], + [ 0.58751830000000, -0.76027237000000], + [ 0.76027237000000, 0.58751830000000], + [ 0.23389032000000, 0], + [ -0.08832942000000, -0.11430184000000], + [ 0, 0], + [ 0.03516384000000, 0] + ]); + + af = np.array([ [a1,a2] ], dtype=object) + + s1 = a1[::-1] + s2 = a2[::-1] + + sf = np.array([ [s1,s2] ], dtype=object) + + + return af, sf Index: trunk/src/src/irls_dn.py =================================================================== diff --git a/trunk/src/src/irls_dn.py b/trunk/src/src/irls_dn.py new file mode 10644 --- /dev/null (revision 0) +++ b/trunk/src/src/irls_dn.py (revision 18) @@ -0,0 +1,80 @@ +''' +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 + 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 + I = sps.eye(M); + + # 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)) + + # Sparse matrix + for i in range (1, N): + for j in range (1,N): + if(i==j): + temp[i,j] = w + + # 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); + + # See if this subproblem is converging + delu = np.linalg.norm(u_new-u)/np.linalg.norm(u); + epschange = delu < (np.sqrt(eps)/100); + + # Make update + u = u_new; + + + eps = eps/10; # decrease eps + # Print info + print 'eps =',eps; + + return u + + + Index: trunk/src/src/irls_dn2.py =================================================================== diff --git a/trunk/src/src/irls_dn2.py b/trunk/src/src/irls_dn2.py new file mode 10644 --- /dev/null (revision 0) +++ b/trunk/src/src/irls_dn2.py (revision 18) @@ -0,0 +1,50 @@ +''' +Created on May 27, 2014 + +@author: Yolian Amaro +''' + +from irls_dn import * + +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 fzero to localize the root even better + + # Line Search + + alpha = 2; # 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; +# +# 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; +# print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid); +# +# # Refinement using fzero +# lambda10 = [lambda1 lambda1*alpha]; # interval with zero-crossing +# f = @(lambda1) np.norm(A*irls_dn(A,b,p,lambda1) - b)^2 - G; +# opts = optimset('fzero'); +# # opts.Display = 'iter'; +# opts.Display = 'none'; +# opts.TolX = 0.01*lambda1; +# lambda1 = fzero(f,lambda10,opts); +# +# u = irls_dn(A,b,p,lambda1); +# +# +# return u; +