nikcleju@0: # -*- coding: utf-8 -*-
nikcleju@0: """
nikcleju@17: Main script for approximate reconstruction tests.
nikcleju@17: Author: Nicolae Cleju
nikcleju@19: 
nikcleju@17: """
nikcleju@17: __author__ = "Nicolae Cleju"
nikcleju@17: __license__ = "GPL"
nikcleju@17: __email__ = "nikcleju@gmail.com"
nikcleju@0: 
nikcleju@0: 
nikcleju@15: import numpy
nikcleju@0: import scipy.io
nikcleju@0: import math
nikcleju@0: import os
nikcleju@0: import time
nikcleju@17: import multiprocessing
nikcleju@17: import sys
nikcleju@0: 
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: proccount = multiprocessing.Value('I', 0) # 'I' = unsigned int, see docs (multiprocessing, array)
nikcleju@15: 
nikcleju@17: # Contains pre-defined simulation parameters
nikcleju@13: import stdparams_approx
nikcleju@17: 
nikcleju@17: # Analysis operator and data generation functions
nikcleju@13: import AnalysisGenerate
nikcleju@0: 
nikcleju@13: # For exceptions
nikcleju@13: import pyCSalgos.BP.l1qec
nikcleju@13: import pyCSalgos.BP.l1qc
nikcleju@13: import pyCSalgos.NESTA.NESTA
nikcleju@21: import pyCSalgos.SL0.EllipseProj
nikcleju@15: 
nikcleju@19: # For plotting with right axes
nikcleju@19: import utils
nikcleju@19: 
nikcleju@17: 
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@0:     currmodule = sys.modules[__name__]
nikcleju@0:     currmodule.proccount = share
nikcleju@15:     currmodule.ntasks = ntasks
nikcleju@21:     currmodule.printLock = printLock
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:   algosN = globalparams['algosN']
nikcleju@15:   algosL = globalparams['algosL']
nikcleju@15:   lambdas = globalparams['lambdas']
nikcleju@15:   
nikcleju@15:   # Process parameters
nikcleju@21:   # noiselevel = norm(noise)/norm(signal) = SNR**(-1/2) = 10**-(SNRdb/20)
nikcleju@21:   noiselevel = 1.0 / (10.0**(SNRdb/20.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@0:           
nikcleju@15:           # Append task params
nikcleju@15:           taskparams.append((algosN,algosL, Omega,y,lambdas,realnoise,M,x0))
nikcleju@0:   
nikcleju@15:   return taskparams
nikcleju@0: 
nikcleju@19: def processResults(params, taskparams, taskresults):
nikcleju@15:   """
nikcleju@15:   Process the raw task results
nikcleju@15:   """
nikcleju@15:   deltas = params['deltas']
nikcleju@15:   rhos = params['rhos']
nikcleju@15:   algosN = params['algosN']
nikcleju@15:   algosL = params['algosL']
nikcleju@15:   lambdas = params['lambdas']
nikcleju@19:   numvects = params['numvects']
nikcleju@0:   
nikcleju@19:   err = dict()
nikcleju@19:   relerr = dict()
nikcleju@19:   mrelerrN = dict()
nikcleju@19:   mrelerrL = dict()
nikcleju@0:   meanmatrix = dict()
nikcleju@0:   elapsed = dict()
nikcleju@19: 
nikcleju@19:   # Prepare storage variables for algorithms non-Lambda
nikcleju@15:   for algo in algosN:
nikcleju@19:     err[algo[1]]     = numpy.zeros(numvects)
nikcleju@19:     relerr[algo[1]]  = numpy.zeros(numvects)
nikcleju@15:     meanmatrix[algo[1]]   = numpy.zeros((rhos.size, deltas.size))
nikcleju@0:     elapsed[algo[1]] = 0
nikcleju@19: 
nikcleju@19:   # Prepare storage variables for algorithms with Lambda    
nikcleju@15:   for algo in algosL:
nikcleju@19:     err[algo[1]]     = numpy.zeros((lambdas.size, numvects))
nikcleju@19:     relerr[algo[1]]  = numpy.zeros((lambdas.size, numvects))
nikcleju@15:     meanmatrix[algo[1]]   = numpy.zeros((lambdas.size, rhos.size, deltas.size))
nikcleju@15:     elapsed[algo[1]] = numpy.zeros(lambdas.size)
nikcleju@0: 
nikcleju@15:   # Process results  
nikcleju@0:   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@19:       algosN,algosL,Omega,y,lambdas,realnoise,M,x0 = taskparams[idx]
nikcleju@19:       xrec,addelapsed = taskresults[idx]
nikcleju@0:       idx = idx+1
nikcleju@19:       
nikcleju@19:       # Optionally compare not with original signal x0, but with solution of 
nikcleju@19:       #  another algorithm (e.g. GAP, NESTA etc)
nikcleju@19:       if 'reference_signal' in params:
nikcleju@19:         xref = xrec[params['reference_signal']]
nikcleju@19:       else:
nikcleju@19:         xref = x0
nikcleju@19: 
nikcleju@19:       # Compute errors and relative errors
nikcleju@19:       for iy in numpy.arange(y.shape[1]):
nikcleju@19:         for algofunc,algoname in algosN: 
nikcleju@19:           err[algoname][iy]    = numpy.linalg.norm(xref[:,iy] - xrec[algoname][:,iy])
nikcleju@19:           relerr[algoname][iy] = err[algoname][iy] / numpy.linalg.norm(xref[:,iy])
nikcleju@19:       for ilbd,lbd in zip(numpy.arange(lambdas.size),lambdas):
nikcleju@19:         for iy in numpy.arange(y.shape[1]):
nikcleju@19:           for algofunc,algoname in algosL:
nikcleju@19:             err[algoname][ilbd,iy]    = numpy.linalg.norm(xref[:,iy] - xrec[algoname][ilbd,:,iy])
nikcleju@19:             relerr[algoname][ilbd,iy] = err[algoname][ilbd,iy] / numpy.linalg.norm(xref[:,iy])
nikcleju@19: 
nikcleju@19:       # Compute average relative errors and prepare matrix to plot
nikcleju@19:       for algofunc,algoname in algosN: 
nikcleju@19:         mrelerrN[algoname] = numpy.mean(relerr[algoname])
nikcleju@19:         meanmatrix[algoname][irho,idelta] = 1 - mrelerrN[algoname]
nikcleju@19:         if meanmatrix[algoname][irho,idelta] < 0 or math.isnan(meanmatrix[algoname][irho,idelta]):
nikcleju@19:           meanmatrix[algoname][irho,idelta] = 0
nikcleju@19:         elapsed[algoname] = elapsed[algoname] + addelapsed[algoname]
nikcleju@19:       for algofunc, algoname in algosL:
nikcleju@15:         for ilbd in numpy.arange(lambdas.size):
nikcleju@19:           mrelerrL[algoname] = numpy.mean(relerr[algoname],1)
nikcleju@19:           meanmatrix[algoname][ilbd,irho,idelta] = 1 - mrelerrL[algoname][ilbd]
nikcleju@19:           if meanmatrix[algoname][ilbd,irho,idelta] < 0 or math.isnan(meanmatrix[algoname][ilbd,irho,idelta]):
nikcleju@19:             meanmatrix[algoname][ilbd,irho,idelta] = 0
nikcleju@19:           elapsed[algoname][ilbd] = elapsed[algoname][ilbd] + addelapsed[algoname][ilbd]
nikcleju@0: 
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:   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['lambdas'] = params['lambdas']
nikcleju@15:   tosave['saveplotbase'] = params['saveplotbase']
nikcleju@15:   tosave['saveplotexts'] = params['saveplotexts']  
nikcleju@15: 
nikcleju@15:   # Save algo names as cell array
nikcleju@15:   obj_arr = numpy.zeros((len(params['algosN']),), dtype=numpy.object)
nikcleju@15:   idx = 0
nikcleju@15:   for algotuple in params['algosN']:
nikcleju@15:     obj_arr[idx] = algotuple[1]
nikcleju@15:     idx = idx+1
nikcleju@15:   tosave['algosNnames'] = obj_arr
nikcleju@15:   
nikcleju@15:   obj_arr = numpy.zeros((len(params['algosL']),), dtype=numpy.object)
nikcleju@15:   idx = 0
nikcleju@15:   for algotuple in params['algosL']:
nikcleju@15:     obj_arr[idx] = algotuple[1]
nikcleju@15:     idx = idx+1
nikcleju@15:   tosave['algosLnames'] = obj_arr
nikcleju@15:   
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['lambdas'] = mdict['lambdas']
nikcleju@15:   params['saveplotbase'] = mdict['saveplotbase'][0]
nikcleju@15:   params['saveplotexts'] = mdict['saveplotexts']
nikcleju@15:   
nikcleju@15:   algonames = mdict['algosNnames'][:,0]
nikcleju@15:   params['algosNnames'] = []
nikcleju@15:   for algoname in algonames:
nikcleju@15:     params['algosNnames'].append(algoname[0])
nikcleju@15: 
nikcleju@15:   algonames = mdict['algosLnames'][:,0]
nikcleju@15:   params['algosLnames'] = []
nikcleju@15:   for algoname in algonames:
nikcleju@15:     params['algosLnames'].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:   algosNnames = params['algosNnames']
nikcleju@15:   algosLnames = params['algosLnames']
nikcleju@15:   lambdas = params['lambdas']
nikcleju@15:   
nikcleju@15:   for algoname in algosNnames:
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:   for algoname in algosLnames:
nikcleju@15:     for ilbd in numpy.arange(lambdas.size):
nikcleju@15:       plt.figure()
nikcleju@15:       plt.imshow(meanmatrix[algoname][ilbd], cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower')
nikcleju@15:       for ext in saveplotexts:
nikcleju@15:         plt.savefig(saveplotbase + algoname + ('_lbd%.0e' % lambdas[ilbd]) + '.' + ext, bbox_inches='tight')
nikcleju@15: 
nikcleju@19: def plotProveEll1(savedataname):
nikcleju@19:   """
nikcleju@19:   Plot ...
nikcleju@19:   The files are saved in the current folder.
nikcleju@19:   """
nikcleju@19:   params, procresults = loadSim(savedataname)
nikcleju@19:   meanmatrix = procresults['meanmatrix']
nikcleju@19:   saveplotbase = params['saveplotbase']  
nikcleju@19:   saveplotexts = params['saveplotexts']
nikcleju@19:   algosNnames = params['algosNnames']
nikcleju@19:   algosLnames = params['algosLnames']
nikcleju@19:   lambdas = params['lambdas']
nikcleju@19:   
nikcleju@19:   toplot = numpy.zeros((len(lambdas),len(algosNnames) + len(algosLnames)))
nikcleju@19:   
nikcleju@19:   idxcol = 0
nikcleju@19:   for algoname in algosNnames:
nikcleju@19:     toplot[:,idxcol] = (1 - meanmatrix[algoname][0,0]) * numpy.ones(len(lambdas))
nikcleju@19:     idxcol = idxcol + 1
nikcleju@19:   for algoname in algosLnames:
nikcleju@19:     for ilbd in numpy.arange(len(lambdas)):
nikcleju@19:       toplot[ilbd,idxcol] = 1 - meanmatrix[algoname][ilbd][0,0]
nikcleju@19:     idxcol = idxcol + 1
nikcleju@19:   
nikcleju@19:   plt.figure()
nikcleju@19:   plt.plot(toplot)
nikcleju@19:   for ext in saveplotexts:
nikcleju@19:     plt.savefig(saveplotbase + '.' + ext, bbox_inches='tight')
nikcleju@19: 
nikcleju@15: #==========================
nikcleju@15: # Main function
nikcleju@15: #==========================
nikcleju@15: def run(params):
nikcleju@17:   """
nikcleju@17:   Run simulation with given parameters
nikcleju@17:   """
nikcleju@17:   
nikcleju@15:   print "This is analysis recovery ABS approximation script by Nic"
nikcleju@17:   print "Running simulation"
nikcleju@15: 
nikcleju@15:   algosN = params['algosN']
nikcleju@15:   algosL = params['algosL']
nikcleju@15:   d = params['d']
nikcleju@15:   sigma = params['sigma']
nikcleju@15:   deltas = params['deltas']
nikcleju@15:   rhos = params['rhos']
nikcleju@15:   lambdas = params['lambdas']
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 algosN],[algotuple[1] for algotuple in algosL]
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@21:     pool.close()
nikcleju@21:     pool.join()
nikcleju@15:   else:
nikcleju@15:     for taskparam in taskparams:
nikcleju@15:       taskresults.append(run_once_tuple(taskparam))
nikcleju@15: 
nikcleju@15:   # Process results
nikcleju@19:   procresults = processResults(params, taskparams, taskresults)
nikcleju@15: 
nikcleju@0:   # Save
nikcleju@15:   saveSim(params, procresults)
nikcleju@15:       
nikcleju@0:   print "Finished."
nikcleju@0:   
nikcleju@0: 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@0:   results = run_once(*t)
nikcleju@17:   
nikcleju@17:   if currmodule.printLock:
nikcleju@17:     currmodule.printLock.acquire()  
nikcleju@17:     
nikcleju@17:     currmodule.proccount.value = currmodule.proccount.value + 1
nikcleju@17:     print "================================"
nikcleju@17:     print "Finished task",currmodule.proccount.value,"of",currmodule.ntasks
nikcleju@17:     print "================================"
nikcleju@17:     
nikcleju@17:     currmodule.printLock.release()
nikcleju@17:     
nikcleju@0:   return results
nikcleju@0: 
nikcleju@0: def run_once(algosN,algosL,Omega,y,lambdas,realnoise,M,x0):
nikcleju@17:   """
nikcleju@17:   Run single task (i.e. task function)
nikcleju@17:   """
nikcleju@0:   
nikcleju@0:   d = Omega.shape[1]  
nikcleju@0:   
nikcleju@0:   xrec = dict()
nikcleju@19: #  err = dict()
nikcleju@19: #  relerr = dict()
nikcleju@0:   elapsed = dict()
nikcleju@0: 
nikcleju@0:   # Prepare storage variables for algorithms non-Lambda
nikcleju@15:   for algo in algosN:
nikcleju@15:     xrec[algo[1]]    = numpy.zeros((d, y.shape[1]))
nikcleju@0:     elapsed[algo[1]] = 0
nikcleju@0:   # Prepare storage variables for algorithms with Lambda    
nikcleju@15:   for algo in algosL:
nikcleju@15:     xrec[algo[1]]    = numpy.zeros((lambdas.size, d, y.shape[1]))
nikcleju@15:     elapsed[algo[1]] = numpy.zeros(lambdas.size)
nikcleju@0:   
nikcleju@0:   # Run algorithms non-Lambda
nikcleju@15:   for iy in numpy.arange(y.shape[1]):
nikcleju@0:     for algofunc,strname in algosN:
nikcleju@15:       epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
nikcleju@0:       try:
nikcleju@0:         timestart = time.time()
nikcleju@0:         xrec[strname][:,iy] = algofunc(y[:,iy],M,Omega,epsilon)
nikcleju@0:         elapsed[strname] = elapsed[strname] + (time.time() - timestart)
nikcleju@0:       except pyCSalgos.BP.l1qec.l1qecInputValueError as e:
nikcleju@0:         print "Caught exception when running algorithm",strname," :",e.message
nikcleju@0:       except pyCSalgos.NESTA.NESTA.NestaError as e:
nikcleju@0:         print "Caught exception when running algorithm",strname," :",e.message
nikcleju@21:       except pyCSalgos.SL0.EllipseProj.EllipseProjDaiError as e:
nikcleju@21:         print "Caught exception when running algorithm",strname," :",e.message
nikcleju@0: 
nikcleju@0:   # Run algorithms with Lambda
nikcleju@15:   for ilbd,lbd in zip(numpy.arange(lambdas.size),lambdas):
nikcleju@15:     for iy in numpy.arange(y.shape[1]):
nikcleju@15:       D = numpy.linalg.pinv(Omega)
nikcleju@15:       #U,S,Vt = numpy.linalg.svd(D)
nikcleju@0:       for algofunc,strname in algosL:
nikcleju@15:         epsilon = 1.1 * numpy.linalg.norm(realnoise[:,iy])
nikcleju@0:         try:
nikcleju@0:           timestart = time.time()
nikcleju@13:           #gamma = algofunc(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd)
nikcleju@19:           xrec[strname][ilbd][:,iy] = algofunc(y[:,iy],M,Omega,epsilon,lbd)
nikcleju@0:           elapsed[strname][ilbd] = elapsed[strname][ilbd] + (time.time() - timestart)
nikcleju@0:         except pyCSalgos.BP.l1qc.l1qcInputValueError as e:
nikcleju@0:           print "Caught exception when running algorithm",strname," :",e.message
nikcleju@21:         except pyCSalgos.SL0.EllipseProj.EllipseProjDaiError as e:
nikcleju@21:           print "Caught exception when running algorithm",strname," :",e.message
nikcleju@21:           
nikcleju@0: 
nikcleju@19:   return xrec, elapsed
nikcleju@19: 
nikcleju@19: def generateFigProveEll1():
nikcleju@19:   """
nikcleju@19:   Generates figure xxx (prove theorem IV.2 in the ell_1 case).
nikcleju@19:   Figures are saved in the current folder.
nikcleju@19:   """
nikcleju@19:   #paramsl1prove['reference_signal'] = nesta[1] # 'NESTA'
nikcleju@19:   run(stdparams_approx.paramsl1prove)
nikcleju@21:   #plotProveEll1(stdparams_approx.paramsl1prove['savedataname'])
nikcleju@21:   utils.replot_ProveEll1(stdparams_approx.paramsl1prove['savedataname'],
nikcleju@21:                      algonames = None, # will read them from mat file
nikcleju@21:                      doshow=False,
nikcleju@21:                      dosave=True,
nikcleju@21:                      saveplotbase=stdparams_approx.paramsl1prove['saveplotbase'],
nikcleju@21:                      saveplotexts=stdparams_approx.paramsl1prove['saveplotexts'])  
nikcleju@0: 
nikcleju@17: def generateFig():
nikcleju@17:   """
nikcleju@17:   Generates figures. 
nikcleju@17:   Figures are saved in the current folder.
nikcleju@17:   """
nikcleju@21:   #stdparams_approx.params1['ncpus'] = 1
nikcleju@21:   run(stdparams_approx.params1)
nikcleju@21:   utils.replot_approx(stdparams_approx.params1['savedataname'],
nikcleju@21:                      algonames = None, # will read them from mat file
nikcleju@21:                      doshow=False,
nikcleju@21:                      dosave=True,
nikcleju@21:                      saveplotbase=stdparams_approx.params1['saveplotbase'],
nikcleju@21:                      saveplotexts=stdparams_approx.params1['saveplotexts'])
nikcleju@21: 
nikcleju@21:   #stdparams_approx.params2['ncpus'] = 1
nikcleju@21:   run(stdparams_approx.params2)
nikcleju@21:   utils.replot_approx(stdparams_approx.params2['savedataname'],
nikcleju@21:                      algonames = None, # will read them from mat file
nikcleju@21:                      doshow=False,
nikcleju@21:                      dosave=True,
nikcleju@21:                      saveplotbase=stdparams_approx.params2['saveplotbase'],
nikcleju@21:                      saveplotexts=stdparams_approx.params2['saveplotexts'])
nikcleju@21:                      
nikcleju@21:   #stdparams_approx.params3['ncpus'] = 1
nikcleju@21:   run(stdparams_approx.params3)
nikcleju@21:   utils.replot_approx(stdparams_approx.params3['savedataname'],
nikcleju@21:                      algonames = None, # will read them from mat file
nikcleju@21:                      doshow=False,
nikcleju@21:                      dosave=True,
nikcleju@21:                      saveplotbase=stdparams_approx.params3['saveplotbase'],
nikcleju@21:                      saveplotexts=stdparams_approx.params3['saveplotexts'])
nikcleju@21: 
nikcleju@21:   #stdparams_approx.params4['ncpus'] = 1
nikcleju@21:   run(stdparams_approx.params4)
nikcleju@21:   utils.replot_approx(stdparams_approx.params4['savedataname'],
nikcleju@21:                      algonames = None, # will read them from mat file
nikcleju@21:                      doshow=False,
nikcleju@21:                      dosave=True,
nikcleju@21:                      saveplotbase=stdparams_approx.params4['saveplotbase'],
nikcleju@21:                      saveplotexts=stdparams_approx.params4['saveplotexts'])                     
nikcleju@21: 
nikcleju@0: # Script main
nikcleju@0: if __name__ == "__main__":
nikcleju@17: 
nikcleju@19:   #stdparams_approx.paramsl1prove['ncpus'] = 1
nikcleju@19:   #generateFigProveEll1()
nikcleju@19: 
nikcleju@19:   #generateFig()
nikcleju@15:   
nikcleju@17:   # Run test parameters
nikcleju@17:   #stdparams_approx.paramstest['ncpus'] = 1
nikcleju@17:   #run(stdparams_approx.paramstest)
nikcleju@17:   #plot(stdparams_approx.paramstest['savedataname'])
nikcleju@21:   #utils.replot_approx(stdparams_approx.paramstest['savedataname'], algonames = None,doshow=False,dosave=True,saveplotbase=stdparams_approx.paramstest['saveplotbase'],saveplotexts=stdparams_approx.paramstest['saveplotexts'])
nikcleju@21:   
nikcleju@21:   #stdparams_approx.params5['ncpus'] = 3
nikcleju@21:   #run(stdparams_approx.params5)
nikcleju@21:   #utils.replot_approx(stdparams_approx.params5['savedataname'],
nikcleju@21:   #                   algonames = None, # will read them from mat file
nikcleju@21:   #                   doshow=False,
nikcleju@21:   #                   dosave=True,
nikcleju@21:   #                   saveplotbase=stdparams_approx.params5['saveplotbase'],
nikcleju@21:   #                   saveplotexts=stdparams_approx.params5['saveplotexts'])  
nikcleju@21: 
nikcleju@21: #  stdparams_approx.params3sl0['ncpus'] = 1
nikcleju@21: #  run(stdparams_approx.params3sl0)
nikcleju@21: #  utils.replot_approx(stdparams_approx.params3sl0['savedataname'],
nikcleju@21: #                     algonames = None, # will read them from mat file
nikcleju@21: #                     doshow=False,
nikcleju@21: #                     dosave=True,
nikcleju@21: #                     saveplotbase=stdparams_approx.params3sl0['saveplotbase'],
nikcleju@21: #                     saveplotexts=stdparams_approx.params3sl0['saveplotexts'])  
nikcleju@21: 
nikcleju@21:   #stdparams_approx.params3['ncpus'] = 1
nikcleju@21:   run(stdparams_approx.params3)
nikcleju@21:   utils.replot_approx(stdparams_approx.params3['savedataname'],
nikcleju@19:                      algonames = None, # will read them from mat file
nikcleju@19:                      doshow=False,
nikcleju@19:                      dosave=True,
nikcleju@21:                      saveplotbase=stdparams_approx.params3['saveplotbase'],
nikcleju@21:                      saveplotexts=stdparams_approx.params3['saveplotexts'])  
nikcleju@21:                      
nikcleju@21:   #stdparams_approx.params4['ncpus'] = 1
nikcleju@21:   run(stdparams_approx.params4)
nikcleju@21:   utils.replot_approx(stdparams_approx.params4['savedataname'],
nikcleju@21:                      algonames = None, # will read them from mat file
nikcleju@21:                      doshow=False,
nikcleju@21:                      dosave=True,
nikcleju@21:                      saveplotbase=stdparams_approx.params4['saveplotbase'],
nikcleju@21:                      saveplotexts=stdparams_approx.params4['saveplotexts'])