Add more util tests and remove more pdb

Author: mjhough

dfols/params.py

@@ -113,7 +113,7 @@ def __init__(self, n, npt, maxfun, objfun_has_noise=False):
         self.params["dykstra.d_tol"] = 1e-10
         self.params["dykstra.max_iters"] = 100
         # Matrix rank algorithm
-        self.params["matrix_rank.r_tol"] = 1e-16
+        self.params["matrix_rank.r_tol"] = 1e-15
 
         self.params_changed = {}
         for p in self.params:

dfols/tests/test_util.py

@@ -165,4 +165,135 @@ def runTest(self):
             self.assertTrue(np.linalg.norm(dirns[i, :]) <= delta + 1e-10, "Unconstrained: dirn %i too long" % i)
             self.assertTrue(np.all(dirns[i, :] >= lower), "Direction %i below lower bound" % i)
             self.assertTrue(np.all(dirns[i, :] <= upper), "Direction %i above upper bound" % i)
-        # self.assertTrue(False, "bad")
+        # self.assertTrue(False, "bad")
+
+# Trivial case of full rank
+class TestMatrixRankQR1(unittest.TestCase):
+    def runTest(self):
+        mr_tol = 1e-15
+        A = np.array([
+            [1,0,0,0],
+            [0,1,0,0],
+            [0,0,1,0],
+            [0,0,0,1]])
+        rank, D = qr_rank(A,mr_tol)
+        self.assertTrue(np.all(D > mr_tol), "Incorrect diagonal matrix output")
+        self.assertTrue(rank == 4, "Incorrect rank output")
+
+# Full rank but QR has negative entries for diag(R)
+class TestMatrixRankQR2(unittest.TestCase):
+    def runTest(self):
+        mr_tol = 1e-15
+        A = np.array([
+            [1,2,3,4],
+            [0,6,7,8],
+            [-1,-2,-2,-1],
+            [4,2,2,1]])
+        rank, D = qr_rank(A,mr_tol)
+        self.assertTrue(np.all(D > mr_tol), "Incorrect diagonal matrix output")
+        self.assertTrue(rank == 4, "Incorrect rank output")
+
+
+# Full rank but QR has negative entries for diag(R)
+class TestMatrixRankQR3(unittest.TestCase):
+    def runTest(self):
+        mr_tol = 1e-15
+        A = np.array([
+            [1,2,3,4],
+            [2,6,4,8],
+            [-1,-2,-3,-4],
+            [1,3,2,4]])
+        rank, D = qr_rank(A,mr_tol)
+        self.assertTrue(np.all(D[0:2] > mr_tol), "Incorrect diagonal matrix output (rows 1,2)")
+        self.assertTrue(np.all(D[2:4] <= mr_tol), "Incorrect diagonal matrix output (rows 3,4)")
+        self.assertTrue(rank == 2, "Incorrect rank output")
+
+
+class TestDykstraBoxInt(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([0,0])
+        lower = np.array([-0.01, -0.1])
+        upper = np.array([0.01, 0.5])
+        boxproj = lambda x: pbox(x,lower,upper)
+        P = [boxproj]
+        xproj = dykstra(P,x0)
+        self.assertTrue(np.all(xproj == x0), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxExt(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([-2,5])
+        lower = np.array([-1, -1])
+        upper = np.array([0.5, 0.9])
+        boxproj = lambda x: pbox(x,lower,upper)
+        P = [boxproj]
+        xproj = dykstra(P,x0)
+        xtrue = np.array([-1,0.9])
+        self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+class TestDykstraBallInt(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([0,0])
+        ballproj = lambda x: pball(x,x0+1,2)
+        P = [ballproj]
+        xproj = dykstra(P,x0)
+        self.assertTrue(np.all(xproj == x0), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBallExt(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([-3,5])
+        ballproj = lambda x: pball(x,np.array([-0.5,1]),1)
+        P = [ballproj]
+        xproj = dykstra(P,x0)
+        xtrue = np.array([-1.02999894, 1.8479983])
+        self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxBallInt(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([0.72,1.1])
+        lower = np.array([0.7, -2.0])
+        upper = np.array([1.0, 2])
+        boxproj = lambda x: pbox(x,lower,upper)
+        ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+        P = [boxproj,ballproj]
+        xproj = dykstra(P,x0)
+        self.assertTrue(np.all(xproj == x0), "Incorrect point returned by Dykstra")
+
+class TestDykstraBoxBallExt1(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([0,4])
+        lower = np.array([0.7, -2.0])
+        upper = np.array([1.0, 2])
+        boxproj = lambda x: pbox(x,lower,upper)
+        ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+        P = [boxproj,ballproj]
+        xproj = dykstra(P,x0)
+        xtrue = np.array([0.6940582, 1.1576116])
+        self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxBallExt2(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([0.8,-3])
+        lower = np.array([0.7, -2.0])
+        upper = np.array([1.0, 2])
+        boxproj = lambda x: pbox(x,lower,upper)
+        ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+        P = [boxproj,ballproj]
+        xproj = dykstra(P,x0)
+        xtrue = np.array([0.68976232, 0.8372417])
+        self.assertTrue(np.allclose(xproj, xtrue), "Incorrect point returned by Dykstra")
+
+
+class TestDykstraBoxBallBdry(unittest.TestCase):
+    def runTest(self):
+        x0 = np.array([0.7,0.85])
+        lower = np.array([0.7, -2.0])
+        upper = np.array([1.0, 2])
+        boxproj = lambda x: pbox(x,lower,upper)
+        ballproj = lambda x: pball(x,np.array([0.5,1]),0.25)
+        P = [boxproj,ballproj]
+        xproj = dykstra(P,x0)
+        self.assertTrue(np.allclose(xproj, x0), "Incorrect point returned by Dykstra")

dfols/util.py

@@ -30,8 +30,6 @@
 import scipy.linalg as LA
 import sys
 
-import pdb
-
 
 __all__ = ['sumsq', 'eval_least_squares_objective', 'model_value', 'random_orthog_directions_within_bounds',
            'random_directions_within_bounds', 'apply_scaling', 'remove_scaling', 'pbox', 'pball', 'dykstra', 'qr_rank']
@@ -211,7 +209,7 @@ def remove_scaling(x_scaled, scaling_changes):
     return shift + x_scaled * scale
 
 
-def dykstra(P,x0,max_iter=100,tol=1.0e-10):
+def dykstra(P,x0,max_iter=100,tol=1e-10):
     x = x0.copy()
     p = len(P)
     y = np.zeros((p,x0.shape[0]))
@@ -226,7 +224,6 @@ def dykstra(P,x0,max_iter=100,tol=1.0e-10):
             x = P[i](prev_x - y[i,:])
 
             # Update increment
-            #  pdb.set_trace()
             prev_y = y[i,:].copy()
             y[i,:] = x - (prev_x - prev_y)
 
@@ -251,7 +248,7 @@ def pbox(x,l,u):
 can be given by the number of nonzero diagonal elements of
 R in the QR factorization of A.
 '''
-def qr_rank(A,tol=1e-25):
+def qr_rank(A,tol=1e-15):
     m,n = A.shape
     assert m == n, "Input matrix must be square"
     Q,R = LA.qr(A)