comparison test_exact_v2.py @ 15:a27cfe83fe12

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

+# -*- 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'])