nikcleju@14: # -*- coding: utf-8 -*- nikcleju@14: """ nikcleju@17: Old version, ignore it nikcleju@14: nikcleju@17: Author: Nicolae Cleju nikcleju@14: """ nikcleju@14: nikcleju@14: import numpy nikcleju@14: import scipy.io nikcleju@14: import math nikcleju@14: import os nikcleju@14: import time nikcleju@14: nikcleju@14: import multiprocessing nikcleju@14: import sys nikcleju@14: _currmodule = sys.modules[__name__] nikcleju@14: # Lock for printing in a thread-safe way nikcleju@14: _printLock = None nikcleju@14: nikcleju@14: import stdparams_exact nikcleju@14: import AnalysisGenerate nikcleju@14: nikcleju@14: # For exceptions nikcleju@14: import pyCSalgos.BP.l1eq_pd nikcleju@14: import pyCSalgos.NESTA.NESTA nikcleju@14: nikcleju@14: nikcleju@14: def _initProcess(share, njobs, printLock): nikcleju@14: """ nikcleju@14: Pool initializer function (multiprocessing) nikcleju@14: Needed to pass the shared variable to the worker processes nikcleju@14: The variables must be global in the module in order to be seen later in run_once_tuple() nikcleju@14: see http://stackoverflow.com/questions/1675766/how-to-combine-pool-map-with-array-shared-memory-in-python-multiprocessing nikcleju@14: """ nikcleju@14: currmodule = sys.modules[__name__] nikcleju@14: currmodule.proccount = share nikcleju@14: currmodule.njobs = njobs nikcleju@14: currmodule._printLock = printLock nikcleju@14: nikcleju@14: #========================== nikcleju@14: # Interface run functions nikcleju@14: #========================== nikcleju@14: def run(std=stdparams_exact.std1,ncpus=None): nikcleju@14: nikcleju@14: algos,d,sigma,deltas,rhos,numvects,SNRdb,dosavedata,savedataname,doshowplot,dosaveplot,saveplotbase,saveplotexts = std() nikcleju@14: run_multi(algos, d,sigma,deltas,rhos,numvects,SNRdb,dosavedata=dosavedata,savedataname=savedataname,\ nikcleju@14: ncpus=ncpus,\ nikcleju@14: doshowplot=doshowplot,dosaveplot=dosaveplot,saveplotbase=saveplotbase,saveplotexts=saveplotexts) nikcleju@14: nikcleju@14: #========================== nikcleju@14: # Main functions nikcleju@14: #========================== nikcleju@14: def run_multi(algos, d, sigma, deltas, rhos, numvects, SNRdb, ncpus=None,\ nikcleju@14: doshowplot=False, dosaveplot=False, saveplotbase=None, saveplotexts=None,\ nikcleju@14: dosavedata=False, savedataname=None): nikcleju@14: nikcleju@14: print "This is analysis recovery ABS approximation script by Nic" nikcleju@14: print "Running phase transition ( run_multi() )" nikcleju@14: nikcleju@14: if ncpus is None: nikcleju@14: print " Running in parallel with default",multiprocessing.cpu_count(),"threads using \"multiprocessing\" package" nikcleju@14: if multiprocessing.cpu_count() == 1: nikcleju@14: doparallel = False nikcleju@14: else: nikcleju@14: doparallel = True nikcleju@14: elif ncpus > 1: nikcleju@14: print " Running in parallel with",ncpus,"threads using \"multiprocessing\" package" nikcleju@14: doparallel = True nikcleju@14: elif ncpus == 1: nikcleju@14: print "Running single thread" nikcleju@14: doparallel = False nikcleju@14: else: nikcleju@14: print "Wrong number of threads, exiting" nikcleju@14: return nikcleju@14: nikcleju@14: if dosaveplot or doshowplot: nikcleju@14: try: nikcleju@14: import matplotlib nikcleju@14: if doshowplot or os.name == 'nt': nikcleju@14: print "Importing matplotlib with default (GUI) backend... ", nikcleju@14: else: nikcleju@14: print "Importing matplotlib with \"Cairo\" backend... ", nikcleju@14: matplotlib.use('Cairo') nikcleju@14: import matplotlib.pyplot as plt nikcleju@14: import matplotlib.cm as cm nikcleju@14: import matplotlib.colors as mcolors nikcleju@14: print "OK" nikcleju@14: except: nikcleju@14: print "FAIL" nikcleju@14: print "Importing matplotlib.pyplot failed. No figures at all" nikcleju@14: print "Try selecting a different backend" nikcleju@14: doshowplot = False nikcleju@14: dosaveplot = False nikcleju@14: nikcleju@14: # Print summary of parameters nikcleju@14: print "Parameters:" nikcleju@14: if doshowplot: nikcleju@14: print " Showing figures" nikcleju@14: else: nikcleju@14: print " Not showing figures" nikcleju@14: if dosaveplot: nikcleju@14: print " Saving figures as "+saveplotbase+"* with extensions ",saveplotexts nikcleju@14: else: nikcleju@14: print " Not saving figures" nikcleju@14: print " Running algorithms",[algotuple[1] for algotuple in algos] nikcleju@14: nikcleju@14: nalgos = len(algos) nikcleju@14: nikcleju@14: meanmatrix = dict() nikcleju@14: elapsed = dict() nikcleju@14: for i,algo in zip(numpy.arange(nalgos),algos): nikcleju@14: meanmatrix[algo[1]] = numpy.zeros((rhos.size, deltas.size)) nikcleju@14: elapsed[algo[1]] = 0 nikcleju@14: nikcleju@14: # Prepare parameters nikcleju@14: jobparams = [] nikcleju@14: print " (delta, rho) pairs to be run:" nikcleju@14: for idelta,delta in zip(numpy.arange(deltas.size),deltas): nikcleju@14: for irho,rho in zip(numpy.arange(rhos.size),rhos): nikcleju@14: nikcleju@14: # Generate data and operator nikcleju@14: Omega,x0,y,M,realnoise = generateData(d,sigma,delta,rho,numvects,SNRdb) nikcleju@14: nikcleju@14: #Save the parameters, and run after nikcleju@14: print " delta = ",delta," rho = ",rho nikcleju@14: jobparams.append((algos,Omega,y,M,x0)) nikcleju@14: nikcleju@14: print "End of parameters" nikcleju@14: nikcleju@14: _currmodule.njobs = len(jobparams) nikcleju@14: # Thread-safe variable to store number of finished jobs nikcleju@14: _currmodule.proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array) nikcleju@14: nikcleju@14: # Run nikcleju@14: jobresults = [] nikcleju@14: nikcleju@14: if doparallel: nikcleju@14: _currmodule._printLock = multiprocessing.Lock() nikcleju@14: pool = multiprocessing.Pool(ncpus,initializer=_initProcess,initargs=(_currmodule.proccount,_currmodule.njobs,_currmodule._printLock)) nikcleju@14: jobresults = pool.map(run_once_tuple, jobparams) nikcleju@14: else: nikcleju@14: for jobparam in jobparams: nikcleju@14: jobresults.append(run_once_tuple(jobparam)) nikcleju@14: nikcleju@14: # Read results nikcleju@14: idx = 0 nikcleju@14: for idelta,delta in zip(numpy.arange(deltas.size),deltas): nikcleju@14: for irho,rho in zip(numpy.arange(rhos.size),rhos): nikcleju@14: mrelerr,addelapsed = jobresults[idx] nikcleju@14: idx = idx+1 nikcleju@14: for algotuple in algos: nikcleju@14: meanmatrix[algotuple[1]][irho,idelta] = mrelerr[algotuple[1]] nikcleju@14: if meanmatrix[algotuple[1]][irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][irho,idelta]): nikcleju@14: meanmatrix[algotuple[1]][irho,idelta] = 0 nikcleju@14: elapsed[algotuple[1]] = elapsed[algotuple[1]] + addelapsed[algotuple[1]] nikcleju@14: nikcleju@14: # Save nikcleju@14: if dosavedata: nikcleju@14: tosave = dict() nikcleju@14: tosave['meanmatrix'] = meanmatrix nikcleju@14: tosave['elapsed'] = elapsed nikcleju@14: tosave['d'] = d nikcleju@14: tosave['sigma'] = sigma nikcleju@14: tosave['deltas'] = deltas nikcleju@14: tosave['rhos'] = rhos nikcleju@14: tosave['numvects'] = numvects nikcleju@14: tosave['SNRdb'] = SNRdb nikcleju@14: # Save algo names as cell array nikcleju@14: obj_arr = numpy.zeros((len(algos),), dtype=numpy.object) nikcleju@14: idx = 0 nikcleju@14: for algotuple in algos: nikcleju@14: obj_arr[idx] = algotuple[1] nikcleju@14: idx = idx+1 nikcleju@14: tosave['algonames'] = obj_arr nikcleju@14: try: nikcleju@14: scipy.io.savemat(savedataname, tosave) nikcleju@14: except: nikcleju@14: print "Save error" nikcleju@14: # Show nikcleju@14: if doshowplot or dosaveplot: nikcleju@14: for algotuple in algos: nikcleju@14: algoname = algotuple[1] nikcleju@14: plt.figure() nikcleju@14: plt.imshow(meanmatrix[algoname], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower') nikcleju@14: if dosaveplot: nikcleju@14: for ext in saveplotexts: nikcleju@14: plt.savefig(saveplotbase + algoname + '.' + ext, bbox_inches='tight') nikcleju@14: if doshowplot: nikcleju@14: plt.show() nikcleju@14: nikcleju@14: print "Finished." nikcleju@14: nikcleju@14: def run_once_tuple(t): nikcleju@14: results = run_once(*t) nikcleju@14: nikcleju@14: if _currmodule._printLock: nikcleju@14: _currmodule._printLock.acquire() nikcleju@14: nikcleju@14: _currmodule.proccount.value = _currmodule.proccount.value + 1 nikcleju@14: print "================================" nikcleju@14: print "Finished job",_currmodule.proccount.value,"/",_currmodule.njobs,"jobs remaining",_currmodule.njobs - _currmodule.proccount.value,"/",_currmodule.njobs nikcleju@14: print "================================" nikcleju@14: nikcleju@14: _currmodule._printLock.release() nikcleju@14: nikcleju@14: return results nikcleju@14: nikcleju@14: def run_once(algos,Omega,y,M,x0): nikcleju@14: nikcleju@14: d = Omega.shape[1] nikcleju@14: nikcleju@14: nalgos = len(algos) nikcleju@14: nikcleju@14: xrec = dict() nikcleju@14: err = dict() nikcleju@14: relerr = dict() nikcleju@14: elapsed = dict() nikcleju@14: nikcleju@14: # Prepare storage variables for algorithms nikcleju@14: for i,algo in zip(numpy.arange(nalgos),algos): nikcleju@14: xrec[algo[1]] = numpy.zeros((d, y.shape[1])) nikcleju@14: err[algo[1]] = numpy.zeros(y.shape[1]) nikcleju@14: relerr[algo[1]] = numpy.zeros(y.shape[1]) nikcleju@14: elapsed[algo[1]] = 0 nikcleju@14: nikcleju@14: # Run algorithms nikcleju@14: for iy in numpy.arange(y.shape[1]): nikcleju@14: for algofunc,strname in algos: nikcleju@14: try: nikcleju@14: timestart = time.time() nikcleju@14: xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega) nikcleju@14: elapsed[strname] = elapsed[strname] + (time.time() - timestart) nikcleju@14: except pyCSalgos.BP.l1eq_pd.l1eqNotImplementedError as e: nikcleju@14: if _currmodule._printLock: nikcleju@14: _currmodule._printLock.acquire() nikcleju@14: print "Caught exception when running algorithm",strname," :",e.message nikcleju@14: _currmodule._printLock.release() nikcleju@14: err[strname][iy] = numpy.linalg.norm(x0[:,iy] - xrec[strname][:,iy]) nikcleju@14: relerr[strname][iy] = err[strname][iy] / numpy.linalg.norm(x0[:,iy]) nikcleju@14: for algofunc,strname in algos: nikcleju@14: if _currmodule._printLock: nikcleju@14: _currmodule._printLock.acquire() nikcleju@14: print strname,' : avg relative error = ',numpy.mean(relerr[strname]) nikcleju@14: _currmodule._printLock.release() nikcleju@14: nikcleju@14: # Prepare results nikcleju@14: #mrelerr = dict() nikcleju@14: #for algotuple in algos: nikcleju@14: # mrelerr[algotuple[1]] = numpy.mean(relerr[algotuple[1]]) nikcleju@14: #return mrelerr,elapsed nikcleju@14: nikcleju@14: # Should return number of reconstructions with error < threshold, not average error nikcleju@14: exactthr = 1e-6 nikcleju@14: mrelerr = dict() nikcleju@14: for algotuple in algos: nikcleju@14: mrelerr[algotuple[1]] = float(numpy.count_nonzero(relerr[algotuple[1]] < exactthr)) / y.shape[1] nikcleju@14: return mrelerr,elapsed nikcleju@14: nikcleju@14: nikcleju@14: def generateData(d,sigma,delta,rho,numvects,SNRdb): nikcleju@14: nikcleju@14: # Process parameters nikcleju@14: noiselevel = 1.0 / (10.0**(SNRdb/10.0)); nikcleju@14: p = round(sigma*d); nikcleju@14: m = round(delta*d); nikcleju@14: l = round(d - rho*m); nikcleju@14: nikcleju@14: # Generate Omega and data based on parameters nikcleju@14: Omega = AnalysisGenerate.Generate_Analysis_Operator(d, p); nikcleju@14: # Optionally make Omega more coherent nikcleju@14: #U,S,Vt = numpy.linalg.svd(Omega); nikcleju@14: #Sdnew = S * (1+numpy.arange(S.size)) # Make D coherent, not Omega! nikcleju@14: #Snew = numpy.vstack((numpy.diag(Sdnew), numpy.zeros((Omega.shape[0] - Omega.shape[1], Omega.shape[1])))) nikcleju@14: #Omega = numpy.dot(U , numpy.dot(Snew,Vt)) nikcleju@14: nikcleju@14: # Generate data nikcleju@14: x0,y,M,Lambda,realnoise = AnalysisGenerate.Generate_Data_Known_Omega(Omega, d,p,m,l,noiselevel, numvects,'l0'); nikcleju@14: nikcleju@14: return Omega,x0,y,M,realnoise nikcleju@14: nikcleju@14: nikcleju@14: def testMatlab(): nikcleju@14: mdict = scipy.io.loadmat("E:\\CS\\Ale mele\\Analysis_ExactRec\\temp.mat") nikcleju@14: algos = stdparams_exact.std1()[0] nikcleju@14: res = run_once(algos, mdict['Omega'].byteswap().newbyteorder(),mdict['y'],mdict['M'],mdict['x0']) nikcleju@14: nikcleju@14: def generateFig(): nikcleju@14: run(stdparams_exact.std1) nikcleju@14: nikcleju@14: # Script main nikcleju@14: if __name__ == "__main__": nikcleju@14: #import cProfile nikcleju@14: #cProfile.run('mainrun()', 'profile') nikcleju@14: #run_mp(stdparams_exact.stdtest) nikcleju@14: #runsingleexampledebug() nikcleju@14: nikcleju@14: run(stdparams_exact.std1, ncpus=3) nikcleju@14: #testMatlab() nikcleju@14: #run(stdparams_exact.stdtest, ncpus=1)