Benchmark llvm backend

Author: bjodah

benchmarks/Lambdify.py

@@ -1,7 +1,9 @@
+#!/usr/bin/env python
 from time import clock
 import numpy as np
 import sympy as sp
 import symengine as se
+import warnings
 
 # Real-life example (ion speciation problem in water chemistry)
 
@@ -10,22 +12,63 @@
 args = np.concatenate((x, p))
 exp = sp.exp
 exprs = [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, p[0] - p[1] + 2*p[10] + 2*p[11] - p[12] - 2*p[13] + p[2] + 2*p[5] + 2*p[6] + 2*p[7] + 2*p[8] + 2*p[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]), -p[0] - p[1] - 15*p[10] - 2*p[11] - 3*p[12] - 4*p[13] - 4*p[2] - 3*p[3] - 2*p[4] - 3*p[6] - 6*p[7] - 9*p[8] - 12*p[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*p[10] - p[2] - p[3] - p[6] - 2*p[7] - 3*p[8] - 4*p[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]), -p[1] - 2*p[11] - 3*p[12] - 4*p[13] - p[4] + exp(x[1]) + 2*exp(x[11]) + 3*exp(x[12]) + 4*exp(x[13]) + exp(x[4]), -p[10] - 2*p[11] - p[12] - p[13] - p[5] - p[6] - p[7] - p[8] - p[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])]
-lmb_symengine = se.Lambdify(args, exprs)
-lmb_sympy = sp.lambdify(args, exprs)
+lmb_sp = sp.lambdify(args, exprs)
+lbm_se = se.Lambdify(args, exprs)
+lbm_se_llvm = se.Lambdify(args, exprs, backend='llvm')
 
 inp = np.ones(28)
 
-tim_symengine = clock()
-res_symengine = np.empty(len(exprs))
-for i in range(500):
-    res_symengine = lmb_symengine(inp)
-    #lmb_symengine.unsafe_real_real(inp, res_symengine)
-tim_symengine = clock() - tim_symengine
-
 tim_sympy = clock()
 for i in range(500):
-    res_sympy = lmb_sympy(*inp)
+    res_sympy = lmb_sp(*inp)
 tim_sympy = clock() - tim_sympy
 
-print('symengine speed-up factor (higher is better) vs sympy: %12.5g' %
-      (tim_sympy/tim_symengine))
+tim_se = clock()
+res_se = np.empty(len(exprs))
+for i in range(500):
+    res_se = lbm_se(inp)
+tim_se = clock() - tim_se
+
+tim_se_llvm = clock()
+res_se_llvm = np.empty(len(exprs))
+for i in range(500):
+    res_se_llvm = lbm_se_llvm(inp)
+tim_se_llvm = clock() - tim_se_llvm
+
+
+print('SymEngine (lambda double) speed-up factor (higher is better) vs sympy: %12.5g' %
+      (tim_sympy/tim_se))
+
+print('symengine (LLVM)          speed-up factor (higher is better) vs sympy: %12.5g' %
+      (tim_sympy/tim_se_llvm))
+
+import itertools
+from functools import reduce
+from operator import mul
+
+def ManualLLVM(inputs, *outputs):
+    outputs_ravel = list(itertools.chain(*outputs))
+    cb = se.Lambdify(inputs, outputs_ravel, backend="llvm")
+    def func(*args):
+        result = []
+        n = np.empty(len(outputs_ravel))
+        t = cb.unsafe_real(np.concatenate([arg.ravel() for arg in args]), n)
+        start = 0
+        for output in outputs:
+            elems = reduce(mul, output.shape)
+            result.append(n[start:start+elems].reshape(output.shape))
+            start += elems
+        return result
+    return func
+
+lbm_se_llvm_manual = ManualLLVM(args, np.array(exprs))
+tim_se_llvm_manual = clock()
+res_se_llvm_manual = np.empty(len(exprs))
+for i in range(500):
+    res_se_llvm_manual = lbm_se_llvm_manual(inp)
+tim_se_llvm_manual = clock() - tim_se_llvm_manual
+print('symengine (ManualLLVM)    speed-up factor (higher is better) vs sympy: %12.5g' %
+      (tim_sympy/tim_se_llvm_manual))
+
+if tim_se_llvm_manual < tim_se_llvm:
+    warnings.warn("Cython code for Lambdify.__call__ is slow.")

symengine/lib/symengine_wrapper.pyx

@@ -3386,7 +3386,7 @@ def Lambdify(args, *exprs, bool real=True, backend="lambda"):
     return LambdaDouble(args, *exprs, real=real)
 
 
-def LambdifyCSE(args, *exprs, real=True, cse=None, concatenate=None):
+def LambdifyCSE(args, *exprs, cse=None, concatenate=None, **kwargs):
     """
     Analogous with Lambdify but performs common subexpression elimination
     internally. See docstring of Lambdify.
@@ -3395,13 +3395,14 @@ def LambdifyCSE(args, *exprs, real=True, cse=None, concatenate=None):
     ----------
     args: iterable of symbols
     exprs: iterable of expressions (with symbols from args)
-    real: bool (default: True)
     cse: callback (default: None)
         defaults to sympy.cse (see SymPy documentation)
     concatenate: callback (default: numpy.concatenate)
         Examples when not using numpy:
         ``lambda tup: tup[0]+list(tup[1])``
         ``lambda tup: tup[0]+array.array('d', tup[1])``
+    \*\*kwargs: Keyword arguments passed onto Lambdify
+
     """
     if cse is None:
         from sympy import cse
@@ -3424,16 +3425,16 @@ def LambdifyCSE(args, *exprs, real=True, cse=None, concatenate=None):
             else:
                 raise NotImplementedError("n-dimensional output not yet supported.")
             n_taken += size
-        lmb = Lambdify(tuple(args) + cse_symbs, *new_exprs, real=real)
-        cse_lambda = Lambdify(args, cse_exprs, real=real)
-        def cb(inp, out=None, **kwargs):
-            cse_vals = cse_lambda(inp, **kwargs)
+        lmb = Lambdify(tuple(args) + cse_symbs, *new_exprs, **kwargs)
+        cse_lambda = Lambdify(args, cse_exprs, **kwargs)
+        def cb(inp, out=None, **kw):
+            cse_vals = cse_lambda(inp, **kw)
             new_inp = concatenate((inp, cse_vals))
-            return lmb(new_inp, out, **kwargs)
+            return lmb(new_inp, out, **kw)
 
         return cb
     else:
-        return Lambdify(args, *exprs, real=real)
+        return Lambdify(args, *exprs, **kwargs)
 
 
 def has_symbol(obj, symbol=None):