bayesian/bayesian_optimization/bayesopt_water_accumulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from sklearn import preprocessing\n",
    "import numpy as np\n",
    "import GPy\n",
    "from scipy import stats\n",
    "from scipy import optimize\n",
    "import matplotlib.pyplot as plt\n",
    "from tqdm import tqdm\n",
    "from functools import partial\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from anastruct.basic import converge\n",
    "from anastruct.material.profile import HEA, IPE\n",
    "from anastruct.fem.system import SystemElements, Vertex\n",
    "from anastruct.material.units import to_kNm2, to_kN\n",
    "import anastruct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "# anastruct version == anastruct==1.1.0a1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bayesian Optimization for water accumulation analysis\n",
    "\n",
    "See these posts:\n",
    "\n",
    "* [water accumulation](https://www.ritchievink.com/blog/2017/08/23/a-nonlinear-water-accumulation-analysis-in-python/)\n",
    "* [bayesian optimization](https://www.ritchievink.com/blog/2019/08/25/algorithm-breakdown-bayesian-optimization/)\n",
    "\n",
    "For optimization purposes we make the `structure` function take the following inputs:\n",
    "\n",
    "* mp_span_1\n",
    "* mp_span_2\n",
    "* mp_span_3\n",
    "* k1\n",
    "* k2\n",
    "\n",
    "These parameters will be optimized within the following boundaries:\n",
    "\n",
    "```python\n",
    "{'mp_span_1': {'range': [70, 91], 'type': 'float'},\n",
    " 'mp_span_2': {'range': [240, 312], 'type': 'float'},\n",
    " 'mp_span_3': {'range': [25, 32], 'type': 'float'},\n",
    " 'k1': {'range': [2000, 2600], 'type': 'float'},\n",
    " 'k2': {'range': [3000, 3900], 'type': 'float'}}\n",
    "```\n",
    "\n",
    "The stiffness of the rotation is assumed to be linearly increasing with `mp`. See the `structure` function;\n",
    "\n",
    "```python\n",
    "    first = dict(spring={1: 9e3 / 70 * mp_span_1}, mp={1: mp_span_1})\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# constants\n",
    "E = 2.1e5 / 1.3\n",
    "b = 5 * 1.09\n",
    "q_water = 10\n",
    "\n",
    "# axes height levels\n",
    "h_1 = 0\n",
    "h_2 = 0.258\n",
    "h_3 = 0.046\n",
    "h_4 = 0.274\n",
    "h_5 = 0.032\n",
    "h_6 = 0.15\n",
    "\n",
    "# beam spans\n",
    "span_1 = span_2 = 21.9\n",
    "span_3 = 8.9\n",
    "\n",
    "# Vertices at the axes\n",
    "p1 = Vertex(0, h_1)\n",
    "p2 = Vertex(span_1 * 0.5, h_2)\n",
    "p3 = Vertex(span_1, h_3)\n",
    "p4 = Vertex(span_1 + span_2 * 0.5, h_4)\n",
    "p5 = Vertex(span_1 + span_2, h_5)\n",
    "p6 = Vertex(span_1 + span_2 + span_3, h_6)\n",
    "\n",
    "\n",
    "def structure(\n",
    "    mp_span_1=70,\n",
    "    mp_span_2=240,\n",
    "    mp_span_3=25,\n",
    "    k1=2e3,\n",
    "    k2=3e3,\n",
    "):\n",
    "    \"\"\"\n",
    "    Build the structure from left to right, starting at axis 1.\n",
    "\n",
    "    variables:\n",
    "    EA = Young's modulus * Area\n",
    "    EI = Young's modulus * moment of Inertia\n",
    "    g = Weight [kN/ m]\n",
    "    elements = reference of the element id's that were created\n",
    "    dl = c.t.c distance different nodes.\n",
    "    \"\"\"\n",
    "\n",
    "    dl = 0.2\n",
    "\n",
    "    ## SPAN 1 AND 2\n",
    "\n",
    "    # The elements between axis 1 and 3 are an IPE 450 member.\n",
    "    EA = to_kN(E * IPE[450][\"A\"])  # Y\n",
    "    EI = to_kNm2(E * IPE[450][\"Iy\"])\n",
    "    g = IPE[450][\"G\"] / 100\n",
    "\n",
    "    # New system.\n",
    "    ss = SystemElements(mesh=3)\n",
    "\n",
    "    # span 1\n",
    "    first = dict(spring={1: 9e3 / 70 * mp_span_1}, mp={1: mp_span_1})\n",
    "\n",
    "    elements = ss.add_multiple_elements(\n",
    "        location=[p1, p2], dl=dl, first=first, EA=EA, EI=EI, g=g\n",
    "    )\n",
    "    elements += ss.add_multiple_elements(location=p3, dl=dl, EA=EA, EI=EI, g=g)\n",
    "\n",
    "    # span 2\n",
    "    first = dict(spring={1: 40e3 / 240 * mp_span_2}, mp={1: mp_span_2})\n",
    "    elements += ss.add_multiple_elements(\n",
    "        location=p4, dl=dl, first=first, EA=EA, EI=EI, g=g\n",
    "    )\n",
    "    elements += ss.add_multiple_elements(location=p5, dl=dl, EA=EA, EI=EI, g=g)\n",
    "\n",
    "    ## SPAN 3\n",
    "\n",
    "    # span 3\n",
    "    # different IPE\n",
    "    g = IPE[240][\"G\"] / 100\n",
    "    EA = to_kN(E * IPE[240][\"A\"])\n",
    "    EI = to_kNm2(E * IPE[240][\"Iy\"])\n",
    "    first = dict(spring={1: 15e3 / 25 * mp_span_3}, mp={1: mp_span_3})\n",
    "\n",
    "    elements += ss.add_multiple_elements(\n",
    "        location=p6, first=first, dl=dl, EA=EA, EI=EI, g=g\n",
    "    )\n",
    "\n",
    "    # Add a dead load of -2 kN/m to all elements.\n",
    "    ss.q_load(-2, elements, direction=\"y\")\n",
    "\n",
    "    ## COLUMNS\n",
    "\n",
    "    # column height\n",
    "    h = 7.2\n",
    "\n",
    "    # left column\n",
    "    EA = to_kN(E * IPE[220][\"A\"])\n",
    "    EI = to_kNm2(E * HEA[220][\"Iy\"])\n",
    "    left = ss.add_element([[0, 0], [0, -h]], EA=EA, EI=EI)\n",
    "\n",
    "    # right column\n",
    "    EA = to_kN(E * IPE[180][\"A\"])\n",
    "    EI = to_kNm2(E * HEA[180][\"Iy\"])\n",
    "    right = ss.add_element([p6, Vertex(p6.x, -h)], EA=EA, EI=EI)\n",
    "\n",
    "    ## SUPPORTS\n",
    "\n",
    "    # node ids for the support\n",
    "    id_left = max(ss.element_map[left].node_map.keys())\n",
    "    id_top_right = min(ss.element_map[right].node_map.keys())\n",
    "    id_btm_right = max(ss.element_map[right].node_map.keys())\n",
    "\n",
    "    # Add supports. The location of the supports is defined with the nodes id.\n",
    "    ss.add_support_hinged((id_left, id_btm_right))\n",
    "\n",
    "    # Retrieve the node ids at axis 2 and 3\n",
    "    id_p3 = ss.find_node_id(p3)\n",
    "    id_p5 = ss.find_node_id(p5)\n",
    "\n",
    "    ss.add_support_roll(id_top_right, direction=1)\n",
    "\n",
    "    # Add translational spring supports at axes 2 and 3\n",
    "    ss.add_support_spring(id_p3, translation=2, k=k1, roll=True)\n",
    "    ss.add_support_spring(id_p5, translation=2, k=k2, roll=True)\n",
    "    return ss\n",
    "\n",
    "\n",
    "def water_load(ss, water_height, deflection=None):\n",
    "    \"\"\"\n",
    "    :param ss: (SystemElements) object.\n",
    "    :param water_height: (flt) Water level.\n",
    "    :param deflection: (array) Computed deflection.\n",
    "    :return (flt) The cubic meters of water on the structure\n",
    "    \"\"\"\n",
    "\n",
    "    # The horizontal distance between the nodes.\n",
    "    dl = np.diff(ss.nodes_range(\"x\"))\n",
    "\n",
    "    if deflection is None:\n",
    "        deflection = np.zeros(len(ss.node_map))\n",
    "\n",
    "    # Height of the nodes\n",
    "    y = np.array(ss.nodes_range(\"y\"))\n",
    "\n",
    "    # An array with point loads.\n",
    "    # cubic meters * weight water\n",
    "    force_water = (water_height - y[:-3] - deflection[:-3]) * q_water * b * dl[:-2]\n",
    "\n",
    "    cubics = 0\n",
    "    n = force_water.shape[0]\n",
    "    for k in ss.node_map:\n",
    "        if k > n:\n",
    "            break\n",
    "        point_load = force_water[k - 1]\n",
    "\n",
    "        if point_load > 0:\n",
    "            ss.point_load(k, Fx=0, Fy=-point_load)\n",
    "            cubics += point_load / q_water\n",
    "\n",
    "    return cubics\n",
    "\n",
    "\n",
    "def det_water_height(c, deflection=None, hyperparams={}):\n",
    "    \"\"\"\n",
    "    :param c: (flt) Cubic meters.\n",
    "    :param deflection: (array) Node deflection values.\n",
    "    :return (SystemElement, flt) The structure and the redistributed water level is returned.\n",
    "    \"\"\"\n",
    "    wh = 0.1\n",
    "\n",
    "    while True:\n",
    "        ss = structure(**hyperparams)\n",
    "        cubics = water_load(ss, wh, deflection)\n",
    "\n",
    "        factor = converge(cubics, c)\n",
    "        if 0.9999 <= factor <= 1.0001:\n",
    "            return ss, wh\n",
    "\n",
    "        wh *= factor\n",
    "\n",
    "\n",
    "def determine_max_water_cubic(hyper_params, return_ss=False):\n",
    "    cubics = [0]\n",
    "    water_heights = [0]\n",
    "\n",
    "    deflection = None\n",
    "    max_water_level = 0\n",
    "\n",
    "    # Iterate from 8 m3 to 15 m3 of water.\n",
    "    for cubic in reversed(\n",
    "        (14 - np.logspace(0, 1, 35)[:-10])\n",
    "    ):  # This loop computes the results per m3 of storaged water.\n",
    "        wh = 0.05\n",
    "        lastwh = 0.2\n",
    "\n",
    "        print(f\"Starting analysis of {cubic} m3\")\n",
    "\n",
    "        c = 1\n",
    "        for _ in range(100):  # This loop redistributes the water until the water level converges.\n",
    "\n",
    "            # redistribute the water\n",
    "            ss, wh = det_water_height(cubic, deflection, hyper_params)\n",
    "\n",
    "            # Do a non linear calculation!!\n",
    "            ss.solve(max_iter=100, verbosity=1)\n",
    "            deflection = ss.get_node_result_range(\"uy\")\n",
    "\n",
    "            # Some breaking conditions\n",
    "            if min(deflection) < -1:\n",
    "                print(min(deflection), \"Breaking due to exceeding max deflection\")\n",
    "                break\n",
    "            if 0.9999 < lastwh / wh < 1.001:\n",
    "                print(f\"Convergence in {c} iterations.\")\n",
    "                cubics.append(cubic)\n",
    "                water_heights.append(wh)\n",
    "                break\n",
    "\n",
    "            lastwh = wh\n",
    "            c += 1\n",
    "\n",
    "        if wh > max_water_level:\n",
    "            max_water_level = wh\n",
    "        else:\n",
    "            print(\"Breaking. Water level isn't rising.\")\n",
    "            break\n",
    "    \n",
    "    if return_ss:\n",
    "        return ss, cubics, water_heights\n",
    "    \n",
    "    return max_water_level\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def select_hyperparams_random(param_ranges):\n",
    "    \"\"\"\n",
    "    Select hyperparameters at random.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    param_ranges : dict\n",
    "        Named parameter ranges.\n",
    "\n",
    "        Example:\n",
    "        \n",
    "        {\n",
    "            'foo': {\n",
    "                'range': [1, 10],\n",
    "                'type': 'float'\n",
    "            }\n",
    "            'bar': {\n",
    "                'range': [10, 1000],\n",
    "                'type': 'int'\n",
    "            }\n",
    "        }\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    selection : dict\n",
    "        Randomly selected hyperparameters within given boundaries.\n",
    "        \n",
    "        Example:\n",
    "        {'foo': 4.213, 'bar': 935}\n",
    "\n",
    "    \"\"\"\n",
    "    selection = {}\n",
    "    for k in param_ranges:\n",
    "        val = np.random.choice(\n",
    "            np.linspace(*param_ranges[k]['range'], num=100)\n",
    "        )\n",
    "\n",
    "        dtype = param_ranges[k]['type']\n",
    "        if dtype is 'int':\n",
    "            val = int(val)\n",
    "        selection[k] = val\n",
    "    return selection\n",
    "\n",
    "\n",
    "def expected_improvement(f, y_current, x_proposed):\n",
    "    \"\"\"\n",
    "    Return E(max(f_proposed - f_current), 0)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "\n",
    "    f : GP predict function\n",
    "    y_current : float\n",
    "        Current best evaluation f(x+)\n",
    "    x_proposed : np.array\n",
    "        Proposal parameters. Shape: (1, 1)\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    expected_improvement : float\n",
    "        E(max(f_proposed - f_current), 0)\n",
    "    \"\"\"\n",
    "    mu, var = f(x_proposed)\n",
    "    std = var ** 0.5\n",
    "    delta = mu - y_current\n",
    "\n",
    "    # x / inf = 0\n",
    "    std[std == 0] = np.inf\n",
    "    z = delta / std\n",
    "    unit_norm = stats.norm()\n",
    "    return delta * unit_norm.cdf(z) + std * unit_norm.pdf(z)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BayesOpt:\n",
    "    def __init__(\n",
    "        self, param_ranges, f, random_trials=5, optimization_trials=20, kernel=None\n",
    "    ):\n",
    "        \"\"\"\n",
    "        Parameters\n",
    "        ----------\n",
    "\n",
    "        param_ranges : dict\n",
    "        f : function\n",
    "            black box function to evaluate\n",
    "        random_trials : int\n",
    "            Number of random trials to run before optimization starts\n",
    "        optimization_trials : int\n",
    "            Number of optimization trials to run.\n",
    "            Together with the random_trials this is the total budget\n",
    "        kernel: GPy.kern.src.kern.Kern\n",
    "            GPy kernel for the Gaussian Process.\n",
    "            If None given, RBF kernel is used\n",
    "        \"\"\"\n",
    "        self.param_ranges = param_ranges\n",
    "        self.f = f\n",
    "        self.random_trials = random_trials\n",
    "        self.optimization_trials = optimization_trials\n",
    "        self.n_trials = random_trials + optimization_trials\n",
    "        self.x = np.zeros((self.n_trials, len(param_ranges)))\n",
    "        self.y = np.zeros((self.n_trials, 1))\n",
    "\n",
    "        if kernel is None:\n",
    "            self.kernel = GPy.kern.RBF(\n",
    "                input_dim=self.x.shape[1], variance=1, lengthscale=1\n",
    "            )\n",
    "        else:\n",
    "            self.kernel = kernel\n",
    "        self.gp = None\n",
    "        self.bounds = np.array([pr[\"range\"] for pr in param_ranges.values()])\n",
    "\n",
    "    @property\n",
    "    def best_params(self):\n",
    "        \"\"\"\n",
    "        Select best parameters.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        best_parameters : dict\n",
    "        \"\"\"\n",
    "        return self._prepare_kwargs(self.x[self.y.argmax()])\n",
    "\n",
    "    def fit(self):\n",
    "        self._random_search()\n",
    "        self._bayesian_search()\n",
    "\n",
    "    def _random_search(self):\n",
    "        \"\"\"\n",
    "        Run the random trials budget\n",
    "        \"\"\"\n",
    "        print(f\"Starting {self.random_trials} random trials...\")\n",
    "        for i in tqdm(range(self.random_trials)):\n",
    "            hp = select_hyperparams_random(self.param_ranges)\n",
    "            self.x[i] = np.array(list(hp.values()))\n",
    "            self.y[i] = self.f(hp)\n",
    "            print(f'HP: {hp}\\t y:{self.y[i]}')\n",
    "\n",
    "    def _bayesian_search(self):\n",
    "        \"\"\"\n",
    "        Run the Bayesian Optimization budget\n",
    "        \"\"\"\n",
    "        print(f\"Starting {self.optimization_trials} optimization trials...\")\n",
    "        for i in tqdm(\n",
    "            range(self.random_trials, self.random_trials + self.optimization_trials)\n",
    "        ):\n",
    "            self.x[i], self.y[i] = self._single_iter()\n",
    "\n",
    "    def _single_iter(self, x=None):\n",
    "        \"\"\"\n",
    "        Fit a GP and retrieve and evaluate a new\n",
    "        parameter proposal.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        out : tuple[np.array[flt], np.array[flt]]\n",
    "            (x, f(x))\n",
    "\n",
    "        \"\"\"\n",
    "        self._fit_gp()\n",
    "        if x is None:\n",
    "            x = self._new_proposal()\n",
    "        y = self.f(self._prepare_kwargs(x))\n",
    "        print(f'HP: {self._prepare_kwargs(x)}\\t y:{y}')\n",
    "        return x, y\n",
    "\n",
    "    def _fit_gp(self, noise_var=0):\n",
    "        \"\"\"\n",
    "        Fit a GP on the currently observed data points.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        noise_var : flt\n",
    "            GPY argmument noise_var\n",
    "        \"\"\"\n",
    "        mask = self.x.sum(axis=1) != 0\n",
    "        self.gp = GPy.models.GPRegression(\n",
    "            self.x[mask],\n",
    "            self.y[mask],\n",
    "            normalizer=True,\n",
    "            kernel=self.kernel,\n",
    "            noise_var=noise_var,\n",
    "        )\n",
    "        self.gp.optimize()\n",
    "\n",
    "    def _new_proposal(self, n=25):\n",
    "        \"\"\"\n",
    "        Get a new parameter proposal by maximizing\n",
    "        the acquisition function.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        n : int\n",
    "            Number of retries.\n",
    "            Each new retry the optimization is\n",
    "            started in another parameter location.\n",
    "            This improves the chance of finding a global optimum.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        proposal : dict\n",
    "            Example:\n",
    "           {'foo': 4.213, 'bar': 935}\n",
    "        \"\"\"\n",
    "\n",
    "        def f(x):\n",
    "            return -expected_improvement(\n",
    "                f=self.gp.predict, y_current=self.y.max(), x_proposed=x[None, :]\n",
    "            )\n",
    "\n",
    "        x0 = np.random.uniform(\n",
    "            low=self.bounds[:, 0], high=self.bounds[:, 1], size=(n, self.x.shape[1])\n",
    "        )\n",
    "        proposal = None\n",
    "        best_ei = np.inf\n",
    "        for x0_ in x0:\n",
    "            res = optimize.minimize(f, x0_, bounds=self.bounds)\n",
    "            if res.success and res.fun < best_ei:\n",
    "                best_ei = res.fun\n",
    "                proposal = res.x\n",
    "            if np.isnan(res.fun):\n",
    "                raise ValueError(\"NaN within bounds\")\n",
    "        return proposal\n",
    "\n",
    "    def _prepare_kwargs(self, x):\n",
    "        \"\"\"\n",
    "        Create a dictionary with named parameters\n",
    "        and the proper python types.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        x : np.array\n",
    "            Example:\n",
    "            [4.213, 935.03]\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        hyperparameters : dict\n",
    "\n",
    "            Example:\n",
    "            {'foo': 4.213, 'bar': 935}\n",
    "        \"\"\"\n",
    "        # create hyper parameter dict\n",
    "        hp = dict(zip(self.param_ranges.keys(), x))\n",
    "        # cast values\n",
    "        for k in self.param_ranges:\n",
    "            if self.param_ranges[k][\"type\"] == \"int\":\n",
    "                hp[k] = int(hp[k])\n",
    "            elif self.param_ranges[k][\"type\"] == \"float\":\n",
    "                hp[k] = float(hp[k])\n",
    "            else:\n",
    "                raise ValueError(\"Parameter type not known\")\n",
    "        return hp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "maximum_factor = 1.3\n",
    "\n",
    "param_ranges = {\n",
    "    'mp_span_1': {\n",
    "        'range': [70, int(70 * maximum_factor)],\n",
    "        'type': 'float'\n",
    "    },\n",
    "    'mp_span_2': {\n",
    "        'range': [240, int(240 * maximum_factor)],\n",
    "        'type': 'float'\n",
    "    },\n",
    "    'mp_span_3': {\n",
    "        'range': [25, int(25 * maximum_factor)],\n",
    "        'type': 'float'\n",
    "    },\n",
    "    'k1': {\n",
    "        'range': [2000, int(2000 * maximum_factor)],\n",
    "        'type': 'float'\n",
    "    },\n",
    "    'k2': {\n",
    "        'range': [3000, int(3000 * maximum_factor)],\n",
    "        'type': 'float'\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bo = BayesOpt(param_ranges, \n",
    "              determine_max_water_cubic, \n",
    "              random_trials=5, \n",
    "              optimization_trials=15, \n",
    "              kernel=GPy.kern.Matern52(input_dim=3) + GPy.kern.Matern32(input_dim=3))\n",
    "bo.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'mp_span_1': 79.73448251348799,\n",
       " 'mp_span_2': 265.11052593961267,\n",
       " 'mp_span_3': 29.237956757987163,\n",
       " 'k1': 2570.259044964325,\n",
       " 'k2': 3443.585144096509}"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bo.best_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'iterations')"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAEWCAYAAACufwpNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhcZZn38e8v6SSQRbLjAIEEWRRFEcL2jmhGhBcchiiDAiKiIqgjjuL4Oqgjg4zOiIrbyKggQUQ2RXEyDhhxADcGSICEJWwRg0kUCF2dkFSTVDp9v3+cp0JRVHdXkjpV1d2/z3XV1afOeld1dd39PM8591FEYGZm1ggjWh2AmZkNHU4qZmbWME4qZmbWME4qZmbWME4qZmbWME4qZmbWME4q1lSS1kvas9VxNJKkkLRXC477QUlPpfd0SrOPXyOeT0n6botjuEnS6a2MYbiTr1OxMknLgV2AXSLimYr59wIHALMiYnlromtfkgLYOyKWNfGYo4BngcMiYkmr4uiLpJnAH4BREdGT0zHOB/aKiHfmsX/bNm6pWLU/AKeUn0jaHxjbunCsDzsDOwAPNmJnkjoasZ9Gabd4rH5OKlbtSuBdFc9PB75fuYKkv5Z0r6RnJa1I/zGWl50k6Q+SXpKeHyvpSUnT0vMtXUWSvifpP1KXxXpJv5P0Uklfk9Ql6WFJr63Y9wu6mdL2n0vTcyStlPQJSU9L+rOkt0h6s6RHJRUkfarWC5Z0aIpxZMW8t0q6L00fIul/Ja1J+/2mpNF97Os2Se+reP5uSb+teP5ySTeneB6R9Pa+fhGS3iPpIUnrJD0u6f1p/j7AI2m1NZJukfTr9HxJei9PSuseJ2lxiv12Sa+u2P9ySf+YXmex1hd5es8/IOmxtI+LJamPeM+X9IP0tBzPmhTP4Wmd96bX1CVpgaQ9qo71IUmPAY+leV9Pn7FnJd0t6Yg0/xjgU8BJaf/l1tqW91/SCEn/JOmJ9Jn4vqSd0rKZ6XinS/qjpGckfboilkMkLUrHfUrSV/r6PVmViPDDDyICYDnwJrIvrFcAI4GVwB5AADPTenOA/cn+KXk18BTwlor9XAV8D5gC/Ak4rmJZkHVZkNZ5BjiI7L/uW8haSu9Kx/4ccGutbSu2/1xFTD3AecAo4ExgNXA1MAF4JfAcWRderdf+e+Coiuc/As5N0wcBhwEdwEzgIeCjfbym24D3VSx7N/DbND0OWAG8J+3rten179dHTH8NvAwQ8AagGzgwLZuZjtvRz/vzWuBp4ND0fp6efsdjKn7fi4EZwI59xBDAz4CJwO7pPT2mj3XPB37QT3xzgWVkn60O4J+A26uOdTMwuRwP8E6yz1EH8A/Ak8AO1cer2MeW9x94bzrensB44CfAlVXxXQrsCLwG2Ai8Ii3/X+C0ND2erJux5X+jg+HhlorVUm6tHEX2BbqqcmFE3BYR90dEb0TcB1xD9qVX9iHgjWR/4P8VET/r51g3RMTdEbEBuAHYEBHfj4jNwHVkX4z12gR8PiI2AdcCU4GvR8S6iHgQWEr25VHLNaRuP0kTgDeneaT47oiInsjGlL5T9XrrdRywPCIuT/u6F/gx8LZaK0fEf0fE7yPzK+AXwBFbcbyzgO9ExJ0RsTkiriD74jysYp1vRMSKiHiun/18ISLWRMQfgVvJxte2xQeAf4uIhyIbZ/lX4IDK1kpaXijHExE/iIjO9H5dBIwB9q3zeKcCX4mIxyNiPfBJ4OSqFtlnI+K5yMallvD852MTsJekqRGxPiLu2MbXPOw4qVgtVwLvIPsv+/vVC1N30a2SVktaS/ZlMbW8PCLWkP2n/yrgogGO9VTF9HM1no/firg7UzIqb1tr/33t72rgBEljgBOAeyLiCci6myT9LHWRPUv2ZTi1j/30Zw/g0NSNtEbSGrIvvpfWWllZ1+EdqatsDVmi25rj7gH8Q9XxZpCdjFG2oo79PFkx3c3W/U6q4/l6RSwFslbYrn3FI+njqbtsbdpmJ+p/D3YBnqh4/gRZi2fninl9vbYzgH2AhyUtlHRcnccc9pxU7EXSl+kfyL7EflJjlauB+cCMiNgJ+DbZlwMAkg4g63q4BvhGA0Pr5oUnDdT8Mt4WEbGU7EvnWLKEenXF4m8BD5OdWfUSsr78muMKQLGfGFcAv4qIiRWP8RHxweqdpOT2Y+DLwM4RMRG4sZ/j1rKCrOVWebyxEXFN5Uvfiv1tjVr7XQG8vyqeHSPi9lrbpfGTTwBvByal92Atz78HA8X+J7JEVrY7WRfpU7VXrwgi4rGIOAWYDlwIXC9p3EDbmZOK9e0M4I0RUayxbAJQiIgNkg4h+xIGQNIOwA/IvnjfA+wq6e8aFNNi4B2SRqaB2m3pgurP1cBHgNeTtbTKJpCdvrte0suBFyWBqhhPkDRW2UkFZ1Qs+xmwj6TTJI1Kj4MlvaLGfkaTdfWsBnokHQscPUD8T5GNH5RdCnwgtSwlaZyykywmDLCfRlgN9FbF823gk5JeCSBpJ0k1u/6SCWRJYDXQIek84CUVy58CZkrq63vsGuAcSbMkjSdrYV4XdZziLOmdkqZFRC+wJs3uHWg7c1KxPqS+/EV9LP474AJJ68gGxn9YsezfgBUR8a2I2Eg20Po5SXs3IKyPAH9D9kd+KvDTBuyzUnls6JaouE4H+DhZ4lxH9kV9XT/7+CpQIvvCu4LspAUAImIdWWI4mey/6CfJ/gseU72TtO7fk723Xen48weI/3zgitS99Pb0+zsT+GbaxzKyLs3cRUQ38HngdymewyLiBrLXe23qRnyArGXYlwXAz4FHyVqRG3hh91g58XdKuqfG9vPIunJ/Tdby3gB8uM6XcAzwoKT1wNeBkwcYd7LEFz+amVnDuKViZmYN46RiZmYN46RiZmYN46RiZmYNM6yLtk2dOjVmzpzZ6jDMzAaVu++++5mImFZr2bBOKjNnzmTRor7OmjUzs1okPdHXMnd/mZlZwzipmJlZwzipmJlZwzipmJlZwzipmJlZwzipmJlZwzipmJlZwwzr61RseFq/sYcrbl/Oxk2bB17ZbIg68hU785oZExu+XycVG3ZuffhpvrTgEQC0NfdRNBtCpr9kBycVs0Z4Zv1GAO79zFFMGje6xdGYDS0eU7Fhp1AsMUKw046jWh2K2ZDjpGLDTmexxKSxoxkxwn1fZo3mpGLDTlexxGR3e5nlwknFhp1OJxWz3Dip2LBTcFIxy42Tig07Tipm+XFSsWFlc2+wprvEFCcVs1w4qdiwsva5TfQGbqmY5cRJxYaVQjG78NEXPZrlw0nFhpXO9SUApowb0+JIzIYmJxUbVrq6s6Ti7i+zfDip2LDSWUwtlfFOKmZ5cFKxYaWQur8mjnXdL7M8OKnYsNJZLDFhTAdjOka2OhSzISnXpCLpGEmPSFom6dway8dIui4tv1PSzDR/tKTLJd0vaYmkORXb3Jb2uTg9plcse7ukpZIelHR1nq/NBqeu7hKT3fVllpvc7qciaSRwMXAUsBJYKGl+RCytWO0MoCsi9pJ0MnAhcBJwJkBE7J+Sxk2SDo6I3rTdqRGxqOp4ewOfBP4yIroqk41Zma+mN8tXni2VQ4BlEfF4RJSAa4G5VevMBa5I09cDR0oSsB9wC0BEPA2sAWYPcLwzgYsjoqtiO7MX6FxfYvJYJxWzvOSZVHYFVlQ8X5nm1VwnInqAtcAUYAlwvKQOSbOAg4AZFdtdnrq+PpOSEMA+wD6SfifpDknH1ApK0lmSFklatHr16u19jTbIuKVilq92HaifR5aEFgFfA24HNqdlp0bE/sAR6XFamt8B7A3MAU4BLpX0ohswR8QlETE7ImZPmzYt1xdh7SUiKHhMxSxXeSaVVbywdbFbmldzHUkdwE5AZ0T0RMQ5EXFARMwFJgKPAkTEqvRzHXA1WTcbZElofkRsiog/pPX3zuWV2aBULG2m1NPrYpJmOcozqSwE9pY0S9Jo4GRgftU684HT0/SJwC0REZLGShoHIOkooCcilqbusKlp/ijgOOCBtP1PyVoppHX2AR7P7dXZoFO+RmWSx1TMcpPb2V8R0SPpbGABMBKYFxEPSroAWBQR84HLgCslLQMKZIkHYDqwQFIvWWum3MU1Js0flfb5S+DStGwBcLSkpWRdZf8vIjrzen02+HSmYpK+mt4sP7klFYCIuBG4sWreeRXTG4C31dhuObBvjflFskH7WscK4GPpYfYiz9f9cjFJs7y060C9WcOVKxT7lGKz/Dip2LBRSMUkffaXWX6cVGzYKBRLjO4YwbjRrvtllhcnFRs2CsXs3vTPXy9rZo3mpGLDRqFY8unEZjlzUrFho7NY8unEZjlzUrFhw3W/zPLnpGLDRpeTilnunFRsWNjYs5l1G3t8jYpZzpxUbFjoKm4CfI2KWd6cVGxY2FL3y91fZrlyUrFhYUtLxXW/zHLlpGLDQrmlMnncqBZHYja0OanYsLCl7pdbKma5clKxYaFQLDFCMHFHt1TM8uSkYsNCuUTLiBGu+2WWJycVGxYKxRKTfOaXWe6cVGxY6PTV9GZN4aRiw0K57L2Z5ctJxYYF1/0yaw4nFRvyenuDrm4nFbNmcFKxIW/Nc5voDZxUzJrAScWGvMKWq+mdVMzy5qRiQ14h1f2a4qvpzXLnpGJDXrmlMsl1v8xy19HXAkmT+9swIgqND8es8TpT3S+3VMzy119L5W5gUfpZ/VhUz84lHSPpEUnLJJ1bY/kYSdel5XdKmpnmj5Z0uaT7JS2RNKdim9vSPhenx/Sqff6tpJA0u54YbegrrM+SilsqZvnrs6USEbO2Z8eSRgIXA0cBK4GFkuZHxNKK1c4AuiJiL0knAxcCJwFnphj2T0njJkkHR0Rv2u7UiHhRYpM0AfgIcOf2xG5DS6G7xIQxHYzpGNnqUMyGvAHHVJR5p6TPpOe7Szqkjn0fAiyLiMcjogRcC8ytWmcucEWavh44UpKA/YBbACLiaWANUE/L41/IEtOGOta1YcJ1v8yap56B+v8ADgfekZ6vI2uBDGRXYEXF85VpXs11IqIHWAtMAZYAx0vqkDQLOAiYUbHd5anr6zMpCSHpQGBGRPx3f0FJOkvSIkmLVq9eXcfLsMGu4KvpzZqmnqRyaER8iPTff0R0AXn/hc4jS0KLgK8BtwOb07JTI2J/4Ij0OE3SCOArwD8MtOOIuCQiZkfE7GnTpuUSvLWXzvWu+2XWLPUklU1pfCQAJE0DevvfBIBVvLB1sVuaV3MdSR3ATkBnRPRExDkRcUBEzAUmAo8CRMSq9HMdcDVZN9sE4FXAbZKWA4cB8z1Yb4BLtJg1UT1J5RvADcB0SZ8Hfgv8ax3bLQT2ljRL0mjgZGB+1TrzgdPT9InALRERksZKGgcg6SigJyKWpu6wqWn+KOA44IGIWBsRUyNiZkTMBO4Ajq81mG/DS0S47L1ZE/V59ldZRFwl6W7gSEDAWyLioTq265F0NrAAGAnMi4gHJV0ALIqI+cBlwJWSlgEFssQDMB1YIKmXrDVzWpo/Js0flfb5S+DS+l+uDTfF0mZKPb1OKmZNMmBSkfQN4NqIqGdw/gUi4kbgxqp551VMbwDeVmO75cC+NeYXyQbtBzrunK2N1Yam8jUqTipmzVFP99fdwD9J+r2kL3ucwgaTQreTilkzDZhUIuKKiHgzcDDwCHChpMdyj8ysAVyh2Ky5tqag5F7Ay4E9gIfzCcessTrXu+6XWTPVc0X9F1PL5ALgfmB2RPxN7pGZNUAhFZOcPN4tFbNmGHCgHvg9cHhEPJN3MGaNVuguMXrkCMaNdt0vs2aop/vrUuAYSefBVtX+Mmu5wvrsGpVUzcfMclZPUrmYrPbXKel5vbW/zFrOdb/Mmque7q9DI+JASfdCVvsrXSFv1vY6iyWmeDzFrGnyrP1l1nJd3SUmjXVSMWuWPGt/mbVceUzFzJojt9pfZq22sWcz6zb2uOy9WRP1mVQkTa54+jRwTeWyiCjkGZjZ9uoqbgJ8jYpZM/XXUrmbbBylfC5mpJ9K03vmGJfZdtty4aPHVMyaps+kEhGzmhmIWaNtSSru/jJrmq2p/WU2qHSmYpI+pdiseZxUbMh6vqXiYpJmzeKkYkNWV7GEBDvtOKrVoZgNG/0mFUkjJbnMvQ1KncXswseRI1z3y6xZ+k0qEbEZeETS7k2Kx6xhXPfLrPnqqf01CXhQ0l1AsTwzIo7PLSqzBuh0UjFrunqSymdyj8IsB13FEi+bNr7VYZgNK/Xco/5XwHJgVJpeCNyTc1xm261QLPlqerMmq+d2wmcC1wPfSbN2BX6aZ1Bm26u3N+jqLrnul1mT1XNK8YeAvwSeBYiIx4DpeQZltr3WPLeJ3vDV9GbNVk9S2RgRpfITSR08XwfMrC25RItZa9STVH4l6VPAjpKOAn4E/Fe+YZltHycVs9aoJ6mcC6wG7gfeD9wYEZ/ONSqz7VRIdb+cVMyaq56k8uGIuDQi3hYRJ0bEpZI+Us/OJR0j6RFJyySdW2P5GEnXpeV3SpqZ5o+WdLmk+yUtkTSnYpvb0j4Xp8f0NP9jkpZKuk/S/0jao54YbWjqTC2VKa77ZdZU9SSV02vMe/dAG6X72l8MHAvsB5wiab+q1c4AuiJiL+CrwIVp/pkAEbE/cBRwkaTKWE+NiAPS4+k0715gdkS8muxstS/W8dpsiOpKSWXSONf9MmumPpOKpFMk/RcwS9L8isetQD13fTwEWBYRj6eB/muBuVXrzAWuSNPXA0dKElkSugUgJY01wOz+DhYRt0ZEd3p6B7BbHTHaENVZLDF+TAdjOka2OhSzYaW/K+pvB/4MTAUuqpi/Drivjn3vCqyoeL4SOLSvdSKiR9JaYAqwBDhe0jXADOCg9POutN3lkjYDPwY+FxHVZ6OdAdxUKyhJZwFnAey+u0uaDVWu+2XWGv3d+fEJ4Ang8OaFs8U84BXAohTD7cDmtOzUiFglaQJZUjkN+H55Q0nvJGvVvKHWjiPiEuASgNmzZ/vU6CHKScWsNeq5ov4wSQslrZdUkrRZ0rN17HsVWeuibLc0r+Y66fqXnYDOiOiJiHPSmMlcYCLwKEBErEo/1wFXk3WzlWN9E/Bp4PiI2FhHjDZEOamYtUY9A/XfBE4BHgN2BN5HNgA/kIXA3pJmSRoNnAzMr1pnPs+fCHAicEtEhKSxksYBpGtjeiJiqaQOSVPT/FHAccAD6flryUrJHF8xeG/DlJOKWWvUU6WYiFgmaWS6v8rlku4FPjnANj2SzgYWACOBeRHxoKQLgEURMR+4DLhS0jKywf+T0+bTgQWSeslaM6el+WPS/FFpn78ELk3LvgSMB36UjfXzR5fnH54igs6i636ZtUI9SaU7tTQWS/oi2eB9XbchjogbgRur5p1XMb0BeFuN7ZYD+9aYXyQbtK91rDfVE5MNfd2lzZR6epnkpGLWdPUkh9PSemeT3aRrBvC3eQZltj1cosWsdeppqewFPB0RzwKfzTkes+32/NX0TipmzVZPS+VdwBJJd0j6kqS/kTQp78DMtpXrfpm1zoAtlYg4HUDSLmRnaF0M7FLPtmatUChuApxUzFphwMSQLiY8AtgfeIbsFOPf5ByX2TZzS8WsdeppbXwN+D3wbeDWdGaWWdvqLJYYPXIE48e4MW3WbAOOqUTEVOC9wA7A5yXdJenK3CMz20aF9dmFj+l6JTNronrKtLwE2B3YA5hJVkqlN9+wzLZdV3fJ16iYtUg9/QO/rXh8MyJW5huS2fbx1fRmrVPP2V+vbkYgZo1SKJaYMWlsq8MwG5bqKrdiNpiUx1TMrPmcVGxIKfX0sm5jj5OKWYs4qdiQ0tXtul9mrdTnmIqkfwf6vDNiRPx9LhGZbYfO9a77ZdZK/bVUFgF3k12fciDZTboeAw4A/BdrbckVis1aq7971F8BIOmDwOsioic9/zYu02JtquDuL7OWqmdMZRLwkorn49M8s7ZTWO+6X2atVM/Fj18A7pV0KyDg9cD5eQZltq0KxRISTBzrpGLWCvVc/Hi5pJuAQ9Osf4yIJ/MNy2zbdBZLTBo7mpEjXPfLrBXqqf0l4E3AayLiP4HRkg7JPTKzbdDVXWLS2FGtDsNs2KpnTOU/gMOBU9LzdWQ36jJrO53rS0wZN6bVYZgNW/UklUMj4kPABoCI6MKnFFubKhRdosWslepJKpskjSRdCClpGi59b22qUCwxebyTilmr1JNUvgHcAEyX9HmyEvj/mmtUZtugtzfo6i4x2Wd+mbVMPWd/XSXpbuBIslOK3xIRD+UemdlWWvvcJnrD16iYtVK9N/F+DHi2vL6k3SPij7lFZbYNOlOJlinu/jJrmQGTiqQPA/8MPAVsJmutBOCbd1lbcd0vs9arZ0zlI8C+EfHKiHh1ROxf790gJR0j6RFJyySdW2P5GEnXpeV3SpqZ5o+WdLmk+yUtkTSnYpvb0j4Xp8f0/vZlw0c5qUzymIpZy9STVFYAa7d2x+mMsYuBY4H9gFMk7Ve12hlAV0TsBXwVuDDNPxMgIvYHjgIuklQZ66kRcUB6PD3AvmyYKLj7y6zl+rufysfS5OPAbZL+G9hYXh4RXxlg34cAyyLi8bS/a4G5wNKKdebyfB2x64Fvpiv49wNuScd5WtIaYDZwVz/Hq7mviOjznjA2tBSKLiZp1mr9tVQmpMcfgZvJLnicUPEYyK5krZyylWlezXVSaf21wBRgCXC8pA5Js4CDgBkV212eur4+k5JQf/t6AUlnSVokadHq1avreBk2WHQWS4wf08GYjpGtDsVs2OrvfiqfbWYgVeYBryC7UdgTwO1kJwlA1vW1StIE4MfAacD3691xRFwCXAIwe/Zst2KGkK5iiUnjXPfLrJXqKSh5s6SJFc8nSVpQx75X8cLWxW5pXs11JHUAOwGdEdETEeekMZO5wETgUYCIWJV+rgOuJutm63NfdcRpQ0RnscRk1/0ya6l6BuqnRcSa8pNU+2t6HdstBPaWNEvSaOBkYH7VOvOB09P0icAtERGSxkoaByDpKKAnIpam7rCpaf4o4Djggf72VUecNkQUiiXfm96sxeq5+HFz5cWOkvYg1QHrT0T0SDobWACMBOZFxIOSLgAWRcR84DLgSknLgAJZ4oEsaS2Q1EvWAjktzR+T5o9K+/wlcGla1te+bJgoFEu8/KUvGXhFM8tNPUnl08BvJf2K7MLHI4Cz6tl5RNwI3Fg177yK6Q3A22pstxzYt8b8Itmgfa1j1dyXDQ8RkbVUfDqxWUvVU/vr55IOBA5Lsz4aEc/kG5bZ1ukubWZjT69PJzZrsXprf20GngZ2APaTRET8Or+wzLaOS7SYtYd6an+9j6xUy27AYrIWy/8Cb8w3NLP6lYtJuuy9WWvVW/vrYOCJiPgr4LXAmv43MWuurnJS8ZiKWUvVk1Q2pEFwJI2JiIepMYhu1kpbyt67+8uspeoZU1mZLn78KXCzpC6yq9zN2obrfpm1h3rO/nprmjxf0q1kV6r/PNeozLZSZ7HEqJFi/Jh6zz0xszzU9Rco6XXA3hFxuaRpZMUb/5BrZGZboatYYvK40TxfX9TMWqGe2l//DPwj8Mk0axTwgzyDMttaBdf9MmsL9QzUvxU4HigCRMSfqK/0vVnTdLrul1lbqCeplFJhxgAoF3o0ayeFYolJTipmLVdPUvmhpO8AEyWdSVbE8bv5hmW2dVyh2Kw91HP215dT+flnya5POS8ibs49MrM6lXp6Wbehx6cTm7WBus7+SknkZgBJIySdGhFX5RqZWZ26ul33y6xd9Nn9Jeklkj4p6ZuSjlbmbOBx4O3NC9Gsf53rnVTM2kV/LZUrgS6y4pHvAz5Fdj+Vt0TE4ibEZlYXt1TM2kd/SWXPiNgfQNJ3gT8Du5frgJm1C9f9Mmsf/Z39tak8ERGbgZVOKNaOCutd98usXfTXUnmNpGfTtIAd03MBERG+Gbi1hUKxhAQTfS8Vs5brM6lExMhmBmK2rQrdJSbuOIqRI1z3y6zV6rn40aytFVIxSTNrPScVG/Q615eY4mKSZm3BScUGvazu16hWh2FmOKnYENDV7bL3Zu3CScUGtd7eoKt7k69RMWsTTio2qK19bhObe8MD9WZtwknFBrXy1fROKmbtIdekIukYSY9IWibp3BrLx0i6Li2/U9LMNH+0pMsl3S9piaQ5NbadL+mBiucHSLpD0mJJiyQdkuNLszbhul9m7SW3pCJpJHAxcCywH3CKpP2qVjsD6IqIvYCvAhem+WcCpNpjRwEXSdoSq6QTgPVV+/oi8NmIOAA4Lz23Ic4Vis3aS54tlUOAZRHxeESUgGuBuVXrzAWuSNPXA0dKElkSugUgIp4G1gCzASSNBz4GfK5qXwGUS8fsBPypoa/G2lKhXExyvJOKWTvIM6nsCqyoeL4yzau5TkT0AGuBKcAS4HhJHZJmAQcBM9I2/wJcBHRX7eujwJckrQC+DHyyVlCSzkrdY4tWr169ra/N2kShmBWTnOS6X2ZtoV0H6ueRJaFFwNeA24HNkg4AXhYRN9TY5oPAORExAzgHuKzWjiPikoiYHRGzp02blk/01jSF4ibGjR7JDqNcqs6sHeSZVFbxfOsCYLc0r+Y6kjrIuq06I6InIs6JiAMiYi4wEXgUOByYLWk58FtgH0m3pX2dDvwkTf+IrPvNhrhCcSOT3fVl1jbyTCoLgb0lzZI0GjgZmF+1znyyZABwInBLRISksZLGAUg6CuiJiKUR8a2I2CUiZgKvAx6NiDlp+z8Bb0jTbwQey+uFWfvoLJaY7K4vs7bR3/1UtktE9KR72i8ARgLzIuJBSRcAiyJiPlkX1ZWSlgEFssQDMB1YIKmXrDVzWh2HPBP4emrxbADOauwrsnZUKJaYPsElWszaRW5JBSAibgRurJp3XsX0BuBtNbZbDuw7wL6XA6+qeP5bsgF9G0a6iiVe/lLfL86sXbTrQL3ZgCKCzmLJpxObtREnFRu0ukub2djT69OJzdqIk4oNWlsufPTV9GZtw0nFBq2Ci0matR0nFRu0tiQVj6mYtQ0nFRu0tpS995iKWdtwUrFBq1z3yy0Vs/bhpGKDVqG4iVEjxYQxuV5uZWZbwUnFBq1CcSOTx40muwYJch8AAAqYSURBVFuCmbUDJxUbtArFkq9RMWszTio2aPlqerP246Rig1ZXscTkcS4madZOnFRs0Ooslnw1vVmbcVKxQanU08u6DT0eUzFrM04qNih1dftqerN25KRig5KLSZq1JycVG5RcTNKsPTmp2KDU6aRi1pacVGxQKqxPdb+cVMzaipOKDUqF7k1I+OwvszbjpGKDUqG4kYk7jmLkCNf9MmsnTio2KBWKJSa568us7Tip2KDUud5X05u1IycVG5S6uksepDdrQ04qNigVXEzSrC05qdig09sbdHVvYvK4Ua0Oxcyq5JpUJB0j6RFJyySdW2P5GEnXpeV3SpqZ5o+WdLmk+yUtkTSnxrbzJT1QNe/Dkh6W9KCkL+b0sqzF1j63ic294ZaKWRvK7ebekkYCFwNHASuBhZLmR8TSitXOALoiYi9JJwMXAicBZwJExP6SpgM3STo4InrTvk8A1lcd76+AucBrImJj2s6GoEK3636Ztas8WyqHAMsi4vGIKAHXkn3pV5oLXJGmrweOVHbD8f2AWwAi4mlgDTAbQNJ44GPA56r29UHgCxGxsWI7G4Jc98usfeXWUgF2BVZUPF8JHNrXOhHRI2ktMAVYAhwv6RpgBnBQ+nkX8C/ARUB31b72AY6Q9HlgA/DxiFhYHZSks4CzAHbfffdtemE/XLiCS3/z+DZta9uvuLEHcFIxa0d5JpXtMQ94BbAIeAK4Hdgs6QDgZRFxTnn8pUIHMBk4DDgY+KGkPSMiKleKiEuASwBmz54dbIOJY0ex987jt2VTa5A5Y0ezz84TWh2GmVXJM6msImtdlO2W5tVaZ6WkDmAnoDMlgnPKK0m6HXgUeAMwW9JystinS7otIuaQtYR+kra9S1IvMBVY3egXdvQrX8rRr3xpo3drZjbo5TmmshDYW9IsSaOBk4H5VevMB05P0ycCt0RESBoraRyApKOAnohYGhHfiohdImIm8Drg0ZRQAH4K/FXaZh9gNPBMfi/PzMyq5dZSSWMkZwMLgJHAvIh4UNIFwKKImA9cBlwpaRlQIEs8ANOBBam1sQo4rY5DzgPmpdOMS8Dp1V1fZmaWLw3n793Zs2fHokWLWh2GmdmgIunuiJhda5mvqDczs4ZxUjEzs4ZxUjEzs4ZxUjEzs4ZxUjEzs4YZ1md/SVpNdsX+tphKe18H4/i2j+Pbfu0eo+PbdntExLRaC4Z1Utkekhb1dUpdO3B828fxbb92j9Hx5cPdX2Zm1jBOKmZm1jBOKtvuklYHMADHt30c3/Zr9xgdXw48pmJmZg3jloqZmTWMk4qZmTWMk8oAJB0j6RFJyySdW2P5GEnXpeV31rgjZZ6xzZB0q6Slkh6U9JEa68yRtFbS4vQ4r1nxpeMvl3R/OvaLSkIr8430/t0n6cAmxrZvxfuyWNKzkj5atU7T3z9J8yQ9nW7jUJ43WdLNkh5LPyf1se3paZ3HJJ1ea50cYvuSpIfT7+8GSRP72Lbfz0LOMZ4vaVXF7/HNfWzb7997jvFdVxHbckmL+9i2Ke/hdokIP/p4kN0H5vfAnmQ3/VoC7Fe1zt8B307TJwPXNTG+vwAOTNMTyO6OWR3fHOBnLXwPlwNT+1n+ZuAmQGS3gr6zhb/rJ8ku6mrp+we8HjgQeKBi3heBc9P0ucCFNbabDDyefk5K05OaENvRQEeavrBWbPV8FnKO8Xzg43V8Bvr9e88rvqrlFwHntfI93J6HWyr9OwRYFhGPR0QJuBaYW7XOXOCKNH09cKQkNSO4iPhzRNyTptcBDwG7NuPYDTQX+H5k7gAmSvqLFsRxJPD7iNjWCgsNExG/JrtpXaXKz9kVwFtqbPp/gZsjohARXcDNwDF5xxYRv4iInvT0DrJbh7dMH+9fPer5e99u/cWXvjveDlzT6OM2i5NK/3YFVlQ8X8mLv7S3rJP+sNYCU5oSXYXU7fZa4M4aiw+XtETSTZJe2dTAIIBfSLpb0lk1ltfzHjfDyfT9h9zK969s54j4c5p+Eti5xjrt8F6+l6zlWctAn4W8nZ266Ob10X3YDu/fEcBTEfFYH8tb/R4OyEllCJA0Hvgx8NGIeLZq8T1kXTqvAf4d+GmTw3tdRBwIHAt8SNLrm3z8AUkaDRwP/KjG4la/fy8SWT9I210LIOnTQA9wVR+rtPKz8C3gZcABwJ/Jupja0Sn030pp+78nJ5X+rQJmVDzfLc2ruY6kDmAnoLMp0WXHHEWWUK6KiJ9UL4+IZyNifZq+ERglaWqz4ouIVenn08ANZF0Mlep5j/N2LHBPRDxVvaDV71+Fp8rdgunn0zXWadl7KendwHHAqSnpvUgdn4XcRMRTEbE5InqBS/s4dks/i+n74wTgur7WaeV7WC8nlf4tBPaWNCv9N3syML9qnflA+SybE4Fb+vqjarTU/3oZ8FBEfKWPdV5aHuORdAjZ77wpSU/SOEkTytNkA7oPVK02H3hXOgvsMGBtRTdPs/T532Er378qlZ+z04H/rLHOAuBoSZNS987RaV6uJB0DfAI4PiK6+1inns9CnjFWjtO9tY9j1/P3nqc3AQ9HxMpaC1v9Htat1WcKtPuD7OykR8nOCvl0mncB2R8QwA5k3SbLgLuAPZsY2+vIukHuAxanx5uBDwAfSOucDTxIdibLHcD/aWJ8e6bjLkkxlN+/yvgEXJze3/uB2U3+/Y4jSxI7Vcxr6ftHluD+DGwi69c/g2yc7n+Ax4BfApPTurOB71Zs+970WVwGvKdJsS0jG4sofwbLZ0PuAtzY32ehie/flenzdR9ZoviL6hjT8xf9vTcjvjT/e+XPXcW6LXkPt+fhMi1mZtYw7v4yM7OGcVIxM7OGcVIxM7OGcVIxM7OGcVIxM7OGcVIx2w6Sbk8/Z0p6R4P3/alaxzJrZz6l2KwBJM0hq4J73FZs0xHPF2KstXx9RIxvRHxmzeKWitl2kLQ+TX4BOCLd5+IcSSPTfUYWpiKG70/rz5H0G0nzgaVp3k9TgcAHy0UCJX0B2DHt76rKY6XqA1+S9EC6t8ZJFfu+TdL1yu5vclVFNYAvKLvvzn2SvtzM98iGl45WB2A2RJxLRUslJYe1EXGwpDHA7yT9Iq17IPCqiPhDev7eiChI2hFYKOnHEXGupLMj4oAaxzqBrDDia4CpaZtfp2WvBV4J/An4HfCXkh4iK03y8ogI9XETLbNGcEvFLB9Hk9U0W0x2O4IpwN5p2V0VCQXg7yWVy8DMqFivL68DromsQOJTwK+Agyv2vTKywomLgZlkt2PYAFwm6QSgZn0us0ZwUjHLh4APR8QB6TErIsotleKWlbKxmDcBh0dWXv9esnpy22pjxfRmsjsy9pBVs72erJLwz7dj/2b9clIxa4x1ZLd0LlsAfDDdmgBJ+6TKstV2AroiolvSy8luqVy2qbx9ld8AJ6Vxm2lkt6e9q6/A0v12doqsdP85ZN1mZrnwmIpZY9wHbE7dWN8Dvk7W9XRPGixfTe1bAP8c+EAa93iErAus7BLgPkn3RMSpFfNvAA4nq1YbwCci4smUlGqZAPynpB3IWlAf27aXaDYwn1JsZmYN4+4vMzNrGCcVMzNrGCcVMzNrGCcVMzNrGCcVMzNrGCcVMzNrGCcVMzNrmP8PMwJoVZqeuu8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.title('Maximum value after n iterations')\n",
    "plt.plot(np.maximum.accumulate(bo.y))\n",
    "plt.ylabel('Reached water level')\n",
    "plt.xlabel('iterations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7ff625609358>"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD4CAYAAADlwTGnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAYxUlEQVR4nO3df4yc1X3v8fenXtt13RQDXnrD2lybYqG6tcqPiUMvCcoNF2xHEUuRae0i4rQobpVauv1F66hqErlXClyai5TWSuMUXIJCcEoTsmpJHFo3uleJcD0GY2Mch43r4N3SeIttoCmVWfjeP+ZsMgyzu896dn6ez0sa7TPnOc8z55l9Zj7znHNmVxGBmZnl58fa3QAzM2sPB4CZWaYcAGZmmXIAmJllygFgZpapvnY3YCYWL14cy5Yta3czzMy6yv79+/8tIvpry7sqAJYtW0a5XG53M8zMuoqk79UrdxeQmVmmHABmZpkqFACS1ko6KmlY0tY666+T9KSkcUnra9Z9TdIZSX9bU75c0t60z12S5jV2KGZmNhPTBoCkOcB2YB2wEtgoaWVNteeBDwIP1dnFPcDtdcrvBu6NiMuA08AdxZttZmaNKnIFsBoYjohjEXEWeBgYrK4QEccj4iDwRu3GEfEPwCvVZZIEvBd4JBU9ANw88+abmdm5KhIAA8CJqvsjqawRFwJnImJ8un1K2iypLKk8NjbW4MOamdmEjh8EjogdEVGKiFJ//1umsZqZ2TkqEgCjwNKq+0tSWSNeBBZJmvgewmzs08zMZqBIAOwDVqRZO/OADcBQIw8alX9C8I/AxIyhTcBXGtmnmZnNzLQBkPrptwC7gSPAFyPisKRtkm4CkPQOSSPArcBnJB2e2F7S/wP+Grhe0oikNWnVHwK/K2mYypjAfbN5YGZmNjV1038EK5VK4T8FYWY2M5L2R0SptrzjB4HNzKw5HABmZplyAJiZZcoBYGaWKQeAmVmmHABmZplyAJiZZcoBYGaWKQeAmVmmHABmZplyAJiZZcoBYGaWKQeAmVmmHABmZplyAJiZZcoBYGaWKQeAmVmmHABmZplyAJiZZcoBYGaWqUIBIGmtpKOShiVtrbP+OklPShqXtL5m3SZJz6Xbpqryb6R9Hki3ixo/HDMzK6pvugqS5gDbgRuAEWCfpKGIeLaq2vPAB4Hfr9n2AuBjQAkIYH/a9nSqcltElBs+CjMzm7EiVwCrgeGIOBYRZ4GHgcHqChFxPCIOAm/UbLsGeDwiTqU3/ceBtbPQbjMza1CRABgATlTdH0llRUy37c7U/fPHklRvB5I2SypLKo+NjRV8WDMzm047B4Fvi4hVwLvT7fZ6lSJiR0SUIqLU39/f0gaamfWyIgEwCiytur8klRUx6bYRMfHzFeAhKl1NZmbWIkUCYB+wQtJySfOADcBQwf3vBm6UdL6k84Ebgd2S+iQtBpA0F3g/8MzMm29mZudq2gCIiHFgC5U38yPAFyPisKRtkm4CkPQOSSPArcBnJB1O254C/oRKiOwDtqWy+VSC4CBwgMpVwWdn/ejMzGxSioh2t6GwUqkU5bJnjZqZzYSk/RFRqi33N4HNzDLlADAzy5QDwMwsUw4AM7NMOQDMzDLlADAzy5QDwMwsUw4AM7NMOQDMzDLlADAzy5QDwMwsUw4AM7NMOQDMzDLlADAzy5QDwMwsUw4AM7NMOQDMzDLlADAzy5QDwMwsU4UCQNJaSUclDUvaWmf9dZKelDQuaX3Nuk2Snku3TVXlV0s6lPb5KUlq/HDMzKyoaQNA0hxgO7AOWAlslLSyptrzwAeBh2q2vQD4GPBOYDXwMUnnp9WfBj4ErEi3ted8FGZmNmNFrgBWA8MRcSwizgIPA4PVFSLieEQcBN6o2XYN8HhEnIqI08DjwFpJbwd+KiKeiIgAPgfc3OjBmJlZcUUCYAA4UXV/JJUVMdm2A2l52n1K2iypLKk8NjZW8GHNzGw6HT8IHBE7IqIUEaX+/v52N8fMrGcUCYBRYGnV/SWprIjJth1Ny+eyTzMzmwVFAmAfsELScknzgA3AUMH97wZulHR+Gvy9EdgdES8AL0u6Js3++QDwlXNov5mZnaNpAyAixoEtVN7MjwBfjIjDkrZJuglA0jskjQC3Ap+RdDhtewr4Eyohsg/YlsoAPgz8JTAMfBf46qwemZmZTUmVSTjdoVQqRblcbnczzMy6iqT9EVGqLe/4QWAzM2sOB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmeprdwPMpvLoU6Pcs/so/3LmVS5etIA711zOzVcOtLtZZj3BAWAd69GnRvnIlw7x6muvAzB65lU+8qVDAA4Bs1ngLiDrWPfsPvrDN/8Jr772OvfsPtqmFpn1FgeAdax/OfPqjMrNbGYKBYCktZKOShqWtLXO+vmSdqX1eyUtS+XzJO2UdEjS05LeU7XNN9I+D6TbRbN0TNYjLl60YEblZjYz0waApDnAdmAdsBLYKGllTbU7gNMRcRlwL3B3Kv8QQESsAm4APimp+jFvi4gr0u1kY4divebONZezYO6cN5UtmDuHO9dc3qYWdZ9Hnxrl2rv2sHzr33HtXXt49KnRdjfJOkiRQeDVwHBEHAOQ9DAwCDxbVWcQ+HhafgT4c0miEhh7ACLipKQzQAn4p1lpvfW0iYHebp4F1M5ZTB5Et+kUCYAB4ETV/RHgnZPViYhxSS8BFwJPAzdJ+gKwFLg6/ZwIgJ2SXgf+BvhfERG1Dy5pM7AZ4JJLLil4WD/iaYTd7eYrB9r6+2rk/Gn3G/BUg+h+DRg0fxro/cDPAmXge8C3gIkz8raIGJX0NioBcDvwudodRMQOYAdAqVR6S0BMpd0vQOtujZ4/7X4D9iB643r9A2SRQeBRKp/aJyxJZXXrSOoDzgNejIjxiPid1Mc/CCwCvgMQEaPp5yvAQ1S6mmaVpxFaIxo9f9r9BuxB9MZMfAAYPfMqwY8+APTSOEqRANgHrJC0XNI8YAMwVFNnCNiUltcDeyIiJP2EpIUAkm4AxiPiWUl9khan8rnA+4FnZuF43qTdL0Drbo2eP+1+A/YgemNy+AA5bQBExDiwBdgNHAG+GBGHJW2TdFOqdh9woaRh4HeBiamiFwFPSjoC/CGVbh6A+cBuSQeBA1SuID47S8f0Q+1+AVp3a/T8afcb8M1XDvCJW1YxsGgBAgYWLeATt6zqqS6MZsrhA2ShMYCIeAx4rKbso1XL/wncWme748BbzvaI+AGVAeGmunPN5W/qwwV/AspNI324jZ4/nTCLqd2D6N3s4kULGK3zZt9LHyB7+m8BdcIL0Nqn0UHc2Th//AbcvXL4AKk6My87VqlUinK53O5mWJe49q49dT/BDSxawDe3vrcNLbJu0yuzgCTtj4hSbXlPXwHMhl45AXLUC324Pv/aq9ev4BwAU/D3CLr7Dajb+3B9/lmz+a+BTiGHaWBT6fZ50O2ehdOo3M8/az5fAUyhF7oQGtHub7I2qtsnAeR+/vWCTr+CdgBModu7EBrVC29A3dyHm/v51+26oQvPXUBT6PYuhEb5i3Ttlfv51+26oQvPVwBT6PYuhEblMA+6k3XC+dfpXRidrBuuoB0A0+jmLoRGdcIbUO7aef51QxdGJ+uGLjwHgE0p5wDMXbdPAmi3briCdgB0OF+CW7t0QxdGJ+uGK2gHQAfzJbi1Uzd0YXS6Tr+C9iygDtYNswisd3kWUu/zFUAH64VLcHdhda9u6MLodc1+/TgAOli3X4K7C6v7dXoXRi9rxevHXUAdrNsvwd2FZXbuWvH68RVAB+v2S/Be6MKyxrgL8Ny14vXjAOhw3XwJ3u1dWNYYdwE2phWvn0JdQJLWSjoqaVjS1jrr50valdbvlbQslc+TtFPSIUlPS3pP1TZXp/JhSZ+SpFk6JusQ3d6FZY1xF2BjWvH6mTYAJM0BtgPrgJXARkkra6rdAZyOiMuAe4G7U/mHACJiFXAD8ElJE4/56bR+RbqtbexQrNPcfOUAn7hlFQOLFiAq/4rxE7es8qe/TLgLsDGteP0U6QJaDQxHxDEASQ8Dg8CzVXUGgY+n5UeAP0+f6FcCewAi4qSkM0BJ0gngpyLiibTPzwE3A19t+Iiso3RzF5Y1xl2AjWv266dIF9AAcKLq/kgqq1snIsaBl4ALgaeBmyT1SVoOXA0sTfVHptmnmXUxdwF2vmYPAt8P/CxQBr4HfAt4fcotakjaDGwGuOSSS2a7fT3PszCsXbp9FlsOigTAKJVP7ROWpLJ6dUYk9QHnAS9GRAC/M1FJ0reA7wCn036m2icAEbED2AFQKpWiQHst8SwMazd3AXa2Il1A+4AVkpZLmgdsAIZq6gwBm9LyemBPRISkn5C0EEDSDcB4RDwbES8AL0u6Jo0VfAD4ymwckP2IZ2GY2VSmvQKIiHFJW4DdwBzg/og4LGkbUI6IIeA+4EFJw8ApKiEBcBGwW9IbVD7h31616w8DfwUsoDL46wHgWeZZGGY2lUJjABHxGPBYTdlHq5b/E7i1znbHgbojPhFRBn5+Bm21GfIsDOt2HsNqLv8toB7mWRjWzSbGsEbPvErwozGsR5+qO1w46T6uvWsPy7f+HdfetWdG2+bAfwqiydr5CcazMKybNfovKT0JYnoOgCbqhBPQszCsWzU6huX/aTw9dwE1kWfhmJ27ycaqio5heRLE9BwATeQT0OzcNTqG1WiA5MAB0EQ+Ac3OXaN/DM2TIKbnMYAmunPN5W8aAwCfgGYz0cgYlidBTM8B0EQ+Ac3ay5MgpuYAaDKfgGbWqTwGYGaWKQeAmVmmHABmZplyAJiZZcoBYGaWKQeAmVmmHABmZplyAJiZZcoBYGaWKQeAmVmmHABmZpkqFACS1ko6KmlY0tY66+dL2pXW75W0LJXPlfSApEOSjkj6SNU2x1P5AUnl2TogMzMrZtoAkDQH2A6sA1YCGyWtrKl2B3A6Ii4D7gXuTuW3AvMjYhVwNfAbE+GQ/PeIuCIiSg0dhZmZzViRK4DVwHBEHIuIs8DDwGBNnUHggbT8CHC9JAEBLJTUBywAzgIvz0rLzcysIUUCYAA4UXV/JJXVrRMR48BLwIVUwuAHwAvA88CfRsSptE0AX5e0X9LmyR5c0mZJZUnlsbGxAs01M7Mimj0IvBp4HbgYWA78nqRL07p3RcRVVLqWfkvSdfV2EBE7IqIUEaX+/v4mN9fMLB9FAmAUWFp1f0kqq1sndfecB7wI/CrwtYh4LSJOAt8ESgARMZp+ngS+TCUszMysRYoEwD5ghaTlkuYBG4ChmjpDwKa0vB7YExFBpdvnvQCSFgLXAN+WtFDS26rKbwSeafRgzMysuGn/JWREjEvaAuwG5gD3R8RhSduAckQMAfcBD0oaBk5RCQmozB7aKekwIGBnRBxM3UBfrowT0wc8FBFfm+2DMzOzyanyQb07lEqlKJf9lQEzs5mQtL/edHt/E9jMLFMOADOzTDkAzMwy5QAwM8uUA8DMLFMOADOzTDkAzMwy5QAwM8uUA8DMLFMOADOzTDkAzMwy5QAwM8uUA8DMLFMOADOzTDkAzMwy5QAwM8uUA8DMLFMOADOzTDkAzMwyVSgAJK2VdFTSsKStddbPl7Qrrd8raVkqnyvpAUmHJB2R9JGi+zQzs+aaNgAkzQG2A+uAlcBGSStrqt0BnI6Iy4B7gbtT+a3A/IhYBVwN/IakZQX3aWZmTVTkCmA1MBwRxyLiLPAwMFhTZxB4IC0/AlwvSUAACyX1AQuAs8DLBfdpZmZNVCQABoATVfdHUlndOhExDrwEXEglDH4AvAA8D/xpRJwquE8AJG2WVJZUHhsbK9BcMzMrotmDwKuB14GLgeXA70m6dCY7iIgdEVGKiFJ/f38z2mhmlqUiATAKLK26vySV1a2TunvOA14EfhX4WkS8FhEngW8CpYL7NDOzJioSAPuAFZKWS5oHbACGauoMAZvS8npgT0QElW6f9wJIWghcA3y74D7NzKyJ+qarEBHjkrYAu4E5wP0RcVjSNqAcEUPAfcCDkoaBU1Te0KEy02enpMOAgJ0RcRCg3j5n+djMzGwKqnxQ7w6lUinK5XK7m2Fm1lUk7Y+IUm25vwlsZpYpB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmSoUAJLWSjoqaVjS1jrr50valdbvlbQsld8m6UDV7Q1JV6R130j7nFh30WwemJmZTW3aAJA0B9gOrANWAhslraypdgdwOiIuA+4F7gaIiM9HxBURcQVwO/DPEXGgarvbJtZHxMlZOB4zMyuoyBXAamA4Io5FxFngYWCwps4g8EBafgS4XpJq6mxM25qZWQcoEgADwImq+yOprG6diBgHXgIurKnzK8AXasp2pu6fP64TGABI2iypLKk8NjZWoLlmZlZESwaBJb0T+I+IeKaq+LaIWAW8O91ur7dtROyIiFJElPr7+1vQWjOzPBQJgFFgadX9Jamsbh1JfcB5wItV6zdQ8+k/IkbTz1eAh6h0NZmZWYsUCYB9wApJyyXNo/JmPlRTZwjYlJbXA3siIgAk/Rjwy1T1/0vqk7Q4Lc8F3g88g5mZtUzfdBUiYlzSFmA3MAe4PyIOS9oGlCNiCLgPeFDSMHCKSkhMuA44ERHHqsrmA7vTm/8c4O+Bz87KEZmZWSFKH9S7QqlUinK53O5mmJl1FUn7I6JUW+5vApuZZcoBYGaWKQeAmVmmHABmZplyAJiZZcoBYGaWKQeAmVmmuup7AJLGgO+d4+aLgX+bxebMNrevMW5fY9y+xnR6+/5rRLzlj6l1VQA0QlK53hchOoXb1xi3rzFuX2M6vX2TcReQmVmmHABmZpnKKQB2tLsB03D7GuP2Ncbta0ynt6+ubMYAzMzszXK6AjAzsyoOADOzTPVcAEhaK+mopGFJW+usny9pV1q/V9KyFrZtqaR/lPSspMOS/medOu+R9JKkA+n20Va1Lz3+cUmH0mO/5Z8vqOJT6fk7KOmqFrbt8qrn5YCklyX9dk2dlj5/ku6XdFLSM1VlF0h6XNJz6ef5k2y7KdV5TtKmenWa1L57JH07/f6+LGnRJNtOeS40sX0flzRa9Tt83yTbTvlab2L7dlW17bikA5Ns2/Tnr2ER0TM3Kv9d7LvApcA84GlgZU2dDwN/kZY3ALta2L63A1el5bcB36nTvvcAf9vG5/A4sHiK9e8DvgoIuAbY28bf9b9S+YJL254/Kv/x7irgmaqy/w1sTctbgbvrbHcBcCz9PD8tn9+i9t0I9KXlu+u1r8i50MT2fRz4/QK//ylf681qX836TwIfbdfz1+it164AVgPDEXEsIs5S+T/EgzV1BoEH0vIjwPWS1IrGRcQLEfFkWn4FOAIMtOKxZ9Eg8LmoeAJYJOntbWjH9cB3I+Jcvxk+KyLi/1L5N6jVqs+xB4Cb62y6Bng8Ik5FxGngcWBtK9oXEV+PiPF09wlgyWw/blGTPH9FFHmtN2yq9qX3jV8GvjDbj9sqvRYAA8CJqvsjvPUN9od10ovgJeDClrSuSup6uhLYW2f1L0p6WtJXJf1cSxsGAXxd0n5Jm+usL/Ict8IGJn/htfP5A/jpiHghLf8r8NN16nTK8/jrVK7o6pnuXGimLamL6v5JutA64fl7N/D9iHhukvXtfP4K6bUA6AqSfhL4G+C3I+LlmtVPUunW+AXgz4BHW9y8d0XEVcA64LckXdfix5+WpHnATcBf11nd7ufvTaLSF9CRc60l/REwDnx+kirtOhc+DfwMcAXwApVulk60kak//Xf8a6nXAmAUWFp1f0kqq1tHUh9wHvBiS1pXecy5VN78Px8RX6pdHxEvR8S/p+XHgLmSFreqfRExmn6eBL5M5VK7WpHnuNnWAU9GxPdrV7T7+Uu+P9Etln6erFOnrc+jpA8C7wduSyH1FgXOhaaIiO9HxOsR8Qbw2Uket93PXx9wC7Brsjrtev5motcCYB+wQtLy9ClxAzBUU2cImJhxsR7YM9kLYLalPsP7gCMR8X8mqfNfJsYkJK2m8jtqSUBJWijpbRPLVAYLn6mpNgR8IM0GugZ4qaq7o1Um/eTVzuevSvU5tgn4Sp06u4EbJZ2fujhuTGVNJ2kt8AfATRHxH5PUKXIuNKt91WNKvzTJ4xZ5rTfT/wC+HREj9Va28/mbkXaPQs/2jcosle9QmSHwR6lsG5WTHeDHqXQdDAP/BFzawra9i0p3wEHgQLq9D/hN4DdTnS3AYSqzGp4A/lsL23dpetynUxsmnr/q9gnYnp7fQ0Cpxb/fhVTe0M+rKmvb80cliF4AXqPSD30HlTGlfwCeA/4euCDVLQF/WbXtr6fzcBj4tRa2b5hK//nEOTgxK+5i4LGpzoUWte/BdG4dpPKm/vba9qX7b3mtt6J9qfyvJs65qrotf/4avflPQZiZZarXuoDMzKwgB4CZWaYcAGZmmXIAmJllygFgZpYpB4CZWaYcAGZmmfr/8iSD6re6AYQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(np.arange(bo.y.shape[0]), bo.y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting analysis of 8.919781953086979 m3\n",
      "Convergence in 5 iterations.\n",
      "Starting analysis of 9.252438621002577 m3\n",
      "Convergence in 3 iterations.\n",
      "Starting analysis of 9.563312669021387 m3\n",
      "Convergence in 3 iterations.\n",
      "Starting analysis of 9.85383044020318 m3\n",
      "Convergence in 3 iterations.\n",
      "Starting analysis of 10.125324879543868 m3\n",
      "Convergence in 3 iterations.\n",
      "Starting analysis of 10.379041649754082 m3\n",
      "Convergence in 4 iterations.\n",
      "Starting analysis of 10.616144846571766 m3\n",
      "Convergence in 39 iterations.\n",
      "Breaking. Water level isn't rising.\n"
     ]
    }
   ],
   "source": [
    "ss, cubics, water_heights = determine_max_water_cubic(bo.best_params, return_ss=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ff62427fc18>]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD7CAYAAABkO19ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU9d3+8feHQAJJlCUBVBZBARGtCkbA3bpibaW/ihWwLW2p2AoubbVF26q1z9OndrOt0IWn2vooq6g1KmrdutkWCUF2kMiWsK+BJIRsn98fM+CYBhlg1jP367q4mDnnZM73KNc9J/d8zxlzd0REJLhaJXsAIiISXwp6EZGAU9CLiAScgl5EJOAU9CIiAaegFxEJuKiC3syGmdlKMyszs4ktrL/EzErNrMHMRjRbN8bMVoX/jInVwEVEJDp2uHn0ZpYFvAdcBVQA84BR7r4sYptewPHA3UCxu88OL+8ElABFgAPzgXPdfVesD0RERFrWOoptBgNl7r4awMxmAMOBg0Hv7mvD65qa/ew1wGvuvjO8/jVgGDD9UDsrLCz0Xr16RX8EIiLC/Pnzt7t755bWRRP03YDyiOcVwJAo993Sz3ZrvpGZjQPGAfTs2ZOSkpIoX15ERADMbN2h1qXEh7HuPsXdi9y9qHPnFt+QRETkKEUT9BuAHhHPu4eXReNYflZERGIgmqCfB/Q1s95mlg2MBIqjfP1XgavNrKOZdQSuDi8TEZEEOWzQu3sDMIFQQC8HZrn7UjN7yMyuBzCz88ysArgR+J2ZLQ3/7E7gB4TeLOYBDx34YFZERBLjsNMrE62oqMj1YayIyJExs/nuXtTSupT4MFZEROJHQS8iEnDRzKMXEZEoNTU5e2rrqdzX8p8BJx7PZad1SeiYFPQiIkegqcnZVrWfil01VOzaR/nO0N8Vu/ZRvquGjbv3Ud946M8+u3Vox9sTL0/giBX0IiL/wT0U5mu317BmexVrttewdns1a3eE/tTWf/huL4X52XTrmMvHurVn2Jkn0Dk/h/bt2tAhN5v27doc/PO7v73PtLnrcXfMLGHHo6AXkYxV19DE2h3VvL+1irKtVZRtq+L9bVWs2VZNdV3jwe3aZBk9OuXSuyCPC/sUcnJBLj065tK9Yzu6d8ylXXZWVPs7sX1b9jc0UV3XSH5O4uJXQS8igVe9v4H3t1WxaksVq8KhvnpbFet21tDY9EHN0q1DO07tkk/RyZ3oXZhHr8I8ehfkcVKHtrTOOva5KwV5OQDsqNqvoBcRORp7autDZ+Zbqli1dS+rtobCfcPufQe3aZNl9CrI47QTjuO6s06kT5d8Tu2czymd88jNjm8kFuRnA7C9qo6TC/Liuq9ICnoRSWtNTc5j/1jDY/9Yw+Y9tQeX57RuRZ8u+ZzXqyOjuvSgT5fj6NMln14FuTE5Oz8ahfkfnNEnkoJeRNLWruo67n56IW+s2MrFfQsZc0Ev+nbJp2/XfLp3zCWrVeI+8IzGgTP6HdV1Cd2vgl5E0tL8dbu4fVop26vq+P71Z/CF809O6EyWo9EpLxz0OqMXETk0d+f3f1/Dw6+s4MQObZn9tfM5q3uHZA8rKjmtszgupzXbq3RGLyLSot01oarm9eVbufbME/jRDWfRvl2bZA/riBTkZ7NT1Y2IyH8qXb+L26ctYOve2rSpalpSkJ/DjmpVNyIiBzWvap752gVpU9W0pCAvm/U7axK6TwW9iKSsyKpm2Bkn8PCI9KtqmivIz6F0/e6E7lNBLyIpKbKqefBTAxhzQa+0rGqaK8zPZmf1fpqanFYJmv6poBeRlOIeugDqRy+HZ9V89QLO7pG+VU1zBXnZNDns3ld/cLplvCnoRSRlhKqaRby+fAvXnNGVH484O+2rmuY6RVwdq6AXkYyyYP0uJoSrmgc+NYAvBqSqaa4w74P73fTtmph9KuhFJKkiq5oT2gevqmmuIHxGn8i59Ap6EUmaypp6vvn0Ql5fvoWrB3TlJyPOpn1usKqa5j64303i5tIr6EUkKSKrmvs/OYAvXRjMqqa5jrnZmJHQ2yAo6EUkodydx99ey49eXk6X49ry9Fcv4JwAVzXNZbUyOuVmJ/TGZgp6EUmYypp67p69kNeWbeGqAV35aQZUNS0pyM9mh87oRSRo3i3fzfippWzZU8v3PjmAL2dIVdOSgrzE3u9GQS8icfWfVc35DOzZMdnDSqpO+dks37gnYftT0ItI3FTW1HPP7IX8edkWrjy9Kz+98Sw65CbmIqFUVpiXzXZ19CKS7haW72b8tFI2V9by3etOZ+xFvTO2qmmuID+HPbUN1DU0kd06/t9fq6AXkZhyd/74z7X8cI6qmkM5MJd+V00dXY9vG/f9KehFJGYq99XzrdkLeXWpqpqPUpAXujp2e9V+Bb2IpI+F5buZML2UTbtV1RxO4YGrYxM0xVJBLyLHxN154p9r+e9wVTPrq+czSFXNRzpwv5tETbGM6lMAMxtmZivNrMzMJrawPsfMZobXzzWzXuHlbczsCTNbbGbLzeze2A5fRJKpcl89X3uqlAdfWMal/Trz0h0XKeSjcOD2xClzRm9mWcBk4CqgAphnZsXuvixis7HALnfvY2YjgYeBm4AbgRx3/5iZ5QLLzGy6u6+N9YGISGItqgjNqlFVc+SOb9uaNlmWsPvdRHNGPxgoc/fV7l4HzACGN9tmOPBE+PFs4AoL/R93IM/MWgPtgDogcVcJiEjMuTt/fHsNN/zmnzQ2OjNvPZ+vXHyKQv4ImFno6tgEzaWPpqPvBpRHPK8AhhxqG3dvMLNKoIBQ6A8HNgG5wNfdfWfzHZjZOGAcQM+ePY/wEEQkUfbU1vPt2Yt4eclmrujfhZ999mzNqjlKBfnZCbsnfbw/jB0MNAInAR2Bv5vZ6+6+OnIjd58CTAEoKiryOI9JRI7C4opKxk8rZePufXznE6fzlYtV1RyLgvwctico6KOpbjYAPSKedw8va3GbcE3THtgBjAZecfd6d98KvA0UHeugRSRxDsyqueE3/6ShsYmZt57PLZeoqjlWhXmJu1VxNEE/D+hrZr3NLBsYCRQ326YYGBN+PAJ4090dWA9cDmBmecBQYEUsBi4i8bentp7bppbyQPFSLupbyEt3XMy5J2tWTSwk8lbFh61uwp37BOBVIAt43N2XmtlDQIm7FwOPAU+aWRmwk9CbAYRm6/zBzJYCBvzB3RfF40BEJLYOVDUbdu/jvk/05ysXnUKrVjqLj5VOeTnsq2+kpq6B3Oz4tuhRvbq7zwHmNFt2f8TjWkJTKZv/XFVLy0Ukdbk7T/57Hf/14nIK87OZdetQzj25U7KHFTgFEVfH5nZKgaAXkcywp7aeic8sYs7izVzevws/u/FsOuZpVk08HLgNwvaq/fTolBvXfSnoRQSAJRsquW1qqKq599r+3HKxqpp4OnBjs0T09Ap6kQzn7jz173X84MXlFKiqSZgD1U0i5tIr6EUy2J7aeu59ZjEvLd7Ex0/rzM8+e87B+7BIfB28VXECbmymoBfJUEs2hGbVVOzax8Rr+zNOVU1CtcvOIi87S9WNiMReZFXTKS+bmeOGUtRLVU0ydMpPzEVTCnqRDLK3tp6Jzy7mpUWbuOy0zvxcVU1SFeTlsEMdvYjESmRV8+1h/bn1ElU1yVaYn82G3bVx34+CXiTg3J2n5q7nBy8so1NeNjPGDeU8VTUpoSAvh0UVlXHfj4JeJMD21tZz77OLeXHRJi7t15lHblJVk0oO3KrY3eN6kzgFvUhALd1YyfippZTv2se3hp3GVy85VVVNiinIz6Ghydmzr4H2uW3ith8FvUjAuDvT3lnP919YRqfcbKbfMpTBvVXVpKKDt0Go3q+gF5Ho7K2t577nlvDCwo1c2q8zP//s2RTk5yR7WHIIkbdBOLVz/PajoBcJiKUbK5kwbQHrd9aoqkkTBz4vifdcegW9SJqLrGo65rZRVZNGPqhu4juXXkEvksaq9jdw37OLKV64kUv6deYRVTVppaPO6EXkoyzbuIfx00pZt6Oae645ja9dqqom3bTJakWH3DZxv9+Ngl4kzbg7098p58EXlh6saoacUpDsYclRKsjLZkec72CpoBdJI5FVzcV9C3nkpnMoVFWT1gryc3RGLyIhyzftYfzUUtaqqgmUwvxs3ttSFdd9KOhFUpy7M2NeOQ8WL6V9O1U1QdMpL/63KlbQi6Swqv0NfOe5xTz/rqqaoCrIy2FXTT0NjU20zmoVl30o6EVSVGRVc/fV/bjtsj6qagLowFz6nTV1dDmubVz2oaAXSTHuzsx55TwQrmqm3TKUoapqAuvAdQ87qhT0IhmhOlzV/ElVTcYoOHjRVPxm3ijoRVLE8k2hC6DWbq/mm1f147aP9yFLVU3gHTyjj+NcegW9SJJFVjXHt2vD1K8M5fxTVdVkigMdvc7oRQKqen8D3/3TEp5bsIGL+oSqms7HqarJJMe3bUNWK9MZvUgQrdi8h9umhqqab1zVj/GqajJSq1YWnkuvM3qRwHB3ZpWUc//zoarmqa8M4YJTC5M9LEmigrxstivoRYIhsqq5sE8Bv7hpoKoaoTA/R9WNSBCs3LyX26bOZ/X2ar5+ZT8mXK6qRkIK8rNZv74mbq8f1fW2ZjbMzFaaWZmZTWxhfY6ZzQyvn2tmvSLWnWVm/zKzpWa22Mzic0WASIpyd2bNK2f45H9Qua+BqWOHcOeVfRXyclBBXk5c73dz2DN6M8sCJgNXARXAPDMrdvdlEZuNBXa5ex8zGwk8DNxkZq2Bp4DPu/tCMysA6mN+FCIpqnp/A9/70xKeDVc1j9x0TtyufpT0VZCfTXVdI7X1jbRtkxXz14+muhkMlLn7agAzmwEMByKDfjjwYPjxbGCSmRlwNbDI3RcCuPuOGI1bJOWpqpFoHbw6trqObh3axfz1o6luugHlEc8rwsta3MbdG4BKoADoB7iZvWpmpWb2rZZ2YGbjzKzEzEq2bdt2pMcgknJmlaiqkeh9cL+b+NQ38f4wtjVwEXAeUAO8YWbz3f2NyI3cfQowBaCoqMjjPCaRuKmpa+B7f1rKM6UVXHBqAb8YqapGDq8gzlfHRhP0G4AeEc+7h5e1tE1FuJdvD+wgdPb/N3ffDmBmc4BBwBuIBMx7W/Zy29RS3t9WxV1X9uX2y3UWL9EpzAud0W+P0xl9NNXNPKCvmfU2s2xgJFDcbJtiYEz48QjgTXd34FXgY2aWG34DuJQPd/sigfB0STnXT/oHu2vqeWrsEO66sp9CXqJ28Iy+Okln9O7eYGYTCIV2FvC4uy81s4eAEncvBh4DnjSzMmAnoTcD3H2Xmf2c0JuFA3Pc/aW4HIlIEkRWNeefUsAvR6mqkSOXm51F2zatktvRu/scYE6zZfdHPK4FbjzEzz5FaIqlSKCsClc1ZduquPOKvtxxhaoaOTpmFp5Ln7yOXkSamT2/gu/9aQl5OVk8NXYIF/bRvWrk2Iw8rwcntNc3TIkkXU1dA/c/v5TZ88NVzchz6HK8qho5drdf0Tdur62gF4lSZFVzxxV9uVNVjaQJBb1IFJ6ZX8F3w1XNk18ewkV9VdVI+lDQi3yEfXWN3P/8Ep6eX8HQUzrxq5EDVdVI2lHQixzCqi17GT+tlFVbVdVIelPQi7TgQFWTm53F/315MBf37ZzsIYkcNQW9SIR9dY08ULyEWSUVDOndiV+NGkhXVTWS5hT0ImFlW0OzalZtreL2y/tw5xV9aZ0V1XfziKQ0Bb0I8GxpBd95LlTVPPGlwVzST1WNBIeCXjJaZFUzuHcnHlVVIwGkoJeMVba1ivFTS3lv614mfLwPd12pqkaCSUEvGem5BaGqpm2bLP74pcFcqqpGAkxBLxllX10jDxYvZWZJOYN7hy6AiteNpERShYJeMsaBqmblFlU1klkU9JIR/rRgA/c9t5i2bbJ44suqaiSzKOgl0GrrQ1XNjHnlDO4VugBKVY1kGgW9BNb720JVzYrNexn/8VP5+pX9VNVIRlLQSyA9/+4G7n12cXhWzXlcdlqXZA9JJGkU9BIotfWNfP+FpUx/p5zzenXkV6MGcmL7dskelkhSKeglMCKrmtsuO5VvXKWqRgQU9BIQz7+7gfueXUx261aqakSaUdBLWgtVNcuY/s56VTUih6Cgl7S1elsVt4Wrmq9ddirfVFUj0iIFvaSlyKrmD186j4+rqhE5JAW9pJXa+kYeenEZ0+aup+jkjjw6WlWNyOEo6CVtrN5WxfhpC1i+aQ9fC8+qaaOqRuSwFPSSFooXbuTeZxaFqpovnsfH+6uqEYmWgl5SWm19Iz94cRlT567n3JM78uiogZzUQVWNyJFQ0EvKWrO9mtumlrJ80x5uvfQU7r76NFU1IkdBQS8p6YWFG5n4zCLatG7F418s4vL+XZM9JJG0paCXlBJZ1Qzq2YFJowepqhE5RlH9Hmxmw8xspZmVmdnEFtbnmNnM8Pq5Ztar2fqeZlZlZnfHZtgSRGu2V/OZX/+TqXPXc+slpzDz1vMV8iIxcNgzejPLAiYDVwEVwDwzK3b3ZRGbjQV2uXsfMxsJPAzcFLH+58DLsRu2BM0LCzdy77OLaZ1lPDamiCtOV1UjEivRVDeDgTJ3Xw1gZjOA4UBk0A8HHgw/ng1MMjNzdzezTwNrgOqYjVoCo7a+kf96aRlP/TtU1Tw6ehDddBYvElPRBH03oDzieQUw5FDbuHuDmVUCBWZWC3yb0G8Dh6xtzGwcMA6gZ8+eUQ9e0tva8KyaZZv2cOslp3D3NZpVIxIP8f4w9kHgEXevMrNDbuTuU4ApAEVFRR7nMUkKeHHRRiY+s5isVsbvv1DElQNU1YjESzRBvwHoEfG8e3hZS9tUmFlroD2wg9CZ/wgz+zHQAWgys1p3n3TMI5e0VFvfyH+/tJwn/72OgeFZNapqROIrmqCfB/Q1s96EAn0kMLrZNsXAGOBfwAjgTXd34OIDG5jZg0CVQj5zrd1ezfhppSzduIdxl5zCPapqRBLisEEf7twnAK8CWcDj7r7UzB4CSty9GHgMeNLMyoCdhN4MRA56adEmvv3MIlU1IklgoRPv1FFUVOQlJSXJHobESG19Iz+cs5z/+9c6zunRgUmjB9K9Y26yhyUSOGY2392LWlqnK2MlbtbtCFU1Szbs4ZaLe3PPNf3Jbq2qRiTRFPQSF3MWb+LbsxfRqpXxv18o4ipVNSJJo6CXmNrf0MgPX1rOE6pqRFKGgl5iZt2OaiZMW8DiDZV85aLefGuYqhqRVKCgl5h4efEmvjV7EWYw5fPncvUZJyR7SCISpqCXY7K/oZH/mbOCP/5zLWf36MCkUQPp0UlVjUgqUdDLUVu/o4bx00pV1YikOAW9HBVVNSLpQ0EvR0RVjUj6UdBL1NbvqGHC9FIWVVQy9qLefFtVjUhaUNBLVF5Zsol7Zi8C4HefP5drVNWIpA0FvXykD1U13dszafQgVTUiaUZBL4dUvrOGCdNKWVhRyZcu7MW9156uqkYkDSnopUWvLNnMPbMXAvDbz53LsDNV1YikKwW9fEhdQxP/8/Jy/vD2Ws7q3p7JqmpE0p6CXg5qXtVMvLY/Oa2zkj0sETlGCnoB4NWlm7nn6YU48NvPDWLYmScme0giEiMK+gxX19DEj15eweNvr+Gs7u2ZNGoQPQtU1YgEiYI+g0VWNV+8oBf3fkJVjUgQKegzlKoakcyhoM8wdQ1NPPzKCh77xxo+1i00q0ZVjUiwKegzSPnOGiZMX8DC8t2qakQyiII+Q/x56Wbufnoh7vCbmwdx7cdU1YhkCgV9wDWvaiaNHsjJBXnJHpaIJJCCPsAqdtUwYdoC3i3fzZjzT+a+605XVSOSgRT0AfXasi18c9a7uMOvbx7EJ1TViGQsBX3A1Dc28eNXVvC/f1/Dmd2OZ/LoQapqRDKcgj5AKnbVcPv0BSxYr6pGRD6goA+I15dt4ZtPL6SxyZk8ehDXnaWqRkRCFPRprr6xiZ+8upIpf1vNGSeFqppehapqROQDCvo0tmH3PiZMK2XB+t184fyTue8Tp9O2jaoaEfkwBX2aemP5Fr4xS1WNiBxeVF8AambDzGylmZWZ2cQW1ueY2czw+rlm1iu8/Cozm29mi8N/Xx7b4Wee+sYmfjhnOWOfKKFbh3a8ePtFCnkR+UiHPaM3syxgMnAVUAHMM7Nid18WsdlYYJe79zGzkcDDwE3AduBT7r7RzM4EXgW6xfogMsWG3fu4fVoppet387mhPfnudQNU1YjIYUVT3QwGytx9NYCZzQCGA5FBPxx4MPx4NjDJzMzdF0RssxRoZ2Y57r7/mEeeYd5YHppV09DoPDpqIJ86+6RkD0lE0kQ0Qd8NKI94XgEMOdQ27t5gZpVAAaEz+gNuAEoV8kemvrGJn766kt/9bTUDTjyeyTcPordm1YjIEUjIh7FmdgahOufqQ6wfB4wD6NmzZyKGlBY2hmfVqKoRkWMRTdBvAHpEPO8eXtbSNhVm1hpoD+wAMLPuwHPAF9z9/ZZ24O5TgCkARUVFfiQHEFRvrgjNqqlvaOJXowZyvaoaETlK0QT9PKCvmfUmFOgjgdHNtikGxgD/AkYAb7q7m1kH4CVgoru/HbthB1d9YxM//fNKfvfX1Zx+4vH8WlWNiByjwwZ9uHOfQGjGTBbwuLsvNbOHgBJ3LwYeA540szJgJ6E3A4AJQB/gfjO7P7zsanffGusDCYKNu/dx+/QFzF+3i5uH9OR7n1RVIyLHztxTqykpKirykpKSZA8j4d5asZWvz3qX+oYm/ueGs1TViMgRMbP57l7U0jpdGZtk9Y1N/OzP7/Hbv77P6Scez+TRAzmlc36yhyUiAaKgT6KNu/dxx/QFlKzbxeghPblfVY2IxIGCPkneWrGVb8x6l7qGJn458hyGn6MLhkUkPhT0CRZZ1fQ/4Th+ffMgVTUiElcK+gTaVLmP26eFqppRg3vywKdU1YhI/CnoE+StlVv5xkxVNSKSeAr6OGtobOJnr73Hb/4Sqmom3zyIU1XViEgCKejjaFNlaFbNvLWqakQkeRT0cfKXlVv5xqyF1NY38oubzuHTA1XViEhyKOhjrKGxiZ+/9h6/VlUjIilCQR9DmytruWP6At5Zu5NRg3vwwKfOUFUjIkmnoI+Rv763ja/PfFdVjYikHAX9MWpobOKR199j8lvvc1rXUFXTp4uqGhFJHQr6Y7BlTy23T1/AO2t2MvK8UFXTLltVjYikFgX9UfpbuKqpqWvkkZvO5v8N7J7sIYmItEhBf4QaGpv4xeurmPyXMvp1OY7JNw+kT5fjkj0sEZFDUtAfgciq5qaiHjx4vaoaEUl9CvooRVY1P7vxbG44V1WNiKQHBf1hNDQ28cs3VjHprTL6dsln5s2DVNWISFpR0H+ELXtCF0DNXbOTzxZ15/vXn6mqRkTSjoL+EP6+aht3zVBVIyLpT0HfTGOT88vX3+PRt8ro0zmfGeMG0berqhoRSV8K+ghb99Ryx4wF/Hv1Tm48tzvfH34Gudn6TyQi6U0pFvb3VaFZNdX7G/npjWczQlWNiARExgd986pm+i2qakQkWDI66COrmhHnduchVTUiEkAZm2r/WLWdu2YuoGp/Az8ZcRY3FvVI9pBEROIi44K+scn55RurePTNVfTpnM+0W4bST1WNiARYRgX91j213DnjXf61egc3DOrODz6tqkZEgi9jUu7tsu3cOeNdqvbXq6oRkYwS+KBvbHJ+9cYqfvXmKk7tnM+0W4aoqhGRjBLooN+6t5Y7p4eqms8M6sZ/ffpMVTUiknECm3qqakREQlpFs5GZDTOzlWZWZmYTW1ifY2Yzw+vnmlmviHX3hpevNLNrYjf0ljU2Ob94/T0+99hc2rdrzfPjL1LIi0hGO+wZvZllAZOBq4AKYJ6ZFbv7sojNxgK73L2PmY0EHgZuMrMBwEjgDOAk4HUz6+fujbE+EIBte/dz18wFvF22g88M7MYPPn0meTmB/aVFRCQq0aTgYKDM3VcDmNkMYDgQGfTDgQfDj2cDk8zMwstnuPt+YI2ZlYVf71+xGf4HFldU8uUn5rG3tp4fjziLG8/tTmgIIiKZLZqg7waURzyvAIYcaht3bzCzSqAgvPzfzX62W/MdmNk4YBxAz549ox37h3Tv2I7+JxzHd68bwGknaFaNiMgBUXX08ebuU9y9yN2LOnfufFSv0TEvmyfHDlHIi4g0E03QbwAiP83sHl7W4jZm1hpoD+yI8mdFRCSOogn6eUBfM+ttZtmEPlwtbrZNMTAm/HgE8Ka7e3j5yPCsnN5AX+Cd2AxdRESicdiOPty5TwBeBbKAx919qZk9BJS4ezHwGPBk+MPWnYTeDAhvN4vQB7cNwPh4zbgREZGWWejEO3UUFRV5SUlJsochIpJWzGy+uxe1tC4lPowVEZH4UdCLiAScgl5EJOAU9CIiAZdyH8aa2TZg3TG8RCGwPUbDSVU6xuDIhOPUMSbGye7e4hWnKRf0x8rMSg71yXNQ6BiDIxOOU8eYfKpuREQCTkEvIhJwQQz6KckeQALoGIMjE45Tx5hkgevoRUTkw4J4Ri8iIhEU9CIiAReYoD/cF5inOzPrYWZvmdkyM1tqZncme0zxZGZZZrbAzF5M9ljiwcw6mNlsM1thZsvN7PxkjynWzOzr4X+rS8xsupm1TfaYYsHMHjezrWa2JGJZJzN7zcxWhf/umMwxNheIoI/4AvNrgQHAqPAXkwdJA/BNdx8ADAXGB/AYI90JLE/2IOLol8Ar7t4fOJuAHauZdQPuAIrc/UxCtzgfmdxRxcwfgWHNlk0E3nD3vsAb4ecpIxBBT8QXmLt7HXDgC8wDw903uXtp+PFeQsHwH9+/GwRm1h24Dvh9sscSD2bWHriE0Pc44O517r47uaOKi9ZAu/C3zuUCG5M8nphw978R+t6NSMOBJ8KPnwA+ndBBHUZQgr6lLzAPZAgCmFkvYCAwN7kjiZtfAN8CmpI9kDjpDWwD/hCup35vZnnJHlQsufsG4KfAemATUOnuf07uqOKqq7tvCj/eDHRN5mCaC0rQZwwzyweeAe5y9z3JHk+smdknga3uPj/ZY4mj1sAg4DfuPhCoJsV+1T9W4Y56OKE3tdgRcKEAAAE/SURBVJOAPDP7XHJHlRjhr1FNqXnrQQn6jPgScjNrQyjkp7r7s8keT5xcCFxvZmsJVXCXm9lTyR1SzFUAFe5+4Dey2YSCP0iuBNa4+zZ3rweeBS5I8pjiaYuZnQgQ/ntrksfzIUEJ+mi+wDytmZkR6nSXu/vPkz2eeHH3e929u7v3IvT/8U13D9SZoLtvBsrN7LTwoisIfa9ykKwHhppZbvjf7hUE7APnZoqBMeHHY4DnkziW/3DYLwdPB4f6AvMkDyvWLgQ+Dyw2s3fDy+5z9zlJHJMcvduBqeETk9XAl5I8nphy97lmNhsoJTRjbAEpfpuAaJnZdOAyoNDMKoAHgB8Bs8xsLKHbrH82eSP8T7oFgohIwAWluhERkUNQ0IuIBJyCXkQk4BT0IiIBp6AXEQk4Bb2ISMAp6EVEAu7/A3N+T3mB6JUSAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(cubics, water_heights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3yV1f3H3yd7hySEBLIDAdkrBELYULe4QBkq4ACtq+vX2vHTqvXXVltbtbZCVURlioqiqFUQwgwr7BmyQwjZe9/z++O5sSG5N1ySO5Pzfr3Oi+R5zpPzaRPveb7nu4SUEoVCoVAorhUnWwtQKBQKhWOiNhCFQqFQdAq1gSgUCoWiU6gNRKFQKBSdQm0gCoVCoegULrYWYE169+4to6OjbS1DoVAoHIpDhw4VSSmD217vURtIdHQ0Bw8etLUMhUKhcCiEEFmGrqsjLIVCoVB0CrWBKBQKhaJTqA1EoVAoFJ1CbSAKhUKh6BRqA1EoFApFp1AbiEKhUCg6hdEwXiFEoAnP66SUZWbUo1AoFAoHoaM8kIv6ITqY4wxEmlWRQqGwOTqdZMPBHG4YGkqAt5ut5SjslI42kNNSytEdPSyESDWzHoVCYQd8d7qAZz45zplLlfx+9lBby1HYKR35QBJNeN6UOQqFwsFYkZwOwMZDuVTWNdpYjcJeMbqBSCnrWr4WQgQIIUYIIca0jLZzFApF9+BQVgkHs0q5a0wYVfVNbDyUa2tJCjvlqrWwhBAvAouBC0BL/1sJzLCcLIVCYSuW70jH39OVF28fRlZxDav2ZLIoMRonp47coYqeiClhvPcA/aWU06SU0/VDbR4KRTfkQmEV354u4IHEKLzdXVg8MZrM4hq+P3vZ1tIUdogpG8gJoJelhSgUCtvz9s50XJ2deCAxGoAbh4US6ufBe3sybapLYZ+YsoH8EUgVQnwjhPi8ZVhamEKhsC6FlfV8fDiPu8eEE+zrDoCrsxP3J0ax83wR5wsqbaxQYW+YsoGsAv4M/An4a6uhUCi6Eav2ZNLYrOORyTFXXJ+fEIm7i5OyQhTtMKWhVI2U8nWLK1EoFDajur6JD/Zlcf2QEGKDfa64F+jtxh2jwvjkcB6/vOE6/L1cbaRSYW+YYoHsFEL8UQiR2DaMV9H9+TQ1l19uPIqU8uqTFQ7L+gM5lNc2smxqf4P3FydFU9vYzLoD2VZWprBnTLFAWrLRJ7S6psJ4ewDNOslfvjlHXlktNw/vy7RBfWwtSWEBGpt1vLMrg3HRAYyJDDA4Z3BfPybEBvL+3iwemhSDi7Oqw6owwQJpFbo7XYXx9ix2nLtMXlktbs5OvLb1vLJCuilbjueTV1bL0imGrY8WFk+MIa+slu9OF1hJmcLeMbqBCCFuvdrDpsy5yvM3CiHOCiHShBDPGLj/MyHEKSHEMSHEViFEVKt7zUKII/qhosIswJqUbIJ93fnNzdeRml3G7rRiW0tSmBkpJct3pNM/2JuZ13VsYf5oSAjhAZ6s3J1pHXEKu6cjC+QVIcTo1n6PtgP4v84uLIRwBt4EbgKGAPOFEEPaTEsF4qWUI4CNwMut7tVKKUfpx+zO6lAYJq+slm1nLnNvfATzx0cS6ufBa1vPKSukm7ErrYhT+RUsnRJ71UxzZyfBosRoUjJKOHmx3EoKFfZMRz6QAuDVqzx/vgtrJwBpUsp0ACHEOuB24FTLBCnl963m7wPu68J6imtg/f5sJDAvIQJ3F2cem9af5z4/yd70Yib2721reQozsSI5nWBfd+4YHWbS/HviI3j123Os2pPJy3NGWlidwt7pqJjiNCP+j9bj7i6sHQbktPo+V3/NGA8BX7X63kMIcVAIsU8IcYexh4QQS/XzDhYWFnZBbs+hsVnHugM5TBsYTHiAFwD3jougj687r2/tyjuDwp44kVfOzvNFLEmKxt3F2aRn/L1cuXtsGJuOXKS4qt7CChX2jkOEUggh7gPigVdaXY6SUsYDC4C/CyEMegCllCuklPFSyvjg4GArqHV8tp6+zOXKehaO/8HlhIerM49O7c++9BJS0pUvpDvw753peLs5X/F7NoXFE6NpaNKxdr8K6e3p2HIDyQMiWn0frr92BUKIWcBvgdlSyh9eeaSUefp/04Ht/DfcWNFFVqdk0dffg2mDrtxw5ydE0tvHnTe2pdlImcJc5JbW8MWxfOYnROLveW2JgQP6+DI5rjcf7MuisVlnIYUKR8CWG8gBIE4IESOEcAPmAVdEUwkhRgPL0TaPy62uBwgh3PVf9waSaOU7UXSerOJqdp4vYt64yHax/p5uziybEsuutCIOZZXYSKHCHLyzKwMBPDgp5qpzDfFgUgwFFfV8deKSeYUpHIqrbiBCiLsMjJlCiC5llUkpm4AngG+A08AGKeVJIcQLQoiWqKpXAB/gozbhuoOBg0KIo8D3wJ+klGoDMQNr9+fg7CS4d1yEwfsLJ0QS5O3Ga1uVFeKolNc0sv5ADreN7Ee/Xp6d+hlTBwYT09ub93ZnmFmdwpEwJRP9IbTWtS0RUdOAQ0CMEOIFKeUHnV1cSrkF2NLm2rOtvp5l5Lk9wPDOrqswTH1TMx8dzGHW4D6E+nsYnOPl5sIjU2L501dnSM0uZbSRzGWF/fJhShY1Dc0snRJrfJJOB4dWwpA7wDuo3W0nJ8GixCh+v/kUR3PKGBmhOj70REw5wnIBBksp79ZHXQ1BK2UyHviVJcUprMs3Jwsorm5gwVWcqvdPiCLAy1X5QhyQusZmVu7OZMrAYAb39TM+8dSn8OXPYOdfjE65e2w4Pu4urFRWSI/FlA0kQkrZunbBZf21EqDRMrIUtmD1viwiAj2ZPKDjPA9vdxcenhzLtjOXOZZbZiV1CnPwaWoeRVX1PNqh9dEM2/+sfX1kNTTWGpzm6+HK3Phwvjyez+WKOguoVdg7pmwg24UQXwghFgkhFgGf6a95A+rTo5uQdrmKlIwSFiREmdT7+oHEKPw9XXld+UIcBp1O8u/kdIaF+ZHYv/2x1A+c+ASKzsK4h6GuHE5+anTqosRomnSSD1NUSG9PxJQN5HHgPWCUfrwPPC6lrJZSTregNoUVWZOSjauzYG58uEnzfT1ceWhSDN+dLuBEnipr4Qh8e7qA9KJqlk7pjxBGXhJ0zbDjT9BnKNz0MgTFwcGVRn9mdG9vZgzqw5qULOqbmi2kXGGvmFKNV0opN0opf6ofG6UqiNStqGtsZuOhHG4YGkpvH3eTn1s0MRpfDxf+oXwhDsGK5HTCAzy5eVio8UnH1kNxGkx7BpycIf5ByN0Pl44bfWRJUgxFVQ18cTTfAqoV9oypYbznhRDlQogKIUSlEKLCGuIU1uGLY/lU1DVdc0ayv6crS5Ji+PrkJU7nqz8Je+ZgZgmHskp5uKNeHk0NsP2P0HckDL5NuzZyHji7d2iFJA0IIq6PDyv3ZKhimz0MU46wXkZL5POXUvpJKX2llB2EbygcjTUpWcQGezMhNvCan30wKRofd2WF2DvLk9Pp5eXKPUbyewBI/QDKsmHG/0LLEZdXIAy7C45tgPoqg48JIVicFM2JvAoOZZVaQL3CXjFlAymQUp62uBKFTTh1sYLD2WUsHB9l/Fy8A3p5ubF4YjRbTuRzrqDSAgoVXSXtchXfnirggQlReLkZSf1qrIXkVyBiAgxok34V/yA0VMKJjUbXuHN0GH4eLqpXSA/DlA3koBBivRBifutsdIsrU1iFNfuzcHNx4u4xppXzNsRDk2LwdHVWVoid8vbOdNxdnHhgYrTxSQffhcp8mPG7/1ofLYSPg5BhcOAdMHJE5eXmwvyESL4+eYmLZYbDfhXdD1M2ED+gBrgeuE0/utSJUGEfVNc3sSn1IreO6EsvL7dO/5wAbzceSIxm87GLpF02fMyhsA2XK+v45HAec8aGGw+QqK+Cna9C7DSImdz+vhAwdjFcOgYXDxtd6/7EKKSUfLAvyxzSFQ6AKVFYSwyMB60hTmFZPj96kar6a3eeG+KRyTF4uDjz5vfKCrEn3tudSaNOx8OTO0gcTHkLaoo034cxRtwLrt4dOtPDA7y4fkgoa/dnU9ugQnp7Ah31RP+l/t83hBCvtx3Wk6iwBFJKPtyXxXWhvoyJ7HodoyAfd+5PjOKzI3lkFFWbQaGiq1TVN/HhvixuHBpKTG9vw5Nqy2DP6zDwJgiPN/7DPPxg+Bw48bH2jBGWJEVTVtPIZ0fadWZQdEM6skBaHOcH0Yonth0KB+ZYbjknL1awcHxkp5znhnhkciyuzk7KCrET1u3PpqKuqeOiiXv/oWWbT//N1X9g/BJorNEisoyQEBPI4L5+rNydqUJ6ewAdtbTdrP+yRkq5qvVA84koHJjVKVl4uTmb3AvbFIJ93Vk4PopPU/PIKlZWiC1pbNbx7q4MEqIDjVdMri6Cff+CoXdC3xFX/6H9Rmvj4LtGnelCCJYkRXO2oJK9qnNlt8cUJ/qvTbymcBDKaxvZfDSf20f1w9fj2rrRXY1lU2NxdhL88/sLZv25imvjy2P5XCyvY9nUDqyPXX/TLIppJlgfLcQ/CIWnISfF6JTZI/sR6O2mQnp7AB35QG4SQrwBhLXxf7wHNFlNocLsbErNo7axmQUJXXeetyXEz4MFCZF8fDiXnBJlqNoCKSVv7bjAgD4+TB9kpO9bRT4ceBtGzIPggab/8GF3g7ufZoUYwcPVmQUJkXx3uoDsYvU30J3pyAK5iOb/qONK38fnwA2Wl6awBFJKVqdkMSLcn+Hh/hZZY9nUWJyE4F87lBViC3aeL+LMpUqWTo41Xlk5+RXQNcG0a2zp4+atRWSd3AQ1xtsa3zchCmcheH9v5rX9fIVD0ZEP5Kje3zGgjQ/kEymlqlfgoBzMKuVcQRULx0dabI2+/p7cMy6cjw7mkKeSyqzOiuR0+vi6c/vofoYnlGbC4fdhzAMQEH3tC8Q/CM31WukTI4T6e3DT8L6sP5hDdb06sOiumOIDiRZCbBRCnBJCpLcMiytTWIQ1Kdn4urtw20gjHy5m4rFpAwB4a7uyQqzJibxydqUVsSQpBncXZ8OTdrwMwgmm/E/nFgkZAlFJWma6zni+x5KkaCrrmvjkcG7n1lHYPaZsICuBf6H5Paaj9QP50JKiFJahpLqBL4/nc9eYMOM1kcxEWC9P5oyNYP2BHC6Vq2511mJFcjo+7i4sMGZhFp2Ho2u1ZlF+XXiJGPcwlGVB2ndGp4yJDGBkRC9W7slEp1Mhvd0RUzYQTynlVkBIKbOklL8HbrGsLIUl+PhQLg1Nuqv2PDcXP57WH53eoauwPDklNXx5PJ/5CRH4exqJrvv+/8DFEyb9tGuLDb4NfEJh/4oOpy2ZGE16YTXJ5wu7tp7CLjFlA6kXQjgB54UQTwgh7gR8zLG4EOJGIcRZIUSaEOIZA/d/pj86OyaE2CqEiGp1b5G+T8l5fatdRQdIKVmzP5v4qAAGhfpaZc2IQC/uGhPG2v3Zqme2FXhnVwYCrcGTQS6dgJOfwITHwCe4a4s5u2qJhWnfQbHxF4Sbh/cl2Ned9/Zkdm09hV1iygbyNOAFPAWMBe4HuvyBLYRwBt4EbgKGAPOFEEPaTEsF4qWUI4CNaL1JEEIEAs8B44EE4DkhhJFsKQXA3gvFZBRVs3CC5Zznhnh8+gCadJLlycptZklKqxtYfyCH2aP60a+Xp+FJ378E7v4w8QnzLDp2MTi5aL4QI7i5OHHf+Ci2ny3kQqEqtNndMKWY4gEpZZWUMldfSPEuKeU+M6ydAKRJKdOllA3AOuD2Nmt/L6VsCSTfB7Q07L4B+FZKWaKPCPsWuNEMmrotq1Oy6eXlyk3D+lp13aggb+4YFcbqlCwKK+utunZP4sN9WdQ2NhsvW5J7CM5ugaQnwdNM71q+oTB4Nhz5EBqM53ssGB+Jm7MT7ysrpNvRUSLhZiHE58aGGdYOA3JafZ+rv2aMh4CvrvVZIcRSIcRBIcTBwsKeeQ57ubKOb05eYs6YcDxcjUTmWJDHp/enoUnH2zuVFWIJ6hqbWbU3k2mDgrku1Eiz0G0vglcQjH/UvIsnLNVqaR3/yOiUYF93bhvZj42HcqmoazTv+gqb0lEozl+spuIqCCHuA+KBqdf6rJRyBbACID4+vkeGgnx0MJcmnWS+BXM/OiI22IfZI/vx/t4slk6JJchYXwpFp/j4cC5FVQ3GrY/MXZD+PVz/Erib2f8VOUFrNrX/31peiZHCnEuSovn4cC4bDuR0XFpe4VB0lEi4o6NhhrXzgNYNmsP1165ACDEL+C1aX/b6a3lWAc06ydr92UzsH0T/YLPEPnSKJ2YMoK6pmbd3ZdhMQ3ekWSd5e2cGw8P8SYwNaj9BStj2B/DtC+MeMr8AIbSQ3oLjHdbHGhbmz7joAN7fm0WzCuntNlzVByKEyGidQGjGRMIDQJwQIkYI4QbMQyuT0nrt0cBytM3jcqtb3wDXCyEC9M7z6/XXFG1IPl9Ibmmt8bwAKzGgjy+3jujH+3syKa1usKmW7sS3pwrIKKpm2dRYw2X507ZC9l6Y8gtwNeJc7yoj7tGc81cJ6V08MYbskhq2nbnc4TyF42BKFFY8ME4/JgOvY4ZEQillE/AE2gf/aWCDlPKkEOIFIcRs/bRX0EKGPxJCHGnxvUgpS4AX0TahA8AL+muKNqzel01vHzeuHxJqayk8OWMA1Q3NvLtbWSHmQErJ8uQLRAR6cuNQA79fKTXfR69IGP2A5YS4ecPohXDqM6i8ZHTaDUND6OfvwXt71O+/u2BKFFZxq5Enpfw7ZkoklFJukVIOlFL2l1K+pL/2rJSyZaOYJaUMkVKO0o/ZrZ59V0o5QD+M99nswVwsq2XbmQLuiY/AzcWUdwXLMjDEl5uHh/Le7kzKa5QztasczColNbuMRybH4uJs4Pd7ejPkH4FpvwaXzve8N4lxD2vFGQ+tMjrFxdmJ+xOj2Z1WzNlLlZbVo7AKphxhjWk14oUQj9Kx811hJ6w/kIME5ifY9viqNU/OiKOyvklZIWZg+Y50ArxcmTs2ov1NXbOW99F7oFY919IE9Yf+M+HQSmg2/nIwb1wEHq5OygrpJpjyWvrXVuOPaMmE91hSlKLrNDXrWHcgmylxwUQEetlazg8M7uvHDUNDeHd3hgrp7AJplyv57nQB9ydG4+lmIDT7+EYoPKO1qnWyUuh2wlKozIfTxqP8A7zduHN0GJ+m5ilfWDfAlCOs6a3Gj6SUj0gpz1pDnKLzbD1zmYKKeouWbe8sT86Io7KuiVWqY12n+XdyBu4uTixKNFDXrLkRtv8fhA6Hwbe3v28p4q6HwFitTW4HLJoYTV2jjnUHcjqcp7B/TDnC6iWEeEoI8WrrzoTWEKfoPKtTsgn182DGdUY60tmQYWH+zBrch7d3ZVCprJBr5nJFHZ+m5jE3PtxwTk3qh1rPjxn/C05W9H05OcH4xyD3AOQeNDrtulA/JvYP4oO9mTQ166ynT2F2TPnr2gJEA8e5sjOhwk7JLq5h5/lC5iVEGHau2gFPzoijvLaR9/dm2VqKw7FyTyaNOh0PTzKQkNdYp3UbDE/QLAJrM2q+1vL2KlbI4onRXCyv4z+nCqwkTGEJTPl08ZBS/kxKubJ1Z0KLK1N0mrUHshHAveMMOFfthJERvZg2KJi3d6arjnXXQFV9Ex/uy+KmYaFE9/ZuP+Hgu1CRBzP/12hWuEVx99Uy0k9tgoqLRqfNHBxCRKAnK1UwhUNjygbygRDiESFEXyFEYMuwuDJFp2ho0rHhQA4zB4fQ199CiWNm4qmZcZTWNPLhPmWFmMq6/dlU1jWxdEr/9jfrq2DnXyFmKsRMsb64FhIeAamDA28bneLsJFiUGM2BzFJO5JVbUZzCnJiygTSgJfTt5b/HV8YPOBU25ZuTlyiubrBL53lbxkQGMDmuNyuS06lpUFbI1Whs1vHOrgzGxwQyKqJX+wkpb0FNkeb7sCUB0TDoZji4EhprjU6bGx+Bl5szK1UwhcNiygbyc2CAlDJaShmjH6oamp2yJiWb8ABPpsR1sWGQlXh6ZhzF1Q2sScm2tRS7Z/PRi+SX17FsqoH//GrLYM/rMPBGiBhnfXFtmfBjqC2BYxuMTvH3dGXO2HA2H71IUZUq9e+ImLKBpAHGi/0r7Ia0y1XsTS9mfkIkTk42OP/uBPHRgUzsH8RbO9Kpa2y2tRy7RUrJiuR04vr4MG2ggci6PW9oZdWn/9b64gwRNVELI973L62kihEWTYymoVmnXiAcFFM2kGrgiBBiuQrjtW/W7s/GxUlwT7z9Os8N8fTMOIqq6lm7X32IGCP5fBFnLlWydEps+5eDqkLtg3rondB3hG0EtkUIzQopPA3p241O6x/sw9SBwXywL4uGJhXS62iYsoFsAl4C9qDCeO2WusZmNh7K5YZhoQT7Ola/jfGxQYyPCeStHReUFWKE5TsuEOLnzu2jDPRN2/U3aKqFab+xvrCOGHY3eAdrvpkOWJIUTWFlPV+dyLeSMIW5MCUTfZWhYQ1xCtPZcjyf8tpGh3CeG+LpmXEUVNSz4aDKTm7L8dxy9lwo5sGkmPZFMcvztGinkfMheKBtBBrDxR3iH4JzX0PxBaPTpsQFE9vbm3eVM93hsGU/EIUZWZ2STWxvb8NNhRyAxP5BxEcF8K/tF6hvUlZIa5YnX8DH3cVwR8nkV7SQ2am/sr4wU4h/EJzdIGW50SlOToLFSdEczSkjNbvUiuIUXcVm/UAU5uPMpQoOZZWyYHyk4aZCDoAQgqdnxZFfXsfGQ7m2lmM35JTUsOV4PgvGR+Ln4XrlzdJMSP0Axi6CAAM1sewB3xAYNkcrr1JrfHO4a0w4vu4uKqTXwbBpPxCFeViTko2bixN3jwm3tZQuMWlAb0ZH9uKf319QDlU97+zKwNlJsCQpuv3N5L+AcIbJP7e6rmsi8XForNay5I3g4+7CPeMi2HI8n4KKOiuKU3QF1Q/Ewamub+KTw3ncOrwvAd4WbhpkYYQQPDUzjryyWj45rKyQ0uoG1h/IYfbIsPZVBUoy4OhaGLsY/PrZRJ/JhA7TeoWkLIcm4/keixKjaZZSVSZwIDrTD2QMqh+I3bD56EWq6pts3vPcXEwbGMyIcH/e3J5GYw+v1PrBvixqG5tZOsVA4uBOvfUx6afWF9YZJj4JVQVw/COjUyKDvJh5XQhrUrJVNJ6D0Jl+IEtVPxD7YXVKNoNCfBkbFWBrKWZBCMHTM+PIKallU2qereXYjLrGZlbtyWT6oGAGhfpeebMkHY6shfgl4NfXNgKvldhpEDJcS3jUGX8xeDApmuLqBjYfNV6IUWE/mHKE9X9CiF6tvg8QQvzBsrIUpnAst4zjeeUsnOC4znNDzLiuD0P7+fGP79N6bL+IjYdyKa5uMFw0Mfkv4OzqONYHaImFE5/UuiSmfWd0WmL/IAaF+LJydyaygwx2hX1gyhHWTVLKspZvpJSlwM2Wk6QwlTUp2Xi6OnPHaAPJZQ5Miy8kq7iGz3vgm2izTvL2znRGhvszIbZN4eviC3B0nRYe6xtqG4GdZdhd4Bem1ewyghBaSO+p/AoOZKqQXnvHlA3EWQjxQ2qzEMITMEuqsxDiRiHEWSFEmhDiGQP3pwghDgshmoQQc9rcaxZCHNEP402YuykVdY18duQis0f2ax/e2Q340eAQrgv15R/b0mjW9aw30f+cvERmcQ1Lp/Rvb1m2WB9JT9tGXFdwdoUJj0HmTriYanTaHaPC6OXlqnqFOACmbCCrga1CiIeEEA8B3wJdzkQXQjgDbwI3AUOA+UKIIW2mZQOLgTUGfkStlHKUfszuqh5HY1NqHrWNzSyc0D2c521xctJ8IelF1XxxrOdYIVJKlienExnoxY3D2lgYxRfg2Dotu9vRrI8WxizSOhbuecPoFE83Z+aNi+Sbk5fILVV1XO0ZU5zofwb+AAzWjxellC+bYe0EIE1KmS6lbADWAbe3WTtTSnkM6JkH4UaQUrImJZvhYf6MCDfQF6KbcMPQUAaG+PDGtjR0PcQKOZBZypGcMh6ZHINz26KJO14GZ3eY9BPbiDMHHn5a6PHJTVBqPFz3/sQohBB8oEJ6u0xW2gm+ff+PNDY2mP1nm+JEjwG2Syl/IaX8BZAshIg2w9phQOvCR7n6a6biIYQ4KITYJ4S4w9gkIcRS/byDhYWFndVqVxzOLuXMpcpuE7prDCcnwZMz4ki7XMWWHlJob/mOCwR6uzFnbJuKykXn4fgGGPcQ+Bgo5+5IjH9Uc6p30Dc9rJcnNwwNYd3+HNVsrAvUNzWTu/5/SEr/G+WF5v9vyJQjrI+40gJo1l+zNVFSynhgAfB3IYSBcBWQUq6QUsZLKeODgx2jydLVWL0vGx93F2aPtPMEMjNw8/C+9A/25o2t3d8KOV9QydYzl3kgMQpPN+crbya/olkfSQ5sfbTgH6aVNzn8foflTZYkxVBe28im1J5zhGlu1m9cR1LjHi4OXUbvfuYvd2PKBuKiP2ICQP+1OVKe84DWr1nh+msmIaXM0/+bDmwHRptBk91TWt3AF8fzuXN0GN7u3b8ggLOTFpF1tqCS/5y6ZGs5FmVFcjoerk48kBh95Y2SdC0Bb9xD4NM9XoKY+ORVy5vERwUwLMyP9/ZkqJDeTpB8toBRp16h3LUPA27/tUXWMGUDKRRC/OCkFkLcDhSZYe0DQJwQIkYI4QbMA0yKptLnorjrv+4NJAGnzKDJ7vn4cC4NTbpuf3zVmltH9CO2tzevbU3rth8kBRV1bDqSx9yxEQS2LUmz6+/g5Kp96HYXWsqb7PuX0b7pQggWT4zhXEEVu9OKrSzQsSmqqmfr+jcY4ZSB500vgJuXRdYxZQN5FPiNECJbCJED/ApY1tWFpZRNwBPAN8BpYIOU8qQQ4oWWDUsIMU4IkQvMBZYLIU7qHx8MHBRCHAW+B/4kpez2G0iL83xMZC8G9/WztRyr4ewkeHz6AE7nV/DtqQJby7EIK3dn0qyTPDw55sob5XlwZA2Mvs9xI6+MMflnUF2oVeo1wm0j+9Lbx4339qiQXlORUvLb9ft4rHk1tcEjcBt1r8XWuuoZiJTyAjBBCOGj/77KXItLKbcAW9pce7bV1wfQjrbaPrcHGG4uHcRNAG0AACAASURBVI7C3vRi0ouq+evckbaWYnVuH9WP17ed5/Vt5/nRkJBulXlfWdfI6n1Z3DSsL1FB3lfe3PM6IB0z7+NqRCVBxHjY/boWmeXcPp/J3cWZBQmRvPF9GlnF1e3//1G0Y9WeTK5LX0WoawncugacTLETOodJP1kIcQvwY+BnQohnhRDPXu0ZhflZnZKNv6crt4xwkPpHZsTF2YnHpw/gRF4F35+9bGs5ZmXd/hwq65vaF02sugyHVsGIe+2330dXEAIm/QzKs+H4RqPT7psQhYuTYNUeFdJ7NU7nV/DOV3v4sduXyCG3Q1SiRdczJYz3LeBe4ElAoB0ndcO/ZvumsLKe/5y8xJyx4Xi4Ol/9gW7InaPDCA/w7Fa+kIYmHe/uzmBCbCAjI9rk9Ox9E5rqtA/Z7srAG6DPUK2vu5Eii338PLhleF8+OphDVb0K6TVGbUMzT61N5VeuG3Bz0iFmPW/xNU2xQCZKKR8ASqWUzwOJgJ01X+7+fHQoh8ZmyfyEnuM8b4ur3go5mlPGjnPdI6dn89GL5JfXsaxt0cSaEq3X+dA7ofcA24izBkJovpCis3D2S6PTFifFUFnfxMaDOUbn9HRe2nIK98Jj3Cq3I8Y/CoExV3+oi5iygbSESNQIIfoBjUDPO0OxITqd5jyfEBvIgD4+tpZjU+4eE05YL09e23re4a0QKSUrktMZFOLLtEFtwnP3r4CGKvvvNmgOhtwBAdGw81Uw8jsdFdGL0ZG9WLU3q9vnA3WGb05e4sN9WbzZ+2PwCoIpv7DKuqZsIF/oy7m/AhwGMjFcm0phIZLPF5JbWsvC8erk0M3Ficem9Sc1u8zhQzu3nyvkbEElj0yJvTIooL5SC28ddLMW7trdcXbREiQvHob07UanLZ4YTUZRdbexPs3FpfI6fvXxMR7pfZKoylSY/hvw8LfK2qbUwnpRSlkmpfwYzfdxXetIKYXlWZOSTZC3GzcM7WZhnJ1kbnw4oX4evLb1nENbISt2pBPq59G+osCh96CurGdYHy2MWgA+obDrVaNTbh7elxA/d1buybSeLjunWSf56foj6Bob+KXTGgi+DsYsttr61xTfJaWsl1KWW0qMoj355bVsPXOZufERuLlYLhzPkXB3ceaxaf05kFnK3nTHtEKO5ZaxN72YBydFX/l7bW7UrI/oyRAebzuB1sbFHSY+ARnJkHvQ4BRXZyfuGx9F8rlC0i6bLZvAoVmefIG96cW8P+wIrhWZcP1LmkVnJdQnkp2z/kAOzTrJgh7sPDfEveMi6OPrzutbz9taSqdYnpyOr7tL+6CIEx9DRR5MfMo2wmzJ2MXg0UvzhRhhwfhI3FycVGIhcCSnjFf/c457hnoxMn05DJgFcbOsqkFtIHZMU7OOdftzmDIwmMggy5QicFQ8XJ15dGp/9qWXkOJgVkh2cQ1fHc9nwYRIfFs3A5NSS6oLHgxxP7KdQFvh7gvjl2nRWAWGC0sE+bhz+8h+fHwoj/LaRisLtB+q6pt4el0qIX4evOj/BaKhGq63fqdxUxMJw4QQE/UdAqcIIaZYWpgCvj9byKWKOhb2oLpX18L8hEh6+7jzxrY0W0u5Jt7elY6zk+DBpDZhlmlb4fJJSHpKC2/tiYx/FFy9YedfjE5ZnBRNbWMzGw703JDeZz87QU5JDW/d6It76krNeusz2Oo6TEkk/DOwG/gd8D/6YZ0YsR7O6pQsQvzcmXmdg/d/sBCebs4smxLLrrQiDmWV2FqOSZRUN7DhYA53jAojxM/jypt7XgPfflqp856KVyAkPAwnPoHCcwanDO3nT0JMIKv2Zva4dscAnx3J45PDeTw5I47hp/4Cbt5a5JUNMMUCuQMYJKW8WUp5m370uBay1ianpIYd5wq5d1wkLs7qpNEYCydEEuTtxmtbHcMK+WBvFnWNuvZlSy6mag7kCY+Cizm6JTgwE58CV0+tB4oRlkyMJre0lu9Od8/imsbIKanhd5+eID4qgKeic+Dc11rOh3dvm+gx5ZMpHWhf5UxhUdbuz0YA88ZFXHVuT8bLzYVHpsSSfK6Q1GzjzYnsgdqGZlbtzWTGdX2IC/G98ubu18HNVzuK6Ol499Z6n5zYCEWGXwx+NCSEsF6erNzdc5zpTc06nlqXCgL+Pncozt/8GgJitGM/G2HKBlIDHBFCLBdCvN4yLC2sJ9PQpGPDwVxmXBdCv16etpZj99w/IYoAL1e794VsPJxLSXUDy9paH6WZcGoTxC+xWgKY3TPxKa0DoxErxMXZiQcSo9iXXsLp/Aori7MNr289T2p2GS/dOZzw86u18i83/lELgbYRpmwgnwMvAnuAQ62GwkJ8e6qAoqp65Tw3EW93Fx6eHMu2M5c5lltmazkGadZJ3t6ZzsiIXiTEBF55c+8/QTjDhMdsI84e8emjWSHHN0DxBYNT7h0XgYerE+/tzrSuNhuQkl7MP75PY87YcGb3d4Xtf9TCdgfeaFNdpmSirwI2APuklKtahuWl9VxWp2QR1suTKQO7SftSK/BAYhT+nq68bqe+kG9OXiKruIZlbcuW1JZC6gcwfC74df8e99fExKfA2Q12/tXg7V5ebtw1JpxNR/IoqW4wOKc7UF7TyE/WHyEqyJvnZw+Frc9DYw3c+CebR+uZEoV1G3AE+Fr//SghhEmtZxXXTnphFXsuFLNgfCTOTj00lLMT+Hq48tCkGL47XcCJPPsqliClZPmOC0QFebUvR5P6ofZhoKyP9viGwNglcHSd1hfeAIsnRlPfpGPt/mwri7MOUkqe+eQYhZX1vDZvFN5Fx7S/mQmPQe84W8sz6Qjr90ACUAYgpTwCxHb0gKLzrN2fjYuTYG58u0aMiquwaGI0vh4u/MPOfCEpGSUczS3n4cmxV74U6Jq1qrtRSdB3hO0E2jNJT4OTi1ErZGCIL5MG9OaDvVk0NhvuJ+LIrD+Qw1cnLvGLGwYxop8ffPUr8A6GKb+0tTTAtA2k0UD9q+73m7ID6hqb+ehQLtcPDaGPr8fVH1Bcgb+nK0uSYvj65CW7cqyuSE4nyNuNuWPbvBSc+xrKsrXsa4Vh/PpqkWlH12nBBgZYkhTNpYo6vjl5yarSLE3a5Sqe33yKpAFBLJ0cq/mDcvfDj54HDz9bywNM20BOCiEWAM5CiDghxBtoDnWFmfnqRD5lNY2qbHsXeDApGh93+7FCzhVUsu3MZR5IjG7fSTLlLfALh0G32EacozDpJyCcjNbImj6oD1FBXqzsRs70+iatu6CHqxOv3jMKp4YK+PZZCBsLI+bZWt4PmLKBPAkMBeqBtUAF8BNLiuqprN6XTUxvbxJjg2wtxWHp5eXG4onRbDmRz7mCSlvLYUVyOh6uTtyf2OaloOCkljiY8LBVq6c6JH79YMwDcGQ1lLbvi+7kJFiUGM2hrFK7jcK7Vl75+iyn8it4Zc5IrWLB1heguhBu+Ss42U9isSlRWDVSyt9KKcdJKeP1X9eZY3EhxI1CiLNCiDQhxDMG7k8RQhwWQjQJIea0ubdICHFePxaZQ48tOXupkoNZpSxIiMRJOc+7xEOTYvB0dba5FXKpvI7PjuRxb3wEgd5tsstTloOLB4xx+D9d6zDpZ1qoc/LLBm/PiQ/H2825W4T07jhXyNu7MnggMYpZQ0Ig9xAceAcSlkK/0baWdwWmRGENFEKsEEL8RwixrWV0dWEhhDPwJnATMASYL4QY0mZaNrCYNh0QhRCBwHPAeDQH/3NCiICuarIla1KycHN24u625+SKaybA240HEqPZfOyiTftGrNydQbNO8vDkNjEnNSVwbAOMuEer/aS4Ov5hEP8gHFlrMDvdz8OVufERbD52kcuVZnm/tQlFVfX8fMNRBoX48pubB0NzE3zxNPiGwvTf2lpeO0yxhT4CUrmymOL/mGHtBCBNSpkupWwA1gG3t54gpcyUUh6jvdP+BuBbKWWJlLIU+BawbUZNF6hpaOKTw3ncPDy0/ZuqolM8MjkGDxdn3vzeNlZIZV0ja1KyuXl4XyIC25TiP/w+NNXatASFQzL5Z1rW9Y4/Gby9aGI0jc2SNSmOGdIrpeQXHx2loq6R1+eP1nxm+5fDpeNw05/txnHeGlM2kCYp5b+klPullIdahhnWDgNa12PO1V8z67NCiKVCiINCiIOFhfbZS3nz0YtU1jexcIJynpuLIB937k+M4rMjeWQUVVt9/bX7s6msb2LZlP5X3mhuggNvax0HQ4ZaXZdD49NHi1g7vtFgv5CY3t5MHxTMh/uyqW9qtoHArvHenky2ny3kd7cMZlCoL5TnwraXIO56GGyf9WuNbiBCiED9UdFmIcSPhRB9W67przsEUsoVet9NfHCwfWZ2r0nJJq6PD/FRDn0KZ3c8MjkWV2cnq1shDU063t2VSWJsEMPD29S2OrsFynOU9dFZJj6lNZ7a/n8Gby9JiqGoqp4tx/OtLKxrnLpYwR+3nGHW4D7c3/Ii+dWvQOrg5r/YPOPcGB1ZIIeAg8AitCOr1rWwDDctvjbygNalZsP11yz9rF1xPLeco7nlLBwfeWWJC0WXCfZ1Z+H4KD5NzSOr2HpWyOdHL3Kpoo5lUw3k2x74N/hHwKCbrKanW+EVCImPw+nNWgn8NkyO603/YG9W7s5ESsfoFVLb0MxT61Lp5eXKy3NGap8DZ7bAmS9g2q8gwH5PJoxuIFLKGCllrP7ftsMcmegHgDghRIwQwg2Yh1a40RS+Aa4XQgTonefX6685HGv2Z+Hh6sSdY5Tz3BIsm6plf//ze8MF+cyNlJIVyRe4LtSXqW1rmRWlaaG7YxeBk7PhH6C4OhMeA88A+L69FSKEYHFSDMdyyzls5+X9W3jxy1OkXa7i1XtGaT7Q2lL44qfQZygkPmFreR1iShTWLiHES/qQW9+rzTcVKWUT8ATaB/9pYIOU8qQQ4gUhxGz92uOEELnAXGC5EOKk/tkStArBB/TjBf01h6KyrpHPjlxk9sh++HuqliuWIMTPgwUJkXx8OJeckhqLr7f9bCHnCqpY2rZoIsChlVpZjtH3W1xHt8bDXytxcv4/kJ3S7vbdY8Lw9XBxiMTCr09cYk1KNsumxDIpTt8U6qtntJyPO/4Jzvb9uWCKE/1+4CxwN7BH75D+mzkWl1JukVIOlFL2l1K+pL/2rJTyc/3XB6SU4VJKbyllkJRyaKtn35VSDtCPlebQY202HblITUMzC1TmuUVZNjUWJyH41w7LWyHLky/Q19+D20a2qazbWAdH1mhHV76hhh9WmE7CUq0m1LYX293ycnNh3rgIvjpxifzyWhuIM4388lqe+eQYw8P8+fn1g7SLZ7+CY+tg8s+h3yjbCjQBUxIJM9DCZLcCyYAXYP3u7d0MKSWr92UxtJ8fI9s6WhVmpa+/J/eMC+ejgznklVnuA+VoThn70kt4MCkG17ZtiE9vhtoSLZdB0XXcvLUP2cydkL6j3e0HEqORUvLhvvaZ6/ZAs07y0/VHaGjS8dq8Ubi5OGlHV5t/oh1dTTFHpoTlMeUI6wKwCQgB3gGGSSkdNufCXjicXcaZS5UsHB+lnOdW4LFpAwB4a7vlrJAVyen4ergwL8FAG+JDKyEgGmKmWWz9HsfYJeAXBtv+AG0c5hGBXswaHMKalGzqGu0vpPetHRfYl17C72cPJTbYR7vY+ujKxTHywUw5wnodLSN8PvAUsEgI0b/jRxRXY01KNj7uLswepZoIWYOwXp7MGRvB+gM5XCo3f6ZyVnE1X53IZ+H4KHw92pxbXz4DWbu1qrJ2VMfI4XH10N7Uc/fD+W/b3V6SFENpTSOfH7loA3HGSc0u5dVvz3HriL7/rdDsYEdXLZhyhPWalHIuMAsthPf3wDkL6+rWlNU08MWxi9w+qh8+7qqQnrX48bT+6KTkLQv4Qt7emYGzk2BJUnT7m4feAydXGHWf2dft8Yy+T7Pstr0IuisLVkyIDeS6UF/e3Z1hNyG9lXWNPL3uCKF+Hrx053Dt9KGmBDY/7VBHVy2YcoT1VyFECpACjACeBWzfCsuB+fhwHvVNOlW23cpEBHpx15gw1u7P5nKF+ayQ4qp6PjqUw52jw7TKqa1prIWja2HwbeBjn4msDo2zK0x9Bi4dg1ObrrglhLahn7lUyb50+wjSfO6zk+SW1vDavFFa5KWU8PmT2iZy578c5uiqBVPs6b3AbCnlUCnlI/qe6Ib7SyquipSS1SlZjI7sxZB+9lfbprvz+PQBNOkky5PN9yf8/t4s6hp1LJ1iID3qzJdQV6aVI1dYhhH3QJ8hmhXS3HjFrdtHhRHg5cp7ezJsJO6/bErN45PUPJ6aGUd8tL6Yx6H3tITBWc9B35E21dcZTDnC2iilLLCGmJ5ASkYJ6YXVyvqwEVFB3twxKozVKVkUVtZ3+efVNjTz/t5MZg3uw4A+BtKkUj/UMs9jpnZ5LYURnJxh5rNa3/TD719xy8PVmfkJkXx7qsAqeUDGyC6u4XebThAfFcAT07WADgrPwte/htjpMOFxm2nrCsqjZ2VWp2Tj5+HCrSP62lpKj+Xx6f1paNLx9s6uWyEbD+VQWtPI0rZFEwHKciB9O4xaoJznlmbgjRCZCDv+DA1Xlq25P1GLdPzARiG9jc06nlqXihDw93mjcHF2gqZ62PgQuHnBnW857N9HR8UUY6wppCdQVFXP1yfyuXtsePv2pgqrERvsw+yR/Xh/bxbFVZ23Qpp1kn/vzGBURC/GRRsohHl0HSC1DURhWYSAWc9DVQHs+9cVt/r6e3LTsFDW7c+mpqHJ6tJe++48R3LK+ONdwwkP0Jf2/+55KDgOt//ToRNLO9r2NgIIIbZaSUu356ODuTQ2SxaOj7S1lB7PEzMGUNfUzNu7On82/vWJS2SX1PDoVANlS3Q6OPKhVrY9ILprYhWmETkeBt0Mu1/TnNKtWJIUTUVdEx8ftm7N1X3pxby5PY25Y8O5dYQ+ZP/8d7DvTS2bfpBjp9R1tIE4CSF+AwwUQvys7bCWwO6CTidZuz+b8TGBhs/KFVZlQB9fbh3Rj/f3ZFJa3XDNz0spWZ58geggL340xMAbZPYeKM3UwkwV1mPms9BQBTv/esXlMZEBjAj35z0rhvSW1TTw0/VHiA7y5vez9VWYKi7Cp8s0p/+PXrCKDkvS0QYyD2gGXABfA0NxDexKKyK7pIYFyvqwG56cMYDqhmbe3X3tVsi+9BKO5ZbzyBSt2m87UleDm6/dNgLqtvQZDCPnw/5/az4oPUIIFk+M5kJhNTvPF1lchpSSZz4+TlFVPa/PG423u4sWIfbREi20e+4qcPW0uA5L01E597NSyj8DD0opn287rKixW7A6JYtAbzduHOa4553djYEhvtw8PJT3dmdSXtN49QdasSL5AkHebtxtqAx/faWWkzDsLs1JqrAu036t/bv9j1dcvmVEX3r7uPPenkyLS1h3IIevT17iF9cP+m9TsW+fg5x9cPsbEDzQ4hqsgSmu/z1CiFdb2sLqEwtV9b9roKCiju9OX2ZufDjuLsp5bk88OSOOyvqma7JCzl6q5PuzhSyaGG04GOLkJmisUcdXtqJXBCQ8oiVwXj79w2V3F2cWjo9k25nLFm1znHa5iuc3n2TSgN48MlmfG3TqM73fYxkMu9tia1sbUzaQd4FK4B79qAAcsny6rVh/IIdmnWRBgjq+sjcG9/XjhqEhvLs7g4o606yQFcnpeLo6/7f1aFuOrYfA/hA+zoxKFdfE5J+Dmw9svdLPsHBCJK7OglUWskLqm5p5am0qXm4uvHrPSJychNZIbNPjEBYP1//BIuvaClM2kP5SyueklOn68Txgjo6EPYKmZh1r92czOa43UUHetpajMMCTM+KorGtilQkNiPLLa/n8aB73josgwNtA2YnyPMjcpWVHqyrLtsMrUGs6dXYLZO/74XIfXw9uG9GPjw7mUGniC8O18PLXZzmVX8HLd4+gj58HNNTAhge0kitz33O4UiVXw5QNpFYIManlGyFEEmC/XVrsjO1nC8kvr1Ohu3bMsDB/Zg3uw9u7Mq76obJydyY6CQ9NMpImdWIjIGH4XPMLVVwbEx4DnxD47vdXlHtfnBRNdUMzHx3MNety289e5p1dGSxKjGLWkBBtzS9/DpdPwd3/1o7WuhmmbCCPAW8KITKFEFnAP4BllpXVfVidkkUfX3dmDg6xtRRFBzw1M47y2kbe32s8W7mirpE1KdncPLwvEYFGnOPHPoKwsRCkOh7YHDdvmPYMZO+Fc9/8cHlEeC/GRgWwam8mOp15QnoLK+v5xUdHGRTiy69v1vfbO7wKjq6Bqb+CAbPMso69YUotrCNSypFolXiHSylHSymPWV6a45NbWsP2c4XMGxfRvkOdwq4YEd6L6YOCeXtnOtX1hrOV16RkU1XfxDJDRRNBc9gWHIfh91hQqeKaGH2/5o/a+jzo/ttYavHEaLKKa/j+7OUuL6HTSX7x0VEq65p4ff5oLbAi7xBs+R/oPxOm/rLLa9grppRz9xdCvApsA7apKCzTWbc/BwHcq5znDsGTM+MorWk02Aa1oUnHyt0ZJA0IYliYkT//YxtAOGvhuwr7wNkVZv6vdox0bP0Pl28cFkqonwcrTfB7XY2VezLZca6Q390ymEGhvlBdBOsf0EqU3P22Vuyxm6KisCxEY7OOdQdymD6oD2G9HD9hqCcwJjKAyXG9WZGc3q5m0mdH8iioqDdcNBG00iXHN0LsNPDpY3GtimtgyB3aseK2l6BR6wPj6uzE/YlR7Eor4nxBZad/9MmL5fz5qzPMGhzCfROiNCtn44Naa9p7PtCc+d0YFYVlIb49VUBRVT0LJyjrw5F4emYcxdUNrEnJ/uGaTidZkZzOdaG+TInrbfjBnH1Qng0j7rWSUoXJCKGVDanIhf3Lf7g8PyESdxcnVnYypLe2QQvZ7eXlystzRmj10La9CBk74NZXHao1bWexaRSWEOJGIcRZIUSaEOIZA/fdhRDr9fdThBDR+uvRQohaIcQR/XjLHHrMyZqUbMJ6eTJ1oHobdSTiowNJGhDEWzvSqWvUzsy3n7vM+ctVLDNUNLGFYxvA1Quuu8WKahUmEz0J4m7QamTpCy0Gertxx6gwPjmce82VCABe+OIU6UXV/O3eUQR6u8HpzbDrbzB2SY9JIjVlA3mU/0ZhZWKmKCwhhDPwJnATMASYL4QY0mbaQ0CplHIA8Dfgz63uXZBSjtKPR7uqx5xkFFWzK62IeeMiDNdJUtg1T82Io6iqnrX7NSvkrR3p9PP3+G811bY0N2qlSwbdDO4+VlSquCZm/V4rM7Pr1R8uLU6Kpq5Rx7oD2UYfM8TXJ/JZuz+bpVNiSRrQG4rOw6ePaUdlN/356j+gm9DhBiKEcAIGtYrCGmHGKKwEIE1/LNYArANubzPndmCV/uuNwExh9BXQfli7PxtnJ8G947pf3HdPYHxsEONjAnlrxwVS0ovZn1HCg5NijEfSpe+A2lLlPLd3QobAyAWQshzKtA1jcF8/JsQG8v7eLJqadSb9mPzyWn718XFGhPvz8x8NgvoqWH+fliR4z/vg4m7J/xV2RYcbiJRSB/xS/3WFlLLCjGuHATmtvs/VXzM4R0rZBJQDQfp7MUKIVCHEDiHEZGOLCCGWttTxKiwsNJ96I9Q1NvPRwRyuHxKiZaIqHJKnZ8ZRUFHPox8ewtfDhXkdRdKd+lSrvNt/pvUEKjrH9N+AcNIc6nqWJMWQV1bLd6ev3rm7WSf5ybojNDbreG3eaNycBXz2OBSdgzkrwd9Acc1ujClHWN8JIX4hhIgQQgS2DIsr65h8IFJKORr4GbBGCOFnaKKUcoWUMl5KGR8cHGxxYd+cvERpTaMq2+7gJPYPYlx0AKU1jdw3IQofdxfDE5sb4cyXWmMgV/XCYPf4h8H4R7WQ3kvHAZg1OITwAE/eNSGk960dF0jJKOH52UOJ6e2tdT88tQlmPgexPa/vvSkbyL3A40AycEg/Dpph7Tyg9RlPuP6awTlCCBfAHyiWUtZLKYsBpJSHgAuAXdRHXr0vm6ggL5L6G4nWUTgEQgieuek6hof5syQp2vjEDP3x1ZA7rKZN0UUm/RQ8e2nl1QFnJ8GixGj2Z5Rw8mK50cdSs0t59dtz3DayH3PGhmvJgt8+C4Nu0epu9UBMyUSPMTDMEcZ7AIgTQsQIIdzQGlh93mbO58Ai/ddzgG1SSimECNY74RFCxAJxQLoZNHWJcwWV7M8sYUFCpFaFU+HQjI0KZPOTk+jj24FlcXKTVvW1m5aq6JZ49oLJv4ALWyF9OwD3xEfg6erMe0askMq6Rp5al0qonwd/uGMYor5Caw7lGwq3/6PHFs60WX0NvU/jCeAb4DSwQUp5UgjxghCipY3bO0CQECIN7aiqJdR3CnBMCHEEzbn+qJTyyibINmBNSjZuzk7a24mi+9PcCGe+gEE3qeMrRyPhEfCP1CwInQ5/L1fuHhvGZ0cvUlxV3276s5+dJK+0ltfmjcLfwwU+fwrKc2HOu90+WbAjbFqgSUq5RUo5UErZX0r5kv7as1LKz/Vf10kp50opB0gpE6SU6frrH0sph+pDeMdIKTfb8n8HaElFHx/O5cZhoQT59JwojB5NRrI6vnJUXNxhxu8g/yic/ATQ6mM1NOl+CN9u4dPUXD5NzePpmQOJjw6EQyv1fo//hYgEW6i3G1SFPzOx+dhFKuuaVNn2nsSpluMrFX3lkAyfC6HDtaZTTfUM6OPL5LjefLAvi0Z9SG9WcTX/u+kk46IDeHx6fyg8B1//Wou4m9gz/R6tMbqBCCHGdDSsKdIRWJ2SzYA+PiTE9FxztkfR3Ainv4CBN4KrqnXmkDg5waznoSwLDr4LwINJMRRU1PPViUs0Nut4et0RhIC/zxuNCzrY9Kj2+77jX9rzPRwjsYkA/LWDexKYYWYtDsuJvHKO5pTx3G1DjJe6UHQvMndCbQkMVcdXuzfNVgAAFsZJREFUDs2AmVoBzB0vw6gFTB0YTExvb1buzuDspQqO5JTxjwWjtYKoyX/RIq/mrARf1d8HOthApJTTrSnEkVmzPxt3FyfuGq2c5z2G019ota9U9JXjM+t5WDEVdr+G08xnWZQYxe83nyI1u4x74sO1EjaXjsP2P8HQu1TFgVaY0g/ESwjxOyHECv33cUKIWy0vzTGoqm/is9Q8bhvZD38vV1vLUVgDKeHsV9B/hjq+6g70G6X5Q/b+EyouMic+Al93F2J7e/PcbUOhqR4+fRQ8A+CWjg5meh6mHOKtBBqAifrv84A/WEyRg7EpNY/qhmblPO9JXEyFyouq8m53YsbvQDbD9j/i4+7Cp49PZMOjiXi7u8Cuv0PBCZj9eo8O2TWEqf1AXgYaAaSUNYA66AeklKxOyWZIXz9GRfSytRyFtTi7RaunFHeDrZUozEVANIx7GFI/hMtnGNDHl94+7lCSrpWAH3a3lu+juAJTNpAGIYQnmuMcIUR/oH2mTQ/kSE4Zp/MrWDA+UjnPexJntkBkIngHXX2uwnGY/AstLHvr89r3Ump9zZ3d4PqXOn62h2LKBvJ74GsgQgixGtiKvkJvT2d1Sjbebs7cMbptEWFFt6UkAy6fVMdX3RHvIJj0E83CzNqrNYhK+w5m/Bb8+tpanV3SURgvAFLK/wghDgET0I6unpZSFllcmZ1TXtPI5qMXuXtsuPFKrYrux9kt2r+DbratDoVlGP8Y7H8bvvkNVF2GkOEw7hFbq7JbTInC2gxcD2yXUn6hNg+Njw/nUt+kY0FHfSIU3Y8zW6DPEAiMsbUShSVw84Lpv4aLh7Ue6rf8FZzVC6IxTDnC+gswGTglhNgohJgjhOjRleOklKzZn83IiF4MC/O3tRyFtagpgew9yvro7oxcAFGTYOJTEDne1mrsGlOOsHYAO/Tl02cAjwDvAgYbOPUE9meUkHa5ipfnjLC1FIU1OfcNSB1cpzaQbo2zCyz+oseWaL8WTLLN9FFYt6E1lxrDf/uU90hWp2Tj6+HCbSP62VqKwpqc/RJ8+0Lf0bZWorA0avMwiatuIEKIDUACWiTWP4Ad+l7pPZLiqnq+OpHPwvFReLo521qOwlo01cOF77WMZVVET6EATLNA3gHmSymbLS3GEdh4KJfGZqkyz3sa2XuhoQoGquRBhaIFU16ldgK/VrWwQKfTnOcJ0YHEhfjaWo7Cmpz/Vksoi5liayUKhd2gamFdA7svFJFVXMPCCcr66HGc/w9ETwI3b1srUSjsBlUL6xpYvS+bQG83bhwWamspCmtSmsn/t3fnUVJV1x7Hvz/mUaRlkDAICgqNAxGEGPEJCooDoNFEkSSgJJgsTTTG+DB5L6jRxGiehjjkiQbFAREHFPAtsAExalCEiAwCQkgIts0kCCgzvd8f53YoOt1NdXVNXb0/a9Wqqlu36u6zKGr3veecfdjyMXQemOlInMsqXgsrTht37KFgxUau6NmO+nVqUOf52nnw1v2ZjiKzVheE+y7nZzYO57JMPJ3oYzm8FtZZwMhUBpWNpry/noPFxrCaNPPcLKz/vOkjaNer5l7/X10AzTvBMSdkOhLnssoRz0DMrAD4BiFpPAf0MrN5yTi4pEGSVklaI2lMGa/Xl/R89Pp7kjrGvHZbtH2VpJQOjTlYbDy34J/07dyCTi1q0DXw9QtC8lBteP2/oLgGjt7evwf+/mfoMtDnBjhXSrkJRNLpJTfgOKAI+BToEG2rkmhm+8PAhUA+MExSfqndRgHbzKwz8ADw2+i9+cBVQHdgEPBI9Hkp8ebHm/h0+x6urmlDdxc9AfWawkX3QtGHsOzFTEeUfuvehgO7/fKVc2Wo6BJWRWs3GqGsSVX0BtaY2VoASZOBocBHMfsMJZSTB3gReEhh4Y2hwGQz2wv8XdKa6PPmVzGmMs18awH5TXYxML91Kj4+O+3eBsunQo/h0PNaWDQR5twJ3YZA3RpUCm11AdRpEEZgOecOU24CMbP+KT52W2B9zPNPgNKVy/61j5kdkLQdOCba/m6p95a5KIek0cBogA4dKn8GUbx/H2M2/Yx9DfOoWzwUateQH88PJ8OBPdDrmjDz+vy74Kkh8Pb90P/nmY4ufVYXhL4fX/vcuX+T8zUZzGy8mfUys14tW7as9Ptr1a1H3jfu49idy+G1m0PHcq4rLoYF46HdGXDsKWHb8efAyVfA2w/AltWZjS9V9u6Ef7wDf30a3v1fmH4jbP2bD991rhyZLHRfCLSPed4u2lbWPp9IqgM0Az6L873J020wnPOf8OZvoc1p0Oe6lB0qK6yeFdaCPve/Dt9+wa9hTQHM+AmMmJ4bncqbV8HyV8LqcxuXEY1WP6R2fThpUEZCcy7bZTKBvA90kdSJ8ON/FXB1qX2mASMIfRtXAHPNzCRNAyZJuh/4CtAFWJDSaM8ZAxuWhmGtrfKh09kpPVxGvfsIHNU29HfEatoaBtweEsgHz8Dp38lEdFVXfBBWzoD3HoV17wAKa5z3/zm06QEtT4T6R0Gd+mF/n33uXJniqcY7x8zOO9K2yor6NG4AZgG1gQlmtlzSncBCM5tGKOT4dNRJvpWQZIj2m0LocD8AXJ/yYo+1asFlj8Lj58ELI2D0PDg6B0dlbVgWhq0OuB1q1/33108fCUtfhJljoONZkHd8mgOsAjNYMQ3m3g1bVoV/v4G/glO/BU29uoBzlSUr55p+tOpgI+ANoB+HypccBcw0s67pCDCZevXqZQsXLqzah2xZA4/1h+Yd4dpZYQnMXPLq9bDsZfjJcmiUV/Y+n6+HP54FLbrAtTPLTjTZZu2bMHssfPoBtDgxnG10GwK1alBVAecSJGmRmfUqvb2iTvTrgEVA1+i+5PYqYV2QmqlFZ7j8T+Fy1rQf5Van+o4iWDIFThtWfvIAOLo9DH4AChfCvHvSF18idhTBCyPDCLIvt8DQh+GH86H7ZZ48nKuiiobxjgPGSfqRmT2Yxpiy34nnw3n/HeZFtDkVzrox0xElx18eDP0DX//Rkfc9+XJYMxfe+l0YWJA/5MjvSaeDB+D9x8LlqoP7oP8vwhrXNWkOi3MpFs+a6A9K+jrQMXZ/M3sqhXFlv743Q9ESmH17GOp6QlXnVWbYl1tg4YTQH5DXKb73XPw/sHklTL0Omh8XEkk2+GQhzLgpnCV2HgAX3Ve9+mqcqyaOOA9E0tPA74C+wBnR7d+uhdU4Urgc0rIbvHANbP17piOqmvkPh4mDfW+O/z11G8BVk6BhHjw3LPSNZNLubTD9Jnh8QEiI33wShr/oycO5FCm3E/1fO0grgHw70o7VQFI60UvbuhbG9w/DXr9XUD2HfO7eBg+cAl0GhB/dyipaAk9eAo2aw8jXoFm7pIdYoeJiWPxsOBvcvRX6/AD63QYNjkpvHM7lqEQ60UssA3yMY3nyjocr/gSbV4QRTNUxz74zLqz3ffYtib2/zanw3amwa2tIJOk8E/l0MUw4H6bdEMqtj34TBv3Gk4dzaVBRNd7p0YS9FsBHkmZJmlZyS1+I1UDnAXDeL0PxwXfGZTqaytleCO/+MfR9HHty4p/Ttid85xXY9Rk8di78873kxViWXVvhtZ/C+H5hxcBL/wjXzAzJzDmXFhV1ov8ubVHkgrNuCiXP59wRfog7D8h0RPGZ9xuw4jBKqara9YRRBTB5GEy8BAbdA72uTW7Jkz07Qn/N/Idh/5ehrEy/26Dh0ck7hnMuLhUN430znYFUeyWd6ps/hhevDTPVs73zdtPK0HfQ54dhFFUytOoK35sDL40KxSeXT4XB46q+mt+urbDoSfjLH0KfTbchYTJgq25JCds5V3nxjMLaKWlHqdt6SVMlZfkvZJrVawxXPQsIJg+HvV9kOqLymcGs26BeEzj7p8n97EZ5MPylkDiKlsDDfeDVGypfxbe4GNbNh2k/hvvzw9ld214hOV/5tCcP5zIsnmKKvyestzGJUM7kKuAE4K/ABEKZE1cirxN88wl45nJ46XshoWTjjOflU+Fvc+HCe6HxMcn//Fq1oOdI6HIB/Pm+cKbzwTPQvjecdCEc1zfM6m/Y/NB79u2C7etDuZF1f4E1s2FHIdRpCKd+M4yuat09+bE65xISzzDeD83stFLbFptZj7Jey2YpGcZbngWPwf/dAmd8P0xky6bS53t2wENnhOq6338jPQnui81hidyVM0JfUYl6TcPxzWDv9kPbGzQLSab7ZaGcev2mqY/ROVem8obxxnMGskvStwhLykIoq74nelwNx6ymSe/vh9FB8x8KZyVnXp/piA554274YiMMm5S+s6MmLeGcW8NteyFsWBIuae0sguID0T6tQ4XcVvnhVivn1ztzrlqLJ4EMB8YBjxASxrvAtyU1BG5IYWzV38Bfwef/hFm/CJPr8odmOqJQqv29R+GMUWHobSY0axtuJ12YmeM755IinlpYa4HB5bz8dnLDyTG1asE3xsPEIaE/pF4T6FylZVSqZvc2mPqDMCJq4J2Zi8M5lxPKTSCSbjWzeyU9SBmXqszsxymNLFfUbQhXPx+SyOThMPyFzKxmaAYzbg6XrkZV05IrzrmsUtFF5hXR/UIOXw+k5Obi1SgPvvtKmGsx6crUz9Iuy4LxsPzlMOmu7enpP75zLucccRTWv3aUGpnZrhTHk1JpHYVVlp0b4ImL4MvNoUpshz7pOe6aOfDsFXDihXDlM9457ZyrlISLKUo6U9JHwMro+WmSHklBjLmv6bEwYjo0bglPDYVVM1N/zE0rQ7n5VvmhP8aTh3MuSeL5Nfk9cAHwGYCZfQj8RyqDymnN2sKo10PJj8lXh/IcqbJpJUwcDHXqw7DnoH6T1B3LOVfjxPXnqJmVrs99MAWx1ByNW4QzkRP6w/Qbw+3A3uQeY+NH8OTFYQLjyBlhfoVzziVRPAlkfbSkrUmqK+kWDnWwJ0RSnqQCSauj++bl7Dci2me1pBEx2+dJWiVpcXRrVZV4MqJ+U7h6SlgBcNGTMOGCcMaQDCtmhM+rXTcs8NTypOR8rnPOxYgngfwAuB5oCxQCPaLnVTEGmGNmXYA50fPDSMoDxgJ9gN7A2FKJZriZ9Yhum6oYT2bUqg0DxoaO7W3r4NGzQ92oRM9G9u+Ggl/C88PDXI9Rr0OLLsmN2TnnIkdMIGa2xcyGm1lrM2tlZt82s8+qeNyhwMTo8UTg0jL2uQAoMLOtZrYNKAAGVfG42anbYLh+AXS9GObeBX/4KiycEH8iMYNlL8NDvcOCVj2vgWtn+WUr51xKVTSRsMwJhCWqOJGwtZkVRY83AK3L2KctENv38km0rcQTkg4CLwF3lbdmu6TRwGiADh2y+Ae1ScuwHnnPkTD3bpjxE5h9B3S/FLoOhjanhX1K7NsFm1fCxzPhw8nw+TpofQpcOh06+RgH51zqVVTKJHbCxB2Ey0lxkzSbstdSP2zpOzMzSZUtyjjczAolNSUkkO8AT5W1o5mNB8ZDmAdSyeOk3/H9oNM5sHYeLJ4ES6YcGqnVMA/qNAiPdxYR8rvCe877Zahcm42l451zOamiFQlLLjEh6abY5/Ews3LXdJW0UVIbMyuS1AYoqw+jkMPXGmkHzIs+uzC63ylpEqGPpMwEUi1JYYTWCf1h35dQuAg2LA3Va4v3h0tWzdqHBZXa94Gj2mQ6YudcDRRPNV5Iftn2acAI4J7o/tUy9pkF/Dqm4/x84DZJdYCjzWyLpLrAJcDsJMeXPeo1Dpek/LKUcy7LZGpa8j3AQEmrgQHRcyT1kvQ4gJltBX4FvB/d7oy21QdmSVoCLCacqTyW/iY451zNVm4tLEk7OXTm0QgoqYMlQtfFUakPL7kyXgvLOeeqoUqvSGhmvoaoc865cnllPeeccwnxBOKccy4hnkCcc84lxBOIc865hHgCcc45l5C4l7TNBZI2A+sSfHsLYEsSw8lWNaGdNaGN4O3MNZls53Fm1rL0xhqVQKpC0sKyxkHnmprQzprQRvB25ppsbKdfwnLOOZcQTyDOOecS4gkkfuMzHUCa1IR21oQ2grcz12RdO70PxDnnXEL8DMQ551xCPIE455xLiCeQI5A0SNIqSWskjcl0PMkiaYKkTZKWxWzLk1QgaXV037yiz6gOJLWX9IakjyQtl3RjtD2n2iqpgaQFkj6M2nlHtL2TpPei7+/zkuplOtaqklRb0geSZkTPc7GN/5C0VNJiSQujbVn3nfUEUgFJtYGHgQuBfGCYpPzMRpU0TwKDSm0bA8wxsy7AnOh5dXcA+KmZ5QNfA66P/g1zra17gXPN7DSgBzBI0teA3wIPmFlnYBswKoMxJsuNwIqY57nYRoD+ZtYjZu5H1n1nPYFUrDewxszWmtk+YDIwNMMxJYWZ/RnYWmrzUGBi9HgicGlag0oBMysys79Gj3cSfnjakmNtteCL6Gnd6GbAucCL0fZq305J7YCLgcej5yLH2liBrPvOegKpWFtgfczzT6Jtuaq1mRVFjzcArTMZTLJJ6gh8FXiPHGxrdGlnMbAJKAD+BnxuZgeiXXLh+/t74FagOHp+DLnXRgjJ/3VJiySNjrZl3Xe23BUJXc1mZiYpZ8Z4S2oCvATcZGY7wh+uQa601cwOAj0kHQ1MBbpmOKSkknQJsMnMFknql+l4UqyvmRVKagUUSFoZ+2K2fGf9DKRihUD7mOftom25aqOkNgDR/aYMx5MUkuoSksezZvZytDkn2wpgZp8DbwBnAkdLKvlDsbp/f88Chkj6B+Fy8rnAOHKrjQCYWWF0v4nwx0BvsvA76wmkYu8DXaJRHvWAq4BpGY4plaYBI6LHI4BXMxhLUkTXyP8ErDCz+2Neyqm2SmoZnXkgqSEwkNDf8wZwRbRbtW6nmd1mZu3MrCPh/+JcMxtODrURQFJjSU1LHgPnA8vIwu+sz0Q/AkkXEa671gYmmNndGQ4pKSQ9B/QjlIjeCIwFXgGmAB0IZe+/ZWalO9qrFUl9gbeApRy6bv5zQj9IzrRV0qmEjtXahD8Mp5jZnZKOJ/y1ngd8AHzbzPZmLtLkiC5h3WJml+RaG6P2TI2e1gEmmdndko4hy76znkCcc84lxC9hOeecS4gnEOeccwnxBOKccy4hnkCcc84lxBOIc865hHgCca4SJP0iqna7JKqU2ifFxzsYHecrlXjP2VH14WVH3tu5xPkwXufiJOlM4H6gn5ntldQCqGdmn6bwmF+YWZME3tcRmGFmJyc9KOcifgbiXPzaAFtKJqmZ2ZaS5BGt33BvtIbDAkmdo+2Do7UqPpA0W1LraPvt0Zos8yStlfTjeAKQ9IWk+6KzoNmSesd8xpAUtdu5MnkCcS5+rwPtJX0s6RFJ55R6fbuZnQI8RKheAPA28DUz+yphtvStMft3BS4g1DkaG9XsOpLGhBIe3YGdwF2EsiWXAXcm2C7nEuIJxLk4Rett9ARGA5uB5yWNjNnluZj7M6PH7YBZkpYCPwO6x+z/mpntNbMthMJ48ZTn3gfMjB4vBd40s/3R446VbZNzVeEJxLlKMLODZjbPzMYCNwCXx75cxuMHgYeiM5PrgAYx+8TWazpIfMsr7LdDHZfFJZ9hZsVxvt+5pPEE4lycJJ0kqUvMph6EonYlroy5nx89bsah8uIjcC6H+F8szsWvCfBgVDb9ALCGcDmrRHNJSwhnBcOibbcDL0jaBswFOqUvXOdSy4fxOpcE0SJHvaL+jGR+rg/jdVnLL2E5l912JDKREJgOJDWZOVean4E455xLiJ+BOOecS4gnEOeccwnxBOKccy4hnkCcc84lxBOIc865hPw/Lq9xuAp8CGoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ss.nodes_range('x')[:-2], ss.nodes_range('y')[:-2])\n",
    "plt.plot(ss.nodes_range('x')[:-2], [a + b for a, b in zip(ss.nodes_range('y')[:-2], ss.get_node_result_range(\"uy\")[:-2])])\n",
    "\n",
    "plt.ylabel(\"Height level roof when accumulating [m]\")\n",
    "plt.xlabel(\"Span [m]\")\n",
    "plt.show()"
   ]
  },
  {
   "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": 2
}