{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ADM1F: Steady State\n",
    "\n",
    "Here we run the steady state case and comparing it with the Matlab results. Make sure to compile `build/adm1f.cxx`. \n",
    "\n",
    "Author: Elchin Jafarov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Steady State Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import subprocess\n",
    "import sklearn.metrics as sklm\n",
    "import xlrd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# navigate to simulations folder\n",
    "os.chdir('../../simulations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/elchin/project/ADM1F_WM/build/adm1f\r\n"
     ]
    }
   ],
   "source": [
    "# check the path to the executable \n",
    "!echo $ADM1F_EXE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vliq [m3] is: 3400.000000\n",
      "Vgas [m3] is: 300.000000\n",
      "Reading parameters in file: params.dat\n",
      "Reading influent values in file: influent.dat\n",
      "Reading initial condition values in file: ic.dat\n",
      "Running as steady state problem.\n",
      "Solving.\n",
      "Done!\n"
     ]
    }
   ],
   "source": [
    "# running the executable in the cell\n",
    "!$ADM1F_EXE -steady"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove the output files\n",
    "!sh clean.sh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<subprocess.Popen at 0x7f84a4974040>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# or run using subprocess \n",
    "subprocess.Popen('$ADM1F_EXE -ts_monitor -steady', shell=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the run is successful then `indicator-***.out` should be saved in the `simulations` folder. Here take the last time step saved in the last indicator file (`indicator-062.out`). The accending numeration of the output files corresponds to the time iterations taken towards the steady state condition.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Comparison of the ADM1F (C++) with ADM1 (Matlab)\n",
    "\n",
    "The ADM1F runs much faster than the corresponding Matlab version. The main difference between C++ and the Matlab versions of the model is that ADM1F uses optimized solvers from the PETCS package to solve the corresponding mass balance equations. The ADM1F allows usage of the diffrent solvers. The ADM1(Matlab) is using ode45 nonstiff differential equation solver that cannot be changed. Below we benchmark ADM1F(C++) outputs with the ADM1(Matlab). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the output produced by ADM1(Matlab) from the xls file\n",
    "wb = xlrd.open_workbook('../docs/jupyter_notebook/out_sludge.xls')\n",
    "sheet = wb.sheet_by_index(1)\n",
    "results_matlab = [sheet.cell_value(4,i) for i in range(66)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "load the last file from the steady runs and compare it with the Matlab output. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_c=np.loadtxt('indicator-062.out', skiprows=2, unpack=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Matlab</th>\n",
       "      <th>C++</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10.116210</td>\n",
       "      <td>10.11620</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.529384</td>\n",
       "      <td>4.52938</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>83.095031</td>\n",
       "      <td>83.09500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>9.171951</td>\n",
       "      <td>9.17149</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>11.824647</td>\n",
       "      <td>11.82410</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      Matlab       C++\n",
       "0  10.116210  10.11620\n",
       "1   4.529384   4.52938\n",
       "2  83.095031  83.09500\n",
       "3   9.171951   9.17149\n",
       "4  11.824647  11.82410"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res = pd.DataFrame({\n",
    "    \"Matlab\": np.asarray(results_matlab),\n",
    "    \"C++\": np.asarray(results_c[:-1])})\n",
    "res.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAEbCAYAAAAMKCkgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABG90lEQVR4nO2dd5wkZZn4v0/35NnZHNllWXZZcmZJCpIFFEURFQ8VPZHjwk89c4YzB06PMx4ZUURJJ4rAESXvsoAL7MKyETbn3cmp+/n9UVXdNT0dZ6Z6+p15vp/PfKa7qrrqrZ6a93mfLKqKYRiGYURJbLgHYBiGYYx8TNgYhmEYkWPCxjAMw4gcEzaGYRhG5JiwMQzDMCLHhI1hGIYROSZsDKcRkVNFRP2f+4fx+pNL+MxjIvLzIo4J7mvB4EdaGqFrd5b72oUIja11uMdiFI8JG2OkcAjwoeCNiFzpT0gPZx4oIv/s73ullAuIyFoR+fwQjLUYLgCOK9O1cvFJYJ/wBvG4VESeEZEWEWkWkRdE5IsiMrZM45oBfKZM1zKGCBM2xrAiIjVDdKqtqrorY9tm4CQRmZOx/RPAm0N03UhQ1Z3AtlI/JyJVIiJDNIzdqrolY9stwM+AvwJnAIcD3wBOwxOQ2cZ0pYjcVOxFfa3uY7n2q+pmYE+x5zMqAxM2RlnxJ5JfichVIrINeEpEDhaRe/2V8lYR+b2ITA995jARedhfRbeKyBIROa2Iy+0A7gU+HjrX4cCBwB0Z45onIn8Skc0i0uav1s8Ljxtvlf/jwIyT4/4m+eNfLyIdIrJURD6e5dAqEblaRHb5Pz8WkZL/H/2J/BUR+ZiIrAK6gEYRGSci1/jfZ4uI/C1sjvP33+Lv7xSR1SLymQLX+gBwMXCxqn5bVRep6lpVvVdVzwX+t9TxG6MHEzbGcPBhQICTgU8BjwOv4JmNzgTGAH8KTb63Apv8/UcCVwLF+hKuBy4JnesTwB+BlozjxgD3AWcBRwB3AneJyIH+/guA9cC38Mw4M3Jcrw54ATgPz7R3NfA/InJGxnEX4/3/nQj8E3AZAzcN7Qv8A/B+f+xdeEJ2pj+Oo/C+40dEJBj3d4DD/P0HAP8IbChwnYuB11X1rmw7VXX3AMdvjAKqhnsAxqhkjap+DkBEvgUsUdUvBTtF5KPATmABsAhPo7hKVV/zD1lZwrXuB6qBM0TkcTxBdz6eUEuhqkuAJaFN3xWRdwEXAt9R1Z0ikgBafDNOVlR1A/Dj0KZrROR0PH9S2H+0CfiUesUJXxOR/YHPAj8p4d4CaoCPBCYv/3pHAlNUtcM/5hv+/XwE+BHed/qCqi7y979RxHXmA68VPMowsmCajTEcPB96fQzwNt881upHGK3z983zf/8EuE5EHhGRr4W0jYKoagK4GW/l/h5gm6o+mXmciDSKyI9EZJlv1mrFE3azS7kxEYn7Y3xJRHb457kgy3me1b5VcJ8BZuZzsoe/IxH5dWjX+gzfyjFAA7At43s9lPR3+ivgg75J8ioROaWY2yviGETk4ozrfhW4OGP8F4eO/3XG8ScDv844vqS/g1F5mGZjDAdtodcxPJNPtiivLQCqeqWI/A44FzgbuEJELlfVG4q83o3AS8Ac/3U2rgLO8cexAmgHfoOnNZTC54HPAZ8GXgZage8BU0s8TzaODL1uDr1uyzguhvfdnZzlHM0AqnqfiOyD952eAdwrIrerajb/UsDrwEFFjPMeYGHo/afwTHpfCm0LC8dv4n3/Ab/DN2OGtm0s4rpGBWPCxhhuXgA+ALyhqj25DlLVFXhC4L9F5FfApUBRwkZVV4jIIuAtwHtzHHYS8BtVvRNAROrwtIDXQ8d0A/EClzsJ+LOq3uKfR4D9gd0Zxx0vIhLSbk4ANqpqMzlQ1WLNhy8A04Ckqq7Oc77teNFlt4jIfcDvfSHeleMjtwK3icgF2fw2IjJeVXeragshn5iI7ATG5hq/qm4FtoaO78CLLizFXGpUOGZGM4abXwDjgD+IyPEiMldEzvQjqZpEpF5EfiFe8uQcETkeb0JfVuJ1zgWm5vG3vA68V0SOFpHDgN/iOfvDrAVOFpGZkjuJ83U8/9BJvrnv53gO/Ez2Av5LRA4QkQuBLwA/Le2WcvIQ8BRekMW5IrKviJwoIv8hIieD5ysTkfeIyHwROQjP1Lc6j6ABL7DiD8DvROQbInKsiOwjIueIyL14ZkrDyIoJG2NYUdWNwFuBJJ4zfymeAOryfxLABOAmYDlwN55/47MlXqc9Sx5OmM/ira6fwItKe9Z/HeabwN7AKnLnwHwHL6jhPrwIsDY8s1Amv8PTkhYC1+JFzQ2JsPG1pXcAj/jnXo4nKA4gbY7qAr6LFxTxFNAEvKuI834Iz0R4HvAonqnw+8Df8ExfhpEVsU6dhsuIyKl4k94U3yw0YhAvGXUNcKyqLi7ztRV4v6reUfDgYUC8pM+fq+qY4R6LURym2RgjhbUicvdwD2Ko8H0oS4d5GLeISMUJcD9i7dcFDzQqCtNsDKcRkXq8SCeANlXdNJzjGSpEZCZQ779dV8CXEsX19/Nf5g0yGA4qeWxGbkzYGIZhGJFjZjTDMAwjckZ1ns3kyZN1zpw5wz0MwzAMp3j++ee3q+qUUj4zqoXNnDlzWLy4rEE+hmEYziMixdTS64OZ0QzDMIzIMWFjGIZhRI4JG8MwDCNyTNgYhmEYkWPCxjAMw4gcEzaGYRhG5JiwMQzDMCLHhI1hGMYoQFW58p6l3LrwzWG5/qhO6jQMwxgtrNvZwU1Pr2VKUy3/cPzssl/fNBvDMIxRwM72bgB2tXUzHAWYTdgYhmGMAnb5wqY3qbR09Zb9+iZsDMMwRgF72ntSr3e39eQ5MhpM2BiGYYwCdvuaDaS1nHJiwsYwDGMUsCuk2ZiwMQzDMCJhT0fIjNZuZjTDMAwjAnaZGc0wDMOImt19zGim2RiGYRgRsLuPGc00G8MwDCMC+kajmWZjGIZhREAfM1qbaTaGYRjGEJNIKs2dFvpsGIZhREhzRw/hcmgW+mwYhmEMOUFwwMTGGsA0G8MwDCMCAuEya0I9VTGhvTtBV2+irGMwYWMYhjHCCYpwTmioYXxDNVB+U5oJG8MwjBFOoNmMb6hmQsPwmNJM2BiGYYxwdoc0m5SwKXObARM2hmEYI5wgQGBcfXXKjGaajWEYhjGk7DYzmmEYhhE1YTPa+EYLEDAMwzAiINBixoU1mzKXrDFhYxiGMcIJGqd5AQKBz8Y0G8MwDGMISYU+16c1m3K3GTBhYxiGMcIJ/DPjG6qZMEwla0zYGIZhjGB6E0laOnsRgbF11SkzmgUIGIZhGEPGnlCOTSwmjPfNaDtNszEMwzCGit2h4ADw/DbgCaFEUnN+bqgxYWMYhjGCCQIBxvlCpioeo6muClWvz025qChhIyLniMhyEVkpIl/Osr9WRP7g718oInP87XNEpENE/u7//LrsgzcMw6hAwsEBAcNRRaCqbFcqgIjEgV8AZwHrgedE5B5VXRY67BPALlXdT0QuAn4IfNDft0pVjyznmA3DMCqdcPWAgAkN1by5s7y5NpWk2RwHrFTV1araDdwGnJ9xzPnAzf7rO4AzRETKOEbDMAyn2JVhRgNS4c/lzLWpJGEzE1gXer/e35b1GFXtBfYAk/x9+4rIiyLyNxE5OddFROQyEVksIou3bds2dKM3DMOoQPZ0ZNNsAjPa6NRsBsMmYLaqHgV8FrhVRMZmO1BVr1HVBaq6YMqUKWUdpGEYRrkJN04LSHfrHJ2azQZg79D7Wf62rMeISBUwDtihql2qugNAVZ8HVgH7Rz5iwzCMCidfgMDOMhbjrCRh8xwwX0T2FZEa4CLgnoxj7gEu8V9fCDyiqioiU/wAA0RkLjAfWF2mcRuGYVQsgRltfEaAAJTXjFYx0Wiq2isi/wY8AMSBG1R1qYh8C1isqvcA1wO3iMhKYCeeQAJ4G/AtEekBksDlqrqz/HdhGIZRWYSLcAaMH4ZinBUjbABU9a/AXzO2fTP0uhN4f5bP3QncGfkADcMwHCN76HP582wqyYxmGIZhDDGBsBkX9tkMQ7dOEzaGYRgjlJ5EktauXmICTbVpQ5ZpNoZhGMaQEQ4OiMXS+e/hPBvV8hTjNGFjGIYxQtmdJTgAoL4mTm1VjO7eJB09ibKMxYSNYRjGCCVbjk1AuasImLAxDMMYoexq759jExAIoF1lSuw0YWMYhjFCyWVGg/IHCZiwMQzDGKFkqx4QMLHRzGiGYRjGEJCtCGdAuYtxmrAxDMMYoaSrB+Qxo7WZZmMYhmEMgnT1gDwBAqbZGIZhGINhd4cnSPJpNmZGMwzDMAZFKs+mvr9mE9RHswABwzAMY1DkS+ocb6HPhhE9L6/fwz1LNg73MIxhYk9HDzc9tYbWrt7hHkqk7M4TjVbuPJuK6mdjGOVAVfmnWxazcU8ncyY1cPis8cM9JKPM3LboTb5/32us2tbGt99z6HAPJxK6e5O0dSeoigljavtP9RMDn41FoxlGNCzd2MzGPZ0APLFi+zCPxhgOdvolWu5ZspGu3vIUoiw3QXDA+IZqRKTf/qa6KmICLV299CSSkY/HhI0x6nj41a2p10+s2DaMIzGGi6DS8Z6OHh59bWuBo90kFfacpVQNQCwmofbQ0Ws3JmyMUcfDr21JvX7hjd20d49su73Rn47utDZz5wsbhnEk0bE7TxHOgHJWETBh4zidPQnaRriTcyjZ0tzJS+v3UFsV44BpTXQnkjy3dtdwD2vAqOqId3JHQbiHy6OvbWVHa9cwjiYaAgGSLccmoJxtBkzYOEwiqbz7509y+n8+VrYy4a7ziG8yOWm/yZx+0FQAnnTYlPbLx1Zx+JUP8PRK8z2VQqcvbOqr4/QmlT+PwMjEtBktt2YzoYxVBEzYOMyzq3fw+pZWtjR3cfXDK4Z7OE7w8KueCe2Mg6Zx8n6TAXhy5Y7hHNKA6UkkufGpNSQVbn5m7XAPxynafTPau4/YC4C7Xhx5prR81QMCUrk2ZVismrBxmLtD/yC/ffYNVm9rHcbRVD6dPQme9DWAMw6aytH7TKC2Ksarm5rZ1uKeGeVvy7exvdWbJB59bVvZyo6MBAIz2ruO2IumuipeWr+HFVtahnlUQ8uuPAmdAeVsM2DCxlE6uhPc9/ImwDMJ9SaV79/32jCPqrJ5etV2OnuSHDZzHNPG1lFXHee4fSem9rnGHc+vB6A6LnQnktzrPw9GYYIAgfEN1Zx3+Ayg9ECBl9fv4fon15BI6pCPbyjIV4QzwAIEjII8+OoW2roTHLn3eH7ygSNoqInz4LItPLPKTZNQOXjID3k+w/fVgCeoAZ5yzOexs62bh1/bQkzgc28/AID/HYGmoKhI+Wxq4lxw9CzA+/6KFRydPQku/c1zfPsvy7jj+XWRjXMw7CnCjFbOKgImbBzl7he8Ve17j5rJ1LF1XH7KPAC+c+8ykhW60hpOVJVHAmFz4LTU9rcGfpsV21F153u75+8b6Ekop+w/hY+csA/11XGeW7uLdTvbh3toTtARChBYsM8EZk9sYHNzZ9GLtZufXsuWZs/0+j+Pr65I7SboU5OtCGdAOkDAzGhGFra3dvH4iu1UxSRlAvjkyXOZPraOpRubR6Szc7As3djM5uZOpo2t5dCZY1PbD54xlomNNWzc08ma7W3DOMLSuN03oV14zN401lZx9iGeADXtpjgCM1p9dRwR4YKjZwJwl7+Iy8eejh5++dgqABpr4qze1saDyzYPaBy3LnyT2xdHoxnt7ijssxlfxjYDJmwc5M9LNpJIeqvaSWNqAc8c8IWzPXPKVQ8st0TFDIKqAacfOK1P6Y5YTHjLvEkAqeCBSmfZxmaWbmxmXH11yiT4nqO8yfLuFzc4paENF509XnmW+po4ABcc5ZnS7ntlc8G8pWsfX82ejh6O33di6n/uV4+tKvl737C7g6/e/TJfvuvlSP5f8xXhDLA8GyMvwer1vf5qLOC9R83k0Jlj2dzcybWPrxmOoVUsQdWAM0P+moCTQqY0F7jTX32/+4i9qKv2JsuT9pvM5DG1rN7exkvr9wzn8Cqe3kSS7kQSEait8qbA2ZMaOHbOBDp6Etz/Sm4tZVtLF9c/6f1vffGcA/ngsbOZ2FjDkvV7SvaXBjlfiaSyYsvQR5IWU0FgQiUGCIjIv4jIUhFpF5G5/rYvi8gHohuekcmqba0sWb+HptoqzjxoWp99sZjw9XceDMCv/7aKLc2dwzHEiiOoGlBXHUv5aMKcNN/b9szqHfSWoSDhYOhJJFOLjQuPmZXaXhWPpXJG7jZTWl46e32txjehBbzPDxTIZ0r7+SMr6OhJcOZB0zhmnwnU18T5+FvmAPCrv60qaRzhmmzLNw9t2HVnT4KOngTVcaHR196yMT6k2UStERclbETkM8DXgWuAcPnQDcC/Df2wjFwEE825h01PrWrDnDB3Em8/eBodPQk+c9vfzZxG36oB2b6zWRMamDOpgZbOXl7aUNlawWPLt7GjrZv5U8dw+KxxffYFfoc/L9lYliq+rhL214R5x+EzqK2K8fSqHXznL8v6LTzW7Wzn1kVvIkLKfAbw0RPn0FgT54kV23m5SK2yozvRJwLytSEWNns60tUDslV8DqipijGmtopEUmnujHauKFazuRz4pKpeDYRH9AJwyJCPyshKMqmpVWtgo8/GN847mMljanlm9Q4uuWERLZ3l6VdRTkqJ/gmqBpx+4LScxwTazVMVbkoLwmwvPGZWv0nkkL3Gst/UMexo6y6pmvXzb+zi/lc2jRpfTxD2nLnwGFtXzbfOP4SqmHDdk2u4+LqFfZJ9f/rQ6/QklPceOZMDpjelto9rqOYfjp8NeBaFYnh61Xa6epNUx72/4fItzYO6p0wCE1q+sOeAXLk2T6/czqOvbR2y56JYYbMP8EqW7T1A/ZCMxCjI82/uYv2uDmaMq+OEfSflPG7viQ384Z9OYPrYOp5bu4sPX7+IPWXqMx41a7e3cdlvFjP/a3/lpB8+wuW3PM/PH1nBo8u3srWlv9kws2pALgK/zRMVHCSwo7WLh1/dSkw8/1wmIpLafveLhWt9bW3u5NO3vcj7fvU0l//2BT56wyI27u4Y8nFXGh2hHJtMPnjsbH5/2QlMaapl4ZqdnPezJ3j+jZ0s39zC3S9uoDou/PtZ+/f73CdOmkt1XPjrK5uKimp82Ne2g7/XUJvRdhURHBCQGSSwbmc7l968mH+4biEfv+k5Lv/t82wfgkKlxXbqXA0cDbyRsf0dwLJBj8JHRM4BrgbiwHWq+oOM/bXAb4BjgB3AB1V1rb/vK8AngATwKVV9YKjGVSnc5Wc4n3/kTGKx3KoxwLwpY7j98hP50LXPsmTdbj507bPc8onjUtFrrtHc2cPPH1nJjU+toSfhrbTW7+pg/a4O7l+adujWxGMgEBMQ3+IbrhqQixPnTiYm8OKbu2jr6qUxS2fD4eaeJRvpTSqnHTCFqTnu5fwj9+LHDyzn/5ZupqWzh6a6/pNNbyLJLc++wU/+73VaunqprYpR75uBzv7p41zx7kN439Ez+2hOyaTy7Ood3LroTdbt6uDCY2bxgQWzqK3K7Q+oVHKZ0QKOnTORe//fSfzbrS+yaO1OPvg/z7LPpAZU4UPHzWbviQ39PjN9XB0XHDWLPyxexzWPr+L7Fxye8/rhnK+PnDCH+17ezPbWbra1dDGlaWj+P4spwhkQCKStzZ384tGV/OyRFXT2JGms8XxaDyzdwnNrd/Ht8w/lnX6qxUAo9j/qKuDnItKA57M5UUQ+AnwR+McBXz2EiMSBXwBnAeuB50TkHlUNC7NPALtUdT8RuQj4IfBBETkYuAjPpLcX8JCI7K+qkbXgU1XauhP0JpIkkur9qNKbUCaPqc26asqkN5FkxdZWxjdUM3lMLdXx3IpmV2+Ce1/yVqvZVrXZ2HtiA7dffiIXX7uQZZuaueiaZ/ndpcfnnKjKgapnG66vjlNTVVixTiSV2557k5/83+vsaOtGBN5/zCz+/az9ae3qZenGPbyyoZlXNuxh2cZmWrKErYp4k0Q+xjVUc9is8SxZt5tFa3Zy2oG5taCB0JNI0t6VoDuRpMf/6e71IqLmTRmT164ecEcotyYXsyY0cPy+E1m4Zif3v7KZ9y9IH5tIKi++uYtv/mkpyzZ5ZpvTD5zKle86hPqaOF+562UeenULn799Cfe/spnvXXAoMRHueH49ty16k7U70gmjS9bt5pePruSfT53HBxbsndUXVqmEEzpzMXVsHb/75PH84L7XuP7JNaza1kZ9dZx/O32/nJ+57JS5/PH5ddz5/AY+c+b+ORc3yzb1zfk6YHoTi9/YxfLNLQMSNttauhhXX93n/6mY6gEBgWbz6dv+3qdm3NffeRA9iSRfuvMlnlq5g3+99QXue2UG3zp/YG20ixI2qnqjiFQB3wMagFuAjXgaxB8GdOX+HAesVNXVACJyG3A+fTWn84Er/dd34AlA8bffpqpdwBoRWemf75l8F3xzZzv//Nvn+21XBUX93977zp4Eezp6aO7s8X539JDLbTChoZonvnR61r7fYb5z76vc9PRavPuFSY01TG2qY+rYWuqr4ySSSlKVpEJzRw/Nnb0cPGNsH3txIWaMq+e2fzqBD1+3kNe3tHLO1U8we2IDTXVVNNZUMaauKjXO7kSSnl5/MkwqvYlk6jsIvpfgXWDGDX8FddUx6quraKiJU18Tp746Tk8iyeY9nWzc08HmPZ1s2tNJlx8NNKGhmqlNdUxpqmVqUy1NdVW0dPXS0tlLa2cvLV09bG3uYqtvNz92zgS+ed4hHBZyjO8/rYn3HhWMT1PnVoWkKgrERYoS/iftN4kl63bzg/teS03sAUrxdutk0tPEdrd7z8ru9m7aunOvez51xnw+m8U0E2bt9jaWbmxmbF1VXnMgeIuRhWt28sP7l3PjU2tTz2tYEM8cX88V7zqYsw5O5x1d+9FjuOuFDVz556U89OoWFv7nDjp7EilNcsa4Oj547N7sM6mBXz+2muVbWvjmn5byy0dX8c+nzuOYfSawraWLrS2dqb/bzvZuSvjqsrLPpAY+//YDCmrzDy7bwg1PruHqi47Mu6AKJtS6As9EdTzGN847mKNmj+fHDyzn0pPnMrUp93nnTRnDOYdM575XNnP9k2v46jsOynrcI6mcr6mISErYvLa5OeU7LJaVW1s4+7+eYPrYOj5/9v6cf4Rn9SimCGdAIJA6ehLMndLIt88/tE/U5m8/cTy/W/gm3/vrq/zlpU08u3pgJbGKthWo6rXAtSIyGYip6lD3Up0JhFNp1wPH5zpGVXtFZA8wyd/+bMZnsy7/ReQy4DKAmun7cV+emPpCNNTEqY7HiMfE+xFhe2sXu9p7WL+rnQOnj837+RVbPTttU10VbV29bG/tZntrN8vy1FP84LG5V7W5mNpUx22XncjHblzES+v3pPqvDxeNNXE6ehLsau9hV3sPywtU2505vp6vvuMg3nHY9LwagIgMaoV9xkHT+MWjq1i+paXgmEolHhMaauLUVsWojns/Xb0JtjR38XoR9vptvs18/rSmgvd47mEz+O69r7K9taufrX1iYw0XHbs3/+/0+f0EsIjwvmNmceK8SXzpzpd4YsV2YuLlJn3ouNmcesBU4v6Ef/4RM3lg6WaufngFr21u4Yp7lpbydZTMWQdP46jZE/Ie85tn1vLM6h08tWo77z1qVs7jOlNmtOJc1ucdvhfnHb5XUcf+86nzuO+VzfzmmbV87C1z2Gt8f5d24K8JAlYO9BePA/HbrN7WRiKpbNjdwb//YQnXPr6Gr7zjwKJybALOOng6j6/YzoXHzOLSk/ftZxoVET58wj68bf4UvnDHEhau2VnyOKFIYSMihwBxVX1JVbeHth8O9GaYuioaVb0GL4Sb+YccoVdffHTGfk/LELzfIIh4kSvj6qtTP2PrqqjKYvY6/xdPsWTd7lS/jHwEx9z08WM5YtZ4drR1s7W5iy3NnXQnksQEYuIJslhMGFNbxTEF/uFyMbGxhrv/5a2s2tbqaQ5dnvbQ1tVLS1cvAlRXxaiJS2oyrIpJanIP5njve5HU62BfUj1TX3t3go5uL8a/vTtBXIQZ4+uYMc77mT6uPhVqubOtO7US3tbSRavvKxlbV0VTXTVNvtY1e2JD1u96qDl69gRuv/zEnPlJQmFTF3jfx9i6asY3eM/K+IZqxtRW9ROUjy3fysdufI62IsLTg26sDUVoaOPqq/nrp0/mzZ3t/rPqjWNMXVVKWORjr/H1/OYfj+PZ1TuZM7mBGeP6T5ixmHDuYTM4+5Dp/N+yLVz/5GqaO3qZOrbW11TrmNpUy6QxNUVdMxe3PPMGC9fsZOXW1oLCZvU2zzHf1pX/f68YM9pAOXzWeN552AzufXkTP7r/Nf7roqP67N/W0sWS9bupqYrx1v28IJ8D/EXpQBY4wRyy39QxtHf1smxTMx+5flEqt6YYzeak+ZN59POnFjxu9qQGfv/JE1i0dicn/rDkoRat2VyD5095KWP7wXh5NieVful+bADCy/ZZ/rZsx6z3zXrj8AIFivlsP8bVV/OOwwbu8MpG8EduL/DAh4+pr/YE17SxdUwbW8dhjCvwyYERjwn7TyveBBcl8ZgwpcmbmA4mvwZYTo6dM7Fs1wqCEIpZmARO7WKEDXj+umyO7GIREU6clzviMSAWE845dDrnHDp9wNfKx8qtrSxcs5PVBSK8OnsSbNzjRdIVyi3LF402FHz53AN58NUt/O/fN3LJW+b0EZKPLd+KKrxl3iQaary//wH+/+TrW1pIJLUk4RwsVI6dM4Er3nUINz29ll88upIWP2cmXxHOgRCLCSfMLfxcZP1skccdDizKsv054LABXTn7ueaLyL4iUoPn8L8n45h7gEv81xcCj6gXBH4PcJGI1IrIvsD8HOONnOABKiaZsr3HO6ax1h3nqjF0BCvrYoRN4PNprKm8KLkomTtlDEDBxoBrtrel/IgFNZvu7Hk2Q8XeExu49KR9AfjWX5b1yVMJEozPCAWgjGuoZsa4Ojp7kryxo7RisOnIuirqquNcfso8Hv/CaXzy5H05Ye5Ejp9bvsVTIYoVNgnIutyeAEXaFgqgqr14WtIDwKvAH1V1qYh8S0Te7R92PTDJDwD4LPBl/7NLgT/iBRPcD/xrlJFo+QhWnkWZ0bqiXWEZlU1asyliYeIf0zDKFiZzJzcCaRNZLlaFhFGgueSiM0IzWsC/nLYfk8fU8uKbu7lniRdF2t2b5PHXvWTbzGjHAwbotwkEa3jBOqGxhq+982Buu+xEJldQqkOxwuZvwNf88GQAfDPW14DHh2owqvpXVd1fVeep6nf9bd9U1Xv8152q+n5V3U9Vjwsi1/x93/U/d4Cq3jdUYyqV4I9elB3eP2a0rVYNj8DkWmglDunFS8Moe1bmTvGEzdodbXnr1oWFUVuBqs1R+mwCxtRW8YWzvQjDH973Gh3dCRat2Ulbd4IDpzcxa0JfE2cQTFRq2ZrUIsSB56LYEX4ReBJYKSJP+ttOAsYAb4tiYK4S/NE7Cmg2iaSmy5w7lKNgDB2BRttRjGZTQoDASKKhpoq9xtWxcU8n63d1MMfXdDIJm9kKWRU6uvu2F4iKC4/Zm5uffoNlm5q59onVqaz+07PkcA00Iq29RF/ecFKUZqOqy/H8NrcCE/2f3wFHqOqr0Q3PPRqKXK2GV1eF8geMkUnKv9eTKNhddbT6bCDkt9me228TDiAoVrOJOhE1HhO+cZ5Xhf1Xj63iry97OQ3Z8qRSZrQSI9Laut1ZhBQdT6qqm1T1a6r6Tv/n66pauADTKCM9geR/4AP114IDRi/xmFBXHfOShnvzL06CFexo9O8FprRcfhtV7bOvEnw2ASfOm8Q5h0yno8fLqZrQUM2Re/cP4Z43ZQxVMWHtjraCVpEw7SmfTeUvQkrpZ9MgIm8RkfeIyAXhnygH6BoNRYY+W3CAAeHoxULCZvQuToIggVU5ItKCHK2AgppNmQX3V95xoFezDzgtlBgbpqYqxtwpjah6IdDF0h5xGPdQUmw/mzPxinA+CdyFVyom+Lk9stE5SMqMVsAOb8EBBhS/OAnMsi44goeaeVM9M9qqHJpNsL3YSND2Mmo2APtMakzVVAs3vMskldxZgt8m8OW5MI8Uq9lcDdwLzFLVWMZP5YvUMhKos4VU4VKT9IyRSTBJFFqcdPS4Y5sfatK5NtmFTeDLOWQvb7Iu9F12Rpxnk41PnTGfpf9xNm/J0ik2IAgSKCUirc2heaRYYTMH+Lb5aApTn9JsCqxUR2koq9GX+iJX46NZs5kxto666hjbW7tSHSjDBELo0JleKmDBhd4wmZ4K+VVSEWklNFJrH4EBAk8BBxQ8ykitVNsL2o3deUiM6Ah8MIUSO0ezzyYWE/adnLuSQLDt0L08YTOctdEGw0ASO4NFigsBAsWO8NfAVSKyF/AyXofOFKr6wlAPzFWKtRunV6qV9cAb5SXQVApNkKl8iurKn1SiYN6URl7d1MzqbW39CnIGYc8pzaYnkbfGWKWasGeOr6eptqqkRmou5V8V++Te4f++Jss+xeusaRAWNsWtVBscWJEY0RE8Lx0FQ+X9CXIUajaQO9emqzfBup3txATmTG6goSbuVR7vSeTsJ9VZpjybUhER9p/exPNFNlJT1VSwgwvm1WLNaPvm+ZkbzdDcpNhKvin114EViREdxWo2bQ5FHUXBvBy5Nm/uaCepXofS2qp4UYu94fLZFMMBqSCBwn6bzh6vuWFtVWxQbRzKRbGdOt+IeiAjhaIdvqlY/9E5eRgejUVMjomk14HU66sUfV+fSmTu5CD8ua9mE4Q9B4mfnvDu9kLJs3TTUNV0BYEi2pKXm1LK1qTSJxyxjpSS1HmuiPxFRJaJyN7+tktF5IzohuceDdXpPJtwafFMOlJ5NpW3ujLKRzE+vpTJtTqet1PpSCZdkLOdRKi0T2BWC4RRoTy3rl5PG6iJx8rSkK9Ugt42xZStaXfM71tsUufFeCX8V+CZzoL2b3G8Ip2GT1U8Rk2VV4Kkqzd3lVqX4uON6Ggowuya9te4sYKNgsbaKqaPraO7N8mGXR2p7aszNJtCZuy0v6byBA2kqz8HjdTy0e5Y7lWx3/gXgU+q6r8D4SXDs8CRQz0o10mXjs9jN7Y8G4PinhWXKvtGSSBQVoWCBIKw57QZLb+mWMn+GiitkZpruVfFCpv5wDNZtrdCBfX0rRCKqXeVcviO0ugiw6O+iJYUbanwVjcmlahICZutIWHjhz3Pm9LXjJYrzy3d2bJy/++KzbdxLfeqWGGzEdg/y/a3AauGbjgjg+Ls8BYgYIQ0mzwBAha56DEvFf7sCZidbd3sbu+hsSbOVD9MOF3+J79mU2lhz2EOKLJsTWoOcST3qlhhcw3w3yLyVv/93iJyCfAj4FeRjMxhAtt6/gnEAgSMYn023rNSqaafcpGukdba5/e8qWNSgRMNBSoydFa4GQ1IdfHc1tqV9zjXNJtiQ59/JCLjgAeBOuBRoAu4SlV/EeH4nCSISMtnGhmtbX6NvpSiBY/WHJuAoNVAEBSQCg4Ide8sZMJOdemsYM1mbJ13Dy2dBSrHO+azKThKEakC3g78BPgucDCeRrRMVXO3zhvFBCsNc/oahWgoJUDAkRVsVMwcX09tVYytLV20dPakcm4CjQeK8NlUaF20MEHlg9bO/kVHw1Rq2Z1cFDSjqWovXg+bJlVtV9XFqrrIBE1ugpVGvo6B6XI1bjwoRjQ0FhFM4lJl3yjxCnKmtZvMhE4owWdTwd9lU52XWVJQs3HMFF+sz2YJsF+UAxlJpFer+SKM3FKBjWgoxowWPCuj3YwG4SCB1n4JnVCEz8aBaLSUZlOgcrxr+VfFCpsrgf/0W0LvLSITwz8Rjs9J0nbj7A9LMqlOqPNG9KQDBPLlZFnoc0Cgxby+pZU3d7QDpLQdKLzQc+H/rqlon41bGm+xT++9/u+78Ko8BwhW9bkfRSeWVcedKKBnREcw6bV3J0gmlViW58GqTaQJhM1jy7fRm1Rmjq/vE1lWMEDAgWi0tLAp1mfjxiKk2FGeFukoRhiBKp8r9NmCA4yAeEyoq47R2ZOkszeRdeIw/16awIz26iavKnLYXwNhH1j+pM5KzrNpDJnRVDVnPTzXfDbFhj7/LeqBjCTS3Tqzr65s8jDCNNZU0dnTTVtXLmFjPpuAsMkM+oY9Q3ihl782WiWb0arjMeqr43T0JGjvTuSs6pxODK/cewlTStXnw0Tk5yJyn4jM8Le9R0SOim54blKozUAqOMCRzF8jWgo5tYPnxZVJJUqa6qpT1QKgb9gzlBL6XJmFOAPG1BUOEnCpJTQUX/X57cBzwEzgdKDe3zUPuCKaoblLQVW+xzQbI02w6Mi1OElXm3BjUomasOkstxktV1KnG4K7mCAB1wIEihXv3wY+q6rvBbpD2x8DjhvqQblOeqWaX7OxycOAwpqNa+aSqAlrMzk1m5wLvcr32QA01RYOEnDNvFqssDkU+GuW7TsBC33OoKHaJg+jeFKJiAV8fK7UwIqaIEigrjrGjLF1ffY1puoS5vfZVHoEVylmtJGm2ezEM6FlcjSwfuiGMzJIPfCFJg9HHhIjWor18bmygo2aeb7pbO7kMf1CxWurYsQEunuT9Cb6Ny9sdyCpE6CptnAVgXSgkRvPRbGjvBX4sYh8AC+vpkpETgGuAm6ManCuEkweucrVtKU0GzceEiNaGos0/Zgm7HHSfpP5+FvncMr+U/rtExEaaqpo7eqlvSfB2IzWz+nv0pEAgRzCJplUZwRnQLGz3deBm4A38BI5l/m/b8UrzmmESJtFcsX6m2ZjpGkoYPpJNdqzxQngtV6/4l2H5NzfUBP3hE1XgrF11X32uZBnA6EAgRxzSGdvur21K4nhOZ9eEZkNrFOPHuBiEfkGnuksBryoqivKNE6nCBy+uVoMpEKfHVF/jWhJt6ToP7EkkkpXbxIRb2IxCtNYWwUtXVmTql3Is4HCAQIumlbzPb1rgCkAIvKIiIxX1dWqeoeq/nEoBY1fY+1BEVnh/56Q47hL/GNW+M3bgu2PichyEfm7/zN1qMY2EILJo63bywDOpKPHLceeES0NeXx8Kbt8dTxnJrnRl3SuTf/v0xWTZCEzmouJ4fmETQsw2X99KlCd+9BB82XgYVWdDzzsv++DX/DzCuB4vHDrKzKE0sWqeqT/szXCsRakKh6jpipGUqGrt7+TMm0WcedBMaIjn8/Gtcq+lUC+8OcOR/wchdoMuJgYnm+kDwGPiMir/vu7RaQ724Gqevogx3E+nkADuBkvf+dLGcecDTyoqjsBRORB4Bzg94O8diQ01MTp7k3S3p3oZx/usAABI0S+wq22MCmdfMU4O3u8xV+l+2wKtRlwMTE832z3EeAf8frYnAIsB9ojGsc0Vd3kv94MTMtyzExgXej9evqGY98oIgngTuA7ms1+BYjIZcBlALNnzx7suHPSWFPF7vYe2rp6mdhY02efawX0jGjJNzm228KkZBpzFMLtTSTpTnj+r9qqyvZ/BQECzSPIZ5NvpFOAX6qqisiRwOdUdfdALyQiDwHTs+z6WviNf72sgiIPF6vqBhFpwhM2HwF+k+1AVb0GuAZgwYIFpV6naBryhD+bacQIk6+NeDpL3BYmxdKQoxBup2/SrnfA/9VUIKkzMBFWuu8pTL7Zbg0wA9hK3x42A0JVz8y1T0S2iMgMVd3kF/nM5nPZQNrUBjALz9yGqm7wf7eIyK14Pp2swqZc5Ost71rmrxEtgdaSbWHS5ljiXiWQy2fjir8G0j6b3AEC7i1Cig0QOIVoAwTuAYLoskuAP2U55gHg7SIywQ8MeDvwgIhUichkABGpBs4DXolwrEWRzzTiWgE9I1oa8yxMUg2yHJggK4Xgfy8zb6nTkbpokPbZ5AwQcNA6UmyAgBBtgMAPgD+KyCfwEkc/ACAiC4DLVfVSVd0pIt/Gqz4N8C1/WyOe0KnG6xj6EHDtIMczaPI5fTscqc9klIeiFiYOOYKHm1zRfa6EPUPh2mhBCwWXFiEVESCgqjuAM7JsXwxcGnp/A3BDxjFtwDFRjGsw5Ostn3buufOgGNFRzMLEJUfwcJMrb8klM9qYmrSwSSS1X5UAF/2+OUeqqh3ALwCGIkBgtJE2jeRJ1HPoQTGiI1+LgVQ+hS1MiiYVnJMhvDscqR4AEIsJY2q9Gm9t3b39yu64WMy3qPg/VT3NBE1p1OdQ5ZNJdeqhN6InX4uB1MLENJuiSQXn5DCj1TkyQTflqSLQ5mCQUb7aaP8NfEVV2/zXOVHVTw35yBwnV8fAzt4Eqm4V0DOiJVh0dPQkSCa1T9n8dOtfdyaV4Sbn/17KjFbZOTYB+YIEUoEjDi1C8o30MNIRaIflOS6yXBWXydWjpN3Bh8SIllhMqK+O09GToKMn0aenvIv5FMNNQ468JdcsCukggf6JnanKEg4tQvL5bE7L9toojlwRMe1mgzey0FjrCZu27t4+wsbFTPHhJld0n0vRaJDOtWnOotm4uGh1Q590kFwRMelSNe48JEb01OdwalsCcOnkDH12pJdNQNBmIJvPJu3Lc+NeoAhhIyL1InKFiLwkIq0i0iIiS0Tk6yJSX45Buki6XE2GZtPt1urKKA+5ggQsQKB00mkH2ZM6XTGj5StZ46Jmk3ekIlIFPILXMO1+4F68BM+DgW8C54rIKaqau1H2KKXQ5OGSrdWInlwlVtKZ4va8FEuuigzO+WzyNFBrc3AeKSQWL8NL6jxaVZeGd4jIocCjwCeBX0UzPHcpbBZxZ0ViRE9OP4OZXUsmV3BOR3eyz/5KJ199tMD368q9QGEz2oXAdzMFDYCqvgJ8H3h/FANznZRmkxkg4KCt1YienJqNBZSUTE08RlVM6E0q3aHmha4FCIxJtRnIVw3cnUVIIWFzCJ4ZLRcPAYcO3XBGDums8IwAgS7TbIz+NOYIKLHFSemISFbh7ZzPJkcDtYSjieGFhM0EYFue/duA8UM2mhFEoTLnNnkYYXJVnEgnddripBRSwju02HOpNhqkAwQyfTZhQRNzKDG8kLCJA/mc/0n/GCODXA2crEunkY3GLH6G3kSSrl43OktWGinhHdIK2h0rV5Or8rOrQUaFlksC/FZEunLsrx3i8YwYwvWZVDXVGbDDwWqtRvRk68HSHqr4XOmdJSuNbCVrOp3TbLIHCLQ7aoovNNqbizjHsHbErFSq4zFq4jG6/dVpkEjWZjZ4IwsNWVbiZnIdONmKcbrm58hVG83VOSSvsFHVj5drICORhto43e1J2rsTKWHj6qrEiJZUImKoNbR1dB04gc8mbMZ2LRptbOCzycwXcnQRYobgCAm66IWdvlZ+xMhGYzYfg+VkDZismo1jZrQxOQIE2hwNGjFhEyHZyma4qgIb0ZKeHEM+G2svMGCyNVALQp9dqY1WXx0nHhM6e5L0JNL5QsGCxBWhGWDCJkKylc3ocHRVYkRLoL10ZFmY1JtmUzLZAi5cM6OJSMpv05ZF43VtDjFhEyHZSta0OabKG+Uh0F7CZp/2VHsBe1ZKJfg+Ay1ANZ0IWedQGHm2IAFXE33d+dYdpDHb6ioVI+/WqsSIlvrq/g5tq/g8cDI1m67eJKp+KZu4O9NeOrEzLWzMZ2P0I1tWuIu9w43oSa3EeyyYZChI+2y87zPtr3FrystWRcB8NkY/siWWtVs4q5GFbBUnUsLGAgRKJtOq4Jq/JmBMlvporgaOmLCJkMxe6KqayqMw04gRJluobru1Fxgw6UK43nfoWthzQKqKQFc264hbz4UJmwjJDL/s7PHsxrVVMeIOFdAzoieYBDt7kiSSClh7gcGQ2byww7Gw54BsbQYsQMDoR6aTst2CA4wcxGLpsvjBxBi0FHdtBVsJZEaCdjpqRku1hu50P9nXhE2EZDop2x1V5Y3ykFkfLViVu2abrwQymxemunQ69r/XlKU1tKtVn03YREimk9LFvuFG+chsDR1MKq5NkJVAZvNC14pwBmT12ThqXjVhEyGZTkpX1V+jPGQGCbiaKV4JpH02vmbjWC+bgGxJnR2OziMmbCIk3a3TX6k6uiIxykNjRi09y8kaOPUZ/3tBL5sGxzSbMVmTOt2MUjRhEyGZuROWEW7ko//ixJ6XgRJuyx4uVeNsgEBX2Gfj5r2YsImQXGYRW6ka2cgMELDnZeBUx2PUVMVIqleqxlmfTa3ns8lWG801368JmwjJrORrAQJGPjIDSixUfnCEq64H/4Ou5dmkNRvvWUgklc6eJCJQV+XWvVSEsBGRiSLyoIis8H9PyHHc/SKyW0T+krF9XxFZKCIrReQPIlJTnpHnJ7OSr6uOPaM81GeEypvPZnCEo/tczbPJ9NmEIxRjjiWGV4SwAb4MPKyq84GH/ffZ+DHwkSzbfwj8VFX3A3YBn4hklCXSUJ3h8LUAASMPgQbT1p2gN5GkuzdJTLyKE0bphH1grprRUrXROt1fsFbKU3w+cLP/+mbgPdkOUtWHgZbwNhER4HTgjkKfLzfhiBivLpo5fI3chH02QQ29xpoqvEfcKJWG2nRip6u10eqq49TEY3QnknT1JkLtBdy6D6gcYTNNVTf5rzcD00r47CRgt6oGHrT1wMxcB4vIZSKyWEQWb9u2bWCjLZKaqhjVcSGRVLp6k+lmWA4+KEb0hFfiwbPimtmnkmhMCe+Es3k20NeU1uZoewGAsi2xReQhYHqWXV8Lv1FVFRGNahyqeg1wDcCCBQsiu05AQ00Vezp66OhOpNv8OvigGNETrqVnwQGDJxwN6qpmA16QwM62blo7e51O9C3biFX1zFz7RGSLiMxQ1U0iMgPYWsKpdwDjRaTK125mARsGOdwho7Emzp6Onj4PvIsPihE9jaGKExb2PHjC0aCu+mygbxUBVys+Q+WY0e4BLvFfXwL8qdgPqqoCjwIXDuTzUROuPtvmaDKWUR7qQwElbdZkb9CEo0HTSZ2VMuUVT6pbZ1dPWrNx0O9bKd/8D4CzRGQFcKb/HhFZICLXBQeJyBPA7cAZIrJeRM72d30J+KyIrMTz4Vxf1tHnIRxh1OFomQmjPPTRbKzJ3qAJV/BwNc8GYIyf2Nka8tm4uAipiCdZVXcAZ2TZvhi4NPT+5ByfXw0cF9kAB0E4wshCn418NIQaflkwyeAJ+2w6HTajNYUCBAINzcVW4RUhbEYy4cSy1INiwsbIQkMfk2sQTGL/ogMlq8/Gwf+9cBUBl81o7o3YMcKrq0AFtgABIxvhhl8dDudTVAp9fDYOR6ONCTVQ6+r1m8A5KDRt1ouYPrkTFiBg5CHc8KvNKoQPmrDPprPHm6Rd9NkEDdRaunrpTXjZGqbZGP1oCDVxSoUtOvjAG9ETLovfYaHPgyZI6mzu7KU74RWvdLH0Tzip0wu+NZ+NkYVAld/d3kNSvYe9Ku7eA29ET11VHBHo7EmmCi+asBk4gQVhZ1uX97467mTpn6aM+mjg5nNhs17EBJrN9tYu/717D4lRHmIxSWm9wfNi/r2BE3x3O9q6ATf9NZA9QMBF86oJm4gJhMu2lkDYuPeQGOWj3hYnQ0bw3e1o9YSNi/4a6Bsg0O5wrp4Jm4gJHnibPIxiCMyutjgZPMGEHDQeczUwJxUg0NnrdBUSe5IjJm1G81ZXDWYWMfKQ+bw0OjipVAqZCzvXzWgtnb0k/QABF0PibeaLmMyVqk0eRj6CCXJPRw/g5gq2Ushc2Ln6XWa2hgY3zWjujdgxggzw7oQX529mNCMfmc+HBQgMnExNxlXNJngGWrt6CYLpXBSc5rOJmEx112zwRj4yV6y2OBk48Zj0ETCuCpvqeIy66hiJpKY0Xhc1GxM2EZM5WdjkYeSj//Pi3qRSSYS/Txe1gYAgSEAVRKCu2r2p270RO0bmZGGTh5GPzMxwW5wMjvD36WroM6QTO8HTalxMTjVhEzGZ6q6LUSRG+Qg/L/GYOFlepZIIf5+umtEgHSQA7mpo9iRHTOaD4eqDYpSH8PPR4Gh5lUqirxnN3eluTF1Ys3FzDnH323eEmqoY1fH0hOGiY88oH+Hnw8Vii5VGwwjRbMaEzGiumuJN2JSBPhExjq5KjPIQFjC2MBk8Yc3GaZ+NHyAA7vrxTNiUgcYM555h5GKkRE9VCuH/PZe/zz6ajaO5VyZsykB4AjHTiJGPsInEFiaDp4/wdlizGWs+G6MYwhOINU4z8mE+m6FlpAibMRaNZhRD+IG38iNGPvpEozk6qVQS4YVencPfZ9hn46rGa8KmDJgd3iiWcB6Wq1FHlUT4+3Ras6l1X+M1YVMGGixAwCiSxhr3bfOVxIgJfa4Lm+LdnENM2JSBRgsQMIqkvo8W7OakUkmMFKtCnwABR+cQEzZlwAIEjGIxzWZoGTGaTW04z8bNRYgJmzIQrK5qqmJUxe0rN3JTVx1L9SxxNZ+ikmgcKYU4w2Y0RxchNvOVgSACzVaqRiFEJKX9ujqpVBJ9NBuHv88xJmyMYqhPTR62UjUKE2g0rk4qlUSfcjUOV9AeEzavOqrxuvvtO0SgytvkYRRDoAFb5OLgCb7DmrjbJuxYTFLhz65qaO5++w4RRBWZsDGKwZ6XoSPwdYR9Hq4yJmWOd/Ne3By1YwQrVVdXJEZ5mTelkeWbm5k9qWG4h+I8Expr+Po7D2Lq2LrhHsqgecu8STy7egd7T6wf7qEMCFHV4R7DsLFgwQJdvHhx5NfZ1dbNR29YxD8cP5sPHTc78usZbtPVm2BHazd7jXdzUjGiQVVJqtfBdbgRkedVdUEpnzHNpgxMaKzhz//vpOEehuEItVVxEzRGP0SE+PDLmQFTET4bEZkoIg+KyAr/94Qcx90vIrtF5C8Z228SkTUi8nf/58iyDNwwDMMoiooQNsCXgYdVdT7wsP8+Gz8GPpJj3xdU9Uj/5+8RjNEwDMMYIJUibM4HbvZf3wy8J9tBqvow0FKmMRmGYRhDRKUIm2mqusl/vRmYNoBzfFdEXhKRn4pIba6DROQyEVksIou3bds2oMEahmEYpVE2YSMiD4nIK1l+zg8fp154XKkhcl8BDgSOBSYCX8p1oKpeo6oLVHXBlClTSr0NwzAMYwCULRpNVc/MtU9EtojIDFXdJCIzgK0lnjvQirpE5Ebg84MYqmEYhjHEVIoZ7R7gEv/1JcCfSvmwL6AQEcHz97wylIMzDMMwBkelCJsfAGeJyArgTP89IrJARK4LDhKRJ4DbgTNEZL2InO3v+p2IvAy8DEwGvlPW0RuGYRh5GdUVBESkBVg+3OOIkMnA9uEeRESM5HsDuz/XGen3d4CqNpXygdFeQWB5qSUXXEJEFo/U+xvJ9wZ2f64zGu6v1M9UihnNMAzDGMGYsDEMwzAiZ7QLm2uGewARM5LvbyTfG9j9uY7dXwajOkDAMAzDKA+jXbMxDMMwyoAJG8MwDCNyRqWwEZFzRGS5iKwUkVztDJxBRG4Qka0i8kpoW1E9glxARPYWkUdFZJmILBWRT/vbR8Q9ikidiCwSkSX+/f2Hv31fEVnoP6d/EJGa4R7rQBGRuIi8GPSiGmH3tlZEXvZ7aS32t42IZxNARMaLyB0i8pqIvCoiJw7k/kadsBGROPAL4FzgYOBDInLw8I5q0NwEnJOxrdgeQS7QC3xOVQ8GTgD+1f+bjZR77AJOV9UjgCOBc0TkBOCHwE9VdT9gF/CJ4RvioPk08Gro/Ui6N4DT/F5aQW7NSHk2Aa4G7lfVA4Ej8P6Opd+fqo6qH+BE4IHQ+68AXxnucQ3Bfc0BXgm9Xw7M8F/PwEtgHfZxDtG9/gk4ayTeI9AAvAAcj5eBXuVv7/PcuvQDzPInpNOBvwAyUu7NH/9aYHLGthHxbALjgDX4wWSDub9Rp9kAM4F1offr/W0jjaHoEVRxiMgc4ChgISPoHn0z09/xKp4/CKwCdqtqr3+Iy8/pfwFfBJL++0mMnHsDryXK/4nI8yJymb9tpDyb+wLbgBt9M+h1ItLIAO5vNAqbUYd6yw/nY9xFZAxwJ/AZVW0O73P9HlU1oapH4mkBx+H1Z3IeETkP2Kqqzw/3WCLkJFU9Gs80/68i8rbwTsefzSrgaOBXqnoU0EaGyazY+xuNwmYDsHfo/Sx/20hjS6j1Qsk9gioNEanGEzS/U9W7/M0j6h4BVHU38CieaWm8iAT1C119Tt8KvFtE1gK34ZnSrmZk3BsAqrrB/70VuBtvsTBSns31wHpVXei/vwNP+JR8f6NR2DwHzPejYWqAi/D66Yw0BtUjqJLw+xRdD7yqqj8J7RoR9ygiU0RkvP+6Hs8f9Sqe0LnQP8zJ+1PVr6jqLFWdg/e/9oiqXswIuDcAEWkUkabgNfB2vH5aI+LZVNXNwDoROcDfdAawjAHc36isICAi78CzI8eBG1T1u8M7osEhIr8HTsUra74FuAL4X+CPwGzgDeADqrpzmIY4KETkJOAJvH5Fgd3/q3h+G+fvUUQOB27Gex5jwB9V9VsiMhdPG5gIvAh8WFW7hm+kg0NETgU+r6rnjZR78+/jbv9tFXCrqn5XRCYxAp5NABE5ErgOqAFWAx/Hf04p4f5GpbAxDMMwystoNKMZhmEYZcaEjWEYhhE5JmwMwzCMyDFhYxiGYUSOCRvDMAwjckzYGEYWROSmoEJxpSAi5/tVdntF5KbhHo9hlIIJG6Pi8Cd6FZFvZGw/1d8+ebjGNsxcj1dFYR+8KsqRICIfE5HWCM+vInJh4SONkYQJG6NS6QS+ICJThnsgQ4lfdmcgnxuPV8DyAVXdoKp7hnRghhExJmyMSuVRvNLt38h1QDZNR0Tm+NsWZBxzrl+Vt0NEnhCRWSJyit+wrFVE/uJnfWde4+sissU/5ka/nEywT0TkiyKyyj/vyyLy4Sxj+ZCIPCIiHcA/5biXCSJys4js8s/1kIgcEtwDXs8XgEf8c55a6nn8/f20lvD36J/3RqDR36YicqV/3FoRuVJEfut/H5tF5PMZ5+qntfif+3zw2t98u3/sWn/73iLyJxHZKSLt4jXquijbPRpuYsLGqFSSeNVlLxeReUNwvv8APoPXJ2YC8Afgm8BleKV+DgGuzPjMKXjNos4A3odX9+qHof3fwWv69a94jfi+D/yPiLwz4zzfB37pH/O/OcZ3kz+28/EKObYD9/vC7Wl/fPjjmOFvK/U8xfA03vfU7l9nBnBVaP9n8eq2HY1XFul7InJBkecGONb//Un/3MH7X+L18jkN714/A+wu4bxGhVNV+BDDGB5U9a8i8hTwXbwijoPhG6r6BICI/Br4GXCMqr7gb7uZdGHIgATwcVVtBV4RkS8B14vIV/z9nwXeHpwXWCMix+EJn3tD5/mZqt6Ra2AiMh94N3CKqj7ub/sI8CZwsapeJyJBVd2dfnHEks+DV98qL6raLSJ7vJdZr7MwVEvwdRE5Fu97uCvLsdnOv01EwOtnEz7/PsCdqrrEf7+mmPMZ7mDCxqh0vgQ8IyI/HuR5Xgq93uL/fjlj29TMz/iCJuAZvGKE84BaoA5PawgXGKzGM/+FWVxgbAfhaXLPBBtUdY+IvIynDRXLUJ0nH89keV+KZpOLq4Ffi8g5eF097x7hPXBGHWZGMyoaVV2EF4H1oyy7gwrQEtqWywHfEz6tf+7MbaX8PwTHvgs4MvRzCJ65LUxbCefNZKgq5QbnSdL3+4Lc39lAr1Py+VX1eryukDcC+wNPB74iY2RgwsZwga8CJwPnZGzf5v+eEdp25BBe9zDxepQEnAB047VsXgZ0Afuo6sqMnzdKvM6reP+LJwYbRGQscJh/naE8zzagwd8ecGTGebrx2h1k44Qs718Nvd9G6O8hItPo+/cBT/D3O7+qrlfVa1T1A6T9acYIwYSNUfGo6krgGvrnlqwE1gFXisj+IvJ24OtDeOkq4AYROUREzgJ+AFyrqm2q2oLnOL9KRP5RRPYTkSNF5HJJ96EvClVdgdd86n9E5GQROQz4LdAM3DrE51mIp2l93x/z+4B/yTjVWqBORM7yI9QaQvtOEJGviMh8Efkk8FHgp6H9j+C1Rl4gIkfhBSx0Zjn/GSIyXUQmAIjI1SJyjojMFa9/yjmUJmiNCseEjeEK3wJ6wxt8M9hFwFxgCV7E2VeH8Jp/A5bihWHfjTeRfjG0/xt4EWyf9497EC9abCDO7Y8Di/A6IC7Ci8w6R1U7hvI8foOri/G6gb6Mpz30CS9X1aeBXwO/x9NUwvf8E+BwvIZn3wG+mRH88Dm8BluP4bUQvo7+LYM/hxd1ts4/D3hz0c/wBMyDeD60SzBGDNY8zTCMovBzYn6uqlcVOtYwMjHNxjAMw4gcEzaGYRhG5JgZzTAMw4gc02wMwzCMyDFhYxiGYUSOCRvDMAwjckzYGIZhGJFjwsYwDMOInP8PSm2LPYDgCxQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(res['Matlab']-res['C++'],linewidth=2)\n",
    "plt.xlabel('Number of outputs',fontsize=14)\n",
    "plt.ylabel('Difference ',fontsize=14);\n",
    "plt.title('res[Matlab]-res[C++]',fontsize=14)\n",
    "plt.xlim([0,60]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MAE: 0.0136\n",
      "RMSE: 0.0014\n",
      "R2 Score: 1.0\n"
     ]
    }
   ],
   "source": [
    "rmse=sklm.mean_squared_error(res['Matlab'],res['C++'])\n",
    "mae=sklm.mean_absolute_error(res['Matlab'],res['C++'])\n",
    "r2=sklm.r2_score(res['Matlab'],res['C++'])\n",
    "print('MAE:',round(mae,4))\n",
    "print('RMSE:',round(rmse,4))\n",
    "print('R2 Score:',round(r2,4))"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}