Mercurial > hg > absrec
changeset 19:2837cfeaf353
Fixed plot functions in utils.py
Started working for test_approx.py
| author | Nic Cleju <nikcleju@gmail.com> |
|---|---|
| date | Thu, 05 Apr 2012 13:59:22 +0300 |
| parents | 4a967f4f18a0 |
| children | eccc7a5b9ee3 |
| files | stdparams_approx.py test_approx.py test_exact.py utils.py |
| diffstat | 4 files changed, 209 insertions(+), 70 deletions(-) [+] |
line wrap: on
line diff
--- a/stdparams_approx.py Thu Apr 05 11:01:22 2012 +0300 +++ b/stdparams_approx.py Thu Apr 05 13:59:22 2012 +0300 @@ -34,6 +34,22 @@ paramstest['saveplotexts'] = ('png','pdf','eps') +# Test parameters +paramsl1prove = dict() +paramsl1prove['algosN'] = nesta, # tuple of algorithms not depending on lambda +paramsl1prove['algosL'] = lambda_bp, # tuple of algorithms depending on lambda (ABS-lambda) +paramsl1prove['d'] = 50.0 +paramsl1prove['sigma'] = 2.0 +paramsl1prove['deltas'] = numpy.array([0.8]) # m = 0.8*50 = 40 +paramsl1prove['rhos'] = numpy.array([0.15]) # l = d - rho*m = 50 - 0.15*40 = 44 +paramsl1prove['numvects'] = 100; # Number of vectors to generate +paramsl1prove['SNRdb'] = 40.; # This is norm(signal)/norm(noise), so power, not energy +paramsl1prove['lambdas'] = numpy.array([0.00001, 0.0001, 0.001, 0.01, 0.1, 1, 10, 100, 1000, 10000]) +paramsl1prove['savedataname'] = 'prove_approx_ell1.mat' +paramsl1prove['saveplotbase'] = 'prove_approx_ell1' +paramsl1prove['saveplotexts'] = ('png','pdf','eps') + + # Standard parameters 1 # All algorithms, 100 vectors # d = 50, sigma = 2, delta and rho full resolution (0.05 step), lambdas = 0, 1e-4, 1e-2, 1, 100, 10000
--- a/test_approx.py Thu Apr 05 11:01:22 2012 +0300 +++ b/test_approx.py Thu Apr 05 13:59:22 2012 +0300 @@ -2,6 +2,7 @@ """ Main script for approximate reconstruction tests. Author: Nicolae Cleju + """ __author__ = "Nicolae Cleju" __license__ = "GPL" @@ -48,6 +49,9 @@ import pyCSalgos.BP.l1qc import pyCSalgos.NESTA.NESTA +# For plotting with right axes +import utils + def initProcess(share, ntasks, printLock): @@ -102,7 +106,7 @@ return taskparams -def processResults(params, taskresults): +def processResults(params, taskparams, taskresults): """ Process the raw task results """ @@ -111,13 +115,26 @@ algosN = params['algosN'] algosL = params['algosL'] lambdas = params['lambdas'] + numvects = params['numvects'] + err = dict() + relerr = dict() + mrelerrN = dict() + mrelerrL = dict() meanmatrix = dict() elapsed = dict() + + # Prepare storage variables for algorithms non-Lambda for algo in algosN: + err[algo[1]] = numpy.zeros(numvects) + relerr[algo[1]] = numpy.zeros(numvects) meanmatrix[algo[1]] = numpy.zeros((rhos.size, deltas.size)) elapsed[algo[1]] = 0 + + # Prepare storage variables for algorithms with Lambda for algo in algosL: + err[algo[1]] = numpy.zeros((lambdas.size, numvects)) + relerr[algo[1]] = numpy.zeros((lambdas.size, numvects)) meanmatrix[algo[1]] = numpy.zeros((lambdas.size, rhos.size, deltas.size)) elapsed[algo[1]] = numpy.zeros(lambdas.size) @@ -125,21 +142,43 @@ idx = 0 for idelta,delta in zip(numpy.arange(deltas.size),deltas): for irho,rho in zip(numpy.arange(rhos.size),rhos): - mrelerrN,mrelerrL,addelapsed = taskresults[idx] + algosN,algosL,Omega,y,lambdas,realnoise,M,x0 = taskparams[idx] + xrec,addelapsed = taskresults[idx] idx = idx+1 - for algotuple in algosN: - meanmatrix[algotuple[1]][irho,idelta] = 1 - mrelerrN[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]] - for algotuple in algosL: + + # Optionally compare not with original signal x0, but with solution of + # another algorithm (e.g. GAP, NESTA etc) + if 'reference_signal' in params: + xref = xrec[params['reference_signal']] + else: + xref = x0 + + # Compute errors and relative errors + for iy in numpy.arange(y.shape[1]): + for algofunc,algoname in algosN: + err[algoname][iy] = numpy.linalg.norm(xref[:,iy] - xrec[algoname][:,iy]) + relerr[algoname][iy] = err[algoname][iy] / numpy.linalg.norm(xref[:,iy]) + for ilbd,lbd in zip(numpy.arange(lambdas.size),lambdas): + for iy in numpy.arange(y.shape[1]): + for algofunc,algoname in algosL: + err[algoname][ilbd,iy] = numpy.linalg.norm(xref[:,iy] - xrec[algoname][ilbd,:,iy]) + relerr[algoname][ilbd,iy] = err[algoname][ilbd,iy] / numpy.linalg.norm(xref[:,iy]) + + # Compute average relative errors and prepare matrix to plot + for algofunc,algoname in algosN: + mrelerrN[algoname] = numpy.mean(relerr[algoname]) + meanmatrix[algoname][irho,idelta] = 1 - mrelerrN[algoname] + if meanmatrix[algoname][irho,idelta] < 0 or math.isnan(meanmatrix[algoname][irho,idelta]): + meanmatrix[algoname][irho,idelta] = 0 + elapsed[algoname] = elapsed[algoname] + addelapsed[algoname] + for algofunc, algoname in algosL: for ilbd in numpy.arange(lambdas.size): - meanmatrix[algotuple[1]][ilbd,irho,idelta] = 1 - mrelerrL[algotuple[1]][ilbd] - if meanmatrix[algotuple[1]][ilbd,irho,idelta] < 0 or math.isnan(meanmatrix[algotuple[1]][ilbd,irho,idelta]): - meanmatrix[algotuple[1]][ilbd,irho,idelta] = 0 - elapsed[algotuple[1]][ilbd] = elapsed[algotuple[1]][ilbd] + addelapsed[algotuple[1]][ilbd] + mrelerrL[algoname] = numpy.mean(relerr[algoname],1) + meanmatrix[algoname][ilbd,irho,idelta] = 1 - mrelerrL[algoname][ilbd] + if meanmatrix[algoname][ilbd,irho,idelta] < 0 or math.isnan(meanmatrix[algoname][ilbd,irho,idelta]): + meanmatrix[algoname][ilbd,irho,idelta] = 0 + elapsed[algoname][ilbd] = elapsed[algoname][ilbd] + addelapsed[algoname][ilbd] - procresults = dict() procresults['meanmatrix'] = meanmatrix procresults['elapsed'] = elapsed @@ -240,6 +279,35 @@ for ext in saveplotexts: plt.savefig(saveplotbase + algoname + ('_lbd%.0e' % lambdas[ilbd]) + '.' + ext, bbox_inches='tight') +def plotProveEll1(savedataname): + """ + Plot ... + The files are saved in the current folder. + """ + params, procresults = loadSim(savedataname) + meanmatrix = procresults['meanmatrix'] + saveplotbase = params['saveplotbase'] + saveplotexts = params['saveplotexts'] + algosNnames = params['algosNnames'] + algosLnames = params['algosLnames'] + lambdas = params['lambdas'] + + toplot = numpy.zeros((len(lambdas),len(algosNnames) + len(algosLnames))) + + idxcol = 0 + for algoname in algosNnames: + toplot[:,idxcol] = (1 - meanmatrix[algoname][0,0]) * numpy.ones(len(lambdas)) + idxcol = idxcol + 1 + for algoname in algosLnames: + for ilbd in numpy.arange(len(lambdas)): + toplot[ilbd,idxcol] = 1 - meanmatrix[algoname][ilbd][0,0] + idxcol = idxcol + 1 + + plt.figure() + plt.plot(toplot) + for ext in saveplotexts: + plt.savefig(saveplotbase + '.' + ext, bbox_inches='tight') + #========================== # Main function #========================== @@ -305,7 +373,7 @@ taskresults.append(run_once_tuple(taskparam)) # Process results - procresults = processResults(params, taskresults) + procresults = processResults(params, taskparams, taskresults) # Save saveSim(params, procresults) @@ -341,21 +409,17 @@ d = Omega.shape[1] xrec = dict() - err = dict() - relerr = dict() +# err = dict() +# relerr = dict() elapsed = dict() # Prepare storage variables for algorithms non-Lambda for algo in algosN: 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 # Prepare storage variables for algorithms with Lambda for algo in algosL: xrec[algo[1]] = numpy.zeros((lambdas.size, d, y.shape[1])) - err[algo[1]] = numpy.zeros((lambdas.size, y.shape[1])) - relerr[algo[1]] = numpy.zeros((lambdas.size, y.shape[1])) elapsed[algo[1]] = numpy.zeros(lambdas.size) # Run algorithms non-Lambda @@ -370,10 +434,6 @@ print "Caught exception when running algorithm",strname," :",e.message except pyCSalgos.NESTA.NESTA.NestaError as e: print "Caught exception when running algorithm",strname," :",e.message - 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 algosN: - print strname,' : avg relative error = ',numpy.mean(relerr[strname]) # Run algorithms with Lambda for ilbd,lbd in zip(numpy.arange(lambdas.size),lambdas): @@ -385,26 +445,22 @@ try: timestart = time.time() #gamma = algofunc(y[:,iy],M,Omega,D,U,S,Vt,epsilon,lbd) - gamma = algofunc(y[:,iy],M,Omega,epsilon,lbd) + xrec[strname][ilbd][:,iy] = algofunc(y[:,iy],M,Omega,epsilon,lbd) elapsed[strname][ilbd] = elapsed[strname][ilbd] + (time.time() - timestart) except pyCSalgos.BP.l1qc.l1qcInputValueError as e: print "Caught exception when running algorithm",strname," :",e.message - xrec[strname][ilbd,:,iy] = numpy.dot(D,gamma) - err[strname][ilbd,iy] = numpy.linalg.norm(x0[:,iy] - xrec[strname][ilbd,:,iy]) - relerr[strname][ilbd,iy] = err[strname][ilbd,iy] / numpy.linalg.norm(x0[:,iy]) - print 'Lambda = ',lbd,' :' - for algofunc,strname in algosL: - print ' ',strname,' : avg relative error = ',numpy.mean(relerr[strname][ilbd,:]) - # Prepare results - mrelerrN = dict() - for algotuple in algosN: - mrelerrN[algotuple[1]] = numpy.mean(relerr[algotuple[1]]) - mrelerrL = dict() - for algotuple in algosL: - mrelerrL[algotuple[1]] = numpy.mean(relerr[algotuple[1]],1) - - return mrelerrN,mrelerrL,elapsed + return xrec, elapsed + +def generateFigProveEll1(): + """ + Generates figure xxx (prove theorem IV.2 in the ell_1 case). + Figures are saved in the current folder. + """ + #paramsl1prove['reference_signal'] = nesta[1] # 'NESTA' + run(stdparams_approx.paramsl1prove) + plotProveEll1(stdparams_approx.paramsl1prove['savedataname']) + def generateFig(): """ @@ -417,9 +473,18 @@ # Script main if __name__ == "__main__": - generateFig() + #stdparams_approx.paramsl1prove['ncpus'] = 1 + #generateFigProveEll1() + + #generateFig() # Run test parameters #stdparams_approx.paramstest['ncpus'] = 1 #run(stdparams_approx.paramstest) #plot(stdparams_approx.paramstest['savedataname']) + utils.replot_approx(stdparams_approx.paramstest['savedataname'], + algonames = None, # will read them from mat file + doshow=False, + dosave=True, + saveplotbase=stdparams_approx.paramstest['saveplotbase'], + saveplotexts=stdparams_approx.paramstest['saveplotexts'])
--- a/test_exact.py Thu Apr 05 11:01:22 2012 +0300 +++ b/test_exact.py Thu Apr 05 13:59:22 2012 +0300 @@ -47,6 +47,9 @@ import pyCSalgos.BP.l1eq_pd import pyCSalgos.NESTA.NESTA +# For plotting with right axes +import utils + def initProcess(share, ntasks, printLock): @@ -371,8 +374,14 @@ Generates figures from paper "Analysis-based sparse reconstruction with synthesis-based solvers". The figures are saved in the current folder. """ - run(stdparams_exact.params1) - plot(stdparams_exact.params1['savedataname']) + #run(stdparams_exact.params1) + #plot(stdparams_exact.params1['savedataname']) + utils.replot_exact(stdparams_exact.params1['savedataname'], + algonames = None, # will read them from mat file + doshow=False, + dosave=True, + saveplotbase=stdparams_exact.params1['saveplotbase'], + saveplotexts=stdparams_exact.params1['saveplotexts']) # Script main if __name__ == "__main__":
--- a/utils.py Thu Apr 05 11:01:22 2012 +0300 +++ b/utils.py Thu Apr 05 13:59:22 2012 +0300 @@ -14,25 +14,58 @@ # Sample call #utils.loadshowmatrices_multipixels('H:\\CS\\Python\\Results\\pt_std1\\approx_pt_std1.mat', dosave=True, saveplotbase='approx_pt_std1_',saveplotexts=('png','eps','pdf')) -def loadshowmatrices_multipixels(filename, algonames = None, doshow=True, dosave=False, saveplotbase=None, saveplotexts=None): +def replot_exact(filename, algonames = None, doshow=True, dosave=False, saveplotbase=None, saveplotexts=None): + """ + Replot exact recovery results from mat file, with better axis ticks and + other custom tweaked options. + """ + + mdict = scipy.io.loadmat(filename) + + if algonames == None: + if 'algonames' in mdict: + algonames = mdict['algonames'] + else: + print "No algonames given, and couldn't find them in mat file." + print "Exiting." + return + + loadshowmatrices_multipixels(filename, algonames, [], [], algonames, [], doshow, dosave, saveplotbase, saveplotexts) + +def replot_approx(filename, algonames = None, doshow=True, dosave=False, saveplotbase=None, saveplotexts=None): + """ + Replot exact recovery results from mat file, with better axis ticks and + other custom tweaked options. + """ + + mdict = scipy.io.loadmat(filename) + + if algonames == None: + if 'algosNnames' in mdict and 'algosLnames' in mdict: + algonames = numpy.vstack((mdict['algosNnames'], mdict['algosLnames'])) + else: + print "No algonames given, and couldn't find them in mat file." + print "Exiting." + return + + if dosave: + lambdas = mdict['lambdas'] + threshs = [(0.85,2,0),(0.8,2,0.4),(0.5,2,1)] + withticks = ['GAP'] + withnoaxes = [algoname[0][0] for algoname in algonames if algoname not in withticks] + #withnoaxes.remove('GAP') + loadshowmatrices_multipixels(filename, algonames, lambdas, threshs, withticks, withnoaxes, doshow, dosave, saveplotbase, saveplotexts) + +def loadshowmatrices_multipixels(filename, algonames, lambdas, threshs = [], withticks = [], withnoaxes = [], doshow=True, dosave=False, saveplotbase=None, saveplotexts=None): if dosave and (saveplotbase is None or saveplotexts is None): print('Error: please specify name and extensions for saving') raise Exception('Name or extensions for saving not specified') mdict = scipy.io.loadmat(filename) - - if dosave: - lambdas = mdict['lambdas'] - - if algonames == None: - algonames = mdict['algonames'] - -# thresh = 0.90 - N = 10 -# border = 2 -# bordercolor = 0 # black - + + N = 10 # one data box = NxN + for algonameobj in algonames: algoname = algonameobj[0][0] print algoname @@ -44,10 +77,14 @@ for i in numpy.arange(rows): for j in numpy.arange(cols): bigmatrix[i*N:i*N+N,j*N:j*N+N] = mdict['meanmatrix'][algoname][0,0][i,j] - bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,10,0.95,2,0) - #bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,10,0.9, 2,0.2) - bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,10,0.8, 2,0.4) - bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,10,0.5, 2,1) + + for thrval,width,color in threshs: + bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,N,thrval,width,color) + + #bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,10,0.95,2,0) + #bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,10,0.8, 2,0.4) + #bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0],bigmatrix,10,0.5, 2,1) + # # Mark 95% border # if mdict['meanmatrix'][algoname][0,0][i,j] > thresh: # # Top border @@ -64,11 +101,11 @@ # bigmatrix[i*N:i*N+N,j*N+N-border:j*N+N] = bordercolor plt.figure() - #plt.imshow(mdict['meanmatrix'][algoname][0,0], cmap=cm.gray, interpolation='nearest',origin='lower') plt.imshow(bigmatrix, cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower') - if algoname == 'GAP': + + if algoname in withticks: int_setticks() - else: + if algoname in withnoaxes: plt.gca().get_xaxis().set_visible(False) plt.gca().get_yaxis().set_visible(False) @@ -87,10 +124,13 @@ for i in numpy.arange(rows): for j in numpy.arange(cols): bigmatrix[i*N:i*N+N,j*N:j*N+N] = matrix[i,j] - bigmatrix = int_drawseparation(matrix,bigmatrix,10,0.95,2,0) - #bigmatrix = int_drawseparation(matrix,bigmatrix,10,0.9, 2,0.2) - bigmatrix = int_drawseparation(matrix,bigmatrix,10,0.8, 2,0.4) - bigmatrix = int_drawseparation(matrix,bigmatrix,10,0.5, 2,1) + + for thrval,width,color in threshs: + bigmatrix = int_drawseparation(mdict['meanmatrix'][algoname][0,0][ilbd],bigmatrix,N,thrval,width,color) + + #bigmatrix = int_drawseparation(matrix,bigmatrix,10,0.95,2,0) + #bigmatrix = int_drawseparation(matrix,bigmatrix,10,0.8, 2,0.4) + #bigmatrix = int_drawseparation(matrix,bigmatrix,10,0.5, 2,1) # # Mark 95% border # if matrix[i,j] > thresh: # # Top border @@ -109,13 +149,21 @@ plt.figure() #plt.imshow(matrix, cmap=cm.gray, interpolation='nearest',origin='lower') plt.imshow(bigmatrix, cmap=cm.gray, norm=mcolors.Normalize(0,1), interpolation='nearest',origin='lower') - plt.gca().get_xaxis().set_visible(False) - plt.gca().get_yaxis().set_visible(False) + + #plt.gca().get_xaxis().set_visible(False) + #plt.gca().get_yaxis().set_visible(False) #int_setticks() + if algoname in withticks: + int_setticks() + if algoname in withnoaxes: + plt.gca().get_xaxis().set_visible(False) + plt.gca().get_yaxis().set_visible(False) + if dosave: for ext in saveplotexts: plt.savefig(saveplotbase + algoname + ('_lbd%.0e' % lambdas[ilbd]) + '.' + ext, bbox_inches='tight') - ilbd = ilbd + 1 + + ilbd = ilbd + 1 if doshow: plt.show() print "Finished." @@ -154,6 +202,7 @@ bigmatrix[i*N:i*N+N,j*N+N-border:j*N+N] = bordercolor return bigmatrix + def int_setticks():