Evaluar una expresión matemática (función) para un gran número de valores de entrada rápidamente

Las siguientes preguntas

  • Evaluando una expresión matemática en una cadena
  • Ecuación de análisis en Python
  • Forma segura de analizar la fórmula matemática proporcionada por el usuario en Python
  • Evaluar ecuaciones matemáticas de entrada de usuario insegura en Python

y sus respectivas respuestas me hicieron pensar cómo podría analizar una sola expresión matemática (en términos generales, en la línea de esta respuesta https://stackoverflow.com/a/594294/1672565 ) dada por un usuario (más o menos confiable) de manera eficiente para valores de entrada de 20k a 30k que provienen de una base de datos. Implementé un punto de referencia rápido y sucio para poder comparar diferentes soluciones.

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

# Solución # 1: eval [sí, totalmente inseguro]

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

# Solución # 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') 

# Solución # 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') 

# Solución # 2c: sympy – lambdify con numexpr [y 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') 

# Solución # 3a: asteval [basado en ast] – con cadena mágica ( 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') 

# Solución # 3b (M Newville): asteval [basado en ast] – parse & run ( 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') 

# Solución # 3c (M Newville): asteval [basado en ast]: analizar y ejecutar con 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') 

# Solución # 4: simpleeval [basado en 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') 

# Solución # 5 numexpr [y 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') 

En mi antigua máquina de prueba (Python 3.4, Linux 3.11 x86_64, dos núcleos, 1.8GHz) obtengo los siguientes resultados:

 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 

Lo que sobresale es la increíble velocidad de evaluación , aunque no quiero usar esto en la vida real. La segunda mejor solución parece ser numexpr , que depende de numpy , una dependencia que me gustaría evitar, aunque esto no es un requisito difícil. La siguiente mejor cosa es simpleeval , que se construye alrededor de ast . Aeval , otra solución basada en ast, sufre el hecho de que primero tengo que convertir cada valor de entrada flotante en una cadena, alrededor de la cual no pude encontrar una manera. Sympy fue inicialmente mi favorita porque ofrece la solución más flexible y aparentemente más segura, pero terminó siendo la última con una distancia impresionante de la segunda a la última solución.

Actualización 1 : Hay un enfoque mucho más rápido usando sympy . Ver solución 2b. Es casi tan bueno como numexpr , aunque no estoy seguro de si Sympy realmente lo usa internamente.

Actualización 2 : las implementaciones de Sympy ahora usan sympify en lugar de simplificar (según lo recomendado por su desarrollador líder, asmeurer – gracias). No está usando numexpr a menos que se le pida explícitamente que lo haga (vea la solución 2c). También agregué dos soluciones significativamente más rápidas basadas en asteval (gracias a M Newville).


¿Qué opciones tengo para acelerar aún más cualquiera de las soluciones relativamente más seguras? ¿Hay otros enfoques seguros (-ish) usando ast directamente, por ejemplo?

Dado que ha preguntado acerca de la evolución, hay una manera de usarlo y obtener resultados más rápidos:

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

Es decir, primero puede analizar (“precomstackr”) la función de entrada del usuario y luego insertar cada nuevo valor de x en la tabla de símbolos y usar Interpreter.run() para evaluar la expresión comstackda para ese valor. En tu escala, creo que esto te llevará a cerca de 0,5 segundos.

Si estás dispuesto a usar numpy , una solución híbrida:

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

debería ser mucho más rápido y comparable en tiempo de ejecución al uso de numexpr .

Si le está pasando una cadena a sympy.simplify (que no se recomienda su uso; se recomienda usar sympify explícitamente), eso va a usar sympy.sympify para convertirla en una expresión SymPy, que usa eval internamente.

CPython (y pypy) utilizan un lenguaje de stack muy simple para ejecutar funciones, y es bastante fácil escribir el código de bytes por ti mismo, usando el módulo 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 pn not in consts: consts.append(pn) idx = consts.index(pn) 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 

Esto tiene la clara ventaja de generar esencialmente la misma función que eval , y se escala casi exactamente igual que compile + eval (el paso de compile es ligeramente más lento que el de eval , y eval precomputa todo lo que pueda ( 1+1+x se comstack como 2+x ).

A modo de comparación, eval finaliza su prueba de 20k en 0.0125 segundos y realiza la makefunction en 0.014 segundos. Aumentando el número de iteraciones a 2,000,000, eval finaliza en 1.23 segundos y hace que la makefunction termine en 1.32 segundos.

Una nota interesante, pypy reconoce que eval y makefunction producen esencialmente la misma función, por lo que el calentamiento de JIT para el primero acelera el segundo.

No soy un codificador de Python, por lo que no puedo proporcionar el código de Python. Pero creo que puedo proporcionar un esquema simple que minimiza sus dependencias y aún se ejecuta bastante rápido.

La clave aquí es construir algo que esté cerca de eval sin ser eval. Entonces, lo que quiere hacer es “comstackr” la ecuación del usuario en algo que se pueda evaluar rápidamente. OP ha mostrado una serie de soluciones.

Aquí hay otro basado en la evaluación de la ecuación como pulido inverso .

En aras de la discusión, suponga que puede convertir la ecuación en RPN (notación de pulido inverso). Esto significa que los operandos vienen antes que los operadores, por ejemplo, para la fórmula del usuario:

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

obtienes lectura RPN equivalente de izquierda a derecha:

  x 2 ** y 2 ** + sqrt 

De hecho, podemos tratar los “operandos” (por ejemplo, variables y constantes) como operadores que toman cero operandos. Ahora todo en RPN es un operador.

Si tratamos cada elemento de operador como un token (supongamos un pequeño entero único escrito como ” RPNelemento ” para cada uno) y los almacenamos en una matriz “RPN”, podemos evaluar dicha fórmula usando una stack de pushdown bastante rápido:

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

Puede alinear las operaciones de push y pop para acelerarlo un poco. Si el RPN suministrado está bien formado, este código es perfectamente seguro.

Ahora, ¿cómo conseguir el RPN? Respuesta: cree un pequeño analizador de descenso recursivo, cuyas acciones agreguen operadores RPN a la matriz RPN. Vea mi respuesta SO para saber cómo construir fácilmente un analizador de descenso recursivo para las ecuaciones típicas.

Tendrá que organizar para poner las constantes encontradas en el análisis en K1, K2, … si no son valores especiales que ocurren comúnmente (como he mostrado para “0” y “1”; puede agregar más si es útil ).

Esta solución debe ser de unos cientos de líneas a lo sumo, y tiene cero dependencias en otros paquetes.

(Expertos en Python: siéntase libre de editar el código para hacerlo Pythonesque).

He usado la biblioteca ExprTK de C ++ en el pasado con gran éxito. Aquí hay una prueba de velocidad de referencia entre otros analizadores de C ++ (por ejemplo, Muparser, MathExpr, ATMSP, etc.) y ExprTK sale a la cabeza.

Hay un envoltorio de Python para ExprTK llamado cexprtk que he usado y que he encontrado que es muy rápido. Puede comstackr la expresión matemática solo una vez y luego evaluar esta expresión serializada tantas veces como sea necesario. Aquí hay un código de ejemplo simple que usa cexprtk con 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') 

En mi máquina (Linux, dual core, 2.5GHz), para una demo de 20000 esto se completa en 0.0202 segundos.

Para una longitud de demostración de 2,000,000 cexprtk finaliza en 1.23 segundos.