Mercurial > hg > absrec

--- a/ABSlambda.py	Tue Mar 20 17:18:23 2012 +0200
+++ b/ABSlambda.py	Tue Apr 03 10:18:11 2012 +0300
@@ -21,7 +21,7 @@
   aggD = numpy.vstack((aggDupper, lbd * aggDlower))
   aggy = numpy.concatenate((y, numpy.zeros(N-n)))

-  return pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,A_pinv,true_s)
+  return numpy.dot(D, pyCSalgos.SL0.SL0_approx.SL0_approx(aggD,aggy,epsilon,sigma_min,sigma_decrease_factor,mu_0,L,A_pinv,true_s))

 def bp(y,M,Omega,epsilon,lbd, x0, lbtol=1e-3, mu=10, cgtol=1e-8, cgmaxiter=200, verbose=False):

@@ -33,7 +33,7 @@
   aggD = numpy.vstack((aggDupper, lbd * aggDlower))
   aggy = numpy.concatenate((y, numpy.zeros(N-n)))

-  return pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon, lbtol, mu, cgtol, cgmaxiter, verbose)
+  return numpy.dot(D, pyCSalgos.BP.l1qc.l1qc_logbarrier(x0,aggD,aggD.T,aggy,epsilon, lbtol, mu, cgtol, cgmaxiter, verbose))

 def ompeps(y,M,Omega,epsilon,lbd):

@@ -48,7 +48,7 @@
   opts = dict()
   opts['stopCrit'] = 'mse'
   opts['stopTol'] = epsilon**2 / aggy.size
-  return pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0]
+  return numpy.dot(D, pyCSalgos.OMP.omp_QR.greed_omp_qr(aggy,aggD,aggD.shape[1],opts)[0])

 def tst_recom(y,M,Omega,epsilon,lbd, nsweep=300, xinitial=None, ro=None):

@@ -61,5 +61,5 @@
   aggy = numpy.concatenate((y, numpy.zeros(N-n)))

   tol = epsilon / numpy.linalg.norm(aggy)
-  return pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep, tol, xinitial, ro)
+  return numpy.dot(D, pyCSalgos.RecomTST.RecommendedTST.RecommendedTST(aggD, aggy, nsweep, tol, xinitial, ro))

\ No newline at end of file
--- a/runbatch.py	Tue Mar 20 17:18:23 2012 +0200
+++ b/runbatch.py	Tue Apr 03 10:18:11 2012 +0300
@@ -18,6 +18,7 @@
 printLock = None

 import stdparams_exact
+import stdparams_approx
 import AnalysisGenerate

 # For exceptions
@@ -29,11 +30,10 @@
     Class to run batch
     """

-    def __init__(self, xs, ys, prerunfunc, paramfunc, runfunc, postrunfunc):
-        self.prerunfunc = prerunfunc
+    def __init__(self, paramfunc, immedResultsFunc, processResultsFunc):
         self.paramfunc = paramfunc
-        self.runfunc = runfunc
-        self.postrunfunc = postrunfunc
+        self.immedResultsFunc = immedResultsFunc
+        self.processResultsFunc = processResultsFunc

     def initProcess(self, share, njobs, printLock):
         """
@@ -46,56 +46,13 @@
         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):
+    def run(self, globalparams):

       print "This is RunBatch.run() by Nic"

-      ncpus = params['ncpus']
-      savedataname = params['savedataname']
+      ncpus = globalparams['ncpus']
+      savedataname = globalparams['savedataname']

       if ncpus is None:
           print "  Running in parallel with default",multiprocessing.cpu_count(),"threads using \"multiprocessing\" package"
@@ -114,7 +71,8 @@
           return

       # Prepare parameters
-      taskparams = generateTaskParams(params)
+      #taskparams = generateTaskParams(params)
+      taskparams, refparams = self.paramfunc(globalparams)

       # Store global variables
       currmodule.ntasks = len(taskparams)
@@ -124,25 +82,233 @@
       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)
+        taskresults = pool.map(self.run_once_tuple, globalparams,taskparams)
       else:
-        for taskparam in taskparams:
-          taskresults.append(run_once_tuple(taskparam))
+        for taskparam,refparam in zip(taskparams, refparams):
+          #taskresults.append(self.run_once_tuple(globalparams,taskparam))
+          taskresults.append(self.run_once(globalparams,taskparam,refparam))

       # Process results
       procresults = processResults(params, taskresults)
+      #procresults = []
+      #for taskresult in taskresults
+      #    procresults.append(self.processResultsFunc(globalparams, taskresult))

       # Save
-      saveSim(params, procresults)
+      saveSim(globalparams, procresults)

       print "Finished."

-    def run_once_tuple(self,t):
-      results = run_once(*t)
+    def run_once_tuple(self,globalparams,taskparam):
+      results = self.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
+      return results
+
+    def run_once(self, globalparams, algoparams, refparam):
+
+      d = globalparams['d']
+
+      #xrec = dict()
+      #err = dict()
+      #relerr = dict()
+      #elapsed = dict()
+
+      # Prepare storage variables for algorithms non-Lambda
+      #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 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)
+
+      taskresult = []
+      rawresults = None
+      for algoparam in algoparams:
+        try:
+          #timestart = time.time()
+          #xrec[strname][:,iy] = algoparams[0](*algoparams[1:])
+          rawresults = algoparam[0](*algoparam[1:])
+          #elapsed[strname] = elapsed[strname] + (time.time() - timestart)
+        except pyCSalgos.BP.l1qec.l1qecInputValueError as e:
+          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
+          print "Caught error! :",e.message
+        #err[strname][iy]    = numpy.linalg.norm(x0[:,iy] - xrec[strname][:,iy])
+        #relerr[strname][iy] = err[strname][iy] / numpy.linalg.norm(x0[:,iy])
+        #taskresult.append(immediateProcessResults(rawresults))
+        taskresult.append(self.immedResultsFunc(rawresults, refparam))
+
+      return taskresult
+
+
+#      # Run algorithms with Lambda
+#      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 * numpy.linalg.norm(realnoise[:,iy])
+#            try:
+#              timestart = time.time()
+#              #gamma = algofunc(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd)
+#              gamma = algofunc(y[:,iy],M,Omega,epsilon,lbd)
+#              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] = 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 = ',numpy.mean(relerr[strname][ilbd,:])
+
+#      # Prepare results
+#      mrelerrN = dict()
+#      for algotuple in algosN:
+#        mrelerrN[algotuple[1]] = numpy.mean(relerr[algotuple[1]])
+#      mrelerrL = dict()
+#      for algotuple in algosL:
+#        mrelerrL[algotuple[1]] = numpy.mean(relerr[algotuple[1]],1)
+#
+#      return mrelerrN,mrelerrL,elapsed
+
+
+
+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);
+          D = numpy.linalg.pinv(Omega)
+
+          # Generate data
+          x0,y,M,Lambda,realnoise = AnalysisGenerate.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0')
+
+          boxparams = []
+          refparams = []
+          # Prepare algos and params to run
+          for iy in numpy.arange(y.shape[1]):
+            for algofunc,strname in algosN:
+              epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
+              boxparams.append((algofunc,y[:,iy],M,Omega,epsilon))
+              refparams.append(x0[:,iy])
+          for ilbd,lbd in zip(numpy.arange(lambdas.size),lambdas):
+            for iy in numpy.arange(y.shape[1]):
+              for algofunc,strname in algosL:
+                epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
+                boxparams.append((algofunc,y[:,iy],M,Omega,epsilon,lbd))
+                refparams.append(x0[:,iy])
+
+          # algoparams = algo function + parameters it requires
+          taskparams.append(boxparams)
+
+  return taskparams, refparams
+
+
+def immedResults(xrec,x0):
+    if xrec == None:
+        return None
+    err = numpy.linalg.norm(x0 - xrec)
+    relerr = err / numpy.linalg.norm(x0)
+
+    return err,relerr
+
+def processResults(globalparams, taskresults):
+  deltas = globalparams['deltas']
+  rhos = globalparams['rhos']
+  algosN = globalparams['algosN']
+  algosL = globalparams['algosL']
+  lambdas = globalparams['lambdas']
+
+  meanmatrix = dict()
+  elapsed = dict()
+  for algo in algosN:
+    meanmatrix[algo[1]]   = numpy.zeros((rhos.size, deltas.size))
+    elapsed[algo[1]] = 0
+  for algo in algosL:
+    meanmatrix[algo[1]]   = numpy.zeros((lambdas.size, rhos.size, deltas.size))
+    elapsed[algo[1]] = numpy.zeros(lambdas.size)
+
+  # To fix this
+  idx = 0
+  for idelta,delta in zip(numpy.arange(deltas.size),deltas):
+    for irho,rho in zip(numpy.arange(rhos.size),rhos):
+        immedResults = taskresults[idx]
+        for iy in numpy.arange(y.shape[1]):
+          for algofunc,strname in algosN:
+            epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
+            boxparams.append((algofunc,y[:,iy],M,Omega,epsilon))
+            refparams.append(x0[:,iy])
+        for ilbd,lbd in zip(numpy.arange(lambdas.size),lambdas):
+          for iy in numpy.arange(y.shape[1]):
+            for algofunc,strname in algosL:
+              epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
+              boxparams.append((algofunc,y[:,iy],M,Omega,epsilon,lbd))
+              refparams.append(x0[:,iy])
+
+
+  # Process results
+  idx = 0
+  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]
+      mrelerrN,mrelerrL = taskresults[idx]
+      idx = idx+1
+      for algotuple in algosN:
+        meanmatrix[algotuple[1]][irho,idelta] = 1 - mrelerrN[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]]
+      for algotuple in algosL:
+        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
+
+# Script main
+if __name__ == "__main__":
+
+  stdparams_approx.paramstest['ncpus'] = 1
+  rb = RunBatch(generateTaskParams, immedResults, processResults)
+  rb.run(stdparams_approx.paramstest)
+
+  print "Finished"
\ No newline at end of file
--- a/test_exact.py	Tue Mar 20 17:18:23 2012 +0200
+++ b/test_exact.py	Tue Apr 03 10:18:11 2012 +0300
@@ -11,11 +11,28 @@
 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__]
+currmodule = sys.modules[__name__]
 # Lock for printing in a thread-safe way
-_printLock = None
+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
@@ -25,7 +42,7 @@
 import pyCSalgos.NESTA.NESTA


-def _initProcess(share, njobs, printLock):
+def initProcess(share, ntasks, printLock):
     """
     Pool initializer function (multiprocessing)
     Needed to pass the shared variable to the worker processes
@@ -34,29 +51,175 @@
     """
     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']
+  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()
+
 #==========================
-# Interface run functions
+# Main function
 #==========================
-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):
+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:
@@ -73,139 +236,53 @@
       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))
+  taskparams = generateTaskParams(params)

-  print "End of parameters"
+  # 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))

-  _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]]
+  # 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
-    # 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()
+  saveSim(params, procresults)

   print "Finished."

 def run_once_tuple(t):
   results = run_once(*t)

-  if _currmodule._printLock:
-    _currmodule._printLock.acquire()
+  if currmodule.printLock:
+    currmodule.printLock.acquire()

-    _currmodule.proccount.value = _currmodule.proccount.value + 1
+    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 "Finished task",currmodule.proccount.value,"/",currmodule.ntasks,"tasks remaining",currmodule.ntasks - currmodule.proccount.value,"/",currmodule.ntasks
     print "================================"

-    _currmodule._printLock.release()
+    currmodule.printLock.release()

   return results

 def run_once(algos,Omega,y,M,x0):
+  """
+  Run single task (task function)
+  """

   d = Omega.shape[1]

@@ -231,17 +308,17 @@
         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()
+        if currmodule.printLock:
+          currmodule.printLock.acquire()
           print "Caught exception when running algorithm",strname," :",e.message
-          _currmodule._printLock.release()
+          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()
+    if currmodule.printLock:
+      currmodule.printLock.acquire()
       print strname,' : avg relative error = ',numpy.mean(relerr[strname])
-      _currmodule._printLock.release()
+      currmodule.printLock.release()

   # Prepare results
   #mrelerr = dict()
@@ -256,28 +333,6 @@
     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 testMatlab():
     mdict = scipy.io.loadmat("E:\\CS\\Ale mele\\Analysis_ExactRec\\temp.mat")
@@ -294,6 +349,6 @@
   #run_mp(stdparams_exact.stdtest)
   #runsingleexampledebug()

-  run(stdparams_exact.std1, ncpus=3)
-  #testMatlab()
-  #run(stdparams_exact.stdtest, ncpus=1)
+  stdparams_exact.paramstest['ncpus'] = 1
+  run(stdparams_exact.paramstest)
+  plot(stdparams_exact.paramstest['savedataname'])
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test_exact_old.py	Tue Apr 03 10:18:11 2012 +0300
@@ -0,0 +1,299 @@
+# -*- 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 _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 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)
--- a/test_exact_v2.py	Tue Mar 20 17:18:23 2012 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,354 +0,0 @@
-# -*- 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'])