Set use_numpy at init of Lambdify (6x performance gain)

Author: bjodah

benchmarks/Lambdify.py

@@ -18,17 +18,20 @@
 
 inp = np.ones(28)
 
+lmb_sp(*inp)
 tim_sympy = clock()
 for i in range(500):
     res_sympy = lmb_sp(*inp)
 tim_sympy = clock() - tim_sympy
 
+lbm_se(inp)
 tim_se = clock()
 res_se = np.empty(len(exprs))
 for i in range(500):
     res_se = lbm_se(inp)
 tim_se = clock() - tim_se
 
+lbm_se_llvm(inp)
 tim_se_llvm = clock()
 res_se_llvm = np.empty(len(exprs))
 for i in range(500):
@@ -62,6 +65,7 @@ def func(*args):
     return func
 
 lbm_se_llvm_manual = ManualLLVM(args, np.array(exprs))
+lbm_se_llvm_manual(inp)
 tim_se_llvm_manual = clock()
 res_se_llvm_manual = np.empty(len(exprs))
 for i in range(500):
@@ -72,3 +76,16 @@ def func(*args):
 
 if tim_se_llvm_manual < tim_se_llvm:
     warnings.warn("Cython code for Lambdify.__call__ is slow.")
+
+import setuptools
+import pyximport
+pyximport.install()
+from Lambdify_reference import _benchmark_reference_for_Lambdify as lmb_ref
+
+lmb_ref(inp)
+tim_ref = clock()
+for i in range(500):
+    res_ref = lmb_ref(inp)
+tim_ref = clock() - tim_ref
+print('Hard-coded Cython code    speed-up factor (higher is better) vs sympy: %12.5g' %
+      (tim_sympy/tim_ref))

benchmarks/Lambdify_reference.pyx

@@ -0,0 +1,10 @@
+# Benchmark reference:
+cimport numpy as cnp
+import numpy as np
+from libc.math cimport exp as exp_
+
+def _benchmark_reference_for_Lambdify(cnp.ndarray[cnp.float64_t] x):
+    cdef cnp.ndarray[cnp.float64_t] out = np.empty(14)
+    cdef double * data = <double *>out.data
+    data[:] = [x[0] + x[1] - x[4] + 36.252574322669, x[0] - x[2] + x[3] + 21.3219379611249, x[3] + x[5] - x[6] + 9.9011158998744, 2*x[3] + x[5] - x[7] + 18.190422234653, 3*x[3] + x[5] - x[8] + 24.8679190043357, 4*x[3] + x[5] - x[9] + 29.9336062089226, -x[10] + 5*x[3] + x[5] + 28.5520551531262, 2*x[0] + x[11] - 2*x[4] - 2*x[5] + 32.4401680272417, 3*x[1] - x[12] + x[5] + 34.9992934135095, 4*x[1] - x[13] + x[5] + 37.0716199972041, x[14+0] - x[14+1] + 2*x[14+10] + 2*x[14+11] - x[14+12] - 2*x[14+13] + x[14+2] + 2*x[14+5] + 2*x[14+6] + 2*x[14+7] + 2*x[14+8] + 2*x[14+9] - exp_(x[0]) + exp_(x[1]) - 2*exp_(x[10]) - 2*exp_(x[11]) + exp_(x[12]) + 2*exp_(x[13]) - exp_(x[2]) - 2*exp_(x[5]) - 2*exp_(x[6]) - 2*exp_(x[7]) - 2*exp_(x[8]) - 2*exp_(x[9]), -x[14+0] - x[14+1] - 15*x[14+10] - 2*x[14+11] - 3*x[14+12] - 4*x[14+13] - 4*x[14+2] - 3*x[14+3] - 2*x[14+4] - 3*x[14+6] - 6*x[14+7] - 9*x[14+8] - 12*x[14+9] + exp_(x[0]) + exp_(x[1]) + 15*exp_(x[10]) + 2*exp_(x[11]) + 3*exp_(x[12]) + 4*exp_(x[13]) + 4*exp_(x[2]) + 3*exp_(x[3]) + 2*exp_(x[4]) + 3*exp_(x[6]) + 6*exp_(x[7]) + 9*exp_(x[8]) + 12*exp_(x[9]), -5*x[14+10] - x[14+2] - x[14+3] - x[14+6] - 2*x[14+7] - 3*x[14+8] - 4*x[14+9] + 5*exp_(x[10]) + exp_(x[2]) + exp_(x[3]) + exp_(x[6]) + 2*exp_(x[7]) + 3*exp_(x[8]) + 4*exp_(x[9]), -x[14+1] - 2*x[14+11] - 3*x[14+12] - 4*x[14+13] - x[14+4] + exp_(x[1]) + 2*exp_(x[11]) + 3*exp_(x[12]) + 4*exp_(x[13]) + exp_(x[4]), -x[14+10] - 2*x[14+11] - x[14+12] - x[14+13] - x[14+5] - x[14+6] - x[14+7] - x[14+8] - x[14+9] + exp_(x[10]) + 2*exp_(x[11]) + exp_(x[12]) + exp_(x[13]) + exp_(x[5]) + exp_(x[6]) + exp_(x[7]) + exp_(x[8]) + exp_(x[9])]
+    return out

symengine/lib/symengine_wrapper.pyx

@@ -3092,6 +3092,8 @@ cdef class _Lambdify(object):
         the shape of exprs is preserved
     real : bool
         Whether datatype is ``double`` (``double complex`` otherwise).
+    use_numpy: bool (default: None)
+        None -> use numpy if available
 
     Returns
     -------
@@ -3115,9 +3117,22 @@ cdef class _Lambdify(object):
     cdef int *accum_out_sizes
     cdef readonly bool real
     cdef readonly int n_exprs
+    cdef readonly bint use_numpy
+    cdef object _np
 
-    def __cinit__(self, args, *exprs, bool real=True):
+    def __cinit__(self, args, *exprs, bool real=True, use_numpy=None):
         self.real = real
+        if use_numpy is None:
+            try:
+                import numpy as np
+            except ImportError:
+                use_numpy = False  # we will use cython.view.array instead
+            else:
+                use_numpy = True
+        if use_numpy is True:
+            import numpy as np
+            self._np = np
+        self.use_numpy = use_numpy
         self.out_shapes = [get_shape(expr) for expr in exprs]
         self.n_exprs = len(exprs)
         self.args_size = _size(args)
@@ -3136,7 +3151,7 @@ cdef class _Lambdify(object):
         free(self.out_sizes)
         free(self.accum_out_sizes)
 
-    def __init__(self, args, *exprs, bool real=True):
+    def __init__(self, args, *exprs, bool real=True, use_numpy=None):
         cdef:
             Basic e_
             size_t ri, ci, nr, nc
@@ -3197,7 +3212,7 @@ cdef class _Lambdify(object):
             raise ValueError("Size of out incompatible with number of exprs.")
         self.unsafe_complex(inp, out)
 
-    def __call__(self, inp, out=None, use_numpy=None):
+    def __call__(self, inp, out=None):
         """
         Parameters
         ----------
@@ -3208,26 +3223,26 @@ cdef class _Lambdify(object):
             If ``None``: an output container will be allocated (NumPy ndarray or
             cython.view.array). If ``len(exprs) > 0`` output is found in the corresponding
             order. Note that ``out`` is not reshaped.
-        use_numpy: bool (default: None)
-            None -> use numpy if available
 
         Returns
         -------
         If ``len(exprs) == 1``: ``numpy.ndarray`` or ``cython.view.array``, otherwise
         a tuple of such.
 
         """
-        cdef cython.view.array tmp
-        cdef double[::1] real_out_view, real_inp_view
-        cdef double complex[::1] cmplx_out_view, cmplx_inp_view
-        cdef size_t nbroadcast = 1
+        cdef:
+            cython.view.array tmp
+            double[::1] real_out_view, real_inp_view
+            double complex[::1] cmplx_out_view, cmplx_inp_view
+            size_t nbroadcast = 1
+            long inp_size
 
         try:
             inp_shape = getattr(inp, 'shape', (len(inp),))
         except TypeError:
             inp = tuple(inp)
             inp_shape = (len(inp),)
-        inp_size = long(reduce(mul, inp_shape))
+        inp_size = reduce(mul, inp_shape)
         if inp_size % self.args_size != 0:
             raise ValueError("Broadcasting failed")
         nbroadcast = inp_size // self.args_size
@@ -3236,36 +3251,27 @@ cdef class _Lambdify(object):
         new_out_shapes = [inp_shape[:-1] + out_shape for out_shape in self.out_shapes]
         new_out_sizes = [nbroadcast*self.out_sizes[i] for i in range(self.n_exprs)]
         new_tot_out_size = nbroadcast * self.tot_out_size
-        if use_numpy is None:
-            try:
-                import numpy as np
-            except ImportError:
-                use_numpy = False  # we will use cython.view.array instead
-            else:
-                use_numpy = True
-        elif use_numpy is True:
-            import numpy as np
 
-        if use_numpy:
-            numpy_dtype = np.float64 if self.real else np.complex128
+        if self.use_numpy:
+            numpy_dtype = self._np.float64 if self.real else self._np.complex128
             if isinstance(inp, DenseMatrixBase):
-                arr = np.empty(inp_size, dtype=numpy_dtype)
+                arr = self._np.empty(inp_size, dtype=numpy_dtype)
                 if self.real:
                     inp.dump_real(arr)
                 else:
                     inp.dump_complex(arr)
                 inp = arr
             else:
-                inp = np.ascontiguousarray(inp, dtype=numpy_dtype)
+                inp = self._np.ascontiguousarray(inp, dtype=numpy_dtype)
             if inp.ndim > 1:
                 inp = inp.ravel()
         else:
             inp = with_buffer(inp, self.real)
 
         if out is None:
             # allocate output container
-            if use_numpy:
-                out = np.empty(new_tot_out_size, dtype=numpy_dtype)
+            if self.use_numpy:
+                out = self._np.empty(new_tot_out_size, dtype=numpy_dtype)
             else:
                 if self.real:
                     out = cython.view.array((new_tot_out_size,),
@@ -3276,11 +3282,11 @@ cdef class _Lambdify(object):
             reshape_outs = len(new_out_shapes[0]) > 1
         else:
             reshape_outs = False
-            if use_numpy:
+            if self.use_numpy:
                 try:
                     out_dtype = out.dtype
                 except AttributeError:
-                    out = np.asarray(out)
+                    out = self._np.asarray(out)
                     out_dtype = out.dtype
                 if out_dtype != numpy_dtype:
                     raise TypeError("Output array is of incorrect type")
@@ -3310,7 +3316,7 @@ cdef class _Lambdify(object):
                 self.unsafe_complex(complex_inp_view[idx*self.args_size:(idx+1)*self.args_size],
                                     complex_out_view[idx*self.tot_out_size:(idx+1)*self.tot_out_size])
 
-        if use_numpy and reshape_outs:
+        if self.use_numpy and reshape_outs:
             out = out.reshape((nbroadcast, self.tot_out_size))
             result = [out[:, self.accum_out_sizes[idx]:self.accum_out_sizes[idx+1]].reshape(new_out_shapes[idx])
                       for idx in range(self.n_exprs)]
@@ -3375,15 +3381,15 @@ IF HAVE_SYMENGINE_LLVM:
             self.lambda_double[0].call(&out[0], &inp[0])
 
 
-def Lambdify(args, *exprs, bool real=True, backend="lambda"):
+def Lambdify(args, *exprs, bool real=True, backend="lambda", use_numpy=None):
     if backend == "llvm":
         IF HAVE_SYMENGINE_LLVM:
-            return LLVMDouble(args, *exprs, real=real)
+            return LLVMDouble(args, *exprs, real=real, use_numpy=use_numpy)
         ELSE:
             raise ValueError("""llvm backend is chosen, but symengine is not compiled
                                 with llvm support.""")
 
-    return LambdaDouble(args, *exprs, real=real)
+    return LambdaDouble(args, *exprs, real=real, use_numpy=use_numpy)
 
 
 def LambdifyCSE(args, *exprs, cse=None, concatenate=None, **kwargs):

symengine/tests/test_lambdify.py

@@ -183,9 +183,9 @@ def test_array_out_no_numpy():
     if sys.version_info[0] < 3:
         return  # requires Py3
     args, exprs, inp, check = _get_array()
-    lmb = se.Lambdify(args, exprs)
+    lmb = se.Lambdify(args, exprs, use_numpy=False)
     out1 = array.array('d', [0]*len(exprs))
-    out2 = lmb(inp, out1, use_numpy=False)
+    out2 = lmb(inp, out1)
     # Ensure buffer points to still data point:
     assert out1.buffer_info() == out2.buffer_info()
     assert out1 is out2
@@ -198,10 +198,10 @@ def test_array_out_no_numpy():
 
 def test_memview_out():
     args, exprs, inp, check = _get_array()
-    lmb = se.Lambdify(args, exprs)
-    cy_arr1 = lmb(inp, use_numpy=False)
+    lmb = se.Lambdify(args, exprs, use_numpy=False)
+    cy_arr1 = lmb(inp)
     check(cy_arr1)
-    cy_arr2 = lmb(inp, cy_arr1, use_numpy=False)
+    cy_arr2 = lmb(inp, cy_arr1)
     check(cy_arr2)
     assert cy_arr2[0] != -1
     cy_arr1[0] = -1
@@ -249,17 +249,17 @@ def _get_cse_exprs():
 def test_cse_list_input():
     args, exprs, inp, ref = _get_cse_exprs()
     lmb = se.LambdifyCSE(args, exprs, concatenate=lambda tup:
-                         tup[0]+list(tup[1]))
-    out = lmb(inp, use_numpy=False)
+                         tup[0]+list(tup[1]), use_numpy=False)
+    out = lmb(inp)
     assert allclose(out, ref)
 
 
 def test_cse_array_input():
     args, exprs, inp, ref = _get_cse_exprs()
     inp = array.array('d', inp)
     lmb = se.LambdifyCSE(args, exprs, concatenate=lambda tup:
-                         tup[0]+array.array('d', tup[1]))
-    out = lmb(inp, use_numpy=False)
+                         tup[0]+array.array('d', tup[1]), use_numpy=False)
+    out = lmb(inp)
     assert allclose(out, ref)
 
 
@@ -299,12 +299,12 @@ def test_broadcast_fortran():
         check(A[i, ...], inp[i, :])
 
 
-def _get_1_to_2by3_matrix():
+def _get_1_to_2by3_matrix(use_numpy=None):
     x = se.symbols('x')
     args = x,
     exprs = se.DenseMatrix(2, 3, [x+1, x+2, x+3,
                                   1/x, 1/(x*x), 1/(x**3.0)])
-    L = se.Lambdify(args, exprs)
+    L = se.Lambdify(args, exprs, use_numpy=use_numpy)
 
     def check(A, inp):
         X, = inp
@@ -318,9 +318,9 @@ def check(A, inp):
 
 
 def _test_2dim_Matrix(use_numpy):
-    L, check = _get_1_to_2by3_matrix()
+    L, check = _get_1_to_2by3_matrix(use_numpy=use_numpy)
     inp = [7]
-    check(L(inp, use_numpy=use_numpy), inp)
+    check(L(inp), inp)
 
 
 def test_2dim_Matrix():
@@ -335,9 +335,9 @@ def test_2dim_Matrix_numpy():
 
 
 def _test_2dim_Matrix_broadcast(use_numpy):
-    L, check = _get_1_to_2by3_matrix()
+    L, check = _get_1_to_2by3_matrix(use_numpy=use_numpy)
     inp = range(1, 5)
-    out = L(inp, use_numpy=use_numpy)
+    out = L(inp)
     for i in range(len(inp)):
         check(out[i, ...], (inp[i],))
 
@@ -443,10 +443,10 @@ def ravelled(A):
         return L
 
 
-def _get_2_to_2by2_list(real=True):
+def _get_2_to_2by2_list(real=True, use_numpy=None):
     args = x, y = se.symbols('x y')
     exprs = [[x + y*y, y*y], [x*y*y, se.sqrt(x)+y*y]]
-    L = se.Lambdify(args, exprs, real=real)
+    L = se.Lambdify(args, exprs, real=real, use_numpy=use_numpy)
 
     def check(A, inp):
         X, Y = inp
@@ -461,19 +461,19 @@ def check(A, inp):
 
 
 def test_2_to_2by2_list():
-    L, check = _get_2_to_2by2_list()
+    L, check = _get_2_to_2by2_list(use_numpy=False)
     inp = [13, 17]
-    A = L(inp, use_numpy=False)
+    A = L(inp)
     check(A, inp)
 
 
 # @pytest.mark.skipif(not HAVE_NUMPY, reason='requires numpy')
 def test_2_to_2by2_numpy():
     if not HAVE_NUMPY:  # nosetests work-around
         return
-    L, check = _get_2_to_2by2_list()
+    L, check = _get_2_to_2by2_list(use_numpy=True)
     inp = [13, 17]
-    A = L(inp, use_numpy=True)
+    A = L(inp)
     check(A, inp)
 
 
@@ -502,9 +502,9 @@ def test_unsafe_complex():
 
 def test_itertools_chain():
     args, exprs, inp, check = _get_array()
-    L = se.Lambdify(args, exprs)
+    L = se.Lambdify(args, exprs, use_numpy=False)
     inp = itertools.chain([inp[0]], (inp[1],), [inp[2]])
-    A = L(inp, use_numpy=False)
+    A = L(inp)
     check(A)