Mercurial > hg > pycsalgos

changeset 52:768b57e446ab

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Created Analysis.py and working
author nikcleju
date Thu, 08 Dec 2011 09:05:04 +0000
parents eb4c66488ddf
children 991794ff9544
files pyCSalgos/Analysis.py pyCSalgos/GAP/GAP.py scripts/ABSapprox.py scripts/algos.py scripts/stdparams.py
diffstat 5 files changed, 211 insertions(+), 200 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/pyCSalgos/Analysis.py	Thu Dec 08 09:05:04 2011 +0000
@@ -0,0 +1,127 @@
+# -*- 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/GAP/GAP.py	Wed Dec 07 12:44:19 2011 +0000
+++ b/pyCSalgos/GAP/GAP.py	Thu Dec 08 09:05:04 2011 +0000
@@ -5,134 +5,13 @@
 @author: ncleju
 """
 
-#from numpy import *
-#from scipy import *
+
 import numpy
 import numpy.linalg
 import scipy as sp
-#import scipy.stats #from scipy import stats
-#import scipy.linalg #from scipy import lnalg
+
 import math
 
-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
-
-#####################
 
 #function [xhat, arepr, lagmult] = ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, ilagmult, params)
 def ArgminOperL2Constrained(y, M, MH, Omega, OmegaH, Lambdahat, xinit, ilagmult, params):
--- a/scripts/ABSapprox.py	Wed Dec 07 12:44:19 2011 +0000
+++ b/scripts/ABSapprox.py	Thu Dec 08 09:05:04 2011 +0000
@@ -11,16 +11,7 @@
 import os
 
 import stdparams
-
-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
-
+import pyCSalgos.Analysis
 
 #==========================
 # Pool initializer function (multiprocessing)
@@ -279,6 +270,8 @@
   
   return mrelerrN,mrelerrL
 
+
+
 def generateData(d,sigma,delta,rho,numvects,SNRdb):
 
   # Process parameters
@@ -288,7 +281,7 @@
   l = round(d - rho*m);
   
   # Generate Omega and data based on parameters
-  Omega = pyCSalgos.GAP.GAP.Generate_Analysis_Operator(d, p);
+  Omega = pyCSalgos.Analysis.Generate_Analysis_Operator(d, p);
   # Optionally make Omega more coherent
   U,S,Vt = np.linalg.svd(Omega);
   Sdnew = S * (1+np.arange(S.size)) # Make D coherent, not Omega!
@@ -296,10 +289,11 @@
   Omega = np.dot(U , np.dot(Snew,Vt))
 
   # Generate data  
-  x0,y,M,Lambda,realnoise = pyCSalgos.GAP.GAP.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');
+  x0,y,M,Lambda,realnoise = pyCSalgos.Analysis.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');
   
   return Omega,x0,y,M,realnoise
-  
+
+
 def runsingleexampledebug():
   d = 50.0;
   sigma = 2.0
--- a/scripts/algos.py	Wed Dec 07 12:44:19 2011 +0000
+++ b/scripts/algos.py	Thu Dec 08 09:05:04 2011 +0000
@@ -5,6 +5,16 @@
 @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
 #==========================
@@ -14,26 +24,26 @@
                    "stopping_coefficient_size" : 1e-4,\
                    "l2solver" : 'pseudoinverse',\
                    "noise_level": epsilon}
-  return pyCSalgos.GAP.GAP.GAP(y,M,M.T,Omega,Omega.T,gapparams,np.zeros(Omega.shape[1]))[0]
+  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 = np.linalg.pinv(Omega)
-  U,S,Vt = np.linalg.svd(D)
-  Aeps = np.dot(M,D)
+  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 = np.zeros(N)
-  return np.dot(D , pyCSalgos.BP.l1qec.l1qec_logbarrier(x0,Aeps,Aeps.T,y,epsilon,Aexact,Aexact.T,np.zeros(N-n)))
+  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 = np.linalg.pinv(Omega)
-  U,S,Vt = np.linalg.svd(D)
-  Aeps = np.dot(M,D)
+  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
   
@@ -41,14 +51,14 @@
   sigma_decrease_factor = 0.5
   mu_0 = 2
   L = 10
-  return np.dot(D , pyCSalgos.SL0.SL0_approx.SL0_approx_analysis(Aeps,Aexact,y,epsilon,sigmamin,sigma_decrease_factor,mu_0,L))
+  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 = np.linalg.svd(M, full_matrices = True)
+  U,S,V = numpy.linalg.svd(M, full_matrices = True)
   V = V.T         # Make like Matlab
   m,n = M.shape   # Make like Matlab
-  S = np.hstack((np.diag(S), np.zeros((m,n-m))))  
+  S = numpy.hstack((numpy.diag(S), numpy.zeros((m,n-m))))  
 
   opt_muf = 1e-3
   optsUSV = {'U':U, 'S':S, 'V':V}
@@ -59,12 +69,12 @@
 def run_sl0(y,M,Omega,D,U,S,Vt,epsilon,lbd):
   
   N,n = Omega.shape
-  #D = np.linalg.pinv(Omega)
-  #U,S,Vt = np.linalg.svd(D)
-  aggDupper = np.dot(M,D)
+  #D = numpy.linalg.pinv(Omega)
+  #U,S,Vt = numpy.linalg.svd(D)
+  aggDupper = numpy.dot(M,D)
   aggDlower = Vt[-(N-n):,:]
-  aggD = np.concatenate((aggDupper, lbd * aggDlower))
-  aggy = np.concatenate((y, np.zeros(N-n)))
+  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
+  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
   
   sigmamin = 0.001
   sigma_decrease_factor = 0.5
@@ -75,25 +85,25 @@
 def run_bp(y,M,Omega,D,U,S,Vt,epsilon,lbd):
   
   N,n = Omega.shape
-  #D = np.linalg.pinv(Omega)
-  #U,S,Vt = np.linalg.svd(D)
-  aggDupper = np.dot(M,D)
+  #D = numpy.linalg.pinv(Omega)
+  #U,S,Vt = numpy.linalg.svd(D)
+  aggDupper = numpy.dot(M,D)
   aggDlower = Vt[-(N-n):,:]
-  aggD = np.concatenate((aggDupper, lbd * aggDlower))
-  aggy = np.concatenate((y, np.zeros(N-n)))
+  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
+  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
 
-  x0 = np.zeros(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 = np.linalg.pinv(Omega)
-  #U,S,Vt = np.linalg.svd(D)
-  aggDupper = np.dot(M,D)
+  #D = numpy.linalg.pinv(Omega)
+  #U,S,Vt = numpy.linalg.svd(D)
+  aggDupper = numpy.dot(M,D)
   aggDlower = Vt[-(N-n):,:]
-  aggD = np.concatenate((aggDupper, lbd * aggDlower))
-  aggy = np.concatenate((y, np.zeros(N-n)))
+  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
+  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
   
   opts = dict()
   opts['stopCrit'] = 'mse'
@@ -103,15 +113,15 @@
 def run_tst(y,M,Omega,D,U,S,Vt,epsilon,lbd):
   
   N,n = Omega.shape
-  #D = np.linalg.pinv(Omega)
-  #U,S,Vt = np.linalg.svd(D)
-  aggDupper = np.dot(M,D)
+  #D = numpy.linalg.pinv(Omega)
+  #U,S,Vt = numpy.linalg.svd(D)
+  aggDupper = numpy.dot(M,D)
   aggDlower = Vt[-(N-n):,:]
-  aggD = np.concatenate((aggDupper, lbd * aggDlower))
-  aggy = np.concatenate((y, np.zeros(N-n)))
+  aggD = numpy.concatenate((aggDupper, lbd * aggDlower))
+  aggy = numpy.concatenate((y, numpy.zeros(N-n)))
   
   nsweep = 300
-  tol = epsilon / np.linalg.norm(aggy)
+  tol = epsilon / numpy.linalg.norm(aggy)
   return pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep=nsweep, tol=tol)
 
 
--- a/scripts/stdparams.py	Wed Dec 07 12:44:19 2011 +0000
+++ b/scripts/stdparams.py	Thu Dec 08 09:05:04 2011 +0000
@@ -5,6 +5,7 @@
 @author: ncleju
 """
 
+import numpy
 from algos import *
 
 #==========================
@@ -23,16 +24,16 @@
   
   d = 50.0
   sigma = 2.0
-  deltas = np.array([0.05, 0.45, 0.95])
-  rhos = np.array([0.05, 0.45, 0.95])
-  #deltas = np.array([0.95])
-  #deltas = np.arange(0.05,1.,0.05)
-  #rhos = np.array([0.05])
+  deltas = numpy.array([0.05, 0.45, 0.95])
+  rhos = numpy.array([0.05, 0.45, 0.95])
+  #deltas = numpy.array([0.95])
+  #deltas = numpy.arange(0.05,1.,0.05)
+  #rhos = numpy.array([0.05])
   numvects = 10; # Number of vectors to generate
   SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_stdtest.mat'
@@ -56,13 +57,13 @@
   
   d = 50.0;
   sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std1.mat'
@@ -86,13 +87,13 @@
   
   d = 20.0
   sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std2.mat'
@@ -117,13 +118,13 @@
   
   d = 50.0;
   sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std3.mat'
@@ -147,13 +148,13 @@
   
   d = 20.0
   sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std4.mat'
@@ -177,13 +178,13 @@
   
   d = 50.0;
   sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std1nesta.mat'
@@ -207,13 +208,13 @@
   
   d = 20.0
   sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 20.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std2nesta.mat'
@@ -237,13 +238,13 @@
   
   d = 50.0;
   sigma = 2.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std3nesta.mat'
@@ -267,13 +268,13 @@
   
   d = 20.0
   sigma = 10.0
-  deltas = np.arange(0.05,1.,0.05)
-  rhos = np.arange(0.05,1.,0.05)
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
   numvects = 100; # Number of vectors to generate
   SNRdb = 10.;    # This is norm(signal)/norm(noise), so power, not energy
   # Values for lambda
   #lambdas = [0 10.^linspace(-5, 4, 10)];
-  lambdas = np.array([0., 0.0001, 0.01, 1, 100, 10000])
+  lambdas = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
   
   dosavedata = True
   savedataname = 'approx_pt_std4nesta.mat'