Оценка математического выражения (функция) для большого числа входных значений быстрых

Question

следующих вопросов

• Evaluating a mathematical expression in a string
  
• Equation parsing in Python
  
• Safe way to parse user-supplied mathematical formula in Python
  
• Evaluate math equations from unsafe user input in Python

и их соответствующие ответы заставляют меня думать, как я могу разобрать одно математическое выражение (в общих чертах в соответствии с этими ответами https://stackoverflow.com/a/594294/1672565), которое дает (более или менее доверенный) пользователь эффективно для входных значений от 20 до 30 тыс., Поступающих из базы данных. Я реализовал быстрый и грязный тест, чтобы сравнить различные решения.

# 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 
     }) 

# Решение # 1: Eval [да, совершенно небезопасно]

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') 

# Решение # 2a: SymPy - evalf (http://www.sympy.org)

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') 

# Решение №2b: sympy - lambdify (http://www.sympy.org)

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') 

# Решение # 2с: SymPy - lambdify с numexpr [и NumPy] (http://www.sympy.org)

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') 

# Решение # 3a: asteval [на основе аст] - шпагатом магии (http://newville.github.io/asteval/index.html)

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') 

# Решение # 3b (M Newville): аст Eval [на основе аст] - синтаксический & пробег (http://newville.github.io/asteval/index.html)

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') 

# Решение # 3в (M Newville): asteval [на основе аст] - синтаксический & бег с NumPy (http://newville.github.io/asteval/index.html)

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') 

# Решение # 4: simpleeval [на основе AST] (https://github.com/danthedeckie/simpleeval)

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') 

# Решение # 5 numexpr [и NumPy] (https://github.com/pydata/numexpr)

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') 

На моей старой тестовой машине (Python 3.4, Linux x86_64 3,11, два ядра, 1.8GHz) Я получаю следующие результаты:

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 

Что торчит это невероятная скорость Eval, хотя я не хочу использовать это в реальной жизни. Вторым лучшим решением является numexpr, что зависит от numpy - зависимости, которого я бы хотел избежать, хотя это не является жестким требованием. Следующее лучшее - simpleeval, который строится вокруг ast. aeval, другое решение, основанное на астрах, страдает тем, что мне нужно преобразовать каждое значение входного значения float в строку сначала, вокруг которой я не мог найти способ. sympy изначально был моим любимым, потому что он предлагает наиболее гибкое и, по-видимому, безопасное решение, но в итоге он оказался последним с впечатляющим расстоянием до второго до последнего решения.

Обновление 1: Существует более быстрый подход, используя sympy. См. Решение 2b. Это почти так же хорошо, как numexpr, хотя я не уверен, действительно ли sympy фактически использует его внутри.

Update 2: В SymPy реализации теперь используют sympify вместо упростить (в соответствии с рекомендациями его ведущим разработчиком, asmeurer - спасибо). Он не использует numexpr, если это явно не предлагается сделать (см. Решение 2c). Я также добавил два значительно более быстрых решения на основе asteval (спасибо M Newville).

---

Какие у меня есть возможности для ускорения любых относительно безопасных решений? Существуют ли другие, безопасные (-ий) подходы, например, с использованием ast?

Answer 1

Поскольку вы спросили о asteval, там это способ использовать его и получить более быстрые результаты:

aeval = Interpreter() 
time_start = time.time() 
expr = aeval.parse(userinput_function) 
for item in database_xy: 
    aeval.symtable['x'] = item['x'] 
    item['y_aeval'] = aeval.run(expr) 
time_end = time.time() 

То есть, вы можете сначала разобрать («предварительной компиляции») функции входа пользователя, а затем вставьте каждое новое значение x в таблицу символов и используйте Interpreter.run() для оценки скомпилированного выражения для этого значения. По вашему масштабу, я думаю, что это приблизит вас к 0,5 секундам.

Если вы готовы использовать numpy, гибридное решение:

aeval = Interpreter() 
time_start = time.time() 
expr = aeval.parse(userinput_function) 
x = numpy.array([item['x'] for item in database_xy]) 
aeval.symtable['x'] = x 
y = aeval.run(expr) 
time_end = time.time() 

должен быть гораздо быстрее, и сопоставим по времени выполнения с использованием numexpr.

Answer 2

Если вы передаете строку sympy.simplify (что не рекомендуется использование, рекомендуется использовать sympify явно), что собирается использовать sympy.sympify, чтобы преобразовать его в выражение SymPy, который использует eval внутренне.

Answer 3

CPython (и pypy) использует очень простой язык стека для выполнения функций, и довольно просто написать байт-код самостоятельно, используя модуль ast.

import sys 
PY3 = sys.version_info.major > 2 
import ast 
from ast import parse 
import types 
from dis import opmap 

ops = { 
    ast.Mult: opmap['BINARY_MULTIPLY'], 
    ast.Add: opmap['BINARY_ADD'], 
    ast.Sub: opmap['BINARY_SUBTRACT'], 
    ast.Div: opmap['BINARY_TRUE_DIVIDE'], 
    ast.Pow: opmap['BINARY_POWER'], 
} 
LOAD_CONST = opmap['LOAD_CONST'] 
RETURN_VALUE = opmap['RETURN_VALUE'] 
LOAD_FAST = opmap['LOAD_FAST'] 
def process(consts, bytecode, p, stackSize=0): 
    if isinstance(p, ast.Expr): 
     return process(consts, bytecode, p.value, stackSize) 
    if isinstance(p, ast.BinOp): 
     szl = process(consts, bytecode, p.left, stackSize) 
     szr = process(consts, bytecode, p.right, stackSize) 
     if type(p.op) in ops: 
      bytecode.append(ops[type(p.op)]) 
     else: 
      print(p.op) 
      raise Exception("unspported opcode") 
     return max(szl, szr) + stackSize + 1 
    if isinstance(p, ast.Num): 
     if p.n not in consts: 
      consts.append(p.n) 
     idx = consts.index(p.n) 
     bytecode.append(LOAD_CONST) 
     bytecode.append(idx % 256) 
     bytecode.append(idx // 256) 
     return stackSize + 1 
    if isinstance(p, ast.Name): 
     bytecode.append(LOAD_FAST) 
     bytecode.append(0) 
     bytecode.append(0) 
     return stackSize + 1 
    raise Exception("unsupported token") 

def makefunction(inp): 
    def f(x): 
     pass 

    if PY3: 
     oldcode = f.__code__ 
     kwonly = oldcode.co_kwonlyargcount 
    else: 
     oldcode = f.func_code 
    stack_size = 0 
    consts = [None] 
    bytecode = [] 
    p = ast.parse(inp).body[0] 
    stack_size = process(consts, bytecode, p, stack_size) 
    bytecode.append(RETURN_VALUE) 
    bytecode = bytes(bytearray(bytecode)) 
    consts = tuple(consts) 
    if PY3: 
     code = types.CodeType(oldcode.co_argcount, oldcode.co_kwonlyargcount, oldcode.co_nlocals, stack_size, oldcode.co_flags, bytecode, consts, oldcode.co_names, oldcode.co_varnames, oldcode.co_filename, 'f', oldcode.co_firstlineno, b'') 
     f.__code__ = code 
    else: 
     code = types.CodeType(oldcode.co_argcount, oldcode.co_nlocals, stack_size, oldcode.co_flags, bytecode, consts, oldcode.co_names, oldcode.co_varnames, oldcode.co_filename, 'f', oldcode.co_firstlineno, '') 
     f.func_code = code 
    return f 

Это имеет явное преимущество генерации по существу ту же функцию, что и eval, и весы почти точно так же, как compile + eval (compile шаг немного медленнее, чем eval «s, и eval будет предвычисления все, что он может (1+1+x компилируется в 2+x).

для сравнения, eval заканчивает свой 20k тест на 0,0125 секунды, а makefunction заканчивается в 0,014 секунд. Увеличение числа итераций 2000000, eval заканчивается за 1.23 секунд, а makefunction заканчивается в 1,32 секунды.

Интересно отметить, что pypy признает, что eval и makefunction производят по существу ту же функцию, поэтому прогрев JIT для первого ускоряет второй.

Answer 4

Я не кодер Python, поэтому я не могу предоставить код Python. Но я думаю, что могу предоставить простую схему, которая имитирует ваши зависимости и все еще работает довольно быстро.

Ключевым моментом здесь является создание чего-то, что близко к eval, не будучи eval. Итак, что вы хотите сделать, это «скомпилировать» уравнение пользователя во что-то, что можно быстро оценить. ОП продемонстрировал ряд решений.

Вот еще один, основанный на оценке уравнения как Reverse Polish.

Для обсуждения предположим, что вы можете преобразовать уравнение в RPN (обратная полировка). Это означает, что операнды прийти до операторов, например, для формулы пользователя:

 sqrt(x**2 + y**2) 

вы получите RPN эквивалент чтения слева направо:

  x 2 ** y 2 ** + sqrt 

На самом деле, мы можем рассматривать «операнды», (например, переменные и константы) в качестве операторов, которые принимают нулевые операнды. Теперь каждый в RPN является оператором.

Если рассматривать каждый элемент оператора в качестве маркеров (предположит, что единственное малое число записывается в виде «RPNelement» ниже для каждого) и хранить их в массиве «RPN», можно оценить такую ​​формулу, используя стек магазинную довольно быстро:

 stack = {}; // make the stack empty 
     do i=1,len(RPN),1 
      case RPN[i]: 
       "0": push(stack,0); 
       "1": push(stack,1); 
       "+": push(stack,pop(stack)+pop(stack));break; 
       "-": push(stack,pop(stack)-pop(stack));break; 
       "**": push(stack,power(pop(stack),pop(stack)));break; 
       "x": push(stack,x);break; 
       "y": push(stack,y);break; 
       "K1": push(stack,K1);break; 
       ... // as many K1s as you have typical constants in a formula 
      endcase 
     enddo 
     answer=pop(stack); 

Вы можете встроить операции для push и pop для ускорения бит. Если поставляемый RPN хорошо сформирован, этот код совершенно безопасен.

Теперь, как получить RPN? Ответ: создайте небольшой рекурсивный синтаксический анализатор, чьи действия присоединяют RPN-операторы к массиву RPN. См. my SO answer for how to build a recursive descent parser easily для типичных уравнений.

Вам нужно организовать, чтобы положить константы, встречающиеся при разборе в K1, K2, ... если они не являются особыми, обычно встречающимися значениями (как я показал для «0» и «1», вы можете добавьте больше, если полезно).

Это решение должно составлять не более нескольких сотен строк и иметь нулевые зависимости от других пакетов.

(эксперты Python: не стесняйтесь редактировать код, чтобы сделать его Pythonesque).

Answer 5

Я использовал библиотеку C++ ExprTK в прошлом с большим успехом. Here - это тест скорости теста среди других сиквелов C++ (например, Muparser, MathExpr, ATMSP и т. Д.), А ExprTK - сверху.

Существует оболочка Python для ExprTK, называемая cexprtk, которую я использовал и нашел очень быстро. Вы можете скомпилировать математическое выражение только один раз, а затем оценить это сериализованное выражение столько раз, сколько требуется. Вот простой пример кода с использованием cexprtk с userinput_function:

import cexprtk 
import time 

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) 

time_start = time.time() 
x = 1 

st = cexprtk.Symbol_Table({"x":x}, add_constants = True) # Setup the symbol table 
Expr = cexprtk.Expression(userinput_function, st) # Apply the symbol table to the userinput_function 

for x in range(0,demo_len,1): 
    st.variables['x'] = x # Update the symbol table with the new x value 
    Expr() # evaluate expression 
time_end = time.time() 

print('1 cexprtk: ' + str(round(time_end - time_start, 4)) + ' seconds') 

На моей машине (Linux, двухъядерный, 2.5GHz), для демо длиной 20000 это завершает в 0.0202 секунд.

Для демонстрационной длины 2 000 000 cexprtk заканчивается в 1,23 секунды.