diff test_exact.py @ 16:23e9b536ba71

Renamed test_exact.py and test_approx.py. Fixed in ABSlambda.py multiplication with D. Still not finished with runbatch.py, probably I'll stop it.

--- 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'])