changeset 15:a27cfe83fe12

Changing, changing, trying to get a common framework for batch jobs

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/runbatch.py	Tue Mar 20 17:18:23 2012 +0200
@@ -0,0 +1,148 @@
+# -*- 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
+
+class RunBatch():
+    """
+    Class to run batch
+    """
+    
+    def __init__(self, xs, ys, prerunfunc, paramfunc, runfunc, postrunfunc):
+        self.prerunfunc = prerunfunc
+        self.paramfunc = paramfunc
+        self.runfunc = runfunc
+        self.postrunfunc = postrunfunc
+        
+    def initProcess(self, 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
+     
+    def generateTaskParams(self,globalparams):
+      """
+      Generate a list of task parameters
+      """
+      taskparams = []
+      SNRdb = globalparams['SNRdb'] 
+      sigma = globalparams['sigma'] 
+      d = globalparams['d'] 
+      deltas = globalparams['deltas'] 
+      rhos = globalparams['rhos'] 
+      numvects = globalparams['numvects']
+      algosN = globalparams['algosN']
+      algosL = globalparams['algosL']
+      lambdas = globalparams['lambdas']
+      
+      # Process parameters
+      noiselevel = 1.0 / (10.0**(SNRdb/10.0));
+      
+      for delta in deltas:
+          for rho in rhos:
+              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')
+              
+              # Append task params
+              algoparams = []
+              for algo in algosN:
+                  algoparams.append(algo[0],algo[1],Omega,y,realnoise,M,x0)
+              taskparams.append((algosN,algosL, Omega,y,lambdas,realnoise,M,x0))
+      
+      return taskparams     
+       
+    def run(self, params):
+
+      print "This is RunBatch.run() by Nic"
+   
+      ncpus = params['ncpus']
+      savedataname = params['savedataname']
+      
+      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  
+      
+      # Prepare parameters
+      taskparams = generateTaskParams(params)
+      
+      # Store global variables
+      currmodule.ntasks = len(taskparams)
+      
+      # Run
+      taskresults = []
+      if doparallel:
+        currmodule.printLock = multiprocessing.Lock()
+        pool = multiprocessing.Pool(ncpus,initializer=initProcess,initargs=(currmodule.proccount,currmodule.ntasks,currmodule.printLock))
+        taskresults = pool.map(run_once_tuple, taskparams)
+      else:
+        for taskparam in taskparams:
+          taskresults.append(run_once_tuple(taskparam))
+    
+      # Process results
+      procresults = processResults(params, taskresults)
+    
+      # Save
+      saveSim(params, procresults)
+          
+      print "Finished."
+
+    def run_once_tuple(self,t):
+      results = run_once(*t)
+      import sys
+      currmodule = sys.modules[__name__]  
+      currmodule.proccount.value = currmodule.proccount.value + 1
+      print "================================"
+      print "Finished task",currmodule.proccount.value,"of",currmodule.ntasks
+      print "================================"
+      return results
\ No newline at end of file

--- a/stdparams_approx.py	Fri Mar 16 13:42:31 2012 +0200
+++ b/stdparams_approx.py	Tue Mar 20 17:18:23 2012 +0200
@@ -16,34 +16,25 @@
 # 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
-  algosN = nesta,      # tuple of algorithms not depending on lambda
+paramstest = dict()
+paramstest['algosN'] = nesta,      # tuple of algorithms not depending on lambda
   #algosL = sl0,bp    # tuple of algorithms depending on lambda (our ABS approach)
-  algosL = lambda_sl0,
-  
-  d = 50.0
-  sigma = 2.0
-  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 = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
-  
-  dosavedata = True
-  savedataname = 'approx_pt_stdtest.mat'
-  doshowplot = False
-  dosaveplot = True
-  saveplotbase = 'approx_pt_stdtest_'
-  saveplotexts = ('png','pdf','eps')
-
-  return algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,\
-          doshowplot,dosaveplot,saveplotbase,saveplotexts   
+paramstest['algosL'] = lambda_sl0,
+paramstest['d'] = 50.0
+paramstest['sigma'] = 2.0
+paramstest['deltas'] = numpy.array([0.05, 0.45, 0.95])
+paramstest['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])
+paramstest['numvects'] = 10; # Number of vectors to generate
+paramstest['SNRdb'] = 20.;    # This is norm(signal)/norm(noise), so power, not energy
+# Values for lambda
+#lambdas = [0 10.^linspace(-5, 4, 10)];
+paramstest['lambdas'] = numpy.array([0., 0.0001, 0.01, 1, 100, 10000])
+paramstest['savedataname'] = 'approx_pt_stdtest.mat'
+paramstest['saveplotbase'] = 'approx_pt_stdtest_'
+paramstest['saveplotexts'] = ('png','pdf','eps')
 
 
 # Standard parameters 1

--- a/stdparams_exact.py	Fri Mar 16 13:42:31 2012 +0200
+++ b/stdparams_exact.py	Tue Mar 20 17:18:23 2012 +0200
@@ -8,140 +8,51 @@
 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   
-
+paramstest = dict()
+paramstest['algos'] = exact_gap,exact_sl0,exact_bp,exact_ompeps,exact_tst       # tuple of algorithms
+#paramstest['algos'] = exact_bp_cvxopt,       # tuple of algorithms
+paramstest['d'] = 50.0
+paramstest['sigma'] = 1.2
+paramstest['deltas'] = numpy.array([0.05, 0.45, 0.95])
+paramstest['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])
+paramstest['numvects'] = 10; # Number of vectors to generate
+paramstest['SNRdb'] = 100.;    # This is norm(signal)/norm(noise), so power, not energy
+paramstest['savedataname'] = 'exact_pt_stdtest.mat'
+paramstest['saveplotbase'] = 'exact_pt_stdtest_'
+paramstest['saveplotexts'] = ('png','pdf','eps')
 
 # 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')
+params1 = dict()
+params1['algos'] = exact_gap,exact_sl0,exact_bp_cvxopt,exact_ompeps,exact_tst       # tuple of algorithms
+params1['d'] = 50.0;
+params1['sigma'] = 1.2
+params1['deltas'] = numpy.arange(0.05,1.,0.05)
+params1['rhos'] = numpy.arange(0.05,1.,0.05)
+params1['numvects'] = 100; # Number of vectors to generate
+params1['SNRdb'] = 100.;    # This is norm(signal)/norm(noise), so power, not energy
+params1['savedataname'] = 'exact_pt_std1.mat'
+params1['saveplotbase'] = 'exact_pt_std1_'
+params1['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
+params2 = dict()
+params2['algos'] = exact_gap,exact_sl0,exact_bp_cvxopt,exact_ompeps,exact_tst       # tuple of algorithms
+params2['d'] = 20.0
+params2['sigma'] = 10.0
+params2['deltas'] = numpy.arange(0.05,1.,0.05)
+params2['rhos'] = numpy.arange(0.05,1.,0.05)
+params2['numvects'] = 100; # Number of vectors to generate
+params2['SNRdb'] = 100.;    # This is norm(signal)/norm(noise), so power, not energy
+params2['savedataname'] = 'exact_pt_std2.mat'
+params2['saveplotbase'] = 'exact_pt_std2_'
+params2['saveplotexts'] = ('png','pdf','eps')
\ No newline at end of file

--- a/test_approx.py	Fri Mar 16 13:42:31 2012 +0200
+++ b/test_approx.py	Tue Mar 20 17:18:23 2012 +0200
@@ -5,12 +5,34 @@
 @author: Nic
 """
 
-import numpy as np
+import numpy
 import scipy.io
 import math
 import os
 import time
 
+try:
+  import matplotlib
+  if 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
+except:
+  print "FAIL"
+  print "Importing matplotlib.pyplot failed. No figures at all"
+  print "Try selecting a different backend"
+
+import multiprocessing
+import sys
+currmodule = sys.modules[__name__]
+# Lock for printing in a thread-safe way
+printLock = None
+proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array)
+
 import stdparams_approx
 import AnalysisGenerate
 
@@ -18,136 +40,83 @@
 import pyCSalgos.BP.l1qec
 import pyCSalgos.BP.l1qc
 import pyCSalgos.NESTA.NESTA
-#==========================
-# 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
-#==========================
-def initProcess(share, njobs):
-    import sys
+
+def initProcess(share, ntasks, 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.ntasks = ntasks
+    currmodule._printLock = printLock
+
+def generateTaskParams(globalparams):
+  """
+  Generate a list of task parameters
+  """
+  taskparams = []
+  SNRdb = globalparams['SNRdb'] 
+  sigma = globalparams['sigma'] 
+  d = globalparams['d'] 
+  deltas = globalparams['deltas'] 
+  rhos = globalparams['rhos'] 
+  numvects = globalparams['numvects']
+  algosN = globalparams['algosN']
+  algosL = globalparams['algosL']
+  lambdas = globalparams['lambdas']
+  
+  # Process parameters
+  noiselevel = 1.0 / (10.0**(SNRdb/10.0));
+  
+  for delta in deltas:
+      for rho in rhos:
+          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')
           
-#==========================
-# Interface run functions
-#==========================
-def run_mp(std=stdparams.std2,ncpus=None):
+          # Append task params
+          taskparams.append((algosN,algosL, Omega,y,lambdas,realnoise,M,x0))
   
-  algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = std()
-  run_multi(algosN, algosL, d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\
-            doparallel=True, ncpus=ncpus,\
-            doshowplot=doshowplot,dosaveplot=dosaveplot,saveplotbase=saveplotbase,saveplotexts=saveplotexts)
+  return taskparams
 
-def run(std=stdparams.std2):
-  algosN,algosL,d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = std()
-  run_multi(algosN, algosL, d,sigma,deltas,rhos,lambdas,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\
-            doparallel=False,\
-            doshowplot=doshowplot,dosaveplot=dosaveplot,saveplotbase=saveplotbase,saveplotexts=saveplotexts)  
-#==========================
-# Main functions  
-#==========================
-def run_multi(algosN, algosL, d, sigma, deltas, rhos, lambdas, numvects, SNRdb,
-            doparallel=False, 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() )"
-
-  # Not only for parallel  
-  #if doparallel:
-  import multiprocessing
-  # Shared value holding the number of finished processes
-  # Add it as global of the module
-  import sys
-  currmodule = sys.modules[__name__]
-  currmodule.proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array)
-    
-  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 doparallel:
-    if ncpus is None:
-      print "  Running in parallel with default threads using \"multiprocessing\" package"
-    else:
-      print "  Running in parallel with",ncpus,"threads using \"multiprocessing\" package"
-  else:
-    print "Running single thread"
-  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 algosN],[algotuple[1] for algotuple in algosL]
-  
-  nalgosN = len(algosN)  
-  nalgosL = len(algosL)
+def processResults(params, taskresults):
+  """
+  Process the raw task results
+  """
+  deltas = params['deltas']
+  rhos = params['rhos']
+  algosN = params['algosN']
+  algosL = params['algosL']
+  lambdas = params['lambdas']
   
   meanmatrix = dict()
   elapsed = dict()
-  for i,algo in zip(np.arange(nalgosN),algosN):
-    meanmatrix[algo[1]]   = np.zeros((rhos.size, deltas.size))
+  for algo in algosN:
+    meanmatrix[algo[1]]   = numpy.zeros((rhos.size, deltas.size))
     elapsed[algo[1]] = 0
-  for i,algo in zip(np.arange(nalgosL),algosL):
-    meanmatrix[algo[1]]   = np.zeros((lambdas.size, rhos.size, deltas.size))
-    elapsed[algo[1]] = np.zeros(lambdas.size)
-  
-  # Prepare parameters
-  jobparams = []
-  print "  (delta, rho) pairs to be run:"
-  for idelta,delta in zip(np.arange(deltas.size),deltas):
-    for irho,rho in zip(np.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((algosN,algosL, Omega,y,lambdas,realnoise,M,x0))
-  
-  # Not only for parallel
-  #if doparallel:
-  currmodule.njobs = deltas.size * rhos.size  
-  print "End of parameters"
-  
-  # Run
-  jobresults = []
-  
-  if doparallel:
-    pool = multiprocessing.Pool(ncpus,initializer=initProcess,initargs=(currmodule.proccount,currmodule.njobs))
-    jobresults = pool.map(run_once_tuple, jobparams)
-  else:
-    for jobparam in jobparams:
-      jobresults.append(run_once_tuple(jobparam))
+  for algo in algosL:
+    meanmatrix[algo[1]]   = numpy.zeros((lambdas.size, rhos.size, deltas.size))
+    elapsed[algo[1]] = numpy.zeros(lambdas.size)
 
-  # Read results
+  # Process results  
   idx = 0
-  for idelta,delta in zip(np.arange(deltas.size),deltas):
-    for irho,rho in zip(np.arange(rhos.size),rhos):
-      mrelerrN,mrelerrL,addelapsed = jobresults[idx]
+  for idelta,delta in zip(numpy.arange(deltas.size),deltas):
+    for irho,rho in zip(numpy.arange(rhos.size),rhos):
+      mrelerrN,mrelerrL,addelapsed = taskresults[idx]
       idx = idx+1
       for algotuple in algosN: 
         meanmatrix[algotuple[1]][irho,idelta] = 1 - mrelerrN[algotuple[1]]
@@ -155,55 +124,174 @@
           meanmatrix[algotuple[1]][irho,idelta] = 0
         elapsed[algotuple[1]] = elapsed[algotuple[1]] + addelapsed[algotuple[1]]
       for algotuple in algosL:
-        for ilbd in np.arange(lambdas.size):
+        for ilbd in numpy.arange(lambdas.size):
           meanmatrix[algotuple[1]][ilbd,irho,idelta] = 1 - mrelerrL[algotuple[1]][ilbd]
           if meanmatrix[algotuple[1]][ilbd,irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][ilbd,irho,idelta]):
             meanmatrix[algotuple[1]][ilbd,irho,idelta] = 0
           elapsed[algotuple[1]][ilbd] = elapsed[algotuple[1]][ilbd] + addelapsed[algotuple[1]][ilbd]
 
+        
+  procresults = dict()
+  procresults['meanmatrix'] = meanmatrix
+  procresults['elapsed'] = elapsed
+  return procresults
+
+def saveSim(params, procresults):
+  """
+  Save simulation to mat file
+  """
+  tosave = dict()
+  tosave['meanmatrix'] = procresults['meanmatrix']
+  tosave['elapsed'] = procresults['elapsed']
+  tosave['d'] = params['d']
+  tosave['sigma'] = params['sigma']
+  tosave['deltas'] = params['deltas']
+  tosave['rhos'] = params['rhos']
+  tosave['numvects'] = params['numvects']
+  tosave['SNRdb'] = params['SNRdb']
+  tosave['lambdas'] = params['lambdas']
+  tosave['saveplotbase'] = params['saveplotbase']
+  tosave['saveplotexts'] = params['saveplotexts']  
+
+  # Save algo names as cell array
+  obj_arr = numpy.zeros((len(params['algosN']),), dtype=numpy.object)
+  idx = 0
+  for algotuple in params['algosN']:
+    obj_arr[idx] = algotuple[1]
+    idx = idx+1
+  tosave['algosNnames'] = obj_arr
+  
+  obj_arr = numpy.zeros((len(params['algosL']),), dtype=numpy.object)
+  idx = 0
+  for algotuple in params['algosL']:
+    obj_arr[idx] = algotuple[1]
+    idx = idx+1
+  tosave['algosLnames'] = obj_arr
+  
+  try:
+    scipy.io.savemat(params['savedataname'], tosave)
+  except:
+    print "Save error"
+
+def loadSim(savedataname):
+  """
+  Load simulation from mat file
+  """
+  mdict = scipy.io.loadmat(savedataname)
+  
+  params = dict()
+  procresults = dict()
+  
+  procresults['meanmatrix'] = mdict['meanmatrix'][0,0]
+  procresults['elapsed'] = mdict['elapsed']
+  params['d'] = mdict['d']
+  params['sigma'] = mdict['sigma']
+  params['deltas'] = mdict['deltas']
+  params['rhos'] = mdict['rhos']
+  params['numvects'] = mdict['numvects']
+  params['SNRdb'] = mdict['SNRdb']
+  params['lambdas'] = mdict['lambdas']
+  params['saveplotbase'] = mdict['saveplotbase'][0]
+  params['saveplotexts'] = mdict['saveplotexts']
+  
+  algonames = mdict['algosNnames'][:,0]
+  params['algosNnames'] = []
+  for algoname in algonames:
+    params['algosNnames'].append(algoname[0])
+
+  algonames = mdict['algosLnames'][:,0]
+  params['algosLnames'] = []
+  for algoname in algonames:
+    params['algosLnames'].append(algoname[0])
+  
+  return params, procresults
+
+def plot(savedataname):
+  """
+  Plot results
+  """
+  params, procresults = loadSim(savedataname)
+  meanmatrix = procresults['meanmatrix']
+  saveplotbase = params['saveplotbase']  
+  saveplotexts = params['saveplotexts']
+  algosNnames = params['algosNnames']
+  algosLnames = params['algosLnames']
+  lambdas = params['lambdas']
+  
+  for algoname in algosNnames:
+    plt.figure()
+    plt.imshow(meanmatrix[algoname], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower')
+    for ext in saveplotexts:
+      plt.savefig(saveplotbase + algoname + '.' + ext, bbox_inches='tight')
+  for algoname in algosLnames:
+    for ilbd in numpy.arange(lambdas.size):
+      plt.figure()
+      plt.imshow(meanmatrix[algoname][ilbd], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower')
+      for ext in saveplotexts:
+        plt.savefig(saveplotbase + algoname + ('_lbd%.0e' % lambdas[ilbd]) + '.' + ext, bbox_inches='tight')
+
+#==========================
+# Main function
+#==========================
+def run(params):
+
+  print "This is analysis recovery ABS approximation script by Nic"
+  print "Running phase transition ( run_multi() )"
+
+  algosN = params['algosN']
+  algosL = params['algosL']
+  d = params['d']
+  sigma = params['sigma']
+  deltas = params['deltas']
+  rhos = params['rhos']
+  lambdas = params['lambdas']
+  numvects = params['numvects']
+  SNRdb = params['SNRdb']
+  ncpus = params['ncpus']
+  savedataname = params['savedataname']
+  
+  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  
+  
+  # Print summary of parameters
+  print "Parameters:"
+  print "  Running algorithms",[algotuple[1] for algotuple in algosN],[algotuple[1] for algotuple in algosL]
+  
+  # Prepare parameters
+  taskparams = generateTaskParams(params)
+  
+  # Store global variables
+  currmodule.ntasks = len(taskparams)
+  
+  # Run
+  taskresults = []
+  if doparallel:
+    currmodule.printLock = multiprocessing.Lock()
+    pool = multiprocessing.Pool(ncpus,initializer=initProcess,initargs=(currmodule.proccount,currmodule.ntasks,currmodule.printLock))
+    taskresults = pool.map(run_once_tuple, taskparams)
+  else:
+    for taskparam in taskparams:
+      taskresults.append(run_once_tuple(taskparam))
+
+  # Process results
+  procresults = processResults(params, taskresults)
+
   # 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
-    tosave['lambdas'] = lambdas
-    # Save algo names as cell array
-    obj_arr = np.zeros((len(algosN)+len(algosL),), dtype=np.object)
-    idx = 0
-    for algotuple in algosN+algosL:
-      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 algosN:
-      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')
-    for algotuple in algosL:
-      algoname = algotuple[1]
-      for ilbd in np.arange(lambdas.size):
-        plt.figure()
-        plt.imshow(meanmatrix[algoname][ilbd], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower')
-        if dosaveplot:
-          for ext in saveplotexts:
-            plt.savefig(saveplotbase + algoname + ('_lbd%.0e' % lambdas[ilbd]) + '.' + ext, bbox_inches='tight')
-    if doshowplot:
-      plt.show()
-    
+  saveSim(params, procresults)
+      
   print "Finished."
   
 def run_once_tuple(t):
@@ -212,7 +300,7 @@
   currmodule = sys.modules[__name__]  
   currmodule.proccount.value = currmodule.proccount.value + 1
   print "================================"
-  print "Finished job",currmodule.proccount.value,"of",currmodule.njobs
+  print "Finished task",currmodule.proccount.value,"of",currmodule.ntasks
   print "================================"
   return results
 
@@ -220,32 +308,28 @@
   
   d = Omega.shape[1]  
   
-  nalgosN = len(algosN)  
-  nalgosL = len(algosL)
-  
-  
   xrec = dict()
   err = dict()
   relerr = dict()
   elapsed = dict()
 
   # Prepare storage variables for algorithms non-Lambda
-  for i,algo in zip(np.arange(nalgosN),algosN):
-    xrec[algo[1]]    = np.zeros((d, y.shape[1]))
-    err[algo[1]]     = np.zeros(y.shape[1])
-    relerr[algo[1]]  = np.zeros(y.shape[1])
+  for algo in algosN:
+    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
   # Prepare storage variables for algorithms with Lambda    
-  for i,algo in zip(np.arange(nalgosL),algosL):
-    xrec[algo[1]]    = np.zeros((lambdas.size, d, y.shape[1]))
-    err[algo[1]]     = np.zeros((lambdas.size, y.shape[1]))
-    relerr[algo[1]]  = np.zeros((lambdas.size, y.shape[1]))
-    elapsed[algo[1]] = np.zeros(lambdas.size)
+  for algo in algosL:
+    xrec[algo[1]]    = numpy.zeros((lambdas.size, d, y.shape[1]))
+    err[algo[1]]     = numpy.zeros((lambdas.size, y.shape[1]))
+    relerr[algo[1]]  = numpy.zeros((lambdas.size, y.shape[1]))
+    elapsed[algo[1]] = numpy.zeros(lambdas.size)
   
   # Run algorithms non-Lambda
-  for iy in np.arange(y.shape[1]):
+  for iy in numpy.arange(y.shape[1]):
     for algofunc,strname in algosN:
-      epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
+      epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
       try:
         timestart = time.time()
         xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega,epsilon)
@@ -254,18 +338,18 @@
         print "Caught exception when running algorithm",strname," :",e.message
       except pyCSalgos.NESTA.NESTA.NestaError as e:
         print "Caught exception when running algorithm",strname," :",e.message
-      err[strname][iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][:,iy])
-      relerr[strname][iy] = err[strname][iy] / np.linalg.norm(x0[:,iy])
+      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 algosN:
-    print strname,' : avg relative error = ',np.mean(relerr[strname])  
+    print strname,' : avg relative error = ',numpy.mean(relerr[strname])  
 
   # Run algorithms with Lambda
-  for ilbd,lbd in zip(np.arange(lambdas.size),lambdas):
-    for iy in np.arange(y.shape[1]):
-      D = np.linalg.pinv(Omega)
-      #U,S,Vt = np.linalg.svd(D)
+  for ilbd,lbd in zip(numpy.arange(lambdas.size),lambdas):
+    for iy in numpy.arange(y.shape[1]):
+      D = numpy.linalg.pinv(Omega)
+      #U,S,Vt = numpy.linalg.svd(D)
       for algofunc,strname in algosL:
-        epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
+        epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
         try:
           timestart = time.time()
           #gamma = algofunc(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd)
@@ -273,47 +357,24 @@
           elapsed[strname][ilbd] = elapsed[strname][ilbd] + (time.time() - timestart)
         except pyCSalgos.BP.l1qc.l1qcInputValueError as e:
           print "Caught exception when running algorithm",strname," :",e.message
-        xrec[strname][ilbd,:,iy] = np.dot(D,gamma)
-        err[strname][ilbd,iy]    = np.linalg.norm(x0[:,iy] - xrec[strname][ilbd,:,iy])
-        relerr[strname][ilbd,iy] = err[strname][ilbd,iy] / np.linalg.norm(x0[:,iy])
+        xrec[strname][ilbd,:,iy] = numpy.dot(D,gamma)
+        err[strname][ilbd,iy]    = numpy.linalg.norm(x0[:,iy] - xrec[strname][ilbd,:,iy])
+        relerr[strname][ilbd,iy] = err[strname][ilbd,iy] / numpy.linalg.norm(x0[:,iy])
     print 'Lambda = ',lbd,' :'
     for algofunc,strname in algosL:
-      print '   ',strname,' : avg relative error = ',np.mean(relerr[strname][ilbd,:])
+      print '   ',strname,' : avg relative error = ',numpy.mean(relerr[strname][ilbd,:])
 
   # Prepare results
   mrelerrN = dict()
   for algotuple in algosN:
-    mrelerrN[algotuple[1]] = np.mean(relerr[algotuple[1]])
+    mrelerrN[algotuple[1]] = numpy.mean(relerr[algotuple[1]])
   mrelerrL = dict()
   for algotuple in algosL:
-    mrelerrL[algotuple[1]] = np.mean(relerr[algotuple[1]],1)
+    mrelerrL[algotuple[1]] = numpy.mean(relerr[algotuple[1]],1)
   
   return mrelerrN,mrelerrL,elapsed
 
 
-
-def generateData(d,sigma,delta,rho,numvects,SNRdb):
-
-  # Process parameters
-  noiselevel = 1.0 / (10.0**(SNRdb/20.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 = 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');
-  
-  return Omega,x0,y,M,realnoise
-
-
 def runsingleexampledebug():
   d = 50.0;
   sigma = 2.0
@@ -324,19 +385,19 @@
   lbd = 10000
   
   Omega,x0,y,M,realnoise = generateData(d,sigma,delta,rho,numvects,SNRdb)
-  D = np.linalg.pinv(Omega)
-  U,S,Vt = np.linalg.svd(D)
+  D = numpy.linalg.pinv(Omega)
+  U,S,Vt = numpy.linalg.svd(D)
  
-  xrec   = np.zeros((d, y.shape[1]))
-  err    = np.zeros((y.shape[1]))
-  relerr = np.zeros((y.shape[1]))  
+  xrec   = numpy.zeros((d, y.shape[1]))
+  err    = numpy.zeros((y.shape[1]))
+  relerr = numpy.zeros((y.shape[1]))  
  
-  for iy in np.arange(y.shape[1]):
-    epsilon = 1.1 * np.linalg.norm(realnoise[:,iy])
+  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] = np.dot(D,gamma)
-    err[iy]    = np.linalg.norm(x0[:,iy] - xrec[:,iy])
-    relerr[iy] = err[iy] / np.linalg.norm(x0[:,iy])    
+    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()"
   
@@ -344,5 +405,8 @@
 if __name__ == "__main__":
   #import cProfile
   #cProfile.run('mainrun()', 'profile')    
-  run(stdparams.stdtest)
   #runsingleexampledebug()
+  
+  stdparams_approx.paramstest['ncpus'] = 1
+  run(stdparams_approx.paramstest)
+  plot(stdparams_approx.paramstest['savedataname'])

--- a/test_exact.py	Fri Mar 16 13:42:31 2012 +0200
+++ b/test_exact.py	Tue Mar 20 17:18:23 2012 +0200
@@ -24,16 +24,6 @@
 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):
     """
@@ -141,7 +131,6 @@
   # Thread-safe variable to store number of finished jobs
   _currmodule.proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array)
   
-  
   # Run
   jobresults = []
   
@@ -289,33 +278,6 @@
   
   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")

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test_exact_v2.py	Tue Mar 20 17:18:23 2012 +0200
@@ -0,0 +1,354 @@
+# -*- 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
+
+try:
+  import matplotlib
+  if 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
+except:
+  print "FAIL"
+  print "Importing matplotlib.pyplot failed. No figures at all"
+  print "Try selecting a different backend"
+
+import multiprocessing
+import sys
+currmodule = sys.modules[__name__]
+# Lock for printing in a thread-safe way
+printLock = None
+# Thread-safe variable to store number of finished tasks
+currmodule.proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array)
+
+import stdparams_exact
+import AnalysisGenerate
+
+# For exceptions
+import pyCSalgos.BP.l1eq_pd
+import pyCSalgos.NESTA.NESTA
+
+
+def initProcess(share, ntasks, 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.ntasks = ntasks
+    currmodule._printLock = printLock
+
+
+def generateTaskParams(globalparams):
+  """
+  Generate a list of task parameters
+  """
+  taskparams = []
+  SNRdb = globalparams['SNRdb'] 
+  sigma = globalparams['sigma'] 
+  d = globalparams['d'] 
+  deltas = globalparams['deltas'] 
+  rhos = globalparams['rhos'] 
+  numvects = globalparams['numvects']
+  algos = globalparams['algos']
+  
+  # Process parameters
+  noiselevel = 1.0 / (10.0**(SNRdb/10.0));
+  
+  for delta in deltas:
+      for rho in rhos:
+          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')
+          
+          # Append task params
+          taskparams.append((algos,Omega,y,M,x0))
+  
+  return taskparams
+
+def processResults(params, taskresults):
+  """
+  Process the raw task results
+  """
+  deltas = params['deltas']
+  rhos = params['rhos']
+  algos = params['algos']
+  
+  # Init results
+  meanmatrix = dict()
+  elapsed = dict()
+  for algo in algos:
+    meanmatrix[algo[1]] = numpy.zeros((rhos.size, deltas.size))
+    elapsed[algo[1]] = 0  
+
+  # Process 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 = taskresults[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]]    
+        
+  procresults = dict()
+  procresults['meanmatrix'] = meanmatrix
+  procresults['elapsed'] = elapsed
+  return procresults
+         
+def saveSim(params, procresults):
+  """
+  Save simulation to mat file
+  """
+  #tosaveparams = ['d','sigma','deltas','rhos','numvects','SNRdb']
+  #tosaveprocresults = ['meanmatrix','elapsed']
+
+  tosave = dict()
+  tosave['meanmatrix'] = procresults['meanmatrix']
+  tosave['elapsed'] = procresults['elapsed']
+  tosave['d'] = params['d']
+  tosave['sigma'] = params['sigma']
+  tosave['deltas'] = params['deltas']
+  tosave['rhos'] = params['rhos']
+  tosave['numvects'] = params['numvects']
+  tosave['SNRdb'] = params['SNRdb']
+  tosave['saveplotbase'] = params['saveplotbase']
+  tosave['saveplotexts'] = params['saveplotexts']
+  # Save algo names as cell array
+  obj_arr = numpy.zeros((len(params['algos']),), dtype=numpy.object)
+  idx = 0
+  for algotuple in params['algos']:
+    obj_arr[idx] = algotuple[1]
+    idx = idx+1
+  tosave['algonames'] = obj_arr
+  try:
+    scipy.io.savemat(params['savedataname'], tosave)
+  except:
+    print "Save error"
+    
+def loadSim(savedataname):
+  """
+  Load simulation from mat file
+  """
+  mdict = scipy.io.loadmat(savedataname)
+  
+  params = dict()
+  procresults = dict()
+  
+  procresults['meanmatrix'] = mdict['meanmatrix'][0,0]
+  procresults['elapsed'] = mdict['elapsed']
+  params['d'] = mdict['d']
+  params['sigma'] = mdict['sigma']
+  params['deltas'] = mdict['deltas']
+  params['rhos'] = mdict['rhos']
+  params['numvects'] = mdict['numvects']
+  params['SNRdb'] = mdict['SNRdb']
+  params['saveplotbase'] = mdict['saveplotbase'][0]
+  params['saveplotexts'] = mdict['saveplotexts']
+  
+  algonames = mdict['algonames'][:,0]
+  params['algonames'] = []
+  for algoname in algonames:
+    params['algonames'].append(algoname[0])
+  
+  return params, procresults
+
+def plot(savedataname):
+  """
+  Plot results
+  """
+  params, procresults = loadSim(savedataname)
+  meanmatrix = procresults['meanmatrix']
+  saveplotbase = params['saveplotbase']  
+  saveplotexts = params['saveplotexts']
+  algonames = params['algonames']
+  
+  for algoname in algonames:
+    plt.figure()
+    plt.imshow(meanmatrix[algoname], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower')
+    for ext in saveplotexts:
+      plt.savefig(saveplotbase + algoname + '.' + ext, bbox_inches='tight')
+    #plt.show()         
+      
+#==========================
+# Main function
+#==========================
+def run(params):
+  """
+  Run with given parameters
+  """
+
+  print "This is analysis recovery ABS approximation script by Nic"
+  print "Running phase transition ( run_multi() )"
+
+  algos = params['algos']
+  d = params['d']
+  sigma = params['sigma']
+  deltas = params['deltas']
+  rhos = params['rhos']
+  numvects = params['numvects']
+  SNRdb = params['SNRdb']
+  ncpus = params['ncpus']
+  savedataname = params['savedataname']
+
+  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
+      
+  # Print summary of parameters
+  print "Parameters:"
+  print "  Running algorithms",[algotuple[1] for algotuple in algos]
+  
+  # Prepare parameters
+  taskparams = generateTaskParams(params)
+
+  # Store global variables
+  currmodule.ntasks = len(taskparams)
+ 
+  # Run
+  taskresults = []
+  if doparallel:
+    currmodule.printLock = multiprocessing.Lock()
+    pool = multiprocessing.Pool(ncpus,initializer=initProcess,initargs=(currmodule.proccount,currmodule.ntasks,currmodule.printLock))
+    taskresults = pool.map(run_once_tuple, taskparams)
+  else:
+    for taskparam in taskparams:
+      taskresults.append(run_once_tuple(taskparam))
+
+  # Process results
+  procresults = processResults(params, taskresults)
+
+  # Save
+  saveSim(params, procresults)
+    
+  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 task",currmodule.proccount.value,"/",currmodule.ntasks,"tasks remaining",currmodule.ntasks - currmodule.proccount.value,"/",currmodule.ntasks
+    print "================================"
+
+    currmodule.printLock.release()
+
+  return results
+
+def run_once(algos,Omega,y,M,x0):
+  """
+  Run single task (task function)
+  """
+  
+  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 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()
+  
+  stdparams_exact.paramstest['ncpus'] = 1
+  run(stdparams_exact.paramstest)
+  plot(stdparams_exact.paramstest['savedataname'])