**Motivation**

Take a look at the following picture.

Given are the red, blue, and green curve. I would like to find at each point on the `x`

axis the dominating curve. This is shown as the black graph in the picture. From the properties of the red, green, and blue curve (increasing and constant after a while) this boils down to find the dominating curve on the very right hand side and then move towards the left hand side finding all intersection points and update the dominating curve.

This outlined problem should be solved `T`

times. There is one final twist in the problem. The blue, green, and red curve of the next iteration are constructed via the dominating solution from the previous iteration plus some varying parameters. As an example in the picture above: The solution is the black function. This function is used to generate a new blue, green, and red curve. Then the problem start again to find the dominating one for these new curves etc.

**Question in a nutshell**

In each iteration I start at the fixed very right hand side and evaluate all three functions to see which is the dominating one. This evaluations are taking longer and longer over iteration. *My feeling is that I don't pass optimally the old dominating function to construct the new blue, green, and red curve. Reason: I got in an earlier version a maximum recursion depth error.* Other parts of the code where the value of the current dominating function (which is essential either the green, red, or blue curve) is required are also taking longer and longer with iteration.

For 5 iterations just evaluating the functions on one point on the very right hand side grows:

The results were produced via

```
test = A(5, 120000, 100000)
```

And then running

```
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
```

I would like to know why is this the case and if one can program it more efficiently.

**Detailed Code explanation**

I quickly summarize the most important functions. The complete code can be found further below. If there are any other questions regarding the code I'm more than happy to elaborate / clarify.

Method

`u`

: For the recurring task of generating a new batch of the green, red, and blue curve above we need the old dominating one.`u`

is the initialization to be used in the very first iteration.Method

`_function_template`

: The function generates versions of the green, blue, and red curve by using different parameters. It returns a function of a single input.Method

`eval`

: This is the core function to generate the blue, green, and red versions every time. It takes three varying parameters each iteration:`vfunction`

which is the dominating function from the previous step,`m`

, and`s`

which are two parameters (flaots) affecting the shape of the resulting curve. The other parameters are the same in each iteration. In the code there are sample values for`m`

and`s`

for each iteration. For the more geeky ones: It's to approximate an integral where`m`

and`s`

are the expected mean and standard deviation of the underlying normal distribution. The approximation is done via Gauss-Hermite nodes / weights.Method

`find_all_intersections`

: This is the core method finding in each iteration the dominating one. It constructs a dominating function via a piece wise concatenation of the blue, green, and red curve. This is achieved via the function`piecewise`

.

Here is the complete code

```
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
```

`eval`

/`vfunction`

increases in complexity on each iteration, as if requires re-evaluation of all the underlying and preceding functions.6more comments