simple-pv-simulator/main.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import glob\n",
    "import shutil\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "from matplotlib.ticker import FuncFormatter\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import seaborn as sns\n",
    "import json\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "\n",
    "def clear_folder_make_ess_pv(folder_path):\n",
    "    if os.path.isdir(folder_path):\n",
    "        shutil.rmtree(folder_path)\n",
    "    os.makedirs(folder_path)\n",
    "    os.makedirs(os.path.join(folder_path,'ess'))\n",
    "    os.makedirs(os.path.join(folder_path,'pv'))\n",
    "\n",
    "folder_path = 'plots'\n",
    "clear_folder_make_ess_pv(folder_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from EnergySystem import EnergySystem\n",
    "from config import pv_config, grid_config, ess_config\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Version 0.0.5\n"
     ]
    }
   ],
   "source": [
    "import json\n",
    "\n",
    "print(\"Version 0.0.5\")\n",
    "\n",
    "with open('config.json', 'r') as f:\n",
    "    js_data = json.load(f)\n",
    "\n",
    "\n",
    "    \n",
    "\n",
    "time_interval = js_data[\"time_interval\"][\"numerator\"] / js_data[\"time_interval\"][\"denominator\"]\n",
    "# print(time_interval)\n",
    "\n",
    "pv_loss = js_data[\"pv\"][\"loss\"]\n",
    "pv_cost_per_kW = js_data[\"pv\"][\"cost_per_kW\"]\n",
    "pv_lifetime = js_data[\"pv\"][\"lifetime\"]\n",
    "\n",
    "ess_loss = js_data[\"ess\"][\"loss\"]\n",
    "ess_cost_per_kW = js_data[\"ess\"][\"cost_per_kW\"]\n",
    "ess_lifetime = js_data[\"ess\"][\"lifetime\"]\n",
    "\n",
    "grid_loss = js_data[\"grid\"][\"loss\"]\n",
    "sell_price = js_data[\"grid\"][\"sell_price\"] #kWh\n",
    "grid_capacity = js_data[\"grid\"][\"capacity\"] #kWh\n",
    "\n",
    "pv_begin = js_data[\"pv_capacities\"][\"begin\"]\n",
    "pv_end = js_data[\"pv_capacities\"][\"end\"]\n",
    "pv_groups = js_data[\"pv_capacities\"][\"groups\"]\n",
    "\n",
    "ess_begin = js_data[\"ess_capacities\"][\"begin\"]\n",
    "ess_end = js_data[\"ess_capacities\"][\"end\"]\n",
    "ess_groups = js_data[\"ess_capacities\"][\"groups\"]\n",
    "\n",
    "annot_unmet = js_data[\"annotated\"][\"unmet_prob\"]\n",
    "annot_benefit = js_data[\"annotated\"][\"benefit\"]\n",
    "annot_cost = js_data[\"annotated\"][\"cost\"]\n",
    "annot_roi = js_data[\"annotated\"][\"roi\"]\n",
    "\n",
    "title_unmet = js_data[\"plot_title\"][\"unmet_prob\"]\n",
    "title_cost = js_data[\"plot_title\"][\"cost\"]\n",
    "title_benefit = js_data[\"plot_title\"][\"benefit\"]\n",
    "title_roi = js_data[\"plot_title\"][\"roi\"]\n",
    "\n",
    "\n",
    "figure_size = (js_data[\"figure_size\"][\"length\"], js_data[\"figure_size\"][\"height\"])\n",
    "\n",
    "data = pd.read_csv('combined_data.csv')\n",
    "\n",
    "granularity = js_data[\"time_interval\"][\"numerator\"]\n",
    "\n",
    "months_days = [31,28,31,30,31,30,31,31,30,31,30,31]\n",
    "def get_month_coe(num, granularity):\n",
    "    return 60 / granularity * 24 * months_days[num]\n",
    "\n",
    "months_index = [get_month_coe(num, granularity) for num in range(12)]\n",
    "months_data = []\n",
    "for i in range(1,12):\n",
    "    months_index[i] += months_index[i-1]\n",
    "for i in range(12):\n",
    "    start = 0 if i == 0 else months_index[i-1]\n",
    "    end = months_index[i]\n",
    "    months_data.append(data.iloc[int(start):int(end)])\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "pv_capacities = np.linspace(pv_begin, pv_end, pv_groups)\n",
    "ess_capacities = np.linspace(ess_begin, ess_end, ess_groups)\n",
    "# results = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "# affords = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "# costs = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "# overload_cnt = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "hour_demand = []\n",
    "for index, row in data.iterrows():\n",
    "    time = row['time']\n",
    "    demand = row['demand']\n",
    "    if time.endswith('00'):\n",
    "        hour_demand.append(demand)\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(hour_demand)\n",
    "plt.ylabel('Demand Power / kW')\n",
    "plt.savefig('plots/demand.png')\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_results(results, filename, title_benefit, annot_benefit=False, figure_size=(10, 10)):\n",
    "    df=results\n",
    "    df = df.astype(float)\n",
    "    df.index = df.index / 1000\n",
    "    df.index = df.index.map(int)\n",
    "    df.columns = df.columns / 1000\n",
    "    df.columns = df.columns.map(int)\n",
    "    min_value = df.min().min()\n",
    "    max_value = df.max().max()\n",
    "    max_scale = max(abs(min_value/1000), abs(max_value/1000))\n",
    "\n",
    "    df[df.columns[-1] + 1] = df.iloc[:, -1] \n",
    "    new_Data = pd.DataFrame(index=[df.index[-1] + 1], columns=df.columns)\n",
    "    for i in df.columns:\n",
    "        new_Data[i] = df[i].iloc[-1]\n",
    "    df = pd.concat([df, new_Data])\n",
    "\n",
    "    X, Y = np.meshgrid(np.arange(df.shape[1]), np.arange(df.shape[0]))\n",
    "\n",
    "    def fmt(x,pos):\n",
    "        return '{:.0f}'.format(x/1000)\n",
    "\n",
    "    cmap = sns.color_palette(\"coolwarm\", as_cmap=True)\n",
    "    plt.figure(figsize=figure_size)\n",
    "    ax = sns.heatmap(df/1000, fmt=\".1f\", cmap=cmap, vmin=-max_scale, vmax=max_scale, annot=annot_benefit)\n",
    "    CS = ax.contour(X, Y, df,  colors='black', alpha=0.5)\n",
    "    ax.clabel(CS, inline=True, fontsize=10, fmt=FuncFormatter(fmt))\n",
    "    plt.title(title_benefit)\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.xlim(0, df.shape[1] - 1)\n",
    "    plt.ylim(0, df.shape[0] - 1)\n",
    "    plt.xlabel('ESS Capacity (MWh)')\n",
    "    plt.ylabel('PV Capacity (MW)')\n",
    "    plt.savefig(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_roi(costs, results, filename, title_roi, days=365, annot_roi=False, figure_size=(10, 10)):\n",
    "    costs = costs.astype(float)\n",
    "    costs = costs / 365 \n",
    "    costs = costs * days\n",
    "\n",
    "    df = results\n",
    "    df = costs / df\n",
    "    if 0 in df.index and 0 in df.columns:\n",
    "        df.loc[0,0] = 100\n",
    "    df[df > 80] = 100\n",
    "    # print(df)\n",
    "\n",
    "    df = df.astype(float)\n",
    "    df.index = df.index / 1000\n",
    "    df.index = df.index.map(int)\n",
    "    df.columns = df.columns / 1000\n",
    "    df.columns = df.columns.map(int)\n",
    "    min_value = df.min().min()\n",
    "    max_value = df.max().max()\n",
    "    # print(max_value)\n",
    "    max_scale = max(abs(min_value), abs(max_value))\n",
    "\n",
    "    df[df.columns[-1] + 1] = df.iloc[:, -1] \n",
    "    new_Data = pd.DataFrame(index=[df.index[-1] + 1], columns=df.columns)\n",
    "    for i in df.columns:\n",
    "        new_Data[i] = df[i].iloc[-1]\n",
    "    df = pd.concat([df, new_Data])\n",
    "\n",
    "    X, Y = np.meshgrid(np.arange(df.shape[1]), np.arange(df.shape[0]))\n",
    "\n",
    "    def fmt(x,pos):\n",
    "        return '{:.0f}'.format(x)\n",
    "\n",
    "    cmap = sns.color_palette(\"Greys\", as_cmap=True)\n",
    "    plt.figure(figsize=figure_size)\n",
    "    ax = sns.heatmap(df, fmt=\".1f\", cmap=cmap, vmin=0, vmax=100, annot=annot_benefit)\n",
    "    CS = ax.contour(X, Y, df,  colors='black', alpha=0.5)\n",
    "    ax.clabel(CS, inline=True, fontsize=10, fmt=FuncFormatter(fmt))\n",
    "    plt.title(title_roi)\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.xlim(0, df.shape[1] - 1)\n",
    "    plt.ylim(0, df.shape[0] - 1)\n",
    "    plt.xlabel('ESS Capacity (MWh)')\n",
    "    plt.ylabel('PV Capacity (MW)')\n",
    "    plt.savefig(filename)\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_cost(costs, filename, title_cost, annot_cost=False, figure_size=(10, 10)):\n",
    "    df = costs\n",
    "    df = df.astype(int)\n",
    "    df.index = df.index / 1000\n",
    "    df.index = df.index.map(int)\n",
    "    df.columns = df.columns / 1000\n",
    "    df.columns = df.columns.map(int)\n",
    "\n",
    "    df[df.columns[-1] + 1] = df.iloc[:, -1] \n",
    "    new_Data = pd.DataFrame(index=[df.index[-1] + 1], columns=df.columns)\n",
    "    for i in df.columns:\n",
    "        new_Data[i] = df[i].iloc[-1]\n",
    "    df = pd.concat([df, new_Data])\n",
    "    X, Y = np.meshgrid(np.arange(df.shape[1]), np.arange(df.shape[0]))\n",
    "\n",
    "    def fmt(x, pos):\n",
    "        return '{:.0f}'.format(x / 1000000)\n",
    "\n",
    "    plt.figure(figsize=figure_size)\n",
    "    ax = sns.heatmap(df/1000000,  fmt=\".1f\", cmap='viridis', annot=annot_cost)\n",
    "    CS = ax.contour(X, Y, df,  colors='black', alpha=0.5)\n",
    "    ax.clabel(CS, inline=True, fontsize=10, fmt=FuncFormatter(fmt))\n",
    "    plt.title(title_cost)\n",
    "    plt.gca().invert_yaxis()\n",
    "    plt.xlim(0, df.shape[1] - 1)\n",
    "    plt.ylim(0, df.shape[0] - 1)\n",
    "    plt.xlabel('ESS Capacity (MWh)')\n",
    "    plt.ylabel('PV Capacity (MW)')\n",
    "    plt.savefig(filename)\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_overload(overload_cnt, filename, title_unmet, annot_unmet=False, figure_size=(10, 10), days=365, granularity=15):\n",
    "    df = overload_cnt\n",
    "    # print(days, granularity)\n",
    "    coef = 60 / granularity * days * 24\n",
    "    # print(coef)\n",
    "    # print(df)\n",
    "    df = ( coef - df) / coef\n",
    "    # print(df)\n",
    "\n",
    "    df = df.astype(float)\n",
    "    df.index = df.index / 1000\n",
    "    df.index = df.index.map(int)\n",
    "    df.columns = df.columns / 1000\n",
    "    df.columns = df.columns.map(int)\n",
    "\n",
    "\n",
    "    df[df.columns[-1] + 1] = df.iloc[:, -1] \n",
    "    new_Data = pd.DataFrame(index=[df.index[-1] + 1], columns=df.columns)\n",
    "    for i in df.columns:\n",
    "        new_Data[i] = df[i].iloc[-1]\n",
    "    # print(new_Data)\n",
    "    df = pd.concat([df, new_Data])\n",
    "\n",
    "\n",
    "    plt.figure(figsize=figure_size)\n",
    "    cmap = LinearSegmentedColormap.from_list(\"\", [\"white\", \"blue\"])\n",
    "    ax = sns.heatmap(df, fmt=\".00%\", cmap=cmap, vmin=0, vmax=1, annot=annot_unmet)\n",
    "\n",
    "    cbar = ax.collections[0].colorbar\n",
    "    cbar.set_ticks([0, 0.10, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1])\n",
    "    cbar.set_ticklabels(['0%', '10%', '20%', '30%', '40%', '50%', '60%', '70%', '80%', '90%', '100%'])\n",
    "    cbar.ax.yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x:.0%}'))\n",
    "    X, Y = np.meshgrid(np.arange(df.shape[1]), np.arange(df.shape[0]))\n",
    "\n",
    "    def fmt(x, pos):\n",
    "        return '{:.0f}%'.format(x * 100)\n",
    "    CS = ax.contour(X, Y, df,  colors='black', alpha=0.5)\n",
    "\n",
    "    ax.clabel(CS, inline=True, fontsize=10, fmt=FuncFormatter(fmt))\n",
    "\n",
    "    plt.xlim(0, df.shape[1] - 1)\n",
    "    plt.ylim(0, df.shape[0] - 1)\n",
    "    plt.title(title_unmet)\n",
    "    plt.xlabel('ESS Capacity (MWh)')\n",
    "    plt.ylabel('PV Capacity (MW)')\n",
    "    plt.savefig(filename)\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cal_profit(es: EnergySystem, saved_money, days):\n",
    "    profit = saved_money - es.ess.get_cost_per_year() / 365 * days - es.pv.get_cost_per_year() / 365 * days\n",
    "    return profit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_data(pv_capacity, pv_cost_per_kW, pv_lifetime, pv_loss, ess_capacity, ess_cost_per_kW, ess_lifetime, ess_loss, grid_capacity, grid_loss, sell_price, time_interval, data, days, storage=0):\n",
    "    pv = pv_config(capacity=pv_capacity, \n",
    "                    cost_per_kW=pv_cost_per_kW,\n",
    "                    lifetime=pv_lifetime, \n",
    "                    loss=pv_loss)\n",
    "    ess = ess_config(capacity=ess_capacity, \n",
    "                        cost_per_kW=ess_cost_per_kW, \n",
    "                        lifetime=ess_lifetime, \n",
    "                        loss=ess_loss,\n",
    "                        charge_power=ess_capacity,\n",
    "                        discharge_power=ess_capacity,\n",
    "                        storage=storage)\n",
    "    grid = grid_config(capacity=grid_capacity, \n",
    "                        grid_loss=grid_loss,\n",
    "                        sell_price= sell_price)\n",
    "    energySystem = EnergySystem(pv_type=pv, \n",
    "                                ess_type=ess, \n",
    "                                grid_type= grid)\n",
    "    (benefit, netto_benefit, gen_energy) = energySystem.simulate(data, time_interval)\n",
    "    results = cal_profit(energySystem, benefit, days)\n",
    "    overload_cnt = energySystem.overload_cnt\n",
    "    costs = energySystem.ess.capacity * energySystem.ess.cost_per_kW + energySystem.pv.capacity * energySystem.pv.cost_per_kW\n",
    "    return (results, \n",
    "            overload_cnt,\n",
    "            costs, \n",
    "            netto_benefit, \n",
    "            gen_energy, \n",
    "            energySystem.generated,\n",
    "            energySystem.ess_rest,\n",
    "            energySystem.factory_demand,\n",
    "            energySystem.buy_price_kWh,\n",
    "            energySystem.sell_price_kWh,\n",
    "            energySystem.pv_generated_kWh,\n",
    "            energySystem.grid_need_power_kW,\n",
    "            energySystem.time)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "  0%|          | 0/12 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_1064872/4078438238.py:64: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise an error in a future version of pandas. Value '2390.212877551013' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  months_ess_rest.loc[pv_capacity, ess_capacity] = ess_rest\n",
      "/tmp/ipykernel_1064872/4078438238.py:64: FutureWarning: Setting an item of incompatible dtype is deprecated and will raise an error in a future version of pandas. Value '2390.212877551013' has dtype incompatible with int64, please explicitly cast to a compatible dtype first.\n",
      "  months_ess_rest.loc[pv_capacity, ess_capacity] = ess_rest\n",
      "100%|██████████| 12/12 [00:19<00:00,  1.63s/it]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAD+3UlEQVR4nOzdeXwU9f0/8NfM3kc29x3IxRESTgMIcirIIaIo3lYErVqB9qt8bf35bYtHa6lHrfVEq6K2KlasF4qIyqGA3KDcdwKB3Hc2e878/ghZshIgUbKfJfN6Ph55PGB2snntZHcy7/lckqqqKoiIiIiIiASQRQcgIiIiIiLtYkFCRERERETCsCAhIiIiIiJhWJAQEREREZEwLEiIiIiIiEgYFiRERERERCQMCxIiIiIiIhKGBQkREREREQnDgoSIiIiIiIRhQdIJHD58GJIk4fXXXxcdpc1WrFgBSZKwYsWKdn/vQw89BEmSzn2on6GkpATXXHMNYmNjIUkSnn766Z/1Gs934fg7ovNLZ/1MhcP5+rLLLsMdd9wR+P/rr78OSZKwceNGYZmodTt37oRer8f27dtFRyHqUCxIzgPNfyxa+/p//+//iY4XtqZPnx50rBwOB/r164e//e1vcLvd5/Rn3XvvvVi6dCkeeOAB/Otf/8KECRNa3e/tt9/G008/3ebnzcjICHoNZrMZ3bt3x29/+1tUVlaeo/Sd3+jRo4OOY0xMDAYNGoTXXnsNiqKc05/1xhtvoHfv3rBarejSpQumTZuGY8eOtfn7z/R5lyQJ33333TnNeyYlJSW47777kJOTA6vVCpvNhvz8fPz5z39GdXV1h/7stn6mwlV7P+uhsnr1anzxxRe4//77z9lz/vhc++NzVrOzFT6XX345MjIygrb9+PkcDgdGjRqFTz/9tN05p06dissuu6zd3ydSbm4uJk2ahLlz54qOQtSh9KIDUNs98sgjyMzMDNrWu3dvpKeno7GxEQaDQVCy9hs5ciQaGxthNBo79OeYTCa88sorAIDq6mq8//77uO+++7BhwwYsXLjwnP2cr7/+GldeeSXuu+++wLYePXqc8hrffvttbN++Hffcc0+bn7t///743//9XwCAy+XCpk2b8PTTT2PlypVYv379OXsNnV1aWhrmzZsHACgrK8Obb76J22+/HXv37sVf//rXc/IzPvjgA0yfPh2jRo3C7NmzUVpaikWLFmHv3r1ISUlp13O19nkHgG7dup2TrGezYcMGXHbZZaivr8cvfvEL5OfnAwA2btyIv/71r1i1ahW++OKLDvv5bf1MhavTfdZFn6+feOIJjBkz5py/j1qea1vS6XQ/+7kvvfRSTJs2DaqqoqCgAC+++CImT56MJUuWYPz48W16Dq/Xi2XLlgXOAeeTX/3qV7jssstw4MABZGdni45D1CFYkJxHJk6ciIEDB7b6WMu7UOcDWZZDklmv1+MXv/hF4P8zZ87EhRdeiHfffRdPPfVUqxeJqqrC5XLBYrG0+eeUlpYiKioqaNu5eo2pqalBr+GXv/wl7HY7nnzySezbtw/du3f/2T9DCyIjI4OO41133YWePXviueeew5/+9KdzcoG4cOFCxMTE4PPPPw/87ufOnQuPx9Pu5zrT572jVVdX46qrroJOp8OWLVuQk5MT9Pijjz6Kf/7znx2aoSM/Uz+F0+mE1Wr92c/z41aDUCotLcWnn36K+fPnn/Pn/vG59lzq0aNH0HNPnToVubm5+Mc//tHmguSbb75BXV0dJk2a1CEZzzWfzwdFUWA0GjF27FhER0fjjTfewCOPPCI6GlGHYJetTqC1PsnTp0+H3W5HUVERpkyZArvdjvj4eNx3333w+/1B3//kk0/ioosuQmxsLCwWC/Lz87Fo0aJTfo4kSZg9ezY+/PBD9O7dGyaTCXl5efj8889P2beoqAi33347UlJSYDKZkJmZibvvvjtwYdZaX/BvvvkG1157Lbp27QqTyYQuXbrg3nvvRWNj47k5UGi6oBk9ejSApuMGNHWLuvzyy7F06VIMHDgQFosFL730EgDg4MGDuPbaaxETEwOr1YohQ4YEdRVo7oKgqiqef/75QLeC1l7j6NGj8emnn6KgoCCw34+7J7RVUlISgKaLgGbff/89pk+fjqysLJjNZiQlJeG2225DRUVF0PfW1dXhnnvuQUZGBkwmExISEnDppZdi8+bNQfutW7cOEyZMQGRkJKxWK0aNGoXVq1efkuXbb7/FoEGDYDabkZ2dHTh2bfXee+8hPz8fFosFcXFx+MUvfoGioqKgfdrzfm6r5t9nQ0MDysrK8OCDD8JgMKCsrOyUfe+8805ERUXB5XKd8TllWYbP5zvlrnBH3dFfuHAh8vPzERERAYfDgT59+uAf//hH4HGv14uHH34Y3bt3h9lsRmxsLIYPH45ly5ad8XlfeuklFBUV4amnnjqlGAGAxMRE/OEPfwja9sILLyAvLw8mkwkpKSmYNWvWKd26Ro8ejd69e2Pnzp24+OKLYbVakZqaiscffzywT3s+U82ef/55ZGVlwWKxYPDgwfjmm28wevTowGe95fM2f+6btfaczTk3bdqEkSNHwmq14v/+7/8AAB999BEmTZoUOLdlZ2fjT3/6U9D78Eyf9dONIfn6668xYsQI2Gw2REVF4corr8SuXbuC9mkem7V//35Mnz4dUVFRiIyMxIwZM+B0Ok/5Pf3Yp59+Cp/Ph7Fjx55136qqKgwePBhpaWnYs2fPWfcPpV69eiEuLg4HDhxo8/d8+umnyM3NPe059+DBg5AkCX//+99PeWzNmjWQJAnvvPNOYFtRURFuu+02JCYmBv4Wvvbaa0Hf5/F4MHfuXOTn5yMyMhI2mw0jRozA8uXLg/Zrfk88+eSTePrpp5GdnQ2TyYSdO3cCAAwGA0aPHo2PPvqoza+X6HzDFpLzSE1NDcrLy4O2xcXFnXZ/v9+P8ePH48ILL8STTz6JL7/8En/729+QnZ2Nu+++O7DfP/7xD1xxxRW4+eab4fF4sHDhQlx77bVYvHjxKXeTvv32W/z3v//FzJkzERERgWeeeQZTp05FYWEhYmNjAQDHjh3D4MGDUV1djTvvvBM5OTkoKirCokWL4HQ6T3tx9t5778HpdOLuu+9GbGws1q9fj2effRZHjx7Fe++991MP2yma/4g15wWAPXv24MYbb8Rdd92FO+64Az179kRJSQkuuugiOJ1O/OY3v0FsbCzeeOMNXHHFFVi0aBGuuuoqjBw5Ev/6179wyy23BLoVnM7vf/971NTU4OjRo4E/ena7/ax5vV5v4PfucrmwZcsWPPXUUxg5cmRQl55ly5bh4MGDmDFjBpKSkrBjxw68/PLL2LFjB7777rvARd2vfvUrLFq0CLNnz0Zubi4qKirw7bffYteuXbjgggsANF0cTZw4Efn5+XjwwQchyzIWLFiASy65BN988w0GDx4MAPjhhx8wbtw4xMfH46GHHoLP58ODDz6IxMTENv0uXn/9dcyYMQODBg3CvHnzUFJSgn/84x9YvXo1tmzZEnSHvK3v5/Y4ePAgdDodoqKicMstt+CRRx7Bu+++i9mzZwf28Xg8WLRoEaZOnXrWO9szZszAwoULMXfu3J/dNaS1z7skSYH37bJly3DjjTdizJgxeOyxxwAAu3btwurVq/E///M/AJouYOfNm4df/vKXGDx4MGpra7Fx40Zs3rwZl1566Wl/9scffwyLxYJrrrmmTVkfeughPPzwwxg7dizuvvtu7NmzBy+++CI2bNiA1atXB7U+VVVVYcKECbj66qtx3XXXYdGiRbj//vvRp08fTJw4sV2fKQB48cUXMXv2bIwYMQL33nsvDh8+jClTpiA6OhppaWltyt+aiooKTJw4ETfccAN+8YtfBN7Tr7/+Oux2O+bMmQO73Y6vv/4ac+fORW1tLZ544gkA7f+sf/nll5g4cSKysrLw0EMPobGxEc8++yyGDRuGzZs3n3IRfd111yEzMxPz5s3D5s2b8corryAhISHwPjidNWvWIDY2Funp6Wfcr7y8HJdeeikqKyuxcuXKNncT+vH7FWg
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3wUdf4/8Nds77vZTe+FUEI3dKQJUkRsqGcF0VOPol/1vPPnFdudctbz7O1EvVPxxI4iIoJU6b2XECC9b+rW+f0RsmZJ2wDJ7Cav5+ORhzI7u/vKZMu859MEURRFEBERERERSUAmdQAiIiIiIuq+WJAQEREREZFkWJAQEREREZFkWJAQEREREZFkWJAQEREREZFkWJAQEREREZFkWJAQEREREZFkWJAQEREREZFkWJAQEREREZFkWJB0ASdOnIAgCHjvvfekjhKw1atXQxAErF69ut33feyxxyAIwoUPdR4KCgpw7bXXwmazQRAEvPjii+f1O4a6YPwbUWjpqu+pYPi8vuyyy3DnnXf6/v3ee+9BEARs3bpVskzUvP3790OhUGDv3r1SRyHqUCxIQkDDl0VzP//v//0/qeMFrdtuu83vWJlMJgwcOBDPP/88HA7HBX2u+++/H8uXL8fDDz+M//znP5g6dWqz+3300Ud48cUXA37c5ORkv99Bo9EgPT0df/jDH1BaWnqB0nd948eP9zuOVqsVQ4cOxbvvvguv13tBn+v9999Hv379oNPpkJCQgFmzZiE3Nzfg+7f2fhcEAb/88ssFzduagoICPPjgg+jduzd0Oh30ej0yMzPx97//HeXl5R363IG+p4JVe9/rnWX9+vX44Ycf8NBDD12wxzz7s/bsz6wGbRU+l19+OZKTk/22nf14JpMJ48aNw7ffftvunDNnzsRll13W7vtJKSMjA9OnT8cjjzwidRSiDqWQOgAF7oknnkBKSorftn79+iEpKQm1tbVQKpUSJWu/sWPHora2FiqVqkOfR61W45133gEAlJeX47PPPsODDz6ILVu2YPHixRfseX766SdceeWVePDBB33bevbs2eR3/Oijj7B3717cd999AT/2oEGD8Pvf/x4AUFdXh23btuHFF1/Ezz//jM2bN1+w36Gri4+Px8KFCwEARUVF+OCDD3DHHXfg8OHD+Mc//nFBnuOLL77AbbfdhnHjxmHBggUoLCzEkiVLcPjwYcTGxrbrsZp7vwNAjx49LkjWtmzZsgWXXXYZqqqqcMsttyAzMxMAsHXrVvzjH//AmjVr8MMPP3TY8wf6ngpWLb3Xpf68fvbZZzFx4sQL/jpq/FnbmFwuP+/HvvTSSzFr1iyIoojs7Gy8/vrrmDFjBpYtW4YpU6YE9BgulwsrVqzwfQaEkt/97ne47LLLcOzYMaSlpUkdh6hDsCAJIdOmTcOQIUOava3xVahQIJPJOiWzQqHALbfc4vv3vHnzMHz4cHzyySd44YUXmj1JFEURdXV10Gq1AT9PYWEhLBaL37YL9TvGxcX5/Q6//e1vYTAY8Nxzz+HIkSNIT08/7+foDsxms99xvPvuu9GrVy+88sor+Nvf/nZBThAXL14Mq9WK77//3ve3f+SRR+B0Otv9WK293ztaeXk5rr76asjlcuzYsQO9e/f2u/3JJ5/E22+/3aEZOvI9dS5qamqg0+nO+3HObjXoTIWFhfj222/xxhtvXPDHPvuz9kLq2bOn32PPnDkTGRkZ+Ne//hVwQbJ27VpUVlZi+vTpHZLxQnO73fB6vVCpVJg0aRLCwsLw/vvv44knnpA6GlGHYJetLqC5Psm33XYbDAYDcnJycNVVV8FgMCAiIgIPPvggPB6P3/2fe+45jBo1CjabDVqtFpmZmViyZEmT5xEEAQsWLMCXX36Jfv36Qa1Wo2/fvvj++++b7JuTk4M77rgDsbGxUKvVSElJwdy5c30nZs31BV+7di2uu+46JCYmQq1WIyEhAffffz9qa2svzIFC/QnN+PHjAdQfN6C+W9Tll1+O5cuXY8iQIdBqtXjzzTcBAMePH8d1110Hq9UKnU6HESNG+HUVaOiCIIoiXn31VV+3guZ+x/Hjx+Pbb79Fdna2b7+zuycEKjo6GkD9SUCD3bt347bbbkNqaio0Gg2io6Nx++23o6SkxO++lZWVuO+++5CcnAy1Wo3IyEhceuml2L59u99+mzZtwtSpU2E2m6HT6TBu3DisX7++SZZ169Zh6NCh0Gg0SEtL8x27QH366afIzMyEVqtFeHg4brnlFuTk5Pjt057Xc6Aa/p7V1dUoKirCo48+CqVSiaKioib73nXXXbBYLKirq2v1MWUyGdxud5Orwh11RX/x4sXIzMyE0WiEyWRC//798a9//ct3u8vlwuOPP4709HRoNBrYbDZcfPHFWLFiRauP++abbyInJwcvvPBCk2IEAKKiovCXv/zFb9trr72Gvn37Qq1WIzY2FvPnz2/SrWv8+PHo168f9u/fjwkTJkCn0yEuLg7PPPOMb5/2vKcavPrqq0hNTYVWq8WwYcOwdu1ajB8/3vdeb/y4De/7Bs09ZkPObdu2YezYsdDpdPjTn/4EAPjqq68wffp032dbWloa/va3v/m9Dlt7r7c0huSnn37CmDFjoNfrYbFYcOWVV+LAgQN++zSMzTp69Chuu+02WCwWmM1mzJkzBzU1NU3+Tmf79ttv4Xa7MWnSpDb3LSsrw7BhwxAfH49Dhw61uX9n6tOnD8LDw3Hs2LGA7/Ptt98iIyOjxc/c48ePQxAE/POf/2xy24YNGyAIAj7++GPftpycHNx+++2IioryfRe+++67fvdzOp145JFHkJmZCbPZDL1ejzFjxmDVqlV++zW8Jp577jm8+OKLSEtLg1qtxv79+wEASqUS48ePx1dffRXw70sUathCEkIqKipQXFzsty08PLzF/T0eD6ZMmYLhw4fjueeew48//ojnn38eaWlpmDt3rm+/f/3rX7jiiitw8803w+l0YvHixbjuuuuwdOnSJleT1q1bh88//xzz5s2D0WjESy+9hJkzZ+LkyZOw2WwAgNzcXAwbNgzl5eW466670Lt3b+Tk5GDJkiWoqalp8eTs008/RU1NDebOnQubzYbNmzfj5ZdfxunTp/Hpp5+e62FrouFLrCEvABw6dAg33ngj7r77btx5553o1asXCgoKMGrUKNTU1ODee++FzWbD+++/jyuuuAJLlizB1VdfjbFjx+I///kPbr31Vl+3gpb8+c9/RkVFBU6fPu370jMYDG3mdblcvr97XV0dduzYgRdeeAFjx47169KzYsUKHD9+HHPmzEF0dDT27duHt956C/v27cMvv/ziO6n73e9+hyVLlmDBggXIyMhASUkJ1q1bhwMHDuCiiy4CUH9yNG3aNGRmZuLRRx+FTCbDokWLcMkll2Dt2rUYNmwYAGDPnj2YPHkyIiIi8Nhjj8HtduPRRx9FVFRUQH+L9957D3PmzMHQoUOxcOFCFBQU4F//+hfWr1+PHTt2+F0hD/T13B7Hjx+HXC6HxWLBrbfeiieeeAKffPIJFixY4NvH6XRiyZIlmDlzZptXtufMmYPFixfjkUceOe+uIc293wVB8L1uV6xYgRtvvBETJ07E008/DQA4cOAA1q9fj//7v/8DUH8Cu3DhQvz2t7/FsGHDYLfbsXXrVmzfvh2XXnppi8/99ddfQ6vV4tprrw0o62OPPYbHH38ckyZNwty5c3Ho0CG8/vrr2LJlC9avX+/X+lRWVoapU6fimmuuwfXXX48lS5bgoYceQv/+/TFt2rR2vacA4PXXX8eCBQswZswY3H///Thx4gSuuuoqhIWFIT4+PqD8zSkpKcG0adNwww034JZbbvG9pt977z0YDAY88MADMBgM+Omnn/DII4/Abrfj2WefBdD+9/qPP/6IadOmITU1FY899hhqa2vx8ssvY/To0di+fXuTk+jrr78eKSkpWLhwIbZv34533nkHkZGRvtdBSzZs2ACbzYakpKRW9ysuLsall16K0tJS/PzzzwF3Ezr79QrUF+Mmkymg+weqoqICZWVl7eq+9N133+Hyyy9
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3xV9f0/8NfdM/fe3Ju9SQKEPcIUmSJTFEWtG9RqK1qrtv32Z4erwzpqW61WrYqjroqjigNBWSIbUWaYCRCyb5Kb3OTu8/sj5JILCSRK7ueQ+3o+Hnk84Nxzb145d533+SyFJEkSiIiIiIiIBFCKDkBERERERLGLBQkREREREQnDgoSIiIiIiIRhQUJERERERMKwICEiIiIiImFYkBARERERkTAsSIiIiIiISBgWJEREREREJAwLEiIiIiIiEoYFSQ9QXFwMhUKBl19+WXSUTlu5ciUUCgVWrlzZ5fs+8MADUCgUZz/UD1BRUYHLL78cDocDCoUCf//733/Q33iuk+NzROeWnvqeksPn9axZs3DLLbeE///yyy9DoVBg8+bNwjJR+3bt2gW1Wo0dO3aIjkLUrViQnANavyza+/l//+//iY4nWwsWLIg4VhaLBUOGDMFf//pXeL3es/q77r77bixduhT33nsvXnvtNcyYMaPd/d544w38/e9/7/Tj5uTkRPwNer0evXv3xq9+9Ss4nc6zlL7nmzRpUsRxtNvtGDlyJF566SWEQqGz+rteeeUVDBw4EEajEZmZmbjhhhtw7NixTt//dO93hUKB9evXn9W8p1NRUYFf/vKXKCgogNFohMlkQmFhIf74xz+irq6uW393Z99TctXV93q0rF27Fp9//jl+/etfn7XHPPmz9uTPrFZnKnwuuugi5OTkRGw7+fEsFgsmTpyIjz/+uMs5582bh1mzZnX5fiL1798fs2fPxn333Sc6ClG3UosOQJ330EMPoVevXhHbBg4ciOzsbDQ3N0Oj0QhK1nUTJkxAc3MztFptt/4enU6HF154AQBQV1eHd999F7/85S+xadMmvPXWW2ft93z55Ze45JJL8Mtf/jK8rU+fPqf8jW+88QZ27NiBu+66q9OPPXToUPziF78AAHg8HmzZsgV///vfsWrVKmzcuPGs/Q09XUZGBh5++GEAQFVVFV599VXcfPPN2Lt3L/7yl7+cld/x/vvvY8GCBZg4cSLuuOMOVFZWYvHixdi7dy/S0tK69Fjtvd8BID8//6xkPZNNmzZh1qxZaGxsxHXXXYfCwkIAwObNm/GXv/wFq1evxueff95tv7+z7ym56ui9Lvrz+rHHHsMFF1xw1l9HbT9r21KpVD/4sS+88ELccMMNkCQJJSUl+Ne//oU5c+bg008/xfTp0zv1GH6/H8uWLQt/BpxLfvrTn2LWrFk4cOAA8vLyRMch6hYsSM4hM2fOxIgRI9q9re1VqHOBUqmMSma1Wo3rrrsu/P+FCxdi9OjRePvtt/HEE0+0e5IoSRI8Hg8MBkOnf09lZSVsNlvEtrP1N6anp0f8DT/+8Y9hNpvx+OOPY9++fejdu/cP/h2xwGq1RhzHn/zkJ+jbty/++c9/4g9/+MNZOUF86623YLfb8dlnn4Wf+/vuuw8+n6/Lj3W693t3q6urw6WXXgqVSoVvvvkGBQUFEbf/6U9/wr///e9uzdCd76nvo6mpCUaj8Qc/zsmtBtFUWVmJjz/+GM8+++xZf+yTP2vPpj59+kQ89rx589C/f3/84x//6HRBsmbNGjQ0NGD27NndkvFsCwQCCIVC0Gq1mDp1KuLj4/HKK6/goYceEh2NqFuwy1YP0F6f5AULFsBsNqO0tBRz586F2WxGYmIifvnLXyIYDEbc//HHH8d5550Hh8MBg8GAwsJCLF68+JTfo1AocMcdd+CDDz7AwIEDodPpMGDAAHz22Wen7FtaWoqbb74ZaWlp0Ol06NWrF2677bbwiVl7fcHXrFmDK664AllZWdDpdMjMzMTdd9+N5ubms3Og0HJCM2nSJAAtxw1o6RZ10UUXYenSpRgxYgQMBgOee+45AMDBgwdxxRVXwG63w2g0YsyYMRFdBVq7IEiShKeffjrcraC9v3HSpEn4+OOPUVJSEt7v5O4JnZWSkgKg5SSg1XfffYcFCxYgNzcXer0eKSkpuOmmm1BTUxNx34aGBtx1113IycmBTqdDUlISLrzwQmzdujVivw0bNmDGjBmwWq0wGo2YOHEi1q5de0qWr776CiNHjoRer0deXl742HXWO++8g8LCQhgMBiQkJOC6665DaWlpxD5deT13Vuvz6Xa7UVVVhfvvvx8ajQZVVVWn7HvrrbfCZrPB4/Gc9jGVSiUCgcApV4W764r+W2+9hcLCQsTFxcFisWDQoEH4xz/+Eb7d7/fjwQcfRO/evaHX6+FwOHD++edj2bJlp33c5557DqWlpXjiiSdOKUYAIDk5Gb/73e8itj3zzDMYMGAAdDod0tLScPvtt5/SrWvSpEkYOHAgdu3ahcmTJ8NoNCI9PR2PPvpoeJ+uvKdaPf3008jNzYXBYMCoUaOwZs0aTJo0Kfxeb/u4re/7Vu09ZmvOLVu2YMKECTAajfjNb34DAPjf//6H2bNnhz/b8vLy8Ic//CHidXi693pHY0i+/PJLjB8/HiaTCTabDZdccgl2794dsU/r2Kz9+/djwYIFsNlssFqtuPHGG9HU1HTK83Syjz/+GIFAAFOnTj3jvrW1tRg1ahQyMjJQVFR0xv2jqV+/fkhISMCBAwc6fZ+PP/4Y/fv37/Az9+DBg1AoFPjb3/52ym1ff/01FAoF3nzzzfC20tJS3HTTTUhOTg5/F7700ksR9/P5fLjvvvtQWFgIq9UKk8mE8ePHY8WKFRH7tb4mHn/8cfz9739HXl4edDoddu3aBQDQaDSYNGkS/ve//3X67yU617CF5BxSX1+P6urqiG0JCQkd7h8MBjF9+nSMHj0ajz/+OJYvX46//vWvyMvLw2233Rbe7x//+AcuvvhiXHvttfD5fHjrrbdwxRVXYMmSJadcTfrqq6/w3nvvYeHChYiLi8OTTz6JefPm4fDhw3A4HACAY8eOYdSoUairq8Ott96KgoIClJaWYvHixWhqaurw5Oydd95BU1MTbrvtNjgcDmzcuBFPPfUUjh49infeeef7HrZTtH6JteYFgKKiIlx99dX4yU9+gltuuQV9+/ZFRUUFzjvvPDQ1NeHOO++Ew+HAK6+8gosvvhiLFy/GpZdeigkTJuC1117D9ddfH+5W0JHf/va3qK+vx9GjR8Nfemaz+Yx5/X5/+Hn3eDz45ptv8MQTT2DChAkRXXqWLVuGgwcP4sYbb0RKSgp27tyJ559/Hjt37sT69evDJ3U//elPsXjxYtxxxx3o378/ampq8NVXX2H37t0YPnw4gJaTo5kzZ6KwsBD3338/lEolFi1ahClTpmDNmjUYNWoUAGD79u2YNm0aEhMT8cADDyAQCOD+++9HcnJyp56Ll19+GTfeeCNGjhyJhx9+GBUVFfjHP/6BtWvX4ptvvom4Qt7Z13NXHDx4ECqVCjabDddffz0eeughvP3227jjjjvC+/h8PixevBjz5s0745XtG2+8EW+99Rbuu+++H9w1pL33u0KhCL9uly1bhquvvhoXXHABHnnkEQDA7t27sXbtWvz85z8H0HIC+/DDD+PHP/4xRo0aBZfLhc2bN2Pr1q248MILO/zdH374IQwGAy6//PJOZX3ggQfw4IMPYurUqbjttttQVFSEf/3rX9i0aRPWrl0b0fpUW1uLGTNm4LLLLsOVV16JxYsX49e//jUGDRqEmTNnduk9BQD/+te/cMcdd2D8+PG4++67UVxcjLlz5yI+Ph4ZGRmdyt+empoazJw5E1dddRWuu+668Gv65Zdfhtlsxj333AOz2Ywvv/wS9913H1wuFx577DEAXX+vL1++HDNnzkRubi4eeOABNDc346mnnsK4ceOwdevWU06ir7zySvTq1QsPP/wwtm7dihdeeAFJSUnh10FHvv76azgcDmRnZ592v+rqalx44YVwOp1YtWpVp7s
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3wUdf4/8Nf2ks1ukt30nlBCQgkGIkWaIlWwYPek6Kkn4H3V886fV2x3ylnPs3csZzuxnooICFKl9w4JCem972bb/P6IWVmSkJ7J7L6ej0ceyuzs7iuTLfOeT5MJgiCAiIiIiIhIBHKxAxARERERkf9iQUJERERERKJhQUJERERERKJhQUJERERERKJhQUJERERERKJhQUJERERERKJhQUJERERERKJhQUJERERERKJhQUJERERERKJhQeIDTp8+DZlMhnfeeUfsKB22fv16yGQyrF+/vtP3ffjhhyGTyXo+VDcUFxfj6quvhtlshkwmw3PPPdet31Hq+uPfiKTFV99T/eHzetasWbjttts8/37nnXcgk8mwc+dO0TJR6w4fPgylUomDBw+KHYWoV7EgkYDmL4vWfv7f//t/YsfrtxYuXOh1rIxGI0aMGIFnnnkGjY2NPfpc99xzD1atWoUHHngA77//PmbMmNHqfh9++CGee+65Dj9uQkKC1++g1WoxcOBA/PGPf0RFRUUPpfd9kydP9jqOISEhGD16NN5++2243e4efa53330XQ4cOhV6vR2xsLObPn4+CgoIO3/9873eZTIaff/65R/OeT3FxMe677z6kpKRAr9cjICAAGRkZ+Mc//oGqqqpefe6Ovqf6q86+1/vK5s2b8cMPP+D+++/vscc897P23M+sZu0VPpdddhkSEhK8tp37eEajEZMmTcK3337b6Zzz5s3DrFmzOn0/MaWmpmL27Nl48MEHxY5C1KuUYgegjnv00UeRmJjotW3o0KGIj4+H1WqFSqUSKVnnTZw4EVarFWq1ulefR6PR4M033wQAVFVV4bPPPsN9992HHTt24OOPP+6x5/nxxx9x+eWX47777vNsGzRoUIvf8cMPP8TBgwdx9913d/ix09PT8Yc//AEAYLPZsGvXLjz33HP46aefsH379h77HXxdTEwMli1bBgAoLS3Fe++9h1tvvRXHjx/HP//5zx55ji+++AILFy7EpEmTsHTpUpSUlGDFihU4fvw4oqKiOvVYrb3fAWDAgAE9krU9O3bswKxZs1BXV4ff/OY3yMjIAADs3LkT//znP7Fhwwb88MMPvfb8HX1P9VdtvdfF/rx+6qmncMkll/T46+jsz9qzKRSKbj/2pZdeivnz50MQBOTk5OCVV17BnDlzsHLlSkyfPr1Dj+FwOLB69WrPZ4CU/O53v8OsWbNw6tQpJCcnix2HqFewIJGQmTNnYtSoUa3edvZVKCmQy+V9klmpVOI3v/mN59+LFy/GhRdeiE8++QTPPvtsqyeJgiDAZrNBp9N1+HlKSkoQFBTkta2nfsfo6Giv3+G3v/0tDAYDnn76aZw4cQIDBw7s9nP4A5PJ5HUc77jjDgwePBgvvvgi/v73v/fICeLHH3+MkJAQfP/9956//YMPPgi73d7pxzrf+723VVVV4corr4RCocCePXuQkpLidftjjz2GN954o1cz9OZ7qisaGhqg1+u7/Tjnthr0pZKSEnz77bd49dVXe/yxz/2s7UmDBg3yeux58+YhNTUV//73vztckGzcuBG1tbWYPXt2r2TsaU6nE263G2q1GlOnTkVwcDDeffddPProo2JHI+oV7LLlA1rrk7xw4UIYDAbk5+fjiiuugMFgQGhoKO677z64XC6v+z/99NMYN24czGYzdDodMjIysGLFihbPI5PJsHTpUnz55ZcYOnQoNBoN0tLS8P3337fYNz8/H7feeiuioqKg0WiQmJiIO++803Ni1lpf8I0bN+Kaa65BXFwcNBoNYmNjcc8998BqtfbMgULTCc3kyZMBNB03oKlb1GWXXYZVq1Zh1KhR0Ol0eO211wAAWVlZuOaaaxASEgK9Xo8xY8Z4dRVo7oIgCAJeeuklT7eC1n7HyZMn49tvv0VOTo5nv3O7J3RUREQEgKaTgGb79+/HwoULkZSUBK1Wi4iICNxyyy0oLy/3um9tbS3uvvtuJCQkQKPRICwsDJdeeil2797ttd+2bdswY8YMmEwm6PV6TJo0CZs3b26RZdOmTRg9ejS0Wi2Sk5M9x66jPv30U2RkZECn08FiseA3v/kN8vPzvfbpzOu5o5r/nvX19SgtLcVDDz0ElUqF0tLSFvvefvvtCAoKgs1mO+9jyuVyOJ3OFleFe+uK/scff4yMjAwEBgbCaDRi2LBh+Pe//+253eFw4JFHHsHAgQOh1WphNptx0UUXYfXq1ed93Ndeew35+fl49tlnWxQjABAeHo6//vWvXttefvllpKWlQaPRICoqCkuWLGnRrWvy5MkYOnQoDh8+jClTpkCv1yM6OhpPPvmkZ5/OvKeavfTSS0hKSoJOp0NmZiY2btyIyZMne97rZz9u8/u+WWuP2Zxz165dmDhxIvR6Pf785z8DAL766ivMnj3b89mWnJyMv//9716vw/O919saQ/Ljjz9iwoQJCAgIQFBQEC6//HIcOXLEa5/msVknT57EwoULERQUBJPJhEWLFqGhoaHF3+lc3377LZxOJ6ZOndruvpWVlcjMzERMTAyOHTvW7v59aciQIbBYLDh16lSH7/Ptt98iNTW1zc/crKwsyGQy/Otf/2px25YtWyCTyfDRRx95tuXn5+OWW25BeHi457vw7bff9rqf3W7Hgw8+iIyMDJhMJgQEBGDChAlYt26d137Nr4mnn34azz33HJKTk6HRaHD48GEAgEqlwuTJk/HVV191+Pclkhq2kEhIdXU1ysrKvLZZLJY293e5XJg+fTouvPBCPP3001izZg2eeeYZJCcn48477/Ts9+9//xtz587FTTfdBLvdjo8//hjXXHMNvvnmmxZXkzZt2oTPP/8cixcvRmBgIJ5//nnMmzcPubm5MJvNAICCggJkZmaiqqoKt99+O1JSUpCfn48VK1agoaGhzZOzTz/9FA0NDbjzzjthNpuxfft2vPDCC8jLy8Onn37a1cPWQvOXWHNeADh27BhuuOEG3HHHHbjtttswePBgFBcXY9y4cWhoaMDvf/97mM1mvPvuu5g7dy5WrFiBK6+8EhMnTsT777+Pm2++2dOtoC1/+ctfUF1djby8PM+XnsFgaDevw+Hw/N1tNhv27NmDZ599FhMnTvTq0rN69WpkZWVh0aJFiIiIwKFDh/D666/j0KFD+Pnnnz0ndb/73e+wYsUKLF26FKmpqSgvL8emTZtw5MgRXHDBBQCaTo5mzpyJjIwMPPTQQ5DL5Vi+fDkuvvhibNy4EZmZmQCAAwcOYNq0aQgNDcXDDz8Mp9OJhx56COHh4R36W7zzzjtYtGgRRo8ejWXLlqG4uBj//ve/sXnzZuzZs8frCnlHX8+dkZWVBYVCgaCgINx888149NFH8cknn2Dp0qWefex2O1asWIF58+a1e2V70aJF+Pjjj/Hggw92u2tIa+93mUzmed2uXr0aN9xwAy655BI88cQTAIAjR45g8+bN+L//+z8ATSewy5Ytw29/+1tkZmaipqYGO3fuxO7du3HppZe2+dxff/01dDodrr766g5lffjhh/HII49g6tSpuPPOO3Hs2DG88sor2LFjBzZv3uzV+lRZWYkZM2bgqquuwrXXXosVK1bg/vvvx7BhwzBz5sxOvacA4JVXXsHSpUsxYcIE3HPPPTh9+jSuuOIKBAcHIyYmpkP5W1NeXo6ZM2fi+uuvx29+8xvPa/qdd96BwWDAvffeC4PBgB9//BEPPvggampq8NRTTwHo/Ht9zZo1mDlzJpKSkvDwww/DarXihRdewPjx47F79+4WJ9HXXnstEhMTsWzZMuzevRtvvvkmwsLCPK+DtmzZsgVmsxnx8fHn3a+srAyXXnopKioq8NNPP3W4m9C5r1egqRg3Go0dun9HVVdXo7KyslPdl7777jt
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAD630lEQVR4nOzdd3xT9f4/8Fd20qRJZzrpZJSyZW9EZIqiXL2iVxC96hXQq3y91x93OK8Xr+u6t+K4bty4GLJlLxlldkH3Tlf2+f1RG1vbQtI2OUn6evrI4yEnJ+mrpxnnfT5LIgiCACIiIiIiIhFIxQ5AREREREQ9FwsSIiIiIiISDQsSIiIiIiISDQsSIiIiIiISDQsSIiIiIiISDQsSIiIiIiISDQsSIiIiIiISDQsSIiIiIiISDQsSIiIiIiISDQuSIJCbmwuJRIK33npL7Chu27RpEyQSCTZt2uTxYx944AFIJJLuD9UFJSUl+N3vfofIyEhIJBI8/fTTXfodA50//o0osATre8ofPq9nz56NW265xfXvt956CxKJBHv37hUtE7Xv2LFjkMvlOHLkiNhRiLyKBUkAaP6yaO/2//7f/xM7nt+68cYbWx0rvV6PIUOG4Mknn4TFYunWn3X33Xfjhx9+wIoVK/Duu+9i5syZ7e73/vvv4+mnn3b7eVNSUlr9Dmq1Gn369MFf/vIXVFZWdlP64DdlypRWxzEiIgIjR47Em2++CafT2a0/6+2338bAgQMREhKCXr16YeHChSgsLHT78ed7v0skEuzcubNb855PSUkJ7rnnHmRkZCAkJARarRbDhw/Hv/71L1RXV3v1Z7v7nvJXnr7XfWX79u1Yu3Yt7r333m57zt9+1v72M6vZhQqfyy67DCkpKa22/fb59Ho9Jk+ejG+++cbjnPPnz8fs2bM9fpyYMjMzMWfOHNx3331iRyHyKrnYAch9Dz30EFJTU1ttGzhwIJKTk9HY2AiFQiFSMs9NmjQJjY2NUCqVXv05KpUKr7/+OgCguroan376Ke655x7s2bMHH374Ybf9nB9//BFXXHEF7rnnHte2vn37tvkd33//fRw5cgR33XWX2889dOhQ/N///R8AwGw2Y9++fXj66aexefNm7N69u9t+h2CXmJiIlStXAgDKysrwzjvv4Oabb8bJkyfx6KOPdsvP+Pzzz3HjjTdi8uTJWLZsGUpLS7F69WqcPHkS8fHxHj1Xe+93AOjdu3e3ZL2QPXv2YPbs2airq8Mf/vAHDB8+HACwd+9ePProo9iyZQvWrl3rtZ/v7nvKX3X0Xhf78/rxxx/HJZdc0u2vo5aftS3JZLIuP/ell16KhQsXQhAE5OXl4aWXXsLcuXPx3XffYcaMGW49h81mw7p161yfAYHkT3/6E2bPno0zZ84gPT1d7DhEXsGCJIDMmjULI0aMaPe+llehAoFUKvVJZrlcjj/84Q+ufy9ZsgSjR4/GRx99hKeeeqrdk0RBEGA2m6HRaNz+OaWlpQgLC2u1rbt+x4SEhFa/wx//+EfodDo88cQTOHXqFPr06dPln9ETGAyGVsfxtttuQ79+/fD888/j4Ycf7pYTxA8//BARERH4/vvvXX/7++67D1ar1ePnOt/73duqq6tx5ZVXQiaT4cCBA8jIyGh1/yOPPILXXnvNqxm8+Z7qjIaGBoSEhHT5eX7bauBLpaWl+Oabb/Dyyy93+3P/9rO2O/Xt27fVc8+fPx+ZmZl45pln3C5Itm7ditraWsyZM8crGbub3W6H0+mEUqnEtGnTEB4ejrfffhsPPfSQ2NGIvIJdtoJAe32Sb7zxRuh0OhQUFGDevHnQ6XSIjo7GPffcA4fD0erxTzzxBMaNG4fIyEhoNBoMHz4cq1evbvNzJBIJli1bhi+++AIDBw6ESqXCgAED8P3337fZt6CgADfffDPi4+OhUqmQmpqK22+/3XVi1l5f8K1bt+Lqq69GUlISVCoVevXqhbvvvhuNjY3dc6DQdEIzZcoUAE3HDWjqFnXZZZfhhx9+wIgRI6DRaPDKK68AALKzs3H11VcjIiICISEhGDNmTKuuAs1dEARBwAsvvODqVtDe7zhlyhR88803yMvLc+332+4J7oqNjQXQdBLQ7Oeff8aNN96ItLQ0qNVqxMbG4qabbkJFRUWrx9bW1uKuu+5CSkoKVCoVjEYjLr30Uuzfv7/Vfrt27cLMmTNhMBgQEhKCyZMnY/v27W2ybNu2DSNHjoRarUZ6errr2Lnrk08+wfDhw6HRaBAVFYU//OEPKCgoaLWPJ69ndzX/Pevr61FWVob7778fCoUCZWVlbfa99dZbERYWBrPZfN7nlEqlsNvtba4Ke+uK/ocffojhw4cjNDQUer0egwYNwjPPPOO632az4cEHH0SfPn2gVqsRGRmJCRMmYN26ded93ldeeQUFBQV46qmn2hQjABATE4N//OMfrba9+OKLGDBgAFQqFeLj47F06dI23bqmTJmCgQMH4tixY7j44osREhKChIQEPPbYY659PHlPNXvhhReQlpYGjUaDUaNGYevWrZgyZYrrvd7yeZvf983ae87mnPv27cOkSZMQEhKCv/3tbwCAL7/8EnPmzHF9tqWnp+Phhx9u9To833u9ozEkP/74IyZOnAitVouwsDBcccUVyMrKarVP89is06dP48Ybb0RYWBgMBgMWL16MhoaGNn+n3/rmm29gt9sxbdq0C+5bVVWFUaNGITExESdOnLjg/r7Uv39/REVF4cyZM24/5ptvvkFmZmaHn7nZ2dmQSCT473//2+a+n376CRKJBB988IFrW0FBAW666SbExMS4vgvffPPNVo+zWq247777MHz4cBgMBmi1WkycOBEbN25stV/za+KJJ57A008/jfT0dKhUKhw7dgwAoFAoMGXKFHz55Zdu/75EgYYtJAGkpqYG5eXlrbZFRUV1uL/D4cCMGTMwevRoPPHEE1i/fj2efPJJpKen4/bbb3ft98wzz+Dyyy/H9ddfD6vVig8//BBXX3011qxZ0+Zq0rZt2/DZZ59hyZIlCA0NxbPPPov58+cjPz8fkZGRAIDCwkKMGjUK1dXVuPXWW5GRkYGCggKsXr0aDQ0NHZ6cffLJJ2hoaMDtt9+OyMhI7N69G8899xzOnTuHTz75pLOHrY3mL7HmvABw4sQJLFiwALfddhtuueUW9OvXDyUlJRg3bhwaGhpw5513IjIyEm+//TYuv/xyrF69GldeeSUmTZqEd999FzfccIOrW0FH/v73v6Ompgbnzp1zfenpdLoL5rXZbK6/u9lsxoEDB/DUU09h0qRJrbr0rFu3DtnZ2Vi8eDFiY2Nx9OhRvPrqqzh69Ch27tzpOqn705/+hNWrV2PZsmXIzMxERUUFtm3bhqysLFx00UUAmk6OZs2aheHDh+P++++HVCrFqlWrMHXqVGzduhWjRo0CABw+fBjTp09HdHQ0HnjgAdjtdtx///2IiYlx62/x1ltvYfHixRg5ciRWrlyJkpISPPPMM9i+fTsOHDjQ6gq5u69nT2RnZ0MmkyEsLAw33HADHnroIXz00UdYtmyZax+r1YrVq1dj/vz5F7yyvXjxYnz44Ye47777utw1pL33u0Qicb1u161bhwULFuCSSy7Bf/7zHwBAVlYWtm/fjj//+c8Amk5gV65ciT/+8Y8YNWoUTCYT9u7di/379+PSSy/t8Gd/9dVX0Gg0+N3vfudW1gceeAAPPvggpk2bhttvvx0nTpzASy+9hD179mD79u2tWp+qqqowc+ZMXHXVVbjmmmuwevVq3HvvvRg0aBBmzZrl0XsKAF566SUsW7YMEydOxN13343c3FzMmzcP4eHhSExMdCt/eyoqKjBr1ixce+21+MMf/uB6Tb/11lvQ6XRYvnw5dDodfvzxR9x3330wmUx4/PHHAXj+Xl+/fj1mzZqFtLQ0PPDAA2hsbMRzzz2H8ePHY//+/W1Ooq+55hqkpqZi5cqV2L9/P15//XUYjUbX66AjP/30EyIjI5GcnHze/crLy3HppZeisrISmzdvdrub0G9fr0BTMa7X6916vLtqampQVVXlUfe
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hUZd4+8Ht6zaT3BqGFJDQDCNKRKqIo9kLRRVfAXfXdfX397S6WXddddV3XrqtiWRVX7AUVkSK9Sw+BhATSe51+zu+PkCEjCSSQ5JnJ3J/ryqWcOZncOZk5c77naQpZlmUQEREREREJoBQdgIiIiIiIAhcLEiIiIiIiEoYFCRERERERCcOChIiIiIiIhGFBQkREREREwrAgISIiIiIiYViQEBERERGRMCxIiIiIiIhIGBYkREREREQkDAuSHuDEiRNQKBR46623REdpt3Xr1kGhUGDdunUd/t5HHnkECoWi80NdhJKSElx33XUIDw+HQqHAs88+e1G/o7/zxb8R+Zee+p7yhfP1FVdcgUWLFnn+/dZbb0GhUGDnzp3CMlHrDh06BLVajQMHDoiOQtSlWJD4geYPi9a+/u///k90PJ+1YMECr2NlsVgwZMgQ/OMf/4Ddbu/Un3X//ffju+++w0MPPYR3330XM2bMaHW/999/H88++2y7n7dXr15ev4Ner0e/fv3w+9//HpWVlZ2UvuebOHGi13EMCwvDiBEj8Oabb0KSpE79WW+//TYyMjJgNBqRmJiIefPmobCwsN3ff673u0KhwNatWzs177mUlJTgd7/7HVJTU2E0GmEymZCZmYm//OUvqK6u7tKf3d73lK/q6Hu9u2zatAnff/89HnzwwU57zl+ea395zmp2vsLnyiuvRK9evby2/fL5LBYLJkyYgK+//rrDOefOnYsrrriiw98nUlpaGmbNmoVly5aJjkLUpdSiA1D7PfbYY+jdu7fXtoyMDCQnJ8NqtUKj0QhK1nHjx4+H1WqFVqvt0p+j0+nw+uuvAwCqq6vx8ccf43e/+x127NiBFStWdNrP+fHHH3H11Vfjd7/7nWdb//79z/od33//fRw4cAD33Xdfu5976NCh+J//+R8AgM1mw65du/Dss89i/fr12L59e6f9Dj1dQkICnnjiCQBAWVkZ3nnnHdx55504evQo/va3v3XKz/j000+xYMECTJgwAUuXLkVpaSlWrlyJo0ePIi4urkPP1dr7HQD69u3bKVnPZ8eOHbjiiitQX1+P2267DZmZmQCAnTt34m9/+xs2bNiA77//vst+fnvfU76qrfe66PP1U089hcsvv7zTX0ctz7UtqVSqi37uqVOnYt68eZBlGXl5eXj55Zcxe/ZsrFq1CtOnT2/XczidTqxevdpzDvAnv/71r3HFFVfg+PHj6NOnj+g4RF2CBYkfmTlzJoYPH97qYy3vQvkDpVLZLZnVajVuu+02z78XL16MSy+9FB9++CGeeeaZVi8SZVmGzWaDwWBo988pLS1FSEiI17bO+h3j4+O9fodf/epXMJvNePrpp5GdnY1+/fpd9M8IBMHBwV7H8e6778aAAQPwwgsv4M9//nOnXCCuWLECYWFh+Pbbbz1/+2XLlsHhcHT4uc71fu9q1dXVuOaaa6BSqbBnzx6kpqZ6Pf7444/j3//+d5dm6Mr31IVobGyE0Wi86Of5ZatBdyotLcXXX3+NV155pdOf+5fn2s7Uv39/r+eeO3cu0tLS8K9//avdBclPP/2Euro6zJo1q0sydjaXywVJkqDVajFlyhSEhobi7bffxmOPPSY6GlGXYJetHqC1PskLFiyA2WxGQUEB5syZA7PZjMjISPzud7+D2+32+v6nn34al112GcLDw2EwGJCZmYmVK1ee9XMUCgWWLl2Kzz77DBkZGdDpdEhPT8e333571r4FBQW48847ERcXB51Oh969e+Oee+7xXJi11hf8p59+wvXXX4+kpCTodDokJibi/vvvh9Vq7ZwDhaYLmokTJwJoOm5AU7eoK6+8Et999x2GDx8Og8GAV199FQCQk5OD66+/HmFhYTAajRg1apRXV4HmLgiyLOPFF1/0dCto7XecOHEivv76a+Tl5Xn2+2X3hPaKiYkB0HQR0Gzfvn1YsGABUlJSoNfrERMTgzvuuAMVFRVe31tXV4f77rsPvXr1gk6nQ1RUFKZOnYrdu3d77bdt2zbMmDEDwcHBMBqNmDBhAjZt2nRWlo0bN2LEiBHQ6/Xo06eP59i110cffYTMzEwYDAZERETgtttuQ0FBgdc+HXk9t1fz37OhoQFlZWV4+OGHodFoUFZWdta+d911F0JCQmCz2c75nEqlEi6X66y7wl11R3/FihXIzMxEUFAQLBYLBg0ahH/961+ex51OJx599FH069cPer0e4eHhGDt2LFavXn3O53311VdRUFCAZ5555qxiBACio6Pxxz/+0WvbSy+9hPT0dOh0OsTFxWHJkiVndeuaOHEiMjIycOjQIUyaNAlGoxHx8fF48sknPft05D3V7MUXX0RKSgoMBgNGjhyJn376CRMnTvS811s+b/P7vllrz9mcc9euXRg/fjyMRiP+3//7fwCAzz//HLNmzfKc2/r06YM///nPXq/Dc73X2xpD8uOPP2LcuHEwmUwICQnB1VdfjcOHD3vt0zw269ixY1iwYAFCQkIQHByMhQsXorGx8ay/0y99/fXXcLlcmDJlynn3raqqwsiRI5GQkICsrKzz7t+dBg4ciIiICBw/frzd3/P1118jLS2tzXNuTk4OFAoF/vnPf5712ObNm6FQKPDBBx94thUUFOCOO+5AdHS057PwzTff9Po+h8OBZcuWITMzE8HBwTCZTBg3bhzWrl3rtV/za+Lpp5/Gs88+iz59+kCn0+HQoUMAAI1Gg4kTJ+Lzzz9v9+9L5G/YQuJHampqUF5e7rUtIiKizf3dbjemT5+OSy+9FE8//TR++OEH/OMf/0CfPn1wzz33ePb717/+hauuugq33norHA4HVqxYgeuvvx5fffXVWXeTNm7ciE8++QSLFy9GUFAQnnvuOcydOxf5+fkIDw8HABQWFmLkyJGorq7GXXfdhdTUVBQUFGDlypVobGxs8+Lso48+QmNjI+655x6Eh4dj+/bteP7553Hq1Cl89NFHF3rYztL8IdacFwCysrJw88034+6778aiRYswYMAAlJSU4LLLLkNjYyN+85vfIDw8HG+//TauuuoqrFy5Etdccw3Gjx+Pd999F7fffrunW0Fb/vCHP6CmpganTp3yfOiZzebz5nU6nZ6/u81mw549e/DMM89g/PjxXl16Vq9ejZycHCxcuBAxMTE4ePAgXnvtNRw8eBBbt271XNT9+te/xsqVK7F06VKkpaWhoqICGzduxOHDh3HJJZcAaLo4mjlzJjIzM/Hwww9DqVRi+fLlmDx5Mn766SeMHDkSALB//35MmzYNkZGReOSRR+ByufDwww8jOjq6XX+Lt956CwsXLsSIESPwxBNPoKSkBP/617+wadMm7Nmzx+sOeXtfzx2Rk5MDlUqFkJAQ3H777Xjsscfw4YcfYunSpZ59HA4HVq5ciblz5573zvbChQuxYsUKLFu27KK7hrT2flcoFJ7X7erVq3HzzTfj8ssvx9///ncAwOHDh7Fp0yb89re/BdB0AfvEE0/gV7/6FUaOHIna2lrs3LkTu3fvxtSpU9v82V988QUMBgOuu+66dmV95JFH8Oijj2LKlCm45557kJWVhZdffhk7duzApk2bvFqfqqqqMGPGDFx77bW44YYbsHLlSjz44IMYNGgQZs6c2aH3FAC8/PLLWLp0KcaNG4f7778fJ06cwJw5cxAaGoqEhIR25W9NRUUFZs6ciZtuugm33Xab5zX91ltvwWw244EHHoDZbMaPP/6IZcuWoba2Fk899RSAjr/Xf/jhB8ycORMpKSl45JFHYLVa8fzzz2PMmDHYvXv3WRfRN9xwA3r37o0nnngCu3fvxuuvv46oqCjP66AtmzdvRnh4OJKTk8+5X3l5OaZOnYrKykqsX7++3d2Efvl6BZqKcYvF0q7vb6+
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hT9/k28FtbsmXZlm15b+MNGMwmbMIIIYsmzWiAJE3SAM0voSPN2zarTdLMpplkktFssicBwt57Y4P33ntonvcPx6oNNkjG9rGk+3NdvjBHR/LtY43znO+SCIIggIiIiIiISARSsQMQEREREZHnYkFCRERERESiYUFCRERERESiYUFCRERERESiYUFCRERERESiYUFCRERERESiYUFCRERERESiYUFCRERERESiYUFCRERERESiYUHiBvLz8yGRSPD222+LHcVhmzZtgkQiwaZNm5y+70MPPQSJRNL/oS5CRUUFfvWrXyEgIAASiQTPPffcRf2Orm4o/o3Itbjra2oovF9fdtlluP322+3/f/vttyGRSLBv3z7RMlHPTpw4AblcjmPHjokdhWhAsSBxAZ0fFj19/eUvfxE73pC1dOnSbsdKp9Nh5MiReOaZZ2A0Gvv1Z917771Yu3Yt7r//frz33nuYN29ej/t98MEHeO655xx+3JiYmG6/g1qtxrBhw/CnP/0JtbW1/ZTe/U2fPr3bcdTr9Rg7dizeeust2Gy2fv1Z77zzDtLT0+Hl5YXIyEgsXrwYpaWlDt//fK93iUSCXbt29Wve86moqMAf//hHJCcnw8vLC97e3sjMzMQ///lP1NfXD+jPdvQ1NVQ5+1ofLNu3b8dPP/2E++67r98e8+z32rPfszpdqPC5/PLLERMT023b2Y+n0+kwbdo0fPfdd07nXLRoES677DKn7yem1NRULFiwAA888IDYUYgGlFzsAOS4Rx55BLGxsd22paenIzo6Gm1tbVAoFCIlc97UqVPR1tYGpVI5oD9HpVLhjTfeAADU19fjs88+wx//+Efs3bsXH330Ub/9nJ9//hlXXnkl/vjHP9q3JSYmnvM7fvDBBzh27Bjuuecehx87IyMDf/jDHwAA7e3t2L9/P5577jls3rwZe/bs6bffwd1FRETg8ccfBwBUVVXh3XffxW233Ybs7Gz861//6pef8cUXX2Dp0qWYNm0aVqxYgcrKSqxZswbZ2dkICwtz6rF6er0DQEJCQr9kvZC9e/fisssuQ3NzM37zm98gMzMTALBv3z7861//wpYtW/DTTz8N2M939DU1VPX2Whf7/fqpp57CrFmz+v151PW9tiuZTHbRj33ppZdi8eLFEAQBBQUFeOWVV7Bw4UL88MMPmDt3rkOPYTabsW7dOvt7gCv53e9+h8suuww5OTmIj48XOw7RgGBB4kLmz5+PMWPG9Hhb16tQrkAqlQ5KZrlcjt/85jf2/y9btgzjx4/Hxx9/jGeffbbHk0RBENDe3g6NRuPwz6msrISfn1+3bf31O4aHh3f7HX77299Cq9Xi6aefxunTpzFs2LCL/hmewNfXt9txvPPOO5GUlIQXX3wR//jHP/rlBPGjjz6CXq/Hjz/+aP/bP/DAAzCZTE4/1vle7wOtvr4eV199NWQyGQ4ePIjk5ORutz/66KN4/fXXBzTDQL6m+qK1tRVeXl4X/ThntxoMpsrKSnz33XdYtWpVvz/22e+1/SkxMbHbYy9atAipqan4z3/+43BBsnXrVjQ1NWHBggUDkrG/WSwW2Gw2KJVKzJ49G/7+/njnnXfwyCOPiB2NaECwy5Yb6KlP8tKlS6HValFSUoKrrroKWq0WQUFB+OMf/wir1drt/k8//TQmTZqEgIAAaDQaZGZmYs2aNef8HIlEghUrVuDLL79Eeno6VCoV0tLS8OOPP56zb0lJCW677TaEhYVBpVIhNjYWd911l/3ErKe+4Fu3bsW1116LqKgoqFQqREZG4t5770VbW1v/HCh0nNBMnz4dQMdxAzq6RV1++eVYu3YtxowZA41Gg1dffRUAkJubi2uvvRZ6vR5eXl6YMGFCt64CnV0QBEHASy+9ZO9W0NPvOH36dHz33XcoKCiw73d29wRHhYSEAOg4Ceh05MgRLF26FHFxcVCr1QgJCcGtt96KmpqabvdtamrCPffcg5iYGKhUKhgMBlx66aU4cOBAt/12796NefPmwdfXF15eXpg2bRq2b99+TpZt27Zh7NixUKvViI+Ptx87R3366afIzMyERqNBYGAgfvOb36CkpKTbPs48nx3V+fdsaWlBVVUVHnzwQSgUClRVVZ2z7x133AE/Pz+0t7ef9zGlUiksFss5V4UH6or+Rx99hMzMTPj4+ECn02H48OH4z3/+Y7/dbDbj4YcfxrBhw6BWqxEQEIBLLrkE69atO+/jvvrqqygpKcGzzz57TjECAMHBwfjb3/7WbdvLL7+MtLQ0qFQqhIWFYfny5ed065o+fTrS09Nx4sQJzJgxA15eXggPD8eTTz5p38eZ11Snl156CXFxcdBoNBg3bhy2bt2K6dOn21/rXR+383XfqafH7My5f/9+TJ06FV5eXvh//+//AQC++uorLFiwwP7eFh8fj3/84x/dnofne633Nobk559/xpQpU+Dt7Q0/Pz9ceeWVOHnyZLd9OsdmnTlzBkuXLoWfnx98fX1xyy23oLW19Zy/09m+++47WCwWzJ49+4L71tXVYdy4cYiIiEBWVtYF9x9MKSkpCAwMRE5OjsP3+e6775Camtrre25ubi4kEgn+/e9/n3Pbjh07IJFI8OGHH9q3lZSU4NZbb0VwcLD9s/Ctt97qdj+TyYQHHngAmZmZ8PX1hbe3N6ZMmYKNGzd226/zOfH000/jueeeQ3x8PFQqFU6cOAEAUCgUmD59Or766iuHf18iV8MWEhfS0NCA6urqbtsCAwN73d9qtWLu3LkYP348nn76aaxfvx7PPPMM4uPjcdddd9n3+89//oMrrrgCN910E0wmEz766CNce+21+Pbbb8+5mrRt2zZ8/vnnWLZsGXx8fPD8889j0aJFKCwsREBAAACgtLQU48aNQ319Pe644w4kJyejpKQEa9asQWtra68nZ59++ilaW1tx1113ISAgAHv27MELL7yA4uJifPrpp309bOfo/BDrzAsAWVlZuOGGG3DnnXfi9ttvR1JSEioqKjBp0iS0trbi7rvvRkBAAN555x1cccUVWLNmDa6++mpMnToV7733Hm6++WZ7t4Le/PWvf0VDQwOKi4vtH3parfaCec1ms/3v3t7ejoMHD+LZZ5/F1KlTu3XpWbduHXJzc3HLLbcgJCQEx48fx2uvvYbjx49j165d9pO63/3ud1izZg1WrFiB1NRU1NTUYNu2bTh58iRGjx4NoOPkaP78+cjMzMSDDz4IqVSK1atXY+bMmdi6dSvGjRsHADh69CjmzJmDoKAgPPTQQ7BYLHjwwQcRHBzs0N/i7bffxi233IKxY8fi8ccfR0VFBf7zn/9g+/btOHjwYLcr5I4+n52Rm5sLmUwGPz8/3HzzzXjkkUfw8ccfY8WKFfZ9TCYT1qxZg0WLFl3wyvYtt9yCjz76CA888MBFdw3p6fUukUjsz9t169bhhhtuwKxZs/DEE08AAE6ePInt27fj//7v/wB0nMA+/vjj+O1vf4tx48ahsbER+/btw4EDB3DppZf2+rO//vpraDQa/OpXv3Io60MPPYSHH34Ys2fPxl133YWsrCy88sor2Lt3L7Zv396t9amurg7z5s3DNddcg+uuuw5r1qzBfffdh+HDh2P+/PlOvaYA4JVXXsGKFSswZcoU3HvvvcjPz8dVV10Ff39/REREOJS/JzU1NZg/fz6uv/56/OY3v7E/p99++21otVqsXLkSWq0WP//8Mx544AE0NjbiqaeeAuD8a339+vWYP38+4uLi8NBDD6GtrQ0vvPACJk+ejAMHDpxzEn3dddchNjYWjz/+OA4cOIA33ngDBoPB/jzozY4dOxAQEIDo6Ojz7lddXY1LL70UtbW12Lx5s8PdhM5+vgIdxbhOp3Po/o5qaGhAXV2
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3wUdf4/8Nds32xNNr0ntAChGar0IkVEUdQTPRH00BM4f8p559crtruTO9t56ilWLGfHjiIqKEjv0kIPCaT3nq3z+yNkZSWB3WSTySSv5+OxjwfMzu6+MtndzHs+TRBFUQQREREREZEEFFIHICIiIiKi7osFCRERERERSYYFCRERERERSYYFCRERERERSYYFCRERERERSYYFCRERERERSYYFCRERERERSYYFCRERERERSYYFCRERERERSYYFSRdw6tQpCIKA119/Xeoofvvhhx8gCAJ++OGHgB/70EMPQRCE4Idqg8LCQlx77bWw2WwQBAFPP/10m35GueuMvyOSl676meoM39eXX345Fi5c6P3/66+/DkEQsHPnTskyUfMOHToElUqFAwcOSB2FqF2xIJGBpj8Wzd3+7//+T+p4ndb8+fN9jpXZbMagQYPw5JNPwm63B/W17rnnHqxZswb3338/3nrrLUyfPr3Z/d555x08/fTTfj9vcnKyz8+g0+nQq1cv/OEPf0BZWVmQ0nd9EyZM8DmOYWFhGDZsGF577TV4PJ6gvtYbb7yB9PR0hISEICEhAfPmzUNeXp7fj7/Q510QBGzdujWoeS+ksLAQ9957L9LS0hASEgKDwYCMjAz8/e9/R0VFRbu+tr+fqc4q0M96R9m0aRO++eYb3HfffUF7zl9+1/7yO6vJxQqfK664AsnJyT7bfvl8ZrMZ48ePx5dffhlwzjlz5uDyyy8P+HFS6tevH2bOnIkHHnhA6ihE7UoldQDy3yOPPIKUlBSfbenp6UhKSkJ9fT3UarVEyQI3btw41NfXQ6PRtOvraLVavPLKKwCAiooKfPTRR7j33nuxY8cOvPfee0F7nXXr1uGqq67Cvffe693Wu3fv837Gd955BwcOHMDdd9/t93MPHjwYv//97wEADQ0N2LVrF55++mmsX78e27dvD9rP0NXFx8dj2bJlAIDi4mK8+eabuO2223D06FH885//DMprfPLJJ5g/fz7Gjx+PJUuWoKioCCtXrsTRo0cRGxsb0HM193kHgJ49ewYl68Xs2LEDl19+OWpqavDrX/8aGRkZAICdO3fin//8JzZs2IBvvvmm3V7f389UZ9XSZ13q7+vHH38ckydPDvr76Nzv2nMplco2P/dll12GefPmQRRFZGdn44UXXsCsWbOwevVqTJs2za/ncDqd+Pbbb73fAXLy29/+FpdffjlOnDiBHj16SB2HqF2wIJGRGTNmYOjQoc3ed+5VKDlQKBQdklmlUuHXv/619/+LFi3CiBEj8P777+Opp55q9iRRFEU0NDRAr9f7/TpFRUWwWq0+24L1M8bFxfn8DL/5zW9gNBrxxBNP4NixY+jVq1ebX6M7sFgsPsfxjjvuQJ8+ffDcc8/hb3/7W1BOEN977z2EhYXh66+/9v7uH3jgATgcjoCf60Kf9/ZWUVGBq6++GkqlEnv27EFaWprP/f/4xz/w8ssvt2uG9vxMtUZdXR1CQkLa/Dy/bDXoSEVFRfjyyy+xfPnyoD/3L79rg6l3794+zz1nzhz069cP//nPf/wuSH788UdUV1dj5syZ7ZIx2FwuFzweDzQaDaZMmYLQ0FC88cYbeOSRR6SORtQu2GWrC2iuT/L8+fNhNBqRm5uL2bNnw2g0IiIiAvfeey/cbrfP45944glceumlsNls0Ov1yMjIwMqVK897HUEQsGTJEnz66adIT0+HVqtF//798fXXX5+3b25uLm677TbExsZCq9UiJSUFd955p/fErLm+4D/++COuu+46JCYmQqvVIiEhAffccw/q6+uDc6DQeEIzYcIEAI3HDWjsFnXFFVdgzZo1GDp0KPR6PV588UUAwMmTJ3HdddchLCwMISEhGDlypE9XgaYuCKIo4r///a+3W0FzP+OECRPw5ZdfIjs727vfL7sn+Cs6OhpA40lAk3379mH+/PlITU2FTqdDdHQ0br31VpSWlvo8trq6GnfffTeSk5Oh1WoRGRmJyy67DLt37/bZb9u2bZg+fTosFgtCQkIwfvx4bNq06bwsGzduxLBhw6DT6dCjRw/vsfPXhx9+iIyMDOj1eoSHh+PXv/41cnNzffYJ5P3sr6bfZ21tLYqLi/Hggw9CrVajuLj4vH1vv/12WK1WNDQ0XPA5FQoFXC7XeVeF2+uK/nvvvYeMjAyYTCaYzWYMGDAA//nPf7z3O51OPPzww+jVqxd0Oh1sNhvGjBmDb7/99oLP++KLLyI3NxdPPfXUecUIAERFReEvf/mLz7bnn38e/fv3h1arRWxsLBYvXnxet64JEyYgPT0dhw4dwsSJExESEoK4uDg89thj3n0C+Uw1+e9//4vU1FTo9XoMHz4cP/74IyZMmOD9rJ/7vE2f+ybNPWdTzl27dmHcuHEICQnBn/70JwDAZ599hpkzZ3q/23r06IG//e1vPu/DC33WWxpDsm7dOowdOxYGgwFWqxVXXXUVMjMzffZpGpt1/PhxzJ8/H1arFRaLBQsWLEBdXd15v6df+vLLL+FyuTBlypSL7lteXo7hw4cjPj4eR44cuej+Halv374IDw/HiRMn/H7Ml19+iX79+rX4nXvy5EkIgoB///vf5923efNmCIKAd99917stNzcXt956K6Kiorx/C1977TWfxzkcDjzwwAPIyMiAxWKBwWDA2LFj8f333/vs1/SeeOKJJ/D000+jR48e0Gq1OHToEABArVZjwoQJ+Oyzz/z+eYnkhi0kMlJZWYmSkhKfbeHh4S3u73a7MW3aNIwYMQJPPPEEvvvuOzz55JPo0aMH7rzzTu9+//nPf3DllVfipptugsPhwHvvvYfrrrsOq1atOu9q0saNG/Hxxx9j0aJFMJlMeOaZZzBnzhzk5OTAZrMBAPLy8jB8+HBUVFTg9ttvR1paGnJzc7Fy5UrU1dW1eHL24Ycfoq6uDnfeeSdsNhu2b9+OZ599FmfOnMGHH37Y2sN2nqY/Yk15AeDIkSOYO3cu7rjjDixcuBB9+vRBYWEhLr30UtTV1eGuu+6CzWbDG2+8gSuvvBIrV67E1VdfjXHjxuGtt97CzTff7O1W0JI///nPqKysxJkzZ7x/9IxG40XzOp1O7++9oaEBe/bswVNPPYVx48b5dOn59ttvcfLkSSxYsADR0dE4ePAgXnrpJRw8eBBbt271ntT99re/xcqVK7FkyRL069cPpaWl2LhxIzIzM3HJJZcAaDw5mjFjBjIyMvDggw9CoVBgxYoVmDRpEn788UcMHz4cALB//35MnToVEREReOihh+ByufDggw8iKirKr9/F66+/jgULFmDYsGFYtmwZCgsL8Z///AebNm3Cnj17fK6Q+/t+DsTJkyehVCphtVpx880345FHHsH777+PJUuWePdxOBxYuXIl5syZc9Er2wsWLMB7772HBx54oM1dQ5r7vAuC4H3ffvvtt5g7dy4mT56Mf/3rXwCAzMxMbNq0Cf/v//0/AI0nsMuWLcNvfvMbDB8+HFVVVdi5cyd2796Nyy67rMXX/vzzz6HX63Httdf6lfWhhx7Cww8/jClTpuDOO+/EkSNH8MILL2DHjh3YtGmTT+tTeXk5pk+fjmuuuQbXX389Vq5cifvuuw8DBgzAjBkzAvpMAcALL7yAJUuWYOzYsbjnnntw6tQpzJ49G6GhoYiPj/crf3NKS0sxY8YM3HDDDfj1r3/tfU+//vrrMBqNWLp0KYxGI9atW4cHHngAVVVVePzxxwEE/ln/7rvvMGPGDKSmpuKhhx5CfX09nn32WYwePRq7d+8+7yT6+uuvR0pKCpYtW4bdu3fjlVdeQWRkpPd90JLNmzfDZrMhKSnpgvuVlJTgsssuQ1lZGdavX+93N6Ffvl+BxmLcbDb79Xh/VVZWory8PKD
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3hV9f0H8PfdM/cmudkkhCSMEDZhI0uQIaIo1ToRtWoFtErbn7XD1VrrqLVa90BtVaw46kZFmbKn7L0C2Xvceb6/P0IuuSSBEJN8D7nv1/PkeeDcc2/eOXedz/kujRBCgIiIiIiISAKt7ABERERERBS+WJAQEREREZE0LEiIiIiIiEgaFiRERERERCQNCxIiIiIiIpKGBQkREREREUnDgoSIiIiIiKRhQUJERERERNKwICEiIiIiImlYkHQAhw4dgkajwRtvvCE7SrMtWbIEGo0GS5YsOef7Pvjgg9BoNK0f6ifIy8vDz372M7hcLmg0Gjz99NM/6W8836nxOaLzS0d9T6nh8/riiy/GrbfeGvz/G2+8AY1Gg/Xr10vLRI3bsWMH9Ho9tm3bJjsKUZtiQXIeqPuyaOznd7/7nex4qjVr1qyQY+VwONCvXz/8/e9/h8fjadXfdc8992DRokW477778O9//xuTJ09udL933nkHTz/9dLMft0uXLiF/g9lsRrdu3fDb3/4WxcXFrZS+4xs7dmzIcYyOjsbgwYPx+uuvQ1GUVv1db775Jnr37g2r1YqUlBTMnDkTx48fb/b9z/R+12g0WL16davmPZO8vDz85je/QWZmJqxWK2w2G7Kzs/GXv/wFpaWlbfq7m/ueUqtzfa+3l5UrV+Lrr7/Gvffe22qPefpn7emfWXXOVvhccskl6NKlS8i20x/P4XBgzJgx+Pzzz88554wZM3DxxRef8/1kysrKwtSpU3H//ffLjkLUpvSyA1DzPfzww0hLSwvZ1rt3b6SmpqKmpgYGg0FSsnM3evRo1NTUwGg0tunvMZlMePXVVwEApaWl+OCDD/Cb3/wG69atw4IFC1rt93z33Xe47LLL8Jvf/Ca4rXv37g3+xnfeeQfbtm3D3Xff3ezH7t+/P379618DANxuNzZs2ICnn34aS5cuxdq1a1vtb+jokpOT8eijjwIACgoK8NZbb+GWW27Bnj178Le//a1VfsdHH32EWbNmYcyYMZg7dy7y8/OxcOFC7NmzB0lJSef0WI293wGga9eurZL1bNatW4eLL74YlZWVuP7665GdnQ0AWL9+Pf72t79h2bJl+Prrr9vs9zf3PaVWTb3XZX9eP/HEExg/fnyrv47qf9bWp9PpfvJjX3TRRZg5cyaEEDh8+DBeeOEFTJs2DV9++SUmTZrUrMfw+Xz45ptvgp8B55Nf/vKXuPjii7F//35kZGTIjkPUJliQnEemTJmCQYMGNXpb/atQ5wOtVtsumfV6Pa6//vrg/2fPno2hQ4fivffew1NPPdXoSaIQAm63GxaLpdm/Jz8/H5GRkSHbWutv7NSpU8jf8Itf/AJ2ux1PPvkk9u7di27duv3k3xEOnE5nyHG8/fbb0aNHD/zrX//Cn//851Y5QVywYAGio6Px1VdfBZ/7+++/H16v95wf60zv97ZWWlqKyy+/HDqdDps2bUJmZmbI7Y888gheeeWVNs3Qlu+plqiurobVav3Jj3N6q0F7ys/Px+eff44XX3yx1R/79M/a1tS9e/eQx54xYwaysrLwz3/+s9kFyfLly1FRUYGpU6e2ScbW5vf7oSgKjEYjJkyYgKioKLz55pt4+OGHZUcjahPsstUBNNYnedasWbDb7cjJycH06dNht9sRGxuL3/zmNwgEAiH3f/LJJzFixAi4XC5YLBZkZ2dj4cKFDX6PRqPB3Llz8fHHH6N3794wmUzo1asXvvrqqwb75uTk4JZbbkFSUhJMJhPS0tJwxx13BE/MGusLvnz5clx55ZXo3LkzTCYTUlJScM8996CmpqZ1DhRqT2jGjh0LoPa4AbXdoi655BIsWrQIgwYNgsViwUsvvQQAOHDgAK688kpER0fDarVi2LBhIV0F6rogCCHw3HPPBbsVNPY3jh07Fp9//jkOHz4c3O/07gnNlZCQAKD2JKDO1q1bMWvWLKSnp8NsNiMhIQE333wzioqKQu5bUVGBu+++G126dIHJZEJcXBwuuugibNy4MWS/NWvWYPLkyXA6nbBarRgzZgxWrlzZIMuKFSswePBgmM1mZGRkBI9dc73//vvIzs6GxWJBTEwMrr/+euTk5ITscy6v5+aqez6rqqpQUFCABx54AAaDAQUFBQ32ve222xAZGQm3233Gx9RqtfD7/Q2uCrfVFf0FCxYgOzsbERERcDgc6NOnD/75z38Gb/f5fHjooYfQrVs3mM1muFwuXHDBBfjmm2/O+LgvvfQScnJy8NRTTzUoRgAgPj4ef/zjH0O2Pf/88+jVqxdMJhOSkpIwZ86cBt26xo4di969e2PHjh0YN24crFYrOnXqhMcffzy4z7m8p+o899xzSE9Ph8ViwZAhQ7B8+XKMHTs2+F6v/7h17/s6jT1mXc4NGzZg9OjRsFqt+P3vfw8A+N///oepU6cGP9syMjLw5z//OeR1eKb3elNjSL777juMGjUKNpsNkZGRuOyyy7Bz586QferGZu3btw+zZs1CZGQknE4nbrrpJlRXVzd4nk73+eefw+/3Y8KECWfdt6SkBEOGDEFycjJ279591v3bU8+ePRETE4P9+/c3+z6ff/45srKymvzMPXDgADQaDf7xj380uO2HH36ARqPBu+++G9yWk5ODm2++GfHx8cHvwtdffz3kfl6vF/fffz+ys7PhdDphs9kwatQofP/99yH71b0mnnzySTz99NPIyMiAyWTCjh07AAAGgwFjx47F//73v2b/vUTnG7aQnEfKyspQWFgYsi0mJqbJ/QOBACZNmoShQ4fiySefxLfffou///3vyMjIwB133BHc75///CcuvfRSXHfddfB6vViwYAGuvPJKfPbZZw2uJq1YsQIffvghZs+ejYiICDzzzDOYMWMGjhw5ApfLBQA4fvw4hgwZgtLSUtx2223IzMxETk4OFi5ciOrq6iZPzt5//31UV1fjjjvugMvlwtq1a/Hss8/i2LFjeP/991t62Bqo+xKrywsAu3fvxjXXXIPbb78dt956K3r06IG8vDyMGDEC1dXVuOuuu+ByufDmm2/i0ksvxcKFC3H55Zdj9OjR+Pe//40bbrgh2K2gKX/4wx9QVlaGY8eOBb/07Hb7WfP6fL7g8+52u7Fp0yY89dRTGD16dEiXnm+++QYHDhzATTfdhISEBGzfvh0vv/wytm/fjtWrVwdP6n75y19i4cKFmDt3LrKyslBUVIQVK1Zg586dGDhwIIDak6MpU6YgOzsbDzzwALRaLebPn48LL7wQy5cvx5AhQwAAP/74IyZOnIjY2Fg8+OCD8Pv9eOCBBxAfH9+s5+KNN97ATTfdhMGDB+PRRx9FXl4e/vnPf2LlypXYtGlTyBXy5r6ez8WBAweg0+kQGRmJG264AQ8//DDee+89zJ07N7iP1+vFwoULMWPGjLNe2b7pppuwYMEC3H///T+5a0hj73eNRhN83X7zzTe45pprMH78eDz22GMAgJ07d2LlypX41a9+BaD2BPbRRx/FL37xCwwZMgTl5eVYv349Nm7ciIsuuqjJ3/3JJ5/AYrHgZz/7WbOyPvjgg3jooYcwYcIE3HHHHdi9ezdeeOEFrFu3DitXrgxpfSopKcHkyZNxxRVX4KqrrsLChQtx7733ok+fPpgyZco5vacA4IUXXsDcuXMxatQo3HPPPTh06BCmT5+OqKgoJCcnNyt/Y4qKijBlyhRcffXVuP7664Ov6TfeeAN2ux3z5s2D3W7Hd999h/vvvx/l5eV44oknAJz7e/3bb7/FlClTkJ6ejgcffBA1NTV49tlnMXLkSGzcuLHBSfRVV12FtLQ0PProo9i4cSNeffVVxMXFBV8HTfnhhx/gcrmQmpp6xv0KCwtx0UUXobi4GEuXLm12N6HTX69AbTHucDi
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3RUZf4G8OdOb5lJ7z2UEEqAUEWqSBE76lpRdHVXYF1l3XXd365tC7u2da1YUXft2FFERIoC0lF6CwTSeybJZPr9/REyJhIgQZJ3Jvf5nDPnwJ2byTM3U+73vk2SZVkGERERERGRACrRAYiIiIiISLlYkBARERERkTAsSIiIiIiISBgWJEREREREJAwLEiIiIiIiEoYFCRERERERCcOChIiIiIiIhGFBQkREREREwrAgISIiIiIiYViQ9ABHjhyBJEl49dVXRUfpsFWrVkGSJKxatarTP/vAAw9AkqSzH+pnKCsrwxVXXIGoqChIkoQnnnjiZz3HUBeMfyMKLT31PRUMn9cXXHABbr311sD/X331VUiShM2bNwvLRO3bvXs3NBoNdu7cKToKUZdiQRICWr4s2rv98Y9/FB0vaN10001tjpXVakVubi4ee+wxuFyus/q77rrrLixbtgz33nsv/vvf/2LatGnt7vfmm2/iiSee6PDjpqent3kOBoMBvXv3xu9//3tUV1efpfQ934QJE9ocx8jISAwfPhyvvPIK/H7/Wf1dr732GgYMGACTyYSUlBTMmjULxcXFHf75U73fJUnCd999d1bznkpZWRnuvvtuZGdnw2QywWw2Iy8vD3/7299QW1vbpb+7o++pYNXZ93p3Wbt2Lb788kvcc889Z+0xf/pZ+9PPrBanK3wuvPBCpKent9n208ezWq0YP348Pvvss07nnDlzJi644IJO/5xIOTk5mDFjBu677z7RUYi6lEZ0AOq4hx56CBkZGW22DRgwAGlpaWhqaoJWqxWUrPPGjRuHpqYm6HS6Lv09er0eL730EgCgtrYW77//Pu6++25s2rQJb7/99ln7PV9//TUuueQS3H333YFtffr0OeE5vvnmm9i5cyfuvPPODj/24MGD8bvf/Q4A4HQ6sWXLFjzxxBNYvXo1Nm7ceNaeQ0+XnJyMBQsWAAAqKirw+uuv45ZbbsH+/fvxz3/+86z8jg8//BA33XQTxo8fj3nz5qG8vByLFy/G/v37kZiY2KnHau/9DgC9evU6K1lPZ9OmTbjgggvQ0NCA66+/Hnl5eQCAzZs345///CfWrFmDL7/8sst+f0ffU8HqZO910Z/XjzzyCM4777yz/jpq/Vnbmlqt/tmPff7552PWrFmQZRkFBQV47rnncNFFF2Hp0qWYOnVqhx7D4/Fg+fLlgc+AUPLrX/8aF1xwAQ4dOoSsrCzRcYi6BAuSEDJ9+nQMGzas3ftaX4UKBSqVqlsyazQaXH/99YH/z5kzByNHjsQ777yDxx9/vN2TRFmW4XQ6YTQaO/x7ysvLER4e3mbb2XqOSUlJbZ7DL3/5S1gsFjz66KM4cOAAevfu/bN/hxLYbLY2x/FXv/oV+vbti6effhp//etfz8oJ4ttvv43IyEh88cUXgb/9fffdB7fb3enHOtX7vavV1tbisssug1qtxrZt25Cdnd3m/r///e948cUXuzRDV76nzoTD4YDJZPrZj/PTVoPuVF5ejs8++wwLFy4864/908/as6lPnz5tHnvmzJnIycnBf/7znw4XJN988w3q6+sxY8aMLsl4tnm9Xvj9fuh0OkyePBkRERF47bXX8NBDD4mORtQl2GWrB2ivT/JNN90Ei8WCoqIiXHrppbBYLIiJicHdd98Nn8/X5ucfffRRnHPOOYiKioLRaEReXh4WL158wu+RJAnz5s3DRx99hAEDBkCv16N///744osvTti3qKgIt9xyCxITE6HX65GRkYHbb789cGLWXl/wb775BldeeSVSU1Oh1+uRkpKCu+66C01NTWfnQKH5hGbChAkAmo8b0Nwt6sILL8SyZcswbNgwGI1GPP/88wCA/Px8XHnllYiMjITJZMKoUaPadBVo6YIgyzKeeeaZQLeC9p7jhAkT8Nlnn6GgoCCw30+7J3RUfHw8gOaTgBY//PADbrrpJmRmZsJgMCA+Ph4333wzqqqq2vxsfX097rzzTqSnp0Ov1yM2Nhbnn38+tm7d2ma/DRs2YNq0abDZbDCZTBg/fjzWrl17QpZvv/0Ww4cPh8FgQFZWVuDYddR7772HvLw8GI1GREdH4/rrr0dRUVGbfTrzeu6olr9nY2MjKioqcP/990Or1aKiouKEfW+77TaEh4fD6XSe8jFVKhW8Xu8JV4W76or+22+/jby8PISFhcFqtWLgwIH4z3/+E7jf4/HgwQcfRO/evWEwGBAVFYVzzz0Xy5cvP+XjPv/88ygqKsLjjz9+QjECAHFxcfjzn//cZtuzzz6L/v37Q6/XIzExEXPnzj2hW9eECRMwYMAA7N69GxMnToTJZEJSUhIefvjhwD6deU+1eOaZZ5CZmQmj0YgRI0bgm2++wYQJEwLv9daP2/K+b9HeY7bk3LJlC8aNGweTyYQ//elPAICPP/4YM2bMCHy2ZWVl4a9//Wub1+Gp3usnG0Py9ddfY+zYsTCbzQgPD8cll1yCPXv2tNmnZWzWwYMHcdNNNyE8PBw2mw2zZ8+Gw+E44e/0U5999hm8Xi8mT5582n1ramowYsQIJCcnY9++fafdvzv169cP0dHROHToUId/5rPPPkNOTs5JP3Pz8/MhSRL+/e9/n3DfunXrIEkS3nrrrcC2oqIi3HzzzYiLiwt8F77yyittfs7tduO+++5DXl4ebDYbzGYzxo4di5UrV7bZr+U18eijj+KJJ55AVlYW9Ho9du/eDQDQarWYMGECPv744w4/X6JQwxaSEFJXV4fKyso226Kjo0+6v8/nw9SpUzFy5Eg8+uij+Oqrr/DYY48hKysLt99+e2C///znP7j44otx3XXXwe124+2338aVV16JJUuWnHA16dtvv8UHH3yAOXPmICwsDE8++SRmzpyJo0ePIioqCgBQXFyMESNGoLa2Frfddhuys7NRVFSExYsXw+FwnPTk7L333oPD4cDtt9+OqKgobNy4EU899RQKCwvx3nvvnelhO0HLl1hLXgDYt28frrnmGvzqV7/Crbfeir59+6KsrAznnHMOHA4H7rjjDkRFReG1117DxRdfjMWLF+Oyyy7DuHHj8N///hc33HBDoFvByfzf//0f6urqUFhYGPjSs1gsp83r8XgCf3en04lt27bh8ccfx7hx49p06Vm+fDny8/Mxe/ZsxMfHY9euXXjhhRewa9cufPfdd4GTul//+tdYvHgx5s2bh5ycHFRVVeHbb7/Fnj17MHToUADNJ0fTp09HXl4e7r//fqhUKixatAiTJk3CN998gxEjRgAAduzYgSlTpiAmJgYPPPAAvF4v7r//fsTFxXXob/Hqq69i9uzZGD58OBYsWICysjL85z//wdq1a7Ft27Y2V8g7+nrujPz8fKjVaoSHh+OGG27AQw89hHfeeQfz5s0L7ON2u7F48WLMnDnztFe2Z8+ejbfffhv33Xffz+4a0t77XZKkwOt2+fLluOaaa3DeeefhX//6FwBgz549WLt2LX77298CaD6BXbBgAX75y19ixIgRsNvt2Lx5M7Zu3Yrzzz//pL/7k08+gdFoxBVXXNGhrA888AAefPBBTJ48Gbfffjv27duH5557Dps2bcLatWvbtD7V1NRg2rRpuPzyy3HVVVdh8eLFuOeeezBw4EBMnz69U+8pAHjuuecwb948jB07FnfddReOHDmCSy+9FBEREUhOTu5Q/vZUVVVh+vTpuPrqq3H99dcHXtOvvvoqLBYL5s+fD4vFgq+//hr33Xcf7HY7HnnkEQCdf69/9dVXmD59OjIzM/HAAw+gqakJTz31FMaMGYOtW7eecBJ91VVXISMjAwsWLMDWrVvx0ksvITY2NvA6OJl169YhKioKaWlpp9yvsrIS559/Pqqrq7F69eoOdxP66esVaC7GrVZrh36
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAADjqElEQVR4nOzdd3hUVf4G8PdOn8lk0iuBhNA7GKpKUZFqb2tFwVVXYFVYt/hbF8uuy9pYu+iqWNa2YgVFRHqTjkAICSUkpPc+mXp/f4RMkQQSTObMZN7P8+RR7rR3Tu5M7veec8+RZFmWQUREREREJIBCdAAiIiIiIgpeLEiIiIiIiEgYFiRERERERCQMCxIiIiIiIhKGBQkREREREQnDgoSIiIiIiIRhQUJERERERMKwICEiIiIiImFYkBARERERkTAsSLqAkydPQpIkvPvuu6KjtNmGDRsgSRI2bNjQ7sc+/vjjkCSp40P9CsXFxbjhhhsQFRUFSZLwwgsv/Kr3GOj88XdEgaWrfqb84ft6xowZuOeee1z/fvfddyFJEnbv3i0sE7Xs8OHDUKlUOHTokOgoRJ2KBUkAaP5j0dLPX/7yF9Hx/NZdd93l1VYmkwnDhg3D888/D4vF0qGvtWDBAqxevRqPPPIIPvjgA0ybNq3F+3300Ud44YUX2vy8KSkpXu9Bp9OhT58++OMf/4iKiooOSt/1TZo0yasdIyMjMWrUKLzzzjtwOp0d+lrvvfceBg8eDIPBgO7du2PWrFkoKCho8+PP9nmXJAk//fRTh+Y9m+LiYjz88MPo378/DAYDQkJCkJaWhn/84x+oqqrq1Ndu62fKX7X3s+4rW7duxQ8//IA///nPHfacv/yu/eV3VrNzFT5XXHEFUlJSvLb98vlMJhMmTpyIb7/9tt05r7/+esyYMaPdjxNp4MCBmDlzJhYtWiQ6ClGnUokOQG335JNPomfPnl7bBg8ejOTkZJjNZqjVakHJ2m/ChAkwm83QaDSd+jparRZvvfUWAKCqqgqff/45Hn74YezatQuffPJJh73OunXrcPXVV+Phhx92bevbt+8Z7/Gjjz7CoUOH8NBDD7X5uYcPH44//OEPAIDGxkbs2bMHL7zwAjZu3IidO3d22Hvo6pKSkrB48WIAQGlpKd5//33cfffdyMrKwr/+9a8OeY0vv/wSd911FyZOnIj58+ejpKQEy5cvR1ZWFhITE9v1XC193gGgd+/eHZL1XHbt2oUZM2agrq4Ot99+O9LS0gAAu3fvxr/+9S9s2rQJP/zwQ6e9fls/U/6qtc+66O/rZ599FpdddlmH70ee37WelErlr37uyy+/HLNmzYIsy8jJycHrr7+OK6+8EqtWrcLUqVPb9Bw2mw1r1qxxfQcEkt/97neYMWMGjh8/jl69eomOQ9QpWJAEkOnTp2PkyJEt3uZ5FioQKBQKn2RWqVS4/fbbXf+eO3cuxowZg08//RRLlixp8SBRlmU0NjZCr9e3+XVKSkoQHh7uta2j3mO3bt283sNvf/tbGI1GPPfcczh69Cj69Onzq18jGISFhXm143333Yd+/frhlVdewd///vcOOUD85JNPEBkZie+//971u1+0aBGsVmu7n+tsn/fOVlVVhWuvvRZKpRL79u1D//79vW5/6qmn8J///KdTM3TmZ+p8NDQ0wGAw/Orn+WWvgS+VlJTg22+/xdKlSzv8uX/5XduR+vbt6/Xc119/PQYOHIgXX3yxzQXJ5s2bUVtbi5kzZ3ZKxo5mt9vhdDqh0WgwefJkRERE4L333sOTTz4pOhpRp+CQrS6gpTHJd911F4xGI/Lz83HNNdfAaDQiJiYGDz/8MBwOh9fjn3vuOVx44YWIioqCXq9HWloali9ffsbrSJKE+fPn46uvvsLgwYOh1WoxaNAgfP/992fcNz8/H3fffTcSExOh1WrRs2dP3H///a4Ds5bGgm/evBk33ngjevToAa1Wi+7du2PBggUwm80d01BoOqCZNGkSgKZ2A5qGRV1xxRVYvXo1Ro4cCb1ejzfeeAMAcOLECdx4442IjIyEwWDA2LFjvYYKNA9BkGUZr776qmtYQUvvcdKkSfj222+Rk5Pjut8vhye0VXx8PICmg4BmBw4cwF133YXU1FTodDrEx8djzpw5KC8v93psbW0tHnroIaSkpECr1SI2NhaXX3459u7d63W/HTt2YNq0aQgLC4PBYMDEiROxdevWM7Js2bIFo0aNgk6nQ69evVxt11afffYZ0tLSoNfrER0djdtvvx35+fle92nP/txWzb/P+vp6lJaW4rHHHoNarUZpaekZ97333nsRHh6OxsbGsz6nQqGA3W4/46xwZ53R/+STT5CWlobQ0FCYTCYMGTIEL774out2m82GJ554An369IFOp0NUVBQuvvhirFmz5qzP+8YbbyA/Px9Lliw5oxgBgLi4ODz66KNe21577TUMGjQIWq0WiYmJmDdv3hnDuiZNmoTBgwfj8OHDuOSSS2AwGNCtWzc888wzrvu05zPV7NVXX0Vqair0ej1Gjx6NzZs3Y9KkSa7PuufzNn/um7X0nM059+zZgwkTJsBgMOD//u//AABff/01Zs6c6fpu69WrF/7+97977Ydn+6y3dg3JunXrMH78eISEhCA8PBxXX301MjIyvO7TfG3WsWPHcNdddyE8PBxhYWGYPXs2Ghoazvg9/dK3334Lu92OyZMnn/O+lZWVGD16NJKSkpCZmXnO+/vSgAEDEB0djePHj7f5Md9++y0GDhzY6nfuiRMnIEkS/v3vf59x27Zt2yBJEj7++GPXtvz8fMyZMwdxcXGuv4XvvPOO1+OsVisWLVqEtLQ0hIWFISQkBOPHj8f69eu97te8Tzz33HN44YUX0KtXL2i1Whw+fBgAoFarMWnSJHz99ddtfr9EgYY9JAGkuroaZWVlXtuio6Nbvb/D4cDUqVMxZswYPPfcc/jxxx/x/PPPo1evXrj//vtd93vxxRdx1VVX4bbbboPVasUnn3yCG2+8EStXrjzjbNKWLVvwxRdfYO7cuQgNDcVLL72E66+/Hrm5uYiKigIAFBQUYPTo0aiqqsK9996L/v37Iz8/H8uXL0dDQ0OrB2efffYZGhoacP/99yMqKgo7d+7Eyy+/jLy8PHz22Wfn22xnaP4j1pwXADIzM3HLLbfgvvvuwz333IN+/fqhuLgYF154IRoaGvDAAw8gKioK7733Hq666iosX74c1157LSZMmIAPPvgAd9xxh2tYQWv++te/orq6Gnl5ea4/ekaj8Zx5bTab6/fe2NiIffv2YcmSJZgwYYLXkJ41a9bgxIkTmD17NuLj45Geno4333wT6enp+Omnn1wHdb/73e+wfPlyzJ8/HwMHDkR5eTm2bNmCjIwMXHDBBQCaDo6mT5+OtLQ0PPbYY1AoFFi2bBkuvfRSbN68GaNHjwYAHDx4EFOmTEFMTAwef/xx2O12PPbYY4iLi2vT7+Ldd9/F7NmzMWrUKCxevBjFxcV48cUXsXXrVuzbt8/rDHlb9+f2OHHiBJRKJcLDw3HHHXfgySefxKeffor58+e77mO1WrF8+XJcf/315zyzPXv2bHzyySdYtGjRrx4a0tLnXZIk1367Zs0a3HLLLbjsssvw9NNPAwAyMjKwdetWPPjggwCaDmAXL16M3/72txg9ejRqamqwe/du7N27F5dffnmrr/3NN99Ar9fjhhtuaFPWxx9/HE888QQmT56M+++/H5mZmXj99dexa9cubN261av3qbKyEtOmTcN1112Hm266CcuXL8ef//xnDBkyBNOnT2/XZwoAXn/9dcyfPx/jx4/HggULcPLkSVxzzTWIiIhAUlJSm/K3pLy8HNOnT8fNN9+M22+/3bVPv/vuuzAajVi4cCGMRiPWrVuHRYsWoaamBs8++yyA9n/Wf/zxR0yfPh2pqal4/PHHYTab8fLLL+Oiiy7C3r17zziIvummm9CzZ08sXrwYe/fuxVtvvYXY2FjXftCabdu2ISoqCsnJyWe9X1lZGS6//HJUVFRg48aNbR4m9Mv9FWgqxk0mU5se31bV1dWorKx
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAyQAAAMKCAYAAAB0rIBCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAADY4ElEQVR4nOzdd3hUVf4G8PdOT+8NQgst9BJAECkKUkQRRV0rCK4NsODquu5v17brsouKvResiIoVFektgNTQQ+8ljfQ29fz+CBlmIIFEmTnJnPfzPHnEmcnknZM7k/u9p2lCCAEiIiIiIiIJdLIDEBERERGRuliQEBERERGRNCxIiIiIiIhIGhYkREREREQkDQsSIiIiIiKShgUJERERERFJw4KEiIiIiIikYUFCRERERETSsCAhIiIiIiJpWJAEgEOHDkHTNHz00Ueyo9TZsmXLoGkali1bVu/vffrpp6Fp2sUP9QdkZ2fjhhtuQExMDDRNw8svv/yHXmNj1xB/R9S4BOp7qiF8Xl911VW4++673f//0UcfQdM0bNiwQVomqtnOnTthMBiwfft22VGIfIoFSSNQ/ceipq+//e1vsuM1WHfeeadXW4WHh6Nbt2548cUXYbVaL+rPmjp1KubPn48nnngCn376KUaMGFHj42bNmoWXX365zs/bsmVLr9dgsVjQtm1bPPbYY8jPz79I6QPf4MGDvdoxOjoavXv3xocffgiXy3VRf9bHH3+Mzp07Izg4GM2aNcO4ceNw4sSJOn//+d7vmqbht99+u6h5zyc7OxuPPvooUlNTERwcjJCQEKSlpeHf//43CgsLffqz6/qeaqjq+173l1WrVmHBggV4/PHHL9pznv1Ze/ZnVrULFT5XX301WrZs6XXb2c8XHh6OQYMG4eeff653zrFjx+Kqq66q9/fJ1LFjR4waNQpPPvmk7ChEPmWQHYDq7tlnn0WrVq28buvcuTNatGiBiooKGI1GScnqb+DAgaioqIDJZPLpzzGbzXj//fcBAIWFhfjmm2/w6KOPYv369Zg9e/ZF+zlLlizBtddei0cffdR9W7t27c55jbNmzcL27dvx8MMP1/m5u3fvjr/85S8AgMrKSmzcuBEvv/wyli9fjnXr1l201xDokpOTMW3aNABAbm4uPvnkE9x1113Ys2cP/vvf/16Un/Hdd9/hzjvvxKBBgzBlyhTk5ORgzpw52LNnD5o0aVKv56rp/Q4Abdq0uShZL2T9+vW46qqrUFpaittvvx1paWkAgA0bNuC///0vVqxYgQULFvjs59f1PdVQ1fZel/15/fzzz2PIkCEX/Tjy/Kz1pNfr//BzX3nllRg3bhyEEDh8+DDeeustXHPNNZg3bx6GDx9ep+ew2+1YuHCh+zOgMbnvvvtw1VVXYf/+/WjdurXsOEQ+wYKkERk5ciR69epV432eV6EaA51O55fMBoMBt99+u/v/J02ahEsuuQRffvklZsyYUeNJohAClZWVCAoKqvPPycnJQWRkpNdtF+s1Nm3a1Os1/PnPf0ZoaCheeOEF7N27F23btv3DP0MFERERXu147733on379nj99dfxr3/966KcIM6ePRvR0dH49ddf3b/7J598Ejabrd7Pdb73u68VFhbiuuuug16vR0ZGBlJTU73uf+655/Dee+/5NIMv31O/R3l5OYKDg//w85zda+BPOTk5+Pnnn/H2229f9Oc++7P2YmrXrp3Xc48dOxYdO3bEK6+8UueCZOXKlSgpKcGoUaN8kvFiczgccLlcMJlMGDp0KKKiovDxxx/j2WeflR2NyCc4ZCsA1DQm+c4770RoaCiOHz+OMWPGIDQ0FHFxcXj00UfhdDq9vv+FF17ApZdeipiYGAQFBSEtLQ1z5sw55+domoYpU6bg+++/R+fOnWE2m9GpUyf8+uuv5zz2+PHjuOuuu9CkSROYzWa0atUK999/v/vErKax4CtXrsSNN96I5s2bw2w2o1mzZpg6dSoqKiouTkOh6oRm8ODBAKraDagaFnX11Vdj/vz56NWrF4KCgvDOO+8AAA4cOIAbb7wR0dHRCA4ORt++fb2GClQPQRBC4I033nAPK6jpNQ4ePBg///wzDh8+7H7c2cMT6ioxMRFA1UlAta1bt+LOO+9ESkoKLBYLEhMTMXHiRJw6dcrre0tKSvDwww+jZcuWMJvNiI+Px5VXXolNmzZ5PW7t2rUYMWIEIiIiEBwcjEGDBmHVqlXnZElPT0fv3r1hsVjQunVrd9vV1ddff420tDQEBQUhNjYWt99+O44fP+71mPocz3VV/fssKytDbm4unnrqKRiNRuTm5p7z2HvuuQeRkZGorKw873PqdDo4HI5zrgr76or+7NmzkZaWhrCwMISHh6NLly545ZVX3Pfb7XY888wzaNu2LSwWC2JiYnDZZZdh4cKF533ed955B8ePH8eMGTPOKUYAICEhAf/4xz+8bnvzzTfRqVMnmM1mNGnSBJMnTz5nWNfgwYPRuXNn7Ny5E5dffjmCg4PRtGlTTJ8+3f2Y+rynqr3xxhtISUlBUFAQ+vTpg5UrV2Lw4MHu97rn81a/76vV9JzVOTdu3IiBAwciODgYf//73wEAP/zwA0aNGuX+bGvdujX+9a9/eR2H53uv1zaHZMmSJRgwYABCQkIQGRmJa6+9FpmZmV6PqZ6btW/fPtx5552IjIxEREQEJkyYgPLy8nN+T2f7+eef4XA4MHTo0As+tqCgAH369EFycjJ27959wcf7U4cOHRAbG4v9+/fX+Xt+/vlndOzYsdbP3AMHDkDTNLz00kvn3Ld69WpomoYvvvjCfdvx48cxceJEJCQkuP8Wfvjhh17fZ7PZ8OSTTyItLQ0REREICQnBgAEDsHTpUq/HVR8TL7zwAl5++WW0bt0aZrMZO3fuBAAYjUYMHjwYP/zwQ51fL1Fjwx6SRqSoqAh5eXlet8XGxtb6eKfTieHDh+OSSy7BCy+8gEWLFuHFF19E69atcf/997sf98orr2D06NG47bbbYLPZMHv2bNx444346aefzrmalJ6ejm+//RaTJk1CWFgYXn31VYwdOxZHjhxBTEwMAODEiRPo06cPCgsLcc899yA1NRXHjx/HnDlzUF5eXuvJ2ddff43y8nLcf//9iImJwbp16/Daa6/h2LFj+Prrr39vs52j+o9YdV4A2L17N2655Rbce++9uPvuu9G+fXtkZ2fj0ksvRXl5OR588EHExMTg448/xujRozFnzhxcd911GDhwID799FPccccd7mEFtfm///s/FBUV4dixY+4/eqGhoRfMa7fb3b/3yspKZGRkYMaMGRg4cKDXkJ6FCxfiwIEDmDBhAhITE7Fjxw68++672LFjB3777Tf3Sd19992HOXPmYMqUKejYsSNOnTqF9PR0ZGZmomfPngCqTo5GjhyJtLQ0PPXUU9DpdJg5cyauuOIKrFy5En369AEAbNu2DcOGDUNcXByefvppOBwOPPXUU0hISKjT7+Kjjz7ChAkT0Lt3b0ybNg3Z2dl45ZVXsGrVKmRkZHhdIa/r8VwfBw4cgF6vR2RkJO644w48++yz+PLLLzFlyhT3Y2w2G+bMmYOxY8de8Mr2hAkTMHv2bDz55JN/eGhITe93TdPcx+3ChQtxyy23YMiQIfjf//4HAMjMzMSqVavw0EMPAag6gZ02bRr+/Oc/o0+fPiguLsaGDRuwadMmXHnllbX+7B9//BFBQUG44YYb6pT16aefxjPPPIOhQ4fi/vvvx+7du/HWW29h/fr1WLVqlVfvU0FBAUaMGIHrr78eN910E+bMmYPHH38cXbp0wciRI+v1ngKAt956C1OmTMGAAQMwdepUHDp0CGPGjEFUVBSSk5PrlL8mp06dwsiRI3HzzTfj9ttvdx/TH330EUJDQ/HII48gNDQUS5YswZNPPoni4mI8//zzAOr/Xl+0aBFGjhyJlJQUPP3006ioqMBrr72G/v37Y9OmTeecRN90001o1aoVpk2bhk2bNuH9999HfHy8+ziozerVqxETE4MWLVqc93F5eXm48sorkZ+fj+XLl9d5mNDZxytQVYyHh4f
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy0AAAMKCAYAAACIuauIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3wUdf4/8Nds77vZ9EYKLRB6qCJNkCKgKOpZsZ3eKd6d+vXOn1dsVzx7u1PPU1HvbCeeeoKIKFJFOkrvSSC9b9r2+f0RsiYSSFDCZ9h5PR+PPJTZyea1s23e82mSLMsyiIiIiIiIFEojOgAREREREdHJsGghIiIiIiJFY9FCRERERESKxqKFiIiIiIgUjUULEREREREpGosWIiIiIiJSNBYtRERERESkaCxaiIiIiIhI0Vi0EBERERGRorFoiQL5+fmQJAmvvfaa6ChdtmLFCkiShBUrVpzy7z7wwAOQJOn0h/oRysrKcOmllyI2NhaSJOHpp5/+UY/xbKfE54jOLtH6nlLC5/UFF1yAm2++OfLv1157DZIkYdOmTcIyUcd27doFnU6HHTt2iI5CJByLlrNA6xdKRz//7//9P9HxFOv6669vd6wcDgcGDx6MJ554Aj6f77T+rTvvvBNLly7Fvffei3/961+YPn16h/u99dZbePrpp7t8v5mZme0eg8lkQu/evfHrX/8a1dXVpyl99Js4cWK74+h2uzFixAi8+uqrCIfDp/Vvvf766xgwYAAsFgvS09Mxb948FBcXd/n3T/Z+lyQJX3/99WnNezJlZWW4++67kZOTA4vFAqvViry8PPzpT39CbW1tt/7trr6nlOpU3+tnytq1a/HZZ5/hnnvuOW33+f3P2u9/ZrXqrDiaNWsWMjMz2237/v05HA5MmDABixcvPuWcc+fOxQUXXHDKvydS//79MXPmTNx3332ioxAJpxMdgLruoYceQlZWVrttAwYMQEZGBpqbm6HX6wUlO3Xjx49Hc3MzDAZDt/4do9GIl19+GQBQW1uL999/H3fffTc2btyId95557T9neXLl+Oiiy7C3XffHdnWp0+f4x7jW2+9hR07duCOO+7o8n0PGTIE//d//wcA8Hq92Lx5M55++mmsXLkSGzZsOG2PIdqlpaXh4YcfBgBUVFTgjTfewE033YR9+/bhr3/962n5Gx988AGuv/56TJgwAbfffjvKy8uxcOFC7Nu3DykpKad0Xx293wGgV69epyVrZzZu3IgLLrgADQ0NuOaaa5CXlwcA2LRpE/76179i1apV+Oyzz7rt73f1PaVUJ3qvi/68fuyxxzB58uTT/jpq+1nbllar/dH3ff7552PevHmQZRkFBQV44YUXMHv2bCxZsgTTpk3r0n0EAgEsW7Ys8hlwNvn5z3+OCy64AAcPHkTPnj1FxyEShkXLWWTGjBkYPnx4h7e1vZp1NtBoNGcks06nwzXXXBP592233YZRo0bh3XffxZNPPtnhiaQsy/B6vTCbzV3+O+Xl5XC5XO22na7HmJqa2u4x/PSnP4XNZsPjjz+O/fv3o3fv3j/6b6iB0+lsdxx/9rOfoW/fvvjb3/6GP/7xj6flJPKdd96B2+3Gp59+Gnnu77vvPvj9/lO+r5O937tbbW0tLr74Ymi1WmzduhU5OTntbv/zn/+Mf/7zn92aoTvfUz9EU1MTLBbLj76f77c+nEnl5eVYvHgxXnzxxdN+39//rD2d+vTp0+6+586di/79++OZZ57pctGyevVq1NfXY+bMmd2S8XQLBoMIh8MwGAyYMmUKYmJi8Prrr+Ohhx4SHY1IGHYPiwId9ZG+/vrrYbPZUFRUhDlz5sBmsyE+Ph533303QqFQu99//PHHcc455yA2NhZmsxl5eXlYuHDhcX9HkiTcfvvt+PDDDzFgwAAYjUbk5ubi008/PW7foqIi3HTTTUhJSYHRaERWVhZuvfXWyMlbR33TV69ejcsuuww9evSA0WhEeno67rzzTjQ3N5+eA4WWk56JEycCaDluQEsXrFmzZmHp0qUYPnw4zGYz/vGPfwAADh06hMsuuwxutxsWiwWjR49u1y2htbuDLMv4+9//HunC0NFjnDhxIhYvXoyCgoLIft/vCtFVSUlJAFpOFFp9++23uP7665GdnQ2TyYSkpCTceOONqKqqave79fX1uOOOO5CZmQmj0YiEhAScf/752LJlS7v91q9fj+nTp8PpdMJisWDChAlYu3btcVnWrFmDESNGwGQyoWfPnpFj11Xvvfce8vLyYDabERcXh2uuuQZFRUXt9jmV13NXtT6fjY2NqKiowP333w+9Xo+Kiorj9r3lllvgcrng9XpPep8ajQbBYPC4q8vd1TLwzjvvIC8vD3a7HQ6HAwMHDsQzzzwTuT0QCODBBx9E7969YTKZEBsbi3PPPRfLli076f3+4x//QFFREZ588snjChYASExMxO9///t2255//nnk5ubCaDQiJSUF8+fPP64L2cSJEzFgwADs2rULkyZNgsViQWpqKh599NHIPqfynmr197//HdnZ2TCbzRg5ciRWr16NiRMnRt7rbe+39X3fqqP7bM25efNmjB8/HhaLBb/97W8BAB999BFmzpwZ+Wzr2bMn/vjHP7Z7HZ7svX6iMS3Lly/HuHHjYLVa4XK5cNFFF2H37t3t9mkdK3bgwAFcf/31cLlccDqduOGGG9DU1HTc8/R9ixcvRjAYxJQpUzrdt6amBiNHjkRaWhr27t3b6f5nUr9+/RAXF4eDBw92+XcWL16M/v37n/Az99ChQ5AkCU899dRxt3311VeQJAlvv/12ZFtRURFuvPFGJCYmRr4LX3311Xa/5/f7cd999yEvLw9OpxNWqxXjxo3Dl19+2W6/1tfE448/jqeffho9e/aE0WjErl27AAB6vR4TJ07ERx991OXHSxSN2NJyFqmrq0NlZWW7bXFxcSfcPxQKYdq0aRg1ahQef/xxfP7553jiiSfQs2dP3HrrrZH9nnnmGVx44YW4+uqr4ff78c477+Cyyy7DokWLjrsqtWbNGvz3v//FbbfdBrvdjmeffRZz585FYWEhYmNjAQDFxcUYOXIkamtrccsttyAnJwdFRUVYuHAhmpqaTngC995776GpqQm33norYmNjsWHDBjz33HM4evQo3nvvvR962I7T+kXXmhcA9u7diyuvvBI/+9nPcPPNN6Nv374oKyvDOeecg6amJvzyl79EbGwsXn/9dVx44YVYuHAhLr74YowfPx7/+te/cO2110a6MJzI7373O9TV1eHo0aORL0abzdZp3kAgEHnevV4vtm7diieffBLjx49v131o2bJlOHToEG644QYkJSVh586deOmll7Bz5058/fXXkRO/n//851i4cCFuv/129O/fH1VVVVizZg12796NYcOGAWg5gZoxYwby8vJw//33Q6PRYMGCBTjvvPOwevVqjBw5EgCwfft2TJ06FfHx8XjggQcQDAZx//33IzExsUvPxWuvvYYbbrgBI0aMwMMPP4yysjI888wzWLt2LbZu3druSntXX8+n4tChQ9BqtXC5XLj22mvx0EMP4d1338Xtt98e2cfv92PhwoWYO3dup1fIb7jhBrzzzju47777fnQ3lI7e75IkRV63y5Ytw5VXXonJkyfjkUceAQDs3r0ba9euxa9+9SsALSe5Dz/8MH76059i5MiR8Hg82LRpE7Zs2YLzzz//hH/7f//7H8xmMy699NIuZX3ggQfw4IMPYsqUKbj11luxd+9evPDCC9i4cSPWrl3brhWrpqYG06dPxyWXXILLL78cCxcuxD333IOBAwdixowZp/SeAoAXXngBt99+O8aNG4c777wT+fn5mDNnDmJiYpCWltal/B2pqqrCjBkzcMUVV+Caa66JvKZfe+012Gw23HXXXbDZbFi+fDnuu+8+eDwePPbYYwBO/b3++eefY8aMGcjOzsYDDzyA5uZmPPfccxg7diy2bNly3In25ZdfjqysLDz88MPYsmULXn75ZSQkJEReByfy1VdfITY2FhkZGSfdr7KyEueffz6qq6uxcuXKLndJ+v7rFWg
      "text/plain": [
       "<Figure size 1000x900 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "months_results = []\n",
    "months_costs = []\n",
    "months_overload = []\n",
    "months_nettos = []\n",
    "months_gen_energy = []\n",
    "months_gen_energy2 = []\n",
    "months_ess_rest = pd.DataFrame(30, index=pv_capacities, columns= ess_capacities)\n",
    "months_csv_data = {}\n",
    "for index, month_data in tqdm(enumerate(months_data), total=len(months_data), position=0, leave= True):\n",
    "    results = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "    costs = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "    overload_cnt = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "    nettos = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "    gen_energies = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "    gen_energies2 = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "    factory_demands = {}\n",
    "    buy_prices= {}\n",
    "    sell_prices = {}\n",
    "    pv_generates = {}\n",
    "    grid_need_powers = {}\n",
    "    times = {}\n",
    "    for pv_capacity in tqdm(pv_capacities, total=len(pv_capacities), desc=f'generating pv for month {index + 1}',position=1, leave=False):\n",
    "        factory_demands[pv_capacity] = {}\n",
    "        buy_prices[pv_capacity] = {}\n",
    "        sell_prices[pv_capacity] = {}\n",
    "        pv_generates[pv_capacity] = {}\n",
    "        grid_need_powers[pv_capacity] = {}\n",
    "        times[pv_capacity] = {}\n",
    "        for ess_capacity in ess_capacities:\n",
    "            (result, \n",
    "             overload, \n",
    "             cost, \n",
    "             netto,\n",
    "             gen_energy,\n",
    "             gen_energy2,\n",
    "             ess_rest,\n",
    "             factory_demand,\n",
    "             buy_price,\n",
    "             sell_price,\n",
    "             pv_generate,\n",
    "             grid_need_power,\n",
    "             time) = generate_data(pv_capacity=pv_capacity,\n",
    "                        pv_cost_per_kW=pv_cost_per_kW,\n",
    "                        pv_lifetime=pv_lifetime, \n",
    "                        pv_loss=pv_loss, \n",
    "                        ess_capacity=ess_capacity, \n",
    "                        ess_cost_per_kW=ess_cost_per_kW, \n",
    "                        ess_lifetime=ess_lifetime, \n",
    "                        ess_loss=ess_loss, \n",
    "                        grid_capacity=grid_capacity, \n",
    "                        grid_loss=grid_loss, \n",
    "                        sell_price=sell_price, \n",
    "                        time_interval=time_interval, \n",
    "                        data=month_data, \n",
    "                        days=months_days[index],\n",
    "                        storage=months_ess_rest.loc[pv_capacity, ess_capacity])\n",
    "            results.loc[pv_capacity,ess_capacity] = result\n",
    "            overload_cnt.loc[pv_capacity,ess_capacity] = overload\n",
    "            costs.loc[pv_capacity,ess_capacity] = cost\n",
    "            nettos.loc[pv_capacity,ess_capacity] = netto\n",
    "            gen_energies.loc[pv_capacity, ess_capacity] = gen_energy\n",
    "            gen_energies2.loc[pv_capacity, ess_capacity] = gen_energy2\n",
    "            months_ess_rest.loc[pv_capacity, ess_capacity] = ess_rest\n",
    "            factory_demands[pv_capacity][ess_capacity] = factory_demand\n",
    "            buy_prices[pv_capacity][ess_capacity] = buy_price\n",
    "            sell_prices[pv_capacity][ess_capacity] = sell_price\n",
    "            pv_generates[pv_capacity][ess_capacity] = pv_generate\n",
    "            grid_need_powers[pv_capacity][ess_capacity] = grid_need_power\n",
    "            times[pv_capacity][ess_capacity] = time\n",
    "    months_csv_data[index] = {\"factory_demand\": factory_demands, \"buy_price\": buy_prices, \"sell_price\": sell_prices, \"pv_generate\": pv_generates, \"grid_need_power\": grid_need_powers, \"time\": times}\n",
    "    months_results.append(results)\n",
    "    months_costs.append(costs)\n",
    "    months_overload.append(overload_cnt)\n",
    "    months_nettos.append(nettos)\n",
    "    months_gen_energy.append(gen_energies)\n",
    "    months_gen_energy2.append(gen_energies2)\n",
    "    draw_results(results=results, \n",
    "                    filename=f'plots/pv-{pv_capacity}-ess-{ess_capacity}-month-{index+1}-benefit.png',\n",
    "                    title_benefit=title_benefit,\n",
    "                    annot_benefit=annot_benefit,\n",
    "                    figure_size=figure_size)\n",
    "    draw_overload(overload_cnt=overload_cnt, \n",
    "                    filename=f'plots/pv-{pv_capacity}-ess-{ess_capacity}-month-{index+1}-unmet.png',\n",
    "                    title_unmet=title_unmet,\n",
    "                    annot_unmet=annot_unmet,\n",
    "                    figure_size=figure_size,\n",
    "                    days=months_days[index],\n",
    "                    granularity=granularity)\n",
    "annual_result = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "annual_costs = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "annual_overload = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "annual_nettos = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "annual_gen = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "annual_gen2 = pd.DataFrame(index=pv_capacities, columns= ess_capacities)\n",
    "\n",
    "\n",
    "# get the yearly results\n",
    "for pv_capacity in pv_capacities:\n",
    "    for ess_capacity in ess_capacities:\n",
    "        results = 0\n",
    "        costs = 0\n",
    "        overload_cnt = 0\n",
    "        nettos = 0\n",
    "        gen = 0\n",
    "        gen2 = 0\n",
    "        for index, month_data in enumerate(months_data):\n",
    "            results += months_results[index].loc[pv_capacity,ess_capacity]\n",
    "            costs += months_costs[index].loc[pv_capacity,ess_capacity]\n",
    "            overload_cnt += months_overload[index].loc[pv_capacity, ess_capacity]\n",
    "            nettos += months_nettos[index].loc[pv_capacity, ess_capacity]\n",
    "            gen += months_gen_energy[index].loc[pv_capacity, ess_capacity]\n",
    "            gen2 += months_gen_energy[index].loc[pv_capacity, ess_capacity]\n",
    "        annual_result.loc[pv_capacity, ess_capacity] = results\n",
    "        annual_costs.loc[pv_capacity, ess_capacity] = costs\n",
    "        annual_overload.loc[pv_capacity, ess_capacity] = overload_cnt\n",
    "        annual_nettos.loc[pv_capacity, ess_capacity] = nettos\n",
    "        annual_gen.loc[pv_capacity, ess_capacity] = gen\n",
    "        annual_gen2.loc[pv_capacity, ess_capacity] = gen2\n",
    "\n",
    "draw_cost(costs=annual_costs,\n",
    "          filename='plots/annual_cost.png',\n",
    "        title_cost=title_cost,\n",
    "        annot_cost=annot_cost,\n",
    "        figure_size=figure_size)\n",
    "draw_results(results=annual_result,\n",
    "                filename='plots/annual_benefit.png',\n",
    "                title_benefit=title_benefit,\n",
    "                annot_benefit=annot_benefit,\n",
    "                figure_size=figure_size)\n",
    "draw_overload(overload_cnt=annual_overload,\n",
    "                filename='plots/annual_unmet.png',\n",
    "                title_unmet=title_unmet,\n",
    "                annot_unmet=annot_unmet,\n",
    "                figure_size=figure_size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "def collapse_months_csv_data(months_csv_data, column_name,pv_capacies, ess_capacities):\n",
    "    data = {}\n",
    "    for pv_capacity in pv_capacities:\n",
    "        data[pv_capacity] = {}\n",
    "        for ess_capacity in ess_capacities:\n",
    "            annual_data = []\n",
    "            for index, month_data in enumerate(months_data):\n",
    "                annual_data.extend(months_csv_data[index][column_name][pv_capacity][ess_capacity])\n",
    "                # months_csv_data[index][column_name][pv_capacity][ess_capacity] = months_csv_data[index][column_name][pv_capacity][ess_capacity].tolist()\n",
    "\n",
    "            data[pv_capacity][ess_capacity] = annual_data\n",
    "    return data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [],
   "source": [
    "annual_pv_gen = collapse_months_csv_data(months_csv_data, \"pv_generate\", pv_capacities, ess_capacities)\n",
    "annual_time = collapse_months_csv_data(months_csv_data, \"time\", pv_capacities, ess_capacities)\n",
    "annual_buy_price = collapse_months_csv_data(months_csv_data, \"buy_price\",pv_capacities, ess_capacities)\n",
    "annual_sell_price = collapse_months_csv_data(months_csv_data, \"sell_price\", pv_capacities, ess_capacities)\n",
    "annual_factory_demand = collapse_months_csv_data(months_csv_data, \"factory_demand\", pv_capacities, ess_capacities)\n",
    "annual_grid_need_power = collapse_months_csv_data(months_csv_data, \"grid_need_power\", pv_capacities, ess_capacities)\n",
    "\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "for pv_capacity in pv_capacities:\n",
    "    for ess_capacity in ess_capacities:\n",
    "        with open(f'data/annual_data-pv-{pv_capacity}-ess-{ess_capacity}.csv', 'w') as f:\n",
    "            f.write(\"date, time,pv_generate (kW),factory_demand (kW),buy_price (USD/MWh),sell_price (USD/MWh),grid_need_power (kW)\\n\")\n",
    "            start_date = datetime(2023, 1, 1, 0, 0, 0)\n",
    "            for i in range(len(annual_time[pv_capacity][ess_capacity])):\n",
    "                current_date = start_date + timedelta(hours=i)\n",
    "                formate_date = current_date.strftime(\"%Y-%m-%d\")\n",
    "                f.write(f\"{formate_date},{annual_time[pv_capacity][ess_capacity][i]},{int(annual_pv_gen[pv_capacity][ess_capacity][i])},{int(annual_factory_demand[pv_capacity][ess_capacity][i])},{int(annual_buy_price[pv_capacity][ess_capacity][i]*1000)},{int(annual_sell_price[pv_capacity][ess_capacity][i]*1000)},{int(annual_grid_need_power[pv_capacity][ess_capacity][i])} \\n\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_data(data, filename):\n",
    "    data.to_csv(filename+'.csv')\n",
    "    data.to_json(filename + '.json')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.isdir('data'):\n",
    "    os.makedirs('data')\n",
    "\n",
    "save_data(annual_result, f'data/{pv_begin}-{pv_end}-{pv_groups}-{ess_begin}-{ess_end}-{ess_groups}-results')\n",
    "save_data(annual_costs, f'data/{pv_begin}-{pv_end}-{pv_groups}-{ess_begin}-{ess_end}-{ess_groups}-costs')\n",
    "save_data(annual_overload, f'data/{pv_begin}-{pv_end}-{pv_groups}-{ess_begin}-{ess_end}-{ess_groups}-overload_cnt')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzEAAANXCAYAAAAFO9HHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOzdd3gc9Z0/8PfM9r7aVa+WLTe5F1wwxhgIpoeSEAIhlFQw4QKXcqRdwiXhflxCEhICl1wCCQmEkhAgVGNsg3HFNrh3WbLV21Ztn/n9sdbawjYuWP6Odt6v59GD2R2v3xpJq/nMt3wkVVVVEBERERERDRGy6ABEREREREQng0UMERERERENKSxiiIiIiIhoSGERQ0REREREQwqLGCIiIiIiGlJYxBARERER0ZDCIoaIiIiIiIYUFjFERERERDSksIghIiIiIqIhhUUMERERERENKSxiiIjOoBUrVuCHP/whAoHAoP0bP/3pT/HPf/5z0F6fiIhINBYxRERn0IoVK/CjH/2IRQwREdHHwCKGiIiIiIiGFBYxRERnyA9/+EN885vfBADU1tZCkiRIkoR9+/YBAP7yl79g2rRpsNls8Pl8uP7667F///4Br7Fr1y5ce+21KC0thdVqRWVlJa6//noEg0EAgCRJiEaj+NOf/pR7/VtuueVMfppERESDzig6ABGRXlxzzTXYuXMnnnrqKfziF79AYWEhAKCoqAg/+clP8P3vfx/XXXcdvvjFL6KzsxO//vWvce6552LDhg3wer1IJpNYsGABEokEvva1r6G0tBTNzc3417/+hUAgAI/HgyeeeAJf/OIXMWPGDHz5y18GAIwYMULkp01ERHTaSaqqqqJDEBHpxc9+9jN885vfRENDA4YNGwYAaGxsxIgRI3DffffhO9/5Tu7YzZs3Y8qUKfjRj36E73znO3j//fcxZcoUPPvss/jUpz51zH/D6XTiU5/6FB5//PFB/myIiIjE4HQyIiLB/vGPf0BRFFx33XXo6urKfZSWlmLkyJFYsmQJAMDj8QAAXn/9dfT19YmMTEREJBSnkxERCbZr1y6oqoqRI0ce9XmTyQQgu47mnnvuwYMPPoi//vWvmDt3Lq688kp87nOfyxU4REREesAihohIMEVRIEkSXn31VRgMhiOedzqduT///Oc/xy233IIXXngBb7zxBu666y7cf//9WLVqFSorK89kbCIiImFYxBARnUGSJB3x2IgRI6CqKmprazFq1KjjvsaECRMwYcIEfO9738OKFSswZ84cPProo/jxj398zH+DiIgon3BNDBHRGeRwOABgQLPLa665BgaDAT/60Y/w4b1WVFVFd3c3ACAUCiGdTg94fsKECZBlGYlEYsC/MZjNNImIiETjSAwR0Rk0bdo0AMB3v/tdXH/99TCZTLjiiivw4x//GPfeey/27duHq666Ci6XCw0NDXj++efx5S9/Gd/4xjfw1ltv4c4778SnP/1pjBo1Cul0Gk888QQMBgOuvfbaAf/Gm2++iQcffBDl5eWora3FzJkzRX3KREREpx23WCYiOsN+/OMf49FHH0VraysURcltt/yPf/wDv/jFL7BhwwYAQFVVFS644ALcddddGDVqFBoaGvDjH/8Yy5YtQ3NzM+x2OyZNmoTvfve7uOCCC3Kvv2PHDnz5y1/G2rVrEYvFcPPNN3O7ZSIiyissYoiIiIiIaEjhmhgiIiIiIhpSWMQQEREREdGQwiKGiIiIiIiGFBYxREREREQ0pLCIISIiIiKiIYVFDBERERERDSksYoiIiIiIaEgxig4wGK665YcAgMtv/A6MJrPYMBqyYtFf0NG8G1PPuRrVdZNEx9GMWF8Yrz/zc0iSjE/e/APRcYTavis+4P/XvfUIYpFuTJhzEzz+GkGpxNuy+in0tu/ByClXoKRKvz87eza9htaG91A1ai5qxswTHUeY5j2r0LDlTRRVjsfoqVcBALatbxQbimgIeuGR0aIjHNPLJu1muyy1Q3QETcjLIoaI6HQaOflKqJk0jGab6ChERESEPC1izv/kHQAAg9EkOIm2jD9rAVKT5sHp8omOoilmiw0zz/+s6BiaZLa6oChpyHJevlWcMLPFITqCJlQMn4miinGwWN2iowhVUj0ZvtJRMBg40k9EJEpeXpm4C4pFR9Akt7dIdARNMhiMKKvW7rCxSBPO/pzoCKQhVkcBrI4C0TGEM5qsMJqsomMQ0SCSTJLoCHQceVnEEBGdTm2N69EX6UZReT1cBRWi4xAREeleXhYx299fCgAYNWEuZINBbBgNaWnchr5IAMUVdRyVOUw6nUJr4zYAQNWIiYLTkBZ1t+1Ab/seONzFui5iAl37EA22wVVQAbevSnQcYULdTehp3w2HpwRFFeNExyEi0qW8LmLqxp3NIuYw+3auy+5OZrGziDlMKhnHunf+AUmSWcR8yLa1zyIRC6Fu4qVwestExyHBulu353Yn03MREw604MDuFSiqHM8ihihPyUZOJ9O6vCxiiOj06At3IRbpRiaTFB2FiIiIKIfNLomIiIiIaEjhSAwRERER0WEkE+/zax2/QkRERERENKSwiCEiIiIioiGF08mIiI6jbsKlyNSnYLY6RUchIqIzgLuTaV9eFjHzLvsSAEA25OWnd8rqp16AunFnw8XtlQcwma2Ydu61kCS+YX2Y0WyDyWKHJOt7q3KL3SM6giaUDZsOf+lo3Z+PosoJ8BbVwmi0io5CRKRbeXmVX1Ck32Z0H8XrZ5+PozEaTagaPkF0DE2adM4toiOQhthdhbC7CkXHEM5sccBscYiOQUSka3lZxBARnU4d+zciFu2Bv3Q0m34SEemAZOLsDK3LyyJm16Z3AQAj6mdBNuh7Gszh2vbvRKwvhKLSWjg9ftFxNCOdTqGjeTcAoLxmrOA0pEWdLVvQ274HVkeBrouYUHcTouFOOD2lcBXod8Q71HsAwc59sLuK4C8bLToOEZEu5eXuZFvWLcKWdYugKBnRUTRl7/Y1+GDlv9DTeUB0FE1JJeNYs+RprF36rOgomrNj/T+x8d0/IxpqFx2FNKCzZSv2bHwVPe27RUcRKtxzAI3bl6KrdZvoKEREupWXIzFEdHpEAq2IRbqRTsVFRyEiIjpjuDuZ9uXlSAwREREREeUvFjFERERERDSkcDoZEREREdFhuDuZ9nEkhoiIiIiIhhSOxBARHcfwcRchPToBq80rOgoREREhT4uYcy6+BQAgG/Ly0ztloyfNQ+3os+DxlYqOoikmkwWTz75SdAxNMhhMMBjNkKDvYXWbk32VAKC0Zgq8RcNhc/hERxGqsGIcXAUVMJntoqMQ0SDh7mTal5dX+YWlw0RH0CR/cZXoCJpkNJkxbNRU0TE0afK8L4qOQBricJfA4S4RHUM4i9UFi9UlOgYRka7lZRFDRHQ6dbZsRSLai4KSOl7EExERaUBeFjF7t60BAAwbNQ2ywSA4jXZ0tOxFIhaBr7gKDleB6Diakcmk0dW2DxIkFFeMEB2HNKhj/wfobd8Dk9Wp6yIm3NuMWLQHdlcRnB79TksNB1oQ6m6CzVkIX0md6DhENAgkA6eTaV1e7k62cfUr2Lj6FShKRnQUTdm9ZQXWvfMPdLc3iY6iKclEDCsX/QUr3/yr6Cias+v9f2HLqqcQDXWIjkIa0HFgE3aufwHdrTtERxEq1N2Ehi1vorN5s+goRES6lZdFDBGdHqGe/ejt2IN0KiY6ChEREVFOXk4nIyIiIiI6VTKnk2keR2KIiIiIiGhIYRFDRERERERDCqeTEREREREdRpI5nUzrWMQQER3HsLEXoLLubNgcftFRiIiICHlaxMy64AYAgGzIy0/vlI2acA6qRkyCr6hSdBRNMZksmDDjEkgS77rQ0TncxaIjaEJx1US4CyphdxeJjiKUv2wMHO4SmCwO0VGIiHQrL6/yS6tGiY6gSYWlw0RH0CSjyYwR9TNFx9CkaeffLjoCaYjLWw6Xt1x0DOGsdi+sdq/oGEQ0iCQDl41rXV4WMUREp1N3204k+gLwFtXC7tL3KAQREZEW5GWZ2bT7fTTtfh9KJiM6iqZ0tze
      "text/plain": [
       "<Figure size 1000x1000 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "draw_results(annual_result, 'plots/test.png', 'test', False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_roi(annual_costs, annual_nettos, 'plots/annual_roi.png',  title_roi, 365, annot_benefit, figure_size)\n"
   ]
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}