changeset 14:f2eb027ed101

test_exact.py working. Added bp_cvxopt(). Commented the code that made the operator Omega more coherent.

--- a/ABSexact.py	Mon Mar 12 17:04:00 2012 +0200
+++ b/ABSexact.py	Fri Mar 16 13:42:31 2012 +0200
@@ -7,11 +7,12 @@
 
 import numpy
 import pyCSalgos.BP.l1eq_pd
+import pyCSalgos.BP.cvxopt_lp
 import pyCSalgos.OMP.omp_QR
 import pyCSalgos.SL0.SL0
 import pyCSalgos.TST.RecommendedTST
 
-def bp(y,M,Omega,x0, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False):
+def bp(y,M,Omega,x0, pdtol=1e-3, pdmaxiter=50, cgtol=1e-8, cgmaxiter=200, verbose=False):
 
   N,n = Omega.shape
   D = numpy.linalg.pinv(Omega)
@@ -22,7 +23,21 @@
   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
 
-  return numpy.dot(D , pyCSalgos.BP.l1eq_pd.l1eq_pd(x0,Atilde,Atilde.T,ytilde, lbtol, mu, cgtol, cgmaxiter, verbose))
+  return numpy.dot(D , pyCSalgos.BP.l1eq_pd.l1eq_pd(x0,Atilde,Atilde.T,ytilde, pdtol, pdmaxiter, cgtol, cgmaxiter, verbose))
+
+def bp_cvxopt(y,M,Omega):
+
+  N,n = Omega.shape
+  D = numpy.linalg.pinv(Omega)
+  U,S,Vt = numpy.linalg.svd(D)
+  Aextra = Vt[-(N-n):,:]
+  
+  # Create aggregate problem
+  Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
+  ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
+
+  return numpy.dot(D , pyCSalgos.BP.cvxopt_lp.cvxopt_lp(ytilde, Atilde))
+
 
 def ompeps(y,M,Omega,epsilon):
   
@@ -38,7 +53,7 @@
   opts = dict()
   opts['stopCrit'] = 'mse'
   opts['stopTol'] = epsilon
-  return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0])
+  return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0])
 
 def ompk(y,M,Omega,k):
   
@@ -53,7 +68,7 @@
   
   opts = dict()
   opts['stopTol'] = k
-  return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0])
+  return numpy.dot(D , pyCSalgos.OMP.omp_QR.greed_omp_qr(ytilde,Atilde,Atilde.shape[1],opts)[0])
 
 def sl0(y,M,Omega, sigma_min, sigma_decrease_factor=0.5, mu_0=2, L=3, true_s=None):
   
@@ -79,5 +94,5 @@
   Atilde = numpy.vstack((numpy.dot(M,D), Aextra))
   ytilde = numpy.concatenate((y,numpy.zeros(N-n)))
   
-  return numpy.dot(D, pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(Atilde, ytilde, nsweep, tol, xinitial, ro))
+  return numpy.dot(D, pyCSalgos.TST.RecommendedTST.RecommendedTST(Atilde, ytilde, nsweep, tol, xinitial, ro))
   
\ No newline at end of file

--- a/algos.py	Mon Mar 12 17:04:00 2012 +0200
+++ b/algos.py	Fri Mar 16 13:42:31 2012 +0200
@@ -39,7 +39,15 @@
   Wrapper for BP algorithm for exact analysis recovery
   Algorithm implementation is l1eq_pd() from l1-magic toolbox
   """
-  return ABSexact.bp(y,M,Omega,numpy.zeros(Omega.shape[0]))
+  return ABSexact.bp(y,M,Omega,numpy.zeros(Omega.shape[0]), pdtol=1e-5, pdmaxiter = 100)
+
+def run_exact_bp_cvxopt(y,M,Omega):
+  """
+  Wrapper for BP algorithm for exact analysis recovery
+  Algorithm implementation is using cvxopt linear programming
+  """
+  return ABSexact.bp_cvxopt(y,M,Omega)
+
 
 def run_exact_ompeps(y,M,Omega):
   """
@@ -156,22 +164,23 @@
 ### Algorithm tuples
 ## Exact recovery
 exact_gap = (run_exact_gap, 'GAP')
-exact_bp = (run_exact_bp, 'ABS-exact: BP')
-exact_ompeps = (run_exact_ompeps, 'ABS-exact: OMP-eps')
-exact_sl0 = (run_exact_sl0, 'ABS-exact: SL0')
-exact_tst = (run_exact_tst, 'ABS-exact: TST')
+exact_bp = (run_exact_bp, 'ABSexact_BP')
+exact_bp_cvxopt = (run_exact_bp_cvxopt, 'ABSexact_BP_cvxopt')
+exact_ompeps = (run_exact_ompeps, 'ABSexact_OMPeps')
+exact_sl0 = (run_exact_sl0, 'ABSexact_SL0')
+exact_tst = (run_exact_tst, 'ABSexact_TST')
 ## Approximate recovery
 # Native
 gap = (run_gap, 'GAP')
 nesta = (run_nesta, 'NESTA')
 # ABS-mixed
-mixed_sl0 = (run_mixed_sl0, 'ABS-mixed: SL0')
-mixed_bp = (run_mixed_bp, 'ABS-mixed: BP')
+mixed_sl0 = (run_mixed_sl0, 'ABSmixed_SL0')
+mixed_bp = (run_mixed_bp, 'ABSmixed_BP')
 # ABS-lambda
-lambda_sl0 = (run_lambda_sl0, 'ABS-lambda: SL0')
-lambda_bp = (run_lambda_bp, 'ABS-lambda: BP')
-lambda_ompeps = (run_lambda_ompeps, 'ABS-lambda: OMP-eps')
-lambda_tst = (run_lambda_tst, 'ABS-lambda: TST')
+lambda_sl0 = (run_lambda_sl0, 'ABSlambda_SL0')
+lambda_bp = (run_lambda_bp, 'ABSlambda_BP')
+lambda_ompeps = (run_lambda_ompeps, 'ABSlambda_OMPeps')
+lambda_tst = (run_lambda_tst, 'ABSlambda_TST')
 
 ##==========================
 ## Algorithm functions

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stdparams_exact.py	Fri Mar 16 13:42:31 2012 +0200
@@ -0,0 +1,147 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Dec 07 14:04:40 2011
+
+@author: ncleju
+"""
+
+import numpy
+from algos import *
+
+#==========================
+# Standard parameters
+#==========================
+# Standard parameters for quick testing
+# Algorithms: GAP, SL0 and BP
+# d=50, sigma = 2, delta and rho only 3 x 3, lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+# Useful for short testing 
+def stdtest():
+  # Define which algorithms to run
+  #algos = exact_gap,exact_sl0,exact_bp,exact_ompeps,exact_tst       # tuple of algorithms
+  algos = exact_bp_cvxopt,       # tuple of algorithms
+  
+  d = 50.0
+  sigma = 1.2
+  deltas = numpy.array([0.05, 0.45, 0.95])
+  rhos = numpy.array([0.05, 0.45, 0.95])
+  #deltas = numpy.array([0.6])
+  #deltas = numpy.arange(0.05,1.,0.05)
+  #rhos = numpy.array([0.05])
+  numvects = 10; # Number of vectors to generate
+  SNRdb = 100.;    # This is norm(signal)/norm(noise), so power, not energy
+  
+  dosavedata = True
+  savedataname = 'exact_pt_stdtest.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'exact_pt_stdtest_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algos,d,sigma,deltas,rhos,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts   
+
+
+# Standard parameters 1
+# All algorithms, 100 vectors
+# d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+def std1():
+  # Define which algorithms to run
+  algos = exact_gap,exact_sl0,exact_bp_cvxopt,exact_ompeps,exact_tst       # tuple of algorithms
+  
+  d = 50.0;
+  sigma = 1.2
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
+  numvects = 100; # Number of vectors to generate
+  SNRdb = 100.;    # This is norm(signal)/norm(noise), so power, not energy
+  
+  dosavedata = True
+  savedataname = 'exact_pt_std1.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'exact_pt_std1_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algos,d,sigma,deltas,rhos,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+          
+         
+# Standard parameters 2
+# All algorithms, 100 vectors
+# d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+# Do save data, do save plots, don't show plots
+def std2():
+  # Define which algorithms to run
+  algos = exact_gap,exact_sl0,exact_bp_cvxopt,exact_ompeps,exact_tst       # tuple of algorithms
+  
+  d = 20.0
+  sigma = 10.0
+  deltas = numpy.arange(0.05,1.,0.05)
+  rhos = numpy.arange(0.05,1.,0.05)
+  numvects = 100; # Number of vectors to generate
+  SNRdb = 100.;    # This is norm(signal)/norm(noise), so power, not energy
+  
+  dosavedata = True
+  savedataname = 'exact_pt_std2.mat'
+  doshowplot = False
+  dosaveplot = True
+  saveplotbase = 'exact_pt_std2_'
+  saveplotexts = ('png','pdf','eps')
+
+  return algos,d,sigma,deltas,rhos,numvects,SNRdb,dosavedata,savedataname,\
+          doshowplot,dosaveplot,saveplotbase,saveplotexts
+  
+  
+#  # Standard parameters 3
+## All algorithms, 100 vectors
+## d=50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+## Do save data, do save plots, don't show plots
+## IDENTICAL with 1 but with 10dB SNR noise
+#def std3():
+#  # Define which algorithms to run
+#  algos = exact_gap,exact_sl0,exact_bp,exact_ompeps,exact_tst       # tuple of algorithms
+#  
+#  d = 50.0;
+#  sigma = 2.0
+#  deltas = numpy.arange(0.05,1.,0.05)
+#  rhos = numpy.arange(0.05,1.,0.05)
+#  numvects = 100; # Number of vectors to generate
+#  SNRdb = 100.;    # This is norm(signal)/norm(noise), so power, not energy
+#  
+#  dosavedata = True
+#  savedataname = 'exact_pt_std3.mat'
+#  doshowplot = False
+#  dosaveplot = True
+#  saveplotbase = 'exact_pt_std3_'
+#  saveplotexts = ('png','pdf','eps')
+#
+#  return algos,d,sigma,deltas,rhos,numvects,SNRdb,dosavedata,savedataname,\
+#          doshowplot,dosaveplot,saveplotbase,saveplotexts
+          
+## Standard parameters 4
+## All algorithms, 100 vectors
+## d=20, sigma = 10, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
+## Do save data, do save plots, don't show plots
+## Identical to 2 but with 10dB SNR noise
+#def std4():
+#  # Define which algorithms to run
+#  algos = exact_gap,exact_sl0,exact_bp,exact_ompeps,exact_tst       # tuple of algorithms
+#  
+#  d = 20.0
+#  sigma = 10.0
+#  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
+#  
+#  dosavedata = True
+#  savedataname = 'exact_pt_std4.mat'
+#  doshowplot = False
+#  dosaveplot = True
+#  saveplotbase = 'exact_pt_std4_'
+#  saveplotexts = ('png','pdf','eps')
+#
+#  return algos,d,sigma,deltas,rhos,numvects,SNRdb,dosavedata,savedataname,\
+#          doshowplot,dosaveplot,saveplotbase,saveplotexts

--- a/test_approx.py	Mon Mar 12 17:04:00 2012 +0200
+++ b/test_approx.py	Fri Mar 16 13:42:31 2012 +0200
@@ -137,7 +137,7 @@
   jobresults = []
   
   if doparallel:
-    pool = multiprocessing.Pool(None,initializer=initProcess,initargs=(currmodule.proccount,currmodule.njobs))
+    pool = multiprocessing.Pool(ncpus,initializer=initProcess,initargs=(currmodule.proccount,currmodule.njobs))
     jobresults = pool.map(run_once_tuple, jobparams)
   else:
     for jobparam in jobparams:
@@ -303,10 +303,10 @@
   # Generate Omega and data based on parameters
   Omega = AnalysisGenerate.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!
-  Snew = np.vstack((np.diag(Sdnew), np.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1]))))
-  Omega = np.dot(U , np.dot(Snew,Vt))
+  #U,S,Vt = np.linalg.svd(Omega);
+  #Sdnew = S * (1+np.arange(S.size)) # Make D coherent, not Omega!
+  #Snew = np.vstack((np.diag(Sdnew), np.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1]))))
+  #Omega = np.dot(U , np.dot(Snew,Vt))
 
   # Generate data  
   x0,y,M,Lambda,realnoise = AnalysisGenerate.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0');

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test_exact.py	Fri Mar 16 13:42:31 2012 +0200
@@ -0,0 +1,337 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sat Nov 05 18:08:40 2011
+
+@author: Nic
+"""
+
+import numpy
+import scipy.io
+import math
+import os
+import time
+
+import multiprocessing
+import sys
+_currmodule = sys.modules[__name__]
+# Lock for printing in a thread-safe way
+_printLock = None
+
+import stdparams_exact
+import AnalysisGenerate
+
+# For exceptions
+import pyCSalgos.BP.l1eq_pd
+import pyCSalgos.NESTA.NESTA
+
+#def printsafe(*t):
+#  """
+#  Thread-safe version of print().
+#  Acquires lock before printing, releases it after.
+#  """  
+#  if _currmodule._printLock:
+#    _currmodule._printLock.acquire()
+#    print t
+#    _currmodule._printLock.release()
+
+
+def _initProcess(share, njobs, printLock):
+    """
+    Pool initializer function (multiprocessing)
+    Needed to pass the shared variable to the worker processes
+    The variables must be global in the module in order to be seen later in run_once_tuple()
+    see http://stackoverflow.com/questions/1675766/how-to-combine-pool-map-with-array-shared-memory-in-python-multiprocessing
+    """
+    currmodule = sys.modules[__name__]
+    currmodule.proccount = share
+    currmodule.njobs = njobs
+    currmodule._printLock = printLock
+          
+#==========================
+# Interface run functions
+#==========================
+def run(std=stdparams_exact.std1,ncpus=None):
+  
+  algos,d,sigma,deltas,rhos,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = std()
+  run_multi(algos, d,sigma,deltas,rhos,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\
+            ncpus=ncpus,\
+            doshowplot=doshowplot,dosaveplot=dosaveplot,saveplotbase=saveplotbase,saveplotexts=saveplotexts)
+
+#==========================
+# Main functions  
+#==========================
+def run_multi(algos, d, sigma, deltas, rhos, numvects, SNRdb, ncpus=None,\
+            doshowplot=False, dosaveplot=False, saveplotbase=None, saveplotexts=None,\
+            dosavedata=False, savedataname=None):
+
+  print "This is analysis recovery ABS approximation script by Nic"
+  print "Running phase transition ( run_multi() )"
+
+  if ncpus is None:
+      print "  Running in parallel with default",multiprocessing.cpu_count(),"threads using \"multiprocessing\" package"
+      if multiprocessing.cpu_count() == 1:
+          doparallel = False
+      else:
+          doparallel = True
+  elif ncpus > 1:
+      print "  Running in parallel with",ncpus,"threads using \"multiprocessing\" package"
+      doparallel = True
+  elif ncpus == 1:
+      print "Running single thread"
+      doparallel = False
+  else:
+      print "Wrong number of threads, exiting"
+      return
+      
+  if dosaveplot or doshowplot:
+    try:
+      import matplotlib
+      if doshowplot or os.name == 'nt':
+        print "Importing matplotlib with default (GUI) backend... ",
+      else:
+        print "Importing matplotlib with \"Cairo\" backend... ",
+        matplotlib.use('Cairo')
+      import matplotlib.pyplot as plt
+      import matplotlib.cm as cm
+      import matplotlib.colors as mcolors
+      print "OK"        
+    except:
+      print "FAIL"
+      print "Importing matplotlib.pyplot failed. No figures at all"
+      print "Try selecting a different backend"
+      doshowplot = False
+      dosaveplot = False
+  
+  # Print summary of parameters
+  print "Parameters:"
+  if doshowplot:
+    print "  Showing figures"
+  else:
+    print "  Not showing figures"
+  if dosaveplot:
+    print "  Saving figures as "+saveplotbase+"* with extensions ",saveplotexts
+  else:
+    print "  Not saving figures"
+  print "  Running algorithms",[algotuple[1] for algotuple in algos]
+  
+  nalgos = len(algos)  
+  
+  meanmatrix = dict()
+  elapsed = dict()
+  for i,algo in zip(numpy.arange(nalgos),algos):
+    meanmatrix[algo[1]]   = numpy.zeros((rhos.size, deltas.size))
+    elapsed[algo[1]] = 0
+  
+  # Prepare parameters
+  jobparams = []
+  print "  (delta, rho) pairs to be run:"
+  for idelta,delta in zip(numpy.arange(deltas.size),deltas):
+    for irho,rho in zip(numpy.arange(rhos.size),rhos):
+      
+      # Generate data and operator
+      Omega,x0,y,M,realnoise = generateData(d,sigma,delta,rho,numvects,SNRdb)
+      
+      #Save the parameters, and run after
+      print "    delta = ",delta," rho = ",rho
+      jobparams.append((algos,Omega,y,M,x0))
+
+  print "End of parameters"  
+
+  _currmodule.njobs = len(jobparams)
+  # Thread-safe variable to store number of finished jobs
+  _currmodule.proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array)
+  
+  
+  # Run
+  jobresults = []
+  
+  if doparallel:
+    _currmodule._printLock = multiprocessing.Lock()
+    pool = multiprocessing.Pool(ncpus,initializer=_initProcess,initargs=(_currmodule.proccount,_currmodule.njobs,_currmodule._printLock))
+    jobresults = pool.map(run_once_tuple, jobparams)
+  else:
+    for jobparam in jobparams:
+      jobresults.append(run_once_tuple(jobparam))
+
+  # Read results
+  idx = 0
+  for idelta,delta in zip(numpy.arange(deltas.size),deltas):
+    for irho,rho in zip(numpy.arange(rhos.size),rhos):
+      mrelerr,addelapsed = jobresults[idx]
+      idx = idx+1
+      for algotuple in algos: 
+        meanmatrix[algotuple[1]][irho,idelta] = mrelerr[algotuple[1]]
+        if meanmatrix[algotuple[1]][irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][irho,idelta]):
+          meanmatrix[algotuple[1]][irho,idelta] = 0
+        elapsed[algotuple[1]] = elapsed[algotuple[1]] + addelapsed[algotuple[1]]
+
+  # Save
+  if dosavedata:
+    tosave = dict()
+    tosave['meanmatrix'] = meanmatrix
+    tosave['elapsed'] = elapsed
+    tosave['d'] = d
+    tosave['sigma'] = sigma
+    tosave['deltas'] = deltas
+    tosave['rhos'] = rhos
+    tosave['numvects'] = numvects
+    tosave['SNRdb'] = SNRdb
+    # Save algo names as cell array
+    obj_arr = numpy.zeros((len(algos),), dtype=numpy.object)
+    idx = 0
+    for algotuple in algos:
+      obj_arr[idx] = algotuple[1]
+      idx = idx+1
+    tosave['algonames'] = obj_arr
+    try:
+      scipy.io.savemat(savedataname, tosave)
+    except:
+      print "Save error"
+  # Show
+  if doshowplot or dosaveplot:
+    for algotuple in algos:
+      algoname = algotuple[1]
+      plt.figure()
+      plt.imshow(meanmatrix[algoname], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower')
+      if dosaveplot:
+        for ext in saveplotexts:
+          plt.savefig(saveplotbase + algoname + '.' + ext, bbox_inches='tight')
+    if doshowplot:
+      plt.show()
+    
+  print "Finished."
+  
+def run_once_tuple(t):
+  results = run_once(*t)
+  
+  if _currmodule._printLock:
+    _currmodule._printLock.acquire()  
+    
+    _currmodule.proccount.value = _currmodule.proccount.value + 1
+    print "================================"
+    print "Finished job",_currmodule.proccount.value,"/",_currmodule.njobs,"jobs remaining",_currmodule.njobs - _currmodule.proccount.value,"/",_currmodule.njobs
+    print "================================"
+
+    _currmodule._printLock.release()
+
+  return results
+
+def run_once(algos,Omega,y,M,x0):
+  
+  d = Omega.shape[1]  
+  
+  nalgos = len(algos)  
+ 
+  xrec = dict()
+  err = dict()
+  relerr = dict()
+  elapsed = dict()
+
+  # Prepare storage variables for algorithms
+  for i,algo in zip(numpy.arange(nalgos),algos):
+    xrec[algo[1]]    = numpy.zeros((d, y.shape[1]))
+    err[algo[1]]     = numpy.zeros(y.shape[1])
+    relerr[algo[1]]  = numpy.zeros(y.shape[1])
+    elapsed[algo[1]] = 0
+  
+  # Run algorithms
+  for iy in numpy.arange(y.shape[1]):
+    for algofunc,strname in algos:
+      try:
+        timestart = time.time()
+        xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega)
+        elapsed[strname] = elapsed[strname] + (time.time() - timestart)
+      except pyCSalgos.BP.l1eq_pd.l1eqNotImplementedError as e:
+        if _currmodule._printLock:
+          _currmodule._printLock.acquire()
+          print "Caught exception when running algorithm",strname," :",e.message
+          _currmodule._printLock.release()
+      err[strname][iy]    = numpy.linalg.norm(x0[:,iy] - xrec[strname][:,iy])
+      relerr[strname][iy] = err[strname][iy] / numpy.linalg.norm(x0[:,iy])
+  for algofunc,strname in algos:
+    if _currmodule._printLock:
+      _currmodule._printLock.acquire()
+      print strname,' : avg relative error = ',numpy.mean(relerr[strname])
+      _currmodule._printLock.release()
+
+  # Prepare results
+  #mrelerr = dict()
+  #for algotuple in algos:
+  #  mrelerr[algotuple[1]] = numpy.mean(relerr[algotuple[1]])
+  #return mrelerr,elapsed
+  
+  # Should return number of reconstructions with error < threshold, not average error
+  exactthr = 1e-6
+  mrelerr = dict()
+  for algotuple in algos:
+    mrelerr[algotuple[1]] = float(numpy.count_nonzero(relerr[algotuple[1]] < exactthr)) / y.shape[1]
+  return mrelerr,elapsed
+
+
+def generateData(d,sigma,delta,rho,numvects,SNRdb):
+
+  # Process parameters
+  noiselevel = 1.0 / (10.0**(SNRdb/10.0));
+  p = round(sigma*d);
+  m = round(delta*d);
+  l = round(d - rho*m);
+  
+  # Generate Omega and data based on parameters
+  Omega = AnalysisGenerate.Generate_Analysis_Operator(d, p);
+  # Optionally make Omega more coherent
+  #U,S,Vt = numpy.linalg.svd(Omega);
+  #Sdnew = S * (1+numpy.arange(S.size)) # Make D coherent, not Omega!
+  #Snew = numpy.vstack((numpy.diag(Sdnew), numpy.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1]))))
+  #Omega = numpy.dot(U , numpy.dot(Snew,Vt))
+
+  # Generate data  
+  x0,y,M,Lambda,realnoise = AnalysisGenerate.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
+#  delta = 0.9
+#  rho = 0.05
+#  numvects = 20; # Number of vectors to generate
+#  SNRdb = 7.;    # This is norm(signal)/norm(noise), so power, not energy
+#  lbd = 10000
+#  
+#  Omega,x0,y,M,realnoise = generateData(d,sigma,delta,rho,numvects,SNRdb)
+#  D = numpy.linalg.pinv(Omega)
+#  U,S,Vt = numpy.linalg.svd(D)
+# 
+#  xrec   = numpy.zeros((d, y.shape[1]))
+#  err    = numpy.zeros((y.shape[1]))
+#  relerr = numpy.zeros((y.shape[1]))  
+# 
+#  for iy in numpy.arange(y.shape[1]):
+#    epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
+#    gamma = run_sl0(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd)
+#    xrec[:,iy] = numpy.dot(D,gamma)
+#    err[iy]    = numpy.linalg.norm(x0[:,iy] - xrec[:,iy])
+#    relerr[iy] = err[iy] / numpy.linalg.norm(x0[:,iy])    
+#  
+#  print "Finished runsingleexampledebug()"
+  
+  
+def testMatlab():
+    mdict = scipy.io.loadmat("E:\\CS\\Ale mele\\Analysis_ExactRec\\temp.mat")
+    algos = stdparams_exact.std1()[0]
+    res = run_once(algos, mdict['Omega'].byteswap().newbyteorder(),mdict['y'],mdict['M'],mdict['x0'])
+  
+def generateFig():
+  run(stdparams_exact.std1)
+  
+# Script main
+if __name__ == "__main__":
+  #import cProfile
+  #cProfile.run('mainrun()', 'profile')    
+  #run_mp(stdparams_exact.stdtest)
+  #runsingleexampledebug()
+  
+  run(stdparams_exact.std1, ncpus=3)
+  #testMatlab()
+  #run(stdparams_exact.stdtest, ncpus=1)