Mercurial > hg > pycsalgos

--- a/pyCSalgos/Analysis.py	Mon Feb 27 20:41:11 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,127 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Dec 08 10:56:56 2011
-
-@author: ncleju
-"""
-
-import numpy
-import numpy.linalg
-
-from numpy.random import RandomState
-rng = RandomState()
-
-def Generate_Analysis_Operator(d, p):
-  # generate random tight frame with equal column norms
-  if p == d:
-    T = rng.randn(d,d);
-    [Omega, discard] = numpy.qr(T);
-  else:
-    Omega = rng.randn(p, d);
-    T = numpy.zeros((p, d));
-    tol = 1e-8;
-    max_j = 200;
-    j = 1;
-    while (sum(sum(abs(T-Omega))) > numpy.dot(tol,numpy.dot(p,d)) and j < max_j):
-        j = j + 1;
-        T = Omega;
-        [U, S, Vh] = numpy.linalg.svd(Omega);
-        V = Vh.T
-        #Omega = U * [eye(d); zeros(p-d,d)] * V';
-        Omega2 = numpy.dot(numpy.dot(U, numpy.concatenate((numpy.eye(d), numpy.zeros((p-d,d))))), V.transpose())
-        #Omega = diag(1./sqrt(diag(Omega*Omega')))*Omega;
-        Omega = numpy.dot(numpy.diag(1.0 / numpy.sqrt(numpy.diag(numpy.dot(Omega2,Omega2.transpose())))), Omega2)
-    #end
-    ##disp(j);
-    #end
-  return Omega
-
-
-def Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr):
-  #function [x0,y,M,LambdaMat] = Generate_Data_Known_Omega(Omega, d,p,m,k,noiselevel, numvectors, normstr)
-
-  # Building an analysis problem, which includes the ingredients:
-  #   - Omega - the analysis operator of size p*d
-  #   - M is anunderdetermined measurement matrix of size m*d (m<d)
-  #   - x0 is a vector of length d that satisfies ||Omega*x0||=p-k
-  #   - Lambda is the true location of these k zeros in Omega*x0
-  #   - a measurement vector y0=Mx0 is computed
-  #   - noise contaminated measurement vector y is obtained by
-  #     y = y0 + n where n is an additive gaussian noise with norm(n,2)/norm(y0,2) = noiselevel
-  # Added by Nic:
-  #   - Omega = analysis operator
-  #   - normstr: if 'l0', generate l0 sparse vector (unchanged). If 'l1',
-  #   generate a vector of Laplacian random variables (gamma) and
-  #   pseudoinvert to find x
-
-  # Omega is known as input parameter
-  #Omega=Generate_Analysis_Operator(d, p);
-  # Omega = randn(p,d);
-  # for i = 1:size(Omega,1)
-  #     Omega(i,:) = Omega(i,:) / norm(Omega(i,:));
-  # end
-
-  #Init
-  LambdaMat = numpy.zeros((k,numvectors))
-  x0 = numpy.zeros((d,numvectors))
-  y = numpy.zeros((m,numvectors))
-  realnoise = numpy.zeros((m,numvectors))
-
-  M = rng.randn(m,d);
-
-  #for i=1:numvectors
-  for i in range(0,numvectors):
-
-    # Generate signals
-    #if strcmp(normstr,'l0')
-    if normstr == 'l0':
-        # Unchanged
-
-        #Lambda=randperm(p);
-        Lambda = rng.permutation(int(p));
-        Lambda = numpy.sort(Lambda[0:k]);
-        LambdaMat[:,i] = Lambda; # store for output
-
-        # The signal is drawn at random from the null-space defined by the rows
-        # of the matreix Omega(Lambda,:)
-        [U,D,Vh] = numpy.linalg.svd(Omega[Lambda,:]);
-        V = Vh.T
-        NullSpace = V[:,k:];
-        #print numpy.dot(NullSpace, rng.randn(d-k,1)).shape
-        #print x0[:,i].shape
-        x0[:,i] = numpy.squeeze(numpy.dot(NullSpace, rng.randn(d-k,1)));
-        # Nic: add orthogonality noise
-        #     orthonoiseSNRdb = 6;
-        #     n = randn(p,1);
-        #     #x0(:,i) = x0(:,i) + n / norm(n)^2 * norm(x0(:,i))^2 / 10^(orthonoiseSNRdb/10);
-        #     n = n / norm(n)^2 * norm(Omega * x0(:,i))^2 / 10^(orthonoiseSNRdb/10);
-        #     x0(:,i) = pinv(Omega) * (Omega * x0(:,i) + n);
-
-    #elseif strcmp(normstr, 'l1')
-    elif normstr == 'l1':
-        print('Nic says: not implemented yet')
-        raise Exception('Nic says: not implemented yet')
-        #gamma = laprnd(p,1,0,1);
-        #x0(:,i) = Omega \ gamma;
-    else:
-        #error('normstr must be l0 or l1!');
-        print('Nic says: not implemented yet')
-        raise Exception('Nic says: not implemented yet')
-    #end
-
-    # Acquire measurements
-    y[:,i]  = numpy.dot(M, x0[:,i])
-
-    # Add noise
-    t_norm = numpy.linalg.norm(y[:,i],2)
-    n = numpy.squeeze(rng.randn(m, 1))
-    # In case n i just a number, nuit an array, norm() fails
-    if n.ndim == 0:
-      nnorm = abs(n)
-    else:
-      nnorm = numpy.linalg.norm(n, 2);
-    y[:,i] = y[:,i] + noiselevel * t_norm * n / nnorm
-    realnoise[:,i] = noiselevel * t_norm * n / nnorm
-  #end
-
-  return x0,y,M,LambdaMat,realnoise
--- a/pyCSalgos/RecomTST/RecommendedTST.py	Mon Feb 27 20:41:11 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,152 +0,0 @@
-
-import numpy as np
-import math
-
-#function beta = RecommendedTST(X,Y, nsweep,tol,xinitial,ro)
-def RecommendedTST(X, Y, nsweep=300, tol=0.00001, xinitial=None, ro=None):
-
-  # function beta=RecommendedTST(X,y, nsweep,tol,xinitial,ro)
-  # This function gets the measurement matrix and the measurements and
-  # the number of runs and applies the TST algorithm with optimally tuned parameters
-  # to the problem. For more information you may refer to the paper,
-  # "Optimally tuned iterative reconstruction algorithms for compressed
-  # sensing," by Arian Maleki and David L. Donoho.
-  #           X  : Measurement matrix; We assume that all the columns have
-  #               almost equal $\ell_2$ norms. The tunning has been done on
-  #               matrices with unit column norm.
-  #            y : output vector
-  #       nsweep : number of iterations. The default value is 300.
-  #          tol : if the relative prediction error i.e. ||Y-Ax||_2/ ||Y||_2 <
-  #               tol the algorithm will stop. If not provided the default
-  #               value is zero and tha algorithm will run for nsweep
-  #               iterations. The Default value is 0.00001.
-  #     xinitial : This is an optional parameter. It will be used as an
-  #                initialization of the algorithm. All the results mentioned
-  #                in the paper have used initialization at the zero vector
-  #                which is our default value. For default value you can enter
-  #                0 here.
-  #        ro    : This is a again an optional parameter. If not given the
-  #                algorithm will use the default optimal values. It specifies
-  #                the sparsity level. For the default value you may also used if
-  #                rostar=0;
-  #
-  # Outputs:
-  #      beta : the estimated coeffs.
-  #
-  # References:
-  # For more information about this algorithm and to find the papers about
-  # related algorithms like CoSaMP and SP please refer to the paper mentioned
-  # above and the references of that paper.
-  #
-  #---------------------
-  # Original Matab code:
-  #
-  #colnorm=mean((sum(X.^2)).^(.5));
-  #X=X./colnorm;
-  #Y=Y./colnorm;
-  #[n,p]=size(X);
-  #delta=n/p;
-  #if nargin<3
-  #    nsweep=300;
-  #end
-  #if nargin<4
-  #    tol=0.00001;
-  #end
-  #if nargin<5 | xinitial==0
-  #    xinitial = zeros(p,1);
-  #end
-  #if nargin<6 | ro==0
-  #    ro=0.044417*delta^2+ 0.34142*delta+0.14844;
-  #end
-  #
-  #
-  #k1=floor(ro*n);
-  #k2=floor(ro*n);
-  #
-  #
-  ##initialization
-  #x1=xinitial;
-  #I=[];
-  #
-  #for sweep=1:nsweep
-  #    r=Y-X*x1;
-  #    c=X'*r;
-  #    [csort, i_csort]=sort(abs(c));
-  #    I=union(I,i_csort(end-k2+1:end));
-  #    xt = X(:,I) \ Y;
-  #    [xtsort, i_xtsort]=sort(abs(xt));
-  #
-  #    J=I(i_xtsort(end-k1+1:end));
-  #    x1=zeros(p,1);
-  #    x1(J)=xt(i_xtsort(end-k1+1:end));
-  #    I=J;
-  #    if norm(Y-X*x1)/norm(Y)<tol
-  #        break
-  #    end
-  #
-  #end
-  #
-  #beta=x1;
-  #
-  # End of original Matab code
-  #----------------------------
-
-  #colnorm=mean((sum(X.^2)).^(.5));
-  colnorm = np.mean(np.sqrt((X**2).sum(0)))
-  X = X / colnorm
-  Y = Y / colnorm
-  [n,p] = X.shape
-  delta = float(n) / p
-  #  if nargin<3
-  #      nsweep=300;
-  #  end
-  #  if nargin<4
-  #      tol=0.00001;
-  #  end
-  #  if nargin<5 | xinitial==0
-  #      xinitial = zeros(p,1);
-  #  end
-  #  if nargin<6 | ro==0
-  #      ro=0.044417*delta^2+ 0.34142*delta+0.14844;
-  #  end
-  if xinitial is None:
-    xinitial = np.zeros(p)
-  if ro == None:
-    ro = 0.044417*delta**2 + 0.34142*delta + 0.14844
-
-  k1 = math.floor(ro*n)
-  k2 = math.floor(ro*n)
-
-  #initialization
-  x1 = xinitial.copy()
-  I = []
-
-  for sweep in np.arange(nsweep):
-      r = Y - np.dot(X,x1)
-      c = np.dot(X.T, r)
-      #[csort, i_csort] = np.sort(np.abs(c))
-      i_csort = np.argsort(np.abs(c))
-      #I = numpy.union1d(I , i_csort(end-k2+1:end))
-      I = np.union1d(I , i_csort[-k2:])
-      #xt = X[:,I] \ Y
-      # Make sure X[:,np.int_(I)] is a 2-dimensional matrix even if I has a single value (and therefore yields a column)
-      if I.size is 1:
-        a = np.reshape(X[:,np.int_(I)],(X.shape[0],1))
-      else:
-        a = X[:,np.int_(I)]
-      xt = np.linalg.lstsq(a, Y)[0]
-      #[xtsort, i_xtsort] = np.sort(np.abs(xt))
-      i_xtsort = np.argsort(np.abs(xt))
-
-      J = I[i_xtsort[-k1:]]
-      x1 = np.zeros(p)
-      x1[np.int_(J)] = xt[i_xtsort[-k1:]]
-      I = J.copy()
-      if np.linalg.norm(Y-np.dot(X,x1)) / np.linalg.norm(Y) < tol:
-          break
-      #end
-
-  #end
-
-  return x1.copy()
-
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pyCSalgos/TST/RecommendedTST.py	Fri Mar 09 16:25:31 2012 +0200
@@ -0,0 +1,152 @@
+
+import numpy as np
+import math
+
+#function beta = RecommendedTST(X,Y, nsweep,tol,xinitial,ro)
+def RecommendedTST(X, Y, nsweep=300, tol=0.00001, xinitial=None, ro=None):
+
+  # function beta=RecommendedTST(X,y, nsweep,tol,xinitial,ro)
+  # This function gets the measurement matrix and the measurements and
+  # the number of runs and applies the TST algorithm with optimally tuned parameters
+  # to the problem. For more information you may refer to the paper,
+  # "Optimally tuned iterative reconstruction algorithms for compressed
+  # sensing," by Arian Maleki and David L. Donoho.
+  #           X  : Measurement matrix; We assume that all the columns have
+  #               almost equal $\ell_2$ norms. The tunning has been done on
+  #               matrices with unit column norm.
+  #            y : output vector
+  #       nsweep : number of iterations. The default value is 300.
+  #          tol : if the relative prediction error i.e. ||Y-Ax||_2/ ||Y||_2 <
+  #               tol the algorithm will stop. If not provided the default
+  #               value is zero and tha algorithm will run for nsweep
+  #               iterations. The Default value is 0.00001.
+  #     xinitial : This is an optional parameter. It will be used as an
+  #                initialization of the algorithm. All the results mentioned
+  #                in the paper have used initialization at the zero vector
+  #                which is our default value. For default value you can enter
+  #                0 here.
+  #        ro    : This is a again an optional parameter. If not given the
+  #                algorithm will use the default optimal values. It specifies
+  #                the sparsity level. For the default value you may also used if
+  #                rostar=0;
+  #
+  # Outputs:
+  #      beta : the estimated coeffs.
+  #
+  # References:
+  # For more information about this algorithm and to find the papers about
+  # related algorithms like CoSaMP and SP please refer to the paper mentioned
+  # above and the references of that paper.
+  #
+  #---------------------
+  # Original Matab code:
+  #
+  #colnorm=mean((sum(X.^2)).^(.5));
+  #X=X./colnorm;
+  #Y=Y./colnorm;
+  #[n,p]=size(X);
+  #delta=n/p;
+  #if nargin<3
+  #    nsweep=300;
+  #end
+  #if nargin<4
+  #    tol=0.00001;
+  #end
+  #if nargin<5 | xinitial==0
+  #    xinitial = zeros(p,1);
+  #end
+  #if nargin<6 | ro==0
+  #    ro=0.044417*delta^2+ 0.34142*delta+0.14844;
+  #end
+  #
+  #
+  #k1=floor(ro*n);
+  #k2=floor(ro*n);
+  #
+  #
+  ##initialization
+  #x1=xinitial;
+  #I=[];
+  #
+  #for sweep=1:nsweep
+  #    r=Y-X*x1;
+  #    c=X'*r;
+  #    [csort, i_csort]=sort(abs(c));
+  #    I=union(I,i_csort(end-k2+1:end));
+  #    xt = X(:,I) \ Y;
+  #    [xtsort, i_xtsort]=sort(abs(xt));
+  #
+  #    J=I(i_xtsort(end-k1+1:end));
+  #    x1=zeros(p,1);
+  #    x1(J)=xt(i_xtsort(end-k1+1:end));
+  #    I=J;
+  #    if norm(Y-X*x1)/norm(Y)<tol
+  #        break
+  #    end
+  #
+  #end
+  #
+  #beta=x1;
+  #
+  # End of original Matab code
+  #----------------------------
+
+  #colnorm=mean((sum(X.^2)).^(.5));
+  colnorm = np.mean(np.sqrt((X**2).sum(0)))
+  X = X / colnorm
+  Y = Y / colnorm
+  [n,p] = X.shape
+  delta = float(n) / p
+  #  if nargin<3
+  #      nsweep=300;
+  #  end
+  #  if nargin<4
+  #      tol=0.00001;
+  #  end
+  #  if nargin<5 | xinitial==0
+  #      xinitial = zeros(p,1);
+  #  end
+  #  if nargin<6 | ro==0
+  #      ro=0.044417*delta^2+ 0.34142*delta+0.14844;
+  #  end
+  if xinitial is None:
+    xinitial = np.zeros(p)
+  if ro == None:
+    ro = 0.044417*delta**2 + 0.34142*delta + 0.14844
+
+  k1 = math.floor(ro*n)
+  k2 = math.floor(ro*n)
+
+  #initialization
+  x1 = xinitial.copy()
+  I = []
+
+  for sweep in np.arange(nsweep):
+      r = Y - np.dot(X,x1)
+      c = np.dot(X.T, r)
+      #[csort, i_csort] = np.sort(np.abs(c))
+      i_csort = np.argsort(np.abs(c))
+      #I = numpy.union1d(I , i_csort(end-k2+1:end))
+      I = np.union1d(I , i_csort[-k2:])
+      #xt = X[:,I] \ Y
+      # Make sure X[:,np.int_(I)] is a 2-dimensional matrix even if I has a single value (and therefore yields a column)
+      if I.size is 1:
+        a = np.reshape(X[:,np.int_(I)],(X.shape[0],1))
+      else:
+        a = X[:,np.int_(I)]
+      xt = np.linalg.lstsq(a, Y)[0]
+      #[xtsort, i_xtsort] = np.sort(np.abs(xt))
+      i_xtsort = np.argsort(np.abs(xt))
+
+      J = I[i_xtsort[-k1:]]
+      x1 = np.zeros(p)
+      x1[np.int_(J)] = xt[i_xtsort[-k1:]]
+      I = J.copy()
+      if np.linalg.norm(Y-np.dot(X,x1)) / np.linalg.norm(Y) < tol:
+          break
+      #end
+
+  #end
+
+  return x1.copy()
+
--- a/pyCSalgos/algos.py	Mon Feb 27 20:41:11 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,138 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Dec 07 14:06:13 2011
-
-@author: ncleju
-"""
-
-import numpy
-import pyCSalgos
-import pyCSalgos.GAP.GAP
-import pyCSalgos.BP.l1qc
-import pyCSalgos.BP.l1qec
-import pyCSalgos.SL0.SL0_approx
-import pyCSalgos.OMP.omp_QR
-import pyCSalgos.RecomTST.RecommendedTST
-import pyCSalgos.NESTA.NESTA
-
-#==========================
-# Algorithm functions
-#==========================
-def run_gap(y,M,Omega,epsilon):
-  gapparams = {"num_iteration" : 1000,\
-                   "greedy_level" : 0.9,\
-                   "stopping_coefficient_size" : 1e-4,\
-                   "l2solver" : 'pseudoinverse',\
-                   "noise_level": epsilon}
-  return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,numpy.zeros(Omega.shape[1]))[0]
-
-def run_bp_analysis(y,M,Omega,epsilon):
-
-  N,n = Omega.shape
-  D = numpy.linalg.pinv(Omega)
-  U,S,Vt = numpy.linalg.svd(D)
-  Aeps = numpy.dot(M,D)
-  Aexact = Vt[-(N-n):,:]
-  # We don't ned any aggregate matrices anymore
-
-  x0 = numpy.zeros(N)
-  return numpy.dot(D , pyCSalgos.BP.l1qec.l1qec_logbarrier(x0,Aeps,Aeps.T,y,epsilon,Aexact,Aexact.T,numpy.zeros(N-n)))
-
-def run_sl0_analysis(y,M,Omega,epsilon):
-
-  N,n = Omega.shape
-  D = numpy.linalg.pinv(Omega)
-  U,S,Vt = numpy.linalg.svd(D)
-  Aeps = numpy.dot(M,D)
-  Aexact = Vt[-(N-n):,:]
-  # We don't ned any aggregate matrices anymore
-
-  sigmamin = 0.001
-  sigma_decrease_factor = 0.5
-  mu_0 = 2
-  L = 10
-  return numpy.dot(D , pyCSalgos.SL0.SL0_approx.SL0_approx_analysis(Aeps,Aexact,y,epsilon,sigmamin,sigma_decrease_factor,mu_0,L))
-
-def run_nesta(y,M,Omega,epsilon):
-
-  U,S,V = numpy.linalg.svd(M, full_matrices = True)
-  V = V.T         # Make like Matlab
-  m,n = M.shape   # Make like Matlab
-  S = numpy.hstack((numpy.diag(S), numpy.zeros((m,n-m))))
-
-  opt_muf = 1e-3
-  optsUSV = {'U':U, 'S':S, 'V':V}
-  opts = {'U':Omega, 'Ut':Omega.T.copy(), 'USV':optsUSV, 'TolVar':1e-5, 'Verbose':0}
-  return pyCSalgos.NESTA.NESTA.NESTA(M, None, y, opt_muf, epsilon, opts)[0]
-
-
-def run_sl0(y,M,Omega,D,U,S,Vt,epsilon,lbd):
-
-  N,n = Omega.shape
-  #D = numpy.linalg.pinv(Omega)
-  #U,S,Vt = numpy.linalg.svd(D)
-  aggDupper = numpy.dot(M,D)
-  aggDlower = Vt[-(N-n):,:]
-  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
-  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
-
-  sigmamin = 0.001
-  sigma_decrease_factor = 0.5
-  mu_0 = 2
-  L = 10
-  return pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigmamin,sigma_decrease_factor,mu_0,L)
-
-def run_bp(y,M,Omega,D,U,S,Vt,epsilon,lbd):
-
-  N,n = Omega.shape
-  #D = numpy.linalg.pinv(Omega)
-  #U,S,Vt = numpy.linalg.svd(D)
-  aggDupper = numpy.dot(M,D)
-  aggDlower = Vt[-(N-n):,:]
-  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
-  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
-
-  x0 = numpy.zeros(N)
-  return pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon)
-
-def run_ompeps(y,M,Omega,D,U,S,Vt,epsilon,lbd):
-
-  N,n = Omega.shape
-  #D = numpy.linalg.pinv(Omega)
-  #U,S,Vt = numpy.linalg.svd(D)
-  aggDupper = numpy.dot(M,D)
-  aggDlower = Vt[-(N-n):,:]
-  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
-  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
-
-  opts = dict()
-  opts['stopCrit'] = 'mse'
-  opts['stopTol'] = epsilon**2 / aggy.size
-  return pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0]
-
-def run_tst(y,M,Omega,D,U,S,Vt,epsilon,lbd):
-
-  N,n = Omega.shape
-  #D = numpy.linalg.pinv(Omega)
-  #U,S,Vt = numpy.linalg.svd(D)
-  aggDupper = numpy.dot(M,D)
-  aggDlower = Vt[-(N-n):,:]
-  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
-  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
-
-  nsweep = 300
-  tol = epsilon / numpy.linalg.norm(aggy)
-  return pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep=nsweep, tol=tol)
-
-
-#==========================
-# Define tuples (algorithm function, name)
-#==========================
-gap = (run_gap, 'GAP')
-sl0 = (run_sl0, 'SL0a')
-sl0analysis = (run_sl0_analysis, 'SL0a2')
-bpanalysis = (run_bp_analysis, 'BPa2')
-nesta = (run_nesta, 'NESTA')
-bp  = (run_bp, 'BP')
-ompeps = (run_ompeps, 'OMPeps')
-tst = (run_tst, 'TST')
\ No newline at end of file