python - ¿Cómo pasar la función de manera eficiente?

algorithm python-2.7

Question

Motivación

Echa un vistazo a la siguiente imagen.

Se dan las curvas roja, azul y verde. Me gustaría encontrar en cada punto del eje x la curva dominante. Esto se muestra como el gráfico negro en la imagen. A partir de las propiedades de la curva roja, verde y azul (que aumenta y se mantiene constante después de un tiempo), esto se reduce a encontrar la curva dominante en el lado derecho y luego moverse hacia el lado izquierdo para encontrar todos los puntos de intersección y actualizar el dominio curva.

Este problema esbozado debe ser resuelto T veces. Hay un giro final en el problema. La curva azul, verde y roja de la siguiente iteración se construye a través de la solución dominante de la iteración anterior más algunos parámetros variables. Como un ejemplo en la imagen de arriba: la solución es la función de negro. Esta función se utiliza para generar una nueva curva azul, verde y roja. Luego, el problema comienza de nuevo para encontrar el dominante para estas nuevas curvas, etc.

Pregunta en pocas palabras
En cada iteración comienzo en el lado derecho fijo y evalúo las tres funciones para ver cuál es la dominante. Estas evaluaciones son cada vez más largas en la iteración. Mi sensación es que no paso de manera óptima la antigua función dominante para construir la nueva curva azul, verde y roja. Motivo: obtuve en una versión anterior un error máximo de profundidad de recursión. Otras partes del código donde se requiere el valor de la función de dominio actual (que es esencial ya sea la curva verde, roja o azul) también están tomando más tiempo y más tiempo con la iteración.

Para 5 iteraciones, solo evaluar las funciones en un punto en el lado derecho crece:

Los resultados fueron producidos vía

test = A(5, 120000, 100000) 

Y luego corriendo

test.find_all_intersections()
>>> test.find_all_intersections()
iteration 4
to compute function values it took
0.0102479457855
iteration 3
to compute function values it took
0.0134601593018
iteration 2
to compute function values it took
0.0294270515442
iteration 1
to compute function values it took
0.109843969345
iteration 0
to compute function values it took
0.823768854141

Me gustaría saber por qué este es el caso y si uno puede programarlo de manera más eficiente.

Explicación detallada del código

Resumo rápidamente las funciones más importantes. El código completo se puede encontrar más abajo. Si hay alguna otra pregunta con respecto al código, estoy más que feliz de elaborar / aclarar.

1. Método u : Para la tarea recurrente de generar un nuevo lote de la curva verde, roja y azul de arriba, necesitamos el anterior que domina. u es la inicialización que se utilizará en la primera iteración.

2. Método _function_template : La función genera versiones de la curva verde, azul y roja utilizando diferentes parámetros. Devuelve una función de una sola entrada.

3. Método eval : esta es la función principal para generar las versiones azul, verde y roja cada vez. Cada iteración toma tres parámetros diferentes: vfunction que es la función dominante del paso anterior, m y s que son dos parámetros (puntos) que afectan la forma de la curva resultante. Los otros parámetros son los mismos en cada iteración. En el código hay valores de muestra para m y s para cada iteración. Para los más geek: es aproximar una integral donde m y s son la media esperada y la desviación estándar de la distribución normal subyacente. La aproximación se realiza a través de los nodos / pesos de Gauss-Hermite.

4. Método find_all_intersections : Este es el método central que encuentra en cada iteración la dominante. Construye una función dominante mediante una concatenación por partes de la curva azul, verde y roja. Esto se logra a través de la función por piecewise .

Aquí está el código completo

import numpy as np
import pandas as pd
from scipy.optimize import brentq
import multiprocessing as mp
import pathos as pt
import timeit
import math
class A(object):
def u(self, w):
_w = np.asarray(w).copy()
_w[_w >= 120000] = 120000
_p = np.maximum(0, 100000 - _w)
return _w - 1000*_p**2
def __init__(self, T, upper_bound, lower_bound):
self.T = T
self.upper_bound = upper_bound
self.lower_bound = lower_bound
def _function_template(self, *args):
def _f(x):
return self.evalv(x, *args)
return _f
def evalv(self, w, c, vfunction, g, m, s, gauss_weights, gauss_nodes):
_A = np.tile(1 + m + math.sqrt(2) * s * gauss_nodes, (np.size(w), 1))
_W = (_A.T * w).T
_W = gauss_weights * vfunction(np.ravel(_W)).reshape(np.size(w),
len(gauss_nodes))
evalue = g*1/math.sqrt(math.pi)*np.sum(_W, axis=1)
return c + evalue
def find_all_intersections(self):
# the hermite gauss weights and nodes for integration
# and additional paramters used for function generation
gauss = np.polynomial.hermite.hermgauss(10)
gauss_nodes = gauss[0]
gauss_weights = gauss[1]
r = np.asarray([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
1., 1., 1., 1., 1., 1., 1., 1., 1.])
m = [[0.038063407778193614, 0.08475713587463352, 0.15420895520972322],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.03836174909668277, 0.08543620707856969, 0.15548297423808233],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.038063407778193614, 0.08475713587463352, 0.15420895520972322],
[0.038063407778193614, 0.08475713587463352, 0.15420895520972322],
[0.03836174909668277, 0.08543620707856969, 0.15548297423808233],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.038063407778193614, 0.08475713587463352, 0.15420895520972322],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.03836174909668277, 0.08543620707856969, 0.15548297423808233],
[0.038063407778193614, 0.08475713587463352, 0.15420895520972322],
[0.038063407778193614, 0.08475713587463352, 0.15420895520972322],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.03836174909668277, 0.08543620707856969, 0.15548297423808233],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624],
[0.038212567720998125, 0.08509661835487026, 0.15484578903763624]]
s = [[0.01945441966324046, 0.04690600929081242, 0.200125178687699],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019529101011406914, 0.04708607140891122, 0.20089341636351565],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.01945441966324046, 0.04690600929081242, 0.200125178687699],
[0.01945441966324046, 0.04690600929081242, 0.200125178687699],
[0.019529101011406914, 0.04708607140891122, 0.20089341636351565],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.01945441966324046, 0.04690600929081242, 0.200125178687699],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019529101011406914, 0.04708607140891122, 0.20089341636351565],
[0.01945441966324046, 0.04690600929081242, 0.200125178687699],
[0.01945441966324046, 0.04690600929081242, 0.200125178687699],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019529101011406914, 0.04708607140891122, 0.20089341636351565],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142],
[0.019491796104351332, 0.04699612658674578, 0.20050966545654142]]
self.solution = []
n_cpu = mp.cpu_count()
pool = pt.multiprocessing.ProcessPool(n_cpu)
# this function is used for multiprocessing
def call_f(f, x):
return f(x)
# this function takes differences for getting cross points
def _diff(f_dom, f_other):
def h(x):
return f_dom(x) - f_other(x)
return h
# finds the root of two function
def find_roots(F, u_bound, l_bound):
try:
sol = brentq(F, a=l_bound,
b=u_bound)
if np.absolute(sol - u_bound) > 1:
return sol
else:
return l_bound
except ValueError:
return l_bound
# piecewise function
def piecewise(l_comp, l_f):
def f(x):
_ind_f = np.digitize(x, l_comp) - 1
if np.isscalar(x):
return l_f[_ind_f](x)
else:
return np.asarray([l_f[_ind_f[i]](x[i])
for i in range(0, len(x))]).ravel()
return f
_u = self.u
for t in range(self.T-1, -1, -1):
print(''iteration'' + '' '' + str(t))
l_bound, u_bound = 0.5*self.lower_bound, self.upper_bound
l_ordered_functions = []
l_roots = []
l_solution = []
# build all function variations
l_functions = [self._function_template(0, _u, r[t], m[t][i], s[t][i],
gauss_weights, gauss_nodes)
for i in range(0, len(m[t]))]
# get the best solution for the upper bound on the very
# right hand side of wealth interval
array_functions = np.asarray(l_functions)
start_time = timeit.default_timer()
functions_values = pool.map(call_f, array_functions.tolist(),
len(m[t]) * [u_bound])
elapsed = timeit.default_timer() - start_time
print(''to compute function values it took'')
print(elapsed)
ind = np.argmax(functions_values)
cross_points = len(m[t]) * [u_bound]
l_roots.insert(0, u_bound)
max_m = m[t][ind]
l_solution.insert(0, max_m)
# move from the upper bound twoards the lower bound
# and find the dominating solution by exploring all cross
# points.
test = True
while test:
l_ordered_functions.insert(0, array_functions[ind])
current_max = l_ordered_functions[0]
l_c_max = len(m[t]) * [current_max]
l_u_cross = len(m[t]) * [cross_points[ind]]
# Find new cross points on the smaller interval
diff = pool.map(_diff, l_c_max, array_functions.tolist())
cross_points = pool.map(find_roots, diff,
l_u_cross, len(m[t]) * [l_bound])
# update the solution, cross points and current
# dominating function.
ind = np.argmax(cross_points)
l_roots.insert(0, cross_points[ind])
max_m = m[t][ind]
l_solution.insert(0, max_m)
if cross_points[ind] <= l_bound:
test = False
l_ordered_functions.insert(0, l_functions[0])
l_roots.insert(0, 0)
l_roots[-1] = np.inf
l_comp = l_roots[:]
l_f = l_ordered_functions[:]
# build piecewise function which is used for next
# iteration.
_u = piecewise(l_comp, l_f)
_sol = pd.DataFrame(data=l_solution,
index=np.asarray(l_roots)[0:-1])
self.solution.insert(0, _sol)
return self.solution

Answer 1

Comencemos cambiando el código para generar la iteración actual:

_u = self.u
for t in range(0, self.T):
print(t)
lparams = np.random.randint(self.a, self.b, 6).reshape(3, 2).tolist()
functions = [self._function_template(_u, *lparams[i])
for i in range(0, 3)]
# evaluate functions
pairs = list(itertools.combinations(functions, 2))
fval = [F(diff(*pairs[i]), self.a, self.b) for i in range(0, 3)]
ind = np.sort(np.unique(np.random.randint(self.a, self.b, 10)))
_u = _temp(ind, np.asarray(functions)[ind % 3])

Mirando en la línea causando el comportamiento,

fval = [F(diff(*pairs[i]), self.a, self.b) for i in range(0, 3)]

Las funciones de interés serían F y diff . Siendo este último sencillo, el primero:

def F(f, a, b):
try:
brentq(f, a=a, b=b)
except ValueError:
pass

Hmm, tragando excepciones, veamos que pasa si nosotros:

def F(f, a, b):
brentq(f, a=a, b=b)

Inmediatamente, para la primera función y en la primera iteración, se produce un error:

ValueError: f (a) y f (b) deben tener signos diferentes

En cuanto a la docs este es un requisito previo de la función de búsqueda de raíz brentq . Cambiemos la definición una vez más para monitorear esta condición en cada iteración.

def F(f, a, b):
try:
brentq(f, a=a, b=b)
except ValueError as e:
print(e)

La salida es

i
f(a) and f(b) must have different signs
f(a) and f(b) must have different signs
f(a) and f(b) must have different signs

para i va de 0 a 57. Es decir, la primera vez que la función F hace un trabajo real es para i=58 . Y sigue haciéndolo para valores más altos de i .

Conclusión : toma más tiempo para estos valores más altos, porque:

1. La raíz nunca se calcula para los valores más bajos.
2. el número de cálculos crece linealmente para i>58

Answer 2

Su código es realmente demasiado complejo para explicar su problema: intente algo más simple. A veces tienes que escribir código solo para demostrar el problema.

Estoy tomando una puñalada, basada únicamente en su descripción en lugar de su código (aunque ejecuté el código y lo verifiqué). Aquí está tu problema:

Método eval: esta es la función principal para generar las versiones azul, verde y roja cada vez. Cada iteración toma tres parámetros diferentes: vfunción que es la función dominante del paso anterior, mys, que son dos parámetros (puntos) que afectan la forma de la curva resultante.

El parámetro vfunction es más complejo en cada iteración. Está pasando una función anidada construida sobre iteraciones anteriores, lo que provoca una ejecución recursiva. Cada iteración aumenta la profundidad de la llamada recursiva.

¿Cómo se puede evitar esto? No hay forma fácil o construida. La respuesta más simple es, suponiendo que las entradas a estas funciones sean coherentes, almacenar el resultado funcional (es decir, los números) en lugar de la función en sí. Puede hacer esto siempre que tenga un número finito de entradas conocidas.

Si las entradas a las funciones subyacentes no son consistentes, entonces no hay acceso directo. Necesitas evaluar repetidamente esas funciones subyacentes. Veo que está realizando un poco de empalme de las funciones subyacentes: puede probar si el costo de hacerlo supera el costo de simplemente aprovechar al max cada una de las funciones subyacentes.

La prueba que corrí (10 iteraciones) tomó unos segundos. No veo eso como un problema.