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