{
 "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": "[encoded data, 14772 characters]\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": "[encoded data, 8672 characters]\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": "[encoded data, 11580 characters]\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": "[encoded data, 29868 characters]\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
}
