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