diff --git a/dolark/dolo_improvements.py b/dolark/dolo_improvements.py index a36b11e..d2aa78c 100644 --- a/dolark/dolo_improvements.py +++ b/dolark/dolo_improvements.py @@ -2,6 +2,11 @@ from dolo.numeric.decision_rule import CallableDecisionRule, cat_grids import numpy as np import tqdm +from dolo.numeric.misc import cartesian + + +class Vector: + pass # not really used class Linear: pass @@ -19,6 +24,89 @@ class Chebychev: 'chebychev': Chebychev() } + +from dolark.dolo_improvements import * +from typing import List, Union + +functions = { + '': [lambda u: u, lambda u: u], + 'exp': [np.exp, np.log], + 'log': [np.log, np.exp], + 'exp(exp)': [lambda x: np.exp(np.exp(x)), lambda x: np.log(np.log(x))], + 'log(log)': [lambda x: np.log(np.log(x)), lambda x: np.exp(np.exp(x))], + 'exp(exp(exp))': [lambda x: np.exp(np.exp(np.exp(x))), lambda x: np.log(np.log(np.log(x)))], + 'log(log(log))': [lambda x: np.log(np.log(np.log(x))), lambda x: np.exp(np.exp(np.exp(x)))] +} + +class WarpGrid(Grid): + + base: Grid + warp: List[str] + + def __init__(self, base:Grid, warp: Union[str, List[str]]): + if isinstance(warp, str): + warp = [str] + self.warp = warp + self.base = base + + assert(base.d == len(warp)) + + self.warp_functions = [] + self.warp_ifunctions = [] + for w in self.warp: + if w not in functions: + raise Exception("Unsupported warp function") + else: + f, i_f = functions[w] + self.warp_functions.append(f) + self.warp_ifunctions.append(i_f) + + @property + def nodes(self): + nn = self.base.nodes.copy() + for i in range(nn.shape[1]): + nn[:,i] = self.warp_functions[i](nn[:,i]) + return nn + + @property + def n_nodes(self): + return self.base.n_nodes + + def node(self, i): + nn = self.base.node(i) + return np.array([self.warp_functions[i](nn[i]) for i in range(len(self.warp))]) + + +class ICartesianGrid(Grid): + + axes: List[Vector] + + def __init__(self, axes: List[Vector]): + self.axes = tuple([np.array(e) for e in axes]) + self.d = len(axes) + self.__nodes__ = None + + @property + def n(self): + return tuple([len(e) for e in self.axes]) + + @property + def nodes(self): + if self.__nodes__ is None: + self.__nodes__ = cartesian(self.axes) + + return self.__nodes__ + + @property + def n_nodes(self): + return np.prod([len(e) for e in self.axes]) + + def node(self, i): + # TODO: improve this! + return self.nodes[i,:] + + + # we keep the same user-facing-API class CallableDecisionRule: @@ -180,7 +268,7 @@ def eval_is(itp, exo_grid: CartesianGrid, endo_grid: CartesianGrid, interp_type: @multimethod def get_coefficients(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_type: Linear, x): - return [x[i].copy() for i in range(x.shape[0])] + return [x[i].reshape(tuple(endo_grid.n)+(-1,)).copy() for i in range(x.shape[0])] @multimethod def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_type: Linear, i, s): @@ -196,6 +284,26 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty return eval_linear(gg, coeffs, s) + +# UnstructuredGrid x ICartesian x Linear + +@multimethod +def get_coefficients(itp, exo_grid: UnstructuredGrid, endo_grid: ICartesianGrid, interp_type: Linear, x): + return [x[i].reshape(tuple(endo_grid.n)+(-1,)).copy() for i in range(x.shape[0])] + +@multimethod +def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: ICartesianGrid, interp_type: Linear, i, s): + + from interpolation.splines import eval_linear + assert(s.ndim==2) + + grid = endo_grid # one single CartesianGrid + coeffs = itp.coefficients[i] + gg = endo_grid.axes + + return eval_linear(gg, coeffs, s) + + # UnstructuredGrid x Cartesian x Cubic @multimethod @@ -205,7 +313,7 @@ def get_coefficients(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, d = len(grid.n) # this gg could be stored as a member of itp gg = tuple( [(grid.min[i], grid.max[i], grid.n[i]) for i in range(d)] ) - return [prefilter_cubic(gg, x[i]) for i in range(x.shape[0])] + return [prefilter_cubic(gg, x[i].reshape(tuple(endo_grid.n)+(-1,))) for i in range(x.shape[0])] @multimethod @@ -213,7 +321,6 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty from interpolation.splines import eval_cubic assert(s.ndim==2) - grid = endo_grid # one single CartesianGrid coeffs = itp.coefficients[i] d = len(grid.n) @@ -221,6 +328,30 @@ def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: CartesianGrid, interp_ty return eval_cubic(gg, coeffs, s) + + +# deal with warped Grids +@multimethod +def get_coefficients(itp, exo_grid, endo_grid: WarpGrid, interp_type, x): + return get_coefficients(itp, exo_grid, endo_grid.base, interp_type, x) + +@multimethod +def eval_is(itp, exo_grid, endo_grid: WarpGrid, interp_type, i, s): + base = endo_grid.base + ss = s.copy() + for k,f in enumerate(endo_grid.warp_ifunctions): + ss[:,k] = f(ss[:,k]) + return eval_is(itp, exo_grid, base, interp_type, i, ss) + +@multimethod +def eval_ms(itp, exo_grid, endo_grid: WarpGrid, interp_type, m, s): + base = endo_grid.base + ss = s.copy() + for i,f in enumerate(endo_grid.warp_ifunctions): + ss[:,i] = f(ss[:,i]) + return eval_ms(itp, exo_grid, base, interp_type, m, ss) + + ###### Test if __name__ == "__main__": diff --git a/dolark/equilibrium.py b/dolark/equilibrium.py index 30b113a..aeed43f 100644 --- a/dolark/equilibrium.py +++ b/dolark/equilibrium.py @@ -13,7 +13,7 @@ def __init__(self, aggmodel, m, μ, dr, y): self.m = m self.μ = μ self.dr = dr - self.x = np.concatenate([e[None,:,:] for e in [dr(i,dr.endo_grid.nodes()) for i in range(max(dr.exo_grid.n_nodes(),1))] ], axis=0) + self.x = np.concatenate([e[None,:,:] for e in [dr(i,dr.endo_grid.nodes) for i in range(max(dr.exo_grid.n_nodes,1))] ], axis=0) self.y = y self.c = dr.coefficients @@ -24,12 +24,12 @@ def __init__(self, aggmodel, m, μ, dr, y): def as_df(self): model = self.aggmodel.model eq = self - exg = np.column_stack([range(eq.dr.exo_grid.n_nodes()), eq.dr.exo_grid.nodes()]) - edg = np.column_stack([eq.dr.endo_grid.nodes()]) + exg = np.column_stack([range(eq.dr.exo_grid.n_nodes, eq.dr.exo_grid.n_nodes)]) + edg = np.column_stack([eq.dr.endo_grid.nodes]) N_m = exg.shape[0] N_s = edg.shape[0] ssg = np.concatenate([exg[:,None,:].repeat(N_s, axis=1), edg[None,:,:].repeat(N_m, axis=0)], axis=2).reshape((N_m*N_s,-1)) - x = np.concatenate([eq.dr(i, edg) for i in range(max(eq.dr.exo_grid.n_nodes(),1))], axis=0) + x = np.concatenate([eq.dr(i, edg) for i in range(max(eq.dr.exo_grid.n_nodes,1))], axis=0) import pandas as pd cols = ['i_m'] + model.symbols['exogenous'] + model.symbols['states'] + ['μ'] + model.symbols['controls'] df = pd.DataFrame(np.column_stack([ssg, eq.μ.ravel(), x]), columns=cols) @@ -55,12 +55,12 @@ def equilibrium(hmodel, m0: 'vector', y0: 'vector', p=None, dr0=None, grids=None Π0, μ0 = ergodic_distribution(hmodel.model, dr, exg, edg, dp) - s = edg.nodes() - if exg.n_nodes()==0: + s = edg.nodes + if exg.n_nodes==0: nn = 1 μμ0 = μ0.data[None,:] else: - nn = exg.n_nodes() + nn = exg.n_nodes μμ0 = μ0.data xx0 = np.concatenate([e[None,:,:] for e in [dr(i,s) for i in range(nn)] ], axis=0) @@ -111,8 +111,8 @@ def fun(u): if verbose: print(colored("done", "green")) - # grid_m = model.exogenous.discretize(to='mc', options=[{},{'N':N_mc}]).nodes() - # grid_s = model.get_grid().nodes() + # grid_m = model.exogenous.discretize(to='mc', options=[{},{'N':N_mc}]).nodes + # grid_s = model.get_grid().nodes # y_ss = solution.x # vector of aggregate endogenous variables m_ss = m0 # vector fo aggregate exogenous diff --git a/dolark/model.py b/dolark/model.py index f378e0a..6865b6d 100644 --- a/dolark/model.py +++ b/dolark/model.py @@ -223,8 +223,8 @@ def 𝒜(self, grids, m0: 'n_e', μ0: "n_m.N" , xx0: "n_m.N.n_x", y0: "n_y", p: # this is so sad mi = self.model.calibration['exogenous'][None,:] # not used anyway... else: - mi = exg.nodes() - s = eng.nodes() + mi = exg.nodes + s = eng.nodes res = sum( [μ0[i,:] @ ℰ(mi[i,:],s,xx0[i,:,:],m0,y0,p) for i in range(xx0.shape[0]) ]) return res diff --git a/dolark/perturbation.py b/dolark/perturbation.py index a2e862e..3673145 100644 --- a/dolark/perturbation.py +++ b/dolark/perturbation.py @@ -81,8 +81,8 @@ def F(hmodel, equilibrium, states, controls, states_f, controls_f, p): dr1.set_values(x1) dr0 = time_iteration(hmodel.model, dr0=dr1, verbose=False, maxit=1, dprocess=tmc) - s = dr0.endo_grid.nodes() - n_m = _mc.n_nodes() + s = dr0.endo_grid.nodes + n_m = _mc.n_nodes xx0 = np.concatenate([e[None,:,:] for e in [dr0(i,s) for i in range(n_m)] ], axis=0) res_0 = xx0-x0 diff --git a/examples/bfs_2017.yaml b/examples/bfs_2017.yaml index 1fa1dba..24cb1ae 100644 --- a/examples/bfs_2017.yaml +++ b/examples/bfs_2017.yaml @@ -12,11 +12,11 @@ equations: - m = ((D+r)/Đ)*(m(-1)-c(-1))*(1/ψ) + ξ arbitrage: - - ( β*(c(1)/c)^(-ρ) )*(ψ(1)^(-ρ))*(D+r(1))/Đ - 1 | 0.0<=c<=m + - ( β*(c(1)/c)^(-ρ) )*(ψ(1)^(-ρ))*(D+r(1)) - 1 | 0.0<=c<=m calibration: ρ: 1 - β: 0.99 + β: 0.98 δ: 0.025 D: 1-δ Đ: 1-0.00625 diff --git a/examples/equilibrium.py b/examples/equilibrium.py index b615beb..d208fb7 100644 --- a/examples/equilibrium.py +++ b/examples/equilibrium.py @@ -55,7 +55,7 @@ eqs = find_steady_state(hmodel2) eqs # results is (for now) a list of equilibrium objects -s = eqs[0][1].dr.endo_grid.nodes().ravel() +s = eqs[0][1].dr.endo_grid.nodes.ravel() i=0 for (w,eq) in eqs: dens = eq.μ.sum(axis=0) # \mu is a 2d array: exogenous x endogenous states diff --git a/examples/reiter_example.py b/examples/reiter_example.py index 96462de..1dfb383 100644 --- a/examples/reiter_example.py +++ b/examples/reiter_example.py @@ -23,6 +23,7 @@ aggmodel.features # %% +aggmodel.agent # %% [markdown] # First we can check whether the one-agent sub-part of it works, or whether we will need to find better initial guess. @@ -39,10 +40,13 @@ eq = find_steady_state(aggmodel) eq +# %% +from matplotlib import pyplot as plt + # %% # lot's look at the aggregate equilibrium for i in range(eq.μ.shape[0]): - s = eq.dr.endo_grid.nodes() # grid for states (temporary) + s = eq.dr.endo_grid.nodes # grid for states (temporary) plt.plot(s, eq.μ[i,:]*(eq.μ[i,:].sum()), label=f"y={eq.dr.exo_grid.node(i)[2]: .2f}") plt.plot(s, eq.μ.sum(axis=0), label='total', color='black') plt.grid() diff --git a/experiments/jagman/heterogeneous_agents_example.ipynb b/experiments/jagman/heterogeneous_agents_example.ipynb index f44d873..f202d66 100644 --- a/experiments/jagman/heterogeneous_agents_example.ipynb +++ b/experiments/jagman/heterogeneous_agents_example.ipynb @@ -429,7 +429,7 @@ "source": [ "%matplotlib inline\n", "\n", - "egrid = np.unique(dr.endo_grid.nodes()[:,0])\n", + "egrid = np.unique(dr.endo_grid.nodes[:,0])\n", "kss = model.calibration['states'][1]\n", "\n", "plt.figure(figsize=(15, 7.5))\n", diff --git a/experiments/krusell_smith/krusell_smith.ipynb b/experiments/krusell_smith/krusell_smith.ipynb index 174210f..2e885d4 100644 --- a/experiments/krusell_smith/krusell_smith.ipynb +++ b/experiments/krusell_smith/krusell_smith.ipynb @@ -269,17 +269,17 @@ "source": [ "import matplotlib.pyplot as plt\n", "plt.figure()\n", - "plt.plot(sol.dr.endo_grid.nodes(),sol.dr(0,sol.dr.endo_grid.nodes()))\n", - "plt.plot(sol.dr.endo_grid.nodes(),sol.dr(1,sol.dr.endo_grid.nodes()))\n", - "plt.plot(sol.dr.endo_grid.nodes(),sol.dr.endo_grid.nodes(),'k--')\n", + "plt.plot(sol.dr.endo_grid.nodes,sol.dr(0,sol.dr.endo_grid.nodes))\n", + "plt.plot(sol.dr.endo_grid.nodes,sol.dr(1,sol.dr.endo_grid.nodes))\n", + "plt.plot(sol.dr.endo_grid.nodes,sol.dr.endo_grid.nodes,'k--')\n", "plt.xlabel(\"$a_t$\")\n", "plt.ylabel(\"$a_{t+1}$\")\n", "plt.legend([\"$e = 0$\",\"$e = 1$\",\"$a_{t+1} = a_{t}$\"])\n", "\n", "plt.figure()\n", - "plt.plot(sol.dr.endo_grid.nodes(),sol.dr(n_K*n_z*n_e-2,sol.dr.endo_grid.nodes()))\n", - "plt.plot(sol.dr.endo_grid.nodes(),sol.dr(n_K*n_z*n_e-1,sol.dr.endo_grid.nodes()))\n", - "plt.plot(sol.dr.endo_grid.nodes(),sol.dr.endo_grid.nodes(),'k--')\n", + "plt.plot(sol.dr.endo_grid.nodes,sol.dr(n_K*n_z*n_e-2,sol.dr.endo_grid.nodes))\n", + "plt.plot(sol.dr.endo_grid.nodes,sol.dr(n_K*n_z*n_e-1,sol.dr.endo_grid.nodes))\n", + "plt.plot(sol.dr.endo_grid.nodes,sol.dr.endo_grid.nodes,'k--')\n", "plt.xlabel(\"$a_t$\")\n", "plt.ylabel(\"$a_{t+1}$\")\n", "plt.legend([\"$e = 0$\",\"$e = 1$\",\"$a_{t+1} = a_{t}$\"])" diff --git a/experiments/krusell_smith/krusell_smith_wout_agg_shocks.ipynb b/experiments/krusell_smith/krusell_smith_wout_agg_shocks.ipynb index 82c4444..e3a4b5a 100644 --- a/experiments/krusell_smith/krusell_smith_wout_agg_shocks.ipynb +++ b/experiments/krusell_smith/krusell_smith_wout_agg_shocks.ipynb @@ -248,9 +248,9 @@ ], "source": [ "import matplotlib.pyplot as plt\n", - "plt.plot(sol.endo_grid.nodes(),sol(0,sol.endo_grid.nodes()))\n", - "plt.plot(sol.endo_grid.nodes(),sol(1,sol.endo_grid.nodes()))\n", - "plt.plot(sol.endo_grid.nodes(),sol.endo_grid.nodes(),'k--')\n", + "plt.plot(sol.endo_grid.nodes,sol(0,sol.endo_grid.nodes))\n", + "plt.plot(sol.endo_grid.nodes,sol(1,sol.endo_grid.nodes))\n", + "plt.plot(sol.endo_grid.nodes,sol.endo_grid.nodes,'k--')\n", "plt.xlabel(\"$a_t$\")\n", "plt.ylabel(\"$a_{t+1}$\")\n", "plt.legend([\"$z = 0$\",\"$z = 1$\",\"$a_{t+1} = a_{t}$\"])" @@ -454,9 +454,9 @@ ], "source": [ "#Policy function\n", - "plt.plot(s.endo_grid.nodes(),sol(0,s.endo_grid.nodes()))\n", - "plt.plot(s.endo_grid.nodes(),sol(1,s.endo_grid.nodes()))\n", - "plt.plot(s.endo_grid.nodes(),sol.endo_grid.nodes(),'k--')\n", + "plt.plot(s.endo_grid.nodes,sol(0,s.endo_grid.nodes))\n", + "plt.plot(s.endo_grid.nodes,sol(1,s.endo_grid.nodes))\n", + "plt.plot(s.endo_grid.nodes,sol.endo_grid.nodes,'k--')\n", "plt.xlabel(\"$a_t$\")\n", "plt.ylabel(\"$a_{t+1}$\")\n", "plt.legend([\"$z = 0$\",\"$z = 1$\",\"$a_{t+1} = a_{t}$\"])" @@ -523,7 +523,7 @@ } ], "source": [ - "plt.plot(s.endo_grid.nodes(),hist)\n", + "plt.plot(s.endo_grid.nodes,hist)\n", "plt.xlabel(\"$a$\")\n", "plt.title(\"Asymptotic density of assets\")" ] diff --git a/experiments/new_grids.py b/experiments/new_grids.py new file mode 100644 index 0000000..9c78acd --- /dev/null +++ b/experiments/new_grids.py @@ -0,0 +1,175 @@ +from dolark.dolo_improvements import multimethod, get_coefficients, eval_is, eval_ms, Linear, EmptyGrid, WarpGrid, CartesianGrid, UnstructuredGrid, DecisionRule, Cubic +import numpy as np + +g = CartesianGrid([0.0,0], [1.0,1.0], [10,10]) +wg = WarpGrid(g, ['exp','exp']) + + +from matplotlib import pyplot as plt +nn = wg.nodes + +plt.plot(nn[:,0], nn[:,1],'o') +for i in range(wg.n_nodes): + nnn = wg.node(i) + plt.plot(nnn[0], nnn[1], 'x', color='red') + + +nn = wg.nodes +exo_grid = UnstructuredGrid(np.array([[0.2, 0.5, 0.7]])) + + + +f = lambda x,y: x**2 + y +vals = f(nn[:,0], nn[:,1]).reshape((1,10,10,1)) +vals + + + +fg = CartesianGrid([1.0,1.0], [2.7,2.7], [1000,1000]) +no = fg.nodes +tvals = f(no[:,0], no[:,1]) + + + + +dr = DecisionRule(exo_grid, wg, 'linear') +dr.set_values(vals) + +s = wg.nodes +x = dr.eval_is(0,s) +assert( abs(x.ravel() - vals.ravel()).max()<1e-10 ) + +abs(dr.eval_is(0,no).ravel() - tvals.ravel()).max() + + +# cubic approximation is better + +dr = DecisionRule(exo_grid, wg, 'cubic') +dr.set_values(vals) + +s = wg.nodes +x = dr.eval_is(0,s) +assert( abs(x.ravel() - vals.ravel()).max()<1e-10 ) + +abs(dr.eval_is(0,no).ravel() - tvals.ravel()).max() + + +## let's check extrapolation properties + + +g = WarpGrid( CartesianGrid([0.0], [1.0], [10]), ['exp'] ) +f = lambda x: np.exp(x[:,0])[:,None] +nodes = g.nodes +vals = f(nodes).reshape((1,10,1)) + + +dr = DecisionRule(exo_grid, g, 'linear') +dr.set_values(vals) + + +fg = CartesianGrid([0.2], [3.0], [1000]) +nn = fg.nodes + +tvals = f(nn) +xx = dr.eval_is(0, nn) + +plt.figure(figsize=(15,10)) +plt.plot(nodes.ravel(), vals.ravel(), 'o', label='data') +plt.plot(nn.ravel(), xx.ravel(), label='true') +plt.plot(nn.ravel(), tvals.ravel(), label='warped') +plt.grid() + + + +## let's check exp(exp) + + +g = WarpGrid( CartesianGrid([-5], [0.5], [10]), ['exp(exp)'] ) +f = lambda x: (x[:,0]**2)[:,None] +nodes = g.nodes +vals = f(nodes).reshape((1,10,1)) + + +dr = DecisionRule(exo_grid, g, 'linear') +dr.set_values(vals) + +fg = CartesianGrid([0], [2], [1000]) +nn = fg.nodes + +tvals = f(nn) +xx = dr.eval_is(0, nn) + +plt.figure(figsize=(15,10)) +plt.plot(nodes.ravel(), vals.ravel(), 'o', label='data') +plt.plot(nn.ravel(), xx.ravel(), label='true') +plt.plot(nn.ravel(), tvals.ravel(), label='warped') +# plt.xscale('log') +# plt.yscale('log') +plt.ylim(0,5) +plt.xlim(0,2) +plt.grid() +# sounds like a good idea, but no thanks. + + +N = 10 + +a = np.exp(np.exp( np.linspace(-10,1, N) )) +b = np.linspace(0,1) + +from dolark.dolo_improvements import ICartesianGrid + + +g = ICartesianGrid([a]) +f = lambda x: (x[:,0]**2)[:,None] +nodes = g.nodes +vals = f(nodes).reshape((1,N,1)) + + +dr = DecisionRule(exo_grid, g, 'linear') +dr.set_values(vals) + +fg = CartesianGrid([0], [20], [1000]) +nn = fg.nodes + +tvals = f(nn) +xx = dr.eval_is(0, nn) + +plt.figure(figsize=(15,10)) +plt.plot(nodes.ravel(), vals.ravel(), 'o', label='data') +plt.plot(nn.ravel(), xx.ravel(), label='true') +plt.plot(nn.ravel(), tvals.ravel(), label='warped') +# plt.xscale('log') +# plt.yscale('log') +# plt.ylim(0,5) +# plt.xlim(0,2) +plt.grid() +# sounds like a good idea, but no thanks. + + +# 2d + +g = ICartesianGrid([a]) +f = lambda x: (x[:,0]**2)[:,None] +nodes = g.nodes +vals = f(nodes).reshape((1,N,1)) + + +dr = DecisionRule(exo_grid, g, 'linear') +dr.set_values(vals) + +fg = CartesianGrid([0], [20], [1000]) +nn = fg.nodes + +tvals = f(nn) +xx = dr.eval_is(0, nn) + +plt.figure(figsize=(15,10)) +plt.plot(nodes.ravel(), vals.ravel(), 'o', label='data') +plt.plot(nn.ravel(), xx.ravel(), label='true') +plt.plot(nn.ravel(), tvals.ravel(), label='warped') +# plt.xscale('log') +# plt.yscale('log') +# plt.ylim(0,5) +# plt.xlim(0,2) +plt.grid() +# sounds like a good idea, but no thanks. diff --git a/experiments/new_grids_discussion.ipynb b/experiments/new_grids_discussion.ipynb new file mode 100644 index 0000000..05993fb --- /dev/null +++ b/experiments/new_grids_discussion.ipynb @@ -0,0 +1,479 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Non cartesian discretizations.\n", + "\n", + "This is a notebook meant for internal communication.\n", + "We describe here two ways to define non cartesian grid:\n", + "- by interpolating on a cartesian product of increasing (but non equally spaced) points\n", + "- by performing a monotone mapping of each dimension\n", + "These are currently implemented in `dolark.dolo_improvements`. Among other things,\n", + "this file contains an experimental rewrite of the decision_rule object which could\n", + "supersede the current implementation in dolo." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib import pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "from dolo.compiler.language import language_element, eval_data\n", + "from dolo.compiler.language import eval_data\n", + "import ruamel.yaml as ry\n", + "eval_txt = lambda txt: eval_data(ry.round_trip_load(txt))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Interpolation on non cartesian grids.\n", + "\n", + "Package interpolate.py, can interpolate (multi)linearly on noncartesian grids.\n", + "We need to add these grids to the language and add relevant capability to the decision rule object." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Non cartesian grids.\n", + "\n", + "Proposed syntax is very simple. ICartesianGrid is defined in `dolo_improvements`.\n", + "\n", + "Quation: what should be the name of this kind of grids? It would feel natural to me to call them `CartesianGrid` but this clashes with some previous uses? If nobody comes up with a better name, I would try to provide two constructors `CartesianGrid(,)` and `CartesianGrid(a=, b=, n=)` which would return respectively a `ICartesianGrid` and `UCartesianGrid` (or `IGrid` and `UGrid` as in `interpolation.py`)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "from dolark.dolo_improvements import ICartesianGrid\n", + "\n", + "# IrregularCartesianGrid: argument is a list of numpy arrays.\n", + "icgrid = ICartesianGrid([ \n", + " np.exp(np.linspace(0,1,10)),\n", + " np.linspace(0,1,10)\n", + "])" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'grid': }" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Let's add it to the language\n", + "@language_element\n", + "def IrregularCartesianGrid(*axes):\n", + " return ICartesianGrid(axes)\n", + "\n", + "txt = \"\"\"\n", + "grid: !IrregularCartesianGrid\n", + " - [0.0, 0.01, 0.02, 0.05, .1]\n", + " - [0.1, 0.2]\n", + "\"\"\"\n", + "eval_txt(txt)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Complementary to the syntax above some new ways of constructing vectors could be added to the language (lists are automatically constructed to numpy arrays). For instance we could have\n", + "grid: !IrregularCartesianGrid\n", + " - \n", + " - !linspace\n", + " a: 0\n", + " b: 1\n", + " n: 10\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAARaUlEQVR4nO3df6zdd13H8efLdmMVoxtyNW6ddsSx2OkUPa2/Z4JCO0w2DJuOaWCgWSLuD2OYbECiTCBINaJxiZtGHCZkThxkRrQOYowxoD3dYLNbiqWO7bYoF8c0YxXW7e0f91s4vdztfm/POff0fu7zkdz0fL+f39/z7aun3+8596SqkCS16+tmPQFJ0nQZ9JLUOINekhpn0EtS4wx6SWrc5llPYKkXvvCFtW3btllPQ5LWlf3793++quaWKzvtgn7btm0Mh8NZT0OS1pUkn3m2Mi/dSFLjDHpJapxBL0mNM+glqXEGvSQ1rte7bpLsBn4f2AT8SVW9a0n5rwG/BBwHFoDXV9VnurLXAm/tqr69qm6f0NxP8qH7jrBn70GOPn6Mc8/ewg27LuKVLzlvGkOdFuPOeuyNxmOtaZr2+bVi0CfZBNwCvAyYB/YlubuqHhypdh8wqKonk/wy8G7g55K8APgNYAAUsL9r+4WJrYDFg3TTXQ9w7KmnATjy+DFuuusBgKn+ZZzVuLMee6PxWGua1uL86nPpZidwqKoOV9WXgTuAK0YrVNU/VNWT3ebHga3d413APVX1WBfu9wC7JzLzEXv2HvzKQTrh2FNPs2fvwUkPdVqMO+uxNxqPtaZpLc6vPkF/HvDoyPZ8t+/Z/CLwt6tpm+S6JMMkw4WFhR5TOtnRx4+tav+kzGrcWY+90XisNU1rcX71Cfoss2/ZbytJ8gssXqbZs5q2VXVbVQ2qajA3t+wneJ/TuWdvWdX+SZnVuLMee6PxWGua1uL86hP088D5I9tbgaNLKyX5KeAtwOVV9aXVtB3XDbsuYssZm07at+WMTdyw66JJD3VajDvrsTcaj7WmaS3Orz7vutkHXJjkAuAIcDVwzWiFJC8BbgV2V9XnRor2Au9Mck63/XLgprFnvcSJGxZr/a6IWY0767E3Go+1pmktzq/0+c7YJK8A3sPi2yv/tKrekeRmYFhVdyf5CPA9wGe7Jo9U1eVd29cDb+72v6Oq3vtcYw0Gg/KXmknS6iTZX1WDZctOty8HN+glafWeK+j9ZKwkNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY3rFfRJdic5mORQkhuXKb80yb1Jjie5cknZu5McSPJQkj9IkklNXpK0shWDPskm4BbgMmA78Ook25dUewS4Fnj/krY/AvwocAnw3cAO4CfGnrUkqbfNPersBA5V1WGAJHcAVwAPnqhQVQ93Zc8saVvAWcCZQIAzgP8ae9aSpN76XLo5D3h0ZHu+27eiqvoY8A/AZ7ufvVX10NJ6Sa5LMkwyXFhY6NO1JKmnPkG/3DX16tN5ku8EvgvYyuI/Di9NcunXdFZ1W1UNqmowNzfXp2tJUk99gn4eOH9keytwtGf/PwN8vKqeqKongL8Ffmh1U5QkjaNP0O8DLkxyQZIzgauBu3v2/wjwE0k2JzmDxRuxX3PpRpI0PSsGfVUdB64H9rIY0ndW1YEkNye5HCDJjiTzwFXArUkOdM0/AHwaeAD4JPDJqvrrKaxDkvQsUtXrcvuaGQwGNRwOZz0NSVpXkuyvqsFyZX4yVpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDWuV9An2Z3kYJJDSW5cpvzSJPcmOZ7kyiVl357k75M8lOTBJNsmM3VJUh8rBn2STcAtwGXAduDVSbYvqfYIcC3w/mW6eB+wp6q+C9gJfG6cCUuSVmdzjzo7gUNVdRggyR3AFcCDJypU1cNd2TOjDbt/EDZX1T1dvScmM21JUl99Lt2cBzw6sj3f7evjxcDjSe5Kcl+SPd3/EE6S5LokwyTDhYWFnl1LkvroE/RZZl/17H8z8OPAG4EdwItYvMRzcmdVt1XVoKoGc3NzPbuWJPXRJ+jngfNHtrcCR3v2Pw/cV1WHq+o48CHg+1c3RUnSOPoE/T7gwiQXJDkTuBq4u2f/+4Bzkpx4mf5SRq7tS5Kmb8Wg716JXw/sBR4C7qyqA0luTnI5QJIdSeaBq4Bbkxzo2j7N4mWbjyZ5gMXLQH88naVIkpaTqr6X29fGYDCo4XA462lI0rqSZH9VDZYr85OxktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1rlfQJ9md5GCSQ0luXKb80iT3Jjme5Mplyr8xyZEkfziJSUuS+lsx6JNsAm4BLgO2A69Osn1JtUeAa4H3P0s3vwX846lPU5J0qvq8ot8JHKqqw1X1ZeAO4IrRClX1cFXdDzyztHGSHwC+Ffj7CcxXkrRKfYL+PODRke35bt+Kknwd8LvADSvUuy7JMMlwYWGhT9eSpJ76BH2W2Vc9+38D8OGqevS5KlXVbVU1qKrB3Nxcz64lSX1s7lFnHjh/ZHsrcLRn/z8M/HiSNwDfAJyZ5Imq+pobupKk6egT9PuAC5NcABwBrgau6dN5Vf38icdJrgUGhrwkra0VL91U1XHgemAv8BBwZ1UdSHJzkssBkuxIMg9cBdya5MA0Jy1J6i9VfS+3r43BYFDD4XDW05CkdSXJ/qoaLFfmJ2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS43oFfZLdSQ4mOZTkxmXKL01yb5LjSa4c2f99ST6W5ECS+5P83CQnL0la2YpBn2QTcAtwGbAdeHWS7UuqPQJcC7x/yf4ngddU1cXAbuA9Sc4ed9KSpP4296izEzhUVYcBktwBXAE8eKJCVT3clT0z2rCqPjXy+GiSzwFzwONjz1yS1EufSzfnAY+ObM93+1YlyU7gTODTy5Rdl2SYZLiwsLDariVJz6FP0GeZfbWaQZJ8G/DnwOuq6pml5VV1W1UNqmowNze3mq4lSSvoE/TzwPkj21uBo30HSPKNwN8Ab62qj69uepKkcfUJ+n3AhUkuSHImcDVwd5/Ou/ofBN5XVX956tOUJJ2qFYO+qo4D1wN7gYeAO6vqQJKbk1wOkGRHknngKuDWJAe65j8LXApcm+QT3c/3TWUlkqRlpWpVl9unbjAY1HA4nPU0JGldSbK/qgbLlfnJWElqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGre5T6Uku4HfBzYBf1JV71pSfinwHuAS4Oqq+sBI2WuBt3abb6+q2ycx8aU+dN8R9uw9yNHHj3Hu2Vu4YddFvPIl501jqNNi3FmPvdF4rDVN0z6/Vgz6JJuAW4CXAfPAviR3V9WDI9UeAa4F3rik7QuA3wAGQAH7u7ZfmMz0F33oviPcdNcDHHvqaQCOPH6Mm+56AGCqfxlnNe6sx95oPNaaprU4v/pcutkJHKqqw1X1ZeAO4IrRClX1cFXdDzyzpO0u4J6qeqwL93uA3ROY90n27D34lYN0wrGnnmbP3oOTHuq0GHfWY280HmtN01qcX32C/jzg0ZHt+W5fH73aJrkuyTDJcGFhoWfXX3X08WOr2j8psxp31mNvNB5rTdNanF99gj7L7Kue/fdqW1W3VdWgqgZzc3M9u/6qc8/esqr9kzKrcWc99kbjsdY0rcX51Sfo54HzR7a3Akd79j9O295u2HURW87YdNK+LWds4oZdF016qNNi3FmPvdF4rDVNa3F+9XnXzT7gwiQXAEeAq4Freva/F3hnknO67ZcDN616lis4ccNird8VMatxZz32RuOx1jStxfmVqpWvwiR5BYtvn9wE/GlVvSPJzcCwqu5OsgP4IHAO8H/Af1bVxV3b1wNv7rp6R1W997nGGgwGNRwOT3lBkrQRJdlfVYNly/oE/Voy6CVp9Z4r6P1krCQ1zqCXpMYZ9JLUOINekhp32t2MTbIAfGaMLl4IfH5C01kvNtqaN9p6wTVvFOOs+TuqatlPnJ52QT+uJMNnu/Pcqo225o22XnDNG8W01uylG0lqnEEvSY1rMehvm/UEZmCjrXmjrRdc80YxlTU3d41eknSyFl/RS5JGGPSS1Lh1E/RJdic5mORQkhuXKX9ekr/oyv8lybaRspu6/QeT7FrLeY/jVNec5GVJ9id5oPvzpWs991M1zvPclX97kieSvHFp29PVmOf2JUk+luRA93yftZZzP1VjnNtnJLm9W+tDSSb+a8+npceaL01yb5LjSa5cUvbaJP/e/bx21YNX1Wn/w+KvR/408CLgTOCTwPYldd4A/FH3+GrgL7rH27v6zwMu6PrZNOs1TXnNLwHO7R5/N3Bk1uuZ9ppHyv8K+EvgjbNezxo8z5uB+4Hv7ba/eQOc29cAd3SPvx54GNg26zVNaM3bgEuA9wFXjux/AXC4+/Oc7vE5qxl/vbyiX/ELyrvt27vHHwB+Mkm6/XdU1Zeq6j+AQ11/p7tTXnNV3VdVJ77J6wBwVpLnrcmsxzPO80ySV7L4l+DAGs13EsZZ88uB+6vqkwBV9d9V9TSnv3HWXMDzk2wGtgBfBv53baY9lhXXXFUPV9X9wDNL2u4C7qmqx6rqC8A9wO7VDL5egr7Pl4x/pU5VHQf+h8VXOON8ufksjbPmUa8C7quqL01pnpN0ymtO8nzgTcDb1mCekzTO8/xioJLs7f7L/+trMN9JGGfNHwC+CHwWeAT4nap6bNoTnoBxcmjsDOvzVYKngz5fMv5sdcb5cvNZGmfNi4XJxcBvs/jKbz0YZ81vA36vqp7oXuCvF+OseTPwY8AO4Engo92XT3x0slOcuHHWvBN4GjiXxcsY/5TkI1V1eLJTnLhxcmjsDFsvr+j7fMn4V+p0/637JuCxnm1PR+OsmSRbWfx6x9dU1aenPtvJGGfNPwi8O8nDwK8Cb05y/bQnPAHjntv/WFWfr6ongQ8D3z/1GY9vnDVfA/xdVT1VVZ8D/hlYD78PZ5wcGj/DZn2ToueNjM0sXnu9gK/eyLh4SZ1f4eSbN3d2jy/m5Juxh1kfN6zGWfPZXf1XzXoda7XmJXV+k/VzM3ac5/kc4F4Wb0puBj4C/PSs1zTlNb8JeC+Lr3KfDzwIXDLrNU1izSN1/4yvvRn7H93zfU73+AWrGn/WB2AVB+oVwKdYvHP9lm7fzcDl3eOzWHy3xSHgX4EXjbR9S9fuIHDZrNcy7TUDb2XxOuYnRn6+ZdbrmfbzPNLHugn6cdcM/AKLN5//DXj3rNcy7TUD39DtP9CF/A2zXssE17yDxVfvXwT+Gzgw0vb13bE4BLxutWP7KxAkqXHr5Rq9JOkUGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcf8P6Salxr0WcGMAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "grid = eval_txt(txt)['grid']\n", + "plt.plot(grid.nodes[:, 0], grid.nodes[:,1], 'o')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The Decision Rue object is amended to interpolate on such grids.\n", + "\n", + "This is done in `dolo_improvements.py`. Other combinations of exogenous grids need to be added.\n", + "\n", + "```python\n", + "\n", + "# UnstructuredGrid x ICartesian x Linear\n", + "@multimethod\n", + "def get_coefficients(itp, exo_grid: UnstructuredGrid, endo_grid: ICartesianGrid, interp_type: Linear, x):\n", + " return [x[i].reshape(tuple(endo_grid.n)+(-1,)).copy() for i in range(x.shape[0])]\n", + "\n", + "@multimethod\n", + "def eval_is(itp, exo_grid: UnstructuredGrid, endo_grid: ICartesianGrid, interp_type: Linear, i, s):\n", + "\n", + " from interpolation.splines import eval_linear\n", + " assert(s.ndim==2)\n", + "\n", + " grid = endo_grid # one single CartesianGrid\n", + " coeffs = itp.coefficients[i]\n", + " gg = endo_grid.axes\n", + "\n", + " return eval_linear(gg, coeffs, s)\n", + "\n", + "```\n", + "\n", + "improvements ideas: \n", + "- grid should have a `.numba_repr()` method that would return a numba compatible object. In the piece of code above, it could be called instead of .axes then the multimethod would work equally for Cartesian and NonCartesian grids (it is already the case for the get_coefficients method)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAI/CAYAAADz4aFLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdbVjVZaL3/d9/LUAQTJgEExwDd46apmJWYveogQzunVmT4+yZdO7R2Xunoj1YkjrHmI11AyYZWGKplTa7Y/TSnc6eZhQzw6C5yodwUkwNr/KhRINRFHkQWOf9oi73NGkpT+d6+H7eqGstWF/Pwze/4/9n6RhjBAAAAACwy2U7AAAAAADAOAMAAAAAr8A4AwAAAAAvwDgDAAAAAC/AOAMAAAAAL8A4AwAAAAAvENSeb9alSxcTHx/fnm95Rc6fP6/w8HDbGQGJs7eHs7eHs7eHs7eHs7eL87eHs7fHW89+9+7dFcaY6Es9167jLD4+Xrt27WrPt7wihYWFGjlypO2MgMTZ28PZ28PZ28PZ28PZ28X528PZ2+OtZ+84zpHLPcdtjQAAAADgBRhnAAAAAOAFGGcAAAAA4AXa9WfOLqWhoUHHjx9XXV2dtYbOnTvro48+avf3DQ0NVffu3RUcHNzu7w0AAADAu1gfZ8ePH1enTp0UHx8vx3GsNJw7d06dOnVq1/c0xqiyslLHjx9XQkJCu743AAAAAO9j/bbGuro6XXvttdaGmS2O4+jaa6+1esUQAAAAgPewPs4kBdww+78C9e8NAAAA4Ju8Ypx5kyeeeEI5OTmXfX7jxo3av39/OxYBAAAACASMs6vEOAMAAADQFnxunG0s+Uy3Z29Twpw/6fbsbdpY8lmLv+eiRYvUu3dvjRo1SgcPHpQkrVixQrfccosGDhyocePGqaamRn/5y1/03//938rIyNCgQYN0+PDhS74OAAAAAK6WT42zjSWfae7re/XZmVoZSZ+dqdXc1/e2aKDt3r1b//Vf/6WSkhK9/vrr2rlzpyTp3nvv1c6dO/XXv/5Vffv21UsvvaRhw4Zp7NixWrRokfbs2aN/+qd/uuTrAAAAAOBq+dQ4W1RwULUNTV97rLahSYsKDjb7exYVFWnMmDHq2LGjrrnmGo0dO1aStG/fPv3whz/UTTfdpNdee02lpaWX/PorfR0AAAAAfBvr/8/Z1fj8TO1VPX6lLvWpiZMmTdLGjRs1cOBArVq1SoWFhZf82it9HQAAAAB8G5+6chYbGXZVj1+J4cOH64033lBtba3OnTunP/7xj5K+/I+pu3XrpoaGBr322msXX9+pUyedO3fu4p8v9zoAAAAAuBo+Nc4y0norLNj9tcfCgt3KSOvd7O85ePBg3XvvvRo0aJDGjRunH/7wh5KkJ598UrfddptSU1PVp0+fi6//2c9+pkWLFikxMVGHDx++7OsAAAAA4Gr41G2N9yTGSfryZ88+P1Or2MgwZaT1vvh4c2VkZGjBggXfeHzatGnfeOz222//2kfpT5s27ZKvAwAAAICr4VPjTPpyoLV0jAEAAACAt/Gp2xoBAAAAwF8xzgAAAADACzDOAAAAAMALMM4AAAAAwAswzgAAAAD4FU9TozyeBtsZV41xJqlbt26SpM8//1w/+clPLNcAAAAAaK4PS9fqvleHqPSzl2ynXDWf+yj9thQbG6v169e36Xs0NjYqKIhjBwAAAFpTRcUB5W2ZoY0NJxVtjDoGRdlOumpcOfs7n376qfr37y9JWrVqle69916NHj1avXr10mOPPXbxdVu2bFFSUpIGDx6s8ePHq7q6WpK0YMEC3XLLLerfv7/uv/9+GWMkSSNHjtSvf/1rjRgxQnl5ee3/FwMAAAD8VENDjX735ym6648/0RsXyjU5opf++NOt+qdu/2o77aoxzr7Fnj17tHbtWu3du1dr167VsWPHVFFRoaeeekpbt27VBx98oCFDhmjx4sWSpBkzZmjnzp3at2+famtr9cYbb1z8XmfOnNH27dv16KOP2vrrAAAAAH7l/Q+W66e/G6qnv/iLBrrC9frI5/XIuNcVHnGd7bRm8a776zbNkcr3tu73vO4m6Z+zm/WlKSkp6ty5syTpxhtv1JEjR3TmzBnt379ft99+uyTpwoULSkpKkiS9/fbbevrpp1VTU6O//e1v6tevn+666y5J0r/+q+8tdwAAAMAbnfh8t3K2PawtTWcUJ2nJD36pkbc9Isfl29eevGuceZkOHTpc/L3b7VZjY6OMMUpNTdXvf//7r722rq5O6enp2rVrl77//e/riSeeUF1d3cXnw8PD260bAAAA8Ef1dVVaVTBDK/9WIiNp+vcGadKPnlNomO/9fNmleNc4a+YVrvY0dOhQTZ8+XWVlZbrhhhtUU1Oj48ePKyYmRpLUpUsXVVdXa/369XzyIwAAANAKjMej7Tue1cL9q3TcLaUGRWpWSq5iY4fYTmtV3jXOfEB0dLRWrVqln//856qvr5ckPfXUU/rBD36g//iP/9BNN92k+Ph43XLLLZZLAQAAAN935EiRsgszVKzz6ilHy/unK+nmqbaz2gTjTNKJEyckSfHx8dq3b58kadKkSZo0adLF1/z9h3skJydr586d3/g+Tz31lJ566qlvPF5YWNi6wQAAAICfq6k+peUF6Xr13AGFGCmj6zD9PDVXwcEdbae1GcYZAAAAAK9hPB5tLlqgnMPrdcrtaGxIV81MfV5dovvaTmtzjDMAAAAAXuHgx39S9rvztcupV18nSM/c/JgG9b/Pdla7YZwBAAAAsKqq6qiWbp6mtbVHdI2RHo9N1b3JC+UOCrGd1q4YZwAAAACs8DQ1asO22co7XqAqRxof1kMPpOWrc2S87TQrGGcAAAAA2t2HpWuVuSNLpa4mDXY66NfDfqvePxhjO8sqxhkAAACAdlNRcUB5W2ZoY8NJRRuj7Ovv1b8Mf0KOy2U7zbqAP4EzZ85oxYoVtjMAAAAAv9bQUKPf/XmK7vrjT/TGhXJNjuilP/50q+4cuYBh9pWAP4UzZ85o5cqV33i8qanJQg0AAADgf97/YLl++ruhevqLv2igK1yvj3xej4x7XeER19lO8yoBf1vjnDlz9Mknn2jQoEEKDg5WRESEunXrpj179ujPf/6zxowZc/E/ps7JyVF1dbWeeOIJHT58WNOnT9cXX3yhjh07asWKFerTp4/lvw0AAADgPU58vls52x7WlqYzipOU1+sXumPoLK6UXUbAj7Ps7Gx9+OGH2rNnjwoLC3XnnXdq3759SkhI0KeffnrZr7v//vv1wgsvqFevXnr//feVnp6ubdu2tV84AAAA4KXq66q0essMragskZE0/XuDNOlHzyk0LMp2mlfzqnG2cMdCHfjbgVb9nn2+10ezb519xa+/9dZblZCQ8K2vqa6u1l/+8heNHz/+4mP19fXNbgQAAAD8gfF4tH3Hs1q4f5WOu6XUoEjNSslVbOwQ22k+wavGmTcIDw+/+PugoCB5PJ6Lf66rq5MkeTweRUZGas+ePe3eBwAAAHijI0eKlF2YoWKdV085Wt4/XUk3T7Wd5VO8apxdzRWu1tKpUydVV1df8rmuXbvq1KlTqqysVEREhN544w2NHj1a11xzjRISErRu3TqNHz9exhh9+OGHGjhwYDvXAwAAAHbVVJ/S8oJ0vXrugEKMNCtmmO77Ua6CgzvaTvM5XjXObLj22mt12223qX///goLC1PXrl0vPhccHKzHH39ct912mxISEr72gR+vvfaapk2bpqeeekoNDQ362c9+xjgDAABAwDAejzYXLVDO4fU65XY0NqSrZqY+ry7RfW2n+ayAH2eS9PLLL6tTp06XfO7BBx/Ugw8++I3HExIStHnz5rZOAwAAALzOwY//pOx352uXU6++TpCeufkxDep/n+0sn8c4AwAAAHBFqqqOaunmaVpbe0TXGOnx2FTdm7xQ7qAQ22l+gXEGAAAA4Ft5mhq1Ydts5R0vUJUjjQ/roQfS8tU5Mt52ml9hnAEAAAC4rA9L1ypzR5ZKXU0a7HTQ3KT56tN7rO0sv+QV48wYI8dxbGe0O2OM7QQAAADgkioqDihvywxtbDipaGOUff29+pfhT8hxuWyn+S3r4yw0NFSVlZW69tprA2qgGWNUWVmp0NBQ2ykAAADARQ0NNVrz5kzln3xXdY40udMPNCUtX+ER19lO83vWx1n37t11/PhxffHFF9Ya6urqrIyk0NBQde/evd3fFwAAALiU9z9Yruw9z6vMbXS7K1yzRyxUQvxI21kBw/o4Cw4OVkJCgtWGwsJCJSYmWm0AAAAAbDnx+W7lbHtYW5rOKE5SXq9f6I6hs7iFsZ1ZH2cAAAAA7Kivq9LqLTO0orJERlL69wZp8o+eU2hYlO20gMQ4AwAAAAKNMdr+/mJl71+l424pNShSs1JyFRs7xHZZQGOcAQAAAAHkyJEiLSzMUJHOq6ccLe+frqSbp9rOghhnAAAAQECoqT6l5QXpevXcAYUYaVbMMN33o1wFB3e0nYavMM4AAAAAP2Y8Hm0uWqCcw+t1yu1obEhXzUx9Xl2i+9pOwz9gnAEAAAB+6uDHf1L2u/O1y6lXXydIz9z8mAb1v892Fi6DcQYAAAD4maqqo1q6eZrW1h7RNUZ6PDZV9yYvlDsoxHYavgXjDAAAAPATnqZGbdg2W3nHC1TlSOPDeuiBtHx1joy3nYYrwDgDAAAA/MCHpWuVuSNLpa4mDXY6aG7SfPXpPdZ2Fq4C4wwAAADwYRUVB5S3ZYY2NpxUtDHKuv7HunP4b+W4XLbTcJUYZwAAAIAPamio0dqtj2hpebHqHGlypx9oSlq+wiOus52GZvrOceY4zsuSxkg6ZYzp/9Vj35O0VlK8pE8l/dQYc7rtMgEAAAD8XztKViqrZInK3Ea3u8I1e8RCJcSPtJ2FFrqSa52rJI3+h8fmSHrLGNNL0ltf/RkAAABAGzrx+W49+p8/1L99mKdaGeX1+oWW/eJ/M8z8xHdeOTPGvOM4Tvw/PHy3pJFf/X61pEJJs1uxCwAAAMBX6uuqtHrLDK2oLJGRlB41UJPTnldoWJTtNLSi5v7MWVdjzAlJMsaccBwnphWbAAAAAHxl+3uLlb3/FR13S6lBkZqVkqvY2CG2s9AGHGPMd7/oyytnb/zdz5ydMcZE/t3zp40xl5ztjuPcL+l+SeratevNa9asaYXs1lVdXa2IiAjbGQGJs7eHs7eHs7eHs7eHs7eL87enJWd/9uw+/fHUK3ov+ILiG41+1umfFRdzZysX+i9v/Xd/xx137DbGXHJdN/fK2UnHcbp9ddWsm6RTl3uhMWa5pOWSNGTIEDNy5MhmvmXbKSwslDd2BQLO3h7O3h7O3h7O3h7O3i7O357mnH1N9SktL0jXq+cOKMQtzeoyTPf9KFfBwR3bJtJP+eK/++aOs/+W9EtJ2V/9+odWKwIAAAACkPF4tLlogXIOr9cpt6OxITGambpUXaL72k5DO7mSj9L/vb788I8ujuMclzRfX46y/+U4zr9JOippfFtGAgAAAP7s4Md/Uva787XLqVdfJ0jPDM7QoJsm2M5CO7uST2v8+WWeSmnlFgAAACCgVFUdVX5ButbUfKprjPR4bKruTV4od1CI7TRY0NzbGgEAAAA0k6epURu2zVbe8QJVOdL4sB56IC1fnSPjbafBIsYZAAAA0I4+LF2rzB1ZKnU1abDTQXOT5qtP77G2s+AFGGcAAABAO6ioOKC8LTO0seGkoo1R1vU/1p3DfyvH5bKdBi/BOAMAAADaUENDjdZufURLy4tV50iTO/1AU9LyFR5xne00eBnGGQAAANBGPv9is376uwdV5ja63RWu2SMWKiF+pO0seCnGGQAAANDKTny+WznbHtaWpjOKk5TX6xe6Y+gsbmHEt2KcAQAAAK2kvq5Kq7fM0IrKEhlJ9zk9NPPnryk0LMp2GnwA4wwAAABoBdvfW6zs/a/ouFtKDYrUrJRcHTpUzTDDFWOcAQAAAC1w5EiRFm7PUJE5r55ytLx/upJunipJOnSo0G4cfArjDAAAAGiGmupTWl6QrlfPHVCIkWbFDNN9P8pVcHBH22nwUYwzAAAA4CoYj0ebixYo5/B6nXI7GhsSo5mpS9Uluq/tNPg4xhkAAABwhQ5+/Cdlvztfu5x69XWC9MzgDA26aYLtLPgJxhkAAADwHaqqjiq/IF1raj7VNUaaFztK45KfljsoxHYa/AjjDAAAALgMT1OjNmybrbzjBapypPFhPfRAWr46R8bbToMfYpwBAAAAl/Bh6Vpl7cjSPleTBjsdNDdpvvr0Hms7C36McQYAAAD8nYqKA8rbMkMbG04q2hhlXf9j3Tn8t3JcLttp8HOMMwAAAEBSQ0ON1m59REvLi1XnSJMjemnK6GUKj7jOdhoCBOMMAAAAAW9HyUpllSxRmdvodle4Zo9YqIT4kbazEGAYZwAAAAhYJz7frZxtD2tL0xnFScrr9QvdMXQWtzDCCsYZAAAAAk59XZVWb5mhFZUlMpLSowZqctrzCg2Lsp2GAMY4AwAAQEDZ/t5iZe9/RcfdUmpQpGal5Co2dojtLIBxBgAAgMBw5EiRFm7PUJE5rwQ5Wt4/XUk3T7WdBVzEOAMAAIBfq6k+pRUF6Vp97oBCjDQrZpju+1GugoM72k4DvoZxBgAAAL9kPB5tLlqgnMPrdcrtaGxIjGamLlWX6L6204BLYpwBAADA7xz8+E/Kfne+djn16usE6ZnBGRp00wTbWcC3YpwBAADAb1RVHVV+QbrW1HyqTkaaFztK45KfljsoxHYa8J0YZwAAAPB5nqZGbdg2W3nHC1TlSOPDeuiBtHx1joy3nQZcMcYZAAAAfNqHpWuVtSNb+1yNGux00Nyk+erTe6ztLOCqMc4AAADgkyoqDihvywxtbDipaGOUdf2Pdefw38pxuWynAc3COAMAAIBPaWio0dqtj2hpebHqHGlyRC9NGb1M4RHX2U4DWoRxBgAAAJ+xo2SlskqWqMxtNMzVUXNGPK2E+JG2s4BWwTgDAACA1ys/UaKctx5SQdNpxUnK7TVRyUMzuIURfoVxBgAAAK9VX1el1VtmaGVliTyS0qMGanLa8woNi7KdBrQ6xhkAAAC80vb3Fit7/ys67pZSgyI1KyVXsbFDbGcBbYZxBgAAAK9y5EiRFm7PUJE5rwQ5Wt4/XUk3T7WdBbQ5xhkAAAC8Qk31Ka0oSNfqcwcUYqRZMUm6LzVPwSEdbacB7YJxBgAAAKuMx6PNRQuUc3i9TrkdjQ2J0cOjnld0zI2204B2xTgDAACANQc//pOy352vXU69+jpBemZwhgbdNMF2FmAF4wwAAADtrqrqqPIL0rWm5lN1MtK82FEal/y03EEhttMAaxhnAAAAaDeepkZt2DZbeccLVOVI48N66IG0fHWOjLedBljHOAMAAEC72Fu6Tpk7MrXP1ajBTgfNTZqvPr3H2s4CvAbjDAAAAG2qouKA8rbM0MaGk4o2RlnX/1h3Dv+tHJfLdhrgVRhnAAAAaBMNDTVau/URLS0vVp0jTY7opSmjlyk84jrbaYBXYpwBAACg1e0oWamskiUqcxsNc3XUnBFPKyF+pO0swKsxzgAAANBqyk+UKOeth1TQdFpxknJ7TVTy0AxuYQSuAOMMAAAALVZfV6XVW2ZoZWWJPJLSowZqctrzCg2Lsp0G+AzGGQAAAFpk+3uLtXD/KzrmllKDIjUrJVexsUNsZwE+h3EGAACAZjlypEgLt2eoyJxXghwt75+upJun2s4CfBbjDAAAAFelpvqUVhSka/W5Awox0qyYJN03KlfBHcJtpwE+jXEGAACAK2I8Hm0uWqCcw+t1yu1obEiMHh71nKJj+tlOA/wC4wwAAADf6eDHf1L2u/O1y6lXXydIzwzO0KCbJtjOAvwK4wwAAACXVVV1VPkF6VpT86k6GWle7CiNS35a7qAQ22mA32GcAQAA4Bs8TY3asG228o4XqMqRxof10ANp+eocGW87DfBbjDMAAAB8zd7Sdcrckal9rkYlOh3066T56tN7rO0swO8xzgAAACBJqqw4qLwt07Wh4aSijVHW9T/WncN/K8flsp0GBATGGQAAQIBraKjR2q2PKL+8WLWONDmil6aMXqbwiOtspwEBhXEGAAAQwHaUrFRWyRKVuY2GuTpqzoinlRA/0nYWEJAYZwAAAAGo/ESJct56SAVNpxUnKbfXRCUPzeAWRsAixhkAAEAAqa+r0uotM7SyskQeSelRAzU57XmFhkXZTgMCHuMMAAAgQGx/b7EW7n9Fx9zSqKBIzUp5VnGxt9jOAvAVxhkAAICfO3t2n9JfnaMic14JcrS8f7qSbp5qOwvAP2CcAQAA+Kma6lNaUZCu1ecOKMRIs2KSdN+oXAV3CLedBuASGGcAAAB+xng82ly0QDmH1+uU21FqUyfNvWulomP62U4D8C0YZwAAAD5uY8lnWlRwUJ+fqdWQ75Wqw7Vr9NegBvV1gvTM4AydqYxjmAE+gHEGAADgwzaWfKa5r+9VUFO5/p9uL2tv578p3GM0OfR2PTTuebmDQlRYWGg7E8AVYJwBAAD4sEWbS5UYsUpHo/fqQ5ejAVVR+uvJf9N/XdNDjwSF2M4DcBUYZwAAAD5qb+k6dfneAu0NlW6odSu0fJyK6oZIks6dqbVcB+BqMc4AAAB8TGXFIeW9ma4NF07q2iCPbvx8sN6v+qkk18XXxEaG2QsE0CyMMwAAAB/R2FCrNVtnKr+8WLWONDmil66Pe0K/+fSEpKaLrwsLdisjrbe9UADNwjgDAADwATtKViqrZInK3EbDXB01Z8TTSogfKUlyh3a5+GmNsZFhykjrrXsS4+wGA7hqjDMAAAAvVn6iRDlvPaSCptOKk5Tba6KSh2bIcf3PLYz3JMYxxgA/wDgDAADwQvV1VVq9ZYZWVpbIIyk9aqAmpz2v0LAo22kA2gjjDAAAwMtsf2+xFu5/Rcfc0qigSM1KXqy4uFttZwFoY4wzAAAAL3HkSJEWbs9QkTmvBDl6sd80DRsyzXYWgHbCOAMAALCspvqUVhSka/W5Awox0qyYJN03KlfBHcJtpwFoR4wzAAAAS4zHo81FC5RzeL1OuR2NDYnRw6OeU3RMP9tpACxgnAEAAFhwqGyTsornaZdTr75OkJ4ZnKFBN02wnQXAIsYZAABAOzpbdVRLC9K1tuZTRRhpXuwojUt+Wu6gENtpACxjnAEAALQDT1OjNmybrbzjBapypPFhPfRAWr46R8bbTgPgJRhnAAAAbWxv6Tpl7sjUPlejEp0O+nXSfPXpPdZ2FgAvwzgDAABoI5UVh5T3Zro2XDipaGOUdf2Pdefw38pxuWynAfBCLRpnjuPMlPTvkoykvZImG2PqWiMMAADAVzU21GnN1oeVX16sWkeaHNFLU0YvU3jEdbbTAHixZo8zx3HiJD0o6UZjTK3jOP9L0s8krWqlNgAAAJ+zo2SlskqWqMxtNMzVUbOHZ6tnQrLtLAA+oKW3NQZJCnMcp0FSR0mftzwJAADA95SfKFHOWw+poOm04iTl9pqo5KEZ3MII4Io1e5wZYz5zHCdH0lFJtZK2GGO2tFoZAACAD6ivq9LqLTO0srJEHknpUQM1Oe15hYZF2U4D4GMcY0zzvtBxoiT9l6R/lXRG0jpJ640x//kPr7tf0v2S1LVr15vXrFnTouC2UF1drYiICNsZAYmzt4ezt4ezt4ezt8dfz/5I+QatOf+Wjgc5+mFDqP7lun9XRERv21nf4K/n7ws4e3u89ezvuOOO3caYIZd6riW3NY6S9Ikx5gtJchzndUnDJH1tnBljlktaLklDhgwxI0eObMFbto3CwkJ5Y1cg4Ozt4ezt4ezt4ezt8bezP3KkSAu3Z6jInFeC49KL/aZq2JB021mX5W/n70s4e3t88exbMs6OShrqOE5HfXlbY4qkXa1SBQAA4IVqqk9pRUG6Vp87oBAjzYpJ0n2jchXcIdx2GgA/0JKfOXvfcZz1kj6Q1CipRF9dIQMAAPAnxuPR5qIFyjm8XqfcjsaGxOjhUc8pOqaf7TQAfqRFn9ZojJkvaX4rtQAAAHidQ2WblFU8T7ucevV1gvTM4AwNummC7SwAfqilH6UPAADgl85WHdPSgmlaW/OpIow0L3aUxiU/LXdQiO00AH6KcQYAAPB3PE2N2rBttvKOF6jKkcaHfl8z0vIVGZVgOw2An2OcAQAAfGVv6Tpl7sjUPlejEp0O+nXSfPXpPdZ2FoAAwTgDAAABr7LikPLeTNeGCycVbYyyrv+x7hz+Wzkul+00AAGEcQYAAAJWY0Od1mx9WPnlxap1pMkRvTRl9DKFR1xnOw1AAGKcAQCAgLSzZKUyS5aozG00zNVRs4dnq2dCsu0sAAGMcQYAAAJK+YkS5bz1kAqaTitOUm6viUoemsEtjACsY5wBAICAUF9XpdVbZmhlZYk8ktKjBmpy2vMKDYuynQYAkhhnAAAgAGx/b7EW7n9Fx9zSqKBIzUperLi4W21nAcDXMM4AAIDfOnKkSAu3Z6jInFeCHL3Yb6qGDUm3nQUAl8Q4AwAAfqfm/BdasXmaVp87oBAjzYpJ0n2jchXcIdx2GgBcFuMMAAD4DePxaHPRAuUcXq9Tbkd3hcRo5qjnFB3Tz3YaAHwnxhkAAPALh8o2Kat4nnY59errBOmZwRkadNME21kAcMUYZwAAwKedrTqmpQXTtLbmU0UYaV7sKI1LflruoBDbaQBwVRhnAADAJ3maGrVx2xzlHd+sM440PvT7mpGWr8ioBNtpANAsjDMAAOBz9pauU+aOTO1zNSrR6aAXkh5X3953284CgBZhnAEAAJ9RWXFIeW+ma8OFk4o2RlnX/1h3Dv+tHJfLdhoAtBjjDAAAeL3Ghjqt2fqw8suLVetIkyN6acroZQqPuM52GgC0GsYZAADwajtLXlLmniUqc3k0zNVRs4dnq2dCsu0sAGh1jDMAAOCVyk+UKOeth1TQdFpxRsrtNVHJQzO4hRGA32KcAQAAr1JfV6XVW2ZoZWWJPJLSowZqctrzCg2Lsp0GAG2KcQYAALzG9vcWa+H+V3TMLY0KitSs5MWKi7vVdunI2BsAACAASURBVBYAtAvGGQAAsO7IkSIt3J6hInNeCXL0Yr+pGjYk3XYWALQrxhkAALCmpqZCKzZN1epzBxRipFkxSbpvVK6CO4TbTgOAdsc4AwAA7c54PNpctEA5h9frlNvRXcExmpn6nKJj+tlOAwBrGGcAAKBdHSrbpKziedrl1KuvE6ScxFlKHDDRdhYAWMc4AwAA7eJs1TEtLZimtTWfKsJI82JHaVzy03IHhdhOAwCvwDgDAABtytPUqNLjKzX/kz0640jjQ7+vGWn5ioxKsJ0GAF6FcQYAANrM3tJ1ytyRqX2uRiU6HfRC0uPq2/tu21kA4JUYZwAAoNVVVhxS3pvp2nDhpKKN0YygW3X/z1fKcblspwGA12KcAQCAVtPYUKc1Wx9Wfnmxah1pckQvTRm9TDt3HWCYAcB3YJwBAIBWsbPkJWXuWaIyl0fDXB01e3i2eiYkf/XsAattAOALGGcAAKBFyk+UKOeth1TQdFpxRsrtNVHJQzO4UgYAV4lxBgAAmqW+rkqrt8zQysoSeSSlRw3Q5LSlCg2Lsp0GAD6JcQYAAK7a9vcWa+H+V3TMLY0KitSs5MWKi7vVdhYA+DTGGQAAuGJHjhRp4fYMFZnzSpCjF/tN1bAh6bazAMAvMM4AAMB3qqmp0IpNU7X63AGFGGlWTJLuG5Wr4A7httMAwG8wzgAAwGUZj0cFRQuUc3i9Trod3RUco5mpzyk6pp/tNADwO4wzAABwSYfKNim7eJ52OvXq4wRpUeIsJQ6YaDsLAPwW4wwAAHzN2apjWlowTWtrPlWEkebFpmhc8iK5g0JspwGAX2OcAQAASZKnqVEb356jvGObdcaRxod+XzPS8hUZlWA7DQACAuMMAABob+k6Ze7I1D5XoxKdDnoh6XH17X237SwACCiMMwAAAlhlxSHlvZmuDRdOqosxyuxxj8aMWCDH5bKdBgABh3EGAEAAamyo05qtDyu/vFi1jjQpopempOUrolM322kAELAYZwAABJidJS8pc88Slbk8GubqqNnDs9UzIdl2FgAEPMYZAAABovxEiXLeekgFTacVZ6TcXhOVPDSDWxgBwEswzgAA8HMX6s5q9ZYZWlH5gTyS0qMGaHLaUoWGRdlOAwD8HcYZAAB+bPt7i7Vw/ys65pZSgiKVkbxYcXG32s4CAFwC4wwAAD905EiRFm7PUJE5rwQ5erHfVA0bkm47CwDwLRhnAAD4kZqaCq3YNFWrzx1QiJFmxSTpvlG5Cu4QbjsNAPAdGGcAAPgB4/GooPhJ5ZSt00m3o7uCYzQz9TlFx/SznQYAuEKMMwAAfNyhsk3KLp6nnU69+jhuLUrMUOKAibazAABXiXEGAICPOlt1TEsLpmltzaeKMNK8bikal7JI7qAQ22kAgGZgnAEA4GM8TY3a+PYc5R3brDOOND70+5qRlq/IqATbaQCAFmCcAQDgQ/aWrlPmjkztczUq0emgF5IeV9/ed9vOAgC0AsYZAAA+oLLikPLeTNeGCyfVxRhl9rhHY0YskONy2U4DALQSxhkAAF6ssaFOa7Y+rPzyYtU60qSIXpqSlq+ITt1spwEAWhnjDAAAL7Wz5CVl7lmiMpdHw1wdNXt4tnomJNvOAgC0EcYZAABepvxEiXLeekgFTacVZ6TcXhOVPDSDWxgBwM8xzgAA8BIX6s9pdcF0raj8QB5J6VEDNDltqULDomynAQDaAeMMAAAv8M57zyp7/8s65pZSgjorI/lZxcXdajsLANCOGGcAAFh09GixFhbO0jvmvBLk6MV+UzVsSLrtLACABYwzAAAsqKmp0IpNU7X63AEFG+nRmKGaMCpPwR3CbacBACxhnAEA0MY2lnymRQUH9fmZWnXr3EETfrBFG6r+rJNuR3cFx2hm6nOKjulnOxMAYBnjDACANrSx5DPNfX2vahua1DO0RFGR6/VCdZN6GZcWJT6mxAETbScCALwE4wwAgDa0qOCggppO6YfdVmpv57+p3mM04OSNOtL070ockGY7DwDgRRhnAAC0EU9To653luvYDXv1V5ejgVWR+vDUr/RuU1c5arSdBwDwMowzAADawN7Sdcrckal91zXqhlq3QsvHqahuyMXnYyPDLNYBALwR4wwAgFZUWXFIeW+ma8OFk+pijP6j4ygtK/uRahvMxdeEBbuVkdbbYiUAwBsxzgAAaAWNDXVau3WmlpYXqdaRJkX00pS0fEV06qYeN/zPpzXGRoYpI6237kmMs50MAPAyjDMAAFpoZ8lLytyzRGUuj5JcHTVneLZ6JiRffP6exDjGGADgOzHOAABopvITJcp56yEVNJ1WnJFyb5io5KQMOS6X7TQAgA9inAEAcJUu1J/T6oLpWlH5gTyS0qMGaHLaUoWGRdlOAwD4MMYZAABX4Z33n1V26cs65pZSgjorI/lZxcXdajsLAOAHGGcAAFyBo0eLtbAwQ++YaiXI0Yv9pmrYkHTbWQAAP8I4AwDgW9TUVGjFpqlafe6Ago30aPRQTUjNU3CHcNtpAAA/wzgDAOASjMejguInlVO2Tifdju4KjtHM1OcUHdPPdhoAwE8xzgAA+AeHyjYpu3iedjr16uO4tSgxQ4kDJtrOAgD4OcYZAABfOVt1TPkF07Sm5lNFGGletxSNS1kkd1CI7TQAQABgnAEAAp6nqVEb356jvGObddqRxod21wNpyxQZlWA7DQAQQBhnAICAtrd0nTJ3ZGqfq1GDnBC9kDRffXvfbTsLABCAWjTOHMeJlLRSUn9JRtKvjDH/uzXCAABoS5UVh5T3Zro2XDipLsYos8fdGjPiSTkul+00AECAaumVszxJm40xP3EcJ0RSx1ZoAgCgzTQ21Gnt1plaWl6kWkeaFNFLU9LyFdGpm+00AECAa/Y4cxznGknDJU2SJGPMBUkXWicLAIDWd+KLAo3/z4dU5vIoydVRc4Znq2dCsu0sAAAktezKWU9JX0h6xXGcgZJ2S3rIGHO+VcoAAGgl5SdKlPPWQypoOq1YI+XeMEHJSY9xCyMAwKs4xpjmfaHjDJH0nqTbjTHvO46TJ+msMWbeP7zufkn3S1LXrl1vXrNmTQuTW191dbUiIiJsZwQkzt4ezt4ezr79NDbVavfxfK33fCIj6W5PnJKuf0BBQZx/e+PfvV2cvz2cvT3eevZ33HHHbmPMkEs915Jxdp2k94wx8V/9+YeS5hhj7rzc1wwZMsTs2rWrWe/XlgoLCzVy5EjbGQGJs7eHs7eHs28f77z/rLJLX9Yxt5TiukYZyc/q449rOHtL+HdvF+dvD2dvj7eeveM4lx1nzb6t0RhT7jjOMcdxehtjDkpKkbS/ud8PAIDWcPRosRYWZugdU60EOXqx31QNG5IuSfr440K7cQAAfIuWflrjA5Je++qTGv+PpMktTwIA4OrV1FRoxaapWn3ugIKN9Gj0UE1IzVNwh3DbaQAAXJEWjTNjzB5Jl7wkBwBAezAejwqKn1RO2TqddDu6KzhGM1OfU3RMP9tpAABclZZeOQMAwJpDZZuUXTxPO5169XHcWpSYocQBE21nAQDQLIwzAIDPOVt1TPkF6VpT84kijDSvW4rGpSySOyjEdhoAAM3GOAMA+AxPU6P+8PZc5R7bpNOOND60ux5IW6bIqATbaQAAtBjjDADgE/aWrlPWjkztdTVqkBOiF5Lmq2/vu21nAQDQahhnAACvVllxSHlvpmvDhZPqYowye9ytMSOelONy2U4DAKBVMc4AAF6psaFOa7fO1NLyItU60qSIXpqSlq+ITt1spwEA0CYYZwAAr7Oz5CVl7lmiMpdHSa6OmjM8Wz0Tkm1nAQDQphhnAACvUX6iRDlvPaSCptOKNVLuDROUnPQYtzACAAIC4wwAYN2F+nNaXTBdKyo/kEfStKgB+lXaUoWGRdlOAwCg3TDOAABWvfP+s8oufVnH3FJKUGdlJD+ruLhbbWcBANDuGGcAACuOHi3WwsIMvWOqlSBHL/abqmFD0m1nAQBgDeMMANCuamoqtHLTNK0695GCjfRo9FBNSM1TcIdw22kAAFjFOAMAtAvj8aig+EnllK3TSbeju4JjNDP1OUXH9LOdBgCAV2CcAQDa3KGyTcounqedTr36OG4tSsxQ4oCJtrMAAPAqjDMAQJs5W3VM+QXpWlPziSKMNK9bisalLJI7KMR2GgAAXodxBgBodZ6mRv3h7bnKPbZJpx1pfGh3PZC2TJFRCbbTAADwWowzAECr2lu6Tlk7MrXX1ahBToheSJqvvr3vtp0FAIDXY5wBAFpFZcUh5b2Zrg0XTqqLMcrscbfGjHhSjstlOw0AAJ/AOAMAtEhjQ53Wbp2ppeVFqnWkSRG9NCUtXxGdutlOAwDApzDOAADNtrPkJWXuWaIyl0dJro6aMzxbPROSbWcBAOCTGGcAgKtWXl6iZ7Y+pM1NpxVrpNwbJig56TFuYQQAoAUYZwCAK3ah/pxWF0zXisoP5JE0LfIm/Wp0vkLDomynAQDg8xhnAIAr8s77zyq79GUdc0spQZ2Vkfys4uJutZ0FAIDfYJwBAL7V0aPFWliYoXdMteLl6MV+UzVsSLrtLAAA/A7jDABwSTU1FVq5eZpWnf1IwUZ6NHqoJqTmKbhDuO00AAD8EuMMAPA1xuNRQfGTyilbp5NuR3cFx2hm6nOKjulnOw0AAL/GOAMAXHSobJOyi+dpp1OvPo5bixIzlDhgou0sAAACAuMMAKCzVceUX5CuNTWfKMJI87qlaFzKIrmDQmynAQAQMBhnABDAPE2N+sPbc5V7bJNOO9L40O56IG2ZIqMSbKcBABBwGGcAEKD2lq5T1o5M7XU1apAToheS5qtv77ttZwEAELAYZwAQYCorDinvzXRtuHBSXYxRZo+7NWbEk3JcLttpAAAENMYZAASIxoY6rd06U0vLi1TrSJMiemlKWr4iOnWznQYAAMQ4A4CAsLPkJWXuWaIyl0dJro6aMzxLPRNSbGcBAIC/wzgDAD9WXr5Hz2x9UJubTivWSLk3TFBy0mPcwggAgBdinAGAH7pQf06vFszQ8srd8kiaFnmTfjU6X6FhUbbTAADAZTDOAMDPvPP+s1pY+rKOuqWUoM7KSH5WcXG32s4CAADfgXEGAH7i6NFiLSzM0DumWvFy9OKNUzTslum2swAAwBVinAGAj6upqdDKzdO06uxHCjbSo9FDNSE1T8Edwm2nAQCAq8A4AwAfZTweFRQ/qZyydTrpdnRXcIxmpj6n6Jh+ttMAAEAzMM4AwAdsLPlMiwoO6vMztYqNDNO0m09o22fPaqdTrz6OW4sSM5Q4YKLtTAAA0AKMMwDwchtLPtPc1/eqtqFJnVyV+qfQlcr57G+K8BjNixulcSmL5A4KsZ0JAABaiHEGAF5uUcFB1TXUKylqjY5Hf6g9LkcDqyJ1sjZdP/3VfbbzAABAK2GcAYCX61T/lgYk/FH7Qo1uqHUrtPxeFdXdIsd2GAAAaFWMMwDwUpUVh7Tkzen6LL5c32v0qN+JRL135l8luSRJsZFhdgMBAECrYpwBgJdpbKjT2q0ztbS8SLWOdLfrem06MkHvXeh88TVhwW5lpPW2WAkAAFob4wwAvMjOkpeUuWeJylweJbk6as7wLPVMSNGQf/i0xoy03ronMc52LgAAaEWMMwDwAuXle/TM1ge1uem0Yo2Ue8MEJSc9Jsf15S2M9yTGMcYAAPBzjDMAsOhC/Tm9WjBDyyt3yyNpWuRNmjw6X2FhUbbTAABAO2OcAYAl77z/rBaWvqyjbiklqLMykp9VXNyttrMAAIAljDMAaGdnz5Vq+qtz9Y6pVrwcvXjjFA27ZbrtLAAAYBnjDADaSU1NhVZunqZVZz9SsJEejR6qCal5Cu4QbjsNAAB4AcYZALQx4/GooPhJ5ZSt00m3o5SmTvr1XSsUE9PfdhoAAPAijDMAaEMfl21WVvFvtNOpVx/HrUWJGar6W3eGGQAA+AbGGQC0gbNVx5RfkK41NZ8owkjzuqVoXMoiuYNCVFhYaDsPAAB4IcYZALQiT1Oj/vD2XOUe26TTjjQ+tLseSFumyKgE22kAAMDLMc4AoJXsLV2nrB2Z2utq1CAnRMuGPq4b+9xjOwsAAPgIxhkAtFBlxSEteXO6Xr9Qri7GKLPH3Roz4kk5LpftNAAA4EMYZwDQTI0NdVq7daaWlhep1pEmRfTSlLR8RXTqZjsNAAD4IMYZADTDzpKXlLlnicpcHiW5OmrO8Cz1TEixnQUAAHwY4wwArkJ5+R49s/VBbW46rVgj5d4wQclJj3ELIwAAaDHGGQBcgQv15/RqwQwtr9wtj6RpkTdp8uh8hYVF2U4DAAB+gnEGAN/hnfef1cLSl3XULSUHdVbGHYvVvftttrMAAICfYZwBwGUcPVqshYUZesdUK16OXrxxiobdMt12FgAA8FOMMwD4BzU1FVq5eZpWnf1IwUZ6NHqoJqTmKbhDuO00AADgxxhnAPAV4/GooPhJ5ZSt00m3ozHB0Zo56jnFdO1vOw0AAAQAxhkASPq4bLOyin+jnU69+jhuLUrMUOKAibazAABAAGGcAQhoZ6uOaVlBun5f84kijDSvW4rGpSySOyjEdhoAAAgwjDMAAcnT1Kg/vD1Xucc26bQjjQ/trgfSlikyKsF2GgAACFCMMwABZ2/pOmXtyNReV6MGOSFaNvRx3djnHttZAAAgwDHOAASMyopDWvLmdL1+oVzXyiizx90aM+JJOS6X7TQAAADGGQD/19hQp7VbZ2ppeZFqHemXETdoatoyRXTqZjsNAADgIsYZAL+2s+QlZe1Zoo9dHiW5OmrO8Cz1TEixnQUAAPANjDMAfqm8fI+e2fqgNjedVqyRcm+YoOSkx7iFEQAAeC3GGQC/cqH+nF4tmKHllbvlkTQt8iZNTluqsI7fs50GAADwrRhnAPzGO+8/q4WlL+uoW0oO6qyMOxare/fbbGcBAABcEcYZAJ939Gixni7M0HZTrXg5evHGKRp2y3TbWQAAAFeFcQbAZ9XUVGjl5mladfYjBRvp0eihmpCap+AO4bbTAAAArhrjDIDPMR6PCoqfVE7ZOp10OxoTHK2Zo55TTNf+ttMAAACajXEGwKd8XLZZWcW/0U6nXn0ct55OnKXBA35hOwsAAKDFGGcAfMLZqmNaVpCu39d8oggj/aZbsn6SkiN3UIjtNAAAgFbBOAPg1TxNjfrD23OVe2yTTjvS+NDueiBtmSKjEmynAQAAtKoWjzPHcdySdkn6zBgzpuVJAPClvaXrlLUjU3tdjRrkhGjZ0Md1Y597bGcBAAC0ida4cvaQpI8kXdMK3wsAVFlxSEvenK7XL5TrWmOU2eNujRnxpByXy3YaAABAm2nROHMcp7ukOyX9f5IeaZUiAAGrsaFOa7fO1NLyItU60i/Db9DU0csU0amb7TQAAIA219IrZ7mSHpPUqRVaAASwnSUvKWvPEn3s8ijJ1VFzhmepZ0KK7SwAAIB24xhjmveFjjNG0r8YY9IdxxkpadalfubMcZz7Jd0vSV27dr15zZo1LchtG9XV1YqIiLCdEZA4e3u85exrzv8fbTrxogqDa9St0ejnHYcrvutP/PoWRm85+0DE2dvD2dvF+dvD2dvjrWd/xx137DbGDLnUcy0ZZ1mSfiGpUVKovvyZs9eNMRMv9zVDhgwxu3btatb7taXCwkKNHDnSdkZA4uztsX32F+rP6dWCGVpeuVseSb+KvEmT05YqrOP3rDW1F9tnH8g4e3s4e7s4f3s4e3u89ewdx7nsOGv2bY3GmLmS5n71BiP15ZWzyw4zAIFrY8lnWlRwUJ+fqVVsZJj+3z7F+u+/rdNRt5Qc1FkZdzyj7t2H2s4EAACwiv/nDECb2ljymea+vle1DU3qHnJA3a95Tc9XNaiHkV64cYpuv2WG7UQAAACv0CrjzBhTKKmwNb4XAP+yqOCg1Hhaw7u+pH1R5frEGA38oqeONkzV7bfcaTsPAADAa3DlDECbMR6P4szvFHLDDpUEuTTgbLgOnZyk4sYecmzHAQAAeBnGGYA28XHZZmUV/0YfxdYrvt6lsM/+We/WjLj4fGxkmMU6AAAA78M4A9CqzlYd07KCdP2+5hOFG+mXoUl6pexunW/4n2tlYcFuZaT1tlgJAADgfRhnAFqFp6lRf3h7rnKPbdJpR/pJaHc9kJavqKieuuEfPq0xI6237kmMs50MAADgVRhnAFpsb+k6Ze3I1F5XowY5IVo29HHd2Oeei8/fkxjHGAMAAPgOjDMAzVZZcUhL3pyu1y+U61pjlNnjbo0Z8aQcl8t2GgAAgM9hnAG4ao0NdVq7daaWlhep1pF+GX6Dpo5epohO3WynAQAA+CzGGYCrsnPPy8oqydPHLo+Gujr+/+3de3TV5Z3v8fezdwKEgAlBQG4VRIulWsSiUj1eEBQ6tZbpsee0PWqnTpeKVWeKJ6fSOm3HjnhBvHVRGEcdx5mOF1YFPXOUtE5Fim3xUmgRKAVvIDeFcDUht/2cP4g0WHCUhDw72e/XWl2QX34hn/Xt4y/7k71/z2bqWTdzzNBxqWNJkiR1eJYzSR/Kpk1LmfHMtcxv2saACHce+1XGfebbvoRRkiSpjVjOJH2g+rpdPFR1NfdufZkcMLn8RL4+YSYl3StSR5MkSepULGeSDmrh4ju5dfkDrM3CuUVlVI6dwaBBY1LHkiRJ6pQsZ5L+zNq1i7htQSXPxd0MITB7xOWcccrVqWNJkiR1apYzSfvU1GzhvvmTeXDnSoojTOlzGhefdw/FXUtTR5MkSer0LGeSiLkcVYt+yO1r5rA5G7iguA/fGv8j+vY7IXU0SZKkgmE5kwrc6jXzuWXRDbwQ6jg+ZLntpOs4eeSlqWNJkiQVHMuZVKDq697h1sc+z8M1r1Ma4Yb+53LRuNvJFnVJHU2SJKkgWc6kApNrauSJZ6dy5/qn2J4JXNRtENecP5NeFcNSR5MkSSpoljOpgLyyYg7TFk9jWaaRE3NFzD797xlx/KTUsSRJkoTlTCoIW7f8kXt+/k3m1m2kIsK0j32BHnEcI44/N3U0SZIkNbOcSZ1YY8MeHn3mW8zc9EtqA1xaeixXTpxFj579WbBgQep4kiRJasFyJnVSLy59gJuX3M3qTI4xme5MPetmjhk6LnUsSZIkHYTlTOpkNm1ayoxnrmV+0zYGRLjz2K8y7jPfJmQyqaNJkiTpA1jOpE6ivm4XD1Vdzb1bX6YpwJVlJ3LZhJmUdK9IHU2SJEkfguVM6gQWLr6TW5c/wNosnFtURuXYGQwaNCZ1LEmSJH0EljOpA1u7dhG3LajkubibIcDsEVdwxilXp44lSZKkQ2A5kzqgmpot3Dd/Mg/uXElxhCl9TuPi8+6huGtp6miSJEk6RJYzqQOJuRxVi37I7WvmsDkb+FxxH6aM/xF9+52QOpokSZJayXImdRCr18znlkU38EKo4/iQ5baTruPkkZemjiVJkqQ2YjmT8tzOHeuYVXUVD9e8TmmEG/qfy0Xjbidb1CV1NEmSJLUhy5mUp3JNjTzx7FTuWvc02wJc1G0Q15w/k14Vw1JHkyRJ0mFgOZPy0Csr5jBt8TSWZRoZGbowa8z3GHH8pNSxJEmSdBhZzqQ8Ur11NXf/7Crm1m2kIsJNH7uQC87+IZlMNnU0SZIkHWaWMymBeUvWM71qFRu21zKgvITrxg+h9p3pzNz0S2oDXFp6LFdOnEWPnv1TR5UkSVI7sZxJ7WzekvVMfXwZtQ1NAPSqf4r7f1vF2q4wJtOdqWfdzDFDxyVOKUmSpPZmOZPa2fSqVdQ2NNG36A2OPepBlvXcQ9+GHKO3nMm93/oxIZNJHVGSJEkJWM6kdrZlRzVnHvnP/KH3m6wCRm0dwAvvXMZrsafFTJIkqYBZzqR2tHDxnQwZdj9LiwMn7O7Km5svYWH9cQAMLC9JnE6SJEkpWc6kdrB27SJuW1DJc3E3HwuB4RvG8usdE/Z9vqQ4S+WE4QkTSpIkKTXLmXQY1dRs4b75k3lw50qKI0zpcxoXn3cP/2/Fdja22K2xcsJwJo0amDquJEmSErKcSYdBzOWoWvRDbl8zh83ZwOeK+zBl/I/o2+8EACaNKrWMSZIkaT+WM6mNrV4zn1sW3cALoY7hIcttJ13HySMvTR1LkiRJec5yJrWRnTvWMavqKh6ueZ3SCN/tP5YvjZtBtqhL6miSJEnqACxnUivlmhp54tmp3LXuabYFuKjbIK45fya9KoaljiZJkqQOxHImtcIrK+YwbfE0lmUaGUkXZo35HiOOn5Q6liRJkjogy5l0CKq3rubun13F3LqNVES4afCFXHD2jWSy/iclSZKkQ+MjSekjaGzYw6PPTGHmpoXUBri09FiunDiLHj37p44mSZKkDs5yJn1ILy59gJuX3M3qTI4xme5MPetmjhk6LnUsSZIkdRKWM+m/sGnTUmY8cy3zm7YxIMKdw77KuNO/TchkUkeTJElSJ2I5kw6ivm4XD1Vdzb1bX6YpwJVlJ3LZhJmUdK9IHU2SJEmdkOVMOoCFi+/i1uX3szYL5xaVUTl2BoMGjUkdS5IkSZ2Y5UxqYe3aRdy2oJLn4m6GALNHXMEZp1ydOpYkSZIKgOVMAmpqtnDf/Mk8uHMlxRGm9DmNi8+7h+KupamjSZIkqUBYzlTQYi5H1aJ/4PY1j7E5G/hccR+mjP8RffudkDqaJEmSCozlTAVr9Zr53LLoBl4IdQwPWW476TpOHnlp6liSJEkqUJYzFZydO9Yxq+oqHq55ndII3+0/li+Nm0G2qEvqaJIkSSpgljMVjFxTI088O5W71j7Ntgxc1G0Q15w/k14Vw1JHkyRJkixnKgyvrJjDtMXTWJZpZGTowqwx32PE8ZNSx5IkSZL2sZypU6uuXsPdVZOZW7eRigg3Db6QC86+kUzWpS9JkqT84iNUdRrzlqxnetUqNmyvZVBZlguHrQg19gAAF8lJREFUPMrcPS9SG+CS0mFMnjibHj37p44pSZIkHZDlTJ3CvCXrmfr4Mmobmjih+7PE3vN5qD5wcq4b3z/3Vo4ZOi51REmSJOkDWc7UKUyvWkXP+CqjBj3Isp576NsQGb5+DKvDly1mkiRJ6hAsZ+rw6ut2cUzRDP4w7E1WAaO2DuCFdy7j1diTQF3qeJIkSdKHYjlTh7Zw8V3cuvx+1vaBE3Z35c3Nl7Cw/rh9nx9QXpIwnSRJkvThWc7UIa1bu4hbF1TyXNzNEODaXl/hjjUnU9vQtO+ckuIslROGJ8soSZIkfRSWM3UoNTVbuG/+ZB7cuZLiCFP6nMbF591DcddS+gz+026NA8pLqJwwnEmjBqaOLEmSJH0oljN1CDGXo2rRP3D7msfYnA18rrgPU8b/iL79Tth3zqRRAy1jkiRJ6rAsZ8p7q9fM55ZFN/BCqGN4yHLbSddx8shLU8eSJEmS2pTlTHlr5451zKq6iodrXqc0wnePGsuXxs8gW9QldTRJkiSpzVnOlHdyTY088exU7lr7NNsycFG3QVxz/kx6VQxLHU2SJEk6bCxnyiuvrJjDtMXTWJZpZGTowqwx32PE8ZNSx5IkSZIOO8uZ8kJ19RrurprM3LqNVES4afCFXHD2jWSyLlFJkiQVBh/5KqlcUz0/efoqZm5aSG2AS0qHMXnibHr07J86miRJktSuLGdK5sWlD3DPG3fyahGMyXRn6lk3c8zQcaljSZIkSUlYztTuNm1ayh3P/A1PN1VzFJE7hn2V8adfT8hkUkeTJEmSkrGcqd3U1+3ioaqruXfryzQFuLLsRIb3+Crj/9vnU0eTJEmSkrOcqV0sXHw3ty6/j7VZOLeojMqxMxg0aAwLFixIHU2SJEnKC5Yztal5S9YzvWoVG7bXMqC8hKtP3cXzb93Cc3E3Q4DZI67gjFOuTh1TkiRJyjuWM7WZeUvWM/XxZdQ2NFESdjKs651MX7uJ4hiZ0mcMF593D8VdS1PHlCRJkvKS5UxtZnrVKmobGji1bC5b+i5mSVGGT+0qZfu7V/L1r389dTxJkiQpr1nO1Ga61v6a0UfPYWX3Jo6uyzB4/Wd5vuZsQupgkiRJUgdgOVOr7dyxjllVV1E99HX25CKf2vwJfl19Cbnm5TWgvCRxQkmSJCn/HXI5CyEMBh4CjgJywL0xxrvbKpjyX66pkSee/Q53rX2KbRmYQD9++eYlPF/XZ985JcVZKicMT5hSkiRJ6hha88xZI3BdjPG3IYSewMshhJ/HGFe0UTblsVdWzOHmxdP4faaRkaELPx7zd3zy+L/8s90aKycMZ9KoganjSpIkSXnvkMtZjHEjsLH577tCCCuBgYDlrBOrrl7D3VWTmVu3kYoINw2+kAvOvpFMdu9SmjRqoGVMkiRJOgRtcs9ZCGEIMApY3Bb/nvJPY8MeHn1mCjM3LaQ2wCWlw5g8cTY9evZPHU2SJEnqFEKMsXX/QAg9gOeAm2KMjx/g85cDlwP069fv04888kirvt/hsHv3bnr06JE6Rt7auOVnPLrj//JqEXy6oYhJfb5GedlJbfJvO/t0nH06zj4dZ5+Os0/L+afj7NPJ19mPHTv25Rjj6AN9rlXlLIRQDPwHUBVjvOO/On/06NHxpZdeOuTvd7gsWLCAc845J3WMvLNp01LueOZveLqpmv5NUPnxrzD+9OsJmUybfQ9nn46zT8fZp+Ps03H2aTn/dJx9Ovk6+xDCQctZa3ZrDMD9wMoPU8zUcdTX7eKhqqu5d+vLNAW4suxELpswk5LuFamjSZIkSZ1Wa+45OwO4BFgWQljafOw7McanWh9LqSxcfDe3Lr+PtVk4t6iMyrEzGDRoTOpYkiRJUqfXmt0aFwGhDbMooXXrnufWZ/83z8XdDAFmj7iCM065OnUsSZIkqWC0yW6N6lhavhfZ0PIGzhrwzzyZe5XiCFP6nMbF591DcdfS1DElSZKkgmI5KzDzlqxn6uPLqG1o4NSyuWzps5ifxgxjYy9u+ItZ9O13QuqIkiRJUkGynBWY6VWrOCrzMr2PnsPK7k0cXZdh8PrP8nKXiRYzSZIkKSHLWQHZuWMdw7r9Pb8/ait1ucinNn+CX1dfQo4iQk1t6niSJElSQbOcFYBcUyNPPPsd7lr3FNvKYeTOMpZtvoznm47ad86A8pKECSVJkiRZzjq5V1bM4ebF0/h9ppGRdOGy/t/kllcHUtvUtO+ckuIslROGJ0wpSZIkyXLWSVVXr+HuqsnMrdtIRYSbBl/IBWffSCZbRFnfP+3WOKC8hMoJw5k0amDqyJIkSVJBs5x1Mo0Ne3j0meuYuek5agNcUjqMyRNn06Nn/33nTBo10DImSZIk5RnLWSfy4tJ/5uYld7E6k2NMpoSpZ97MMceMTx1LkiRJ0odgOesENm3+HXf8/Fqebqqmf4zcMeyrjD/9ekImkzqaJEmSpA/JctbBzFvyp/vFBpfBxI/9G/Pql9EU4MqyE7hswo8p6V6ROqYkSZKkj8hy1oHMW7KeqY8vo7ahiZN6Pk3tkQt4uDHwmVjK342/i8GDP5M6oiRJkqRDZDnrQKZXraJXWMnIwf/G8h71DKiHY9edzfKiL1rMJEmSpA7OctZB1NZUM6zLLSzvt4HXiIx6Zyi/2fJX1FNCoDZ1PEmSJEmtZDnLczGX4z9/dSu3/fEnbOwd+NSu7qze9DUWNg7Zd86A8pJ0ASVJkiS1CctZHnvt9V9w88Lr+Q21HBeyfKv315i+ZgS1jU37zikpzlI5YXjClJIkSZLaguUsD+3etZF/rLqKf9u9mpII1x91Fv9z/B0UFXejYuCfdmscUF5C5YThvqG0JEmS1AlYzhJruTV+/7KufPnY/2DurmfYkg18sWt/rj1vJr2P/Pi+8yeNGmgZkyRJkjohy1lCLbfGP67bS/Ts9VP+qSZyfCzintHf5cRPfil1REmSJEntxHKW0PSqVRQ3beKU/vfxu7Lt1OYiJ2w6kbW5yzjxk+enjidJkiSpHVnOEmlqrGdIZjZvHruc32UCn9pRwdLNl/FWrg+BhtTxJEmSJLUzy1kCv1/xGNMWT2N5vyaOqy2i66aL+OWek/d93q3xJUmSpMJjOTuMWm72MaC8hGvP7MYra3/AvPqN9I2Ry0vP58drxlPbEPd9jVvjS5IkSYXJcnaYtNzsI9DAgOyD3LVqBXsyga+XHssVE2dR2rM/g4e5Nb4kSZIky9lhM71qFbUNTQwrWUy3/vP4Q9fIJ2uKiLsuZsr/qdx3nlvjS5IkSQLL2WGzc9c6Th9wH8vKtnNkY46RGz7N8zsuArKpo0mSJEnKQ5azVnr/fWXXnXcsYedsyoY9wYpM4JRtFbz49jd4PdcbgIFu9iFJkiTpACxnrdDyvjKAoj2/5l9f+hZruuU4IRaze91f8ot3R+87380+JEmSJB1MJnWAjuy9+8q6hnc5s98d7BjyMO8UN/KZ6hP5yV+9yF9f8A0GlpcQ2PuM2c1fPNH7yyRJkiQdkM+ctcKG7bV8ovtz5Po/xdIugZN29uSVTd/g501HkckWu9mHJEmSpA/NcvYRvXeP2fadmziz/yyWlG+nX0Nk2Lpx/HL3+YD3lUmSJEn66CxnH8F795h9vOtTdDvmP1laFBi1rYKXNl9JTSwDvK9MkiRJ0qGxnP0XWu7GWJbdwuh+s1latpuB9YFBaz/HwpqzyIZAIPom0pIkSZIOmeXsIOYtWc8PnlzO9toGIHLqEXN5u99vWJYNjKrux+K3r6AulgKQi5HXb/lc2sCSJEmSOjTL2fvsX8qgf/YNhg14gN/1qOfougzF677Iwj2n7vc1A7zHTJIkSVIrWc5aaPm+Zd3DLj7d+ye8WvEaK0PgxHeGsnjLZTTRZb+v8R4zSZIkSW3BctbC9KpV1DbWcNxRD1F7xGqWZjMMr+nKWxsv5lf1fypg2RDIRe8xkyRJktR2LGctbNheS2nvBWzq9Sqn78rwxrbP8tK7Z+93Tklx1jeTliRJktTmLGct9C/vSkOvX3NSbS3V6yezKh6/3+d7dS/m+5//pMVMkiRJUpuznDWbt2Q97zZuJVe8h1O3B26Kf3oZo6VMkiRJ0uFmOQN+taGBf/3PZTRkqykBNjYcDQRLmSRJkqR2k0kdIB/89I8N1DY00btoEwDrGo8GoHuXIouZJEmSpHZhOQO27okAHJF9G4Adjf2AvRuESJIkSVJ7sJwBvbsFALoXVQOwtfEowDeXliRJktR+LGfAf/94MSXFWboU7aQ4Rt5uOso3l5YkSZLUrtwQBDh9QDEjPjGCOb96l15NOfaU9aZy4vHebyZJkiSp3VjOmk0aNZCfvdJEQ0MRj00dlzqOJEmSpALjyxpb2NZYQ0W2W+oYkiRJkgqQ5ayF6thIRVFp6hiSJEmSCpDlrIXqEOnV5YjUMSRJkiQVIMsZ8KsNDZx7yzz2ZAJ/3BCZt2R96kiSJEmSCkzBl7N5S9bz4Cv17Kl5E4Dde0qZ+vgyC5okSZKkdlXw5Wx61Srqc9Az+zYAdY3l1DY0Mb1qVeJkkiRJkgpJwZezDdtrASgtqgbg3caK/Y5LkiRJUnso+HI2oLwEgK5F2wHY1dR3v+OSJEmS1B4KvpxVThhOlwwUZXcCsK2hLyXFWSonDE+cTJIkSVIhKUodILVJowayYuUKXt1aS7dcpEtZPyonDGfSqIGpo0mSJEkqIAVfzgBOH1DMrgi96wLzrz83dRxJkiRJBajgX9b4nurGd+kV7KqSJEmS0rCcsfdNqDfW1ZCrDZxxyy98jzNJkiRJ7a7gy9l7b0K9K5ujuLEr67fX+ibUkiRJktpdwZezvW9CnWN7NpBp2rt9vm9CLUmSJKm9FXw527C9lp6Z7TSGQK6px37HJUmSJKm9FHw5G1BeQlnRZgAaG4/Y77gkSZIktZeCL2eVE4bTq3gLAHsaywF8E2pJkiRJ7c6944HS4q0AvNvYm17di/n+5z/pm1BLkiRJalcF/czZvCXrmfr4MjKZnQDsaOzLnoZc4lSSJEmSClFBl7PpVauobWgim91bzrY39XOnRkmSJElJFHQ5e29HxpB9l9JcjvpYst9xSZIkSWovBV3O3tuRsaloD0c0/vlxSZIkSWovBV3OKicMp6Q4S122ntJcFnCnRkmSJElpFPRuje/tyHjfS430qO/KwPISKicMd6dGSZIkSe2uoJ85g70FbXc2cnxFb56//lyLmSRJkqQkCr6c5Zoa2Z4JVHQtSx1FkiRJUgEr+HK2c+c6mkKgoqR36iiSJEmSCljBl7Pq7a8BUNG9X+IkkiRJkgqZ5WzHOgB69eifOIkkSZKkQmY52/UWABVHDE6cRJIkSVIhs5y9uxmA3uXHJE4iSZIkqZAVfDnbtmcrAGXlRydOIkmSJKmQFXw521q3nbKmHMXF3VNHkSRJklTAWlXOQggTQwirQghrQgjXt1Wo9lRdv4vyGFLHkCRJklTgDrmchRCywEzgs8AI4CshhBFtFay9bGuqpTxmU8eQJEmSVOBa88zZqcCaGONrMcZ64BHgC20Tq/1U5+opo0vqGJIkSZIKXGvK2UBgXYuP32o+1qEcmelCv2yv1DEkSZIkFbgQYzy0LwzhS8CEGOM3mj++BDg1xnjN+867HLgcoF+/fp9+5JFHWpf4MNi9ezc9evRIHaMgOft0nH06zj4dZ5+Os0/L+afj7NPJ19mPHTv25Rjj6AN9rqgV/+5bQMt3bh4EbHj/STHGe4F7AUaPHh3POeecVnzLw2PBggXkY65C4OzTcfbpOPt0nH06zj4t55+Os0+nI86+NS9rfBE4LoQwNITQBfgy8GTbxJIkSZKkwnLIz5zFGBtDCFcDVUAWeCDGuLzNkkmSJElSAWnNyxqJMT4FPNVGWSRJkiSpYLXqTaglSZIkSW3DciZJkiRJecByJkmSJEl5wHImSZIkSXnAciZJkiRJecByJkmSJEl5wHImSZIkSXnAciZJkiRJecByJkmSJEl5wHImSZIkSXnAciZJkiRJecByJkmSJEl5wHImSZIkSXnAciZJkiRJecByJkmSJEl5wHImSZIkSXnAciZJkiRJecByJkmSJEl5wHImSZIkSXkgxBjb75uF8A7wZrt9ww/vSGBL6hAFytmn4+zTcfbpOPt0nH1azj8dZ59Ovs7+6BhjnwN9ol3LWb4KIbwUYxydOkchcvbpOPt0nH06zj4dZ5+W80/H2afTEWfvyxolSZIkKQ9YziRJkiQpD1jO9ro3dYAC5uzTcfbpOPt0nH06zj4t55+Os0+nw83ee84kSZIkKQ/4zJkkSZIk5YGCLmchhIkhhFUhhDUhhOtT5+nMQgiDQwjPhhBWhhCWhxD+pvn4D0II60MIS5v/9xeps3ZWIYQ3QgjLmuf8UvOxihDCz0MIq5v/7JU6Z2cTQhjeYn0vDSHsDCH8rWv/8AghPBBCeDuE8EqLYwdc52Gve5p/Bvw+hHByuuQd30FmPz2E8Ifm+c4NIZQ3Hx8SQqhtsf5np0ve8R1k9ge9xoQQpjav+1UhhAlpUncOB5n9oy3m/kYIYWnzcdd9G/qAx5Yd+ppfsC9rDCFkgT8C5wFvAS8CX4kxrkgarJMKIfQH+scYfxtC6Am8DEwC/gewO8Z4e9KABSCE8AYwOsa4pcWx24DqGOMtzb+g6BVj/HaqjJ1d83VnPXAa8HVc+20uhHAWsBt4KMZ4QvOxA67z5ger1wB/wd7/T+6OMZ6WKntHd5DZnw/8IsbYGEK4FaB59kOA/3jvPLXOQWb/Aw5wjQkhjAAeBk4FBgDPAB+PMTa1a+hO4kCzf9/nZwA7Yow3uu7b1gc8tvwrOvA1v5CfOTsVWBNjfC3GWA88AnwhcaZOK8a4Mcb42+a/7wJWAgPTphJ71/y/NP/9X9h7UdPhMw54Ncb4ZuognVWMcSFQ/b7DB1vnX2DvA6oYY/wNUN78w16H4ECzjzH+LMbY2Pzhb4BB7R6sABxk3R/MF4BHYox1McbXgTXsfUykQ/BBsw8hBPb+Evrhdg1VID7gsWWHvuYXcjkbCKxr8fFbWBbaRfNvjkYBi5sPXd389PIDvqzusIrAz0IIL4cQLm8+1i/GuBH2XuSAvsnSFYYvs/8Padd++zjYOvfnQPu6DHi6xcdDQwhLQgjPhRDOTBWqkzvQNcZ1337OBDbHGFe3OOa6Pwze99iyQ1/zC7mchQMcK8zXeLajEEIP4KfA38YYdwKzgGHAScBGYEbCeJ3dGTHGk4HPAt9sfimG2kkIoQtwITCn+ZBrPz1/DrSTEMJ3gUbgJ82HNgIfizGOAqYA/x5COCJVvk7qYNcY1337+Qr7/0LOdX8YHOCx5UFPPcCxvFv7hVzO3gIGt/h4ELAhUZaCEEIoZu9/PD+JMT4OEGPcHGNsijHmgH/Cl1YcNjHGDc1/vg3MZe+sN7/3lH7zn2+nS9jpfRb4bYxxM7j229nB1rk/B9pBCOFrwAXA/4rNN7o3v6Rua/PfXwZeBT6eLmXn8wHXGNd9OwghFAFfBB5975jrvu0d6LElHfyaX8jl7EXguBDC0ObfaH8ZeDJxpk6r+XXX9wMrY4x3tDje8rW+fwm88v6vVeuFEEqbb5YlhFAKnM/eWT8JfK35tK8BT6RJWBD2+w2qa79dHWydPwlc2ryD1xj23rS/MUXAziqEMBH4NnBhjLGmxfE+zRvkEEI4BjgOeC1Nys7pA64xTwJfDiF0DSEMZe/sX2jvfAVgPPCHGONb7x1w3betgz22pINf84tSB0ileeeoq4EqIAs8EGNcnjhWZ3YGcAmw7L0tZYHvAF8JIZzE3qeV3wCuSBOv0+sHzN17HaMI+PcY4/wQwovAYyGEvwbWAl9KmLHTCiF0Z+/OsC3X922u/bYXQngYOAc4MoTwFvB94BYOvM6fYu+uXWuAGvbuoKlDdJDZTwW6Aj9vvv78JsZ4JXAWcGMIoRFoAq6MMX7YDS30PgeZ/TkHusbEGJeHEB4DVrD3pabfdKfGQ3eg2ccY7+fP7zEG131bO9hjyw59zS/YrfQlSZIkKZ8U8ssaJUmSJClvWM4kSZIkKQ9YziRJkiQpD1jOJEmSJCkPWM4kSZIkKQ9YziRJkiQpD1jOJEmSJCkPWM4kSZIkKQ/8fx4xkLsPagwEAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "from dolark.dolo_improvements import Linear, \\\n", + " EmptyGrid, WarpGrid, CartesianGrid, ICartesianGrid, \\\n", + " UnstructuredGrid, DecisionRule, Cubic\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n", + "\n", + "# placeholder\n", + "exo_grid = UnstructuredGrid(np.array([[0.2, 0.5, 0.7]]))\n", + "\n", + "# construction of the endogenous grid\n", + "N = 20\n", + "a = np.exp( np.linspace(-1,5,N) )\n", + "endo_grid = ICartesianGrid([a])\n", + "# f = lambda x: (x[:,0]**2)[:,None]\n", + "f = lambda x: np.minimum(x[:,0],1+0.05*(1+(x[:,0]-1)*(1-np.exp(1-(x[:,0])))))[:,None]\n", + "nodes = endo_grid.nodes\n", + "vals = f(nodes).reshape((1,N,1))\n", + "\n", + "# Creation of the decision rule\n", + "dr = DecisionRule(exo_grid, endo_grid, 'linear')\n", + "dr.set_values(vals)\n", + "\n", + "# Evaluation on a finer grid\n", + "fg = CartesianGrid([0.01], [200], [1000])\n", + "nn = fg.nodes\n", + "tvals = f(nn)\n", + "xx = dr.eval_is(0, nn)\n", + "\n", + "plt.figure(figsize=(15,10))\n", + "plt.plot(nodes.ravel(), vals.ravel(), 'o', label='data')\n", + "plt.plot(nn.ravel(), xx.ravel(), label='linear')\n", + "plt.plot(nn.ravel(), tvals.ravel(), label='true')\n", + "plt.legend(loc='upper left')\n", + "plt.grid()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Warped interpolation\n", + "\n", + "We need to approximate a function $x=\\varphi(s)$ where $s=\\mathcal{W}^{1}(\\tilde{s})$. The idea is to discretize the values of $s$ so that the values of $\\tilde{s}$ are discretized uniformly. Then we can say that the grid for $s$ is a warped grid.\n", + "\n", + "First syntax we propose defines a grid for $s$ and a tranformation function (or a list of transformation functions if there are many dimensions)." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "from dolark.dolo_improvements import multimethod, get_coefficients, eval_is, eval_ms, Linear, EmptyGrid, WarpGrid, CartesianGrid, UnstructuredGrid, DecisionRule, Cubic\n", + "import numpy as np\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "base = CartesianGrid([0.0,0], [1.0,1.0], [10,10])\n", + "wg = WarpGrid(\n", + " base,\n", + " ['exp',''] # list of transformation functions (so far only id, exp, exp(exp), exp(exp(exp)) are supported)\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The equivalent yaml would be: (not implemented)\n", + "\n", + "```yaml\n", + "!WarpGrid\n", + " base: !CartesianGrid\n", + " a: [0.0, 0.0]\n", + " b: [1.0, 1.0]\n", + " n: [10, 10]\n", + " warp:\n", + " [exp, '']\n", + " \n", + "```\n", + "\n", + "This is a one to one mapping with the python object, but many other syntaxes could be feasible. For instance:\n", + "\n", + "```yaml\n", + "!CartesianGrid\n", + " a: [0.0, 0.0]\n", + " b: [1.0, 1.0]\n", + " n: [10, 10]\n", + " warp: [exp, '']\n", + "``` \n", + "\n", + "However, it is not clear to me that such an option if interesting enough to be added to the default grid constructor.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD4CAYAAAATpHZ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAYOUlEQVR4nO3df8ydZ10G8Oui22TCZBt9Udg6OrUQCwEHJ+PHjI4AWTdl1WRKqxAgSBNj/RHJkqFmmCGJshjBMMSGEByRLRNhNqRYiM6QAMOesrHRzmotG3tXdC/bGCLDUfz6xzndTp8993uetue+7+/33vVJmr3nPE/Pc90nz747e855z0Uzg4iItOEptQOIiMjiaKiLiDREQ11EpCEa6iIiDdFQFxFpyCm1Drx27Vpbv359rcOLiIS0d+/eb5rZUmp7taG+fv16jMfjWocXEQmJ5D2rbdflFxGRhmioi4g0RENdRKQhGuoiIg3RUBcRacjcT7+Q/DCAXwBwv5m9sGc7AbwPwGUAvgvgzWb25UUHBYCbb7sP1+4+gMPfegTPOfN0XHnJ8/GLF5zj7jFLPHaN49Q+pscMnvP0iZCxK2LmlBJrGfKRxo8AeD+A6xPbLwWwYfrnZQD+cvrPhbr5tvvwjk/ciUe+/wMAwH3fegTv+MSdAHDCT0qOxyzx2DWOU/uYHjN4ztMnQsauiJlTSq1l7uUXM/scgAdX2WUzgOtt4lYAZ5J89qICHnXt7gOPPRlHPfL9H+Da3QdcPWaJx65xnNrH9JjBc54+ETJ2RcycUmoti7imfg6Ae2duL0/vewKS20iOSY5XVlaO6yCHv/XIcd1f6zFLPHaN49Q+pscMQ45bK0+fCBm7ImZOKbWWRQx19tzX27xhZjvMbGRmo6Wl5G+59nrOmacf1/21HrPEY9c4Tu1jesww5Li18vSJkLErYuaUUmtZxFBfBrBu5va5AA4v4HGPceUlz8fpp6455r7TT12DKy95vqvHLPHYNY5T+5geM3jO0ydCxq6ImVNKrWUR3/2yE8B2kjdi8gbpw2b2jQU87jGOvpGwyHeOczxmiceucZzax/SYwXOePhEydkXMnFJqLZzXUUryBgAXA1gL4L8AvBPAqQBgZh+cfqTx/QA2YfKRxreY2dxv6hqNRqYv9BIROT4k95rZKLV97it1M9s6Z7sB+M0TyCYiIgum3ygVEWmIhrqISEM01EVEGqKhLiLSEA11EZGGaKiLiDREQ11EpCEa6iIiDdFQFxFpiIa6iEhDNNRFRBqioS4i0hANdRGRhizi+9SLydHEnbPdu1QLeo22dQ8N7x4yeM7TJ0LGroiZU0qsJcxQz9HEnbPdu1RzeI22dQ8N7x4yeM7TJ0LGroiZU0qtJczllxxN3DnbvUs1h9doW/fQ8O4hg+c8fSJk7IqYOaXUWsIM9RxN3DnbvUs1h9doW/fQ8O4hw5Djemq9j5CxK2LmlFJrCTPUczRx52z3LtUcXqNt3UPDu4cMQ47rqfU+QsauiJlTSq0lzFDP0cSds927VHN4jbZ1Dw3vHjJ4ztMnQsauiJlTSq0lzBulOZq4c7Z7l2oOr9G27qHh3UMGz3n6RMjYFTFzSqm1cNIbXd5oNLLxeFzl2CIiUZHca2aj1PYwl19ERGQ+DXURkYZoqIuINERDXUSkIRrqIiIN0VAXEWmIhrqISEM01EVEGqKhLiLSEA11EZGGDBrqJDeRPEDyIMmrerafR/IWkreRvIPkZYuPKiIi88wd6iTXALgOwKUANgLYSnJjZ7c/BHCTmV0AYAuADyw6qIiIzDfklfqFAA6a2SEzexTAjQA2d/YxAD8y/fkZAA4vLqKIiAw1ZKifA+DemdvL0/tm/RGAN5BcBrALwG/1PRDJbSTHJMcrKysnEFdERFYz5PvU2XNf9/t6twL4iJn9GclXAPgoyRea2f8d85fMdgDYAUy+evd4w+Zo4s7Z7l2qBb1G27qHhncPGTzn6RMhY1fEzCkl1jJkqC8DWDdz+1w88fLKWwFsAgAz+yLJpwJYC+D+RYQE8jRx52z3LtUcXqNt3UPDu4cMnvP0iZCxK2LmlFJrGXL5ZQ+ADSTPJ3kaJm+E7uzs83UArwYAkj8F4KkAFnp9JUcTd85271LN4TXa1j00vHvI4DlPnwgZuyJmTim1lrlD3cyOANgOYDeAuzD5lMs+kteQvHy629sBvI3kVwDcAODNtuBKpRxN3DnbvUs1h9doW/fQ8O4hw5Djemq9j5CxK2LmlFJrGfQ5dTPbZWbPM7OfMLN3T++72sx2Tn/eb2YXmdmLzeynzewzC02JPE3cOdu9SzWH12hb99Dw7iHDkON6ar2PkLErYuaUUmsJ8xulOZq4c7Z7l2oOr9G27qHh3UMGz3n6RMjYFTFzSqm1DHmj1IUcTdw5271LNYfXaFv30PDuIYPnPH0iZOyKmDml1Fq44Evfg41GIxuPx1WOLSISFcm9ZjZKbQ9z+UVERObTUBcRaYiGuohIQzTURUQaoqEuItIQDXURkYZoqIuINERDXUSkIRrqIiIN0VAXEWmIhrqISEM01EVEGqKhLiLSEA11EZGGhPk+dSBPE3fOdu9SLeg12tY9NLx7yOA5T58IGbsiZk4psZYwQz1HE3fOdu9SzeE12tY9NLx7yOA5T58IGbsiZk4ptZYwl19yNHHnbPcu1Rxeo23dQ8O7hwye8/SJkLErYuaUUmsJM9RzNHHnbPcu1Rxeo23dQ8O7hwxDjuup9T5Cxq6ImVNKrSXMUM/RxJ2z3btUc3iNtnUPDe8eMgw5rqfW+wgZuyJmTim1ljBDPUcTd85271LN4TXa1j00vHvI4DlPnwgZuyJmTim1ljBvlOZo4s7Z7l2qObxG27qHhncPGTzn6RMhY1fEzCml1kIzW+gDDjUajWw8Hlc5tohIVCT3mtkotT3M5RcREZlPQ11EpCEa6iIiDdFQFxFpiIa6iEhDBg11kptIHiB5kORViX1+heR+kvtIfmyxMUVEZIi5n1MnuQbAdQBeC2AZwB6SO81s/8w+GwC8A8BFZvYQyWflCiwiImlDXqlfCOCgmR0ys0cB3Ahgc2eftwG4zsweAgAzu3+xMUVEZIghQ/0cAPfO3F6e3jfreQCeR/LzJG8luanvgUhuIzkmOV5ZWTmxxCIikjRkqLPnvu6voZ4CYAOAiwFsBfAhkmc+4S+Z7TCzkZmNlpaWjjeriIjMMWSoLwNYN3P7XACHe/b5ezP7vpl9DcABTIa8iIgUNGSo7wGwgeT5JE8DsAXAzs4+NwN4FQCQXIvJ5ZhDiwwqIiLzzR3qZnYEwHYAuwHcBeAmM9tH8hqSl0932w3gAZL7AdwC4EozeyBXaBER6advaRQRCWTetzSG+T51IE8Td85271It6DXa1j00vHvI4DlPnwgZuyJmTimxljBDPUcTd85271LN4TXa1j00vHvI4DlPnwgZuyJmTim1ljDf/ZKjiTtnu3ep5vAabeseGt49ZPCcp0+EjF0RM6eUWkuYoZ6jiTtnu3ep5vAabeseGt49ZBhyXE+t9xEydkXMnFJqLWGGeo4m7pzt3qWaw2u0rXtoePeQYchxPbXeR8jYFTFzSqm1hBnqOZq4c7Z7l2oOr9G27qHh3UMGz3n6RMjYFTFzSqm1hHmjNEcTd85271LN4TXa1j00vHvI4DlPnwgZuyJmTim1Fn1OXUQkkHmfUw9z+UVERObTUBcRaYiGuohIQzTURUQaoqEuItIQDXURkYZoqIuINERDXUSkIRrqIiIN0VAXEWmIhrqISEM01EVEGqKhLiLSEA11EZGGhPk+dSBPE3fOdu9SLeg12tY9NLx7yOA5T58IGbsiZk4psZYwQz1HE3fOdu9SzeE12tY9NLx7yOA5T58IGbsiZk4ptZYwl19yNHHnbPcu1Rxeo23dQ8O7hwye8/SJkLErYuaUUmsJM9RzNHHnbPcu1Rxeo23dQ8O7hwxDjuup9T5Cxq6ImVNKrSXMUM/RxJ2z3btUc3iNtnUPDe8eMgw5rqfW+wgZuyJmTim1ljBDPUcTd85271LN4TXa1j00vHvI4DlPnwgZuyJmTim1ljBvlOZo4s7Z7l2qObxG27qHhncPGTzn6RMhY1fEzCml1kIzW+gDDjUajWw8Hlc5tohIVCT3mtkotT3M5RcREZlv0FAnuYnkAZIHSV61yn5XkDSSyf+KiIhIPnOHOsk1AK4DcCmAjQC2ktzYs98ZAH4bwJcWHVJERIYZ8kr9QgAHzeyQmT0K4EYAm3v2exeA9wD43gLziYjIcRgy1M8BcO/M7eXpfY8heQGAdWb2qdUeiOQ2kmOS45WVleMOKyIiqxsy1Nlz32MfmSH5FAB/DuDt8x7IzHaY2cjMRktLS8NTiojIIEOG+jKAdTO3zwVweOb2GQBeCOCfSd4N4OUAdurNUhGR8oYM9T0ANpA8n+RpALYA2Hl0o5k9bGZrzWy9ma0HcCuAy81MH0IXESls7lA3syMAtgPYDeAuADeZ2T6S15C8PHdAEREZbtDXBJjZLgC7Ovddndj34pOPJSIiJ0K/USoi0hANdRGRhmioi4g0RENdRKQhYb5PHcjTxJ2z3btUC3qNtnUPDe8eMnjO0ydCxq6ImVNKrCXMUM/RxJ2z3btUc3iNtnUPDe8eMnjO0ydCxq6ImVNKrSXM5ZccTdw5271LNYfXaFv30PDuIYPnPH0iZOyKmDml1FrCDPUcTdw5271LNYfXaFv30PDuIcOQ43pqvY+QsSti5pRSawkz1HM0ceds9y7VHF6jbd1Dw7uHDEOO66n1PkLGroiZU0qtJcxQz9HEnbPdu1RzeI22dQ8N7x4yeM7TJ0LGroiZU0qtJcwbpTmauHO2e5dqDq/Rtu6h4d1DBs95+kTI2BUxc0qptdDM5u+VwWg0svFYX+QoInI8SO41s+RXm4e5/CIiIvNpqIuINERDXUSkIRrqIiIN0VAXEWmIhrqISEM01EVEGqKhLiLSEA11EZGGaKiLiDREQ11EpCEa6iIiDdFQFxFpiIa6iEhDwnyfOpCniTtnu3epFvQabeseGt49ZPCcp0+EjF0RM6eUWEuYoZ6jiTtnu3ep5vAabeseGt49ZPCcp0+EjF0RM6eUWkuYyy85mrhztnuXag6v0bbuoeHdQwbPefpEyNgVMXNKqbWEGeo5mrhztnuXag6v0bbuoeHdQ4Yhx/XUeh8hY1fEzCml1jJoqJPcRPIAyYMkr+rZ/nsk95O8g+Q/knzuQlMiTxN3znbvUs3hNdrWPTS8e8gw5LieWu8jZOyKmDml1FrmDnWSawBcB+BSABsBbCW5sbPbbQBGZvYiAB8H8J6FpkSeJu6c7d6lmsNrtK17aHj3kMFznj4RMnZFzJxSai1D3ii9EMBBMzsEACRvBLAZwP6jO5jZLTP73wrgDYsMCeRp4s7Z7l2qObxG27qHhncPGTzn6RMhY1fEzCml1kIzW30H8goAm8zs16e33wjgZWa2PbH/+wH8p5n9cc+2bQC2AcB555330nvuueck44uIPLmQ3Gtmo9T2IdfU2XNf738JSL4BwAjAtX3bzWyHmY3MbLS0tDTg0CIicjyGXH5ZBrBu5va5AA53dyL5GgB/AODnzOx/FxNPRESOx5BX6nsAbCB5PsnTAGwBsHN2B5IXAPgrAJeb2f2LjykiIkPMHepmdgTAdgC7AdwF4CYz20fyGpKXT3e7FsDTAfwtydtJ7kw8nIiIZDToawLMbBeAXZ37rp75+TULziUiIicgzG+UiojIfBrqIiIN0VAXEWmIhrqISEM01EVEGqKhLiLSEA11EZGGaKiLiDREQ11EpCFhiqeBPE3cOdu9S7Wg12hb99Dw7iGD5zx9ImTsipg5pcRawgz1HE3cOdu9SzWH12hb99Dw7iGD5zx9ImTsipg5pdRawlx+ydHEnbPdu1RzeI22dQ8N7x4yeM7TJ0LGroiZU0qtJcxQz9HEnbPdu1RzeI22dQ8N7x4yDDmup9b7CBm7ImZOKbWWMEM9RxN3znbvUs3hNdrWPTS8e8gw5LieWu8jZOyKmDml1FrCDPUcTdw5271LNYfXaFv30PDuIYPnPH0iZOyKmDml1FrCvFGao4k7Z7t3qebwGm3rHhrePWTwnKdPhIxdETOnlFoLzXo7pLMbjUY2Ho+rHFtEJCqSe81slNoe5vKLiIjMp6EuItIQDXURkYZoqIuINERDXUSkIRrqIiIN0VAXEWmIhrqISEM01EVEGqKhLiLSEA11EZGGaKiLiDREQ11EpCEa6iIiDRn0feokNwF4H4A1AD5kZn/S2f5DAK4H8FIADwB4vZndvdioeZq4c7Z7l2pBr9G27qHh3UMGz3n6RMjYFTFzSom1zB3qJNcAuA7AawEsA9hDcqeZ7Z/Z7a0AHjKznyS5BcCfAnj9IoPmaOLO2e5dqjm8Rtu6h4Z3Dxk85+kTIWNXxMwppdYy5PLLhQAOmtkhM3sUwI0ANnf22Qzgr6c/fxzAq0lyYSmRp4k7Z7t3qebwGm3rHhrePWTwnKdPhIxdETOnlFrLkKF+DoB7Z24vT+/r3cfMjgB4GMAzuw9EchvJMcnxysrKcQXN0cSds927VHN4jbZ1Dw3vHjIMOa6n1vsIGbsiZk4ptZYhQ73vFXe3A2/IPjCzHWY2MrPR0tLSkHyPydHEnbPdu1RzeI22dQ8N7x4yDDmup9b7CBm7ImZOKbWWIUN9GcC6mdvnAjic2ofkKQCeAeDBRQQ8KkcTd85271LN4TXa1j00vHvI4DlPnwgZuyJmTim1liGfftkDYAPJ8wHcB2ALgF/t7LMTwJsAfBHAFQD+yRbcaJ2jiTtnu3ep5vAabeseGt49ZPCcp0+EjF0RM6eUWguHzF6SlwF4LyYfafywmb2b5DUAxma2k+RTAXwUwAWYvELfYmaHVnvM0Whk4/H4pBcgIvJkQnKvmY1S2wd9Tt3MdgHY1bnv6pmfvwfgl080pIiILIZ+o1REpCEa6iIiDdFQFxFpiIa6iEhDBn36JcuByRUA95zgX18L4JsLjFNCtMzKm1e0vEC8zK3mfa6ZJX97s9pQPxkkx6t9pMejaJmVN69oeYF4mZ+seXX5RUSkIRrqIiINiTrUd9QOcAKiZVbevKLlBeJlflLmDXlNXURE+kV9pS4iIj001EVEGuJqqJP8MMn7SX41sZ0k/4LkQZJ3kHzJzLY3kfz36Z83Ocr8a9Osd5D8AskXz2y7m+SdJG8nWeQrKwfkvZjkw9NMt5O8embbJpIHps//VU7yXjmT9askf0Dy7Om2Gs/vOpK3kLyL5D6Sv9Ozj5vzeGBeb+fwkMxuzuOBeRd3HpuZmz8AfhbASwB8NbH9MgCfxqRp6eUAvjS9/2wAh6b/PGv681lOMr/yaBYAlx7NPL19N4C1zp7jiwF8quf+NQD+A8CPAzgNwFcAbKydt7Pv6zD5Lv+az++zAbxk+vMZAP6t+zx5Oo8H5vV2Dg/J7OY8HpK3s/9JnceuXqmb2eewemPSZgDX28StAM4k+WwAlwD4rJk9aGYPAfgsgE35E8/PbGZfmGYCgFsxaY6qZsBznDKkgHzhjjPvVgA3ZIwzl5l9w8y+PP35vwHchSd2+ro5j4fkdXgOD3mOU4qfxyeQ96TOY1dDfYBUCfaQcmwP3orJK7SjDMBnSO4lua1Spj6vIPkVkp8m+YLpfa6fY5I/jMkA/LuZu6s+vyTXY1Ic86XOJpfn8Sp5Z7k6h+dkdncez3uOF3EeDyrJcCRVcD2o+Lomkq/C5F+In5m5+yIzO0zyWQA+S/Jfp69Ma/oyJt8t8R1OGq9uBrAB/p/j1wH4vJnNvqqv9vySfDom/2L+rpl9u7u5569UPY/n5D26j6tzeE5md+fxkOcYCziPo71ST5VgDynHrobkiwB8CMBmM3vg6P1mdnj6z/sBfBKT/zWsysy+bWbfmf68C8CpJNfC+XOMSXfuMf/LWuv5JXkqJv/y/o2ZfaJnF1fn8YC87s7heZm9ncdDnuOpkz+Pc75BcCJ/AKxH+k28n8exbzD9y/T+swF8DZM3l86a/ny2k8znATgI4JWd+58G4IyZn78AYJODvD+Gx38p7UIAX58+36dg8sbd+Xj8DaYX1M473f4MTK67P6328zt9rq4H8N5V9nFzHg/M6+ocHpjZzXk8JO8iz2NXl19I3oDJu9ZrSS4DeCeAUwHAzD6ISU/qZZicYN8F8JbptgdJvgvAnulDXWPH/u9LzcxXA3gmgA+QBIAjNvkmth8F8MnpfacA+JiZ/YODvFcA+A2SRwA8gkmJuAE4QnI7gN14vIB8n4O8APBLAD5jZv8z81erPL8ALgLwRgB3krx9et/vYzIYPZ7HQ/K6OocHZvZ0Hg/JCyzoPNbXBIiINCTaNXUREVmFhrqISEM01EVEGqKhLiLSEA11EZGGaKiLiDREQ11EpCH/Dxz+ZwUkcQmtAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(wg.nodes[:,0], wg.nodes[:,1], 'o')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Decision rules are adapted to deal with warped grids\n", + "\n", + "```python\n", + "@multimethod\n", + "def get_coefficients(itp, exo_grid, endo_grid: WarpGrid, interp_type, x):\n", + " return get_coefficients(itp, exo_grid, endo_grid.base, interp_type, x)\n", + "\n", + "@multimethod\n", + "def eval_is(itp, exo_grid, endo_grid: WarpGrid, interp_type, i, s):\n", + " base = endo_grid.base\n", + " ss = s.copy()\n", + " for k,f in enumerate(endo_grid.warp_ifunctions):\n", + " ss[:,k] = f(ss[:,k])\n", + " return eval_is(itp, exo_grid, base, interp_type, i, ss)\n", + "\n", + "@multimethod\n", + "def eval_ms(itp, exo_grid, endo_grid: WarpGrid, interp_type, m, s):\n", + " base = endo_grid.base\n", + " ss = s.copy()\n", + " for i,f in enumerate(endo_grid.warp_ifunctions):\n", + " ss[:,i] = f(ss[:,i])\n", + " return eval_ms(itp, exo_grid, base, interp_type, m, ss)\n", + "```\n", + "\n", + "- pros:\n", + " - this is compatible with any interpolation method already defined (for instance, linear, cubic)\n", + " - faster if many points\n", + "- cons:\n", + " - monotonicity is preserved but not concavity/convexity (it does not look good, when trying to interpolate asymptotically linear functions)\n", + " - extrapolation would need to be defined separately\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAI/CAYAAADz4aFLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdZ3gV5cKF4Wfv9E4ISSihCwnSQxURKQoWRKXjQUAUFBREPXwioqJUaTZQEEEQ6dWGqCAREBEIPZSEXlMI6T3Z8/2I5oiAUpJMQtb9R/bsKWtekCuLmXnHYhgGIiIiIiIiYi6r2QFERERERERE5UxERERERKRIUDkTEREREREpAlTOREREREREigCVMxERERERkSJA5UxERERERKQIsC/Mg5UpU8aoUqVKYR7yuqSkpODm5mZ2jBJJY28ejb15NPbm0dibR2NvLo2/eTT25imqYx8aGnrRMAzfq31XqOWsSpUq7Ny5szAPeV1CQkJo3bq12TFKJI29eTT25tHYm0djbx6Nvbk0/ubR2JunqI69xWI5da3vdFujiIiIiIhIEaByJiIiIiIiUgSonImIiIiIiBQBhfrM2dVkZWVx9uxZ0tPTTcvg5eXFoUOHTDt+cePs7ExAQAAODg5mRxERERERuW2YXs7Onj2Lh4cHVapUwWKxmJIhKSkJDw8PU45d3BiGQWxsLGfPnqVq1apmxxERERERuW2Yfltjeno6Pj4+phUzuTEWiwUfHx9Tr3SKiIiIiNyOTC9ngIpZMaPfLxERERGR/FckyllRMnr0aKZMmXLN79esWcPBgwcLMZGIiIiIiJQEKmc3SOVMREREREQKQrErZ2t2n+PuiT9TdcR33D3xZ9bsPnfL+5w8eTKBgYHcd999HDlyBIDZs2fTpEkT6tevT5cuXUhNTWXr1q18/fXXDB8+nAYNGnDs2LGrriciIiIiInKjilU5W7P7HK+t2s+5+DQM4Fx8Gq+t2n9LBS00NJSVK1eye/duVq1axY4dOwDo3LkzO3bsYO/evdSqVYs5c+bQokULOnXqxOTJk9mzZw/Vq1e/6noiIiIiIiI3qliVs8k/HCEtK+eyZWlZOUz+4chN73Pz5s107NgRV1dXPD096dSpEwAHDhzgnnvuoW7duixcuJCwsLCrbn+964mIiIiIiPwT099zdiPOx6fd0PLrdbXZB/v168eaNWuoX78+8+bNIyQk5KrbXu96IiIiIiIi/6RYXTkrX8rlhpZfj1atWvHtt9+SlpZGUlIS33zzDZD7Yupy5cqRlZXFwoUL89b38PAgKSkp7/O11hMREREREbkRxaqcDe8QiIuD3WXLXBzsGN4h8Kb3GRwcTOfOnWnQoAFdunThnnvuAWDMmDE0a9aM+++/n6CgoLz1e/bsyeTJk2nYsCHHjh275noiIiIiIiI3oljd1vhYwwpA7rNn5+PTKF/KheEdAvOW36zhw4fzzjvvXLF80KBBVyy7++67L5tKf9CgQVddT0RERERE5EYUq3IGuQXtVsuYiIiIiIhIUVOsbmsUERERERG5XamciYiIiIiIFAEqZyIiIiIiIkWAypmIiIiIiEgRoHImIiIiIiK3FcNmA5vN7Bg3TOWsCNu9ezfPPPMMAPPmzeOFF14AYObMmXzxxRcFfvz77ruPuLi4Aj+OiIiIiEh+Sdu7l5Pde+CyZYvZUW6YylkByc7OvuV9jB8/niFDhlyx/LnnnqNPnz63vP9rMQwDm83Gk08+yccff1xgxxERERERyS/ZFy9y/rWRnOzRk+zoaGzu7mZHumEqZ0CvXr1o1KgRtWvX5tNPP81b7u7uziuvvEJwcDDt2rUjJiYGgNatWzNs2DBatGhBnTp12L59OwCjR49m4MCBtG/fnj59+pCens5TTz1F3bp1adiwIRs3bgRg2rRp9O/fH4D9+/dTp04dUlNTL8uUlJTEvn37qF+//hV5R48ezZQpU/KyvPrqqzRt2pSaNWuyefNmAHJychg+fDhNmjShXr16zJo1C4Dk5GTatWtHcHAwdevW5auvvgLg5MmT1KpVi8GDBxMcHMyZM2fo1KkTixcvzrdxFhERERHJb0ZWFpfmz+fYAw+S8O23+Ax4hmpr15IRHGx2tBumcgbMmDGD0NBQdu7cyYcffkhsbCwAKSkpBAcHs2vXLu69917efvvtvG1SUlLYunUrH3/8cV7RAggNDeWrr75i0aJFzJgxA8gtYIsXL6Zv376kp6czbNgwjh49yurVq3nqqaeYNWsWrq6ul2XauXMnderUua782dnZbN++nffffz8v45w5c/Dy8mLHjh3s2LGD2bNnc+LECZydnVm9ejW7du1i48aNvPLKKxiGAcCRI0fo06cPu3fvpnLlynh7e5ORkZE3HiIiIiIiRUnKtm0cf/xxoiZMxKVhQ6p99RV+r7yCnbub2dFuir3ZAS7z/QiI3J+/+yxbFx6c+I+rzJw5k7Vr1wJw5swZIiIi8PHxwWq10qNHDwB69+5N586d87bp1asXAK1atSIxMZH4+HgAOnXqhIuLCwBbtmzJuy0xKCiIypUrEx4eTr169Zg3bx716tXj2Wef5e67774i04ULF/D19b2uU/wzV6NGjTh58iQAP/74I/v27WPFihUAJCQkEBERQUBAACNHjmTTpk1YrVbOnTtHVFQUAJUrV6Z58+aX7dvPz4/z58/j4+NzXVlERERERApa1vnzRL07iaQffsAhIICAj2fg3qYNFovF7Gi3pGiVMxOEhIQQEhLCb7/9hqurK61btyY9Pf2q6/71N/vvv/F/fnZz+19L//OK1NVERETg7u7O+fPnr/q9i4vLNXP8nZOTEwB2dnZ5z7oZhsFHH31Ehw4dLlt33rx5xMTEEBoaioODA1WqVMk7zl+z/yk9PT2vbIqIiIiImMmWkUHsnDnEfjobAN8Xh1K6f3+sf/w8XNwVrXL2L1e4CkJCQgKlSpXC1dWVw4cPs23btrzvbDYbK1asoGfPnixatIiWLVvmfbd06VLatGnDli1b8PLywsvL64p9t2rVioULF9K2bVvCw8M5ffo0gYGBJCQk8OKLL7Jp0yZeeOEFVqxYQdeuXS/btlatWkydOvWmz6tDhw588skntG3bFgcHB8LDw6lQoQIJCQn4+fnh4ODAxo0bOXXq1DX3YRgGkZGRVKlS5aZziIiIiIjcKsMwSN64kajxE8g6exaPBx7A//+G41C+vNnR8lXRKmcmeOCBB5g+fTr16tUjMDDwstv63NzcCAsLo1GjRnh5ebF06dK877y9vWnRogWJiYnMnTv3qvsePHgwzz33HHXr1sXe3p558+bh5OTEoEGDGDx4MDVr1mTOnDm0adOGVq1a4efnl7dtUFAQCQkJJCUl4eHhccPn9cwzz3Dy5EmCg4MxDANfX1/WrFnDf/7zHx555BEaN25MgwYNCAoKuuY+QkNDad68Ofb2Jf6PiYiIiIiYJOPECaLGTyBl82Yc76hOpXmf4/a3R3FuF5Z/uvUuvzVu3NjYuXPnZcsOHTpErVq1Ci3D1VyrALm7u5OcnHzF8tatWzNlyhQaN25coLnee+89PDw88t51VthefPFFOnXqRLt27a74Lr9+30JCQmjduvUt70dunMbePBp782jszaOxN5fG3zwa+5uXk5xC7MxPiJ3/BVYnJ3yHvID3E09gcXC4ru2L6thbLJZQwzCuWiQ0W2MRNmjQoLznycxQp06dqxYzEREREZGCYhgGCd98w/EHHyT2szl4PfII1dd9T+m+fa+7mBVXul/tH1ztqhnktvDC4OzszJNPPlkox7qaAQMGmHZsERERESl50g8dInLsONJCQ3GuU4eAjz7EpUEDs2MVGpUzERERERExVU58PDEffkjckqXYeXlRdsw7lOrSBYu1ZN3op3ImIiIiIiKmMHJyiF++gpj33ycnMRHvJ57Ad8gL2F1lJvSSQOVMREREREQKXequ3USNHUv6wYO4NmmC/6hROAfWNDuWqVTORERERESk0GRFRxMzdSoJX32Nvb8/FaZNxePBB7FYLGZHM13JuomzmNm9e/dNT6Pfr18/VqxYccXynTt3MnTo0H/c9r777iMuLu6mjisiIiIicjVGZiaxcz/n+IMPkbj2e3wGDqT62u/wfOghFbM/qJwVkOzs7Fvex/jx4xkyZEg+pPmfxo0b8+GHH/7jOk8++SQff/xxvh5XREREREqu5F9/5fhjjxM9aRKujRtT7dtv8Hv5JaxubmZHK1JUzoBevXrRqFEjateuzaeffpq33N3dnVdeeYXg4GDatWtHTEwMkPsS6mHDhtGiRQvq1KnD9u3bARg9ejQDBw6kffv29OnTh/T0dJ566inq1q1Lw4YN2bhxIwDTpk2jf//+AOzfv586deqQmpp6WaakpCT27dtH/fr1gdxp/f/cV7169Vi5cmVexj+tWLGCfv365X1ev34999xzDzVr1uTbb78Fcl8D0LFjx3/cZ6dOnVi8eHH+DK6IiIiIlFiZZ89xdsgQzjz9DEZ2NgGffEzFWTNxrFzZ7GhFkp45A2bMmEHlypVJS0ujSZMmdOnSBR8fH1JSUggODmbq1Km88847vP3220yfPh2AlJQUtm7dyqZNm+jfvz8HDhwAIDQ0lC1btuDi4sLUqVOB3AJ2+PBh2rdvT3h4OMOGDaN169asXr2acePGMWvWLFxdXS/LtHPnTurUqZP3ecyYMXh5ebF//36A67rt8OTJk/zyyy8cO3aMNm3acPTo0cu+v9Y+vb29ycjIIDY2Fh8fn5sZUhEREREpwWzp6cTO/ozYzz4DqxXfYcMo/VQ/rE5OZkcr0opUOXt3+7scvnQ4X/cZVDqIV5u++o/rzJw5k7Vr1wJw5swZIiIi8PHxwWq10qNHDwB69+5N586d87bp1asXAK1atSIxMZH4+Hgg96qTi4sLAFu2bMm7LTEoKIjKlSsTHh5OvXr1mDdvHvXq1ePZZ5/l7rvvviLThQsX8PX1zfu8fv16lixZkvfZ29v7X8+9e/fuWK1WatSoQbVq1Th8+PKx/ad9+vn5cf78eZUzEREREbluhmGQtH490RMmknX+PJ4PPYjf8OE4lCtndrRioUiVMzOEhIQQEhLCb7/9hqurK61btyY9Pf2q6/71QcW/P7T452e3v9w3axjGNY8bERGBu7s758+fv+r3Li4ul+UwDOOqD0r+ddnfc18r47/t8899/VkyRURERET+Tcbx40SNHUfK1q041ahBpfnzcWvW1OxYxUqRKmf/doWrICQkJFCqVClcXV05fPgw27Zty/vOZrOxYsUKevbsyaJFi2jZsmXed0uXLqVNmzZs2bIFLy8vvK7yorxWrVqxcOFC2rZtS3h4OKdPnyYwMJCEhARefPFFNm3axAsvvMCKFSvo2rXrZdvWqlUr77ZIgPbt2zN9+nTef/99IPcWRG9vb/z9/Tl06BCBgYGsXr0aDw+PvG2WL19O3759OXHiBMePHycwMPCy87vWPg3DIDIykipVqtza4IqIiIjIbS8nOZmLMz7m0oIFWF1c8H/9dbx79cRiX6SqRrFQ4icEeeCBB8jOzqZevXq88cYbNG/ePO87Nzc3wsLCaNSoET///DNvvvlm3nfe3t60aNGC5557jjlz5lx134MHDyYnJ4e6devSo0cP5s2bh5OTEy+99BKDBw+mZs2azJkzhxEjRhAdHX3ZtkFBQSQkJJCUlATAqFGjiIuLo06dOtSvXz9vcpGJEyfSsWNH2rZtS7m/XS4ODAzk3nvv5cEHH2TmzJk4Oztf9v219hkaGkrz5s2x1/9QIiIiInINhs1G/Jo1HHvwQS7Nm4fXY49Sfd33lH6yt4rZTSrxo+bk5MSqVasuu+L0V2PGjGHMmDFXLO/SpQsTJky4bNno0aMv++zs7My8efOu2Hbu3Ll5v65YseIVE3X8qX///ixdupRnnnkGd3d35s+ff8U6Xbt2veKqG3DV40LuTJOtW7cGuOY+FyxYwODBg6+6vYiIiIhIWlgYUWPGkrZnD8716lHx449xqVvX7FjFXom/claUDRo0CCcTZrSpU6cO7dq1K/TjioiIiEjRlh0Xx4W3RnOyazcyz5yh3LhxVFmyWMUsn5T4K2f/JDk5+arLQ0JCCuX4zs7OPPnkk4VyrL8aMGBAoR9TRERERIouIyeHuKVLifngQ2zJyZTu8yRlnn8eO09Ps6NdIT07nR9P/Uh8RrzZUW6YypmIiIiIiFxTamgokWPGknH4MK7NmlF21Os41ahhdqwrHI07yvLw5Xxz/BuSMpNo7dGaPvQxO9YNUTkTEREREZErZEVFEz1lConffIN9uXJUeP89PDp0uOarmMyQkZPBT6d+YvmR5eyK3oWD1YH7Kt9Ht5rdSD509bvgijKVMxERERERyWNkZnLpiy+4+PEnGNnZ+Ax6jjIDBmB1dTU7Wp5TiadYfmQ5Xx37iviMeCp6VOSVRq/w6B2P4u3sDUDI4RBzQ94ElTMREREREQEgefNmosaNJ/PkSdzbtMH/tRE4VqpkdiwAsmxZbDy9kWXhy/j9wu/YWexoW6kt3Wp2o1m5ZlgtxX+uwxJfzuLj45k7dy4vv/yy2VFEREREREyReeYMURMmkvzzzzhWrkzFT2fh3qqV2bEAOJ98nhXhK1h9dDUX0y5Szq0cQxoO4fE7HsfX1ffqG13Yi3NaZOEGzQcqZ/HxfPbZZ1eUs5ycHOzs7ExKJSIiIiJS8GxpacTOnk3sZ3PA3h7fV16mdN++WB0dTc2VY8th87nNLDuyjC3ntmCxWGhVoRXdArtxd/m7sbNe5ef0lFjYvwx2L4So/QRUeBjoWejZb0WJL2cjRozgxIkTNGjQAAcHB9zd3SlXrhx79uxh7dq1dOzYkQMHDgAwZcoUkpOTGT16NMeOHeP5558nJiYGV1dXZs+eTVBQkMlnIyIiIiLy7wzDIOmHH4ma9C7Z5y/g2bEjfsP/i4O/v6m5olKiWHV0FasiVhGZEomviy8D6w2kS40ulHMvd+UGOdlwdD3s+RKOrANbFpRvCA9N4WRCOQIK/xRuSYkvZxMnTmTfvn3s2bOHkJAQHn74YQ4cOEDVqlU5efLkNbcbOHAgM2fOpEaNGvz+++8MHjyYn3/+ufCCi4iIiIjchIyjR4kcO47UbdtwCgykwpeTcG3c2LQ8NsPGtvPbWBa+jJAzIeQYObQo34IRTUbQqmIrHKwOV24UE55byPYugeQocC0DzZ6FBk+Af20Asgvp3cT5qUiVs8jx48k4dDhf9+lUK4iyI0de9/pNmzalatWq/7hOcnIyW7dupVu3bnnLMjIybjqjiIiIiEhBy0lK4uL06Vz6ciFWd3f83xiFd48eWOzNqQSxabGsObqGFeErOJt8Fm8nb/rU7kO3Gt2o6Fnxyg3SE+DAKtizEM7uAIsd1OwADf6T+1+7q5S4YqZIlbOiwM3NLe/X9vb22Gy2vM/p6ekA2Gw2SpUqxZ49ewo9n4iIiIjIjTBsNhJWryF62jRyLl2iVLdu+L40DHtv78LPYhjsidnD4sOL+enUT2Tbsmns35ihwUNpV6kdjnZ/e9bNZoOTm3ML2cGvITsNfGtB+7FQrwe4+xX6ORSkIlXObuQKV37x8PAgOfnqL6jz9/cnOjqa2NhY3N3d+fbbb3nggQfw9PSkatWqLF++nG7dumEYBvv27aN+/fqFnF5ERERE5NrS9h8gcuwY0vfuw6VBA/xnzcKlTu1Cz5Galcp3J75j6eGlHIk7goeDBz0Ce9C9Zneqlap25Qbxp2HPotxSFn8anLygQS9o2BvKB0MRehF2fipS5cwMPj4+NGvWjDp16uDi4oL/Xx6CdHBw4M0336RZs2ZUrVr1sgk/Fi5cyKBBgxg7dixZWVn07NlT5UxEREREioTsS5eIee894lesxM7Hh3ITJ+DVqRMWa+G+C+xEwgmWHVnGV0e/IikriUDvQN666y0eqvoQrg5/e6l1dgYcWQu7voBjG3OXVWsN7d6CoIfBwaVQs5uhxJczgLlz5+Lh4XHV74YOHcrQoUOvWF61alXWrVtX0NFERERERK6bkZ1N3JKlxHz4IbbUVEr37UuZF57Hzt290DJk27L55ewvLDm8hG0XtmFvtef+yvfTK6gXDXwbYPn7Va/oQ7BrAexdDGmXwKsi3PsqNPwPlCoaL8AuLCpnIiIiIiK3gZTt24kaO46M8HDcWtyF/+uv41S9eqEdPzYtllURq1gWvozIlEj8Xf0Z0nAInWt0poxLmctXzkjKndxj94LcyT2sDrlXx4KfhGpt4GrvMSsBVM5ERERERIqxrMhIoidNJnHtWhzKl6fChx/gcf/9V16hKgCGYbA3Zi+LDy/mx1M/km3Lplm5ZoxoMoJ7K96LvdX+ryvnFrFdX+QWs6wU8A2CDuNzJ/dwK3PtA5UQKmciIiIiIsWQLTOTS5/P4+LMmZCTQ5nBg/EZ8AxWl4J/Nis1K5W1J9ay9MhSDl86jLuDe+4EH4Hdqeb1twk+Ui7mvo9s9wKIOQwOblCnMwT3gYAmt+3kHjfjX8uZxWKZC3QEog3DqPPHstLAUqAKcBLobhhG3M2GMAyjUJq95A/DMMyOICIiIlKiJYWEEDVhAlmnTuN+Xzv8R4zAMSCgwI97KvEUSw4vyZvgo4Z3Dd5o/gYdq3W8fIIPWw4c35h7lezwWrBl5RaxTh9B7cfB6erzPZR013PlbB4wHfjiL8tGABsMw5hosVhG/PH51ZsJ4OzsTGxsLD4+PipoxYBhGMTGxuLs7Gx2FBEREZESJ/PUKaImTCQ5JATHqlWp+NlnuLe8u0CPmWPLYdPZTSw5soSt57dib8md4KNHUA+C/YIv/xk+/jTsXgi7v4TEs+BSGpoOzH2WzK9Wgea8HfxrOTMMY5PFYqnyt8WPAq3/+PV8IISbLGcBAQGcPXuWmJiYm9k8X6Snp6ts3ABnZ2cCCuFfZkREREQkly01lYuzPuXS3LlYHBzwGz6c0k/2xuLo+O8b36SEjATWHF3D4sOLOZd8Dj9XP55v8Dxda3a9fIKP7Ew48t3lU+BXbwMdxkLgQ2DvVGAZbzc3+8yZv2EYFwAMw7hgsVhu+tXcDg4OVK1a9WY3zxchISE0bNjQ1AwiIiIiIn9nGAZJ339P1KTJZEdG4vVoJ3xfeQUHv5v+8ftfHY07yqLDi/j2+LekZacR7BfMS41eom2ltjhYHf63YuwxCJ2X+7Lo1IvgGVBip8DPL5breX7ojytn3/7lmbN4wzBK/eX7OMMwvK+x7UBgIIC/v3+jJUuW5EPs/JWcnIx7Ib77Qf5HY28ejb15NPbm0dibR2NvLo2/eW5l7O3PncNj6TIcw8PJqhhAUo+eZN1RMFPj2wwbYWlhhCSFEJ4ejj32NHZrzL2e9xLg+L87piy2LMpc/J3y53/AO34fBlYulmnKhXIduFS6PliKzhT4RfXPfZs2bUINw2h8te9utpwdAVr/cdWsHBBiGEbgv+2ncePGxs6dO28ke6EICQmhdevWZscokTT25tHYm0djbx6NvXk09ubS+JvnZsY+JzGRmA8/Im7xYuzc3fF9aRilunXDYpf/xScxM5HVEavzbl30d/WnZ1BPutTogrfzX669xB6DXfNznydLvQhelaBRH2j4JHiUzfdc+aGo/rm3WCzXLGc3e1vj10BfYOIf//3qJvcjIiIiIiKAYbORsGoV0dPeIyc+nlI9uuM7dCj23le9Qe2WHIs/xuLDi/n62NfXvnUxOxMOf5t76+KJX3KvigU+CI2fgmptwWrN91wl3fVMpb+Y3Mk/ylgslrPAW+SWsmUWi+Vp4DTQrSBDioiIiIjcztL27SNyzFjS9+/HJTiYsqNex/nOO/P1GDm2HDaf28zCQwvZdmEbjlZHHqr2EE8EPUEtn7/MpHi1q2RtR0GD3uBZLl8zyeWuZ7bGXtf4ql0+ZxERERERKVGyY2OJnjaNhJWrsPf1pfzkSXh27Jivr5hKzExkTUTurItnk8/i5+rHi8Ev0rlGZ0o7l/4jyB8zLu78/PKrZI2eyp150Vp0niW7nd3sbY0iIiIiInKTjKws4hYvJuaj6djS0yn9dH/KDBqMnbtbvh3jePxxFh1edNmti8MaDbv81sVLxyF0PuxZCCkx4FUR2oyChrpKZgaVMxERERGRQpSy7Xeixo0lI+Iobi1b4j9yJE7V8ufVUjbDxq/nfuXLQ1+y9fxWHK2OPFj1QZ6o9QR3+vxxm2R2Jhxanfss2fGQv1wl6wfV2+oqmYlUzkRERERECkHW+fNETZpM0rp1OAQEEDBjOu5t2+bLLYxp2Wl8c+wbvjz0JScSTuDr4suQhkPoWrPr/25djD+dW8h2LYCUaF0lK4JUzkRERERECpAtI4NLc+dycdanYBiUGfICPk8/jdXZ+Zb3HZUSxZIjS1gevpyEjATu9LmTCfdMoEPlDjjYOYDNBhE/wY45EPFD7kY1H4DG/XWVrAhSORMRERERKQCGYeC4bx/Hx40n68wZPNq3x//V/8OhQoVb3ndYbBgLDi7ghxM/kGPk0LZSW56880mC/YJzr8SlxMLuBRD6OcSdBDc/aPly7q2LpSre8vGlYKiciYiIiIjks4wTJ4iaMAHvTZuxVK9OpblzcGvR4pb2mWPLIeRMCF8c/IJd0btwtXelZ1BPnqj1BBU9KoJhwJntsHMOhK2BnAyo3BLavQlBj4C9Yz6dnRQUlTMRERERkXxiS0nh4syZxM6bj9XRkaSuXWny1ptYHBxuep8pWSmsjljNwkMLOZt8lgruFRjeeDiP13gcD0cPyEiGnXNhx1yI2g9OntCob+6ti361/v0AUmSonImIiIiI3CLDMEj89juiJ08mOzoar8cew++Vl9kSFnbTxexc8jkWHVrEqohVJGcl09CvIS83fpk2Fdtgb7WH6MO5V8n2LIbMJPCvCx3fh7rdwMk9n89QCoPKmYiIiIjILUg/fJioseNI3bkT59q1qfDB+7g2bHhT+zIMgz0xe1hwcAEbTm/AipX7q9xPnzv7UKdMndxp8A9+lTvBx6lfwc4RaneGJk9DQBPIx5dXS+FTORMRERERuQk58fHEfPgRcUuWYOfpSRQk3NIAACAASURBVNl33qZUly5Y7G58BsRsWzY/nfqJL8K+4EDsATwdPelXux+9gnpR1q0sxJ+BDWNg1xe50+CXqgz3vwMNeoObTwGcnZhB5UxERERE5AYYOTnEr1hJzHvvkZOYiHevXvgOHYKdl9cN7yslK4VVEatYcHABF1IuUMWzCqOajeKR6o/gau+S+5Lo7S9B+LrcDWp0gCbP/DENvjV/T0xMp3ImIiIiInKdUnfvJmrsONLDwnBp3Iiyo0bhHBR0w/uJSoli0eFFLD+ynKSsJBr5N2Jks5G0CmiFNTMFdn0JO2bDxXBwLQMtX/pjGvxK+X9SUmSonImIiIiI/IvsmBiip04jYc0a7P38KD9lCp4PP5T7TrEbEB4Xzvyw+aw9sRabYeP+yvfT986+1PWtCxcj4PsRsGdR7gQf5YPh8VlQ+3GwdyqgM5OiROVMREREROQajKwsLn25kIvTp2PLzMRnwADKPPcsVje369+HYfB75O/MOzCPX8//iou9Cz0Ce9C7Vm8C3MpBxI+w7k049vMfE3w8Dk2fhYBGBXhmUhSpnImIiIiIXEXK1q1EjhtP5rFjuLW6h7IjR+JYpcp1b59ly2JH8g5mfDuDw5cO4+Psw9CGQ+ke2B2vnBzY/SXs+AziT4FHOWgzKvf9ZO5+BXdSUqSpnImIiIiI/EXWuXNETXyXpJ9+wqFiRQI+/hj3Nq2v+xbG5MxkVkasZMHBBUSlRlHNqxrvtHiHh6s9jGNMOPwwCvYth+w0qHw33P82BHUEu5t/UbXcHlTOREREREQAW3o6sZ/NIXb2bLBY8B32IqWfegqr0/U97xWZEsnCQwtZEb6C5KxkmpRtwuPujzPo/mewHlkLXzyW+24yexeo1w2aDoSydQv4rKQ4UTkTERERkRLNMAySN2wgasJEss6dw+PBB/D/v//DoVy569r+yKUjzA+bz/cnvsfAoH3l9vSt3Zfazn6cWPEW1g8aQNL53JkW7x8DDXuDa+kCPispjlTORERERKTEyjh+nKhx40n59VecatSg0rx5uDVv9q/b/TnJx9z9c/ntwm+42LvQM6gnve/sTYWESNj0IYStompOJlRrAx2nQY32YL3xF1RLyaFyJiIiIiIlTk5yChc//phLX3yB1cUF/5Ej8X6iFxb7f/7xOMeWw4bTG5hzYA4HYw9SxqUMLwa/SLfqnfE68QssfQrO/A6O7tCoH9uN+jR9uHchnZUUdypnIiIiIlJiGIZB4jffED15CtkxMXh16Yzfyy9j7+Pzj9tl5mTy9bGvmRc2j1OJp6jkUYm37nqLR8q3xGnPEph5DySeBe8q8MBEaPAfcPYkNSSkUM5Lbg8qZyIiIiJSIqQfPEjk2HGk7dqFc926BMyYjku9ev+4TVJmEsvDl7Pg4AIupl3kTp87mXrvVNq5VsRu+2xYPgSyUqHKPfDQZKjZQbcuyk1TORMRERGR21p2XBwxH3xA/NJl2Hl7U27sGLw6d8ZitV5zm4tpF1lwcAHLjiwjOSuZ5uWaM6HleJqlpGDZMhOO/pT7wui63aH5c5p1UfKFypmIiIiI3JaMnBzily0j5v0PyElOxrt3b3yHvICdp+c1tzmdeJrPwz7n66Nfk2XL4v7K99M/6Alqn90Lq16Ei0fAzQ/avA6NngJ330I8I7ndqZyJiIiIyG0nddcuIseMJePQIVybNsX/9ddxDqx5zfXDYsOYu38u60+vx85ix6N3PEq/Sg9S+dBamN8F0uKgXH14fBbUfhzsr+/dZyI3QuVMRERERG4bWdHRRE+ZQuLX32BftiwV3puGxwMPYLFYrlj3z+nw5+yfw7YL23B3cKdf7X709qqD757F8OMDgAFBHaH5YKjUHK6yH5H8onImIiIiIsWekZnJpQULuDjjY4ysLHyee5YyAwdidXW9Yt0/p8Ofe2AuYbFh+Dj7MKzhULrbXPDY8TmcGw1OXtB8EDQdCN6VC/+EpERSORMRERGRYi158xaixo8n88QJ3Fu3xv+1EThWvrJQZdmy+O74d8zZP4eTiSep6FGRNxsPp9OlGJzWvweJ56B0dXhoCtTvBU7uJpyNlGQqZyIiIiJSLGWePUvUhIkkb9iAQ+VKVJw1E/d7771ivYycDFZHrObzA59zPuU8Nb1rMrnR/3H/mTDsvnoNMpNzp8J/eBrUaA//MIujSEFSORMRERGRYsWWlkbs7M+I/ewzsLfH9+WXKd2vL1ZHx8vWS81KZdmRZcw/OJ+LaRep51uPkXd0p1X4JiwrX8x9fqx2Z7jreSjfwKSzEfkflTMRERERKRYMwyDpx5+Ienci2ecv4Pnww/gN/y8OZctetl5CRgKLDi9i4aGFJGQk0KxsUyZWfoymYWuxbH8593myu56HZs+BVwWTzkbkSipnIiIiIlLkZRw9SuS4caT+tg2nmjUp/8VE3Jo2vWydP18cvfTIUlKyUri3fEsGOJan/t7VcGkFeFWCDhMg+Elw8jDpTESuTeVMRERERIqsnORkLk6fwaUvv8Tq6or/qFF49+yBxf5/P8ZGpkTy+YHPWRmxksycTNpXaMWAbCcCd62CtEtQoRF0/RxqdQI7/fgrRZf+dIqIiIhIkWPYbCR89TXRU6eSExtLqa5d8X1pGPalS+etczrxNHMOzOHrY1+DAR3Lt6R/YjJVf1sKOVkQ9DDc9YLeTybFhsqZiIiIiBQpaQfCiBozhrS9e3GuX4+Kn3yCS906ed9HxEUwe/9sfjj5A/YWe7r638VTUWcpv+ULsHeB4D65L432qW7iWYjcOJUzERERESkSsuPiiJn2HvErVmBXujTlxo/H67FHsfwxtX3YxTBm7ZvFxjMbcbV3pa9vM/qcPECZYwvBzQ/ajoLGT4Nr6X85kkjRpHImIiIiIqYysrOJW7qUmA8+xJaSQuk+fSjzwvPYeeRO2rEvZh8z985k87nNeDp6MsinCf85ugOviKXgGwSdpkO97mDvZPKZiNwalTMRERERMU3qjh1Ejh1HxpEjuN7VnLKvv47THXcAsCd6DzP3zuTX879SytGTF73q0TP8V9xTD0Clu+DBSVCjg14aLbcNlTMRERERKXRZUVFET5pM4nffYV++HBU++ACP9vdjsVgIjQpl5t6ZbLuwDW9HT15yrUmP8K24ZR2AwIfh7hehUjOzT0Ek36mciYiIiEihsWVmcmnefC7OnAnZ2ZQZPAifAQOwuriwI3IHM/fOZHvkdko7evJfx4p0i9iGK0egfg9oMRR8A80+BZECo3ImIiIiIoUi+ZdfiBo/gcxTp3Bv1w7/Ea/iEBDA9sjtfPLLJ4RGhVLGwZP/w4eu4XtwcXCH5s9D80HgWd7s+CIFTuVMRERERApU5unTRE2YSPLGjThWqULF2Z/i1rIlv134jZnrXmd39G78HDwYkeVKl5NhOLv5Qbu3oNFT4FLK7PgihUblTEREREQKhC01lYuffsqluZ9jsbfHb/h/8e7dm19jtjPz+97si9mHv707I1MtdI4Ow6n0HdDxA6jXAxyczY4vUuhUzkREREQkXxmGQdK6dURNmkz2hQt4PvIIvv99hd8yDzPzp34ciD1AOTs33kjM5LHYgzhWaAzdx+VO9qGZF6UEUzkTERERkXyTHh5O1LjxpP7+O061alF+8iRCy6by8o6hHIw9SAU7V96KS+HR+NM41OgAj8yFyi3AYjE7uojpVM5ERERE5JblJCYSM306cQsXYXV3x//NNznUMoBR+6exP2w/FawuvBObQMekszjU7gw9X4KydcyOLVKkqJyJiIiIyE0zbDYSVq8meuo0cuLiKNW9G6d7tWTs8S/Ys3EP5SyOjL4YR6eUCzg0eCL3HWU+1c2OLVIkqZyJiIiIyE1J27ePyLHjSN+3D5eGDbk08SXeTf2W0N9X4YcDb1y8xOPpNhwaPQV3PQ9eFcyOLFKkqZyJiIiIyA3Jjo0l+r33SFixEjvfMmS9Ppipvnv4/cTblMGOEbGX6Jpph1PTwdDsOXDzMTuySLGgciYiIiIi18XIziZu0WJiPvoIW1oatp6P8ElwLBvjPqV0pJXhl+LonuOC813DoXF/cPY0O7JIsaJyJiIiIiL/KuX37USNHUtGRARGk/os6uDMVznfU+qShZfj4uiBN64t34KGvcHBxey4IsWSypmIiIiIXFPWhQtETZpE0vfroJwf6wbUY67PATyzLLwYH08vB3/c2kyGul3BzsHsuCLFmsqZiIiIiFzBlpHBpc8/5+KsTzFsOYQ+XJ1ptU7ibBfL4PgEertUwaPDOAjqqBdHi+QTlTMRERERuUzSxo1ETZhI1unTnGpYlsnNY0jxPEn/hESe9AjEq+P7UL2dXhwtks9UzkREREQEgMyTJ4mcMIGUXzaRWM6D6T3tOFQlhicSkuhvuZNSj74GVe5RKRMpICpnIiIiIiWcLSWFizNnETtvHll2sLydHd8Hp/BYagpTLXXwe2wkVG5hdkyR257KmYiIiEgJZRgGiWvXEvXuJHKio/m1rpUv7oV7SGW1Uz0qth8FFZuYHVOkxFA5ExERESmB0o8c4cKYMaTvDOVUWQuzn7Sjkncac73qc0ebN6B8Q7MjipQ4KmciIiIiJUhOQgJRH3xA/OIlpDgbLHrASlqNDN72qU+d1qOhbF2zI4qUWCpnIiIiIiWAkZPDpRUruDB1EpakVH5qaCGsWQ4Dy9ahSZt3wK+W2RFFSjyVMxEREZHbnP2x4+ybPAbHY+cJD4CfH7XRvUZ9Xmw9BotfoNnxROQPKmciIiIit6nsixcJG/0KPuu3c8kd1j1scHfTunzYbjzWMjXMjicif6NyJiIiInKbMbKyODxjPBmfL8Uuy+CnplD5vhq89eAU7H1rmh1PRK5B5UxERESkmFuz+xyTfzjC+fg0Hk7fQc+dK/G5mE14VTAerkRw1f/Q8uE+ZscUkX+hciYiIiJSjK3ZfY7XVu2nXPIhphyZz53HM4gqBb8+7k23ITPwKt+QkJAQs2OKyHVQORMREREpxj79LoSXTnxI0z2JAGxu6sTygF6kejbhGb2rTKRYUTkTERERKYYyk2PY8NEzjPgmHN8ECLvDjkVBj7HL/i6wgSU+zeyIInKDVM5EREREihEjLZ6flwwjcenvBJ2EqNIW3m/Tih+8HrlsvfKlXMwJKCI3TeVMREREpDjISGL7ulHsW/EjzUKhtAMk921HXLv/sumbw5CVk7eqi4Mdwzvo/WUixY3KmYiIiEhRlpXGoc0TWPfNcppuhruTIaFVEA3HzsTJz58mgMXBIW+2xvKlXBjeIZDHGlYwO7mI3CCVMxEREZGiKDuTc79P58ufP6PGLwYPnIXkKqUpP/sDajVsfNmqjzWsoDImchtQORMREREpSmw5JO6ez/wt0zC22ei4xyDbzZFSb/2XoB7/wWK1mp1QRAqIypmIiIhIUWCzkRW2imVbxnHwYDqPbga3DHDp/jiVX34VOy8vsxOKSAFTORMRERExk2FghP/I+l/eZPWFBB7+2ULvKKBBbaq+PR7nwJpmJxSRQqJyJiIiImKWk1vYu+ENPok/R4Pf7BlywEJOmVKUn/oGng89iMViMTuhiBQilTMRERGRwnY2lDMb3uSjhEPYH3Zh4BY7nGxWvAc8hd9zz2F1czM7oYiYQOVMREREpLBEhZGw4W0+vfg7+2Ld6bfemfKxNlxataT866NwrFzZ7IQiYiKVMxEREZGCFneSzJ/HsvjUOlZaPenysxsjw21YAypQfuYoPFq3NjuhiBQBKmciIiIiBSXlIsYvk/jh4CJmuHvQ5KAn4363YG/ngO+wQZR+qh9WJyezU4pIEaFyJiIiIpLfMpLgtxns3vkJU9ydcL7kxWsrrHjHZ+P50AP4DR+OQ7lyZqcUkSJG5UxEREQkv2RnQug8zm6exHsuNg5YPXj2ayu1jmfhWKM6ZT8YhVuzpmanFJEi6pbKmcVieQl4BjCA/cBThmGk50cwERERkWLDZoMDK0nZOIbPiGe5hyedt1rov9OGvaszvq//H969emKx17+Li8i13fTfEBaLpQIwFLjTMIw0i8WyDOgJzMunbCIiIiJFm2HA0Q3kbHiLr5JP8JGPD3eGefDBZntcEzMp1bULvi+9hH3p0mYnFZFi4Fb/+cYecLFYLFmAK3D+1iOJiIiIFANnQ2H9W+yI3M4kX38yckozcqkTASdTcK4XSNk3RuFSt67ZKUWkGLnpcmYYxjmLxTIFOA2kAT8ahvFjviUTERERKYouRsCGdzgTsZapvn5s9/Lj6Y1O3BWain1pZ/zGjcTr8cewWK1mJxWRYsZiGMbNbWixeAMrgR5APLAcWGEYxpd/W28gMBDA39+/0ZIlS24pcEFITk7G3d3d7BglksbePBp782jszaOxN8/tMPaOGZeocnIJ7pHrmVW6FIvd3Lh/r4Vem8AxI4fU1q1J6fgwhqur2VGvcDuMf3GlsTdPUR37Nm3ahBqG0fhq393KbY33AScMw4gBsFgsq4AWwGXlzDCMT4FPARo3bmy0LoIvWQwJCaEo5ioJNPbm0dibR2NvHo29eYr12Gckw9YPyd75EatcHJhRtSr+pzOY/rUb3mcScG3WjLKjXsepRg2zk15TsR7/Yk5jb57iOPa3Us5OA80tFosrubc1tgN25ksqEREREbPlZMPuBbBxPFtzEplcsRKxiWkM/cGF2qGp2Jdzxf/90Xh06IDFYjE7rYjcBm7lmbPfLRbLCmAXkA3s5o8rZCIiIiLFlmFAxI/w05uciD/K1PJV2WI48USoPQ//YoedLYnSg56jzIABWIvgLYwiUnzd0myNhmG8BbyVT1lEREREzHV+D/w4ioTTW5hZtjJLKgbQ6KTB3BBPXM7H4d62Lf6vjcCxYkWzk4rIbUhvQhQRERGJPw0bxpCzfxkrffz5qFoNnGMzmPyzHxX2nMOxsi/+n07EvVUrs5OKyG1M5UxERERKrrR42DINts1kp7MjE2vW5URqPIO3e9EiJBqLfRxlXnmZ0n37YnV0NDutiNzmVM5ERESk5MnOhJ1z4Zd3icxMZGq1eqzLiuHBcHhrgwf2Mefw7NgRv+H/xcHf3+y0IlJCqJyJiIhIyWEYcPAr2PA26XEn+LxKXeZaUykfncCsLeXwDjuLU1AQZd/7CNfGV30NkYhIgVE5ExERkZLhXCisG4lxZhvry9VkSmA94hIuMXxvRer+chqrexK+b76Bd/fuWOz1I5KIFD79zSMiIiK3t8TzsP5t2LeECE9/3q3Tku3Jp+ge5sHj692wxp+iVLdu+L40DHtvb7PTikgJpnImIiIit6fMVNj6Efz6PgnYmFGnDctST1L75EXmbiqLW/g5XBo0wH/2KFzq1DY7rYiIypmIiIjcZmw2OLAC1o8mJ/EcK2vezUeWeIyY44zfXYmqm49j5+OI38QJeHXqhMVqNTuxiAigciYiIiK3kzPbYd1rcG4nO8vXZmKlykQkneaZiIrc91MmpJ2idL9+lHl+MHbu7manFRG5jMqZiIiIFH/xp2H9aDiwkkjPskxt2IF18YdoecKeNzb44XjyJG4t7sL/9ddxql7d7LQiIlelciYiIiLFV0YybHkPfptOugXm1X+IOSkRlDp3lI92VMX/twgcypfH78MP8Lj/fiwWi9mJRUSuSeVMREREih+bDfYuhg3vYCRHElLrPt61JhJ1aT/DIqrTdN1pLMZpfJ5/Hp9nnsbq4mJ2YhGRf6VyJiIiIsXLqa2wbgRc2MvpCg2YENSILbH7eeiCP31+KoP1XDju97XDf8QIHAMCzE4rInLdVM5ERESkeEg4Bz+9CQdWkOZZgdnNn2BezO8EnIzg098qUyr0GI5Vq+L/2We4t7zb7LQiIjdM5UxERESKtqx0+G06bJ6KYcthfZMnmJx6lEtnNjPiYHXqrj+B1eECZYYPp/STvbE4OpqdWETkpqiciYiISNFkGHDke/jhNYg7yYnA+5ng5cJv0ZvpfLos3X/ywhoTjtejnfB95RUc/PzMTiwicktUzkRERKToiQnPfa7s2AZSfQOZdU9/vjgfQvUjDsz5NQCPA6dwqlWLsh98hGtwsNlpRUTyhcqZiIiIFB3pifDLu/D7TAwHN35o8QyTE/aQfPQn3txbhaBfTmLnbsV39FuU6tYNi52d2YlFRPKNypmIiIiY78+p8dePhpQYjtXvzASnbLaf/4Enjpal00+uWJJOUKpHd3yHDsXe29vsxCIi+U7lTERERMx1NhS+Hw7nQkkJaMQnjR9j4Zn13HnUgc83+eMacQ6X4GDKjnod5zvvNDutiEiBUTkTERERcyRHw/q3Yc+XGO7+rG09hKnRW8k6tI6xuypSbcsJ7H2d8Zs8Cc+OHbFYLGYnFhEpUCpnIiIiUrhysmHHbNg4HrLSCG/aj/FcYs/xNfQ97E+HDQ5YMs9S+un+lBk0GDt3N7MTi4gUCpUzERERKTynt8F3r0DUAZKqtebjyrVYfGodjc868XmID86nz+PWsiX+I0fiVK2q2WlFRAqVypmIiIgUOIfMeFg9CPYuwvAM4Ju2LzPtwkYse79j0o7yVNxxGocAT/xnTMe9bVvdwigiJZLKmYiIiBQcWw7snEvT7W+BLZOIZv0Za1xkf8RyBuz3496NVqxE4TPkBXyefhqrs7PZiUVETKNyJiIiIgXj7E747mW4sJfoUnVZEnwvX574jhYnHPk8xBPHyAt4tG+P/6v/h0OFCmanFRExncqZiIiI5K+UWNgwGnZ9geFRjvXt/ss7p37EZedXTN3qS7kDkThWr07Zz6fidtddZqcVESkyVM5EREQkf9hssGs+bHgb0hM53aQf4x1S2Xl4Kf1+c6PtdrBzTqbMiFcp/Z//YHFwMDuxiEiRonImIiIit+787txZGM+FklG5BXNrNuezY19xzyGY84sLjpcS8Xr8cfxefgl7X1+z04qIFEkqZyIiInLz0uJgwxjYORfcfNnS7lXGx2zB+utypmzyouzRSzjXrs35Zx6iVv/+ZqcVESnSVM5ERETkxhkG7FsGP4yEtEtENu7DJBcbW8MW8cw2d1pst2HvacP3nbcp1aULpzZvNjuxiEiRp3ImIiIiNyb2GHz7Epz4hawKwSy852k+ObqaVj9n8ulmBxxSkvDu9QS+Q4dg5+VldloRkWJD5UxERESuT3YG/PoBbJoC9k6Etn6JsfF7sPy8iEkhrvidycS1cWP83xiFc2Cg2WlFRIodlTMRERH5dyc2514ti40g9s6OTPP155ewFQzY7EzTPTnY+7ngN+UtPB9+CIvFYnZaEZFiSeVMREREri0lFn4cBXsXkVOqEivbvcxHp9bRauU+PvnVin1WBj4DBlDmuWexurmZnVZEpFhTORMREZErGQbsWQg/vgEZiYQ17cfYnEgsG5YxYaMTPlHZuLW6h7IjR+JYpYrZaUVEbgsqZyIiInK5mPDcWxhPbSGxYlP+n737ju6qPvw//ryf7AEEAgQStgpuRVGsE1xoXagdamsdtWq1Kg5UVAQHLrSAqN3WTrXflqIVBScgQ2ZQZtgr7BEgJJDxub8/pJT2p61C4Cb5PB/n9LQkl+R13s0BnudzczO0wzF8+Nkofjgmic5z4qS0bkreS33I7t7NWxglqQYZZ5Ik6XOV5fDxczBuMGFqJm+dehND1nzC6X8eyZBPIDkpoGmvO2hy3XXE0tKiXitJ9Y5xJkmSYPHoz18t27SYRUdcxOMZIcGHI3hsdDI5m6ppcP555N17LyktW0a9VJLqLeNMkqREVr4ZRj0EM/5IWZP2/PzU6/jgs4+49oM4Ry6Kk3pIW1r89EGyTuoa9VJJqveMM0mSElEYwpzh8Pa9hGUb+fD47zC4ZBGn/uFdnpkKSRlZNH/gdhpfdSVBsv9ckKQDwT9tJUlKNFtXwYh7oGgEK/OP5MnDuhKOnkDfMQENtoU0uvwymt91F8m5uVEvlaSEYpxJkpQo4nGY9lt4vz+V1ZX87vjLGVk0gx/8pZBDVsZJO+pIWvbtS8bRR0e9VJISknEmSVIi2LAA3rwdlk9gSvsT+Wks4KQ/f8KjM0JiOTm0GHAPjS69lCAWi3qpJCUs40ySpPqsqgImDIExz7AxLZNBx55P2UczuXssZO4MaPL979Pstp+Q1LBh1EslKeEZZ5Ik1Vcrp8GbtxFfN5u/djyFEcs3cMWgz2i3NiT9hC607NuX9I4do14pSdrFOJMkqb6p2A4fDoBJP2NeozwGtenK8X9byv2zQ2jelIJBD9DgvPMIgiDqpZKkPRhnkiTVJ4tHw5u3sX3LCl48+GS2jF/OzeNWkhaPkXvTDTS96UZimZlRr5QkfQHjTJKk+mDHVnivL+G0V3iveXvejHWi54tLKNgEaaefQqsH+5Latm3UKyVJ/4VxJklSXbfwfXjzDlaUr2Vo3rEc+fZ6bl2wk7BVC1o/2Z/sM86IeqEk6SswziRJqqvKS+DdB6ko/COvNGvPhjkt+f6f1xFLSSH3zltoet31xFJTo14pSfqKjDNJkuqi+aPgH72YVLmZNyoP5rzfldFsK6T2OIs2fR4ipUWLqBdKkr4m40ySpLqkfDOM7MOGWa/zy9Q2HDS6CVcvLaOqfQFtXniCrBNPjHqhJGkvGWeSJNUV896m+q1e/LWqnDULCrhoagVhRhq5D/Si2VXfJ0j2r3VJqsv8U1ySpNqubBO8cy+z5w3nzY0t6TYmiaPKqkm++Dw63NeX5CZNol4oSaoBxpkkSbXZnDfZ9vZd/GFLnILxzem5qoodh7al3WMDyTzqqKjXSZJqkHEmSVJtVLaJ8O3evDv3HZbPyuGMT0MqGqbR+LF7ybv8uwSxWNQLJUk1zDiTJKm2WfAeS9+8jRELQk78pBGtKkKC717I0fc8TFKDBlGvkyTtJ8aZJEm1xc5t7Bx5P38dP4LciRmcuT5g29Htaf/4YDI7dox6nSRpPzPOJEmqfh6DcwAAIABJREFUDZaOY+JrN7P4kyqOm5dBaZMMGgzsw6EXfosgCKJeJ0k6AIwzSZKiVFnO2rfvZ9TwtzlqajJHhQE7fnAxx9/Zn1hGRtTrJEkHkHEmSVJEqlZM5p0XbyB7TCUnbE5mw/HtOW7AC2S16xD1NElSBIwzSZIOtKoKZr56G4teHUOnxQGbm6aSNqQfp/W4LOplkqQIGWeSJO1nwwuLGTiqiFUl5ZyauYiLF/6STtOqaZsUsPna8zjpzqeIpaVFPVOSFDHjTJKk/Wh4YTF9hs2kuqKMW0te5BsfrKLJNlhwdC5n/PQPNGrVPuqJkqRawjiTJGk/+unIOVxW+idOmTmV9iuhuFnAIydcxIo2PbjYMJMk7cE4kyRpfwhDtk36Pb2mP0XHWVCeBv/3jQ680vxG4kEyQUl51AslSbWMcSZJUg2Lz3+f8S/eS9rH5XQqg/GHN2Roh5vZltx89zX5OT4mX5L074wzSZJqysqpLPpDb5aNXEXLtTGWt05lzt23M2heAeWV1bsvy0hJonePThEOlSTVRsaZJEn7at1cSt/sy8S3ZtBqXgoZWTGW3XYRZ9/8BMlJyeTs8bTG/JwMevfoRM/OBVGvliTVMsaZJEl7a/My4u8/weS3RpIyI5MWlSnMPqcDpz/8Es2btd19Wc/OBcaYJOl/Ms4kSfq6StfB2GdZMvJPrJzWgKYbM5l/SCatHurHt7peHPU6SVIdZZxJkvRV7dgCE4ay/f2fMfXTTJovakiYEzD/3ks5/5p+pCalRr1QklSHGWeSJP0vlTtgyq+If/Qc02eHJM9sSE4I0y/syFl9htIit03UCyVJ9YBxJknSl6mugk9fhdFPsrxoE6s+zaVRSZyZR2TR9sF+fO+4i6JeKEmqR4wzSZL+UxjCvLfgg8fYvmQhM+YW0GRxY7bnhqx46HJ6XvmwtzBKkmqccSZJ0p6WfAzv9ye+dBqfLWtNUmEe6UlVTLi8I+f1HkpBjrcwSpL2D+NMkiSAVTPgg0cJF37Ayg35rJ3SmqytVUzpnEX7Pv344dHewihJ2r+MM0lSYtu4CD58HGYPo7S8CbNnHkrDxVtZ16KaTXd/i+9c/iDpyelRr5QkJQDjTJKUmLathTFPwfTfU12VxqxVXUkevwLSt/LBFR25+I4htGncLuqVkqQEYpxJkhJLxXaY8AKMH0JYuZPllaezadQiUrevYPyJ2Rx8T19uPfIigiCIeqkkKcHsU5wFQZAD/Bo4EgiB68MwnFgTwyRJqlHxalqueheevwlK17AtuzvzPt5A9tJ5LGsdY0vfK7j6gnvJSM6IeqkkKUHt6ytnQ4CRYRh+KwiCVCCzBjZJklRzwhAWvAfvPUyn9XOpatyFWauPJm3MLCqy4b3rjqTnj5+jTUOfwihJitZex1kQBA2B04FrAcIwrAAqamaWJEk1YNUMeK8vLBlL2Kg9kzf3IH34XJIqVvHh6Y04+p5Hub3juVGvlCQJ2LdXzjoA64HfBkFwDDANuCMMw+01skySpL1VsgI+fAw+ex0ymlDS+scs/PPHNCyeyacHJVF523XccM4dpCWlRb1UkqTdgjAM9+43BkEX4BPglDAMJwVBMATYGoZh3/+47kbgRoC8vLzjX3vttX2cXPNKS0vJzs6OekZC8uyj49lHx7Pff5IrS2mz/G+0WvkPwiBgeYNzWT5hE81nLmVNDow5rwNHn3I9TVJyo56acPy6j5bnHx3PPjq19ey7d+8+LQzDLl/0vn2JsxbAJ2EYttv169OA+8MwvODLfk+XLl3CqVOn7tXn259Gjx5Nt27dop6RkDz76Hj20fHs94OqCpj6Mox5Gso3Ez/sWyxe3JyyP79BPIwzpnsuXe95kh1Lqj37iPh1Hy3PPzqefXRq69kHQfClcbbXtzWGYbgmCIIVQRB0CsOwCDgLmLO3H0+SpK8tDGHBuzDqAdi4kLDd6WzKuIClL7xM5vptFB6WTOodN3HLaT8mJSmF0UtGR71YkqQvta9Pa7wN+NOuJzUuBq7b90mSJH0F6+Z+HmWLPoTcg9lx2vPMfXk46dOGsKEpzLvrFK76/pM0y2wW9VJJkr6SfYqzMAxnAF/4kpwkSfvF9o0w+onPb2NMa0D1GY+yaNwaKm5/mnhyyIiL8uje62nuKOga9VJJkr6WfX3lTJKkA6OqAqb8CkY/DRWlhF1+yPptR1LcZyjpJWVMPDaVRnfcSq+u15Mc8683SVLd499ekqTaLQxh/ih490HYuBAOOovy9jcw76dDSZ8zgpUtYPGDZ/H9bz1CboZPYZQk1V3GmSSp9lo75/PvK1v8EeQeQtWFr7Dgb2OJP3E3Fenw7rdbcf5PBnJR3rFRL5UkaZ8ZZ5Kk2mf7RvhoAEz7LaQ1JDz3SdYsTGPdTY+QXLaTj05IJ7/XXfTu/D1iQSzqtZIk1QjjTJJUe1RXwZRff/7Aj52lcMINlDY6nwWPDyB98SoWtgkovv+bXHPxwzRKaxT1WkmSapRxJkmqHZaOg7fvhXWzoUN3Krvcy/yfv0Ls3VspbQDvXN2BS28ayOVND496qSRJ+4VxJkmK1pZiePchmD0MGrUhvOwVVkxYQcmVN0FVFe+ensXBd9zP/Ydf5i2MkqR6zTiTJEWjaidMfAHGPgthHM64ny3h8Szq9QgZqzYx8+AYpT/+Ljee25uslKyo10qStN8ZZ5KkA2/+KBh5P2xaDIdeSMVRtzPnp4NJG/97tuTABz8+hu9c9zRtG7aNeqkkSQeMcSZJOnA2Lvr80fjzR0LuIcS//RqL3pnGjgHXExLnnXObcEKvx+jd4cyol0qSdMAZZ5Kk/a9iO3z8HEwYCkmphGc/yoaSNiy/6WEyN5Qy7YgU0m6/mdtOvZmUpJSo10qSFAnjTJK0/4QhzH0TRj4AW1fC0d+lvMO1zH7qSbJmLGRDM5jf+3SuvGoATTOaRr1WkqRIGWeSpP1j0xJ4uzcsfA/yjqL6vKHMe/0dwod+SJgKI3u24uw7BtKj5bFRL5UkqVYwziRJNatqJ4wf8vltjLFkwnOeoHhFNuuuvYu0rTuZeHwG+Xffwx2dr/DR+JIk7cE4kyTVnMWjYcTdsHEhHN6TbW2uZu6TT9GgqJji/IB19/Tkqksf8tH4kiR9AeNMkrTvtq2BUQ/CrL9C4/ZUXfQ7Pvvz26T1u43qDPjg+4dywa3P0b5xh6iXSpJUaxlnkqS9F6+GKb+GDx+Hqh2Ep97L4iVZbP1RX1LLqxh3ciMO7/0IPzm0R9RLJUmq9YwzSdLeKZ4Gb90Jqz+FDt3Z2OJ7LHxmCA2Xb2RxuySq7rie68/t5aPxJUn6iowzSdLXU14CHzwKU1+G7Dx2dhvEp6++R4MxD1PREMbe3JXLbhxIs8xmUS+VJKlOMc4kSV9NGMKc4fD2vVC2gfjxP2JuUQZVtz9HWnWcj8/O46R7n+GmNidGvVSSpDrJOJMk/W9bVn7+FMb5I6HlMaxqdycrnv0tDdeWMrdTGg17384Np1zro/ElSdoHxpkk6cvFq2Hyr+DDxyCMU3bsPcz462QaTx3K9iYBRb3P47KrHyM7NTvqpZIk1XnGmSTpi62ZBf+4HYqnEW97JoVLWpDy0KtkBCHjL+nAWfcM4sxmHaNeKUlSvWGcSZL+XWU5jHkaJgwlTMthSdMfseGX79Ng8zymH5NN2z4PccOxl0S9UpKkesc4kyT9y+LR8I9esHkJJS0vZtaoleTOGcHmvBjLHr2Sb19+P6lJqVGvlCSpXjLOJElQtglGPQif/pmqrHZM2XQODf4ylbRUmPy9Yzj/jkE0b9gy6pWSJNVrxpkkJbIwhNl/h7d7E5aXMIdzKXt9Hg1LZ1N4Ui5HPfgk1xxyWtQrJUlKCMaZJCWqbWthxF0w7y3WxA5nwcRmNF06i7Wtk4k9ditXnXeLj8aXJOkAMs4kKdGEIXz2OrxzHzu27WDyimPJnbyO5CyYcdPpXPDjZ2iY3ijqlZIkJRzjTJISyZZieOtOwqJRTN/YkWD8dhpXrKPwzFZ848FBfKPgyKgXSpKUsIwzSUoEYQjTfw/vPsTSNbByRjty15Yy/+B0cu+/j++d8l2CIIh6pSRJCc04k6T6bvMy+MftbJ39MVPnt6bl3ArCRpXMu+divnnto6Qlp0W9UJIkYZxJUv0Vj8PU31A9sh8TF2SSPSOP3LCCTy86lO73D+G03DZRL5QkSXswziSpPtq0GN74CXOmzKBkeg65JSFzj2xEh76PccUx50S9TpIkfQHjTJLqkzCEqb9h3bD+zCzMJn9pI3Y2jbH8kR/Q8zv3+Gh8SZJqMeNMkuqLLSup+MuPGf/BbHJnNqBxEsy54gTO6T2YhllNol4nSZL+B+NMkuq6MCSc8Wdm/Ophdk5OpUVpGnNOaM6x/Z7j8oO7RL1OkiR9RcaZJNVlpeso/tUPKXprPi2LU9nQIpnKR2/jsvN/5KPxJUmqY4wzSaqjyif+gQmDHidvZozs9Bjzr+/GuXc8S0ZaVtTTJEnSXjDOJKkOGF5YzMBRRawqKefQhhXcvWYQjT5eR4sdMeae3IKuj7zIia0Oj3qmJEnaB8aZJNVywwuL6TNsJuWV1Xyr9A0umDCWFusCFrdKpulDD/CtbldFPVGSJNUA40ySarmBo4rovH0K3130Oh2LqtmcHfDrMzoxsd0tTOh2btTzJElSDTHOJKkWixd/xh1zenPwjDJSqmD0sdn8vM1NbIm1JNhSGfU8SZJUg4wzSaqNNi9l3i9vZ92b8zlqY8Dstkn84sjLWZBy4u5L8nMyIhwoSZJqmnEmSbXJtrVs/HtfCl8fTcGiJMiJMfaGcxlccg7lVf+6LCMlid49OkW3U5Ik1TjjTJJqg/ISKj96jvGvvkrjwhSakcScSw7jrAd+xumN8sjb42mN+TkZ9O7RiZ6dC6JeLUmSapBxJklRqigjnPRzZv7fi2yblELelhTmHdWIw/oP5PIjTtt9Wc/OBcaYJEn1nHEmSVGoroTCP7DmjaeZNSmkYHkqO5smsfqJG+h56R0EQRD1QkmSdIAZZ5J0IMXjNF87lp2D7mDclG00n5lG4xSY9/1vcO5dg8jKbBT1QkmSFBHjTJIOhDCEhR8Qvt+PsunLKCxsQH5pGnO6tuD4fj/l+A6do14oSZIiZpxJ0v62fBJ88AjLZ01l4YxcWhY3YEVBKtUDenF5j+uiXidJkmoJ40yS9pf18+H9/myf+Q4TF+TRclYuWRkw8bJjuar/b0hPzYx6oSRJqkWMM0mqaaXrYPRTxKe8wqTiJiRPa0nLHSFzz2jLyf2GUl5UbJhJkqT/j3EmSTWlogw+eRHGDWbhelg5o4C81VUsaZ9O8wcf4FunfguA+UXFEQ+VJEm1kXEmSfsqXg2fvgofDmDL+rVMXtSWVjN3kNqgmiV39eScHz5KSlJK1CslSVItZ5xJ0r5Y+AG89zDVq2cxft1BZE/Io0XVDuZ881BOf/B5muW2jnqhJEmqI4wzSdoba2bBe31h0YfMKWvFxoltaba+nAWHNaTDw49zeedzol4oSZLqGONMkr6Oravhw8dhxp9YX9WYGXMPpdXsrcSbBKzo+wMuvOo+YkEs6pWSJKkOMs4k6auo3AETX4CPf0pFRSXjNnQmd8wamgZbmfutzpx13xAaNWgW9UpJklSHGWeS9N+EIcx9E959iHDzcj4Nu1D23jpabl5DUedcjug3kMsO/UbUKyVJUj1gnEnSl1n9GYzsA8vGsSqlI3MKD6Ng/io25iVT+fRPuOTiWwiCIOqVkiSpnjDOJOk/la6Hjx6Hab9jR1Jjxq05kbyxK8lJgaJrTuXcO39KZnqDqFdKkqR6xjiTpH+qqoDJv4AxzxBWlDElPI3wjSUUlK5k7sn5HN9vEF3aHh31SkmSVE8ZZ5IEMH/U57cwblrE0syuLJ64mZZLF7K8VRrh0/dx2VlXR71QkiTVc8aZpMS2aQmMvB/mj2Rb1kF8srIr+eNWkJUZsOjH53HuLU+RmpIW9UpJkpQAjDNJiamyHMYNhnGDiAcpTKg6i/Q/z6PljhXMO6sDpzz8PF3zDop6pSRJSiDGmaTEU/QOvHMflCxjftbprHqnmLxVc1lyUBYtH+rL5d+4JOqFkiQpARlnkhLHpsXwzv2wYBSbMw5h6qLjaDVlISkNYyy751uce10/kpP8Y1GSJEXDf4VIqv8qy2HcIBg3mGpSGFd6Gg2HLSavejtzLzyCbg8OpUnjllGvlCRJCc44k1S/7XEL46yUUyh5p5jmGxax4IgcDnn4CS47pnvUCyVJkgDjTFJ9VbIc3r4X5r/DurSOfDrrKFrNWkJ1kyRW9b+Oi77bmyAIol4pSZK0m3EmqX6proJJP4OPnmBnFYzb2JVmH66gaayUou+cwNn3DSE7q3HUKyVJkv4/xpmk+mPlNHjrDsLVMymsOo4dH2wgv2QF845vxtH9nqVnxxOjXihJkvSljDNJdd+OrfDhYzD5V6yMt2De9I4ULFzDphYpVD/bi54X3OgtjJIkqdYzziTVXWEIc9+Ed+6jfNNaxq0+ghafbCIntZT515/Bubc/S0Z6dtQrJUmSvhLjTFLdVLIc3u5NvGgkk0sOIhifT37pJopObc2J/QbTpfXhUS+UJEn6WowzSXVLvBom/Rw+fJzFJcksndGelivKWdYmndiz93JZtyujXihJkrRXjDNJdce6ufDGT9i6uJBPlrSjYEY5mVkVLL71As758QBSk9OiXihJkrTXjDNJtV9VBYwfTPyjZxi3sglZU1vQsqKconMO5tS+QzmpebuoF0qSJO0z40xS7VY8Hd68jbnzFrOuMI/m6+IsPiSLVn37c9mJF0S9TpIkqcbsc5wFQZAETAWKwzC8cN8nSRJQWQ6jn2TThy8xdXYzWs/PIakRrLj/Ss77wYMkxZKiXihJklSjauKVszuAuUDDGvhYkgTLJlA57BbGTd9Czoxm5MVh3iVH063P8zTOyYt6nSRJ0n6xT3EWBEEr4AJgAHBXjSySlLh2lsL7/fjs7VfZOq0RLUrSWXBkYzr1f4pLjzw96nWSJEn71b6+cjYYuBdoUANbJCWypeNY8/ubmTmhglZLG1GVm8Tax2/gosvvIAiCqNdJkiTtd0EYhnv3G4PgQuCbYRjeEgRBN+CeL/qesyAIbgRuBMjLyzv+tdde24e5+0dpaSnZ2dlRz0hInn10asvZx6p3UlD0MkvGTqH5zDSqkuCzsw6l5fk3kJaaFfW8/aK2nH0i8uyj49lHy/OPjmcfndp69t27d58WhmGXL3rfvsTZk8DVQBWQzuffczYsDMPvf9nv6dKlSzh16tS9+nz70+jRo+nWrVvUMxKSZx+d2nD24YopTH3uh1SNqyCnNGDeCc3p3H8w7Q7qHOmu/a02nH2i8uyj49lHy/OPjmcfndp69kEQfGmc7fVtjWEY9gH67PoE3fj8lbMvDTNJ2q1qJ8v+cAcL/vgRBcUxilukED52J5eef33UyyRJkiLjzzmTtN8NLyxm4KgiVpWUc3bGQi6b+ytaf1pNo7QYi647g3PuHERaakbUMyVJkiJVI3EWhuFoYHRNfCxJ9cvwwmL6DJtJZUU5vdcPpUvharLK4bMuTTjj6Zc5oaBT1BMlSZJqBV85k7RfPTtyLpeUvMbZn02iYC0szo8x4IwL2NDiPK40zCRJknYzziTtH2HIlkmv03vSo3ScF1KSBa+e3p7fN74JgmSCkvKoF0qSJNUqxpmkGle9ZALjnu1F9rhtdKiED49pxEttb2Z7rNnua/Jz/B4zSZKkPRlnkmrO2jnM/nUvNry9nOYbAxYdlMaGa+7mhZnNKK+s3n1ZRkoSvXt4S6MkSdKejDNJ+27TEtYPe5Dpf59Km0VJxBrFKL7/O5z/g74kxZLIOeRfT2vMz8mgd49O9OxcEPVqSZKkWsU4k7T3tq2l4v0n+HjYW+QWptKCJIouOZIzH3yJhg3/dQtjz84FxpgkSdL/YJxJ+vrKSwjHDabwzVcom5xO/pZU5h+VwxGPPEfPw0+Oep0kSVKdZJxJ+uoqy2HSzyl+ZwhzpqbSalkGW5omsf7Jm7m4560EQRD1QkmSpDrLOJP0v8XjMPP/2DHyET4u3EnezAxyk2DB1adw9l0/JTOjYdQLJUmS6jzjTNJ/t+RjwlEPMLlwKfFp2bQqTaWoa0uO6z+Y49ofHfU6SZKkesM4k/TF1s+H9x5m6bQPWfhpLgUrs1mZn0bsibvoee4Pol4nSZJU7xhnkv5d6XoY/STbJv6eCfOb0GpWLg3TAxbfdC7n3PoUqanpUS+UJEmql4wzSZ+rLIeJLxIfO5jxS5NJn9aMVjtgfrcOnPzw85zY8qCoF0qSJNVrxpmU6MIQZg+Ddx9m3vKNrJ7RlBZrqlnaPpOWffty6ck9o14oSZKUEIwzKZGt/gxG3s+m+Z8weX4+bWc3Jq1BnOV3Xc45P+xPcpJ/REiSJB0o/stLSkTbN9Kx6CWqP3qPsSub02hKHgVV1RR983DOeGgojZvkR71QkiQp4RhnUiKproQpv4HRT7ChuJqFhQW02FDNokMbcVD/AfQ89qyoF0qSJCUs40xKFIs+gpH3s3b5QmbMK6BNUSU0htUPX8s3r+hNLBaLeqEkSVJCM86k+q5kOYzsw85Zb/HxylY0ndqcPCqZcs4hXPbkK2RnN4l6oSRJkjDOpPqrqgI+eZFw9DNMW5vOzimtKCiJM79zM47qP5A2q8sNM0mSpFrEOJPqo6XjYMTdrFi6mLmzW9J6cSVbmyex6enbuPjimwiCgOWrR0e9UpIkSXswzqT6pHQ9vNeXsmmv8/GSluQXNiU3pZJFPzids+96jvT07KgXSpIk6UsYZ1J9EI/DtN8Sf/8RJi2NEUxrSZvSkKKTW3FCv8Ec3/aIqBdKkiTpfzDOpLpu9afw1p0smjebpZ82JX9lnBWt0kh+ujc9z7oq6nWSJEn6iowzqa6q2A4fPcHWsT9nwvzmtJ7VhAYZIUtvuYCzbhlAanJa1AslSZL0NRhnUl206CPib9zGx7O2kVnYnFY7YMGZB3PKw89zYl77qNdJkiRpLxhnUl1StgnefYi5o4exZnpjWqzLZGmHLAoe7k/Pky6Mep0kSZL2gXEm1QVhCLP/zsZh9zJlaoy28xuT2jDGynu/zbnX9iUplhT1QkmSJO0j40yq7bYUU/VGL8a+P5XGM9LJr4b5Fx3FGQ8OJScnL+p1kiRJqiHGmVRbhSFMe4UZf3yUrZPSaLk5nUWH59Cx/1NccvQZUa+TJElSDTPOpNpoy0pW/+4GPhu5mDaL0wmbJLHukeu44Dt3EQRB1OskSZK0HxhnUm0Shuyc+DJjXnqavBlJNA+SWPCdrpx132CysnKiXidJkqT9yDiTaolw62qmPHMlVaPW0HpbEgs653LMo0PofMjxUU+TJEnSAWCcSREZXljMwFFFrCop43reouv00bReHrC6eTJb+t3BxRf+KOqJkiRJOoCMMykCwwuL6TNsJo3LV/Dk8hc4YlYFO1MCJl9yLFc8+hvS0jKjnihJkqQDzDiTIvD8O4XcvGowJxUWk10GUw7P4PmDbiSjcUeuMcwkSZISknEmHUjVlSx49QH6vPcWrVbD0hYBL512JmMyzgcgKCmPeKAkSZKiYpxJB0I8Tsm4XzNx6BDazIzTIBP+eOrB/Dn3h4RByu7L8nMyIhwpSZKkKBln0v4UhlTPe4exQx+gwYSdtN4JC0/Pp+yqJxn2wSbCyurdl2akJNG7R6cIx0qSJClKxpm0v6yYzKzf3sWGdzfQYkPAkg7ptO3/BJec+PktjBmN/vm0xnLyczLo3aMTPTsXRDxakiRJUTHOpJq2bi7rhj3A9BGzabsgiZSGMVbfdwU9rnmQpFjS7st6di4wxiRJkrSbcSbVlC3FVLz7KGP/8R65hWm0jCex4OKj6PbgSzRs1DTqdZIkSarljDNpX+0shQnPM/3vv2T7lHQKNqex6MjGHPbIQC4+4pSo10mSJKmOMM6kvRWvhk9fpfiNx5g5KaTtkgy25Sax8fGbuODynxAEQdQLJUmSVIcYZ9LeWDKW8rfuZ+zEtbT4NJW8GCy68mTO6j2IjMyGUa+TJElSHWScSV/HhoWE7z7EpLHjiE9rQJttqSw4oSWdHxlE5w7HRL1OkiRJdZhxJn0VO7bCmKdZ8sHLLChsROsVDVjVMpXkx+7k4vOujXqdJEmS6gHjTPpvwhA++wulb/Vl3Iw4rWbl0DgNlv7oXM667SlSUzOiXihJkqR6wjiTvsyamcRH3M24T4pIm5ZJ67IkFp7Rnm/0e54T8g+Oep0kSZLqGeNM+k/lJfDRE8x7/w8UT88hf00my9tmkv7Qg1xy2mVRr5MkSVI9ZZxJ/xSGMOPPbPrHw0wqjNFmbg5ZWQEre13KWTf0Jzk5JeqFkiRJqseMMwlgwwKq3ridMePmkjM9g1aVsPC8wzi97ws0zs2Pep0kSZISgHGmxFa1E8YN5tPhL7B5ajb5GzJY0rEhHfoP4JLjzo56nSRJkhKIcabEtWwia1/9CdPHbaPdwmySc2KsfegHnHdVb2KxWNTrJEmSlGCMMyWe8hIq3nmI0X9/m+Yz0sgPk1h42fF0v38I2Q1zo14nSZKkBGWcKbHMeYMpv7qXHRMDWpeksejoXI549DkuOrRr1MskSZKU4IwzJYbtG1n+ux8x583ZtF0ao6xpMiVP3sKFl/446mWSJEkSYJwpAZRNfZWxgx8hvzCgeVKMxd87lbPuGUR6RnbU0yRJkqTdjDPVG8MLixk4qohVJeXk52TQ54xmNH6rF7FRq2hbGrDghGZ0eexnHNfuiKinSpIkSf8f40z1wvDCYvoMm0l5ZTUAXde8SspDH9OkOKC4RRKpA3pzcY9rIl4pSZIkfTnjTPXCwFFFlFdWc1j1PK5d8juOnF1JWVrA/53eiT4vvk5KSlpIU6DYAAAYKklEQVTUEyVJkqT/yjhTvVCxeRUPbniJY6ZvJmsHTDk8k+fb38zm5HweNswkSZJUBxhnqtt2ljLnlTt56pOPyV8bsDg/xpNHXUxhxqkAFORkRDxQkiRJ+mqMM9VN1ZVseH8wk3/xMu3nQGZ2wG/OOIa/5lwFQRIAGSlJ9O7RKeKhkiRJ0ldjnKluCUMqZw5jzAv9aPxJNa2rYME57Tn90d9w6tI4E/d4WmPvHp3o2bkg6sWSJEnSV2Kcqe5Y/gnTf30nWz/YTMHGgCWHZHLIY89y8bHdAejZGGNMkiRJdZZxptpv/XxW/bU3M0bMpf3CJLbnJLG+77Wcd+XdxGKxqNdJkiRJNcI4U+1VXsKOdx9l9F/foEVhKgUksejyLnTv8zxZ2Y2jXidJkiTVKONMtU+8mnDa75j0p6eomhij7ZZUFh3blKMeHcQxHbtEvU6SJEnaL4wz1S5Lx7PstTuZ+9EW2i5LYm2zZLYMvJ0LL/pR1MskSZKk/co4U+2wZSWlb97L2BGTaP1ZKs1SYiy9pjtn3fksqemZUa+TJEmS9jvjTJEK4lXExw1h/J+eJ2VyKu23p7LwG6044ZHnOb7NYVHPkyRJkg4Y40zRWTGFnPfv5INPKmm1KpXigjTSBvbhojO/G/UySZIk6YAzznTglW+mZPh9jP/bR7SbnUyjjBgrb7mA7j8eQHJKatTrJEmSpEgYZzpwwpDqT//CmJ89TIPJMdrtSGZm15ac+8wrNGneJup1kiRJUqSMMx0Y29Yw88VrWPf2MvLXxVjePpPs/o/TojzDMJMkSZKAWNQDVM+FIWvfG8zbV3cj+ZXlZJTHWNP7Ss4ZMYXDup4f9TpJkiSp1vCVM9Wo4YXFDBxVxKqSco7NLuH6FYMpmFhK66qARed34vR+v6RhTvOoZ0qSJEm1jnGmGjO8sJg+w2ays7KS67e9zOlj59F8Eyw4OJ3DnxjChUefHvVESZIkqdYyzlRjBo6cxxllI/nm3FEcshjW5cCQ7ifxWdvvM8EwkyRJkv4r40w1onzRBG6beSeHfVZJPIB3ujTnZ/m3UhnLIijZEfU8SZIkqdbb6zgLgqA18HugBRAHfhmG4ZCaGqa6Ib55BROH3ABvL+eYrVB4cDpDO17H6tSDdl+Tn5MR4UJJkiSpbtiXV86qgLvDMJweBEEDYFoQBO+FYTinhrapNtu5jUWv3cuCP31E2+UBa5rFmNLrBzxZfAzlldW7L8tISaJ3j04RDpUkSZLqhr2OszAMVwOrd/3vbUEQzAUKAOOsPotXs23czxn78xdoMwOapQQsu/oUzuw9lNTUDBru8bTG/JwMevfoRM/OBVGvliRJkmq9GvmesyAI2gGdgUk18fFUO8UXfMjYF+8m/eMddNgOC09qQdcBv+T4gkN2X9Ozc4ExJkmSJO2FIAzDffsAQZANjAEGhGE47AvefyNwI0BeXt7xr7322j59vv2htLSU7OzsqGfUWhllxQSTXmTL2E20Wh2womUyJVdcQfNOp+zzx/bso+PZR8ezj45nHx3PPlqef3Q8++jU1rPv3r37tDAMu3zR+/YpzoIgSAHeAkaFYfjT/3V9ly5dwqlTp+7159tfRo8eTbdu3aKeUftUlLFpRD8m/vEN2s1OojQzoPS6i+h2ywCSkmrmQZ+efXQ8++h49tHx7KPj2UfL84+OZx+d2nr2QRB8aZzty9MaA+A3wNyvEmaqW6rmjGD08/eS80mctjuTWHL2wZzS72c0btYq6mmSJElSvbQvL3+cAlwNzAyCYMautz0QhuHb+z5LkSlZwYyf/5BNI5ZTsD5geYdM2j36FBd2OSfqZZIkSVK9ti9PaxwHBDW4RVGKx1nz9gCm/eqPdCiKkd4gxvr7vsc51/QhFotFvU6SJEmq92rmG4dUp1WsmsmHj15L8wlltK6OsfjCIzij7y/IbpQb9TRJkiQpYRhnCWj4rp9FtqZkO7eV/4rjJi+g7WZYcmgWhw8YygVHfCPqiZIkSVLCMc4SzPDCYvoMm8nR28dx79xhdFwSsrYxTL35cq7u9XjU8yRJkqSEZZwlmN+P+ID7Fg2m86flxGMw8oQ8Xmr5E5rvaMLVUY+TJEmSEphxliDiZSWMH3Q9d78xlybb4NNDUnmp4zUsT+kEwKqS8ogXSpIkSYnNOKvvqipY+NeHWfibN2i7AoqbBTzV9UzGZJ//b5fl52RENFCSJEkSGGf1Vxiy5ZPfMW7oQNoWxmmaCiu+35VtPfox+c0iqKzefWlGShK9e3SKcKwkSZIk46weql49i7HP/pDMj7bSrgyWnJTHSQN+TZeCgwFISk5l4KgiVpWUk5+TQe8enejZuSDi1ZIkSVJiM87qkx1bmf3bW1n1+hRarQlY2SqF7CGPcOFpl/7bZT07FxhjkiRJUi1jnNUHYcjGj3/BJ88Pod0saJAVsPonF3Hmj58gKcn/iyVJkqS6wH+513GVG5Yx+rEraTx6M20qYcmZ7Tntsd/QKLdl1NMkSZIkfQ3GWR0zvLB41/eLlXFzxR84cfJntNoAyzqkcdCjz3Jhl7OjnihJkiRpLxhndcjwwmL6DJtJx/Ip3FX0OocvqGZjQ5h67bl8795BxGKxqCdKkiRJ2kvGWR3y8ogx9F46iM4zSonF4cPjGjO04DZyyeNqw0ySJEmq04yzuqCynMkv/Ii7R0yj+WaY0z6ZXx56FUVpRwOwqqQ84oGSJEmS9pVxVpuFIcvfeY5ZL/6G9ougrHHA4O7fYFSjy/7tsvycjIgGSpIkSaopxlkttX3+R4x56k4KJu2kZQyWXnokZZcMYOyIxVBZvfu6jJQkevfoFOFSSZIkSTXBOKtl4uVbmPDc9wneXEj7rbD4mEZ0HvALOh98DABJaVm7ntZYTn5OBr17dPIHSkuSJEn1gHEWsX89Gr+cq4KRnDr9fdotgzXNYpQPvIMLLrrx367v2bnAGJMkSZLqIeMsQv98NH7+ziIGLHuZo2ZWUJEChRcdzbcH/I6U1PSoJ0qSJEk6QIyzCA19Zzq3rh3I8dPWk7Mdph6ayZCDbiIt92CuMswkSZKkhGKcRSEMmfuXvvR+/2+0XQXL82IMPvmbTMrqBkDgo/ElSZKkhGOcHWCbZr/PhAF3076wgibp8KeTO/KnZtcTBv/6v8JH40uSJEmJxzjbj/Z82MfBDaq4ZcMQCkavpv0OWHpyc7Zf/Sx/+2AjoY/GlyRJkhKecbaf/PNhH+WVVXyn7P84d+JkCtbCslbJNOrXnwtOuxyA1EbFPhpfkiRJknG2vwwcVUS7HbO5bvFvOWJ2NVuz4JXTjuDjDjcz4bSzd1/no/ElSZIkgXG2X8R3bueKRY9w4rR1ZJfDhCMbMKT9TyhNyiXYsjPqeZIkSZJqIeNsH+35fWX5ORn0bjaBjNf+xNkrYGmLGE+ddgmFGafsvt6HfUiSJEn6IsbZPvjX95VV06x6LTfNGMpBn+1gZyrMuOxYHot9j7KqYPf1PuxDkiRJ0pcxzvbBwFFFlFdW8b1tf+C8qZ/RdAtMOzSdVzv3YkS/a8j4z1fVfNiHJEmSpC9hnH1Ne97GeFD1Qp5Z+GuOmlvF2sbw03NO472sSwh2/QxpH/YhSZIk6asyzr6Gf97GuKOykpu2vMQZU5bScDt8fExDBrXtRXmsIeD3lUmSJEn6+oyz/2HPV8piQcBBlXO4qei3HL4gzqqmAYNPvoBJmd12X+/3lUmSJEnaG8bZf7HnAz8Iq7lh00ucNWUZGTvhg+OaMrhVL6pi6QAE4PeVSZIkSdprxtmXGF5YzN1/+ZTqMOSoik+5Ye6f6Lgkzoq8gGeOvZSpGSfvvrYgJ4Px958Z4VpJkiRJdZ1x9h+GFxbT/83ZlJRXkhru4K61Q/nG9LXE4vDuic0Z0qIX8Vjq7uu9jVGSJElSTTDO9vDP2xiTKtZyzbbfc/pnq8jfAPPbxPj5kVcyN7UzAElBQDwMvY1RkiRJUo0xzvYwcFQReRuX8Nz458ncCRsbwatndOT3OT+EIAn4/JWyJy87yiCTJEmSVKOMsz2sKinn7lV/J7ka3junBYMz7yS+K8rg81fMDDNJkiRJ+4NxtsvwwmKSCDlheTFz2sNvs35EnH+Fma+YSZIkSdqfjDNgwqpK/vDBTHLKN9OwLM7m/CQ202j3+xtnptDvoiMMM0mSJEn7jXEG/G1+JeWVIQdVrQCgNLkp8PltjM995xijTJIkSdJ+F4t6QG2wcUcIQJtdcbYupRkA8TA0zCRJkiQdEMYZkJseANC8ah0A65LzAMjPyYhskyRJkqTEYpwBl3dMISMliSYVJQAUJ7f2h0tLkiRJOqD8njPg5PwUDj/scLY9vZ3KJKhsdChPnueTGSVJkiQdOMbZLj07FzCKarZlwNg+l0Q9R5IkSVKC8bbGPQSlOynLCqKeIUmSJCkBGWd7SCmtYmeWRyJJkiTpwLNE9pBaHlKVlRL1DEmSJEkJyDj7p3icrDKgoY/PlyRJknTgGWe7vPXxNNIqYVllwClPfcjwwuKoJ0mSJElKIMYZMGFVJX/6x7sAlKZmUVxSTp9hMw00SZIkSQeMcQb8bX4ljXesBGBzaiMAyiurGTiqKMpZkiRJkhKIcQZs3BGSW7UBgA0pzXa/fVVJeVSTJEmSJCUY4wzITQ/IqdwCwJqkFrvfnp/jw0EkSZIkHRjGGXB5xxQaVGwHYEVyWwAyUpLo3aNTlLMkSZIkJZDkqAfUBifnp1ASVLEjBcpiTSjIyaB3j0707FwQ9TRJkiRJCcI42yW1vILtGbDkqYuiniJJkiQpAXlbI58/Sj++tYLtGf6MM0mSJEnRSPg4G15YzCuzKkgvj1OWEfgzziRJkiRFIuHjbOCoIirikFkOO9I/v8vTn3EmSZIk6UBL+DhbVVJOEFaRXQZl6an/9nZJkiRJOlASPs7yczJoEV9DSjVsT03/t7dLkiRJ0oGS8HHWu0cn2lavAGBbagPAn3EmSZIk6cBL+Efp9+xcwNp/rANgc2pjf8aZJEmSpEgk/CtnwwuLWbd2PQA7G+QbZpIkSZIikdBxNrywmD7DZpJevg2ABZVNfIy+JEmSpEgkdJwNHFVEeWU12Tt3ALA2paWP0ZckSZIUiYSOs38+Lj9r507K0qAkaPxvb5ckSZKkAyWh4+yfj8vP3FHJtgwIdz0fxcfoS5IkSTrQEjrOevfoREZKEpk7qtie+fnbfIy+JEmSpCgk9KP0//lUxswP46zPwcfoS5IkSYpMQscZfB5oU8pDVrdKYvz9Z0Y9R5IkSVKCSujbGgHCeJyMcqjKSvhOlSRJkhShhI+z6i1bSAohnp0S9RRJkiRJCcw427QJgLBBesRLJEmSJCWyhI+zqo2fx1lSg8yIl0iSJElKZAkfZ2XrVgKQ0rhhxEskSZIkJbKEj7Py1csASG2cE/ESSZIkSYks4eOMjRsB6NC8fcRDJEmSJCWyfYqzIAjOC4KgKAiChUEQ3F9Tow6k1B1JxDMzOPakW6OeIkmSJCmB7XWcBUGQBLwInA8cDlwZBMHhNTXsQKnaXEK8YSNIbxT1FEmSJEkJbF9eOTsRWBiG4eIwDCuA14BLambWgVO9aRPx7OyoZ0iSJElKcPsSZwXAij1+vXLX2+qUqk2biDdoEPUMSZIkSQkuCMNw735jEHwb6BGG4Q27fn01cGIYhrf9x3U3AjcC5OXlHf/aa6/t2+Ia1njQILYXtKLiO9+OekpCKi0tJdtXLiPh2UfHs4+OZx8dzz5ann90PPvo1Naz7969+7QwDLt80fuS9+HjrgRa7/HrVsCq/7woDMNfAr8E6NKlS9itW7d9+JT7QbdujB49mlq3K0F49tHx7KPj2UfHs4+OZx8tzz86nn106uLZ78ttjVOAQ4IgaB8EQSpwBfBmzcySJEmSpMSy16+chWFYFQTBT4BRQBLwchiGs2tsmSRJkiQlkH25rZEwDN8G3q6hLZIkSZKUsPbph1BLkiRJkmqGcSZJkiRJtYBxJkmSJEm1gHEmSZIkSbWAcSZJkiRJtYBxJkmSJEm1gHEmSZIkSbWAcSZJkiRJtYBxJkmSJEm1gHEmSZIkSbWAcSZJkiRJtYBxJkmSJEm1gHEmSZIkSbWAcSZJkiRJtYBxJkmSJEm1gHEmSZIkSbWAcSZJkiRJtYBxJkmSJEm1gHEmSZIkSbWAcSZJkiRJtUAQhuGB+2RBsB5YdsA+4VfXFNgQ9YgE5dlHx7OPjmcfHc8+Op59tDz/6Hj20amtZ982DMNmX/SOAxpntVUQBFPDMOwS9Y5E5NlHx7OPjmcfHc8+Op59tDz/6Hj20amLZ+9tjZIkSZJUCxhnkiRJklQLGGef+2XUAxKYZx8dzz46nn10PPvoePbR8vyj49lHp86dvd9zJkmSJEm1gK+cSZIkSVItkNBxFgTBeUEQFAVBsDAIgvuj3lOfBUHQOgiCj4IgmBsEwewgCO7Y9fb+QRAUB0EwY9d/vhn11voqCIKlQRDM3HXOU3e9rUkQBO8FQbBg1383jnpnfRMEQac9vr5nBEGwNQiCXn7t7x9BELwcBMG6IAhm7fG2L/w6Dz73/K6/Az4LguC46JbXfV9y9gODIJi363z/HgRBzq63twuCoHyPr/+fR7e87vuSs//SP2OCIOiz6+u+KAiCHtGsrh++5Oxf3+Pc/1879xdiRRnGcfz7oCVkRUElpkUWehFdaIQGogiVZIRbQbFL1PYHStAgupHqovDKJIOuuggFA/9klLREfyyCurKW3QJTo9SsNpeVEioRit1+Xcy7cNzOLLGenTln5ve5OTPvnoWHd5995n1n3nlPRMTXqd1530KTjC07uubXdlljRMwAvgPuAIaAfqBH0uFSA6uoiJgLzJU0GBGXAAPAPcADwBlJL5caYA1ExAngFkm/NrRtAU5L2pxuUFwuaWNZMVZdqju/AMuAR3Hut1xErATOAG9Iuim1Nc3zNFh9CriL7G/yqqRlZcXe6XL6fjXwqaTRiHgJIPX9dcB749+z85PT9y/SpMZExI3AbmApcDXwCbBI0lihQVdEs76f8POtwO+SNjnvW2uSseUjdHDNr/OTs6XAUUnHJf0N7AG6So6psiQNSxpMx38CR4B55UZlZDm/Ix3vICtqNn1uA45J+rHsQKpK0ufA6QnNeXneRTagkqQDwGXpYm9T0KzvJe2XNJpODwDzCw+sBnLyPk8XsEfSX5J+AI6SjYlsCibr+4gIspvQuwsNqiYmGVt2dM2v8+RsHvBzw/kQniwUIt05WgJ8kZo2pMfL272sbloJ2B8RAxHxRGqbI2kYsiIHXFVadPXQzbkXaed+MfLy3NeBYj0GfNBwviAivoqIzyJiRVlBVVyzGuO8L84KYETS9w1tzvtpMGFs2dE1v86Ts2jSVs81ngWKiIuBt4GnJf0BvAbcACwGhoGtJYZXdcsl3QysAdanpRhWkIi4EFgLvJWanPvl83WgIBHxPDAK7ExNw8C1kpYAzwC7IuLSsuKrqLwa47wvTg/n3pBz3k+DJmPL3K82aWu73K/z5GwIuKbhfD5wsqRYaiEiLiD759kp6R0ASSOSxiT9A7yOl1ZMG0kn0+cpYB9ZX4+MP9JPn6fKi7Dy1gCDkkbAuV+wvDz3daAAEdEL3A08qPSie1pS91s6HgCOAYvKi7J6JqkxzvsCRMRM4D7gzfE2533rNRtb0uE1v86Ts35gYUQsSHe0u4G+kmOqrLTuehtwRNIrDe2Na33vBb6Z+Lt2/iJidnpZloiYDawm6+s+oDd9rRd4t5wIa+GcO6jO/ULl5Xkf8HDawetWspf2h8sIsKoi4k5gI7BW0tmG9ivTBjlExPXAQuB4OVFW0yQ1pg/ojohZEbGArO+/LDq+Grgd+FbS0HiD87618saWdHjNn1l2AGVJO0dtAD4CZgDbJR0qOawqWw48BBwc31IWeA7oiYjFZI+VTwBPlhNe5c0B9mV1jJnALkkfRkQ/sDciHgd+Au4vMcbKioiLyHaGbczvLc791ouI3cAq4IqIGAJeADbTPM/fJ9u16yhwlmwHTZuinL5/FpgFfJzqzwFJ64CVwKaIGAXGgHWS/u+GFjZBTt+valZjJB2KiL3AYbKlpuu9U+PUNet7Sdv47zvG4LxvtbyxZUfX/NpupW9mZmZmZtZO6rys0czMzMzMrG14cmZmZmZmZtYGPDkzMzMzMzNrA56cmZmZmZmZtQFPzszMzMzMzNqAJ2dmZmZmZmZtwJMzMzMzMzOzNuDJmZmZmZmZWRv4F0gEZ8PeUZIGAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# place holder\n", + "exo_grid = UnstructuredGrid(np.array([[0.2, 0.5, 0.7]]))\n", + "N=20\n", + "\n", + "g = WarpGrid(\n", + " CartesianGrid([-1], [5], [N]),\n", + " ['exp'] )\n", + "f = lambda x: np.minimum(x[:,0],1+0.05*(1+(x[:,0]-1)*(1-np.exp(1-(x[:,0])))))[:,None]\n", + "nodes = g.nodes\n", + "vals = f(nodes).reshape((1,N,1))\n", + "\n", + "\n", + "dr = DecisionRule(exo_grid, g, 'linear')\n", + "dr.set_values(vals)\n", + "\n", + "dr2 = DecisionRule(exo_grid, g, 'cubic')\n", + "dr2.set_values(vals)\n", + "\n", + "fg = CartesianGrid([0], [200], [1000])\n", + "nn = fg.nodes\n", + "\n", + "tvals = f(nn)\n", + "xx = dr.eval_is(0, nn)\n", + "xx2 = dr2.eval_is(0, nn)\n", + "\n", + "plt.figure(figsize=(15,10))\n", + "plt.plot(nodes.ravel(), vals.ravel(), 'o', label='data')\n", + "plt.plot(nn.ravel(), xx.ravel(), label='approx (linear)')\n", + "plt.plot(nn.ravel(), xx2.ravel(), label='approx (cubic)')\n", + "plt.plot(nn.ravel(), tvals.ravel(), label='true')\n", + "# plt.xscale('log')\n", + "# plt.yscale('log')\n", + "# plt.ylim(0,10)\n", + "# plt.xlim(100,120)\n", + "plt.legend(loc= 'upper left')\n", + "plt.grid()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/experiments/test_hmodel.py b/experiments/test_hmodel.py index 79cc945..5912d9b 100644 --- a/experiments/test_hmodel.py +++ b/experiments/test_hmodel.py @@ -49,7 +49,7 @@ from matplotlib import pyplot as plt for i, (w, eq) in enumerate(eqss): - s =eq.dr.endo_grid.nodes() + s =eq.dr.endo_grid.nodes for j in range(3): plt.plot(s, w*eq.μ[j,:], label=f"{i}" ) plt.legend()