nikcleju@2: """ 
nikcleju@2: #=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#
nikcleju@2: # Bob L. Sturm <bst@create.aau.dk> 20111018
nikcleju@2: # Department of Architecture, Design and Media Technology
nikcleju@2: # Aalborg University Copenhagen
nikcleju@2: # Lautrupvang 15, 2750 Ballerup, Denmark
nikcleju@2: #=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#=#
nikcleju@2: """
nikcleju@2: 
nikcleju@2: import numpy as np
nikcleju@2: from sklearn.utils import check_random_state
nikcleju@2: import time
nikcleju@2: 
nikcleju@2: from omp_sk_bugfix import orthogonal_mp
nikcleju@2: from omp_QR import greed_omp_qr
nikcleju@2: from omp_QR import omp_qr
nikcleju@2: 
nikcleju@2: """
nikcleju@2: Run a problem suite involving sparse vectors in 
nikcleju@2: ambientDimension dimensional space, with a resolution
nikcleju@2: in the phase plane of numGradations x numGradations,
nikcleju@2: and at each indeterminacy and sparsity pair run
nikcleju@2: numTrials independent trials.
nikcleju@2: 
nikcleju@2: Outputs a text file denoting successes at each phase point.
nikcleju@2: For more on phase transitions, see:
nikcleju@2: D. L. Donoho and J. Tanner, "Precise undersampling theorems," 
nikcleju@2: Proc. IEEE, vol. 98, no. 6, pp. 913-924, June 2010.
nikcleju@2: """
nikcleju@2: 
nikcleju@2: def runProblemSuite(ambientDimension,numGradations,numTrials):
nikcleju@2: 
nikcleju@2:   idx = np.arange(ambientDimension)
nikcleju@2:   phaseDelta = np.linspace(0.05,1,numGradations)
nikcleju@2:   phaseRho = np.linspace(0.05,1,numGradations)
nikcleju@2:   success = np.zeros((numGradations, numGradations))
nikcleju@2:   
nikcleju@2:   #Nic: init timers
nikcleju@2:   t1all = 0
nikcleju@2:   t2all = 0
nikcleju@2:   t3all = 0
nikcleju@2:   
nikcleju@2:   deltaCounter = 0
nikcleju@2:   # delta is number of measurements/
nikcleju@2:   for delta in phaseDelta[:17]:
nikcleju@2:     rhoCounter = 0
nikcleju@2:     for rho in phaseRho:
nikcleju@2:       print(deltaCounter,rhoCounter)
nikcleju@2:       numMeasurements = int(delta*ambientDimension)
nikcleju@2:       sparsity = int(rho*numMeasurements)
nikcleju@2:       # how do I set the following to be random each time?
nikcleju@2:       generator = check_random_state(100)
nikcleju@2:       # create unit norm dictionary
nikcleju@2:       D = generator.randn(numMeasurements, ambientDimension)
nikcleju@2:       D /= np.sqrt(np.sum((D ** 2), axis=0))
nikcleju@2:       # compute Gramian (for efficiency)
nikcleju@2:       DTD = np.dot(D.T,D)
nikcleju@2:   
nikcleju@2:       successCounter = 0
nikcleju@2:       trial = numTrials
nikcleju@2:       while trial > 0:
nikcleju@2:         # generate sparse signal with a minimum non-zero value
nikcleju@2:         x = np.zeros((ambientDimension, 1))
nikcleju@2:         idx2 = idx
nikcleju@2:         generator.shuffle(idx2)
nikcleju@2:         idx3 = idx2[:sparsity]
nikcleju@2:         while np.min(np.abs(x[idx3,0])) < 1e-10 :
nikcleju@2:            x[idx3,0] = generator.randn(sparsity)
nikcleju@2:         # sense sparse signal
nikcleju@2:         y = np.dot(D, x)
nikcleju@2:         
nikcleju@2:         # Nic: Use sparsify OMP function (translated from Matlab)
nikcleju@2:         ompopts = dict({'stopCrit':'M', 'stopTol':2*sparsity})
nikcleju@2:         starttime = time.time()                     # start timer
nikcleju@2:         x_r2, errs, times = greed_omp_qr(y.squeeze().copy(), D.copy(), D.shape[1], ompopts)
nikcleju@2:         t2all = t2all + time.time() - starttime     # stop timer
nikcleju@2:         idx_r2 = np.nonzero(x_r2)[0]        
nikcleju@2:         
nikcleju@2:         # run to two times expected sparsity, or tolerance
nikcleju@2:         # why? Often times, OMP can retrieve the correct solution
nikcleju@2:         # when it is run for more than the expected sparsity
nikcleju@2:         #x_r, idx_r = omp_qr(y,D,DTD,2*sparsity,1e-5)
nikcleju@2:         # Nic: adjust tolerance to match with other function
nikcleju@2:         starttime = time.time()                     # start timer
nikcleju@2:         x_r, idx_r = omp_qr(y.copy(),D.copy(),DTD.copy(),2*sparsity,numMeasurements*1e-14/np.vdot(y,y))
nikcleju@2:         t1all = t1all + time.time() - starttime     # stop timer        
nikcleju@2:         
nikcleju@2:         # Nic: test sklearn omp
nikcleju@2:         starttime = time.time()                     # start timer
nikcleju@2:         x_r3 = orthogonal_mp(D.copy(), y.copy(), n_nonzero_coefs=2*sparsity, tol=numMeasurements*1e-14, precompute_gram=False, copy_X=True)
nikcleju@2:         idx_r3 = np.nonzero(x_r3)[0]
nikcleju@2:         t3all = t3all + time.time() - starttime     # stop timer        
nikcleju@2:         
nikcleju@2:         # Nic: compare results
nikcleju@2:         print 'diff1 = ',np.linalg.norm(x_r.squeeze() - x_r2.squeeze())
nikcleju@2:         print 'diff2 = ',np.linalg.norm(x_r.squeeze() - x_r3.squeeze())
nikcleju@2:         print 'diff3 = ',np.linalg.norm(x_r2.squeeze() - x_r3.squeeze())
nikcleju@2:         print "Bob's total time = ", t1all
nikcleju@2:         print "Nic's total time = ", t2all
nikcleju@2:         print "Skl's total time = ", t3all
nikcleju@2:         if np.linalg.norm(x_r.squeeze() - x_r2.squeeze()) > 1e-10 or \
nikcleju@2:            np.linalg.norm(x_r.squeeze() - x_r3.squeeze()) > 1e-10 or \
nikcleju@2:            np.linalg.norm(x_r2.squeeze() - x_r3.squeeze()) > 1e-10:
nikcleju@2:             print "STOP: Different results"
nikcleju@2:             print "Bob's residual: ||y - D x_r ||_2 = ",np.linalg.norm(y.squeeze() - np.dot(D,x_r).squeeze())
nikcleju@2:             print "Nic's residual: ||y - D x_r ||_2 = ",np.linalg.norm(y.squeeze() - np.dot(D,x_r2).squeeze())
nikcleju@2:             print "Skl's residual: ||y - D x_r ||_2 = ",np.linalg.norm(y.squeeze() - np.dot(D,x_r3).squeeze())
nikcleju@2:             raise ValueError("Different results")
nikcleju@2:   
nikcleju@2:         # debais to remove small entries
nikcleju@2:         for nn in idx_r:
nikcleju@2:           if abs(x_r[nn]) < 1e-10:
nikcleju@2:             x_r[nn] = 0
nikcleju@2:   
nikcleju@2:         # exact recovery condition using support
nikcleju@2:         #if sorted(np.flatnonzero(x_r)) == sorted(np.flatnonzero(x)):
nikcleju@2:         #  successCounter += 1
nikcleju@2:         # exact recovery condition using error in solution
nikcleju@2:         error = x - x_r
nikcleju@2:         """ the following is the exact recovery condition in: A. Maleki 
nikcleju@2:               and D. L. Donoho, "Optimally tuned iterative reconstruction 
nikcleju@2:               algorithms for compressed sensing," IEEE J. Selected Topics 
nikcleju@2:               in Signal Process., vol. 4, pp. 330-341, Apr. 2010. """
nikcleju@2:         if np.vdot(error,error) < np.vdot(x,x)*1e-4:
nikcleju@2:           successCounter += 1
nikcleju@2:         trial -= 1
nikcleju@2:   
nikcleju@2:       success[rhoCounter,deltaCounter] = successCounter
nikcleju@2:       if successCounter == 0:
nikcleju@2:         break
nikcleju@2:   
nikcleju@2:       rhoCounter += 1
nikcleju@2:       #np.savetxt('test.txt',success,fmt='#2.1d',delimiter=',')
nikcleju@2:     deltaCounter += 1
nikcleju@2: 
nikcleju@2: if __name__ == '__main__':
nikcleju@2:   print ('Running problem suite')
nikcleju@2:   ambientDimension = 400
nikcleju@2:   numGradations = 30
nikcleju@2:   numTrials = 1
nikcleju@2:   
nikcleju@2:   #import cProfile
nikcleju@2:   #cProfile.run('runProblemSuite(ambientDimension,numGradations,numTrials)','profres')  
nikcleju@2:   runProblemSuite(ambientDimension,numGradations,numTrials) 
nikcleju@2:   print "Done"
nikcleju@2:   
nikcleju@2:   #import pstats
nikcleju@2:   #p = pstats.Stats('D:\Nic\Dev2\profres')
nikcleju@2:   #p.sort_stats('cumulative').print_stats(10)