Improve performance of simple finite differences.

Author: Benjamin Chrétien

src/roboptim/core/__init__.py

@@ -176,6 +176,9 @@ def impl_compute (self, result, x):
     def impl_gradient (self, result, x, functionId):
         gradient (self._fd, result, x, functionId)
 
+    def impl_jacobian (self, result, x):
+        jacobian (self._fd, result, x)
+
 
 class PyProblem(object):
     def __init__(self, cost):

src/wrap.hh

@@ -194,6 +194,13 @@ namespace roboptim
 	{
 	  fd_t::impl_gradient (gradient, argument, functionId);
 	}
+
+        virtual void impl_jacobian (jacobian_ref jacobian,
+                                    const_argument_ref argument)
+          const
+	{
+	  fd_t::impl_jacobian (jacobian, argument);
+	}
       };
 
 

tests/function_pool.py

@@ -15,6 +15,12 @@ def __init__ (self, n):
         self.n = n
         self.data = np.zeros (n)
         self.jac = np.zeros ((n, 2*n))
+        self.compute_counter = 0
+        self.jacobian_counter = 0
+
+    def reset (self):
+        self.compute_counter = 0
+        self.jacobian_counter = 0
 
     def impl_compute (self, result, x):
         print("Engine: impl_compute")
@@ -33,11 +39,14 @@ def jacobian (self, x):
 
     def computeData (self, x):
         print("Engine: computing data")
+        self.compute_counter += 1
+        # For each square function
         for i in range(self.n):
             self.data[i] = x[2*i]**2 + x[2*i+1]**2
 
     def computeJacobian (self, x):
         print("Engine: computing Jacobian")
+        self.jacobian_counter += 1
         self.jac.fill (0.)
         # For each square function
         for i in range(self.n):
@@ -77,6 +86,9 @@ def test_engine(self):
         engine (x)
         np.testing.assert_almost_equal (engine.data,
                 [xi**2 + yi**2 for xi,yi in x.reshape(engine.n, 2) ])
+        assert engine.compute_counter == 1
+        assert engine.jacobian_counter == 0
+        engine.reset ()
 
         np.testing.assert_almost_equal (engine.jac, np.zeros ((engine.n, 2*engine.n)))
         engine.jacobian (x)
@@ -85,46 +97,71 @@ def test_engine(self):
             for j in range(2):
                 jac[i,2*i+j] = 2. * x[2*i+j]
         np.testing.assert_almost_equal (engine.jac, jac)
+        assert engine.compute_counter == 0
+        assert engine.jacobian_counter == 1
+        engine.reset ()
 
     def test_pool(self):
-        engine = Engine (3)
+        n = 3
+        engine = Engine (n)
         np.testing.assert_almost_equal (engine.data, np.zeros (engine.n))
-        functions = [Square (engine, 0.), Square (engine, 1.), Square (engine, 2.)]
+        functions = [Square (engine, float(i)) for i in range (n)]
         print(engine)
+
         pool = roboptim.core.PyFunctionPool (engine, functions, "Dummy pool")
         print(pool)
+
         x = np.array([10., -5., 1., 2., -1., 1.])
+        assert len(x) == 2 * n
+
         res = pool (x)
         np.testing.assert_almost_equal (engine.data,
                 [xi**2 + yi**2 for xi,yi in x.reshape(engine.n, 2) ])
+        assert engine.compute_counter == 1
+        assert engine.jacobian_counter == 0
+        engine.reset ()
+
         pool_jac = pool.jacobian (x)
         jac = np.zeros ((engine.n, 2*engine.n))
         for i in range(engine.n):
             for j in range(2):
                 jac[i,2*i+j] = 2. * x[2*i+j]
         np.testing.assert_almost_equal (pool_jac, jac)
+        assert engine.compute_counter == 0
+        assert engine.jacobian_counter == 1
+        engine.reset ()
 
     def test_pool_fd(self):
-        engine = Engine (3)
+        n = 3
+        engine = Engine (n)
         print(engine)
-        functions = [Square (engine, 0.), Square (engine, 1.), Square (engine, 2.)]
+        functions = [Square (engine, float(i)) for i in range (n)]
         pool = roboptim.core.PyFunctionPool (engine, functions, "Dummy FD pool")
 
         fd_rule = roboptim.core.FiniteDifferenceRule.SIMPLE
         fd_pool = roboptim.core.PyFiniteDifference (pool, rule = fd_rule)
         print(fd_pool)
 
         x = np.array([10., -5., 1., 2., -1., 1.])
+        assert len(x) == 2 * n
+
         res = fd_pool (x)
         np.testing.assert_almost_equal (engine.data,
                 [xi**2 + yi**2 for xi,yi in x.reshape(engine.n, 2) ])
+        assert engine.compute_counter == 1
+        assert engine.jacobian_counter == 0
+        engine.reset ()
 
         fd_pool_jac = fd_pool.jacobian (x)
         jac = np.zeros ((engine.n, 2*engine.n))
         for i in range(engine.n):
             for j in range(2):
                 jac[i,2*i+j] = 2. * x[2*i+j]
         np.testing.assert_almost_equal (fd_pool_jac, jac, 5)
+        print(engine.compute_counter)
+        assert engine.compute_counter == 1 + len(x) # simple rule: (f(x+h)-f(x))/h
+        assert engine.jacobian_counter == 0
+        engine.reset ()
 
 if __name__ == '__main__':
     unittest.main()