title: Evaluating a mathematical expression (function) for a large number of input values fast
body: <p>The following questions</p>
  
  <ul>
  <li><a href="https://stackoverflow.com/questions/2371436/evaluating-a-mathematical-expression-in-a-string">Evaluating a mathematical expression in a string</a> <br></li>
  <li><a href="https://stackoverflow.com/questions/594266/equation-parsing-in-python?lq=1">Equation parsing in Python</a> <br></li>
  <li><a href="https://stackoverflow.com/questions/11951701/safe-way-to-parse-user-supplied-mathematical-formula-in-python">Safe way to parse user-supplied mathematical formula in Python</a> <br></li>
  <li><a href="https://stackoverflow.com/questions/26505420/evaluate-math-equations-from-unsafe-user-input-in-python">Evaluate math equations from unsafe user input in Python</a></li>
  </ul>
  
  <p>and their respective answers made me think how I could parse a single mathematical expression (in general terms along the lines of this answer <a href="https://stackoverflow.com/a/594294/1672565">https://stackoverflow.com/a/594294/1672565</a>) given by a (more or less trusted) user efficiently for 20k to 30k input values coming from a database. I implemented a quick and dirty benchmark so I could compare different solutions.</p>
  
  <pre><code># Runs with Python 3(.4)
  import pprint
  import time
  
  # This is what I have
  userinput_function = '5*(1-(x*0.1))' # String - numbers should be handled as floats
  demo_len = 20000 # Parameter for benchmark (20k to 30k in real life)
  print_results = False
  
  # Some database, represented by an array of dicts (simplified for this example)
  
  database_xy = []
  for a in range(1, demo_len, 1):
      database_xy.append({
          'x':float(a),
          'y_eval':0,
          'y_sympya':0,
          'y_sympyb':0,
          'y_sympyc':0,
          'y_aevala':0,
          'y_aevalb':0,
          'y_aevalc':0,
          'y_numexpr': 0,
          'y_simpleeval':0
          })
  </code></pre>
  
  <p><strong># Solution #1: eval [yep, totally unsafe]</strong></p>
  
  <pre><code>time_start = time.time()
  func = eval("lambda x: " + userinput_function)
  for item in database_xy:
      item['y_eval'] = func(item['x'])
  time_end = time.time()
  if print_results:
      pprint.pprint(database_xy)
  print('1 eval: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #2a: sympy - evalf (<a href="http://www.sympy.org" rel="noreferrer">http://www.sympy.org</a>)</strong></p>
  
  <pre><code>import sympy
  time_start = time.time()
  x = sympy.symbols('x')
  sympy_function = sympy.sympify(userinput_function)
  for item in database_xy:
      item['y_sympya'] = float(sympy_function.evalf(subs={x:item['x']}))
  time_end = time.time()
  if print_results:
      pprint.pprint(database_xy)
  print('2a sympy: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #2b: sympy - lambdify (<a href="http://www.sympy.org" rel="noreferrer">http://www.sympy.org</a>)</strong></p>
  
  <pre><code>from sympy.utilities.lambdify import lambdify
  import sympy
  import numpy
  time_start = time.time()
  sympy_functionb = sympy.sympify(userinput_function)
  func = lambdify(x, sympy_functionb, 'numpy') # returns a numpy-ready function
  xx = numpy.zeros(len(database_xy))
  for index, item in enumerate(database_xy):
      xx[index] = item['x']
  yy = func(xx)
  for index, item in enumerate(database_xy):
      item['y_sympyb'] = yy[index]
  time_end = time.time()
  if print_results:
      pprint.pprint(database_xy)
  print('2b sympy: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #2c: sympy - lambdify with numexpr [and numpy] (<a href="http://www.sympy.org" rel="noreferrer">http://www.sympy.org</a>)</strong></p>
  
  <pre><code>from sympy.utilities.lambdify import lambdify
  import sympy
  import numpy
  import numexpr
  time_start = time.time()
  sympy_functionb = sympy.sympify(userinput_function)
  func = lambdify(x, sympy_functionb, 'numexpr') # returns a numpy-ready function
  xx = numpy.zeros(len(database_xy))
  for index, item in enumerate(database_xy):
      xx[index] = item['x']
  yy = func(xx)
  for index, item in enumerate(database_xy):
      item['y_sympyc'] = yy[index]
  time_end = time.time()
  if print_results:
      pprint.pprint(database_xy)
  print('2c sympy: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #3a: asteval [based on ast] - with string magic (<a href="http://newville.github.io/asteval/index.html" rel="noreferrer">http://newville.github.io/asteval/index.html</a>)</strong></p>
  
  <pre><code>from asteval import Interpreter
  aevala = Interpreter()
  time_start = time.time()
  aevala('def func(x):\n\treturn ' + userinput_function)
  for item in database_xy:
      item['y_aevala'] = aevala('func(' + str(item['x']) + ')')
  time_end = time.time()
  if print_results:
      pprint.pprint(database_xy)
  print('3a aeval: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #3b (M Newville): asteval [based on ast] - parse &amp; run (<a href="http://newville.github.io/asteval/index.html" rel="noreferrer">http://newville.github.io/asteval/index.html</a>)</strong></p>
  
  <pre><code>from asteval import Interpreter
  aevalb = Interpreter()
  time_start = time.time()
  exprb = aevalb.parse(userinput_function)
  for item in database_xy:
      aevalb.symtable['x'] = item['x']
      item['y_aevalb'] = aevalb.run(exprb)
  time_end = time.time()
  print('3b aeval: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #3c (M Newville): asteval [based on ast] - parse &amp; run with numpy (<a href="http://newville.github.io/asteval/index.html" rel="noreferrer">http://newville.github.io/asteval/index.html</a>)</strong></p>
  
  <pre><code>from asteval import Interpreter
  import numpy
  aevalc = Interpreter()
  time_start = time.time()
  exprc = aevalc.parse(userinput_function)
  x = numpy.array([item['x'] for item in database_xy])
  aevalc.symtable['x'] = x
  y = aevalc.run(exprc)
  for index, item in enumerate(database_xy):
      item['y_aevalc'] = y[index]
  time_end = time.time()
  print('3c aeval: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #4: simpleeval [based on ast] (<a href="https://github.com/danthedeckie/simpleeval" rel="noreferrer">https://github.com/danthedeckie/simpleeval</a>)</strong></p>
  
  <pre><code>from simpleeval import simple_eval
  time_start = time.time()
  for item in database_xy:
      item['y_simpleeval'] = simple_eval(userinput_function, names={'x': item['x']})
  time_end = time.time()
  if print_results:
      pprint.pprint(database_xy)
  print('4 simpleeval: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p><strong># Solution #5 numexpr [and numpy] (<a href="https://github.com/pydata/numexpr" rel="noreferrer">https://github.com/pydata/numexpr</a>)</strong></p>
  
  <pre><code>import numpy
  import numexpr
  time_start = time.time()
  x = numpy.zeros(len(database_xy))
  for index, item in enumerate(database_xy):
      x[index] = item['x']
  y = numexpr.evaluate(userinput_function)
  for index, item in enumerate(database_xy):
      item['y_numexpr'] = y[index]
  time_end = time.time()
  if print_results:
      pprint.pprint(database_xy)
  print('5 numexpr: ' + str(round(time_end - time_start, 4)) + ' seconds')
  </code></pre>
  
  <p>On my old test machine (Python 3.4, Linux 3.11 x86_64, two cores, 1.8GHz) I get the following results:</p>
  
  <pre><code>1 eval: 0.0185 seconds
  2a sympy: 10.671 seconds
  2b sympy: 0.0315 seconds
  2c sympy: 0.0348 seconds
  3a aeval: 2.8368 seconds
  3b aeval: 0.5827 seconds
  3c aeval: 0.0246 seconds
  4 simpleeval: 1.2363 seconds
  5 numexpr: 0.0312 seconds
  </code></pre>
  
  <p>What sticks out is the incredible speed of <em>eval</em>, though I do not want to use this in real life. The second best solution seems to be <em>numexpr</em>, which depends on <em>numpy</em> - a dependency I would like to avoid, although this is not a hard requirement. The next best thing is <em>simpleeval</em>, which is build around <em>ast</em>. <em>aeval</em>, another ast-based solution, suffers from the fact that I have to convert every single float input value into a string first, around which I could not find a way. <em>sympy</em> was initially my favorite because it offers the most flexible and apparently safest solution, but it ended up being last with some impressive distance to the second to last solution. </p>
  
  <p><strong>Update 1</strong>: There is a much faster approach using <em>sympy</em>. See solution 2b. It is almost as good as <em>numexpr</em>, though I am not sure whether <em>sympy</em> is actually using it internally.</p>
  
  <p><strong>Update 2</strong>: The <em>sympy</em> implementations now use <em>sympify</em> instead of <em>simplify</em> (as recommended by its lead developer, asmeurer - thanks). It is not using <em>numexpr</em> unless it is explicitly asked to do so (see solution 2c). I also added two significantly faster solutions based on <em>asteval</em> (thanks to M Newville).</p>
  
  <hr>
  
  <p>What options do I have to speed any of the relatively safer solutions up even further? Are there other, safe(-ish) approaches using ast directly for instance?</p>