nikcleju@2: # -*- coding: utf-8 -*-
nikcleju@2: """
nikcleju@2: Created on Fri Oct 21 14:28:08 2011
nikcleju@2: 
nikcleju@2: @author: Nic
nikcleju@3: 
nikcleju@3: Test BP algorithm
nikcleju@2: """
nikcleju@2: 
nikcleju@3: import numpy as np
nikcleju@3: import numpy.linalg
nikcleju@3: import scipy.io
nikcleju@3: import unittest
nikcleju@3: #import sys
nikcleju@3: #sys.path.append("D:\Nic\Dev\pyCSalgos\trunk")
nikcleju@3: #sys.path.append("D:\Nic\Dev\pyCSalgos\trunk\pyCSalgos\BP")
nikcleju@3: #import pyCSalgos
nikcleju@3: from pyCSalgos.BP.l1qc import l1qc_logbarrier
nikcleju@3: #from l1qc import l1qc_logbarrier
nikcleju@3: 
nikcleju@3: class BPresults(unittest.TestCase):
nikcleju@3:   def testResults(self):
nikcleju@4:     mdict = scipy.io.loadmat('BPtestdata.mat')
nikcleju@3:     
nikcleju@3:     # A = system matrix
nikcleju@3:     # Y = matrix with measurements (on columns)
nikcleju@3:     # X0 = matrix with initial solutions (on columns)
nikcleju@3:     # Eps = vector with epsilon
nikcleju@3:     # Xr = matrix with correct solutions (on columns)
nikcleju@3:     for A,Y,X0,Eps,Xr in zip(mdict['cellA'].squeeze(),mdict['cellY'].squeeze(),mdict['cellX0'].squeeze(),mdict['cellEps'].squeeze(),mdict['cellXr'].squeeze()):
nikcleju@3:       for i in np.arange(Y.shape[1]):
nikcleju@3:         xr = l1qc_logbarrier(X0[:,i], A, np.array([]), Y[:,i], Eps.squeeze()[i])
nikcleju@3:         
nikcleju@3:         # check if found solution is the same as the correct cslution
nikcleju@3:         diff = numpy.linalg.norm(xr - Xr[:,i])
nikcleju@3:         err1 = numpy.linalg.norm(Y[:,i] - np.dot(A,xr))
nikcleju@3:         err2 = numpy.linalg.norm(Y[:,i] - np.dot(A,Xr[:,i]))
nikcleju@3:         norm1 = numpy.linalg.norm(xr,1)
nikcleju@3:         norm2 = numpy.linalg.norm(Xr[:,i],1)
nikcleju@3:         print 'diff = ',diff
nikcleju@3:         print 'err1 = ',err1
nikcleju@3:         print 'err2 = ',err2
nikcleju@3:         print 'norm1 = ',norm1
nikcleju@3:         print 'norm2 = ',norm2
nikcleju@3:         
nikcleju@3:         # It seems Matlab's linsolve and scipy solve are slightly different
nikcleju@3:         # Therefore make a more robust condition:
nikcleju@3:         #  OK;    if   solutions are close enough (diff < 1e-6)
nikcleju@3:         #              or
nikcleju@3:         #              (
nikcleju@3:         #               they fulfill the constraint close enough (differr < 1e-6)
nikcleju@3:         #                 and
nikcleju@3:         #               Python solution has l1 norm no more than 1e-6 larger as the reference solution
nikcleju@3:         #                 (i.e. either norm1 < norm2   or   norm1>norm2 not by more than 1e-6)
nikcleju@3:         #              )
nikcleju@3:         #        
nikcleju@3:         #  ERROR: else
nikcleju@3:         differr  = abs((err1 - err2))
nikcleju@3:         diffnorm = norm1 - norm2  # intentionately no abs(), since norm1 < norm2 is good
nikcleju@3:         if diff < 1e-6 or (differr < 1e-6 and (diffnorm < 1e-6)):
nikcleju@3:           isok = True
nikcleju@3:         else:
nikcleju@3:           isok = False
nikcleju@3: 
nikcleju@4:         #if not isok:
nikcleju@4:         #  print "should raise"
nikcleju@4:         #  #self.assertTrue(isok)
nikcleju@4:         self.assertTrue(isok)
nikcleju@3:   
nikcleju@3: if __name__ == "__main__":
nikcleju@3:     unittest.main(verbosity=2)    
nikcleju@3:     #suite = unittest.TestLoader().loadTestsFromTestCase(CompareResults)
nikcleju@3:     #unittest.TextTestRunner(verbosity=2).run(suite)