nikcleju@15: # -*- coding: utf-8 -*-
nikcleju@15: """
nikcleju@17: Main script for exact reconstruction tests.
nikcleju@17: Author: Nicolae Cleju
nikcleju@17: """
nikcleju@17: __author__ = "Nicolae Cleju"
nikcleju@17: __license__ = "GPL"
nikcleju@17: __email__ = "nikcleju@gmail.com"
nikcleju@15: 
nikcleju@15: 
nikcleju@15: import numpy
nikcleju@15: import scipy.io
nikcleju@15: import math
nikcleju@15: import os
nikcleju@15: import time
nikcleju@17: import multiprocessing
nikcleju@17: import sys
nikcleju@15: 
nikcleju@17: # Try to do smart importing of matplotlib
nikcleju@15: try:
nikcleju@15:   import matplotlib
nikcleju@15:   if os.name == 'nt':
nikcleju@15:     print "Importing matplotlib with default (GUI) backend... "
nikcleju@15:   else:
nikcleju@15:     print "Importing matplotlib with \"Cairo\" backend... "
nikcleju@15:     matplotlib.use('Cairo')
nikcleju@15:   import matplotlib.pyplot as plt
nikcleju@15:   import matplotlib.cm as cm
nikcleju@15:   import matplotlib.colors as mcolors
nikcleju@15: except:
nikcleju@15:   print "FAIL"
nikcleju@15:   print "Importing matplotlib.pyplot failed. No figures at all"
nikcleju@15:   print "Try selecting a different backend"
nikcleju@15: 
nikcleju@15: currmodule = sys.modules[__name__]
nikcleju@17: printLock = None  # Lock for printing in a thread-safe way
nikcleju@17:  # Thread-safe variable to store number of finished tasks
nikcleju@15: currmodule.proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array)
nikcleju@15: 
nikcleju@17: # Contains pre-defined simulation parameters
nikcleju@15: import stdparams_exact
nikcleju@17: 
nikcleju@17: # Analysis operator and data generation functions
nikcleju@15: import AnalysisGenerate
nikcleju@15: 
nikcleju@15: # For exceptions
nikcleju@15: import pyCSalgos.BP.l1eq_pd
nikcleju@15: import pyCSalgos.NESTA.NESTA
nikcleju@15: 
nikcleju@19: # For plotting with right axes
nikcleju@19: import utils
nikcleju@19: 
nikcleju@15: 
nikcleju@17: 
nikcleju@15: def initProcess(share, ntasks, printLock):
nikcleju@15:     """
nikcleju@15:     Pool initializer function (multiprocessing)
nikcleju@15:     Needed to pass the shared variable to the worker processes
nikcleju@15:     The variables must be global in the module in order to be seen later in run_once_tuple()
nikcleju@15:     see http://stackoverflow.com/questions/1675766/how-to-combine-pool-map-with-array-shared-memory-in-python-multiprocessing
nikcleju@15:     """
nikcleju@15:     currmodule = sys.modules[__name__]
nikcleju@15:     currmodule.proccount = share
nikcleju@15:     currmodule.ntasks = ntasks
nikcleju@15:     currmodule._printLock = printLock
nikcleju@15: 
nikcleju@15: 
nikcleju@15: def generateTaskParams(globalparams):
nikcleju@15:   """
nikcleju@17:   Generate a list of task parameters (for parallel running)
nikcleju@15:   """
nikcleju@15:   taskparams = []
nikcleju@15:   SNRdb = globalparams['SNRdb'] 
nikcleju@15:   sigma = globalparams['sigma'] 
nikcleju@15:   d = globalparams['d'] 
nikcleju@15:   deltas = globalparams['deltas'] 
nikcleju@15:   rhos = globalparams['rhos'] 
nikcleju@15:   numvects = globalparams['numvects']
nikcleju@15:   algos = globalparams['algos']
nikcleju@15:   
nikcleju@15:   # Process parameters
nikcleju@15:   noiselevel = 1.0 / (10.0**(SNRdb/10.0));
nikcleju@15:   
nikcleju@15:   for delta in deltas:
nikcleju@15:       for rho in rhos:
nikcleju@15:           p = round(sigma*d);
nikcleju@15:           m = round(delta*d);
nikcleju@15:           l = round(d - rho*m);
nikcleju@15:   
nikcleju@15:           # Generate Omega and data based on parameters
nikcleju@15:           Omega = AnalysisGenerate.Generate_Analysis_Operator(d, p);
nikcleju@15:           # Optionally make Omega more coherent
nikcleju@15:           #U,S,Vt = numpy.linalg.svd(Omega);
nikcleju@15:           #Sdnew = S * (1+numpy.arange(S.size)) # Make D coherent, not Omega!
nikcleju@15:           #Snew = numpy.vstack((numpy.diag(Sdnew), numpy.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1]))))
nikcleju@15:           #Omega = numpy.dot(U , numpy.dot(Snew,Vt))
nikcleju@15:         
nikcleju@15:           # Generate data
nikcleju@15:           x0,y,M,Lambda,realnoise = AnalysisGenerate.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0')
nikcleju@15:           
nikcleju@15:           # Append task params
nikcleju@15:           taskparams.append((algos,Omega,y,M,x0))
nikcleju@15:   
nikcleju@15:   return taskparams
nikcleju@15: 
nikcleju@15: def processResults(params, taskresults):
nikcleju@15:   """
nikcleju@15:   Process the raw task results
nikcleju@15:   """
nikcleju@15:   deltas = params['deltas']
nikcleju@15:   rhos = params['rhos']
nikcleju@15:   algos = params['algos']
nikcleju@15:   
nikcleju@15:   # Init results
nikcleju@15:   meanmatrix = dict()
nikcleju@15:   elapsed = dict()
nikcleju@15:   for algo in algos:
nikcleju@15:     meanmatrix[algo[1]] = numpy.zeros((rhos.size, deltas.size))
nikcleju@15:     elapsed[algo[1]] = 0  
nikcleju@15: 
nikcleju@15:   # Process results  
nikcleju@15:   idx = 0
nikcleju@15:   for idelta,delta in zip(numpy.arange(deltas.size),deltas):
nikcleju@15:     for irho,rho in zip(numpy.arange(rhos.size),rhos):
nikcleju@15:       mrelerr,addelapsed = taskresults[idx]
nikcleju@15:       idx = idx+1
nikcleju@15:       for algotuple in algos: 
nikcleju@15:         meanmatrix[algotuple[1]][irho,idelta] = mrelerr[algotuple[1]]
nikcleju@15:         if meanmatrix[algotuple[1]][irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][irho,idelta]):
nikcleju@15:           meanmatrix[algotuple[1]][irho,idelta] = 0
nikcleju@15:         elapsed[algotuple[1]] = elapsed[algotuple[1]] + addelapsed[algotuple[1]]    
nikcleju@15:         
nikcleju@15:   procresults = dict()
nikcleju@15:   procresults['meanmatrix'] = meanmatrix
nikcleju@15:   procresults['elapsed'] = elapsed
nikcleju@15:   return procresults
nikcleju@15:          
nikcleju@15: def saveSim(params, procresults):
nikcleju@15:   """
nikcleju@15:   Save simulation to mat file
nikcleju@15:   """
nikcleju@15:   #tosaveparams = ['d','sigma','deltas','rhos','numvects','SNRdb']
nikcleju@15:   #tosaveprocresults = ['meanmatrix','elapsed']
nikcleju@15: 
nikcleju@15:   tosave = dict()
nikcleju@15:   tosave['meanmatrix'] = procresults['meanmatrix']
nikcleju@15:   tosave['elapsed'] = procresults['elapsed']
nikcleju@15:   tosave['d'] = params['d']
nikcleju@15:   tosave['sigma'] = params['sigma']
nikcleju@15:   tosave['deltas'] = params['deltas']
nikcleju@15:   tosave['rhos'] = params['rhos']
nikcleju@15:   tosave['numvects'] = params['numvects']
nikcleju@15:   tosave['SNRdb'] = params['SNRdb']
nikcleju@15:   tosave['saveplotbase'] = params['saveplotbase']
nikcleju@15:   tosave['saveplotexts'] = params['saveplotexts']
nikcleju@15:   # Save algo names as cell array
nikcleju@15:   obj_arr = numpy.zeros((len(params['algos']),), dtype=numpy.object)
nikcleju@15:   idx = 0
nikcleju@15:   for algotuple in params['algos']:
nikcleju@15:     obj_arr[idx] = algotuple[1]
nikcleju@15:     idx = idx+1
nikcleju@15:   tosave['algonames'] = obj_arr
nikcleju@15:   try:
nikcleju@15:     scipy.io.savemat(params['savedataname'], tosave)
nikcleju@15:   except:
nikcleju@15:     print "Save error"
nikcleju@15:     
nikcleju@15: def loadSim(savedataname):
nikcleju@15:   """
nikcleju@15:   Load simulation from mat file
nikcleju@15:   """
nikcleju@15:   mdict = scipy.io.loadmat(savedataname)
nikcleju@15:   
nikcleju@15:   params = dict()
nikcleju@15:   procresults = dict()
nikcleju@15:   
nikcleju@15:   procresults['meanmatrix'] = mdict['meanmatrix'][0,0]
nikcleju@15:   procresults['elapsed'] = mdict['elapsed']
nikcleju@15:   params['d'] = mdict['d']
nikcleju@15:   params['sigma'] = mdict['sigma']
nikcleju@15:   params['deltas'] = mdict['deltas']
nikcleju@15:   params['rhos'] = mdict['rhos']
nikcleju@15:   params['numvects'] = mdict['numvects']
nikcleju@15:   params['SNRdb'] = mdict['SNRdb']
nikcleju@15:   params['saveplotbase'] = mdict['saveplotbase'][0]
nikcleju@15:   params['saveplotexts'] = mdict['saveplotexts']
nikcleju@15:   
nikcleju@15:   algonames = mdict['algonames'][:,0]
nikcleju@15:   params['algonames'] = []
nikcleju@15:   for algoname in algonames:
nikcleju@15:     params['algonames'].append(algoname[0])
nikcleju@15:   
nikcleju@15:   return params, procresults
nikcleju@15: 
nikcleju@15: def plot(savedataname):
nikcleju@15:   """
nikcleju@17:   Plot results from a mat file.
nikcleju@17:   The files are saved in the current folder.
nikcleju@15:   """
nikcleju@15:   params, procresults = loadSim(savedataname)
nikcleju@15:   meanmatrix = procresults['meanmatrix']
nikcleju@15:   saveplotbase = params['saveplotbase']  
nikcleju@15:   saveplotexts = params['saveplotexts']
nikcleju@15:   algonames = params['algonames']
nikcleju@15:   
nikcleju@15:   for algoname in algonames:
nikcleju@15:     plt.figure()
nikcleju@15:     plt.imshow(meanmatrix[algoname], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower')
nikcleju@15:     for ext in saveplotexts:
nikcleju@15:       plt.savefig(saveplotbase + algoname + '.' + ext, bbox_inches='tight')
nikcleju@15:     #plt.show()         
nikcleju@15:       
nikcleju@15: #==========================
nikcleju@15: # Main function
nikcleju@15: #==========================
nikcleju@15: def run(params):
nikcleju@15:   """
nikcleju@17:   Run simulation with given parameters
nikcleju@15:   """
nikcleju@15: 
nikcleju@17:   print "This is analysis recovery ABS exact script by Nic"
nikcleju@18:   print "Running simulation"
nikcleju@15: 
nikcleju@15:   algos = params['algos']
nikcleju@15:   d = params['d']
nikcleju@15:   sigma = params['sigma']
nikcleju@15:   deltas = params['deltas']
nikcleju@15:   rhos = params['rhos']
nikcleju@15:   numvects = params['numvects']
nikcleju@15:   SNRdb = params['SNRdb']
nikcleju@18:   if 'ncpus' in params:
nikcleju@18:     ncpus = params['ncpus']
nikcleju@18:   else:
nikcleju@18:     ncpus = None
nikcleju@15:   savedataname = params['savedataname']
nikcleju@15: 
nikcleju@15:   if ncpus is None:
nikcleju@15:       print "  Running in parallel with default",multiprocessing.cpu_count(),"threads using \"multiprocessing\" package"
nikcleju@15:       if multiprocessing.cpu_count() == 1:
nikcleju@15:           doparallel = False
nikcleju@15:       else:
nikcleju@15:           doparallel = True
nikcleju@15:   elif ncpus > 1:
nikcleju@15:       print "  Running in parallel with",ncpus,"threads using \"multiprocessing\" package"
nikcleju@15:       doparallel = True
nikcleju@15:   elif ncpus == 1:
nikcleju@15:       print "Running single thread"
nikcleju@15:       doparallel = False
nikcleju@15:   else:
nikcleju@15:       print "Wrong number of threads, exiting"
nikcleju@15:       return
nikcleju@15:       
nikcleju@15:   # Print summary of parameters
nikcleju@15:   print "Parameters:"
nikcleju@15:   print "  Running algorithms",[algotuple[1] for algotuple in algos]
nikcleju@15:   
nikcleju@15:   # Prepare parameters
nikcleju@18:   print "Generating task parameters..."
nikcleju@15:   taskparams = generateTaskParams(params)
nikcleju@15: 
nikcleju@15:   # Store global variables
nikcleju@15:   currmodule.ntasks = len(taskparams)
nikcleju@15:  
nikcleju@15:   # Run
nikcleju@18:   print "Running..."
nikcleju@15:   taskresults = []
nikcleju@15:   if doparallel:
nikcleju@15:     currmodule.printLock = multiprocessing.Lock()
nikcleju@15:     pool = multiprocessing.Pool(ncpus,initializer=initProcess,initargs=(currmodule.proccount,currmodule.ntasks,currmodule.printLock))
nikcleju@15:     taskresults = pool.map(run_once_tuple, taskparams)
nikcleju@15:   else:
nikcleju@15:     for taskparam in taskparams:
nikcleju@15:       taskresults.append(run_once_tuple(taskparam))
nikcleju@15: 
nikcleju@15:   # Process results
nikcleju@15:   procresults = processResults(params, taskresults)
nikcleju@15: 
nikcleju@15:   # Save
nikcleju@15:   saveSim(params, procresults)
nikcleju@15:     
nikcleju@15:   print "Finished."
nikcleju@15:   
nikcleju@15: def run_once_tuple(t):
nikcleju@17:   """
nikcleju@17:   Wrapper for run_once() that explodes the tuple argument t and shows
nikcleju@17:   the number of finished / remaining tasks
nikcleju@17:   """
nikcleju@17:   
nikcleju@17:   # Call run_once() here
nikcleju@15:   results = run_once(*t)
nikcleju@15:   
nikcleju@15:   if currmodule.printLock:
nikcleju@15:     currmodule.printLock.acquire()  
nikcleju@15:     
nikcleju@15:     currmodule.proccount.value = currmodule.proccount.value + 1
nikcleju@15:     print "================================"
nikcleju@15:     print "Finished task",currmodule.proccount.value,"/",currmodule.ntasks,"tasks remaining",currmodule.ntasks - currmodule.proccount.value,"/",currmodule.ntasks
nikcleju@15:     print "================================"
nikcleju@15: 
nikcleju@15:     currmodule.printLock.release()
nikcleju@15: 
nikcleju@15:   return results
nikcleju@15: 
nikcleju@15: def run_once(algos,Omega,y,M,x0):
nikcleju@15:   """
nikcleju@17:   Run single task (i.e. task function)
nikcleju@15:   """
nikcleju@15:   
nikcleju@15:   d = Omega.shape[1]  
nikcleju@15:   
nikcleju@15:   nalgos = len(algos)  
nikcleju@15:  
nikcleju@15:   xrec = dict()
nikcleju@15:   err = dict()
nikcleju@15:   relerr = dict()
nikcleju@15:   elapsed = dict()
nikcleju@15: 
nikcleju@15:   # Prepare storage variables for algorithms
nikcleju@15:   for i,algo in zip(numpy.arange(nalgos),algos):
nikcleju@15:     xrec[algo[1]]    = numpy.zeros((d, y.shape[1]))
nikcleju@15:     err[algo[1]]     = numpy.zeros(y.shape[1])
nikcleju@15:     relerr[algo[1]]  = numpy.zeros(y.shape[1])
nikcleju@15:     elapsed[algo[1]] = 0
nikcleju@15:   
nikcleju@15:   # Run algorithms
nikcleju@15:   for iy in numpy.arange(y.shape[1]):
nikcleju@15:     for algofunc,strname in algos:
nikcleju@15:       try:
nikcleju@15:         timestart = time.time()
nikcleju@15:         xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega)
nikcleju@15:         elapsed[strname] = elapsed[strname] + (time.time() - timestart)
nikcleju@15:       except pyCSalgos.BP.l1eq_pd.l1eqNotImplementedError as e:
nikcleju@15:         if currmodule.printLock:
nikcleju@15:           currmodule.printLock.acquire()
nikcleju@15:           print "Caught exception when running algorithm",strname," :",e.message
nikcleju@15:           currmodule.printLock.release()
nikcleju@18:       except ValueError as e:
nikcleju@18:         if currmodule.printLock:
nikcleju@18:           currmodule.printLock.acquire()
nikcleju@18:           print "Caught ValueError exception when running algorithm",strname," :",e.message
nikcleju@18:           currmodule.printLock.release()
nikcleju@15:       err[strname][iy]    = numpy.linalg.norm(x0[:,iy] - xrec[strname][:,iy])
nikcleju@15:       relerr[strname][iy] = err[strname][iy] / numpy.linalg.norm(x0[:,iy])
nikcleju@15:   for algofunc,strname in algos:
nikcleju@15:     if currmodule.printLock:
nikcleju@15:       currmodule.printLock.acquire()
nikcleju@15:       print strname,' : avg relative error = ',numpy.mean(relerr[strname])
nikcleju@15:       currmodule.printLock.release()
nikcleju@15: 
nikcleju@15:   # Prepare results
nikcleju@15:   #mrelerr = dict()
nikcleju@15:   #for algotuple in algos:
nikcleju@15:   #  mrelerr[algotuple[1]] = numpy.mean(relerr[algotuple[1]])
nikcleju@15:   #return mrelerr,elapsed
nikcleju@15:   
nikcleju@15:   exactthr = 1e-6
nikcleju@15:   mrelerr = dict()
nikcleju@15:   for algotuple in algos:
nikcleju@15:     mrelerr[algotuple[1]] = float(numpy.count_nonzero(relerr[algotuple[1]] < exactthr)) / y.shape[1]
nikcleju@15:   return mrelerr,elapsed
nikcleju@15: 
nikcleju@15:   
nikcleju@15: def testMatlab():
nikcleju@17:     """
nikcleju@17:     For debugging only.
nikcleju@17:     Load parameters from a mat file saved by Matlab.
nikcleju@17:     """
nikcleju@15:     mdict = scipy.io.loadmat("E:\\CS\\Ale mele\\Analysis_ExactRec\\temp.mat")
nikcleju@15:     algos = stdparams_exact.std1()[0]
nikcleju@15:     res = run_once(algos, mdict['Omega'].byteswap().newbyteorder(),mdict['y'],mdict['M'],mdict['x0'])
nikcleju@15:   
nikcleju@15: def generateFig():
nikcleju@17:   """
nikcleju@17:   Generates figures from paper "Analysis-based sparse reconstruction with synthesis-based solvers".
nikcleju@17:   The figures are saved in the current folder.
nikcleju@17:   """
nikcleju@20:   run(stdparams_exact.params1)
nikcleju@19:   #plot(stdparams_exact.params1['savedataname'])
nikcleju@19:   utils.replot_exact(stdparams_exact.params1['savedataname'],
nikcleju@19:                      algonames = None, # will read them from mat file
nikcleju@19:                      doshow=False,
nikcleju@19:                      dosave=True,
nikcleju@19:                      saveplotbase=stdparams_exact.params1['saveplotbase'],
nikcleju@19:                      saveplotexts=stdparams_exact.params1['saveplotexts'])
nikcleju@15:   
nikcleju@15: # Script main
nikcleju@15: if __name__ == "__main__":
nikcleju@15:   
nikcleju@17:   # Set the number of cpus for paraller running (or comment to leave default = max)
nikcleju@17:   #stdparams_exact.paramstest['ncpus'] = 1
nikcleju@18: 
nikcleju@18:   # Run test parameters
nikcleju@18:   #stdparams_exact.paramstest['ncpus'] = 1
nikcleju@18:   #run(stdparams_exact.paramstest)
nikcleju@18:   #plot(stdparams_exact.paramstest['savedataname'])
nikcleju@18: 
nikcleju@18:   #stdparams_exact.params1['ncpus'] = 1
nikcleju@17:   generateFig()
nikcleju@17: