nikcleju@47: # -*- coding: utf-8 -*-
nikcleju@47: """
nikcleju@47: Created on Sun Nov 06 20:53:14 2011
nikcleju@47: 
nikcleju@47: @author: Nic
nikcleju@47: """
nikcleju@47: import numpy as np
nikcleju@47: import numpy.linalg
nikcleju@47: import scipy.io
nikcleju@47: import unittest
nikcleju@47: from pyCSalgos.NESTA.NESTA import NESTA
nikcleju@47: 
nikcleju@47: class NESTAresults(unittest.TestCase):
nikcleju@47:   def testResults(self):
nikcleju@47:     mdict = scipy.io.loadmat('NESTAtestdata.mat')
nikcleju@47:     
nikcleju@47:     # Add [0,0] indices because data is read from mat file as [1,1] arrays
nikcleju@47:     opt_TolVar = mdict['opt_TolVar'][0,0]
nikcleju@47:     opt_Verbose = mdict['opt_Verbose'][0,0]
nikcleju@47:     opt_muf = mdict['opt_muf'][0,0]
nikcleju@47:     numA = mdict['numA'][0,0]
nikcleju@47:     
nikcleju@47:     # Known bad but good:
nikcleju@47:     known = ()
nikcleju@47:       
nikcleju@47:     sumplus  = 0.0
nikcleju@47:     summinus = 0.0
nikcleju@47:     numplus = 0
nikcleju@47:     numminus = 0    
nikcleju@47:     
nikcleju@47:     # A = system matrix
nikcleju@47:     # Y = matrix with measurements (on columns)
nikcleju@47:     # sigmamin = vector with sigma_mincell
nikcleju@47:     for k,A,Y,M,eps,Xr in zip(np.arange(numA),mdict['cellA'].squeeze(),mdict['cellY'].squeeze(),mdict['cellM'].squeeze(),mdict['cellEps'].squeeze(),mdict['cellXr'].squeeze()):
nikcleju@47: 
nikcleju@47:       # Fix numpy error "LapackError: Parameter a has non-native byte order in lapack_lite.dgesdd"
nikcleju@47:       A = A.newbyteorder('=')
nikcleju@47:       Y = Y.newbyteorder('=')
nikcleju@47:       M = M.newbyteorder('=')
nikcleju@47:       eps = eps.newbyteorder('=')
nikcleju@47:       Xr = Xr.newbyteorder('=')
nikcleju@47:       
nikcleju@47:       eps = eps.squeeze()
nikcleju@47:       
nikcleju@47:       U,S,V = numpy.linalg.svd(M, full_matrices = True)
nikcleju@47:       V = V.T         # Make like Matlab
nikcleju@47:       m,n = M.shape   # Make like Matlab
nikcleju@47:       S = numpy.hstack((numpy.diag(S), numpy.zeros((m,n-m))))
nikcleju@47: 
nikcleju@47:       optsUSV = {'U':U, 'S':S, 'V':V}
nikcleju@47:       opts = {'U':A, 'Ut':A.T.copy(), 'USV':optsUSV, 'TolVar':opt_TolVar, 'Verbose':opt_Verbose}
nikcleju@47:       
nikcleju@47:       for i in np.arange(Y.shape[1]):
nikcleju@47:         xr = NESTA(M, None, Y[:,i], opt_muf, eps[i] * numpy.linalg.norm(Y[:,i]), opts)[0]
nikcleju@47:         
nikcleju@47:         # check if found solution is the same as the correct cslution
nikcleju@47:         diff = numpy.linalg.norm(xr - Xr[:,i])
nikcleju@47:         print "k =",k,"i = ",i
nikcleju@47:         if diff < 1e-6:
nikcleju@47:           print "Recovery OK"
nikcleju@47:           isOK = True
nikcleju@47:         else:
nikcleju@47:           if numpy.linalg.norm(xr,1) < numpy.linalg.norm(Xr[:,i],1):
nikcleju@47:             numplus = numplus+1
nikcleju@47:             sumplus = sumplus + numpy.linalg.norm(Xr[:,i],1) - numpy.linalg.norm(xr,1)
nikcleju@47:           else:
nikcleju@47:             numminus = numminus+1
nikcleju@47:             summinus = summinus + numpy.linalg.norm(xr,1) - numpy.linalg.norm(Xr[:,i],1)
nikcleju@47:          
nikcleju@47:           print "Oops"
nikcleju@47:           if (k,i) not in known:
nikcleju@47:             #isOK = False
nikcleju@47:             print "Should stop here"
nikcleju@47:           else:
nikcleju@47:             print "Known bad but good"
nikcleju@47:             isOK = True
nikcleju@47:         # comment / uncomment this
nikcleju@47:         self.assertTrue(isOK)
nikcleju@47:     print 'Finished test'
nikcleju@47:   
nikcleju@47: if __name__ == "__main__":
nikcleju@47:     #import cProfile
nikcleju@47:     #cProfile.run('unittest.main()', 'profres')
nikcleju@47:     unittest.main()    
nikcleju@47:     #suite = unittest.TestLoader().loadTestsFromTestCase(CompareResults)
nikcleju@47:     #unittest.TextTestRunner(verbosity=2).run(suite)