Turizam.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Turizam\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using numpy backend.\n"
     ]
    }
   ],
   "source": [
    "# Standardna Python biblioteka.\n",
    "import copy\n",
    "import functools\n",
    "import math\n",
    "import numbers\n",
    "import warnings\n",
    "\n",
    "# SciPy paketi.\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# Paket TensorLy.\n",
    "from tensorly.decomposition import parafac\n",
    "\n",
    "# Paket StatsModels.\n",
    "from statsmodels.tsa.seasonal import seasonal_decompose\n",
    "from statsmodels.tsa.stattools import adfuller\n",
    "\n",
    "# Paket Hyperopt.\n",
    "import hyperopt as ho\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tlp import *\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Učitavanje podataka\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_df = pd.read_csv(\n",
    "    'data/BS_TU12.csv',\n",
    "    header = 0,\n",
    "    index_col = None,\n",
    "    encoding = 'utf-8'\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "df16 = raw_df[['Županije', 'Mjesec', 'Tip', '2016 Domaći', '2016 Strani']].copy()\n",
    "df17 = raw_df[['Županije', 'Mjesec', 'Tip', '2017 Domaći', '2017 Strani']].copy()\n",
    "df18 = raw_df[['Županije', 'Mjesec', 'Tip', '2018 Domaći', '2018 Strani']].copy()\n",
    "\n",
    "df16.columns = ['Županije', 'Mjesec', 'Tip', 'Domaći', 'Strani']\n",
    "df17.columns = ['Županije', 'Mjesec', 'Tip', 'Domaći', 'Strani']\n",
    "df18.columns = ['Županije', 'Mjesec', 'Tip', 'Domaći', 'Strani']\n",
    "\n",
    "df16.reset_index(drop = True, inplace = True)\n",
    "df17.reset_index(drop = True, inplace = True)\n",
    "df18.reset_index(drop = True, inplace = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.concat((df16, df17, df18), axis = 0, ignore_index = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "županije = np.array(df['Županije'].unique())\n",
    "po_županijama = df.groupby('Županije')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z = np.zeros((4, len(po_županijama.groups), 3 * 12), dtype = int, order = 'F')\n",
    "for i in iter(range(int(županije.size))):\n",
    "    aux_df = po_županijama.get_group(županije[i]).reset_index(drop = True).copy()\n",
    "    Z[0, i, :] = aux_df.loc[0::2, 'Domaći']\n",
    "    Z[1, i, :] = aux_df.loc[0::2, 'Strani']\n",
    "    Z[2, i, :] = aux_df.loc[1::2, 'Domaći']\n",
    "    Z[3, i, :] = aux_df.loc[1::2, 'Strani']\n",
    "    del aux_df\n",
    "try:\n",
    "    del i\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('data/turizam.npy', Z)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "del raw_df\n",
    "del df16\n",
    "del df17\n",
    "del df18\n",
    "del po_županijama\n",
    "del županije\n",
    "del df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predviđanje\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z = np.load('data/turizam.npy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z_mean = Z.mean(axis = -1).copy(order = 'F')\n",
    "Z_std = Z.std(axis = -1, ddof = 1).copy(order = 'F')\n",
    "\n",
    "Z_norm = (Z - Z_mean.reshape((4, 21, 1))) / Z_std.reshape((4, 21, 1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z_log = np.log1p(Z).copy(order = 'F')\n",
    "\n",
    "Z_log_mean = Z_log.mean(axis = -1).copy(order = 'F')\n",
    "Z_log_std = Z_log.std(axis = -1, ddof = 1).copy(order = 'F')\n",
    "\n",
    "Z_log_norm = (\n",
    "    (Z_log - Z_log_mean.reshape((4, 21, 1))) / Z_log_std.reshape((4, 21, 1))\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO19f9BeR3XecyxbEgjJso1qVH2WZUhjh3pa4hgKCZNxAaeUMjWd0RhIYUwgo9ACpZhMLBgCTqa0TqeEMNMW6sYUM5Nguy7EDEybaKhtkn8o2PwMxhQTSZEq2RBsbDJYisXpH+977aurvbtnd8/Z3fvqPjMafd/7ve+9++7Zfe7Z85OYGTNmzJgxY7VwRu0BzJgxY8YMfczkPmPGjBkriJncZ8yYMWMFMZP7jBkzZqwgZnKfMWPGjBXEmSVvtv6sTbxx4zmi957YSFHXPrE+fjy8XidSaP36J6I/c/y47tTT8bj5CmHdcdXLAUiTEZAvpxT5WCNG/tqylcK1BkrJUFtmW8583Pn6o09sVLuHVKYhebrm6vj+w99n5m0x4ylK7hs3noPnX/4W599++JwNKvd47EK9jXBsp5zhdq19L/r6+w9FyQoAsOFg4u6KwOYD6WSqOf8uxMhkDDGyuvJZ3xK9b9/RS7x/T5E1UEbeligtrz6ksnMhJM8OsXKVyNM1ZwfesPdA1I0AUMk49y2b17gjdy0y98GCaHyLVZPgS2/qFEJ3za9kM5f4bpJxjMlLmxRSiR2wmSsNwk0lqVxI9liO/PoYI/gSD+rh3B14w957mPnymPsV1dyB1SP1VK3C9fn+ojm287gpCeZo5zPKwEr+Gw6uzyLemqeJbo/k7rvSiJ2zXBkBhTX3Tdsu4J+56h1q18sh8ZyJK7GwLI571oRuZZKpoQHGan8WZpmpm2QAfdnF7r0ULV5ikrHYnz58+73X2mjuRPQOAL8KgAF8HcCvANgO4BYA5wG4B8DrmdnADXcyUgikpt1P436ShdR9x7FFVEJTH7tHrYdwH7Hyk9pcU+6fY6bpUOvkFStLDQ20Q8oe1JZjh24sUllan8RdCJI7Ee0A8K8APJeZf0xEtwF4DYBXAPggM99CRB8B8CYAH7YaqLWjzoepHAFrErsPvvvXlOtUUVOe/XtLZadB8FPZg2OosTelce5nAngaEZ0J4OkAjgB4CYDbl3+/GcCr9Ie3wEwAYbRK7C2gJWKQaHobDq5vXp7znsyXpTWC5M7MhwH8BwAHsSD1H2JhhnmEmbtg1EMAdrg+T0R7iOhLRPSlJ3781zqjLoiWiGEMrRPBDDl8RNCKPFOIfRX8B7EIfWdreQbJnYjOAXAVgIsA/G0AmwC8XHoDZr6RmS9n5svPfNqmpEFuPsBZE5GzsDTso9YYO/JOQbsKyfV0IwWf+aIVebbykGkdx3YeD8rTUqYSs8zLAPwFM3+Pmf8GwCcB/AKArUszDQCsAThsNMYnkUvyq4wpE7w1WnpAS06Cq0rwq6ZkSU/1IX9DR/LaZC+JljkI4IVE9HQAPwbwUgBfAnAngN1YRMxcA+AOtVEpI9eZUzu2dtfa94KL27dxfAum9sPSNzbNMLraMtREf85qO1elZKQhy9ZkaPXA6eY0V7aiOHci+i0ArwbwBIAvYxEWuQMLYj93+drrmPmY7zpace4lF9QYaiyw2tms2kQyJkdLufWhIcMrn/WtpHA7HzGkyrPF8MhSsgTK7kkpsfdl6ZuLkK/lnpveGR3nPrkkptLEHhKO5YIaEsdwQdW0R+cSSSyxa8eJj11bgrHEmBiSlzyoU8ggRS5nP+DVyZ7EWHZ5jCyHc22l/VrtS2vz0JhczZKYpgit5InShN5HiCxKJkaU0grH5GaxqVLl2JdLquYuQaxsU2WUUxLEp2y5ZDk1Mh+7T4rmDow/sC0i3iZF7lPLjtPE2GIKZaZKUNvuPkT3XWqnrkugbZKJfWDXlJ1kP1rIsvZe7MYgzRzvy7NkeOSkyD0Fq5AdN+V09RxoEkNtGcZASvAt2tjHkCLL1mUm1eJjHtiuuU2Vc7WSvx1SjoQxi6tW2VFtSLPhfGiR4Gs4x1shDQ2Zdigh29wwvZwyzFOBhaMcmKjN3eXMCRF+t5Brxfz2BTilxagVYqWJmHA6LbQiP8nRXmp2swyP1JBPiNintI988GnzpYuHVSf3DiXqvGtiqouxJWKviRbkl+qU88FCviUewPsPbWtCJhpoJeGqKLmf2EjF2+mtgt12VU0ypdESeUzNJJNSDbIPie+rlRNVCnKSDK3QjOYuQY72kONYndpCG6JVYpfIs/Vop9hQyKmRuu++ltFrw3nSlp+viUdsBJSlpp4j46IO1Y07LuAL33wtANmgS9hiY8mjFNGvurZektiHsGi4rNWJKUXDqyHnEs7VMZToxgSkN/qwcKqmOFSrkXurkC66WhlwUyZ0oC6pu+CTo1arPU1tfQyl5K6tcFmTvHWz7BC0iH6S0TKtwTXhoYzJXKKPOdaFPO6Szdf6A6AkfMWo9h29JJscShA7UC4SStO56iP21kyhqZnIvoio4ffXtss3obmnLshaoZBaiRgatrpSGp9mZcmQ3Epq7kC+mSZHYwf0N3XJh3fsHqwREqmlvQN5PVlz9vuBN+xtX3MPLbxh3LsvusZ1rRKE79uMpWppjN0PiCMLn8YnnctW4uelRbe04NvotYgdOFVuOXIJXWv4e2qlz9Y09TFY1hLSRlFyXyfYb0MyD1WsG74/N2SrBGKLD4Wgdaz3bdyxuUwhjtDRPiaqYuy7p0RmaBKMtPZIDFLnRHMfaGvqHSxrtWuY11JQO969qZK/Eg0wluz713Whv/hSSTLHCZu6AGLGGlNcrLT2bf0Arh3+mmtz1z6BhNZA7P3GrtdaXkmI3HdvuRcAcPujl3nfpxUVBcQpeZMwy0jg0iI7aCVBuRafpByndjRNbJ3yHC29dPqzBJalJHIIRkuTlJYYKCUX371S5sty7GPzVsuEk+tbKY2mNHctlHLY1WwIoLmhWi8bK0XLham0ZJprrkpBTlnpFss2a4W4AvrEXrRZBxFdDODW3kvPBvBeAB9fvr4LwH4AVzPzwzE3T0UqIdQ+osfeL1T3G5hO0osLpSpCtuCs09LgS0RHDeWSc09NJ7eWHIdk7SJ7TfOLCyVOakFyZ+b7ATwPAIhoHYDDAD4FYC+AzzHzDUS0d/n7dYZjPe2g6ZRrhdBrYJWKUvVhJdOY68YoWprNczRhGf1Ss+jbGZHvfymAB5j5AICrANy8fP1mAK/SHNgYWo2AsUKIlFrcLDNssfkAN/OwbmUcU4aVPGMdqq8B8Inlz+cz85Hlz0cBnO/6ABHtAbAHANZvOidljCpYZRKUHOd9TuoZZaB1EquZV1BTubI+gcXEsEsDISRZqFZzKnaoEtF6AP8PwN9l5geJ6BFm3tr7+8PM7GVvK4dqyuS0Fqo1hhQysKgqWDJDVYJW7LexkNYOkny/oZytzCkaaLmTlmaFyA7aNnnTwmFEdBWAtzDzLy1/vx/AFcx8hIi2A7iLmS/2XcPVZs+H0rXfY1CqwFjJEgVaYYk+kilFKi07z1NlGtIWXWRv1ZIyJc/CAjny0oya0YRLttZx7q/FUyYZAPg0gGsA3LD8/47QBU5sjNvYKS34XLCIpZbGCksKjLUSJ6s1P2MmoCn4S6bifLXIgJ0iUuu+p2SsTqn0ACAkdyLaBOBKAL/We/kGALcR0ZsAHABwteRafXIOZZu6PlMiLlo7TKm/4IaLat/RS6ITmUKonagUQ+y1fSEliTylnHN/fiTaOnDy/Mdo7znrxqdAlYySkTycS9R3byFUUkTuzPzXAM4bvPZXWETPJKMj7TGSL03qob/Hpm6PkXo/1XlYWTC17kxtQu8gJfaapK5N6JYbWSOUTvPkGrLr125eD6SdvqTlB3z39OG0bLMn1d5LNOkdQ2q6dn+R9YsXDReQSyOQHruty8Vq1Id3ya7VmOcU5JrbNOobhRB6wI7dN9dJ20rjlQ5jWvvtj172JMG7kJOlWkvxKl5+4IWX/kvx+8fs61ZavGb69xAuLV5yzCu5cEqE1rVUy10rld3quK6RfTyc79D8xkTgSAvyWSNGjikmGYts1VjZNt+Jad3jMvIIOU01nXXSzkux8C24Fp0ypWKmQ8f2khUSc9AniZiHdUlHqOVJd4rEngKrvZpTPkSKKmYZrRDHDqVrzfQxhciKMdRMaqppZtOQ2bBGeO0HtkuWrnr8Y7VeYmPmffKLidXPQYocY2q7a5ysx9BEbRlNnNhI1Yl9SskvFhpeS1mqMQTfotw6orAwy2jX3neZa/r3sFgXFiRv8WAee48FYkk9Ry7VHapTw1TioIdoidT7sNLgW5SRZiN0wK3YpFR9TF0bUtmlkLyV/KRau+ShXaKVYs7emDy5x4ZeaWgTksQkjWv7ELtgTofaMiUJfWhzt3C6pdhlU8igVK0a6d4rWV7AGjVDlCdP7h1ie6dqHRldmzZlcZYIe2ylkbU2SpG6jyR8BJ9rXivVqam/b1JOABpoidg1MlJz+CVX5itD7n1YZeXV7uCu5WGP3cTaSCEFX9OHFPK0TFPXRg3tz9LZLZGlNJnPJ0dL2VlHPBXLUJ0iLGy5vsSbEr02O2hqctZEb0ESGj1AQ76T06H2SAsYyrKToXQv1PCBTYHYgfhmHZOBlebRSpq/RZiZ9pyVDnXccHD9Kf98KBVv3qJztyUc23n8yX+xsJ5b10N+19r3TO4rWbMxWEnN3YpULON2a9aSKdWuzZLstYnBpYFb1P0eQytKhCZiZdTKQ3EsSW04vtZi3ouXHwg16zj7gWPJsfC1mnZ0yFmM1mnpMShlhz/dks9arUHSh7YjtVTfgzFY2N2lD/KYIInQPKXUc2+O3KUoWSXydEQq0ZSq496KzLTLCVg3wUitV1Mi2awVTT0XuevBNQ93v+wDbdeWObG+jB1W0rcw5u8p90y9bgnSci2eU1Lp1xY/9xeq5PukanWxc2VVE8h3fReOHNyB/jkzdwwapgtXz4AnsTZOPmPfuVQt/tOl+UhKvfm7E+5TVHPfuOMCvvDN16pfV7rIahx7pSYOy0JMY4spJgVb2ihCglIaZw5yTFPWRbUkhO6Cy5wQS/RAO6emWmYsC1OTS379EsSX7DzStllmw0VrvP36tzVhW7SAhq26VK2VDr4FaOUH0Bh77Bry3TN0rZBcrTTbVDOFV3MfQezDuxWCH0JT4UgJudUqIz3Eey79rA25E9FWAL8P4FIADOCNAO4HcCuAXQD2A7iamR/2XacjdyDcHqyPWFNOKsmmmoy0HZC1qiWmJmkNiaE0KUhNNJLohpKnEdeYakNqhmuV3IfwJUwNkeo/KSFDS5v7hwD8L2beTUTrATwdwLsBfI6ZbyCivQD2ArguasRLaLbuSiHaXDJN1bYtTzAp87D5wFkATv4+U9jEw6SunDFP4fu2gKl00pI85Mf+1oIPIGcMQc2diM4G8BUAz+bem4nofgBXMPMRItoO4C5mvth3rb7NPZZ8SmmzNY7S1nVlck8WY2Sv9XDSIgkJ4YQ2r4YW1gIpxEJaz30VG52HzFi52rzGejAJhSSi5wG4EcA3Afx9APcAeDuAw8y8dfkeAvBw9/sYUkIhc0m9Vhx7TkMAi8qQUzIdtdbkIbe5Q87mLuGfsjJjWspRyx8xhpx96UOqPK3a7J0J4DIAb2PmLxDRh7AwwTwJZmYicq4QItoDYA8ArN90TszYAOTViKmpTXSLw7J4UalqgS7kNk72QbvJQ642LpFlLrFr9EstjdDe1JSjll27xL7sUFumEnI/BOAQM39h+fvtWJD7g0S0vWeWecj1YWa+EQvNH5u2XZA08hSCbyXVeSyN3ao8rAstkABQxrRm6dyKLTEwlROYZQBCi/VigFNlGdMHV4oSfVJ9CJI7Mx8lor8koouZ+X4AL8XCRPNNANcAuGH5/x2WA12ljj2na+VAa5JvKfJEm9gtK3am1PnXJvbasmthT2rLWBot8zYAf7CMlPkugF/BoqLkbUT0JgAHAFytOrIBrI72OfWiY+A6BroWlHY6e4uIac8WW8a3D0vCCDXH1pajRSetTg5WJ7sYs8zYXFnIsGT/1Jo5PZOoLaOl6ZWIhIm15eWYZ0o229WGRfEpCSybdKQ64aYmx1o+MGsnah8aDlVNm/s9N72z7QzV3MJhGiRfy7kj1dwBXbt76RDKHBmVdoBryRGwi65wYSoZ3i1GywBpmruFLGPkmBItU53cz37g2Cnviyn5W0sDLAFpq7Ec5GiNNedesrlX3bzlQkpZhtyHRQmF6XSSpUseVqGQalj3ODvJfAjfe4bEL3XSlUp7jyU/ybhy7be+h0Rs1mzJ2jcarQs7WIYkAvX66w6/1zDJTCKHlIqmpeVbSpY15Dgmw9yHblHNfcvmNX7+5W8Rv9+nwdduzLEqKKG5W8z7rLm7UdpkU8IfMjU55manTlJzP7GRinZZGkIrqcLS5q61kDU2ubTmT4cxGWk0tB5C4ySTglq+kzHUtr9rFheLnS9Nn0kfuftS2/Ga6kCfRA/Vlpx0KYX2Y6+vQQo5ttTUxRTbM1UyvtqnLYk8fUlpVrAiBAmk+3E4xtZLgbSIHDkWN8v89KvTmnVYZ6haIMdJlGvb0wrBK1mjfogWZKiBEr11SyVDTXkfaihNJWzuWmaZqjb3WBONpo3Xwr6n1bDCOlsvtvuOlhY4JWLQagZiKUvNevQaHcM6aDpbaxTrypFlSv2hPsaIfqVt7oBuZuP+Q9tEzRti0IKW4kNKS7XTFVqydK0zresOUbtQVX8cufMn3YsS52XJfZlqntN8UHeYhM09BZLmDbWccmPoxtP9r/3wAerUMkk1y2hXhwTK1DAZyqlmMpNVyYIQatjbrR+ioQd1rpy1lazq5F6iUqC0DZsPFmnrw6f88ChqHZs9htrp7VYoVXogJxROS5aWMkwp4mfxoJYg1ZnayXMoy5S+wjn9e3NQ1Obe78TkQozgLSelZD/MFmN4u7mtqb23mrouQWmZ1tDggbp5Dh20G1L7zCpacvWdyMdk2bzNHYhzJPomMyWrTorU6xw5uMP5euzirqXldHjqvk/Nw5AcJJnGff+KJJNY2u8yd5ONmb1yr5eK3GgayV7oz7sW0YdkGpKnb96kuREptvkhUuzkqkXM+qG0a099p0nZ3Hn96hz3YzZI13gaiNNgNcxJWhiSQ0fcPpLv/jYk+eEcaPZQjblmLMm3eMpKgWsNutZzat0nq8ghq6Q36bhq15yPRXWbex+1bMwp5qCUOtiaWbYdNMleOq/97/HYhRtPmochIQzJIJbYU6MIYiI2rEk7dCLVPkmk1IU5tvNUc9wPn7MheEKzelCXVGL68rAmcFebv/7JQXMtNkXuHTRDulyEmrNw+pmfrg42Zz9wLEhomtB80IUcZWPH//489L97iBhCcsj5biHzXykt3EUWY6YJjbBJaTkA1wOnW9sueQ7XdSypx8xDCmLnbex+VqGrQ5QoK1yF3C3icWPsftoYarKW47COpnC93n2/2CQn1xG+/wBJTaHPfWC3YF7xkYsLId9Dfy7Hirv1TzP908KYjffked5w0u/a1RNLmzzGHiy1TC+SLNVYFCV3Ok6nDFq7IYQlmWsWSUrF2H0s49ctYqU7xJq3fBqjliaoWfpB60Hi2vzdo3Mxf+MP6yHBdzi28zj2H9omWjuPXUhVs6ot0YpjHdA9iYvInYj2A3gMwAkATzDz5UR0LoBbAewCsB/A1cz8sO86647bxt9aE6xFw98SWmTLces5ae85kVeWZFQidV5DplMpl1EbpU56NZOY/iEzf7/3+14An2PmG4ho7/L361RHN4Ix23DpEEKNRZ1LBDWrA5aA75Q2dJpK5nAKROSTqWYNmBj4IrdS1+4UZGGFEiU/cswyVwG4YvnzzQDuQiFyB/zOP22St8psdGU11ooY8sFHFNoPD6vMR+sSsS4HWYzjVisRqUTGdwfpGkwt/6GdoCRFKHM8FrVqN0nJnQH8CRExgP/CzDcCOJ+Zjyz/fhTA+a4PEtEeAHsAYP2mczKH+xRq17ewgoQQYhZLSshmDbRUs78UYuQolWFKaYCU+7ig3azFtQ+moO2nyNXioSwl9xcz82Ei+lsA9hHRSY9DZuYl8Z+C5YPgRmDRIDtlkFpfXKNa3RSRS/DWD4dVrfeuaauNkWGNh3nJE0PLiNXSO7nGRIZJISJ3Zj68/P8hIvoUgBcAeJCItjPzESLaDuCh0HVOrHfHhrtgtVisCD7l+JfawSe2qTVgk36ei1XX1kOnsJTm5K3IrkML9WVKQXKqTmlubcV1wcJhRLQJwBnM/Njy530AfhvASwH8Vc+hei4z/4bvWqHCYVpIIT/f54coVZ87x9ZnceS3hNYCt4ho0oBmdUhfLkIOtJy1JYlcs+1ljKJVOlfiwBv26ndiIqJnA/jU8tczAfwhM7+fiM4DcBuAnQAOYBEK+QPftSzIPbUoVyz6tszhPbXDHlNJR8tWv0pkP0Sr5N/BomlDKqQ9cXOUIWnxsJR7pDbDbg37D22zIXdNbLhojXe98Z1q14spUBTKAMtNpoqtT2MVN1yqTWDpXpwlYVXrX4pUgk/VmKWROinkLpkbTS24f78xch8SeomonRxMhty3X/+2xc+B5IyUXo1DIbhs2jFpvjmZk1JIy9yOoXYqvZVWmXu0r6HtWpsjSlUJTUmckig6Fp3FJBgzb8bMZ+leuEM0T+6bL34Wn/fud5z02lhNjA5SbSH0xPaRvKW5ItU+qdXAe9UaU6TAguhTohtaTq6LOS3ExttbfO/cuvAhpCpdVvuheXLva+6j78lwfmpr7n1oxf/WihoouTBbqkMPtNsAvKUORanN06UKWQn5awdOuCA9fWhXHk0h9+ZK/m7/+cNRWnX/7/tx8mTuQ361tVRSdzlfu4UheehoQ/v6sfOYQrAahOC6r0ZUiGUnMA1oyTsU2pdiPtVESu16lxI4xFiN9Rj7vFVghBTFe6g+69++1fm30BNRUtLUCikx+d2iCi2kEnWdXdCuAz8Gbfmk1mkPmf+ANBNghxQn6JRiv61q66ciJfs8tBd3b7kXtz962ZO/t+J4nYRD1RcK6YtNt3RsShde7GLyCVqSvFRb0w7BIk1ditKkqGHiaKFcdCy00uhblVcoVDJH8dIOj56sWWaRfnvW4NU08ojR6mMW3fafP5w0HhcksbU5tclLmAdysiWl8dPaKBXS5jpB+Orwt6QV50aPuUImU7PCa4cg5qCFQIMmyN2XGm/doq71Xo0pi8T3nTSJv0S1SE35aPcq1UBK/97Qay2eBFqs67Tv6CXiWHjAjrCtChwWD4X8uf/8OtEk1XZO5WS+TiXkTSNMrGaz8iFyS9BqoiTp5ph8fCUuUsoQW7S7HFvnHTFbZJnW0rzH5vzb7712GmYZSQx3rkYTQsgJGG9v5N7vZ5l1hQfy6mnsW/7fNUQGTvVvxGxQzfIPJYi9e2+JBDXXfQHdRJm+HPv38N3T16zatQd9c1siWibUtFpagE8ze1szi3mM93IVpyqauwsxE18yKQUIV1Q8+4FjJ/0+bAwtjXOPqXutUaIgJaU91RdgmT4vvZ+Gqahk9A+QnswmIfW+EqCVDzK8pzZqtK/UUkhyEhObj5YZkrtFx+8Y+J6SPiIYkrkLfYKXkHspk0LKfcZIYSwWOPZeY/cNIaYglWa8u2SMteL6XZDEdXcIyVRzHCHk3l/blq15PamS1L3v7pd9YBpmGZ/dNzeGPdUhe2zn8ZPG4qs7/8PnbPAS/BixuxCKTnBFHvTHnHLdMUgWq7SEquWDWkISw4gVl4nB0lk/XE8SaGZYhuzUEmhmWaacNLXaFHb7KETOsfX3XUj9nhZoqnBYrCalkTwjLa6vWVcmpaD/2DUk0NZgXJughQJi0vr3mmPNlWXpPr8+SGVqsX60YH0Cto60Gnv4vufSz7Ztlhkr+atZP7xWKdlSWbPaR/eUWtup5ZJL26u1YPEw0IC0hnkorE/j+7UW5gjkRSzlkHhKVzYXdm+598mfL9l5pG1yH2aoajeFaKlGeGuL3RdnHEqnj80W7kxWQ8dyhxpdl1YVIYLviD1Hnh1WtaVeC3kPoQdCiuZe1Oa+7njYUTlGCCGkEMbUFmEsxsLi+ji287jIGZriYB6Tpy/c0oWSpxVLaFbd7ObEdc19vZ99kVGxypW010Io5DMFLRCwNvqy2+d53wKfjb6+mNyJaB2ALwE4zMyvJKKLANwC4DwA9wB4PTN7pbnucQ5GmoS0Pg20EKrVoXYNmdR08z5xS2VV8wHcAjnkNKuQPKifvC50TS59xMgwh+RbkFcuYnw/w4q2GojR3N8O4D4AW5a//w6ADzLzLUT0EQBvAvDh1IHkknlMByerdF9AXjlOg9StnJgSjc4nL03zmFbaekulB6xrBUnf76oPpCm7HLlpySvV/h2b9Zqb82HBSSJyJ6I1AP8EwPsBXEtEBOAlAH55+ZabAVyPSHK30M5THJslS7PWTtYC7B3Y2iGamglQUtJIIQUfIbRQSMqF3M5R1mUWUuoxAflOTWlpgxaL+nWQau6/B+A3AGxe/n4egEeY+Ynl74cA7HB9kIj2ANgDABs2bE0faQQ0tI9S5oBQXG3qYtB2VtfGqvtHxlC7xlJLCJUhKI1WH9gdguRORK8E8BAz30NEV8TegJlvBHAjAGzZvPYk4+Rq7S3UBNdYaCE7bGq9eUkugAZcx/pWQwcBucykx3JpXRPL2O5YSPNDUuVoUQEydq91MtE+gQ3HVDte3weJ5v4LAP4pEb0CwEYsbO4fArCViM5cau9rAMQFz6caEVNCa0h1uo0l0+TUXU9BTJPzEKbgTI2xzWr6W0LQqrXfRylnah+5srOoGNmHVhatBYLkzszvAvAuAFhq7r/OzP+ciP47gN1YRMxcA+AOw3GeFsSei1aP8JqEP3Vo+FxqmdxKt7c8HaEp25w49+sA3EJE/wbAlwHcJPlQitY+E3sYrRHBjHTkdkOyxkzsTyG3N4WlPKPInZnvAnDX8ufvAnhBzOdPbIxfFPNC8qN1IpjxFFomAinm/fgUWm86dEa1O7dIyCAAAB2ESURBVAux6K8av6hbNVFo49jO46OnlJY3YgtE1SJ8J84W5Bkrt9x92IoTOhUhC8JjF9Ip/7TQPLl3SCH53ObDHaa8wLQXjCZC8tR6QLciP6l5byb4k9GK/IaIkWeMmViL7JtokB2DmExUQM/2XivGVivcSjJfJbXpEu3Z+siRnyScLiZ8rhuPD75GMsO5q3EKit2HueGRLWUYd2ghIdGHyZF7DLRjbWstsBDBay2cEg+A0qTeR0y2Y2x8tLTUbv/+sWGugFvWrjktRfilI2gkXYtq3d+FnL2ZK8Oi5H6i4MPLMqOxJS2itEYQ6icr/Wxt+GSoVY87lMKeciqTZi3HJrFJWkcC/mi3WG1eG665zN2jqafm3H2p8XAuXvJXA9LiYKV6Up5OyFl0saV+LeGT4b6jl6gQvHZdEkBOGrFy0qjzlFunJhcW+zK1IXxOhy6tU9ekzDI1iaEmoUv6Ompq8JZHegnBW9QDB+JkGEvwMZmQFkf7Gnb3Gvux5j5MIXqL/SnFpMg9FlYk0SIkx/VWwg+lGrzm6SuFFKxT10NoldSBOGJvtV9pDixLSWiVDClrluk160g5BtY81lvY86SQ2mZDWkKOvVwbpe2zLflJYiNmfCRfqkDc8H5WaEE+MZBq8zkVXlNRtIfqls1r/PzL33LSa9ZFxGpp7VaLVLtIUSwpDB1vLdlqU47NtcnEIpzOiugtC8DVloMWYrovhdCX4z03vbPtHqp9aHRequ1Y7VBqYdYmdkC/wYqmJjgVQk9BbWK3Rmu12i0Qq7nnmmeKk7sWOdQyz9RYgFra3VQ3vhZaIQ+r5JdStftj955UwWrJdBYDyxozOQRfNs49oXCYBjS095nU9VHyAd0CYVhmNJau2W8pu+E8acjOIvKpVFmE1LkuanPftO0C/pmr3pF9nZpt9E4nE0wJ1MhYnfKD2ofSMi5ZituydMQYtEk+x6H67fde27bN/cT6/IgNLY0htbGv9dHRaqHEzlvJ9PXSJrYUGYZIQqM5dm40RenImRTZlQpPjiH13VvuPen32x+97JRrhEpJALLopxYbZE8CuXUuYgjfguRjM+CsFkqrWr42YmrN+JKaNIh9qkgheB+x1zafdcQ+RKiUBCAjedd3t9rHkzTLhKClCVon0Fhs/JqFinIwJrOSoawWxJIj4ymYZzTMataErlUrqI+UBLeczGQTswwRbQTweQAblu+/nZnfR0QXYdE/9TwA9wB4PTNn7URX8aLY6JrWiB0oq72Fxu0ijFY19dI5Cq2F4/m+fyixyVKmmjkmJeZ7SMQaZH/ls74VTfCxSl7u+g9q7kREADYx84+I6CwAfwbg7QCuBfBJZr6FiD4C4KvM/GHftVrU3CVZgMNJlmhULZc8aLXfagvae4da5rZctNaBrBWCHyKX4FOKwqWe4DccXG+jufOC/X+0/PWs5T8G8BIAv7x8/WYA1wPwkrsEEq0jlCUptb33F56P5LU3TGqlOS2M2es1NL6x+ZZcd/ie7lpj8z+Vss5aDVckGJuTlDUcax922d8lYcilS3tYRdCEZGxZBdQFkUOViNZhYXr5KQD/CcADAB5h5ieWbzkEYMfIZ/cA2AMAGzZs9d6nWxghkslNhPItthSHR45NMYfoUyN+fMgheN+DNOW6IWfd8PtraImxMkg5ntdAjAM+tD98pr1Ugh9C6uzWnP/dW+4ddaiO3WNsvWhEBeUqlVEOVSLaCuBTAH4TwMeY+aeWr18A4H8y86W+z1ubZWp2+elQozmAxslC2yQTG5ZnXWWwti1dW3NvzfwSQs6+s5JdTngrkJeZGmvqNY9zZ+ZHiOhOAC8CsJWIzlxq72sADsdcKwclkyf6KEEQKaVEU8IiS6Wqj6Fk4bfaxK6BWp19WmiuYuXo7pO3NKa9PyYJxsy9JR7OkmiZbQD+ZknsTwNwJYDfAXAngN1YRMxcA+AOy4F2WGViH95PW9ur7TCdIUcrtdx992iB+LXQimlNU6ZnCN6zHcCdRPQ1AF8EsI+ZPwPgOgDXEtF3sAiHvEltVCOoRey1oPlAOZ2JfRUTiVqQZwtjmAKkD2rt+ZREy3wNwM86Xv8ugBeojmaGGaxjn2f4EXsSi2m8Aqw+0baWgxBjPo1puqKJSZUfGAuV86E/oTlavEWlOi1I7XetdGKKWchatUhaIIecHpxDDGXum9McWYfWTOmOWlYYOlelZprhmootPdDBwgbfXPmBsx84lhXqmFPnQiNUrEYKe8mSA60Ufmsx+iIFWiYjjeqgGsXlrBPRtGVn1QRdgphom+arQvbhKjUQ+puE9GMLGaUQoy9u1ycwK1LJ6bAuCVkcm8+app6cGv0tn8JSIY2YGpNZykPXdS2rKp8tyCg2oiYEaaGxVOWteps9H8m7yFxz4ZRwtpZclCmnkA59kpfOcY7WnwoLmdUw2Who7JoNWqSEL800zpV1KXl0JK1ZXKxWrsoQ1W3uLpJfVVJPKRmbilRtvlXb6ZSjnoZwbf6acexjn9U6lQ1l14IWPoRGcbG+XFtIMqtO7j7NvQ/N416JptkuuGqCa2TBjUF7ga1KREaL5FKqkYNrD2k2+XBdf7jfptArNZbsWwy3LUru6x5nMZm7oE3wQHmNUDP7zYVaGmDMtWueDloglLGwyJTqox3GHN0xc51D8ikyHVvnJWTkI2vXHo0tGFayKccYqmvufUgdpn20YkbQWJCppK61aEq2ZwPKy64FYo9BqjbfmVqktZZcYZUxa0EakixVpKwc3lJTi1UxuJzAhxRUIffcqo4aaLEkQUpdmZY1dS2skr0dsOup6opc6WPMZOL7TAglHtBaDm8L52nLKNsgeyNVJ/ZcomiteFiMdjcFIrdCK1q7FrGnyHIsImYK60IzommsYJjr75JxSVG6bWJTZpkZ8Vh1Ym/xhNVBO+65j9BDWyvbuNV1MeVQyDGUjmCbPLnH2m9TGjz0IW0ikAMLz7tmRMRUUDKb0cJOKyUDLdNI7mkgNA5pEEPpU9ZQrtYNWErZ3idP7h1SIwRieqcOoUn0pRIfWqkv091f2gYxFpoEoanVpZRyLh0mmboupIqWz7laiti1ZJpTmtua5FeG3C2Qq+XHIHXTA20kTFhAo8SAliPudHHC1XaoWsXAx8ovRnvPIXjLvbty5F4jNFJrIaYukhxbXuvau6YZDUiXVYjgY47xsXK2fnhr1ZqJQejBrXEqLvlATtm7s1mmEZQKx8s95k1Ri49JTssxowF5mmENDd5CnlotEFOQso9SImVqnLRiSjo30WZvSmi1MFUMShzvamvrGqglJw1Hm1WcuwTSwl991NhXuafh3CiYXDmHlLQSipikh+oFAD4O4HwADOBGZv4QEZ0L4FYAuwDsB3A1Mz9sN1Q/cuq4twLLY93Ypk6t/FcDLcosFjE1vC2QYorLIXupzKwcqSmnLa1ImRw/mgaCzTqIaDuA7cx8LxFtBnAPgFcBeAOAHzDzDUS0F8A5zHyd71qSZh1SSBfYKhBCLCTkEJqXsWvkVB/UQE15hgioRPGo1M5iKQ+MFFJvJcKpFaSuieFc3P2yD+g362DmIwCOLH9+jIjuA7ADwFUArli+7WYAd2HRNHsUJ9bb2fNaIfHYTWQxbsk1fSWJ9x29BFhbvDYsY2pZelkydyUimMbGceTgjtHPlFp/kvK5zuSqtVOJRiN+XvN7t1hZ0RK+h9nwtHF3wvWj2uwR0S4AnwdwKYCDzLx1+ToBeLj7ffCZPQD2AMCZZ5/zc89+528mDNMNrYVVuzSudWsyIK6+fIf+8TSWGMaQ8p1KO4m1eo52sJYjECdLIF+etZUpzTUhyZEp4TNwyXD3lnsBAJfsPBKtuYvJnYiegcUD5P3M/EkieqRP5kT0MDOf47vGxh0X8IVvvjZmfMnI6UokhYQEYjpNWSX0uBDU3GFjH67xQI5tTKxleqpRxyi2Vkoqydcm9z6s9ndpJ3AfQzm+59LP2pA7EZ0F4DMA/piZf3f52v0ArmDmI0u7/F3MfLHvOj5yrxFrm4qxzR9bqz6241RLG2oMuaGKufeWprZLSE1K8iFSn4ItOYbkXXM8/I61TCwlT5RA+kmqD4kD18TmvjS53ATgvo7Yl/g0gGsA3LD8/46YG9cKydLA+DjcFS9rhJJpaTNTeKB0yBmr+7PT02K14AvVq9XJTIJWx5WD1AelJFrmxQD+FMDXAfxk+fK7AXwBwG0AdgI4gEUo5A9810qNlkkhx1UUshRTSWSqIaOUmvmuz/cxdUegpNH21E+UqaenGNla3uPAG/aaRMv8GYAxyb405mYpaDHkseVqkK2VALYsLRCLvtw0ZTh2LSvS13h4nyz/uFP0UKap4ZmxKFG6GTjZTBKTddpafZlqPVSlTTs0+6aOITbZQLNxQP+altAk87MfOBYlP0C/FZsLLdm2c+WZ24XJEr49GSrr25KMOlg2qR9CIlctWVYrPxBD8toE71pgsdqXVrnfKR7pU7pplZBhTZSUY+0M4pAcJQ7X1mElz5L1n84ochcPYiNMxiCdsNhFJnl/6kKYIrHXxq617zVFFPsPbVOXoy87uDaxhzAF+3sItfellhJUTXOP0f5iKwbmlBIt4UBxfb6kvb0WJNp7TCnYDjXJvmQIYAslmkMY24NWddotoL0v+5DuUQ1ZR2Wo5iIlWibnKWahRWjWlh6z9Vl2ZapNDhZFp2KQQy5jscul5TgVGabIz4L8tQqHlfaj9PHt915rl6GqgS2b1/inX52WoZp7VCmR4amBViMsrGuA1wyLtIClBt/6CW0KETMlHtRSSOQ5CXJ//uVv8b7HwlwzVTtgray/XPKY6nzXRm4MfodS5G8lZ0uzSAlojr+TZQq5N9esw+VgHSP82gkWw00Uc1zuj1kaXaCpaYYaCfQRIosWy/Bq1MaPzZD2zUOM7GLe6/qe2vKQ+LK0MPzuGu3rJIXBhtAK5ZS+f+w75sx5c5p7H2Okvuoae4daGsxUNfecFoW1YJXVaKm9T31faSNWhin9cyetufvMMS2SuYW9dih0LaLSz2g8FX0Z1SoeljJfJeSo9d4xlGjHF5JvacIvGXVjlZk6hLYcq5N7SkKMCyUbWFteu1bo1RA+MndlqErJv2bVyA7WxFDyBBGSr1ZkTchEVbqUREmUkKfFA7oquWsRe0mUPPZrHNlT2+l1m9dFDrFyK93RJ4QappvcB4rEJuuS5XDuNcg+x/eghVrmN6uTXkpXshCqkLsFqcfWJdGs7mdRQzpnEcU4cVxj6H8fH8lroTUtT1ueGoQgcRxLGtTkkr1V8IJvjnKIPEWWpeviuKLiNJoNFSX3ExvJtC5JDMEPi3/lLKDWyCmEse9aK366tfrgLY1lDL71WisRKkeOFsX4gGnIsoNLpjl7srrNXQppOFOM7U/jaKexIF3jqJ00AejaayVlB4YoEUetjdKyBMrZ3SXIKf8rmafS5Qv6Yxp7AJWuGinFZMi9j9gysn3U6GsZc83SR1AXLMggpSqktkYfI7vUNmlaiUgh1NLOY+SoWedds+xHCL567hZytToxV6kt00K7vJKNIcbQUky2ZqNoH1Kdq9YRC1YnsJIoqcHXDE9uKQxyCM3y3315msS5E9FHAbwSwEPMfOnytXMB3ApgF4D9WLTYezh0rRPr3YtibAFY2oB9WXca9z1ycMfo31p4sPQx5q33OVP7mcQaDnJJlq5lUtcUoqDG0I1d6hzvkEP2IQ0+pqZ77NzXfogCfhlKTwu7t9wLALj90csWLzzrWyedGvYj73tKzDIfA/AfAXy899peAJ9j5huIaO/y9+tSB1HSoZYaGthBQ/vZfOCs0Y3RogPIRQouQh+WjpCQvkWqvPSaVspDizJ0oS/XlHXtIngLO3ssYuSaG2EHpJl+SkBkliGiXQA+09Pc7wdwBTMfIaLtAO5i5otD19lw0Rpvv/5tkyls5GuaMIZYjdaa5GMdQCmykRJDqJ5ODWivxdjvZRX/HlsPqP9+a3lamdhyZBkTlmxJ5idp7r05MmmQPYLzmfnI8uejAM6P+XA3ka1rTmPjzNV4JNAq1qVdg8QVfyuJhc8ldqviVRrxxMPrxGAon1yyj103HdH213pInrmybLGFnevEV5LUgXAv11hkR8swMxPR6K4moj0A9gDAuvO2irThFhyuffgWzZDoh9p6SceThAB9i9+1mUNp5mOfc71nWIIgxd9hGTIZIvpamZcSW7WU1IYk5ruG1anSsp6Pb02lmIz6Y9139BIzgh/T2IF0paOoWSalExNg170ntl56Dc0uFpZmFW1I5BoKew3Nq1XmYwghAmu1YNgYrLJPLaEl3xLjDxF6yaqQnwZwDYAblv/fIfnQusdZ1BDbV5RKS6tPEVgrdmItpJK6q3BYyft3yE0Tr917VaPmfAmUJnattnip0VXW60Lb5zWGoOZORJ8AcAWAZwJ4EMD7APwRgNsA7ARwAItQyB+EbqbZiUmrxoVW0olVXe5YtJStmIPc4lStNWK27IubAql/JAaSfdeaXFIQ4yeRyH0l2uylmmXGYN2BqWaqcwpi7ektwcI30QKRWJB6Srbv8POxmAl/gVrmtEk369BGa8WorNF602QfWmzGogEpEaTU+pE4v32fj4UkcsYFSQkC7YghC0ylr20fkyV37SqRU4Zk4Tx2ITWjvadogVOTo4WzNiRDl2+qFZnHyM+qQmQt1FK8JkvuEkyNEFYdlkf7IaZEDjGNGoYEH5pTi2Yd0nvn7L/WCD4n8co6r2cMzZC7b6HELEhNQg8J0yq7MBaSxSOZQ+vkrFRbcUplzxpHfWmMekjuIXmmkmo/QcmFGJnH1IgC0uah//cWiD5lDDFd0bRRldylG12iffQnThq/nhu3ntO5yULQKRrCmAwsNb5cpPRhrVkYrB/at+/oJeL1KZFnVIEuTyEqF/FL7ezSMfjmIacJdQvEPwatsR1I+ExxcteIU48tVuQKd0yppdGHr/jX+Ge6+zCAMpmWQFhjk5gFcsi+tYzjPkKnnRLO3tCRX3rtEJH47uPKwh7a7TXleOWgAmJOxc8U4p/qwyIGRUMhN+64gC9887Wq1wzVgxgiNrW3ZE3zmj4C3+ZqvQZQKrRCR8dqrZTIUpXcJ3Q/TfnG7kcfanQ40hz/GGJKTHS4+2UfaDvOXZvcUwWR06swl+xLEbtF8wktErBun5ebOJIq4xDJW5iGYuRs+dB2afcxDzoJUuSac8LUIPoUmbvuk0LuTThUU7PlXEWQUrWlHKekq6RCqICYtdbafc8Wj5jaFR5dTT1SU7xjiX1sffa/o6W9XxpVYn0ac4Vp9ven5XrM+S4x+zKlx6vUv2IxL0XJfd3xuGSIzQc4+P4hUbgWkWRC+1XlYiofSmrl1Kxjrp18kZI8I722Zi30WG05xwzjI5cSuRYh4sxZA7Fhmr7PWiBE7K1lsZdUtprQ3FMwPAaOkbwLrgnua3sxNct9tXC6z9S2LafAR+rDwmHDOdEg+5TaMUMn3Rh8zkPAX6jONa5QiWWrWvRDxCg21vBFLaUSnOuBHVM2OiayaAytxd/7MOnaMkOkRtEAegV+pPcuhZSwQe1m2bXqkuQQm0/jttCQpdB4ULTQCc2qj6wErZBzjKmqeYeqqyqkRulYF2KTLGpDM7LCavOWzG4E0ouEtaBdlaxFEruuNdZHjHm1JMm3As0H/a6177XvUD2xkczIfAhXjHKp43EKNAlJq3XcEGN+kRSEUujHIi3GNk0t84MVUhOYYq4dun6HGBlLcgNWldD70IyQSr1GWXJfH9d9JwXWCTO5ZJnrjR8iRHypLfek0EhnHxKCL4ROmsa+KmRvrYhIiD6n6Jx2OGQfXcarz8+SEuGiCa2GIZPIUJUgVkO0aMOnmdzUH18KoUptv7G24FQNX9qIOLeEQfd9hp15hptZw4HYgh26NnzrQfoQ14hzl9SU6a8JlyN9LBvdNf+lslUl2rzm/Yra3DdctMbbr39b1kbKabkn2ViSBt4dUloGAuPjDkWd5GbjuiCRRc68+TCW1u4iBWkkTIecfrhase4h1CL6GCUgVq4aESk5iGkaXrsv7Eo5VDtyf/L3Bhr8hurOAKdudgmpdwiRe4hIYgleoz9jTJSIBiGE7q/dwjDmAT4GbfOfJdGHtF/XQ9OV9NRy6WVJ3SirIIsS37G4Q5WIXg7gQwDWAfh9Zr4h53qlIRXK0ObYJ+wxoh9zHA8XWK5JKUdD9V13iLG5cmWDSs02Y3BlHkvGIhlbi8QO6GjIHSTzEzJrxNyvFmIVC4sgi5adw8maOxGtA/BtAFcCOATgiwBey8zfHPvMhovWeNcb3xl1n9KNfC0KieWSQe7GkmaC5nj4LU9hsQ+fUiaZHLlakKXmCUfTpKFFgKltCyURPLk5FUPEmhFdn+/jPZd+tpxZhoheBOB6Zv5Hy9/fBQDM/O/GPuMqHKYVO51L9LnkpJmWD5TTlFKSgGLC56yil7QfeENo291Lml2GddNd6GSr/VCOTfwKoaSDPMZfMCRfCUJk77rm7i33AgAu2XmkKLnvBvByZv7V5e+vB/APmPmtg/ftAbBn+eulAL6RdMM28EwA3689iERMeezAPP7amMdfFxcz8+aYD5iHQjLzjQBuBAAi+lLs06clTHn8Ux47MI+/Nubx1wURfSn2M2dk3O8wgAt6v68tX5sxY8aMGZWRQ+5fBPB3iOgiIloP4DUAPq0zrBkzZsyYkYNkswwzP0FEbwXwx1iEQn6Umf888LEbU+/XCKY8/imPHZjHXxvz+OsievxFk5hmzJgxY0YZ5JhlZsyYMWNGo5jJfcaMGTNWEEXInYheTkT3E9F3iGhviXtqgoj2E9HXiegrKSFJpUFEHyWih4joG73XziWifUT0f5f/n1NzjD6MjP96Ijq8lMFXiOgVNcfoAxFdQER3EtE3iejPiejty9ebl4Fn7JOYfyLaSET/h4i+uhz/by1fv4iIvrDkoFuXQSDNwTP+jxHRX/Tm/3nBizGz6T8snK0PAHg2gPUAvgrgudb3Vf4O+wE8s/Y4Isb7iwAuA/CN3mv/HsDe5c97AfxO7XFGjv96AL9ee2zC8W8HcNny581YlOl47hRk4Bn7JOYfAAF4xvLnswB8AcALAdwG4DXL1z8C4F/UHmvk+D8GYHfMtUpo7i8A8B1m/i4zHwdwC4CrCtz3tAUzfx7ADwYvXwXg5uXPNwN4VdFBRWBk/JMBMx9h5nuXPz8G4D4AOzABGXjGPgnwAj9a/nrW8h8DeAmA25evNzn3gHf80ShB7jsA/GXv90OY0GJZggH8CRHdsyynMEWcz8xHlj8fBXB+zcEk4q1E9LWl2aY5k4YLRLQLwM9ioYFNSgaDsQMTmX8iWkdEXwHwEIB9WFgOHmHmJ5ZvaZqDhuNn5m7+37+c/w8SUbBf6exQleHFzHwZgH8M4C1E9Iu1B5QDXpz5phYD+2EAzwHwPABHAHyg7nDCIKJnAPgfAP41Mz/a/1vrMnCMfTLzz8wnmPl5WGTNvwBAennGChiOn4guBfAuLL7H8wGcC+C60HVKkPvkyxQw8+Hl/w8B+BQWC2ZqeJCItgPA8v+HKo8nCsz84HLR/wTAf0XjMiCis7Agxz9g5k8uX56EDFxjn9r8AwAzPwLgTgAvArCViLqkzUlwUG/8L1+ay5iZjwH4bxDMfwlyn3SZAiLaRESbu58B/BKmWdny0wCuWf58DYA7Ko4lGh0pLvHP0LAMiIgA3ATgPmb+3d6fmpfB2NinMv9EtI2Iti5/fhoW/Sbuw4Ikdy/f1uTcA6Pj/1ZPKSAs/AXB+S+SoboMm/o9PFWm4P3mN1UCET0bC20dWJRr+MPWx09EnwBwBRZlTh8E8D4Af4RFxMBOLJqpX83MTTotR8Z/BRYmAcYieunXevbrpkBELwbwpwC+DuAny5ffjYXtumkZeMb+Wkxg/ono72HhMF2HhfJ6GzP/9nIf34KFSePLAF631IKbgmf8/xvANiyiab4C4M09x6v7WiXIfcaMGTNmlMXsUJ0xY8aMFcRM7jNmzJixgpjJfcaMGTNWEDO5z5gxY8YKYib3GTNmzFhBzOQ+Y8aMGSuImdxnzJgxYwXx/wGIwCXlnBouCQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "U, V = np.meshgrid(\n",
    "    np.arange(Z.shape[2], dtype = int),\n",
    "    np.arange(Z.shape[0] * Z.shape[1], dtype = int),\n",
    "    indexing = 'xy'\n",
    ")\n",
    "W = np.array(\n",
    "    [Z_log_norm[:, :, i].ravel() for i in iter(range(int(Z.shape[2])))],\n",
    "    dtype = float,\n",
    "    order = 'F'\n",
    ").T.copy(order = 'F')\n",
    "try:\n",
    "    del i\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "fig, ax = plt.subplots(nrows = 1, ncols = 1)\n",
    "ax.contourf(U, V, W)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "del U\n",
    "del V\n",
    "del W\n",
    "\n",
    "del fig\n",
    "del ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEZCAYAAADCJLEQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3hUVfrA8e9J772QHkJCCy0h0puCihU7gl0Q66Lr2nZ1f7q6trV3RUXFhsK6VrDQVJCWkIReUkgP6b3PvL8/ZsCIAULaTML5PM88ydy5N/e9BvPec857z1EigqZpmqZZGxtLB6BpmqZpbdEJStM0TbNKOkFpmqZpVkknKE3TNM0q6QSlaZqmWSWdoDRN0zSrpBOUpmmaZpXsLB1AT/Lz85PIyEhLh6Fpmqa1kpSUVCIi/kdvt2iCUkrNBF4CbIF3ROSpoz6/G5gPtADFwI0ikmX+zADsMO+aLSIXnuh8tQ4+NF/wBA52NjjY2uBoZ2P63u737x3tbDk7NpDLR4dhY6O68Go1TdO0tiilstrcfryZJJRS1cAxdxARj04EZAvsB84EcoGtwBwR2d1qn9OBzSJSp5S6FZgmIrPNn9WIiNvJnDMkZpjMf3YpTQYjjS1Gmlp+/9rUYqTJYKSstonMklpGR3jz+MXDGNyvw5eoaZqmtYNSKklEEo7eftwWlIi4mw9+DCgAPgQUcBUQ1MmYxgBpIpJhPsdSYBZwJEGJyNpW+28Cru7MCYM8nfjXrGHH3UdEWJ6Uy5Mr93Ley+uZN6k/d06PwdXxlOoN7VWqGprZcKCEJoPxyDYREMT01XyLJYCXsz1nDA7QrWNN6wXa+1f3QhEZ2er9G0qpVOD/OnHuECCn1ftcYOxx9p8HrGz13kkplYip++8pEfmyrYOUUguABQDh4eEnDEopxeUJYcwYEsjT3+9l0S8ZfJuazyMXxnJWbL8THq/1nMLKBt7bkMnHm7OpaWxp93EzhgTw3BWj8HS278boNE3rrPYmqFql1FXAUkw3onOA2m6L6ihKqauBBGBqq80RIpKnlIoC1iildohI+tHHisgiYBFAQkJCu2fG9XZ14KlLR3DZ6FAe+nInCz5MYsaQAB65MJZQb5dOXpHWGQcOVbPolwy+TMnDYBTOHR7EteMj8XF1QClTEx9MNxsKzNsUSsGqPYd4/Ls9XPjqet68ejRDgnQXbl9RVN2Aq4Od7u3oQ447BnVkJ6UiMRUzTMSUoDYAd4nIwQ6fWKnxwCMicrb5/d8BROTJo/abAbwCTBWRomP8rPeBb0Vk+fHOmZCQIImJiScda7PByPsbDvLCqv2IwMLpMcyf3B97W12l31NEhK0Hy3nr53RW7y3Cyd6G2QlhzJsURbjvyd0wJGWVcdvH26isb+bJS4ZzcVxoN0WtdScRYWdeFav2HGLVnkPsyq/Cz82RV+fGMS7K19LhaSfhWGNQ7UpQ3UEpZYepSGI6kIepSGKuiOxqtU8csByYKSIHWm33BupEpFEp5QdsBGa1LrBoS0cT1GF5FfX86+td/Lj7EDEBbrwwexTDQjw7/PO0EzMYhZ92H+KtX9JJzq7A28We6yZEHmkxdVRRdQN/+SSZzZllXDs+gofOG4qDnb7hsHYNzQY2ppfy055DrNlTRGFVA0pBfLg30wb687/kPLLK6rj37EHcPCUKpfRYY2/QqQSllBoIvAEEisgwpdQITONS/+5kUOcCL2IqM18sIo8rpR4FEkXka6XUKmA4pgINMJeTK6UmAG8BRkwPG78oIu+e6HydTVCHrdp9iH9+tZPK+mbeuHo0Uwf+qXxf66T6JgP/3ZbL4vWZZJTUEubjzILJUVw2OgxnB9suOUeLwch/ftjHol8yiAv34vWr4gnydO6Sn611naLqBtbtLWbVnkP8eqCE+mYDLg62TInxZ/qQAE4fHICfmyMA1Q3N3P/f7azYUciZQwN59vKReqyxF+hsgvoZuBd4S0TizNt2isjxS+KsTFclKIBDVQ1c/95WDhyq5pnLR+huoi5yqKqBJRsP8vHmbCrqmhkR6smCKVHMjO2HXTd1qa7YUcC9y1JxdrDl5TlxTBjg1y3n0dqntrGFLZllrE8rYUNaCXsLqwEI9nRi+pBAZgwNZFyUD452bd+oiAiLNxzkyRV7CPF25o2rRjM0WI81WrPOJqitInKaUiq5VYJKEZFR3RBrt+nKBAWm8uablySxMaOUf5w7mJsm6y6FjtqZV8ni9Zl8sz2fFqNw1tBA5k+OIiHCu0f+m6YVVXPzh0lkltRy/8zBLNDdQz2mxWBke14lGw6U8GtaCcnZ5TQbBAc7G8ZE+jAx2o8pA/0YGuRxUr+TxINl3P7JNirqmvn3RcO4PCGsG69C64zOJqiVwB3AMhGJV0pdBswTkXO6PtTu09UJCqCxxcDdn6fy3fYC5k3qz4PnDtHP2LST0Sis2VvEO+sz2JRRhquDLZcnhHHDxEgifF17PJ6axhbuW57Kih2FnDOsHy/MHoWTfdd0J2p/VFnfzPc7C1i9p4iNGaVUN7SgFMQGezAp2p9J0X4kRHp3+r9/SU0jCz9N5rf0Uq48LYxHLozVv1Mr1NkEFYWpVHsCUA5kAlcdnnaot+iOBAWmP7SPfrub9387yIUjg3nm8hHH7H7QTHfMn27NYfH6TDJLagn2dOL6iZHMPi3c4uMFIsLbv2bwxIq9nDU0kNeviu+2rsVTTWOLgXX7ivkyOY/Ve4toajES4uXMlIF+TIz2Y8IAv04VvhyLwSg8/9M+XlubTmywB29cNfqkKz+17tXZBGUrIgallCtgIyLV3RFkd+uuBAWmP2xv/pzB09/vZWK0L29ePRp3Jz04e7SGZgN3fJLMqj2HGBnmxfxJ/Zk5rJ/Vley/tyGTf32zm0viQ3j2spG6VdxBRqOQmFXO/5LzWLGjgMr6ZvzcHLhgZDAXjQphRKhnj3Wlrtp9iLs/TwHgtavimRyji5usRWcTVDbwPfAZsEYsVZveSd2ZoA5bnpTL/f/dzuB+7rx3w2kEuDt16/l6k8r6Zm76IJGtWWU8ckEs146PsOpxnpdXH+D5n/Zz/YRIHr5gqFXHam0OHKrmf8l5fJWST15FPc72pkmYL4oLYVK0n8VapdmldSz4MJGcsjq+vH0iMYHuFomjr6usb+bDjQe5JD6UYK8TV8Z2NkG5AOcDVwLxwLfAUhFZf7KBW1JPJCiAtfuKuO2jbfi5O7DkxrH09+v58RRrU1TdwHWLt5JWVM1zV4ziwpHBlg7phESEJ1bs4e1fM/nLGdH87axBlg7J6tU3GfjrZyl8v6sQGwWTY/y5OC6EM4cGWs0MD4WVDZz/yno8nOz48o6JeOieji5TXN3Iu+sz+WhTFjWNLfz7omFcPS7ihMd12YO65odkX8I0BtWrBlp6KkEBpORUcOP7W1HAm9eM5rRInx45rzXKKq3lmne3UFLTyJtXj2ZKL3puTER44L87+Cwxh3+cO5gFUwZYOiSrVV7bxLwPtpKcU8Gd02OYOzbcansQtmSWMfftTUwbFMCia0brLtxOyi2vY9EvGXy2NYcmg5Hzhgdx67QBxAa3byKDDs1mftQPmArMBmYCicAV7T32VDQqzIvlt4znhve3MvutjdxxRgwLz4g+5Qbcd+VXct3irRiMRj6eP5a4cG9Lh3RSlFI8cclwappaeGLFXtyd7Jkz5sSTDp9qcsvruG7xFnLK63l9bjznDO/sYgfda0x/Hx46bwiPfLObV9emsXB6jKVD6pXSiqp5Y10GX6XkoRRcEhfKzVOjiPI/qZWQjqldCUopdRBIBj4H7hWRHpsotjeL8nfju4WTefirXby8+gC/HijmxdmjLFJCbQmbMkq56YNE3J3sWLJgPNEBvbO/39ZG8cIVo6htbOEf/9uBm6MdF/SCLsqesrewiusWb6GuycCHN45hbC+ZB++6CZFsz63khVX7GR7iyemDAywdUq+xPbeC19em88PuQhztbLh2fCQ3Tenf5TOxnLCLz7yw4IMi8miXnpl2rajrCCwBRgOlwOzDE9SaJ5edBxiAhSLyw4nO15NdfEf7JjWfB/+3A4NR+NesYVwaH9KnB91/3FXIHZ8mE+7jwpIbx7RroNTa1TcZuG7xFrZll/P2tQn6Dxrmm5Alibg62PH+jaf1ugU+65sMXPrGb+SW1/H1HZOI1OPFx1VS08g9y1JZt68Ydyc7rp8QyfUTIvE1TzXVUcfq4jthf5OIGDAVSHQpc+J7DTgHGArMUUoNPWq3eUC5iEQDLwBPm48diqlgIxZTl+Pr5p9ntS4YGczKu6YwLMSTe5alcsenyVTWNVs6rG7x+dYcbvkoiaFBHiy7eXyfSE4Azg62vHN9AoOD3LnloyQ2Z5RaOiSLWrGjgGvf3UKghxP/vW1Cr0tOYPqdvmUeg7rloyTqmtq/rtipJqO4hkte/42N6aXcN3MQvz1wBn87a1Cnk9PxtHdAZINS6lWl1GSlVPzhVyfPfWRFXRFpwrTW1Kyj9pkFfGD+fjkwXZmaHbMwVRE2ikgmkGb+eVYtxMuZT24ax30zB/HDzkJmvvQLG9P7zh8507Ng6dz33+1MivHn4/lj8e6GBy8tycPJng9uGEOotzPzPkhkR26lpUOyiA9+O8jtn2xjeKgny28ZT0gvvgkJ83HhlTlx7D9UzX3Lt9NLn6LpVklZZVz6xm/UNLbw6YJx3DYtukee82xvghqFqbXyKPCc+fVsJ8/d1oq6IcfaR0RagErAt53HAqYVdZVSiUqpxOLi4k6G3Hm2NorbpkXzxW0TcLK3Ze47m3hq5V6aWownPtiK1TW18LdlqTy1ci8XjgzmnWsTrKasuKv5ujny8fxxeLnYM/ftTTz27W7SimosHVaPEBGe+WEvD3+9i+mDA/l4/li8XHr/TcjkGH/uOXsQ324v4N31mZYOx6p8v7OAuW9vxtPZni9unUB8DxY6tesviIic3t2BdJeOrqjb3UaEevHdwkk89u1u3vw5nQ1pJVwaH0KItwvBXk6Eerng4WzXK8ap0opquO3jJA4U1XDn9BjunB7T58t2+3k68elN43jq+70s2XiQd9dnMqa/D1eNDWfmsH59cqqrZoORv3+xg+VJucwZE85js2L7VFXqrVMHsCO3kidX7mVosIee1R5YvD6Tx77bzagwL965NqFbu/Pa0t4qvkDgCSBYRM4xjwGNb88aTMeRB7SeXjjUvK2tfXLNCxx6YiqWaM+xVs/FwY4nLxnB1IEBPPTlDh755o/rLbo62BLs5UyIt7Ppq5czod7OTB8SiJuVtE6+Ssnj71/swMneliU3jjmlpo8J83HhtbnxlNQ0siwxl0+3ZHPn0hS8Xey5bHQoc8aEd1m5raUVVzdy51LTpKt3zTDdhPSGm6eToZTimctHcuC1DdzxSTLf/GVSr+667AyjUXh8xR7eXZ/J2bGBvDg7rsvWYTsZJzOb+XuYqvlGmpNFsogM7/CJ27ei7u3AcBG5RSl1JXCJiFyhlIoFPsE07hQMrAZizAUdx2TJKr4TERFKa5vIK68nv6KevMOv8nryK01fy81FFUGeTjx8QSxnxwZa7I9EY4uBx77dzUebskmI8ObVufH087TOhzJ7itEobEgv4ZPN2fy0+xAtRmF8lC9zxoZzdmxgr21VbUwvZeHSZKrqT41lK9KLa7jo1Q3093fl85vHn3Kznzc0G7j78xRW7Cjk+gmR/PP8odh2c4+IVa4H1Y4VdZ2AD4E4oAy4UkQyzMc+CNwItAB3icjKE53PmhNUe9Q1tZCaU8m/vtnF3sJqpg8O4F+zYgn17tmZmbNL67j9k23syKvk5ilR3HP2IKub7NXSiqobjrSqcsvr8XNz5I7TBzBnbHivSVRGo/Da2jReWLWfSD9XXr8qvldW6nXEj7sKWfBhEqMjvFk4PYYpMX59rsXYlvLaJm5akkhiVjkPnTeEeZP698h1dzZBrQMuBX4yrwc1DnhaRKZ2eaTdqLcnqMNaDEbe23CQ53/aD8BdM2K4cVL/HkkSP+4q5G/LUlHAc1eM4syhgd1+zt7MaBR+TSvh9bVpbM4sI8TLmb+eOZCL40K6/a60M0prGrnrsxR+PVDCrFHBPH7xcKvpVu4pnyfm8NyP+zhU1cjQIA9unhrFecOD+tS4W2vZpXVc//4WcsvreeGKUZw3oudmA+lsgooHXgGGATsBf+AyEdne1YF2p76SoA7Lq6jn4a92sWrPIQb3c+fxi4cxOqJ75vxrNhh55od9LPolg+Ehnrx+VTxhPnpNnfYSEX45UMIzP+xlZ14V0QFu3HPWQM6O7Wd1d+ZbMsv4y6fbKK9r5pELYpkzJszqYuwpjS0GvkrJ562f00kvriXMx5mbJkdx+egwi4zJdJe1e4u4Z1kqBhHevjahx+cO7fRkseYxo0GAAvaJSK97yrSvJajDftxVyCNf7yK/soE5Y8K5f+agLi39zSyp5d5lqSRmlXPNuAgeOn9Ir+mmsjYiwsqdhTz74z4yimsZGerJfTMHMzHa8hVjRqPw1i8ZPPvjPsK8nXntqvh2T/bZ1xmNwqo9h3jz53S2ZVfg4+rA9RMiuWZcRK9+1q++ycATK/bw4aYsBvdz57Wr4hlggcKezragLge+F5FqpdRDmJbc+LeIbOv6ULtPX01QALWNLbzw037e++0gXs72PHT+EC4a1bnplBqaDbz5czqvr0vH0daGf188jFmj2nzcTDtJLQYjX2zL48VV+8mvbGBitC/3nj2YUWFeFomnvLaJuz9PYe2+Ys4bHsRTlw7XC262QUTYerCcN39OZ83eIlwcbJl9Whg3Tuzf63oUduZVcufSZNKLa5k/qT/3zhxksRvPziao7SIyQik1CXgM00O6/yciY7s+1O7TlxPUYbvyK3nwfztJyalgcD93bp02oEP95usPlPDPr3aSWVLLBSOD+ed5QwjwOLWr9LpDQ7OBjzdn89raNMpqmzhzaCB3nB7NyB5MVFsPlnHnp8mU1DTxz/OHcPU4615I0lrsK6zmrZ/T+To1H4MI0wb6c/W4CKYNCrDq8UWDUVj0SwbP/7QPX1dHnrtipMVb8J1NUMkiEqeUehLYISKftK7o6y1OhQQFpn+AXybn8cbP6aQV1RDm48yCKQO4fHToCUtmi6obePy7PXyVkk+krwuPXTTslHq2yVJqGlt499dM3l2fQVVDCxOjfbl1ajQTo327LVnUNxl49sd9LN6QSZi36Zmu4aG6S+9k5VfUs3RrDku3ZFNU3UiwpxNzxoQz+7Qwq7upy6uo5+7PUticWca5w/vxxMXDrWImkM4mqG8xPat0JqbuvXpgi4iM7OpAu9OpkqAOO9xv/sbP6SRnV+Dn5sANE/tz9bgIPJ3/2H1jMAqfbMnmP9/vpbHZyK3TBnDrtAGn3DMgllbd0Mwnm7N5Z30mxdWNjAj15LZpAzhraL8unZ0jKauMe5ZtJ7OklqvHhfPAOUNOuSq9rtZsMLJ6zyE+3pzNrwdKsLNRnBUbyNVjIxg/oPtuNNrrq5Q8HvpyJ0YrXFGhK5Z8n4mp9XRAKRWE6QHaH7s+1O5zqiWow0SELZllvPFzOuv2FePmaMdVY8O5cVJ/Aj2c2JlXyYNf7iQ1p4KJ0b48NmtYn5kBobdqaDbwxbY83volnazSOqL8Xbll6gAuGhWCg13Hy5wbmg089+M+3lmfSbCnM/+5bITFu3f6osySWj7dks2yxBzK65qJ8nNl7thwLokPxaeHiyoq65t5+KudfJmST3y4Fy/OjiPc17rGy7qiii8emAQIsKG3FUjAqZugWtudX8WbP6fz7fZ87GxsGDfAl/UHivFxdeCf5w/lwpHBVnNXpZlatit2FPD6unT2FFQR5OnE/MlRzBkThovDybV4krLKuXd5KhnFtcwdG84/ztWtpu7W0Gxg5c4CPt6UTWJWOXY2immD/LkoLoQZQwK7vYdiW3Y5f/kkmcKqBu6cHsNt0wZY5XNcnW1B/R9wOfCFedNFwDIR+XcHg/EBPgMigYPAFSJSftQ+o4A3AA9MixI+LiKfmT97H5iKaXZzgOtFJOVE59UJ6nfZpXUs+jWdlTsKOXd4EPecPehP3X6a9RARft5fzOvr0tmSWYaTvQ0jQryIC/ciLtyb+HCvY453NDQbeP6n/bzzawZBns48fekIJsXoVlNP21dYzRfJuXyVnE9hVQNujnbMHNaPi+NCGBfl26WFFUaj8PavGTzzwz6CvJx4+co44npwFvKT1dkEtQ8YKSIN5vfOQIqIDOpgMP8BykTkKaXUA4C3iNx/1D4DATF3KQYDScAQEakwJ6hvRWT5yZxXJyitL0jKKuOb1AJScirYlV9Js8H0/3CIlzOjwr2IC/MiPsKb2GAPduVXce+yVNKLa5kzJpx/nDtYl49bmMEobM4s5cvkPFbuKKS6sYVAD0dmjQrholEhDAly71QvRlltE/csS2XN3iLOGdaPpy4dYfU3n51NUGuBi0WkwvzeC/hCRM7oYDD7gGkiUmAez1p3omSnlErFNHvFAZ2gNM2kodnA7oIqtmWVk5xTQUp2BXkV9QA42NrQYjTSz8OJpy8boasxrVBDs4HVe4r4X3Ie6/YV0WIUBgW6c3lCKJePDsPT5eQSy5bMMhZ+mkxZbRMPnT+Ea3rJIwMdSlBKqVcwjTmFA6cBP5k/moGpiu+SDgZTISJe5u8VpmXdj/ngh1JqDKaVdWNFxGhOUOOBRkwzmT8gIo0nOq9OUNqp4FBVA8nZFSTnlGOrFLdOG6BbTb1AWW0T3+0o4IttuSRnV+Bkb8NFo0K4elwEw0KOX/5vNApv/JzO8z/tJ9Tbmdfmxp/wGGvS0QR1nflbZ8AeU7JqwVRmjoh8cIxDUUqtAvq18dGDwAetE5JSqlxE2uwgPdzCAq4TkU2tthUCDpgWI0wXkUePcfwCYAFAeHj46KysrGOFrGmaZhV25Vfy4cYsvkzJo6HZSHy4F9eOj+Sc4X9eDLOkppG/mif2PX9EEE9e0vtmAelogrIHHse0rEUWpnn4wjGtDfWPjs7H194uPqWUB6bk9MSxuvOUUtOAe0Tk/BOdV7egNE3rTSrrmlmWlMNHm7I4WFqHn5sDs08LY+7YCEK8nPktvYQ7l6ZQWd+7J/btaIJ6AXAD7haRavM2D0xTHdWJyF0dDOYZoLRVkYSPiNx31D4OwErgGxF58ajPgszJTQEvAA0i8sCJzqsTlKZpvdHhZVs+3JjFmr2HADgt0oetB8uI9HXltaviGRLUe9fq6miCOgAMlKN2UkrZAntFJKaDwfgCn2NqjWVhKjMvU0olALeIyHyl1NWYWmq7Wh16vYikKKXWYFryQwEp5mNqTnRenaA0Tevtcsrq+GRLNl+n5DMuypdHZ8Xi2sufZ+togtovIgNP9jNrpROUpmma9TlWgjpR2t2tlLpWRJYc9cOuBvZ2ZYA9ISkpqcY8/tXb+QEllg6ik/rCNYC+Dmujr8O6tPc6ItraeKIWVAim2SPqMT0oC5CAqarvYhHJO6lQLUwpldhWlu5t+sJ19IVrAH0d1kZfh3Xp7HUctwVlTkBjlVJnALHmzStEZHVHT6hpmqZp7dGukTURWQOs6eZYNE3TNO0I65vWtnstsnQAXaQvXEdfuAbQ12Ft9HVYl05dR7uX29A0TdO0nnSqtaA0TdO0XkInKE3TNM0q6QSlaZqmWSWdoDRN0zSrpBOUpmmaZpV0gtI0TdOskk5QmqZpmlXSCUrTNE2zSjpBaZqmaVZJJyhN0zTNKlk0QSmlZiql9iml0sxLvx/9uaNS6jPz55uVUpHm7ZFKqXqlVIr59WZPx65pmqZ1L4utE2xeNv414EwgF9iqlPpaRHa32m0eUC4i0UqpK4Gngdnmz9JFZFSPBq1pmqb1GEsuZD8GSBORDACl1FJgFtA6Qc0CHjF/vxx4VSmlOnpCPz8/iYyM7OjhmqZpWjdISkoqERH/o7dbMkGFADmt3ucCY4+1j4i0KKUqAV/zZ/2VUslAFfCQiPza1kmUUguABQDh4eEkJiZ23RVYqdKaRmYv2kSotzOh3s6EebsQ6u1i+t7HBW8XezqR5zVN07qUUiqrre2WTFCdUQCEi0ipUmo08KVSKlZEqo7eUUQWYV6TJCEh4ZRYW6SxxUi0vxu5FXUkZ1dQWd/8h89dHGxbJS5novzdiA4wvQLcHXXy0jTNKlgyQeUBYa3eh5q3tbVPrlLKDvAESsW0iFUjgIgkKaXSgYFA328etUOwlzNvXjP6yPuqhmbyyuvJKasjt7ye3PJ6cspN32/OLKOmseXIvu5OdgxolbCizd+H+bhga6MTl6ZpPceSCWorEKOU6o8pEV0JzD1qn6+B64CNwGXAGhERpZQ/UCYiBqVUFBADZPRc6L2Lh5M9HkH2DAny+NNnIkJRdSNpRTV/eP28v5jlSblH9nO0s2FkmBfj+vswpr8v8RFeuDj01ga4pmm9gcX+wpjHlO4AfgBsgcUisksp9SiQKCJfA+8CHyql0oAyTEkMYArwqFKqGTACt4hIWc9fRe+nlCLQw4lADycmRvv94bPKumbSimtIL6phb2E1Ww+W8eraNIxr0rCzUYwI9WRMf1/GRvmQEOGNu5O9ha5C07S+6JRa8j0hIUFOhSKJ7lTd0ExiVjlbMsvYnFHK9txKWoyCjYLYYE/G9PchLtyLkaFehHo76/EsTdNOSCmVJCIJf9quE5TWGXVNLSRnV7A5o5TNmWUk51TQ1GIEwMfVgRGhnowI9WKk+au/u6OFI9Y0zdocK0HpQQStU1wc7JgY7Xeke7Cpxci+wmpScytIzalge24lv+w/gNF8HxTi5cyIUE9GhnmREOHNsBBPnOxtLXgFmqZZK92C0rpdbWMLu/KrSM2pIDXXlLSyy+oAcLC1YUSoJwmRpnGs0RHeeLs6WDhiTdN60kl38Sml/lzy1UpbzxxZO52grEdpTSNJWeUkZpWTeLCMHXmVNBtM/xZjAtxIiPQmIcKHMf19CPNxsXC0mqZ1p44kqBxAAAUEA9Xm792AfBEJa/NAK6YTlPVqaDaQmlNxJGElZpVT3WB6PivK35XpgwOYPiSQhAhv7Gz1JPya1pec9BjU4QRkngt/LZ8AACAASURBVCl8hbnsG6XUBcC53RWodmpysrdlbJQvY6NMM1kZjcKBoho2ppewem8RH/yWxdu/ZuLhZMe0QQFMHxLAtIEBeLro0nZN66tOOAallNohIsOP2rZdREZ0a2TdQLegeq+axhbWHyhm1Z4i1u4torS2CVsbRUKEN9OHBHDG4ACi/Nyw0bNdaFqv0+Eyc6XUj8Aa4CPzpquAGSJyZpdH2c10guobjEYhJbeCNXuKWLXnEHsLqwFwd7QjNsSDYcGeDA/1ZFiIJ/19XXXS0jQr15kE5Qf8C9PsDQC/AA+LSEmXR9nNdILqm/Iq6tlwoIQdeZXsyKtkT0EVjeZnsVwdbIkNNiWrYSEexAZ7EubjrKdp0jQroh/URSeoU0WLwUhacQ07civZaU5auwuqaGg2HtnHx9WBEC9n08vb9DXU2/R9qJcLHs52ehYMTeshHX5QVykVDdwNRLbeX0TO6soANa2r2NnaMLifB4P7eXB5gqnY1GAU0otr2FNQRV6FaUb3vPJ60opNE+PWNxv+8DPcHO3o5+lEkPnVz9PZ/NWJYE9n+nk64eGkk5imdaf29HMsxzRp60eA4QT7appVsrVRDAx0Z2Cg+58+ExHK65rJLa8jr7z+SAIrrGygoKqB/YeKKapu5OjOBhcHW/p5OhHi5UywpzPBXs4EezmZv5oSmp4lQ9M6rj0Jyigir3THyZVSM4GXMM1m/o6IPHXU547AEmA0UArMFpGD5s/+DszDlDQXisgP3RGj1vcppfBxdTDPHejV5j7NBiPF1Y0UVDZQUGlOXubv8ysaWFtYRFF145+O83NzIMjT1IUY4edCf19XInxd6e/nSoC7oy7g0LTjaE+C+sq8bPr/MC8SCJ2fSUIpZQu8BpyJabn3rUqpr0Vkd6vd5gHlIhKtlLoSeBqYrZQaimnpjVhMDxGvUkoNFBHdwtO6hb2tzZGWEXi3uU9ji4FDlY3kVdSTf/hV2UB+RT37i6pZs7eIJsPv42BO9jZE+roS4etCpJ8rkb6uDPB3Y1iIhy7i0DTal6Dmm7/+s9U2AcI7ee4xQJqIZAAopZYCs4DWCWoW8Ij5++XAq8rU6T8LWCoijUCmeb2oMZgWNtQ0i3C0syXc14Vw37anZjIYhfyKeg6W1nKwtI6DJbVkldaSXlzL2r3FR5KXjYKBge6MCvMyvcK9iAlw1ysaa6ecEyaobpzSKATIafU+Fxh7rH3MCxxWAr7m7ZuOOjakm+LUtC5ha6MI83EhzMeFyTF//MxgFAoq69l/qJqUnEpScipYubOQpVtN/4u4ONgyPMSTUWFejDS/gj2ddJGG1qe1p4rPGbgTiBCRW81VfTEisrLbo+sC5u7JBQDh4Z1t9Gla97C1UYR6uxDq7cIZgwMBU/HGwdI6UnMqSDG/3ttw8EhLy9PZniFB7gwN8jR9DfYgJsAdBzs9V6HWN7Sni28xsAOYbH6fDywDOpug8oDWrbNQ87a29slVStkBnpiKJdpzLAAisghYBKbnoDoZs6b1GKUU/f1MBRUXxZk6CJpajOwpqGK7+YHk3flVfLol+0iZvJ2NIjrAjaFBHgwN9iA6wI1gL2cCPXRZvNb7tCdBxYjIHKXU5QAiUqe65l/5ViBGKdUfU3K5Eph71D5fA9dhGlu6DFgjIqKU+hr4RCn1PKYiiRhgSxfEpGlWzcHO5kgX32EGo3CwtJbd+VXsLqhiT0EV69NK+CL5j/dsh8vi+3k4HXnGy/S9qcqwv58rzg66LF6zHu1JUE1KKSdMhRGYE0pTZ09sHlO6A/gBU5n5YhHZpZR6FEg0z57+LvChuQiiDFMSw7zf55gKKlqA23UFn3aqsrVRDPB3Y4C/GxeMDD6yvaSmkYziWgqrGiisrKewspHCKlOJ/Kb0Ug5VN2Iw/t6poBQEezozIMCNAf6misIof1ei/d3wd3fUrS+tx7VnLr6ZwAPAUEzdelOBeSKyuvvD61p6qiNN+53BKJTWmJ7tyimvI6O4lvTiGtKLa8gorqWu6fd7PndHO6IC3BgU6EZcuDfx4d5EB7jpykKtS3RoLj5zV14/TK2UCZgWLPxNRIq6K9DupBOUprWPiFBY1UB6Ua05YdWQVlzD7vwqyuuaAVPSGhnmRXy4F3ER3sSFeeHl4mDhyLXeqDOzme8UkWHdFlkP0glK0zrncGVhcnY527LL2ZZVwd7CKg73FEb5uxIf7m0qhw/1YlA/XVWonViHJ4sFUpRScSKS3A1xaZrWi7SuLLwkPhSA2sYWtudWsi27nOTsctbsLWJ5Ui5gKuoYEuTByFBPRoR6MTLUkwH+emFJrX2O2YJSStmZCxl2AYOAdKAWUzefiEh8z4XZNXQLStO6n4iQW15Pam4F23MrSc2pYGdeJbXmMS1XB1uGhXgyMsyLoUEeRPm7EunnioeTvYUj1yylIy2oLUA8cGG3RaVpWp+j1O8zZpw/wlRVaDAKGcU1pOZWsj23gtTcSt5v9dAxgJ+bI/39XMwtNDf6+7kS5e9KuI+LnhX+FHW8BKUARCS9h2LRNK2PsrVRxAS6ExPozmWjTV2DTS1GskprySipJbOklsxi09c1e4spqck9cqxSEObtwsBAtyNLpgwMdCfK31Unrj7ueAnKXyl197E+FJHnuyEeTdNOEQ52NkeS1tGqGpo5aE5c6eby9/2F1azbV0yLuSLD1kYR4evCIPPPGBTozuAgd/r7uuoxrj7ieAnKFnDD3JLSNE3rKR5O9owI9frT+lxNLUYyS2rZf6ia/Yeq2VdYzd7Car7fVXhkQUk3Rztigz0YEerJsBBPhod4EqmTVq90vARVICKP9lgkmqZpJ+BgZ8Ogfu4M6vfHVldDs4G0ItNzWjvyKtmeV8kHG7NoajGNcbk72hEb4sGIUK8jSSvCx0UnLSt3wjEoTdM0a+dkb6oMHBbiyRWnmeaRbjYY2X+omp15lezIq2THUYUZbo52DA3yIDbEg2HBnsSGeBDt74adrX5uy1ocr8zcR0TKejiebqXLzDXt1NbUYkpau/Ir2ZVfxc68SvYUVB+ZDd7RzobBQR7EBpuS1ohQTwb1c8deJ61uddJl5n0tOWmapjnY2RxpaR1mMAqZJTXszKtiV34lO/Oq+DY1n082Z/9+TLAHI80rHI8M9SLC10VPntsDTjjVUV+iW1CaprWHiJBTVs/2vApScypIzTF1Ex5uaXk625uWPQn1ZGSoF8NDPQnQM753WGemOuqOYHyAz4BI4CBwhYiUt7HfdcBD5rf/FpEPzNvXAUFAvfmzs3rrBLaaplkfpRThvi6E+/7+sHGLwciBohpTwsqtICWnktfWph2Zh9DX1YGhwR4MCfI4smBklJ+rHtPqBIu0oJRS/wHKROQppdQDgLeI3H/UPj5AIpCAaS2qJGC0iJSbE9Q9InJSzSHdgtI0rSvVNbUc6RrcU2BaMHJ/Yc2RQgwHOxsGBbofSVinRfowNNjDwlFbH6tqQQGzgGnm7z8A1gH3H7XP2cBPh8fClFI/ATOBT3smRE3TtONzcbBjTH8fxvT3ObKt2WAkvbjGlLDMqxz/uLuQzxJzuGZcBI9d1CcWh+gRlkpQgSJSYP6+EAhsY58QIKfV+1zztsPeU0oZgP9i6v5rsymolFoALAAIDw/vbNyapmnHZW9rw+B+Hgzu58HFcaZtIsKhqkaMp9CYf1fotgSllFqFabHDoz3Y+o2IiFLqZH9rV4lInlLKHVOCugZY0taOIrIIWGSOqVopte8kz2WN/IASSwfRSX3hGkBfh7XR12Fd2nsdEW1t7LYEJSIzjvWZUuqQUipIRAqUUkFAWwUOefzeDQgQiqkrEBHJM3+tVkp9AozhGAnqKPva6ufsbZRSib39OvrCNYC+Dmujr8O6dPY6LFVe8jVwnfn764Cv2tjnB+AspZS3UsobOAv4QSllp5TyA1BK2QPnAzt7IGZN0zStB1kqQT0FnKmUOgDMML9HKZWglHoHjjwo/Biw1fx61LzNEVOi2g6kYGppvd3zl6BpmqZ1J4sUSYhIKTC9je2JwPxW7xcDi4/apxYY3cFTL+rgcdamL1xHX7gG0NdhbfR1WJdOXccpNZOEpmma1nvoR5w1TdM0q6QTlKZpmmaVdILSNE3TrJJOUJqmaZpV0glK0zRNs0o6QWmapmlWSScoTdM0zSrpBKVpmqZZJZ2gNE3TNKukE5SmaZpmlSyaoJRSM5VS+5RSaeal34/+/G6l1G6l1Hal1GqlVESrzwxKqRTz6+uejVzTNE3rbhabi08pZQvsB87EtFruVmCOiOxutc/pwGYRqVNK3QpME5HZ5s9qRMTNAqFrmqZpPcBSS76DaZHBNBHJAFBKLQVmAUcSlIisbbX/JuDqzpzQz89PIiMjO/MjNE3TtC6WlJRUIiL+R2+3ZIIKAXJavc8Fxh5n/3nAylbvnZRSiUAL8JSIfHmiE0ZGRpKYmPin7UajkFZcw7ascnYXVHHW0H5MivFr10Vo1q+2sYXUnAq2ZZdTWd/M7adH4+XiYOmwtC4gIuSW17Mtu5yUnAqGBXty6ehQS4elnSSlVFZb2y2ZoNpNKXU1kABMbbU5QkTylFJRwBql1A4RSW/j2AXAAoDw8HAAKuubScmpYFtW+ZF/2NUNLQDY2Sg+3pzNfy4dof+h90IiwsHSuiO/223ZFewrrMJo7sm2UbBuXzFL5o0hyNPZssFqJ62h2cCOvMo//H6LqxsB0/+7LeabzfvOHoRSysLRap113ASllIo/3ucisq0T584Dwlq9DzVvOzqGGcCDwFQRaWx17jzz1wyl1DogDvhTghKRRZgXzQoaECtnPv8zacU1iIBSMCjQnfNHBBMf7kV8hDcB7o7c8lESf1uWSlltEzdNierEJWo9obCygS+Sc0k6WE5yTgVltU0AuDnaMSrMiztOjyYuwpu4MC92F1SxYEkSl72xkQ9uHEN0gB7GtGYiwoa0UlbtOURydjm78qtoMd9tRPi6MCnaj7hwL+LDvYkJdONf3+zmjXXplNY08sTFw7Gz1YXKvdlxiySUUmuP+SGIiJzR4RMrZYepSGI6psS0FZgrIrta7RMHLAdmisiBVtu9gToRaVRK+QEbgVmtCyza4hI8UK54/CPiw72JD/dmZJgn7k72f9qvscXA3Z+l8t2OAm6eGsUDMwfruzErlFtexxvr0lmWmEuTwcgAf1fT7zbC9PuNDnDD1ubPv7edeZVc/94WDEbhvRvGMCrMywLRa8cjIqzeU8Qra9NIzanA2d6WEaGeR363ceFe+Lk5tnncC6sO8PLqA8wYEsirc+Nwsre1wBVoJ0MplSQiCX/abskVdZVS5wIvArbAYhF5XCn1KJAoIl8rpVYBw4EC8yHZInKhUmoC8BZgxFQq/6KIvHui8yUkJEhbY1BtMRiFh7/eyUebsrl8dChPXqLvxqzFwZJaXl+Xxhfb8lAKLk8I49apAwjzcTmpn3HN4s2U1jTx5tWjmTLwT+OzmgUYjcIPuwp5ZU0auwuqCPNx5rZp0VwSH4KjXfsTzQe/HeSRb3ZxWoQPb1+XgKfzn29ENevR6QSllBoGDAWcDm8TkSVdFmEPOJkEBaa7sZdWH+DFVQeYMSSAV+fG67sxC0orqua1tel8lZKHva0Nc8aEc/PUqA6PJRVVNXDde1tJK6rm2ctHMmtUSBdHrLWXwSh8uz2f19amsf9QDVF+rtx2ejSzRgVj38Ebw29S87n78xQG+Lux5MYxBHg4nfggzSI6laCUUg8D0zAlqBXAOcB6Ebmsi+PsVieboA77cFMW//fVThIivHnn2tPwdNF3Yz1pT0EVr65JY8XOApzsbLl6XDg3TYkiwL3zf3CqGpqZ/0EiWzLLePiCodwwsX8XRKy1V7PByJfJeby+Lp3MkloGBrpx++nRnD8iuM3u2ZP164Fibv4wCV83B5bcOJb+fq5dELXW1TqboHYAI4FkERmplAoEPhKRM7s+1O7T0QQF8N32Au76LJkB/m58cOMYAvXdWLfbU1DF8z/t56fdh3BztOPa8RHMm9Qf3zbGHjqjodnAwk+T+XH3If5yRjR3nzlQjzl2sxaDkf9uy+XVtWnklNUzNMiDhdOjOWtoP2y6IDG1lppTwQ3vb8VGwfs3jGFYiGeX/nztzwoq63l59QEWTBnQrpuCziaoLSIyRimVBJwOVAN7RGRwB2K3mM4kKIANaSUsWJKIt6sDS24cQ5S/rgDrDtmldbywaj9fpuTh5mjHvEn9uWFC/25tubYYjDz4v518lpjDnDHh/PuiYV1yB6/9kYjw/c5CnvlxHxnFtYwM9WTh9BjOGBzQrTcF6cU1XPvuFirrm1l0zWgmROvnHLtDeW0Tb/yczvu/HQSBpy4dziXxJ35cp7MJ6nXgH8CVwN+AGiBFRG442QuwpM4mKIAduaYKMIAP541laLBHV4SmAcXVjby65gCfbMnGRilumNifW6cO6LEuVRHh2R/38dradGbG9uPVuXG6MKYLbUgr4T/f7yU1t5LoADfuPXsQZw0N7LHWamFlA9cu3szBkjpenRvHWbH9euS8p4K6phYWr8/krZ8zqGlq4ZK4UO6aEdPuwqUuq+JTSkUCHiKy/aQOtAJdkaAAMoprmPv2Zpzsbfh24WTcHHvF885Wq6qhmbd/yeDd9Zk0thiZfVoYC8+IoZ+nZbpR3/k1g39/t4e/zhjInTNiLBJDX7Ijt5L//LCXXw+UEOzpxF/PHMgl8aEWaaFW1jVz7XtbSDtUzXcLJxOpx6Q6panFyGdbs3lpdRolNY3MGBLIvWcPYlA/95P6OV1RxRcCRNDq4V4R+eWkorCwrkpQAFsyy7hy0UYujgvluStGdsnPPNU0NBv4cGMWr61Lo6KumfNGBPG3MwdaRdfpXz9L4auUPD6/eTwJkT6WDqdXSi+u4fkf9/PdjgK8Xey5/fRorh4XYfFK2PyKes556VcifV1YfuuEDlcJnsqMRuGb7fk89+N+ssvqGBPpw/3nDGJ0RMf+X+lsF9/TwGxME7kazJtFRC7sUDQW0pUJCuD5n/bz8uoDvHTlKF2ifBJEhGVJubzw034KKhuYHOPHfWcPZnio9QxeVzc0c97L6zEYhRV3TtbP0ZyEouoGXvhpP58n5uJoZ8P8yVHcNLl/mw/FW8rKHQXc+vE2bp02gPtn9qqhdIv7ZX8xT67cy56CKgb3c+f+mYOZNsi/U121nU1Q+4ARraca6o26OkG1GIzMXrSJ/YXVrLhz8kk9KHoqe2rlXt78OZ2RoZ7cP3Ow1Q5YJ2eXc/mbGzl7WD9enROnK/vaobSmkUve+I38inquGhvB7adH4+/etVWXXeXvX2xn6dYcPp431mr/DVqbr1LyuHNpCuE+LvztrIFcMCK4S6ouj5Wg2tu2zQCs5/bHStjZ2vDi7FEA3Lk0mRaD0cIRWb+PN2fx5s/pzB0bzpe3T7TqPwxx4d789cyBfLe9gGVJuZYOx+o1NBuYvySRwsoGli4YzyMXxlptcgL45/lDifJz5a+fpxyZv1E7tk0Zpdy7bDtj+/vw41+nMGtUSJc/EnC09iaoOiBFKfWWUurlw6/uDKy3CPNx4fFLhrMtu4KXVx848QGnsLV7i/jnlzs5fZA/j14Y2ytaJLdMHcD4KF8e+XoXGcU1lg7HahmMwl1LU0jJqeClK0cxOsLb0iGdkIuDHS/PiaO8tpn7lm/HktO+Wbu0omoWLEkk3NeFRdck9Ng4YnsT1NfAY8BvQFKrlwZcODKYy0aH8uraNDZnlFo6HKu0M6+S2z/ZxtBgD16dG99ryrdtbRQvzB6Fg50Ndy5NoalFt5Lb8sSKPXy/q5CHzhvKzGFBlg6n3WKDPbn/nMGs2nOIjzZnWzocq1RU3cD1723Fwc6W967v2Zl02vVXQkQ+AD7l98T0iXmbZvbIhbGE+7hw12cpVNY1Wzocq5JbXscN72/F28WBxdedhmsvK8vv5+nE05eOYEdeJc/9uM/S4Vid9zdk8u76TK6fEMm8Sb1vqqgbJkQydaA///52N/sKqy0djlWpa2ph3vuJlNY0sfj6hB4fZ29XglJKTQMOAK8BrwP7lVJTujGuXsfN0dRdUFzdyANf6O6Cwyrrm7nhva00NBt474bTeu2EnWfH9uOqseG89UsGvx4otnQ4VuPHXYX869vdnDU0kH+eP9TS4XSIjY3i2ctH4u5kx8JPk2loNpz4oFOAwSgs/DSZXfmVvDInjhGhPb8sTXv7WZ4DzhKRqSIyBTgbeKGzJ1dKzVRK7VNKpSmlHmjjc0el1GfmzzebHxI+/Nnfzdv3KaXO7mwsXWFEqBf3nD2IlTsL+WxrzokP6OOaWozc8mESB0treevq0QwMPLmH96zNQ+cNJSbAjbs/T6W0plcXtHaJlJwKFi5NZkSoFy9dGderp4byd3fk2ctHsu9QNU+u2GPpcCxORPjXN7tYtaeIRy6MZcbQQIvE0d4EZS8iR/o2RGQ/nazqU0rZYmqRnYNplvQ5Sqmjb8HmAeUiEo0pIT5tPnYopmmXYoGZwOvmn2dxCyZHMTHal399s5u0olN3UF1EeOC/29mYUcpTl4yw6mq99nJ2sOXlOXFU1jdz7yk+qJ5dWse897fi7+7Iu9cl4OxgFf/7dcq0QQHcOLE/H2zMYtXuQ5YOx6Le+TWTJRuzWDAlimvHR1osjvYmqESl1DtKqWnm1ztAZx8oGgOkiUiGiDQBS4FZR+0zCzg81rUcmK5MpV+zgKUi0igimUCa+edZnI2N4vkrRuFkb8PCT5NpbDk1uwteWHWAL5Lz+OuMgVw6+sSTRfYWQ4I8+Mc5g1mzt4glG7MsHY5FVNQ1cf37W2gxCu/fMKbNlW17q/vPGcSQIA/uXZ7KoaoGS4djESt2FPD4ij2cNzyIByz8EHN7E9StmGaRWGh+7TJv64wQoHU/WK55W5v7iEgLUAn4tvNYAJRSC5RSiUqpxOLinhk7CPRw4pnLRrK7oIr/fH/qDap/npjDy6sPcNnoUBZOj7Z0OF3uugmRnDE4gMdX7GFPQZWlw+lRDc0GFixJIresnrevTWCAFUxL1ZUc7Wx5Zc4o6psN/O3zVIzGU6uVnJRVxl2fpTA6wpvnrhjZ7c85nUh7q/gaReR5EbkEmA+s7i2zSojIIhFJEJEEf/+eW9Z7xtBArh0fwbvrM3nz53Tqm06NltT6AyX844sdTI7x48lLhveKZ51OllKKZy4bgaezPbd/vI0tmWWWDqlHGI3Cvcu3s+VgGc9eMZIx/fvmHIXRAe48fEEs69NKePDLnafMeGNmSS3zP0gkxMuZt6/tuWedjqe9VXzrlFIeSikfTGXmbyulOlskkQeEtXofat7W5j5KKTvAEyht57EW949zhzB1oD9PrdzL5P+s5e1fMqhrarF0WN0mNaeCWz9KIjrAjdeviu/Tk3D6ujnyypw4qhqaueKtjcxZtImN6X33GTgR4d/f7eGb1HzunzmYC0cGWzqkbnXlaWFcNz6CpVuzmfT0Wp5YsYeSPpyoCirruf69LSileO/60/BxdbB0SED75+JLFpE4pdR8IExEHlZKbReRER0+sSnh7AemY0ouW4G5IrKr1T63A8NF5Bal1JXAJSJyhVIqFvgE07hTMLAaiBGR4zZTunouvvbaklnGS6v3syGtFF9XB26aEsU14yJ63fNAx7P1YBk3vLcVb1d7PlswnmAvZ0uH1CPqmwx8siWbN39Op7i6kTH9fbhregzjB/j2mdaj0Sg8/PUuPtyUxfUTInn4gqF95tpOJK2omlfXpPF1aj4OdjZcPTaCBVOjCHDvnY9LtCWnrI6572yivLaZJfPGEB/e87OAdMWS72dhKlh4UES2djZBmX/uucCLgC2wWEQeV0o9CiSKyNdKKSfgQyAOKAOuFJEM87EPAjcCLcBdIrLyROezVII6LPFgGS+tPsCvB0rwdrFn/uQorpsQ2evXk9qQVsL8DxIJ8nTi45vGEuR5aiSn1hqaDXxqTlSHqho5LdKbhdNjmBTt16v/mBuMpmrMZUm5LJgSxd/PGdyrr6ej0otreG1tGl8m52Fva8PcseHcMnUAgb30ub7DMopruOqdzdQ1GVhy4xhGhvX8s07Q+QR1OfBPYL2I3KaUigKeEZFLuz7U7mPpBHXYtuxyXl59gHX7ivFysWf+pP5cOyESDytajqC91u4t4uaPkujv68pH88da9eSgPaGh2cDniTm8sS6dgsoG4sO9WDg9hqkDO7ccgSU0G4zc/Xkq36Tmc+f0GO6aEdPrrqGrHSyp5dW1afwvOQ9bG8Wc08K4ZdqAXnlTtq+wmqve2YyIWHx18C5bUbc3s5YEdVhKjmmC2TV7i/BwsuPmqQO4YWIkLg69o0X1/c4C/vJpMoP6ufPhjWPxtpJ+a2vQ2GJgWWIur69NI7+ygVFhXtw3cxATBvSO58EaWwws/DSZH3Yd4v6Zg7l12gBLh2RVskpreX1tOv/dlouNUswZE8Zfpsf0mpL7nXmVXPPuZuxtbfjkprFEB1j2IfrOtqCcMD00GwscadOKyI1dGWR3s7YEddiO3EpeXLWf1XuLCHB35M4ZMVyREGbVRQZfpeRx9+epjAz15L0bxugF/Y6hqcXI8qRcXllzgILKBqYM9Oe+swcxLMR6Fmc8WkOzgVs+SmLdvmIeuWAo10/sffPr9ZScsjpeW5vGsqRcnA4vzjglyqq77bdll3Pd4i14ONnz8fyxVrHsfWcT1DJgLzAXeBS4CtgjInd2daDdyVoT1GFbD5bx1Mq9JGWVE+Xnyj1nD+KcYf2srlvls63ZPPDFDsb29+HdXjj5qyU0NBtYsvEgr61Np7K+mQtHBnPPWYMI97WuRS5rG1uY/0EimzJLeeLi4cwZE27pkHqF9OIanv1hHyt3FuLr6sAdZ0Qzd2w4jnaWL9VubVNGKfPe34qfuyMfzx9LqLd1/PvrbII6bNABSwAAIABJREFUXMW3XURGKKXsgV9FZFx3BNtdrD1Bgamcd9WeIv7z/V4OFP1/e/cdH3V9P3D89c6ChBBICBlkkABhzxCGMhQVBNyCWPdA0V9btbZ1ddj+2mq11Wq1LhT3goJWq6KAgyVIBlN2FglkEJKQhJD9/v1xF8wPE3LJ3eXuks/z8bjH3X3ve/d9f/iGe9/3MyvcbtXZN77N4g8ff8/0wX156frxnWKKm450/GQtL61N59WNmdQ3KNdOjHWbqqGyqlpufS2ZtEMlPLlgDFeM6zwzgHSUbTmlPLZyD5sziokJ8edXM4dw6RjHrDprr3X7j7LoLcs4p3dvn+xWHTzsTVBbVHWiiKwDfgrkA1tUdYDjQ3UeT0hQjeoblBVpuTy1ej95x6uYlhDKA7OHurRq6KW16fx15V5mDg/nX9eOc7tfh56koKyKp9ccYFlKDt18vLjdxVVDpZU13PjqFnYfKeOZa8Yxd5TnrOnkblSVtfuP8vjn+9iTV8awyCDunz2Ec13YUWbN7gJ++k4aA8MCeWuh+01PZW+Cug1YAYwGXgMCgYdV9UVHB+pMnpSgGjVXNXTXeYNI6MCZwVWVf355gKfXHODi0ZE8dfVYt24f8yTpRyt4ctU+PttpqRr66YxB/GRCTIdWmxZVVHP9K9+RUXSCF69P5Lyhrpm5urNpaFD+u+MIT6zaR07xSSbFh/DLmYOZGB/SoYnq0x153PP+Vkb0C+KNWyfSO8D9OjOZXnx4ZoJq1LRqqKq2gemD+3Lb1HimJTh3nE11XT2/+3AX/07NZf74aB6fN9qjl1VwV9tySnl85V42ZRwjqLsP10yK5aaz4pw+4Hn3kTIWvZVCUUU1r9w4gakJ7lGV3JnU1DXw7nfZPPvVQY6dqGFUVC8WTo3notGRTv2hp6q8sDadv3+xj6T+wbx68wR6uulQFnuvoMKBR4F+qjrHutzFWaq6xPGhOo8nJ6hGxyqqefe7Q7yxKZuiimoGhweycGo8l42NcvjcWQVlVdzxVqpl3Z/zBvGLCwa7RV16Z5aaXcKrGzJZuSsPEeGiUZEsnBrvlAGUn+7I49f/3k4vf19evGE8Y100SLOrOFlTzwdbc1myIZOMoyeICOrOjWf359qJsQ6/qqmsqeP+5Tv4ZEceF4+O5O/zx7h1e7G9CWollqq936rqGOs0RVtVdZTjQ3WezpCgGlXX1fPf7Xks2ZDJnrwy+vTw4/rJ/bl+cn+HDJZNO1TCnW+lUlFdxz8WjGH2SNMm0ZFyiit549sslibnUF5dx4S4YBZOjWfm8Ai7r2AbGpQnV+/jua/TGd8/mBeuT+xUU/e4u4YGSxvVKxsy2HjwGP6+3swfH80tU+IY4IDZ4XOKK1n0Vip788u4/8Kh3HnOALfrCXw6exNUsqpOaOzNZ922TVXHOiFWp+lMCaqRqrIp/RhLNmTy5d5C/Ly9uGxsPxZOi2doRPtGhi9LyeF3H+4iold3Ft84vt2fY9ivvKqWZSm5vLYxk9ySk8SE+HPL2fEsmBDTrg4VZVW1/OL9bXy1t5BrJsbwx0tHmM4uLrQnr4wlGzL5eNsRahsaOH9oGLdOjeesAe2by3FT+jF+9m4atfUNPHPNOGYMCXNC1I5nb4L6BpgHrFbVRBGZDDyuquc4PFIn6owJqqn0oxW8tjGT5am5VNU2MCEumAVJMVw0OtKm2Slq6xt45NM9vP5tFlMHhfKva8e5ZYNqV1RX38Dq3QUs2ZBJSnYJPfy8uXRsPxYkxTA2prdNX2bpRyu4/c0UDh2r5A+XjuD6SbFu/8u6qygsr+LtTdm8/d0hik/UkBAWyNUTYrhiXBR9bOhxp6q8uSmbP32ym7g+Abx8Y5JDrsY6ir0JKhF4FhgJ7AL6AvNVdYejA3Wmzp6gGpWcqGFpSg5Lk3PILDpBYDcfLhkTydUTYhkT3avZL6XiEzX87J00NmUcY+HUeB6aMxQf01PPLW3LKeXtzdl8uiOPk7X1DA4PZEFSDFcmRre4TMLXewu5+72t+Pl48fx1iUwa0KeDozZsUVVbz0fbDvPelhy25ZTi6y3MHB7OgqQYpiX0bbZ6t7qunof/8z1LU3K4YFgYT1091m07Q7SkXQlKRCYAOaqab213ugPLldRuLN3M27VSm3VdqaVAHJAFLFDVktP2GQu8AAQB9cAjqrrU+trrwDlYVtgFuFlVt7V23K6SoBqpKslZJSxNzuHTnUeoqm1gSHhPFlh/mTV+mTX25Cosr+avV4zqVEu0d2blVbV8siOP95Nz2G79Mps1PIIFE2KYOigUby9BVXn+m3SeWLWP4ZFBLL4xiagushSKp9tfUM7S5Bw+3HqY4hM19OvVnfnjo7kqKYaYEMsMEIVlVdz5dipph0q567xB3OuhHZnam6DSgAtUtVhEpgPvA3cBY4Fhqjq/ncH8DShW1cdE5EEgWFUfOG2fwYCq6gER6YdlocRhqlpqTVCfqOrythy3qyWopsqqavnv9iMsS85he+5x/Ly9mDkinDHRvXhq9QF6+fvy0g3jXTbdvmGfffmNX2a5lFTWWr7MkmJIP1rBpzvyuHRMPx6fN9qte3IZzaupa2DNngLeT85h/YGjqMKUQX2YOSycF9amU15Vx5NXjWGOBw+ubm+C2q6qY6yPnwOOquofrc/b3UlCRPYB56pqnohEAt+o6pBW3rMdS7XiAZOg7LM3v+zUL7PSylqS+gfzvOnJ1SlU19WzencBS5Nz2HCwCAEemD2URdPdvyeX0brDpSdZnpLLspQcDpdaOs28fGOSx3dkam+C2gWMVdU6EdkLLFLVdY2vqerIdgZTqqq9rY8FKGl83sL+E7EsljhCVRusCeosoBrLaroPqmqz6zGLyCJgEUBsbOz47Ozs9oTcKVXX1bMnr5zhkUH4+Zj2ps7mcOlJauoaiHeD2aoNx2poUL4/UkZcaIDHtTc1p6UE1VrXrveAtSJSBJwE1ls/bBA/tP+0dMA1QEQzL/226RNVVRFpMUtar7DeAm5S1Qbr5oewzAfoBywGHsAyy/qPqOpi6z4kJSV1nWkzbNDNx9sMzuzETFtT5+XlJYyKdt8lWxyl1V581i7lkcAqVT1h3TYYCFTVtHYd1MYqPhEJAr4BHm2pOk9EzgV+raoX23DccmBfe2J2M6FAkauDsFNnKAOYcrgbUw73Yms5+qtq39M3tjo4RlU3N7Ntv22xtehj4CbgMev9R6fvICJ+wIfAm6cnJxGJtCY3AS7H0vXdFvuau4z0NCKS4unl6AxlAFMOd2PK4V7sLYerGh4eA2aKyAHgAutzRCRJRF6x7rMAmA7cLCLbrLfGThnviMhOYCeWDP2Xjg3fMAzDcDaXLD6jqseA85vZngLcZn38NvB2C+8/z6kBGoZhGC7X1bpuLXZ1AA7SGcrRGcoAphzuxpTDvdhVji61HpRhGIbhObraFZRhGIbhIUyCMgzDMNySSVCGYRiGWzIJyjAMw3BLJkEZhmEYbskkKMMwDMMtmQRlGIZhuCWToAzDMAy3ZBKUYRiG4ZZMgjIMwzDckklQhmEYhltyyWzmrhIaGqpxcXGuDsMwDMNoIjU1tahdCxZ2JnFxcaSkpLg6DMPBTlTX8emOPOaPj8bLS1wdjmEYbSQi2c1tN1V8hsd7cW0696/YwdoDR10dimEYDuTSBCUis0Vkn4gcFJEHm3m9m4gstb7+nYjEWbf3EZGvRaRCRP7V0XEb7qOqtp53vjsEwOc7810cjWEYjuSyBCUi3sBzwBxgOHCNiAw/bbeFQImqDgKeAh63bq8Cfg/8uoPCNdzUf7YepvhEDQP69mDV7nzq6htcHZJhGA5yxgQlIuUiUtbMrVxEyuw89kTgoKpmqGoN8D5w2Wn7XAa8YX28HDhfRERVT6jqBiyJymjGvvxyKqrrXB2GU6kqr27MZHhkEPdfOISSylq+yyx2dViGYTjIGROUqvZU1aBmbj1VNcjOY0cBOU2e51q3NbuPqtYBx4E+bTmIiCwSkRQRSTl6tGu0Uew+UsaFT6/j0x1HXB2KU60/UMT+ggoWTo3nnMFh+Pt6s3JXnqvDMgzDQdpUxSciYSIS23hzVlCOpKqLVTVJVZP69v1RL8ZOaVhkTwaE9mBF6mFXh+JUSzZk0rdnNy4eE4m/nzczhvbli+8LqG9QV4fmMY6WV/OP1fu5672tvLI+g9TsEqrr6l0dlmEANnYzF5FLgSeBfkAh0B/YA4yw49iHgZgmz6Ot25rbJ1dEfIBewDE7jtkliAjzxkfz9y/2cehYJbF9AlwdksMdLCxn7f6j/HLmYLr5eAMwZ2Qkn+3MJzW7hInxIS6O0L0dLKxgyYYMVqQdpra+gfCe3fnvdssVt5+3F8P7BZEYG0xi/94kxgbTr7e/iyM2uiJbx0H9GZgMrFHVcSIyA7jezmMnAwkiEo8lEf0EuPa0fT4GbgI2AfOBr1TV/Dy2wRXjonhi1T5WpOVy78zBrg7H4V7dmIWfjxfXTfrhQn7G0DD8fLxYuSvPJKhmqCpbMot5eX0Ga/YU0s3Hi/njo1k4NZ6BfQMpLKsi7VApWw+VkHaohHe+y+bVjZkARAR1J7F/b0ZF9WZQWCCDwgKJCfbHx9uMVDGcx9YEVauqx0TES0S8VPVrEXnangOrap2I/Bz4AvAGXlXV70XkT0CKqn4MLAHeEpGDQDGWJAaAiGQBQYCfiFwOzFLV3fbE1Jn06+3P2QP78MHWXO45P6FTDWAtOVHDB2m5XDkuij6B3U5tD+zmw/SEvny+K5/fXzTcLcqsqmxKP8aBwgpuPKs/Ih0fU119A59/n8/L6zLYnnuckB5+3HN+Ajec1Z/QJv9+YUHdmT0ygtkjIwCoqWtgT14ZaYdK2HqolLRDJXzWpCu/n7cXcaEBDAoLZGDfwFP3A/r2IMDPfeYAqG9Q9heUk5JVTHJWCanZJVwwLIz/vWykq0MzWmHrX1GpiAQC64B3RKQQOGHvwVX1M+Cz07Y93ORxFXBVC++Ns/f4nd28xGh+uWw7yVnFTBrQpr4lbu3dLYeoqm3glinxP3pt7qgI1uwpYHtuKeNig10QnYWqsv5AEc98eYCU7BLA8oV/+/QBHRZDVW097205xJINmeSWnCQ+tAd/uXwk8xKj8ffzbvX9fj5ejInpzZiY3twyxbLt+Mla0o9WkF5YwUHr/Z68cj7flU/Tpr+pg0J57tpEegX4Oql0LauqrWd7Tikp2SUkZxWTml1CeZWlR2t4UDdCenTjzc3ZXDupP0MienZ4fIbtbE1Ql2Hp0n0vcB2WtqA/OSsowzFmj4zg9//ZxYq03E6ToGrqGnhzUxbTEkKb/XI5f1g4vt7C57vyXZKgVJVv9h3ln18eYFtOKZG9uvPny0aw4WARj32+l8T+vRnfv2OqH//y6W7e3nyIpP7B/P7i4VwwLBxvO68qe/n7WtqmTvu3ra6rJ/tYJQcLK9h9pIzF6zK4evEm3lw4kbCe3e06pq1W7szj5fUZ7Dx8nNp6S7YcHB7IJWP6MSEumKT+IUQH+1NaWcu0v33NP1bv46UbkjokNqN9bEpQqtr0aumNFnc03EqAnw9zRlk6DvzvpSNt+tXs7j7bmUdBWTWPzRvd7Ou9/H05e2AoK3fl8+CcoR1WpaaqrNlTyDNfHmDn4eNE9fbn0StGMW98FN18vLl0bBSXPLuBn7+7lU/vnkZIDz+nxnOypp7/bD3CleOi+MfVY516LIBuPt4MDu/J4PCezB0VyeQBfVj0VgpXvbiJtxdOIibEuR11ausbeOjDnfTy92Xh1AFMiAtmfP9gegf8+N85uIcft08bwFNr9rM9p5QxMb2dGpvRfja1cJ42YLdKROodMFDX6ADzEqOpqK7ji+89fxogVWXJhkwG9u3BOQktDxmYOyqCQ8WVfH/E+X+iDQ3K57vyuOiZDdz+ZgrHT9byt3mj+ea+c7l2UuypHoa9/H15/rpEjlXUcO/SbTQ4uSv8qt35VFTXMT8p2qnHacnUhFDeuW0SpZW1zH/xW/YXlDv1eN9lFFNaWctDc4bx4JyhnD8svNnk1OjWqXEEB/jyxKp9To3LsI9NCarpgF3AH5gHPO/UyAyHmBRvqdZYkZbr6lDslpJdws7Dx7llSvwZO0DMHB6Bt5elms+ZDhSUM/eZ9dz5dhona+t54qoxfPmrc1gwIQbfZnq3jYzqxcOXDGft/qO8sDbdqbEtT80lqrc/k+NdV7U7LjaYZXechSoseGkT23JKnXasz3blEeDnzblDbBvr2LO7Lz89dxDrDxSxKd2MXHFXbe4jqhb/AS50QjyGg3l5CVcmRrPhYBF5x0+6Ohy7LFmfSe8AX+YlnvmqIKSHH5PiQ/hsVx7OGpWgqvz+o13kl1Xx9NVjWX3vdOaPj242MTV13aRYLh3TjydX7XPaF2P+8So2HiziysQol/dkHBLRk+V3nk1Qd1+ue3kz3x4scvgx6huUVd/nM2NIGN19ba/GvuGs/oQHdeOJVfuc9ndi2MfWKr4rm9zmi8hjmHnwPMa8xChU4cOtnjuzRE5xJat253PtxFib2tLmjIok4+gJDhRWOCWejQePsTmjmHvOT+DycVE2jwcSER69chRxoT24+/2tFJY7/r/Rh1sP06BwZSuJvKPE9glg+Z1nER0cwM2vJTu8ujklq5iiippT3eNt1d3Xm7vOSyA1u4Rv9nWNadA8ja1XUJc0uV0IlPPjiV3brL3LbVhfe8i6fZ+ImKu5M+jfpwcT4oJZkZrrsb8UX/82Cy8Rbjwrzqb9LxwRjgisdMISHKrK37/YS79e3bl2Uttn/Ars5sPz1yVSXlXLPe9tc+jUTKrKirRcxvcPJj60h8M+115hQd1ZesdkRkQF8T9vp/LvlJzW32Sjlbvy6ebjxYyhYW1+74KkGGJDAvj7F/uc3i5otJ2tbVC3NLndrqqPqGqhPQe2Z7kN634/wTLV0mzgeevnGS2YlxhN+tETbM897upQ2qy8qpalyTlcNDqSiF62dVkO69mdpP7BTpk8dvXuArbnHueeCxJOdYJoq6ERQfzpspFsyjjGP9fsd1hsO3KPc7CwotVqUFfoHeDH2wsnMWVQKPct38GSDZl2f6alk0o+0wf3JbBb2wcH+/l4ce/MBHbnlbHSyW2WRtu1ttzGsyLyTEs3O4/d7uU2rNvfV9VqVc0EDlo/z2jB3NGRdPPxYkWq53WWWJaSS0V1HQun/nhg7pnMGRnJ3vxyMovsHlN+Sn2D8uSq/cSH9rA7CSxIimH++Gie/fog6/Y7poppRVou3Xy8uGh0pEM+z9F6dPPhlZuSmDMygj9/spsXvrGvs8jWnFLyy6qYO6pt1XtNXTomioSwQJ5cvc+p64ntLyg3E/G2UWtXUClAKtAdSAQOWG9jAXsHctiz3IYt7wW65nIbzQnq7suFIyL4ePsRj/pPUt+gvP5tJhPighkd3bbxKo1tEo68ivpkxxH2FZRz78zBDpmH7s+XjWRwWE9+sXQb+cfta4+qrqvn4+1HmDUigl7+HT+Dg626+Xjzr2sTuXBEOE+v2U/xiZp2f9bnu/Lw9RbOGxre7s/w9hJ+NWsIGUdP8IGT2mkLyqq4ZvFmHlyx0ymf31m1th7UG6r6BjAaOFdVn1XVZ4HzsSQpt9cVl9toybzx0Rw/WctXe+yqne1Qq3cXkFN8klubmdaoNf16+zMmprfD2qFq6xv4x+r9DI3oycWjHHOF4u/nzXPXJVJVW89d76VRa8cv+K/3FlJaWcu8xGZ/q7kVby/h17OGUF3XwDubs9v1GarKyl35TBkUandCvnBEOKOje/HPNQcc/gOurr6Bu97dSmVNPT89d6BDP7uzs/UnYDCWiVkbBVq32aMty21w2nIbtrzXOM3UQaGEB3XzqDFRr27IJDrYn1kj2leFM3dkBDsPHyenuNLuWJan5pJ9rJJfzxri0O7bg8IC+euVo0jOKuGJL9o/cHR56mHCenZj2hkGMbuThPCeTB/clzc3Z1NT1/bEvOtwGbklJ5k70v4fCyKWhHm49CTvb3FcBw6AJ1fvZ0tWMY9eOZKEcDP3X1vYmqAeA7aKyOsi8gaQBjxq57FPLbchIn5YOj18fNo+jcttwP9fbuNj4CfWXn7xQAKwxc54Oj1vL+HycVF8s+8oRRXVrg6nVV/tLWBLVjE3nx3X7jnk5li/vOzt2lxVW88zXx5gXGxvzh/W9t5irblsbBTXTorlpXUZbM5o+/ioYxXVfLOvkCvGRdk9315HunVKHEfLq/mkHas/r9yVh7eXMHN4+6v3mpqWEMqk+BCe/eoglTV1DvnMr/cW8sI36VwzMYYrxrlfxxV3Z2svvteAScCHwAfAWdaqv3aztik1LrexB1jWuNyGdYFEsCy30ce63MYvgQet7/0eWAbsBj4HfqaqntOw4kLzE6Opa1A+2ubey8F/uiOPO95KZVhkED+Z2P7Fm2P7BDA8MojPdtrXDvXOd4fIO17FfRcOcdr8fr+7aBixIQHcv3xHm78gP9p2hLoGZd54z/oSPGdwXwaFBbJkQ2abhkA0Vu9NHhBCsIPmNRQR7rtwCEUV1bzxbfuqHZs6XHqSe5dtY1hkEH+4xJ61Xbuu1nrxDbXeJ2JZTTfHeutn3WYXVf1MVQer6kBVfcS67WHrWlCoapWqXqWqg1R1oqpmNHnvI9b3DVHVlfbG0lUkhPdkdHQvt+7N996WQ9z1XhpjY3rz/qLJ7eo+3NTcURGkHSptdyeEiuo6nv/6IFMG9eHsgaF2xXImAX4+/G3+aA4VV/L3Nlb1rUjLZVRULwZ7WBWSiHDrlHi+P1LGlsxim9+3r8DSO3OOA6r3mkqKC2HGkL68uDadsqradn9OTV0DP383jbp65fnrEts0w4Xxg9auoH5pvX+ymdsTTozLcKJ5idHszitjT577zff74tp0Hvpgp6Vt4tZJDumNNtvOar7XNmRy7EQNv541xO5YWjN5QB9uOqs/r3+bZfMX9t78Mr4/UuYRnSOac8W4KHoH+LZpXNTKnfmIwKwRjqnea+pXs4Zw/GQtr6zLaH3nFjz++V62Hirl8Xmj3WrAtKdprRffIuv9jGZu53VMiIajXTqmH77e4lZXUarKYyv38tjKvVwyph+Lb0hy2PIgg8ICSQgLbFc1X2llDYvXZzBzeHiHrS91/+yhRAf7c//y7Zysab3mekVqLr7ewqVjPTNB+ft5c92kWFbvKSD7mG1j1lbuymNCXIhT1poaGdWLi0ZHsmRDJsfa0Vb7+a58lmzI5Oaz49x2PJqnsHUuvqtEpKf18e9E5AMRGefc0AxnCe7hx3lDw/jPtiNOHZhoq/oG5Tcf7uLFtelcNymWp68ei5+P/WOMmpozKpLkrGKOlrftC+eldRlUVNfxq1mDHRrPmfTo5sPjV44m61hlq8tB1NU38OHWI8wYEub0Naac6caz4vDxEl7/NqvVfdOPVrC/oII5bZx7ry3uvWAwJ2vruef9bW2qaTh0rJL7lm9nTHQvHpo71GnxdRW2fgv8XlXLRWQqcAGWzgsvOi8sw9nmJUZTVFHNugOuHbxcU9fA3e9t5b0th/jZjIH85fKRTumFNmdkBA1qWSfJVoXlVby+MYtLx/RjaERQ629woLMHhXL95Fhe3ZhJanbLVX3rDxRRVFHtcZ0jThce1J2LR/djWXJOq20/jcuotHVy2LYYFBbIwxcPZ1tOKXP+uZ473kph1+EzTxNWVVvPT99NRYB/XZvY7mmwjB/YmqAa6xkuAhar6qfYP5OE4ULnWn9xr0h13fCxypo6bnszhU935vGbuUO570LnrYA7NKIncX0C2rRG1PNfp1NT38C9F3Tc1VNTD84ZRr9e/tz37x1U1TZf1bc8LZfgAF9mDHF81/eOduuUeE7U1LMs+czjkD7bmce42N5E9vJ3ajw3T4lnwwMzuOf8BL5NP8bFz27gtjeS2d7CulaPfLqHXYfLeHLBWKevINxV2JqgDovIS8DVwGci0q0N7/0REQkRkdUicsB632zlvojcZN3ngIjc1GT7IyKSIyLOWUuhC/Dz8eLSMf1YvbvArqlm2ut4ZS03LNnChgNHeXzeKBZNd+4IexFhzqhIvk0/xj3vb+WtzdnszS9rcQbr3JJK3v3uEAuSoolzUSN3YDcfHp83moyiE/xj9Y8nlD1eWcvq3QVcNjbK4VWirjAquhcT40J4bWNWi1XPh45ZVkp2ZvVeU70D/Lh35mA2Pngev5o5mOSsEi57biM3v7aF1OySU/t9vP0Ib23OZtH0AQ4bl2WArf13F2CZNfwJVS0VkUjgPjuO+yDwpao+Zl1m40HggaY7iEgI8AcgCVAgVUQ+VtUS4L/Av7DMC2i004KkGN7clMWMJ77hhsn9ufHs/k5pdAbLrNOHS09y8GgF6YUV/Dsll8yiEzx3bSJzHDRtUGsWTo3ncMlJNqUfOzUOrGd3H5L6B5MUF8KEuBBGR/eiu683z3x5AATuOi+hQ2JrydSEUK6ZGMsr6zOYPTKCxCYdNT7ZeYSauga3nLm8vW6dGs+db6eyendBs38Xn39v6eji6O7lrQnq7std5ydw85Q43tqczSvrM5n3wrdMSwhl/vhofvPBTsb3D+a+C53f07MrEVsHx1nbnxJU9TUR6QsEWmcSb/tBRfZhmdsvz5rsvlHVIaftc411nzusz1+y7vdek30qVDXQ1uMmJSVpSkpKe0LutLYeKuGltRl8sTsfXy8vrhgXxW3T4ts9JUt1XT1ZRZUcLKzgYGEF6Uct9xlFFVTV/vCrOKxnN55cMMYl0/KoKrklJ0nOKiY5q4SUrOJTCxv6egsjo3qxPaeUm8+O5+FLTl8BpuOVV9Vy4VPr8Pfz5tO7p50aU3Pl8xupqK7ji19Md1rVaEerb1DOfeJrwnt2Z/n/nP2j1y9/biN1DQ18ctc0F0T3gxPVdbzzXTaL12VQVFFDcIAo01UiAAAJxUlEQVQvn90zzenVjp2ViKSqatLp2226ghKRxiuZIcBrgC/wNjClnfGEq2pjn998oLlrYptnLDfab1xsMC/eMJ7MohO8uiGTf6fmsDQlh/OGhnH7tAFMHhDS4pefqpJTfJK0QyWnbnvyyv/fAnxRvf0ZFBbIWQP7MLBvIIPCLDdX9jgTEWJCAogJCTi16mzJiRpSs0tIzi4mJauEuNAe/HSGe0zs2bO7L3+dN5qbXt3C02sO8OCcoWQcrSDtUCkPzXFeu50reHsJN58dz58/2c32nFLGxPwwg/2R0pNsyyl1i6uUHt18WDR9IDdMjuPDrYcZ0S/IJCcnsLWK7wpgHJY5+FDVI43dzlsiImuA5iqKf9v0iaqqiDhtKUsRWQQsAoiNbf+UOZ1dfGgP/nz5SO6dOZi3N2fzxrdZXPPyZkZF9eL26QOYOzKCmvoGduQetySj7FK25ZRQVGFpvwrw82ZMdG/umD6AIRE9Gdg3kIF9Ax02lsnZgnv4ccHwcC5w0/aDcwb35eqkGBavS2f2yAjW7C7ASyyDXDubBUnRPLV6P69uzOSfP/lhNEtjB5eOan+yhb+fd7tWVTZsY2uCqmmaSESk1VZjVb2gpddEpEBEIptU8TW3/sNh4Nwmz6OBb2yMt2kci4HFYKnia+v7u5qQHn7cfX4Ci6YP4MOth3l5fQZ3v7eVPwT4UlZVd+rqKD60B9MH9yUxNpjE2GAGhwc6ZH0ko2W/vXgY6w4c5b5/b6eypp5pCX0JC3JOm6Er9ezue6p99KE5w06tovz5rnyGhPdkQF+ba/UND2drglpmbQPqLSK3A7cCr9hx3MZZyh+z3n/UzD5fAI826eE3C3jIjmMabdDd15trJsZydVIMX+0t5L87jhATHEBi/96MjQn26EGhniqouy+PXjmKW15LBuCBOZ13IOgtU+J4/dtM3tiUxQOzh1JYXkVydjF3u7jTitGxbEpQqvqEiMwEyrC0Qz2sqqvtOO5jWJLeQiAbSy9BRCQJuFNVb1PVYhH5M5ZlOQD+pKrF1v3+BlwLBIhILvCKqv7RjniMFnh5iVtXfXU1M4aEcXVSDF/uLWBWJz4nMSEBzBoewbvfHeLu8xL44vsCVGFuB/X4NNyDzb34/t+bRLyAa1T1HceH5DymF5/RGTQ0KCdq6ujZ3X2XdXeE5KxirnpxE3+5fCQrd+WRV1rFl786p1N1CjEsWurF19pyG0Ei8pCI/EtEZonFz4EMrFc9hmF0LC8v6fTJCSCpfzCjonrx0rp0NmcUM3tkhElOXUxrrdpvYanS2wncBnwNXAVcrqqXOTk2wzC6MBFh4dR4copPUt+gpnqvC2qtDWqAqo4CEJFXgDwgVlXbt/KbYRhGG8wdFclfV+7B19uLEf06dsJew/VaS1CnphVW1XoRyTXJyTCMjuLn48XiG5JQMNV7XVBrCWqMiDQuhiKAv/W5YBlja37SGIbhVE1nkzC6lnb14vNUIlIOnHkFOM8QChS5Ogg7dYYygCmHuzHlcC+2lqO/qv5oYk5bB+p2Fvua68roaUQkxdPL0RnKAKYc7saUw73YWw4zN41hGIbhlkyCMgzDMNxSV0tQi10dgIN0hnJ0hjKAKYe7MeVwL3aVo0t1kjAMwzA8R1e7gjIMwzA8RJdIUCIyW0T2ichBEXnQ1fG0l4hkichOEdkmIh4z662IvCoihSKyq8m2EBFZLSIHrPfBZ/oMd9BCOf4oIoet52SbiMx1ZYy2EJEYEflaRHaLyPcico91u8eckzOUwaPOh4h0F5EtIrLdWo7/tW6PF5HvrN9ZS0XErde3OUM5XheRzCbnY2ybPrezV/GJiDewH5iJZdn4ZCwzse92aWDtICJZQJKqetT4CBGZDlQAb6rqSOu2vwHFqvqY9UdDsKo+4Mo4W9NCOf4IVKjqE66MrS2si4RGqmqadWXsVOBy4GY85JycoQwL8KDzIZbpMXqoaoWI+AIbgHuAXwIfqOr7IvIisF1VX3BlrGdyhnLcCXyiqsvb87ld4QpqInBQVTNUtQZ4HzAT3XYgVV0HFJ+2+TLgDevjN7B8ubi1FsrhcVQ1T1XTrI/LgT1AFB50Ts5QBo+iFhXWp77WmwLnAY1f6m59LuCM5bBLV0hQUUBOk+e5eOAfspUCq0QkVUQWuToYO4Wrap71cT7gyavv/VxEdlirAN22Wqw5IhIHjAO+w0PPyWllAA87HyLiLSLbgEJgNZAOlKpqnXUXj/jOOr0cqtp4Ph6xno+nRKRbWz6zKySozmSqqiYCc4CfWaucPJ5a6pk9ta75BWAgMBbLbP9PujYc24lIILAC+IWqljV9zVPOSTNl8Ljzoar1qjoWiMZS4zPUxSG1y+nlEJGRwENYyjMBCAHaVGXcFRLUYSCmyfNo6zaPo6qHrfeFwIdY/pg9VYG1HaGxPaHQxfG0i6oWWP9jNgAv4yHnxNpOsAJ4R1U/sG72qHPSXBk89XwAqGopljX3zgJ6i0jjVHQe9Z3VpByzrVWxqqrVwGu08Xx0hQSVDCRYe8X4AT8BPnZxTG0mIj2sjcGISA9gFrDrzO9yax8DN1kf3wR85MJY2q3xC93qCjzgnFgbtJcAe1T1H01e8phz0lIZPO18iEhfEeltfeyPpTPXHixf8POtu7n1uYAWy7G3yQ8ewdKO1qbz0el78QFYu5o+DXgDr6rqIy4Oqc1EZACWqyawTPL7rqeUQ0TeA87FMrNxAfAH4D/AMiAWyAYWqKpbd0BooRznYqlOUiALuKNJO45bEpGpwHosK2U3WDf/BksbjkeckzOU4Ro86HyIyGgsnSC8sVwwLFPVP1n/v7+PpVpsK3C99SrELZ2hHF8BfbEs0bQNuLNJZ4rWP7crJCjDMAzD83SFKj7DMAzDA5kEZRiGYbglk6AMwzAMt2QSlGEYhuGWTIIyDMMw3JJP67sYhuFsItIH+NL6NAKoB45an1eq6tkuCcwwXMh0MzcMN+OJM6QbhjOYKj7DcHMiUmG9P1dE1orIRyKSISKPich11nV4dorIQOt+fUVkhYgkW29TXFsCw2gfk6AMw7OMwbLGzjDgBmCwqk4EXgHusu7zT+ApVZ0AzLO+Zhgex7RBGYZnSW6cukdE0oFV1u07gRnWxxcAwy3TnwEQJCKBbZlixjDcgUlQhuFZms7H1tDkeQM//H/2AiaralVHBmYYjmaq+Ayj81nFD9V9iMhYF8ZiGO1mEpRhdD53A0nWVUx3Y2mzMgyPY7qZG4ZhGG7JXEEZhmEYbskkKMMwDMMtmQRlGIZhuCWToAzDMAy3ZBKUYRiG4ZZMgjIMwzDckklQhmEYhlsyCcowDMNwS/8HTSgme0jFtmMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cpd = list(parafac(Z_log_norm.astype(float), 1))\n",
    "for i in iter(range(len(cpd))):\n",
    "    cpd[i] = cpd[i].ravel().copy(order = 'F')\n",
    "try:\n",
    "    del i\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "seasonal_decompose(cpd[-1], freq = 12).plot()\n",
    "adfuller(cpd[-1])\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "del cpd\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mse (a, b):\n",
    "    e = np.abs(b - a)\n",
    "\n",
    "    if isinstance(e, np.ndarray):\n",
    "        e = np.sqrt(np.square(e).mean())\n",
    "\n",
    "    return e\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rel_mse (a, b):\n",
    "    e = np.abs((b - a) / np.maximum(np.abs(a), 1))\n",
    "\n",
    "    if isinstance(e, np.ndarray):\n",
    "        e = np.sqrt(np.square(e).mean())\n",
    "\n",
    "    return e\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optim (\n",
    "    space,\n",
    "    Z,\n",
    "    loss,\n",
    "    k = None,\n",
    "    T0 = None,\n",
    "    p = None,\n",
    "    n = None,\n",
    "    theta = None\n",
    "):\n",
    "    \"\"\"\n",
    "    space: rječnik hiperparametara,\n",
    "    Z: kompletni tenzor (s rezovima koje treba predvidjeti),\n",
    "    loss: funkcija za računanje greške,\n",
    "    k: broj komponenti u CP dekompoziciji,\n",
    "    T0: broj rezova koji se uzimaju u obzir kod predikcije,\n",
    "    p: duljina perioda,\n",
    "    n: broj zadnjih rezova u Z koje treba predvidjeti,\n",
    "    theta: oblik (oblika) paramet(a)ra theta u eksponencijalnom izglađivanju.\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    if 'k' in space:\n",
    "        assert k is None\n",
    "    if 'T0' in space:\n",
    "        assert T0 is None\n",
    "    assert (p is None) == ('p' in space)\n",
    "    assert (n is None) == ('n' in space)\n",
    "    assert (theta is None) == ('theta' in space)\n",
    "\n",
    "    if n is None:\n",
    "        n = int(round(space['n']))\n",
    "\n",
    "    if theta is None:\n",
    "        theta = space['theta']\n",
    "    if isinstance(theta, numbers.Complex):\n",
    "        if isinstance(theta, numbers.Integral):\n",
    "            theta = np.zeros(theta, dtype = float, order = 'F')\n",
    "            for i in iter(range(int(theta.size))):\n",
    "                theta[i] = space['theta[{0:d}]'.format(i)]\n",
    "        else:\n",
    "            theta = space['theta']\n",
    "    else:\n",
    "        if not isinstance(theta, tuple):\n",
    "            if isinstance(theta, np.ndarray):\n",
    "                theta = copy.deepcopy((theta, theta, theta))\n",
    "            else:\n",
    "                if hasattr(theta, '__iter__'):\n",
    "                    theta = copy.deepcopy(tuple(theta))\n",
    "                else:\n",
    "                    theta = copy.deepcopy((theta, theta, theta))\n",
    "        else:\n",
    "            theta = copy.deepcopy(theta)\n",
    "        if len(theta) != 3:\n",
    "            theta = copy.deepcopy((theta, theta, theta))\n",
    "        theta = list(theta)\n",
    "        for i in iter(range(3)):\n",
    "            try:\n",
    "                theta[i] = np.zeros(theta[i], dtype = float, order = 'C')\n",
    "                idx = np.array(\n",
    "                    np.meshgrid(\n",
    "                        *tuple(\n",
    "                            np.arange(theta[i].shape[j], dtype = int)\n",
    "                                for j in iter(range(int(theta[i].ndim)))\n",
    "                        )\n",
    "                    )\n",
    "                ).T.reshape((-1, theta[i].ndim)).copy(order = 'C')\n",
    "                for j in iter(range(int(idx.shape[0]))):\n",
    "                    theta[i][tuple(idx[j])] = (\n",
    "                        space[\n",
    "                            'theta[{0:s}]'.format(\n",
    "                                ','.join(\n",
    "                                    str(r)\n",
    "                                        for r in np.concatenate(\n",
    "                                            ([i], idx[j])\n",
    "                                        )\n",
    "                                )\n",
    "                            )\n",
    "                        ]\n",
    "                    )\n",
    "                try:\n",
    "                    del j\n",
    "                except (NameError, UnboundLocalError):\n",
    "                    pass\n",
    "            except (TypeError, ValueError, AttributeError):\n",
    "                assert (\n",
    "                    hasattr(theta[i], '__iter__') or\n",
    "                    hasattr(theta[i], '__array__')\n",
    "                )\n",
    "                theta[i] = np.array(theta[i], copy = True)\n",
    "        try:\n",
    "            del i\n",
    "        except (NameError, UnboundLocalError):\n",
    "            pass\n",
    "        theta = tuple(theta)\n",
    "\n",
    "    Z = np.moveaxis(Z, -1, 0).copy(order = 'F')\n",
    "\n",
    "    X = np.moveaxis(Z[:-n], 0, -1).copy(order = 'F')\n",
    "    Y = np.moveaxis(Z[-n:], 0, -1).copy(order = 'F')\n",
    "    if n == 1:\n",
    "        if Y.ndim == 1:\n",
    "            Y = Y.dtype.type(Y[0])\n",
    "        else:\n",
    "            Y.shape = Y.shape[:-1]\n",
    "\n",
    "    del Z\n",
    "\n",
    "    Y_predicted = None\n",
    "\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.filterwarnings('error')\n",
    "        with np.errstate(divide = 'raise', invalid = 'raise'):\n",
    "            try:\n",
    "                if n == 1:\n",
    "                    Y_predicted = cp_score(\n",
    "                        X,\n",
    "                        k = int(round(space['k'])) if 'k' in space else k,\n",
    "                        T0 = int(round(space['T0'])) if 'T0' in space else T0,\n",
    "                        predict = (\n",
    "                            lambda X : ExponentialSmooth().fit(\n",
    "                                X.T,\n",
    "                                int(round(space['p'])) if 'p' in space else p\n",
    "                            ).predict(1, theta)\n",
    "                        )\n",
    "                    )\n",
    "                else:\n",
    "                    Y_predicted = cp_score(\n",
    "                        X,\n",
    "                        k = int(round(space['k'])) if 'k' in space else k,\n",
    "                        T0 = int(round(space['T0'])) if 'T0' in space else T0,\n",
    "                        predict = (\n",
    "                            lambda X : ExponentialSmooth().fit(\n",
    "                                X.T,\n",
    "                                int(round(space['p'])) if 'p' in space else p\n",
    "                            ).predict(n, theta).T\n",
    "                        )\n",
    "                    )\n",
    "            except (\n",
    "                TypeError,\n",
    "                ValueError,\n",
    "                FloatingPointError,\n",
    "                np.linalg.LinAlgError,\n",
    "                RuntimeWarning\n",
    "            ):\n",
    "                return {'status' : ho.STATUS_FAIL}\n",
    "\n",
    "    return {'loss' : loss(Y, Y_predicted), 'status' : ho.STATUS_OK}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_theta (space, k = None):\n",
    "    space = copy.deepcopy(space)\n",
    "\n",
    "    if 'k' in space:\n",
    "        assert k is None\n",
    "        space.update(\n",
    "            {\n",
    "                'theta' : (\n",
    "                    int(round(space['k'])),\n",
    "                    int(round(space['k'])),\n",
    "                    int(round(space['k']))\n",
    "                )\n",
    "            }\n",
    "        )\n",
    "    elif k is None:\n",
    "        space.update(\n",
    "            {\n",
    "                'theta' : (\n",
    "                    int(round(default_space['k'])),\n",
    "                    int(round(default_space['k'])),\n",
    "                    int(round(default_space['k']))\n",
    "                )\n",
    "            }\n",
    "        )\n",
    "    else:\n",
    "        space.update({'theta' : (k, k, k)})\n",
    "\n",
    "    return space\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_space = {\n",
    "    'k' : 24,\n",
    "    'T0' : 24,\n",
    "    'p' : 12,\n",
    "    'n' : 1\n",
    "}\n",
    "#default.space.update({'theta' : 1})\n",
    "default_space.update(\n",
    "    {'theta' : (default_space['k'], default_space['k'], default_space['k'])}\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "space = dict()\n",
    "#space.update(\n",
    "#   {\n",
    "#       'k' : ho.hp.quniform('k', 1.0, float(min(Z.shape)), 1.0),\n",
    "#       'T0' : ho.hp.quniform('T0', 1.0, float(Z.shape[-1]), 1.0),\n",
    "#       'p' : ho.hp.quniform('p', 1.0, float(Z.shape[-1]), 1.0),\n",
    "#       'n' : ho.hp.quniform('n', 1.0, float(Z.shape[-1]) - 1.0, 1.0)\n",
    "#   }\n",
    "#)\n",
    "#space.update({'k' : ho.hp.quniform('k', 1.0, float(default_space['k']), 1.0)})\n",
    "space.update(\n",
    "    dict(\n",
    "            (\n",
    "                'theta[{I:d},{J:d}]'.format(I = i, J = j),\n",
    "                ho.hp.uniform(\n",
    "                    'theta[{I:d},{J:d}]'.format(I = i, J = j),\n",
    "                    0.0,\n",
    "                    1.0\n",
    "                )\n",
    "            )\n",
    "                for i in iter(range(3))\n",
    "                for j in iter(range(default_space['k']))\n",
    "    )\n",
    ")\n",
    "try:\n",
    "    del i\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "try:\n",
    "    del j\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100%|██████████| 2000/2000 [13:56<00:00,  1.80it/s, best loss: 0.19946346136660378]\n",
      "1: 8\n",
      "100%|██████████| 2000/2000 [12:42<00:00,  1.78it/s, best loss: 0.15349589493829424]\n",
      "2: 13\n",
      "100%|██████████| 2000/2000 [12:30<00:00,  1.83it/s, best loss: 0.2129370908178003]\n",
      "3: 7\n",
      "100%|██████████| 2000/2000 [12:23<00:00,  1.81it/s, best loss: 0.14365699041043944]\n",
      "4: 7\n",
      "100%|██████████| 2000/2000 [12:26<00:00,  1.78it/s, best loss: 0.16219996987664098]\n",
      "5: 6\n",
      "100%|██████████| 2000/2000 [13:02<00:00,  1.73it/s, best loss: 0.12353518433474606]\n",
      "6: 19\n",
      "100%|██████████| 2000/2000 [15:11<00:00,  1.67it/s, best loss: 0.12031037952551654]\n",
      "7: 15\n",
      "100%|██████████| 2000/2000 [14:11<00:00,  1.26it/s, best loss: 0.12828985645020397]\n",
      "8: 17\n",
      "100%|██████████| 2000/2000 [14:36<00:00,  1.69it/s, best loss: 0.1217945024841068]\n",
      "9: 19\n",
      "100%|██████████| 2000/2000 [13:48<00:00,  1.64it/s, best loss: 0.1682958054214247]\n",
      "10: 12\n",
      "100%|██████████| 2000/2000 [13:52<00:00,  1.67it/s, best loss: 0.11591631573362954]\n",
      "11: 17\n",
      "100%|██████████| 2000/2000 [13:42<00:00,  1.69it/s, best loss: 0.14243428019720158]\n",
      "12: 11\n"
     ]
    }
   ],
   "source": [
    "best = [dict() for r in iter(range(12))]\n",
    "try:\n",
    "    del r\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "for r in iter(range(len(best))):\n",
    "    W = Z[:, :, r:r + 25].copy(order = 'F')\n",
    "    W_y = W[:, :, -1].copy(order = 'F')\n",
    "\n",
    "    W_mean = W[:, :, :-1].mean(axis = -1).copy(order = 'F')\n",
    "    W_std = W[:, :, :-1].std(axis = -1, ddof = 1).copy(order = 'F')\n",
    "\n",
    "    W_norm = (W - W_mean.reshape((4, 21, 1))) / W_std.reshape((4, 21, 1))\n",
    "\n",
    "    W_log = np.log1p(W).copy(order = 'F')\n",
    "\n",
    "    W_log_mean = W_log[:, :, :-1].mean(axis = -1).copy(order = 'F')\n",
    "    W_log_std = W_log[:, :, :-1].std(axis = -1, ddof = 1).copy(order = 'F')\n",
    "\n",
    "    W_log_norm = (\n",
    "        (W_log - W_log_mean.reshape((4, 21, 1))) / W_log_std.reshape((4, 21, 1))\n",
    "    )\n",
    "\n",
    "    trials = ho.Trials()\n",
    "    best[r] = ho.fmin(\n",
    "        fn = lambda S : (\n",
    "            optim(\n",
    "                update_theta(S),\n",
    "                Z = W_log_norm,\n",
    "                loss = (\n",
    "                    lambda a, b : rel_mse(\n",
    "                        W_y,\n",
    "                        np.expm1(W_log_std * b + W_log_mean)\n",
    "                    )\n",
    "                ),\n",
    "                k = None if 'k' in space else default_space['k'],\n",
    "                T0 = default_space['T0'],\n",
    "                p = default_space['p'],\n",
    "                n = default_space['n']\n",
    "            )\n",
    "        ),\n",
    "        space = space,\n",
    "        algo = ho.tpe.suggest,\n",
    "        max_evals = 2000,\n",
    "        trials = trials,\n",
    "        show_progressbar = True\n",
    "    )\n",
    "\n",
    "    if 'k' in best[r]:\n",
    "        print('{R:d}: {K:d}'.format(R = r + 1, K = int(round(best[r]['k']))))\n",
    "\n",
    "        for i in iter(range(3)):\n",
    "            for j in iter(\n",
    "                range(int(round(best[r]['k'])), int(round(default_space['k'])))\n",
    "            ):\n",
    "                best[r].pop('theta[{I:d},{J:d}]'.format(I = i, J = j))\n",
    "            try:\n",
    "                del j\n",
    "            except (NameError, UnboundLocalError):\n",
    "                pass\n",
    "        try:\n",
    "            del i\n",
    "        except (NameError, UnboundLocalError):\n",
    "            pass\n",
    "\n",
    "    del W\n",
    "    del W_mean\n",
    "    del W_std\n",
    "    del W_norm\n",
    "    del W_log\n",
    "    del W_log_mean\n",
    "    del W_log_std\n",
    "    del W_log_norm\n",
    "try:\n",
    "    del r\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "best = tuple(best)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "({'k': 8.0,\n",
       "  'theta[0,0]': 0.4412651671162238,\n",
       "  'theta[0,1]': 0.25625054451861956,\n",
       "  'theta[0,2]': 0.696626857755847,\n",
       "  'theta[0,3]': 0.4310122636835551,\n",
       "  'theta[0,4]': 0.4670730931643621,\n",
       "  'theta[0,5]': 0.15332721757561557,\n",
       "  'theta[0,6]': 0.867464233481899,\n",
       "  'theta[0,7]': 0.8920571056289867,\n",
       "  'theta[1,0]': 0.013181023961012617,\n",
       "  'theta[1,1]': 0.755887722631181,\n",
       "  'theta[1,2]': 0.2111965943875706,\n",
       "  'theta[1,3]': 0.46812910567887456,\n",
       "  'theta[1,4]': 0.09106691743082117,\n",
       "  'theta[1,5]': 0.8817479983072335,\n",
       "  'theta[1,6]': 0.4409550118424658,\n",
       "  'theta[1,7]': 0.5229624349218568,\n",
       "  'theta[2,0]': 0.11047840912397291,\n",
       "  'theta[2,1]': 0.7521683490694231,\n",
       "  'theta[2,2]': 0.929273552172044,\n",
       "  'theta[2,3]': 0.6952854314038778,\n",
       "  'theta[2,4]': 0.19085323912118649,\n",
       "  'theta[2,5]': 0.6545089109797336,\n",
       "  'theta[2,6]': 0.5026060212130528,\n",
       "  'theta[2,7]': 0.07851049904796004},\n",
       " {'k': 13.0,\n",
       "  'theta[0,0]': 0.22180515153149571,\n",
       "  'theta[0,10]': 0.024513313978159663,\n",
       "  'theta[0,11]': 0.5745128939422591,\n",
       "  'theta[0,12]': 0.29550212307984136,\n",
       "  'theta[0,1]': 0.8586119520190444,\n",
       "  'theta[0,2]': 0.4458751003057937,\n",
       "  'theta[0,3]': 0.6330847929369652,\n",
       "  'theta[0,4]': 0.28137359698293074,\n",
       "  'theta[0,5]': 0.8885329935572605,\n",
       "  'theta[0,6]': 0.6907188952050565,\n",
       "  'theta[0,7]': 0.22954131892786633,\n",
       "  'theta[0,8]': 0.09417126398288551,\n",
       "  'theta[0,9]': 0.0666879613293355,\n",
       "  'theta[1,0]': 0.744113727051422,\n",
       "  'theta[1,10]': 0.03709067522359689,\n",
       "  'theta[1,11]': 0.14913574637665464,\n",
       "  'theta[1,12]': 0.1701219366016377,\n",
       "  'theta[1,1]': 0.9300928712344235,\n",
       "  'theta[1,2]': 0.5597955679033022,\n",
       "  'theta[1,3]': 0.9580853268787918,\n",
       "  'theta[1,4]': 0.8006711729241831,\n",
       "  'theta[1,5]': 0.4162164710164478,\n",
       "  'theta[1,6]': 0.09701148269987264,\n",
       "  'theta[1,7]': 0.2674974245495425,\n",
       "  'theta[1,8]': 0.1025554604139313,\n",
       "  'theta[1,9]': 0.8360414392116204,\n",
       "  'theta[2,0]': 0.9543288742094421,\n",
       "  'theta[2,10]': 0.6626691284175656,\n",
       "  'theta[2,11]': 0.6941467554041283,\n",
       "  'theta[2,12]': 0.6351617802862595,\n",
       "  'theta[2,1]': 0.0006679384384742082,\n",
       "  'theta[2,2]': 0.9418886961215871,\n",
       "  'theta[2,3]': 0.95979446427616,\n",
       "  'theta[2,4]': 0.9541515218728465,\n",
       "  'theta[2,5]': 0.5079615297431453,\n",
       "  'theta[2,6]': 0.38253448846447147,\n",
       "  'theta[2,7]': 0.706393615349558,\n",
       "  'theta[2,8]': 0.7057651179865015,\n",
       "  'theta[2,9]': 0.2452200918455297},\n",
       " {'k': 7.0,\n",
       "  'theta[0,0]': 0.22348727304085855,\n",
       "  'theta[0,1]': 0.9997320712113306,\n",
       "  'theta[0,2]': 0.8963209393011883,\n",
       "  'theta[0,3]': 0.7824826196212821,\n",
       "  'theta[0,4]': 0.14313990868089738,\n",
       "  'theta[0,5]': 0.9234216043521764,\n",
       "  'theta[0,6]': 0.11870806450326976,\n",
       "  'theta[1,0]': 0.8467130438656199,\n",
       "  'theta[1,1]': 0.9761883240804103,\n",
       "  'theta[1,2]': 0.9286884915296487,\n",
       "  'theta[1,3]': 0.0006452248886488765,\n",
       "  'theta[1,4]': 0.43299616214876585,\n",
       "  'theta[1,5]': 0.6425591421047931,\n",
       "  'theta[1,6]': 0.6918019022526608,\n",
       "  'theta[2,0]': 0.999965241572744,\n",
       "  'theta[2,1]': 0.9340414824772739,\n",
       "  'theta[2,2]': 0.2713278954627223,\n",
       "  'theta[2,3]': 0.7253414235907887,\n",
       "  'theta[2,4]': 0.9999880049600526,\n",
       "  'theta[2,5]': 0.9343232408794568,\n",
       "  'theta[2,6]': 0.8164749024288755},\n",
       " {'k': 7.0,\n",
       "  'theta[0,0]': 0.34989842809114624,\n",
       "  'theta[0,1]': 0.16712754591864032,\n",
       "  'theta[0,2]': 0.3861455274164493,\n",
       "  'theta[0,3]': 0.2266687502641848,\n",
       "  'theta[0,4]': 0.9447193296876999,\n",
       "  'theta[0,5]': 0.4734340465751576,\n",
       "  'theta[0,6]': 0.3496537196825508,\n",
       "  'theta[1,0]': 0.7730883212536702,\n",
       "  'theta[1,1]': 0.42859689326507244,\n",
       "  'theta[1,2]': 0.6876004919127403,\n",
       "  'theta[1,3]': 0.2870251423605282,\n",
       "  'theta[1,4]': 0.7977992562989599,\n",
       "  'theta[1,5]': 0.5812752473018687,\n",
       "  'theta[1,6]': 0.35958114410367675,\n",
       "  'theta[2,0]': 0.978833369718097,\n",
       "  'theta[2,1]': 0.32037834134977144,\n",
       "  'theta[2,2]': 0.3770033822575118,\n",
       "  'theta[2,3]': 0.7141887825010353,\n",
       "  'theta[2,4]': 0.13758233870373238,\n",
       "  'theta[2,5]': 0.8531881651018143,\n",
       "  'theta[2,6]': 0.29147750948275286},\n",
       " {'k': 6.0,\n",
       "  'theta[0,0]': 0.7732847794317451,\n",
       "  'theta[0,1]': 0.4980988042223483,\n",
       "  'theta[0,2]': 0.21425724648695132,\n",
       "  'theta[0,3]': 0.872278086107851,\n",
       "  'theta[0,4]': 0.8085126105843706,\n",
       "  'theta[0,5]': 0.7114110944830613,\n",
       "  'theta[1,0]': 0.8444618262074036,\n",
       "  'theta[1,1]': 0.05856682607433944,\n",
       "  'theta[1,2]': 0.14408663435544,\n",
       "  'theta[1,3]': 0.34080636834544364,\n",
       "  'theta[1,4]': 0.3133727049985727,\n",
       "  'theta[1,5]': 0.9716141510244888,\n",
       "  'theta[2,0]': 0.8717263402530607,\n",
       "  'theta[2,1]': 0.27930157087528695,\n",
       "  'theta[2,2]': 0.8698171166510085,\n",
       "  'theta[2,3]': 0.902821826222007,\n",
       "  'theta[2,4]': 0.798968588682663,\n",
       "  'theta[2,5]': 0.18363784278854456},\n",
       " {'k': 19.0,\n",
       "  'theta[0,0]': 0.48118619839282245,\n",
       "  'theta[0,10]': 0.19603964526143497,\n",
       "  'theta[0,11]': 0.3467037138502268,\n",
       "  'theta[0,12]': 0.8021347177674241,\n",
       "  'theta[0,13]': 0.2517834518266908,\n",
       "  'theta[0,14]': 0.8327124169093987,\n",
       "  'theta[0,15]': 0.7562177601786806,\n",
       "  'theta[0,16]': 0.4812137545130527,\n",
       "  'theta[0,17]': 0.45424387759594714,\n",
       "  'theta[0,18]': 0.05954621468876023,\n",
       "  'theta[0,1]': 0.5417501670992305,\n",
       "  'theta[0,2]': 0.2987254497314378,\n",
       "  'theta[0,3]': 0.8998165584568242,\n",
       "  'theta[0,4]': 0.6217688314180622,\n",
       "  'theta[0,5]': 0.19375368828408224,\n",
       "  'theta[0,6]': 0.7984037460662272,\n",
       "  'theta[0,7]': 0.24017751820690614,\n",
       "  'theta[0,8]': 0.7097614076714265,\n",
       "  'theta[0,9]': 0.41984069744913716,\n",
       "  'theta[1,0]': 0.0618064114879374,\n",
       "  'theta[1,10]': 0.11649301880537984,\n",
       "  'theta[1,11]': 0.37340973042529096,\n",
       "  'theta[1,12]': 0.49613222985832833,\n",
       "  'theta[1,13]': 0.5703172111834645,\n",
       "  'theta[1,14]': 0.6924138206080253,\n",
       "  'theta[1,15]': 0.29316660939288286,\n",
       "  'theta[1,16]': 0.9153747346499229,\n",
       "  'theta[1,17]': 0.5737530480760913,\n",
       "  'theta[1,18]': 0.5536768441420348,\n",
       "  'theta[1,1]': 0.834599978544904,\n",
       "  'theta[1,2]': 0.6791498836705979,\n",
       "  'theta[1,3]': 0.34944425237598875,\n",
       "  'theta[1,4]': 0.12606746200887448,\n",
       "  'theta[1,5]': 0.3721854878683345,\n",
       "  'theta[1,6]': 0.2840171865688801,\n",
       "  'theta[1,7]': 0.2689753945941694,\n",
       "  'theta[1,8]': 0.927476581899783,\n",
       "  'theta[1,9]': 0.6246732336763258,\n",
       "  'theta[2,0]': 0.1271080730766394,\n",
       "  'theta[2,10]': 0.880648185912847,\n",
       "  'theta[2,11]': 0.8691104103654737,\n",
       "  'theta[2,12]': 0.1259319752353343,\n",
       "  'theta[2,13]': 0.43195117365728,\n",
       "  'theta[2,14]': 0.36180527787709493,\n",
       "  'theta[2,15]': 0.9369085605129124,\n",
       "  'theta[2,16]': 0.5191460478572942,\n",
       "  'theta[2,17]': 0.5156700136533531,\n",
       "  'theta[2,18]': 0.25644643554903207,\n",
       "  'theta[2,1]': 0.6348376044688744,\n",
       "  'theta[2,2]': 0.08041917568290774,\n",
       "  'theta[2,3]': 0.8732453466464752,\n",
       "  'theta[2,4]': 0.06486735918388141,\n",
       "  'theta[2,5]': 0.3095635721990667,\n",
       "  'theta[2,6]': 0.054612363680857395,\n",
       "  'theta[2,7]': 0.79156854127722,\n",
       "  'theta[2,8]': 0.9983476525904292,\n",
       "  'theta[2,9]': 0.373894538416793},\n",
       " {'k': 15.0,\n",
       "  'theta[0,0]': 0.5454959025331912,\n",
       "  'theta[0,10]': 0.34983851039724995,\n",
       "  'theta[0,11]': 0.14162615817487928,\n",
       "  'theta[0,12]': 0.7833285216240913,\n",
       "  'theta[0,13]': 0.22691944123516603,\n",
       "  'theta[0,14]': 0.01637345940039637,\n",
       "  'theta[0,1]': 0.42722339173235374,\n",
       "  'theta[0,2]': 0.5747012264233978,\n",
       "  'theta[0,3]': 0.7593776585493467,\n",
       "  'theta[0,4]': 0.20621264089129515,\n",
       "  'theta[0,5]': 0.8772061591790346,\n",
       "  'theta[0,6]': 0.6580378961039908,\n",
       "  'theta[0,7]': 0.14316208067103708,\n",
       "  'theta[0,8]': 0.5372777856172113,\n",
       "  'theta[0,9]': 0.7363594935478099,\n",
       "  'theta[1,0]': 0.8286190158613365,\n",
       "  'theta[1,10]': 0.26988501226814643,\n",
       "  'theta[1,11]': 0.8771081066073603,\n",
       "  'theta[1,12]': 0.1129839983942626,\n",
       "  'theta[1,13]': 0.3809861924396369,\n",
       "  'theta[1,14]': 0.3032523851195197,\n",
       "  'theta[1,1]': 0.0832754449378107,\n",
       "  'theta[1,2]': 0.6970008179190029,\n",
       "  'theta[1,3]': 0.6997085959406911,\n",
       "  'theta[1,4]': 0.8679168745984424,\n",
       "  'theta[1,5]': 0.5610115983899712,\n",
       "  'theta[1,6]': 0.3063502493323466,\n",
       "  'theta[1,7]': 0.6597994615368208,\n",
       "  'theta[1,8]': 0.744017957065209,\n",
       "  'theta[1,9]': 0.885569310595156,\n",
       "  'theta[2,0]': 0.6456278006894965,\n",
       "  'theta[2,10]': 0.14522069122889933,\n",
       "  'theta[2,11]': 0.7594682301240528,\n",
       "  'theta[2,12]': 0.36226971880415537,\n",
       "  'theta[2,13]': 0.3631337795757019,\n",
       "  'theta[2,14]': 0.17578795829302524,\n",
       "  'theta[2,1]': 0.7007771403480932,\n",
       "  'theta[2,2]': 0.904270006844561,\n",
       "  'theta[2,3]': 0.6773478622841239,\n",
       "  'theta[2,4]': 0.27863037252531103,\n",
       "  'theta[2,5]': 0.11960236094200566,\n",
       "  'theta[2,6]': 0.9672579626839379,\n",
       "  'theta[2,7]': 0.43647878359510645,\n",
       "  'theta[2,8]': 0.2582572543427862,\n",
       "  'theta[2,9]': 0.12129475855449345},\n",
       " {'k': 17.0,\n",
       "  'theta[0,0]': 0.528933689093686,\n",
       "  'theta[0,10]': 0.4213566629221239,\n",
       "  'theta[0,11]': 0.8374590468797857,\n",
       "  'theta[0,12]': 0.3766508122131053,\n",
       "  'theta[0,13]': 0.9142691189811577,\n",
       "  'theta[0,14]': 0.6017130088967376,\n",
       "  'theta[0,15]': 0.6370874022721733,\n",
       "  'theta[0,16]': 0.40019881429471316,\n",
       "  'theta[0,1]': 0.5792081691810814,\n",
       "  'theta[0,2]': 0.4381399486720177,\n",
       "  'theta[0,3]': 0.10639654639792359,\n",
       "  'theta[0,4]': 0.9993517154460222,\n",
       "  'theta[0,5]': 0.2834862669038619,\n",
       "  'theta[0,6]': 0.5070334260817009,\n",
       "  'theta[0,7]': 0.30233074134394594,\n",
       "  'theta[0,8]': 0.7451743799253256,\n",
       "  'theta[0,9]': 0.33630627578764233,\n",
       "  'theta[1,0]': 0.235766966351333,\n",
       "  'theta[1,10]': 0.6450406439417268,\n",
       "  'theta[1,11]': 0.6995682057051938,\n",
       "  'theta[1,12]': 0.33505794346138895,\n",
       "  'theta[1,13]': 0.7302912760948068,\n",
       "  'theta[1,14]': 0.6936618859753111,\n",
       "  'theta[1,15]': 0.692099455301729,\n",
       "  'theta[1,16]': 0.8175998492432148,\n",
       "  'theta[1,1]': 0.5620620016577694,\n",
       "  'theta[1,2]': 0.9266538571139955,\n",
       "  'theta[1,3]': 0.6097207185163109,\n",
       "  'theta[1,4]': 0.25278996045197133,\n",
       "  'theta[1,5]': 0.9002517361282543,\n",
       "  'theta[1,6]': 0.26194460087929317,\n",
       "  'theta[1,7]': 0.6862012613275964,\n",
       "  'theta[1,8]': 0.09170077882152289,\n",
       "  'theta[1,9]': 0.12214380273863822,\n",
       "  'theta[2,0]': 0.10207148169001785,\n",
       "  'theta[2,10]': 0.7834195274660993,\n",
       "  'theta[2,11]': 0.905727115693834,\n",
       "  'theta[2,12]': 0.09366937991485079,\n",
       "  'theta[2,13]': 0.6182086594431908,\n",
       "  'theta[2,14]': 0.512921919304206,\n",
       "  'theta[2,15]': 0.6465719249780122,\n",
       "  'theta[2,16]': 0.29485061546896596,\n",
       "  'theta[2,1]': 0.059402011321588434,\n",
       "  'theta[2,2]': 0.30184595352847593,\n",
       "  'theta[2,3]': 0.3690503598831867,\n",
       "  'theta[2,4]': 0.4431147877965523,\n",
       "  'theta[2,5]': 0.4800991385586125,\n",
       "  'theta[2,6]': 0.13141868147777377,\n",
       "  'theta[2,7]': 0.7111277158869778,\n",
       "  'theta[2,8]': 0.4264501120719979,\n",
       "  'theta[2,9]': 0.9122295461032264},\n",
       " {'k': 19.0,\n",
       "  'theta[0,0]': 0.2885578434950399,\n",
       "  'theta[0,10]': 0.5440014028337139,\n",
       "  'theta[0,11]': 0.22432044191815848,\n",
       "  'theta[0,12]': 0.5520668465990193,\n",
       "  'theta[0,13]': 0.8698485397718212,\n",
       "  'theta[0,14]': 0.4280906653712543,\n",
       "  'theta[0,15]': 0.5508334179237764,\n",
       "  'theta[0,16]': 0.38777128512265596,\n",
       "  'theta[0,17]': 0.4780448062733614,\n",
       "  'theta[0,18]': 0.6416059760877831,\n",
       "  'theta[0,1]': 0.7296842133952055,\n",
       "  'theta[0,2]': 0.3372869505294454,\n",
       "  'theta[0,3]': 0.7260725320669253,\n",
       "  'theta[0,4]': 0.5752615244532281,\n",
       "  'theta[0,5]': 0.28893146384267154,\n",
       "  'theta[0,6]': 0.36961673391657446,\n",
       "  'theta[0,7]': 0.8757879835019874,\n",
       "  'theta[0,8]': 0.2577651789775955,\n",
       "  'theta[0,9]': 0.7010491756220956,\n",
       "  'theta[1,0]': 0.8060269451334479,\n",
       "  'theta[1,10]': 0.6132472178398753,\n",
       "  'theta[1,11]': 0.6315169460451718,\n",
       "  'theta[1,12]': 0.8619805867283833,\n",
       "  'theta[1,13]': 0.08648772762232043,\n",
       "  'theta[1,14]': 0.7611390481500678,\n",
       "  'theta[1,15]': 0.3270925885933659,\n",
       "  'theta[1,16]': 0.5488878581924461,\n",
       "  'theta[1,17]': 0.20471808021682486,\n",
       "  'theta[1,18]': 0.3888804108985178,\n",
       "  'theta[1,1]': 0.7472595681372395,\n",
       "  'theta[1,2]': 0.27939503928538273,\n",
       "  'theta[1,3]': 0.5490070765395593,\n",
       "  'theta[1,4]': 0.822325443675692,\n",
       "  'theta[1,5]': 0.9165424236836558,\n",
       "  'theta[1,6]': 0.6395819166733906,\n",
       "  'theta[1,7]': 0.514434856745126,\n",
       "  'theta[1,8]': 0.1092408045453581,\n",
       "  'theta[1,9]': 0.03704277611598575,\n",
       "  'theta[2,0]': 0.2724393088315462,\n",
       "  'theta[2,10]': 0.3938519331240013,\n",
       "  'theta[2,11]': 0.3865886961541708,\n",
       "  'theta[2,12]': 0.8009859654224211,\n",
       "  'theta[2,13]': 0.5799922423796734,\n",
       "  'theta[2,14]': 0.08161578162426052,\n",
       "  'theta[2,15]': 0.2791707629491813,\n",
       "  'theta[2,16]': 0.5379380786695781,\n",
       "  'theta[2,17]': 0.2032252285031952,\n",
       "  'theta[2,18]': 0.7073884874500638,\n",
       "  'theta[2,1]': 0.4190231754201637,\n",
       "  'theta[2,2]': 0.7984491659176735,\n",
       "  'theta[2,3]': 0.33920402771568275,\n",
       "  'theta[2,4]': 0.5237334919393145,\n",
       "  'theta[2,5]': 0.5683325167189308,\n",
       "  'theta[2,6]': 0.49421257012351494,\n",
       "  'theta[2,7]': 0.49936427047939247,\n",
       "  'theta[2,8]': 0.3220507222347033,\n",
       "  'theta[2,9]': 0.8806354022510281},\n",
       " {'k': 12.0,\n",
       "  'theta[0,0]': 0.5009094443258925,\n",
       "  'theta[0,10]': 0.4145365442927219,\n",
       "  'theta[0,11]': 0.9932163657080815,\n",
       "  'theta[0,1]': 0.40724446127477604,\n",
       "  'theta[0,2]': 0.04272675903939681,\n",
       "  'theta[0,3]': 0.3072169464895358,\n",
       "  'theta[0,4]': 0.5873736562648162,\n",
       "  'theta[0,5]': 0.9957173785046145,\n",
       "  'theta[0,6]': 0.0009503718785932713,\n",
       "  'theta[0,7]': 0.27144662093619054,\n",
       "  'theta[0,8]': 0.5564774264821489,\n",
       "  'theta[0,9]': 0.41756093344900747,\n",
       "  'theta[1,0]': 0.6776908782917106,\n",
       "  'theta[1,10]': 0.08096712741974768,\n",
       "  'theta[1,11]': 0.8072682912132297,\n",
       "  'theta[1,1]': 0.9902472246197508,\n",
       "  'theta[1,2]': 0.5172399472000444,\n",
       "  'theta[1,3]': 0.5877186539982489,\n",
       "  'theta[1,4]': 0.4620296931357893,\n",
       "  'theta[1,5]': 0.37092948317860164,\n",
       "  'theta[1,6]': 0.6049756369056202,\n",
       "  'theta[1,7]': 0.849363055894701,\n",
       "  'theta[1,8]': 0.5734601112904432,\n",
       "  'theta[1,9]': 0.030043271923692938,\n",
       "  'theta[2,0]': 0.35775228546924653,\n",
       "  'theta[2,10]': 0.02899671354170663,\n",
       "  'theta[2,11]': 0.3615306273594776,\n",
       "  'theta[2,1]': 0.8845694256620871,\n",
       "  'theta[2,2]': 0.4534364481247878,\n",
       "  'theta[2,3]': 0.9928266101734101,\n",
       "  'theta[2,4]': 0.23514097693497935,\n",
       "  'theta[2,5]': 0.9997636617826378,\n",
       "  'theta[2,6]': 0.6922560291400964,\n",
       "  'theta[2,7]': 0.6948716679265787,\n",
       "  'theta[2,8]': 0.3341949675645253,\n",
       "  'theta[2,9]': 0.6256831033894328},\n",
       " {'k': 17.0,\n",
       "  'theta[0,0]': 0.9489477465630057,\n",
       "  'theta[0,10]': 0.9267770781848268,\n",
       "  'theta[0,11]': 0.1786217143815108,\n",
       "  'theta[0,12]': 0.9616631332404935,\n",
       "  'theta[0,13]': 0.27084464054800755,\n",
       "  'theta[0,14]': 0.5359887231163851,\n",
       "  'theta[0,15]': 0.8174065018579895,\n",
       "  'theta[0,16]': 0.040020199411211724,\n",
       "  'theta[0,1]': 0.7075789874745368,\n",
       "  'theta[0,2]': 0.9748687817077335,\n",
       "  'theta[0,3]': 0.8365340382170593,\n",
       "  'theta[0,4]': 0.08485148126978294,\n",
       "  'theta[0,5]': 0.010213792630017544,\n",
       "  'theta[0,6]': 0.33257814943064623,\n",
       "  'theta[0,7]': 0.2631023437401141,\n",
       "  'theta[0,8]': 0.9194089400406694,\n",
       "  'theta[0,9]': 0.981317286609801,\n",
       "  'theta[1,0]': 0.45449833029198994,\n",
       "  'theta[1,10]': 0.10653820025758994,\n",
       "  'theta[1,11]': 0.3649345813185623,\n",
       "  'theta[1,12]': 0.9704697969692907,\n",
       "  'theta[1,13]': 0.339400836175172,\n",
       "  'theta[1,14]': 0.3220172938299579,\n",
       "  'theta[1,15]': 0.3425599888384405,\n",
       "  'theta[1,16]': 0.8356121197540083,\n",
       "  'theta[1,1]': 0.13652335185991693,\n",
       "  'theta[1,2]': 0.3550542761767105,\n",
       "  'theta[1,3]': 0.4499621719263331,\n",
       "  'theta[1,4]': 0.5739404059191838,\n",
       "  'theta[1,5]': 0.5735370947754022,\n",
       "  'theta[1,6]': 0.9141328053344362,\n",
       "  'theta[1,7]': 0.9269105136428579,\n",
       "  'theta[1,8]': 0.15979332495378337,\n",
       "  'theta[1,9]': 0.1455715320237873,\n",
       "  'theta[2,0]': 0.26886407680569396,\n",
       "  'theta[2,10]': 0.18094224117069646,\n",
       "  'theta[2,11]': 0.47873504142546547,\n",
       "  'theta[2,12]': 0.8828994732048471,\n",
       "  'theta[2,13]': 0.7369362329195503,\n",
       "  'theta[2,14]': 0.43001308563250656,\n",
       "  'theta[2,15]': 0.8551735189700829,\n",
       "  'theta[2,16]': 0.3746475192051309,\n",
       "  'theta[2,1]': 0.5346892928430168,\n",
       "  'theta[2,2]': 0.22984674982259148,\n",
       "  'theta[2,3]': 0.24843779340934868,\n",
       "  'theta[2,4]': 0.5935565351438259,\n",
       "  'theta[2,5]': 0.9400787639191654,\n",
       "  'theta[2,6]': 0.7715596926052253,\n",
       "  'theta[2,7]': 0.14045800211733275,\n",
       "  'theta[2,8]': 0.3435125685648743,\n",
       "  'theta[2,9]': 0.05893657182455881},\n",
       " {'k': 11.0,\n",
       "  'theta[0,0]': 0.2513978561884933,\n",
       "  'theta[0,10]': 0.4581675472029948,\n",
       "  'theta[0,1]': 0.21993194684061876,\n",
       "  'theta[0,2]': 0.7264286969051923,\n",
       "  'theta[0,3]': 0.2728739280294696,\n",
       "  'theta[0,4]': 0.375154182222539,\n",
       "  'theta[0,5]': 0.34788799693072503,\n",
       "  'theta[0,6]': 0.5301308517622717,\n",
       "  'theta[0,7]': 0.8435498416361362,\n",
       "  'theta[0,8]': 0.43118546248612416,\n",
       "  'theta[0,9]': 0.5892439263011554,\n",
       "  'theta[1,0]': 0.01850947997510715,\n",
       "  'theta[1,10]': 0.7747042328765068,\n",
       "  'theta[1,1]': 0.9295870318587716,\n",
       "  'theta[1,2]': 0.14712383391815625,\n",
       "  'theta[1,3]': 0.6872792613875378,\n",
       "  'theta[1,4]': 0.34743511748019856,\n",
       "  'theta[1,5]': 0.18546670145085054,\n",
       "  'theta[1,6]': 0.4791956313165762,\n",
       "  'theta[1,7]': 0.19951408496509476,\n",
       "  'theta[1,8]': 0.01406274259450546,\n",
       "  'theta[1,9]': 0.37991456061125733,\n",
       "  'theta[2,0]': 0.9434866639055607,\n",
       "  'theta[2,10]': 0.8340056141662335,\n",
       "  'theta[2,1]': 0.5423965392126524,\n",
       "  'theta[2,2]': 0.999709576101717,\n",
       "  'theta[2,3]': 0.3662146302379704,\n",
       "  'theta[2,4]': 0.2908702310090892,\n",
       "  'theta[2,5]': 0.6756068348132871,\n",
       "  'theta[2,6]': 0.5602014120308008,\n",
       "  'theta[2,7]': 0.06558722446427595,\n",
       "  'theta[2,8]': 0.1814198821591568,\n",
       "  'theta[2,9]': 0.998093275616343})"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [None for r in range(12)]\n",
    "try:\n",
    "    del r\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "mu = [None for r in range(12)]\n",
    "try:\n",
    "    del r\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "sigma = [None for r in range(12)]\n",
    "try:\n",
    "    del r\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "Y = [None for r in range(12)]\n",
    "try:\n",
    "    del r\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "for r in range(12):\n",
    "    W = Z[:, :, r:r + 24].copy(order = 'F')\n",
    "\n",
    "    W_log = np.log1p(W).copy(order = 'F')\n",
    "\n",
    "    W_log_mean = W_log.mean(axis = -1).copy(order = 'F')\n",
    "    W_log_std = W_log.std(axis = -1, ddof = 1).copy(order = 'F')\n",
    "\n",
    "    W_log_std[np.isnan(W_log_std) | np.isinf(W_log_std) | np.isclose(1, 1 + W_log_std)] = 1\n",
    "\n",
    "    W_log_norm = (\n",
    "        (W_log - W_log_mean.reshape((4, 21, 1))) / W_log_std.reshape((4, 21, 1))\n",
    "    )\n",
    "\n",
    "    X[r] = W_log_norm\n",
    "    mu[r] = W_log_mean\n",
    "    sigma[r] = W_log_std\n",
    "    Y[r] = Z[:, :, r + 24].copy(order = 'F')\n",
    "\n",
    "    del W\n",
    "    del W_log\n",
    "    del W_log_mean\n",
    "    del W_log_std\n",
    "    del W_log_norm\n",
    "try:\n",
    "    del r\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "X = tuple(X)\n",
    "mu = tuple(mu)\n",
    "sigma = tuple(sigma)\n",
    "Y = tuple(Y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array(X, order = 'F')\n",
    "mu = np.array(mu, order = 'F')\n",
    "sigma = np.array(sigma, order = 'F')\n",
    "Y = np.array(Y, order = 'F')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_year (X, mu, sigma, k, T0, p, theta):\n",
    "    Y = np.array(\n",
    "        [\n",
    "            cp_score(\n",
    "                X[i],\n",
    "                k = k,\n",
    "                T0 = T0,\n",
    "                predict = (\n",
    "                    lambda A : ExponentialSmooth().fit(\n",
    "                        A.T,\n",
    "                        p\n",
    "                    ).predict(1, theta)\n",
    "                )\n",
    "            ) for i in range(12)\n",
    "        ],\n",
    "        order = 'F'\n",
    "    )\n",
    "    try:\n",
    "        del i\n",
    "    except (NameError, UnboundLocalError):\n",
    "        pass\n",
    "    return ((sigma * Y) + mu).copy(order = 'F')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optim2 (\n",
    "    space,\n",
    "    X,\n",
    "    mu,\n",
    "    sigma,\n",
    "    Y,\n",
    "    loss,\n",
    "    k = None,\n",
    "    T0 = None,\n",
    "    p = None,\n",
    "    theta = None\n",
    "):\n",
    "    if 'k' in space:\n",
    "        assert k is None\n",
    "    if 'T0' in space:\n",
    "        assert T0 is None\n",
    "    assert (p is None) == ('p' in space)\n",
    "    assert (theta is None) == ('theta' in space)\n",
    "\n",
    "    if theta is None:\n",
    "        theta = space['theta']\n",
    "    if isinstance(theta, numbers.Complex):\n",
    "        if isinstance(theta, numbers.Integral):\n",
    "            theta = np.zeros(theta, dtype = float, order = 'F')\n",
    "            for i in iter(range(int(theta.size))):\n",
    "                theta[i] = space['theta[{0:d}]'.format(i)]\n",
    "        else:\n",
    "            theta = space['theta']\n",
    "    else:\n",
    "        if not isinstance(theta, tuple):\n",
    "            if isinstance(theta, np.ndarray):\n",
    "                theta = copy.deepcopy((theta, theta, theta))\n",
    "            else:\n",
    "                if hasattr(theta, '__iter__'):\n",
    "                    theta = copy.deepcopy(tuple(theta))\n",
    "                else:\n",
    "                    theta = copy.deepcopy((theta, theta, theta))\n",
    "        else:\n",
    "            theta = copy.deepcopy(theta)\n",
    "        if len(theta) != 3:\n",
    "            theta = copy.deepcopy((theta, theta, theta))\n",
    "        theta = list(theta)\n",
    "        for i in iter(range(3)):\n",
    "            try:\n",
    "                theta[i] = np.zeros(theta[i], dtype = float, order = 'C')\n",
    "                idx = np.array(\n",
    "                    np.meshgrid(\n",
    "                        *tuple(\n",
    "                            np.arange(theta[i].shape[j], dtype = int)\n",
    "                                for j in iter(range(int(theta[i].ndim)))\n",
    "                        )\n",
    "                    )\n",
    "                ).T.reshape((-1, theta[i].ndim)).copy(order = 'C')\n",
    "                for j in iter(range(int(idx.shape[0]))):\n",
    "                    theta[i][tuple(idx[j])] = (\n",
    "                        space[\n",
    "                            'theta[{0:s}]'.format(\n",
    "                                ','.join(\n",
    "                                    str(r)\n",
    "                                        for r in np.concatenate(\n",
    "                                            ([i], idx[j])\n",
    "                                        )\n",
    "                                )\n",
    "                            )\n",
    "                        ]\n",
    "                    )\n",
    "                try:\n",
    "                    del j\n",
    "                except (NameError, UnboundLocalError):\n",
    "                    pass\n",
    "            except (TypeError, ValueError, AttributeError):\n",
    "                assert (\n",
    "                    hasattr(theta[i], '__iter__') or\n",
    "                    hasattr(theta[i], '__array__')\n",
    "                )\n",
    "                theta[i] = np.array(theta[i], copy = True)\n",
    "        try:\n",
    "            del i\n",
    "        except (NameError, UnboundLocalError):\n",
    "            pass\n",
    "        theta = tuple(theta)\n",
    "\n",
    "    Y_predicted = None\n",
    "\n",
    "    with warnings.catch_warnings():\n",
    "        warnings.filterwarnings('error')\n",
    "        with np.errstate(divide = 'raise', invalid = 'raise'):\n",
    "            try:\n",
    "                Y_predicted = predict_year(\n",
    "                    X,\n",
    "                    mu,\n",
    "                    sigma,\n",
    "                    int(round(space['k'])) if 'k' in space else k,\n",
    "                    int(round(space['T0'])) if 'T0' in space else T0,\n",
    "                    int(round(space['p'])) if 'p' in space else p,\n",
    "                    theta\n",
    "                )\n",
    "            except (\n",
    "                TypeError,\n",
    "                ValueError,\n",
    "                FloatingPointError,\n",
    "                np.linalg.LinAlgError,\n",
    "                RuntimeWarning\n",
    "            ):\n",
    "                return {'status' : ho.STATUS_FAIL}\n",
    "\n",
    "    return {'loss' : loss(Y, Y_predicted), 'status' : ho.STATUS_OK}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [01:22<00:00,  1.23it/s, best loss: 0.9971455722566905]\n",
      "100%|██████████| 100/100 [01:22<00:00,  1.20it/s, best loss: 0.9971759472798843]\n",
      "100%|██████████| 100/100 [01:21<00:00,  1.21it/s, best loss: 0.997193706582929]\n",
      "100%|██████████| 100/100 [01:23<00:00,  1.17it/s, best loss: 0.9971940069087487]\n",
      "100%|██████████| 100/100 [01:24<00:00,  1.15it/s, best loss: 0.9971576575397831]\n",
      "100%|██████████| 100/100 [01:27<00:00,  1.11it/s, best loss: 0.9971736000917988]\n",
      "100%|██████████| 100/100 [01:25<00:00,  1.14it/s, best loss: 0.9971699718163951]\n",
      "100%|██████████| 100/100 [01:28<00:00,  1.11it/s, best loss: 0.9971671070642018]\n",
      "100%|██████████| 100/100 [01:29<00:00,  1.09it/s, best loss: 0.9971945690449694]\n",
      "100%|██████████| 100/100 [01:32<00:00,  1.05it/s, best loss: 0.9971790521355937]\n",
      "100%|██████████| 100/100 [01:32<00:00,  1.07it/s, best loss: 0.997144403155174]\n",
      "100%|██████████| 100/100 [01:33<00:00,  1.04it/s, best loss: 0.9971788582969737]\n",
      "100%|██████████| 100/100 [01:36<00:00,  1.02it/s, best loss: 0.9971456712865927]\n",
      "100%|██████████| 100/100 [01:39<00:00,  1.02s/it, best loss: 0.9971651311771245]\n"
     ]
    }
   ],
   "source": [
    "best = [dict() for k in range(6, 20)]\n",
    "try:\n",
    "    del k\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "for k in range(6, 20):\n",
    "    trials = ho.Trials()\n",
    "    best[k - 6] = ho.fmin(\n",
    "        fn = lambda S : (\n",
    "            optim2(\n",
    "                update_theta(S, k),\n",
    "                X = X,\n",
    "                mu = mu,\n",
    "                sigma = sigma,\n",
    "                Y = Y,\n",
    "                loss = rel_mse,\n",
    "                k = k,\n",
    "                T0 = default_space['T0'],\n",
    "                p = default_space['p']\n",
    "            )\n",
    "        ),\n",
    "        space = space,\n",
    "        algo = ho.tpe.suggest,\n",
    "        max_evals = 100,\n",
    "        trials = trials,\n",
    "        show_progressbar = True\n",
    "    )\n",
    "try:\n",
    "    del k\n",
    "except (NameError, UnboundLocalError):\n",
    "    pass\n",
    "\n",
    "best = tuple(best)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}