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 20)
+++ b/trunk/src/src/RICS_paper_compare_noise.py	(revision 21)
@@ -310,30 +310,32 @@
         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 (?)
-        
-        print 'FINALLY!'
-        print Psi.shape
-        
-        # REMOVE THIS?--------------------------------
+        Psi = deb4_basis(Nr);
+        
+        # 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 );
-        #         H2 = H2.T.conj();
-        #         
-        #         print 'H2 shape', H2.shape
-        #         print 'Psi shape', Psi.shape
-        #         
-        #         s = irls_dn2(np.dot(H2,Psi),g2,0.5,G);
+        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 );
+        
+        #H2 = H2.T.conj();
+        
+        #Psi = Psi.T.conj(); # to align matrices
+                 
+        ####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
@@ -345,8 +347,6 @@
     
         # # # Scaling and shifting
         # # # Only necessary for capon solution
-     
-    
         f_capon = f_capon/np.max(f_capon)*np.max(fact);
     
         
Index: trunk/src/src/cshift.py
===================================================================
diff --git a/trunk/src/src/cshift.py b/trunk/src/src/cshift.py
--- a/trunk/src/src/cshift.py	(revision 20)
+++ b/trunk/src/src/cshift.py	(revision 21)
@@ -25,8 +25,6 @@
     N = x.size;
     n = np.arange(N);
     n = np.mod(n-m, N);
-
-    print x.shape
     
     y = x[0,n];
     
Index: trunk/src/src/deb4_basis.py
===================================================================
diff --git a/trunk/src/src/deb4_basis.py b/trunk/src/src/deb4_basis.py
--- a/trunk/src/src/deb4_basis.py	(revision 20)
+++ b/trunk/src/src/deb4_basis.py	(revision 21)
@@ -13,8 +13,7 @@
 def deb4_basis(N):
  
     Psi = np.zeros(shape=(N,2*N+1));
-    idx = 1;
- 
+    idx = 0;
     J = 4;
     [Faf, Fsf] = FSfarras();
     [af, sf] = dualfilt1();
@@ -22,7 +21,6 @@
     # 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):
@@ -36,6 +34,6 @@
                 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
\ No newline at end of file
Index: trunk/src/src/idualtree.py
===================================================================
diff --git a/trunk/src/src/idualtree.py b/trunk/src/src/idualtree.py
--- a/trunk/src/src/idualtree.py	(revision 20)
+++ b/trunk/src/src/idualtree.py	(revision 21)
@@ -36,7 +36,7 @@
     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]);
     
Index: trunk/src/src/irls_dn.py
===================================================================
diff --git a/trunk/src/src/irls_dn.py b/trunk/src/src/irls_dn.py
--- a/trunk/src/src/irls_dn.py	(revision 20)
+++ b/trunk/src/src/irls_dn.py	(revision 21)
@@ -18,16 +18,16 @@
     # 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
+    
+    [Q,R] = linalg.qr(A.T.conj(), mode='economic');
+    
+    
     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--------------
@@ -51,30 +51,33 @@
             # Empty temporary N x N matrix
             temp = np.zeros(shape=(N,N))
             
+            k = 0
             # Sparse matrix
-            for i in range (1, N):
-                for j in range (1,N):
+            for i in range (0, 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
  
Index: trunk/src/src/irls_dn2.py
===================================================================
diff --git a/trunk/src/src/irls_dn2.py b/trunk/src/src/irls_dn2.py
--- a/trunk/src/src/irls_dn2.py	(revision 20)
+++ b/trunk/src/src/irls_dn2.py	(revision 21)
@@ -6,6 +6,8 @@
 
 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):
  
@@ -22,23 +24,29 @@
      
     # 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;
-        print 'lambda1 = #2e \t ||A*u-b||^2 = #.1f\n',lambda1,fid;
+        fid = np.linalg.norm(np.dot(A,u)-b)**2;
+        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');
@@ -46,7 +54,20 @@
 #     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);
 
Index: trunk/src/src/sfb.py
===================================================================
diff --git a/trunk/src/src/sfb.py b/trunk/src/src/sfb.py
--- a/trunk/src/src/sfb.py	(revision 20)
+++ b/trunk/src/src/sfb.py	(revision 21)
@@ -29,27 +29,19 @@
     
     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);
@@ -57,10 +49,10 @@
     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);