view test_exact_old.py @ 21:d395461b92ae tip

Lots and lots of modifications. Approximate recovery script working.

# -*- coding: utf-8 -*-
"""
Old version, ignore it

Author: Nicolae Cleju
"""

import numpy
import scipy.io
import math
import os
import time

import multiprocessing
import sys
_currmodule = sys.modules[__name__]
# Lock for printing in a thread-safe way
_printLock = None

import stdparams_exact
import AnalysisGenerate

# For exceptions
import pyCSalgos.BP.l1eq_pd
import pyCSalgos.NESTA.NESTA


def _initProcess(share, njobs, printLock):
    """
    Pool initializer function (multiprocessing)
    Needed to pass the shared variable to the worker processes
    The variables must be global in the module in order to be seen later in run_once_tuple()
    see http://stackoverflow.com/questions/1675766/how-to-combine-pool-map-with-array-shared-memory-in-python-multiprocessing
    """
    currmodule = sys.modules[__name__]
    currmodule.proccount = share
    currmodule.njobs = njobs
    currmodule._printLock = printLock
          
#==========================
# Interface run functions
#==========================
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):

  print "This is analysis recovery ABS approximation script by Nic"
  print "Running phase transition ( run_multi() )"

  if ncpus is None:
      print "  Running in parallel with default",multiprocessing.cpu_count(),"threads using \"multiprocessing\" package"
      if multiprocessing.cpu_count() == 1:
          doparallel = False
      else:
          doparallel = True
  elif ncpus > 1:
      print "  Running in parallel with",ncpus,"threads using \"multiprocessing\" package"
      doparallel = True
  elif ncpus == 1:
      print "Running single thread"
      doparallel = False
  else:
      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))

  print "End of parameters"  

  _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]]

  # 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()
    
  print "Finished."
  
def run_once_tuple(t):
  results = run_once(*t)
  
  if _currmodule._printLock:
    _currmodule._printLock.acquire()  
    
    _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 "================================"

    _currmodule._printLock.release()

  return results

def run_once(algos,Omega,y,M,x0):
  
  d = Omega.shape[1]  
  
  nalgos = len(algos)  
 
  xrec = dict()
  err = dict()
  relerr = dict()
  elapsed = dict()

  # Prepare storage variables for algorithms
  for i,algo in zip(numpy.arange(nalgos),algos):
    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
  
  # Run algorithms
  for iy in numpy.arange(y.shape[1]):
    for algofunc,strname in algos:
      try:
        timestart = time.time()
        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()
          print "Caught exception when running algorithm",strname," :",e.message
          _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()
      print strname,' : avg relative error = ',numpy.mean(relerr[strname])
      _currmodule._printLock.release()

  # Prepare results
  #mrelerr = dict()
  #for algotuple in algos:
  #  mrelerr[algotuple[1]] = numpy.mean(relerr[algotuple[1]])
  #return mrelerr,elapsed
  
  # Should return number of reconstructions with error < threshold, not average error
  exactthr = 1e-6
  mrelerr = dict()
  for algotuple in algos:
    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")
    algos = stdparams_exact.std1()[0]
    res = run_once(algos, mdict['Omega'].byteswap().newbyteorder(),mdict['y'],mdict['M'],mdict['x0'])
  
def generateFig():
  run(stdparams_exact.std1)
  
# Script main
if __name__ == "__main__":
  #import cProfile
  #cProfile.run('mainrun()', 'profile')    
  #run_mp(stdparams_exact.stdtest)
  #runsingleexampledebug()
  
  run(stdparams_exact.std1, ncpus=3)
  #testMatlab()
  #run(stdparams_exact.stdtest, ncpus=1)