{

 "cells": [

  {

   "cell_type": "markdown",

   "metadata": {},

   "source": [

    "## Tensorflow2 version\n",

    "# Practical Deep Learning Guide for Genomic Prediction\n",

    "## A Keras based guide to implement deep learning using tensorflow 2\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": null,

   "metadata": {},

   "outputs": [],

   "source": [

    "## Tensorflow 2 now incorporates keras\n",

    "## Hyper - Parameter optimization with keras tuner\n",

    "\n",

    "### M Perez-Enciso & LM Zingaretti\n",

    "### miguel.perez@uab.es, m.lau.zingaretti@gmail.com\n",

    "\n",

    "### thanks also to I Vourlaki (ibourlaki@gmail.com)\n",

    "\n",

    "### If you find this resource useful, please cite: \n",

    "### Pérez-Enciso M, Zingaretti LM. 2019. A Guide on Deep Learning for Complex Trait Genomic Prediction. Genes, 10, 553.\n",

    "### Zingaretti LM, Gezan SA, Ferrão LFV, Osorio LF, Monfort A, Muñoz PR, Whitaker VM, Pérez-Enciso M. 2020. Exploring Deep Learning for Complex Trait Genomic Prediction in Polyploid Outcrossing Species. Frontiers in Plant Science 11:25\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 1,

   "metadata": {

    "scrolled": true

   },

   "outputs": [],

   "source": [

    "# main modules needed\n",

    "import pandas as pd\n",

    "import numpy as np\n",

    "import seaborn as sns\n",

    "from sklearn import linear_model\n",

    "from sklearn.model_selection import train_test_split\n",

    "from matplotlib import pyplot as plt\n",

    "from scipy import stats\n",

    "from sklearn.decomposition import PCA\n",

    "from sklearn.preprocessing import StandardScaler\n",

    "from sklearn.preprocessing import scale"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 2,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "tensorflow: 2.4.1\n",

      "keras: 2.4.0\n",

      "kerastuner: 1.0.2\n"

     ]

    }

   ],

   "source": [

    "# DL modules\n",

    "# tensorflow\n",

    "import tensorflow as tf\n",

    "print('tensorflow: %s' % tf.__version__)\n",

    "# keras\n",

    "from tensorflow import keras\n",

    "print('keras: %s' % keras.__version__)\n",

    "import kerastuner as kt\n",

    "print('kerastuner: %s' % kt.__version__)"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 3,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "(479, 1279) (479,)\n",

      "(120, 1279) (120,)\n",

      "       min max mean sd\n",

      "Train: -2.41866172921982 3.27892080508434 0.03656243695565241 0.9744612298270398\n",

      "Test: -2.28909755016522 2.52690454198022 -0.14594506084797887 1.088160002689451\n"

     ]

    },

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEICAYAAACktLTqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVtElEQVR4nO3df5BdZZ3n8ffHJCQK0WCIDBA0sRRWtBDGXkFcFWQs0eDAuqi4OhuU2SzWKLKjhSjlwlbJVtiZ9Qe6O8oMCO74AwVnA+KPQRcqWCoYGEQwMDAShqYCxCgBVkEC3/2jD0wTOul039u53U/er6pUnx/POed7LsmHp59z7jmpKiRJbXnGoAuQJPWf4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhrRknynSTLB11HvyW5IMknBl2H2mG4a8oleWjUn8eT/G7U/Lsmsq+qelNVXdhDLXsnGd7Kukryosnue9R+zkzyt73uZxv7vyrJn07V/tWG2YMuQO2rqt2emE6yDvjTqvr+lu2SzK6qzVNczpuB707xMaSBs+eugUlyeJLhJB9Jcg/wxSS7J/lWkg1JftNNLx61zZO91iQnJPlhkr/s2t6R5E3jHPbNwLfHqGV1N/mz7jeKd3TLj05yQ5L7k/woyYGjtvlIkruTPJjk1iRHJjkK+Bjwjm4/P9vKuR+c5Ppu24uAeaPWbfUzSHIW8Brgc93+P9ct/0ySu5I8kOS6JK8Z53NQ4wx3DdofAM8FXgCsYOTv5Be7+ecDvwM+t43tDwFuBfYA/jtwXpKM1TDJHOC1wBVbrquq13aTL6+q3arqoiQHA+cD/wlYCHwBuDTJ3CT7A+8H/nVVzQfeCKyrqu8C/w24qNvPy8eoYxfg/wD/uzv3bwD/blSTrX4GVXU6cDXw/m7/7++2+SlwULe/rwDfSDIP7bQMdw3a48AZVfVIVf2uqjZW1SVV9duqehA4C3jdNra/s6r+uqoeAy4E9gL23Erb1wI/6/a7PVYAX6iqa6rqsW6s/xHgUOAxYC5wQJI5VbWuqv5pO/d7KDAH+HRVPVpVFzMSzgBM4jOgqv62225zVf2Prrb9t7MeNchw16BtqKqHn5hJ8qwkX0hyZ5IHgNXAgiSztrL9PU9MVNVvu8ndttJ2zCGZbXgB8KFuSOb+JPcD+wJ7V9XtwCnAmcB9Sb6WZO/t3O/ewN311Eey3vnExCQ+A5J8OMnaJJu6Op/DyG8z2kkZ7hq0LZ85/SFGepyHVNWzGeltA4w51DJBEw33u4CzqmrBqD/PqqqvAlTVV6rq3zDyP4ECzu62G+852uuBfbYYPnr+qOnxPoOn7L8bXz8VeDuwe1UtADbRn89MM5ThrulmPiNjzPcneS5wRj92mmQpMLeq1m6j2b3AC0fN/zVwUpJDMmLXJMuSzE+yf5LXJ5kLPNzV/Pio/SxJsrV/Xz8GNgMnJ5mT5K3AK0etH+8z2LLO+d3+NgCzk/wX4NnbOE/tBAx3TTefBp4J/Ar4Cf27bXEZ4/fazwQu7IZg3l5Va4D/yMjFzN8AtwMndG3nAiu7Ou8Bngd8tFv3je7nxiTXb3mQqvo98NZuX78G3gF8c1STT7Ptz+AzwHHdnTTnAN/r2vwjI8M7DzPyW4d2YvFNTNoZJPk28LmqmsiwjDRj2XPXzuIq4MpBFyHtKPbcJalB9twlqUHT4tkye+yxRy1ZsmTQZUjSjHLdddf9qqoWjbVuWoT7kiVLWLNmzaDLkKQZJcmdW1vnsIwkNchwl6QGGe6S1KBpMeYuSRP16KOPMjw8zMMPPzx+4xlu3rx5LF68mDlz5mz3Noa7pBlpeHiY+fPns2TJErbyCP8mVBUbN25keHiYpUuXbvd2DstImpEefvhhFi5c2HSwAyRh4cKFE/4NxXCXNGO1HuxPmMx5Gu6S1CDH3CU1Yclpl/d1f+tWLtvm+o0bN3LkkUcCcM899zBr1iwWLRr5sui1117LLrvsstVt16xZw5e+9CXOOeec/hW8hXHDPcn5wNHAfVX1sm7ZXwBvAX4P/BPwnqq6v1v3UeBERt4xeXJVfW9qStfOpN//cCdivH/k2jktXLiQG264AYAzzzyT3XbbjQ9/+MNPrt+8eTOzZ48dsUNDQwwNDU1pfdszLHMBcNQWy64AXlZVBzLygoCPAiQ5ADgeeGm3zf/a1nsfJaklJ5xwAieddBKHHHIIp556Ktdeey2vetWrOPjggznssMO49dZbAbjqqqs4+uijgZH/Mbz3ve/l8MMP54UvfGHfevPj9tyranWSJVss+/tRsz8BjuumjwG+VlWPAHckuZ2R14f9uC/VStI0Nzw8zI9+9CNmzZrFAw88wNVXX83s2bP5/ve/z8c+9jEuueSSp21zyy23cOWVV/Lggw+y//778773vW9C97SPpR9j7u8FLuqm92Ek7J8w3C17miQrgBUAz3/+88dqIkkzztve9jZmzRoZsNi0aRPLly/ntttuIwmPPvromNssW7aMuXPnMnfuXJ73vOdx7733snjx4p7q6OlumSSnM/Ji3i9PdNuqOreqhqpq6ImLEJI00+26665PTn/84x/niCOO4KabbuKyyy7b6r3qc+fOfXJ61qxZbN68uec6Jt1zT3ICIxdaj6x/eZ3T3cC+o5ot7pZJ0k5n06ZN7LPPyODFBRdcsEOPPalwT3IUcCrwuqr67ahVlwJfSfJJYG/gxcC1PVcpSeOYjnc1nXrqqSxfvpxPfOITLFu2Y+sb9x2qSb4KHA7sAdwLnMHI3TFzgY1ds59U1Uld+9MZGYffDJxSVd8Zr4ihoaHyZR3aFm+F1JbWrl3LS17ykkGXscOMdb5JrquqMe+p3J67Zd45xuLzttH+LOCs8fYrSZo6Pn5AkhpkuEtSgwx3SWqQ4S5JDTLcJalBPvJXUhvOfE6f97dpm6t7eeQvjDw8bJddduGwww7rT71bMNwlaRLGe+TveK666ip22223KQt3h2UkqU+uu+46Xve61/GKV7yCN77xjaxfvx6Ac845hwMOOIADDzyQ448/nnXr1vH5z3+eT33qUxx00EFcffXVfa/Fnrsk9UFV8YEPfIBVq1axaNEiLrroIk4//XTOP/98Vq5cyR133MHcuXO5//77WbBgASeddNKEe/sTYbhLUh888sgj3HTTTbzhDW8A4LHHHmOvvfYC4MADD+Rd73oXxx57LMcee+wOqcdwl6Q+qCpe+tKX8uMfP/3dRJdffjmrV6/msssu46yzzuLnP//5lNfjmLsk9cHcuXPZsGHDk+H+6KOPcvPNN/P4449z1113ccQRR3D22WezadMmHnroIebPn8+DDz44ZfXYc5fUhnFuXZxqz3jGM7j44os5+eST2bRpE5s3b+aUU05hv/32493vfjebNm2iqjj55JNZsGABb3nLWzjuuONYtWoVn/3sZ3nNa17T13oMd0nq0Zlnnvnk9OrVq5+2/oc//OHTlu23337ceOONU1aTwzKS1CDDXZIaZLhLmrHGe5NcKyZznoa7pBlp3rx5bNy4sfmAryo2btzIvHnzJrSdF1QlzUiLFy9meHiYDRs2DLqUKTdv3jwWL148oW0Md0kz0pw5c1i6dOmgy5i2HJaRpAYZ7pLUIMNdkhpkuEtSgwx3SWrQuOGe5Pwk9yW5adSy5ya5Islt3c/du+VJck6S25PcmOQPp7J4SdLYtqfnfgFw1BbLTgN+UFUvBn7QzQO8CXhx92cF8Ff9KVOSNBHjhntVrQZ+vcXiY4ALu+kLgWNHLf9SjfgJsCDJXn2qVZK0nSY75r5nVa3vpu8B9uym9wHuGtVuuFv2NElWJFmTZM3O8A0zSdqRer6gWiMPdpjwwx2q6tyqGqqqoUWLFvVahiRplMmG+71PDLd0P+/rlt8N7Duq3eJumSRpB5psuF8KLO+mlwOrRi3/D91dM4cCm0YN30iSdpBxHxyW5KvA4cAeSYaBM4CVwNeTnAjcCby9a/5t4M3A7cBvgfdMQc2SpHGMG+5V9c6trDpyjLYF/FmvRUnTyZLTLh/IcdetXDaQ46oNfkNVkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhrUU7gn+c9Jbk5yU5KvJpmXZGmSa5LcnuSiJLv0q1hJ0vaZdLgn2Qc4GRiqqpcBs4DjgbOBT1XVi4DfACf2o1BJ0vbrdVhmNvDMJLOBZwHrgdcDF3frLwSO7fEYkqQJmnS4V9XdwF8C/8xIqG8CrgPur6rNXbNhYJ+xtk+yIsmaJGs2bNgw2TIkSWPoZVhmd+AYYCmwN7ArcNT2bl9V51bVUFUNLVq0aLJlSJLG0MuwzB8Bd1TVhqp6FPgm8GpgQTdMA7AYuLvHGiVJEzR7/CZb9c/AoUmeBfwOOBJYA1wJHAd8DVgOrOq1SE0fS067fNAlSNoOvYy5X8PIhdPrgZ93+zoX+Ajw50luBxYC5/WhTknSBPTSc6eqzgDO2GLxL4FX9rJfSVJv/IaqJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDWop3BPsiDJxUluSbI2yauSPDfJFUlu637u3q9iJUnbp9ee+2eA71bVvwJeDqwFTgN+UFUvBn7QzUuSdqBJh3uS5wCvBc4DqKrfV9X9wDHAhV2zC4FjeytRkjRRvfTclwIbgC8m+Yckf5NkV2DPqlrftbkH2HOsjZOsSLImyZoNGzb0UIYkaUu9hPts4A+Bv6qqg4H/xxZDMFVVQI21cVWdW1VDVTW0aNGiHsqQJG1pdg/bDgPDVXVNN38xI+F+b5K9qmp9kr2A+3otUk+15LTLB12CpGlu0j33qroHuCvJ/t2iI4FfAJcCy7tly4FVPVUoSZqwXnruAB8AvpxkF+CXwHsY+R/G15OcCNwJvL3HY0iSJqincK+qG4ChMVYd2ct+JUm98RuqktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSg3p9cJikKTKoRzuvW7lsIMdVf9lzl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJalDP4Z5kVpJ/SPKtbn5pkmuS3J7koiS79F6mJGki+tFz/yCwdtT82cCnqupFwG+AE/twDEnSBPQU7kkWA8uAv+nmA7weuLhrciFwbC/HkCRNXK89908DpwKPd/MLgfuranM3PwzsM9aGSVYkWZNkzYYNG3osQ5I02qTDPcnRwH1Vdd1ktq+qc6tqqKqGFi1aNNkyJElj6OU1e68G/jjJm4F5wLOBzwALkszueu+Lgbt7L1OSNBGT7rlX1UeranFVLQGOB/5vVb0LuBI4rmu2HFjVc5WSpAmZivvcPwL8eZLbGRmDP28KjiFJ2oZehmWeVFVXAVd1078EXtmP/U53g3o7vSSNx2+oSlKDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDJh3uSfZNcmWSXyS5OckHu+XPTXJFktu6n7v3r1xJ0vaY3cO2m4EPVdX1SeYD1yW5AjgB+EFVrUxyGnAa8JHeS9WOsm7ev5/yYyx5+CtTfgxpZzbpnntVra+q67vpB4G1wD7AMcCFXbMLgWN7rFGSNEF9GXNPsgQ4GLgG2LOq1ner7gH23Mo2K5KsSbJmw4YN/ShDktTpOdyT7AZcApxSVQ+MXldVBdRY21XVuVU1VFVDixYt6rUMSdIovYy5k2QOI8H+5ar6Zrf43iR7VdX6JHsB9/VapDRZO+L6AXgNQdPPpMM9SYDzgLVV9clRqy4FlgMru5+reqpQTdpRoSvtrHrpub8a+BPg50lu6JZ9jJFQ/3qSE4E7gbf3VKEkacImHe5V9UMgW1l95GT3K2mwlpx2+cCOvW7lsoEduzV+Q1WSGmS4S1KDDHdJapDhLkkN6uk+9+lgkBd/djRvH5S0vey5S1KDZnzPXdpZ+G1bTYQ9d0lqkOEuSQ0y3CWpQYa7JDXIC6p94C2K8u+Apht77pLUoKZ77vamJO2s7LlLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktSgpu9zlzSzDOrlO+tWLhvIcaeSPXdJapA9d0lP0do3u3fWl4/Yc5ekBk1ZuCc5KsmtSW5PctpUHUeS9HRTMiyTZBbwP4E3AMPAT5NcWlW/mIrjSVIvBnUhF6buYu5U9dxfCdxeVb+sqt8DXwOOmaJjSZK2MFUXVPcB7ho1PwwcMrpBkhXAim72oSS39ruI9L6LPYBf9b6baaXFc4I2z8tz6oujd8RBJn1eObun475gaysGdrdMVZ0LnDuo42+PJGuqamjQdfRTi+cEbZ6X5zRzTMfzmqphmbuBfUfNL+6WSZJ2gKkK958CL06yNMkuwPHApVN0LEnSFqZkWKaqNid5P/A9YBZwflXdPBXHmmLTethoklo8J2jzvDynmWPanVeqatA1SJL6zG+oSlKDDHdJapDhPo4kf5HkliQ3Jvm7JAsGXVOvkrwtyc1JHk8yrW7fmqgWH3OR5Pwk9yW5adC19EuSfZNcmeQX3d+9Dw66pl4lmZfk2iQ/687pvw66ptEM9/FdAbysqg4E/hH46IDr6YebgLcCqwddSC9GPebiTcABwDuTHDDYqvriAuCoQRfRZ5uBD1XVAcChwJ818N/qEeD1VfVy4CDgqCSHDrakf2G4j6Oq/r6qNnezP2Hknv0ZrarWVlXfvxE8AE0+5qKqVgO/HnQd/VRV66vq+m76QWAtI99kn7FqxEPd7Jzuz7S5Q8Vwn5j3At8ZdBF60liPuZjRgbEzSLIEOBi4ZsCl9CzJrCQ3APcBV1TVtDknX9YBJPk+8AdjrDq9qlZ1bU5n5FfLL+/I2iZre85J2tGS7AZcApxSVQ8Mup5eVdVjwEHdtbi/S/KyqpoW10oMd6Cq/mhb65OcwMjTh46sGfLFgPHOqRE+5mIGSTKHkWD/clV9c9D19FNV3Z/kSkaulUyLcHdYZhxJjgJOBf64qn476Hr0FD7mYoZIEuA8YG1VfXLQ9fRDkkVP3D2X5JmMvL/iloEWNYrhPr7PAfOBK5LckOTzgy6oV0n+bZJh4FXA5Um+N+iaJqO70P3EYy7WAl+foY+5eIokXwV+DOyfZDjJiYOuqQ9eDfwJ8Pru39ENSd486KJ6tBdwZZIbGeloXFFV3xpwTU/y8QOS1CB77pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNej/A5w2Jcrt6boNAAAAAElFTkSuQmCC\n",

      "text/plain": [

       "<Figure size 432x288 with 1 Axes>"

      ]

     },

     "metadata": {

      "needs_background": "light"

     },

     "output_type": "display_data"

    },

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEICAYAAAC6fYRZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABPNklEQVR4nO2de5xUdd3435+ZnV12EV0ui7gLCFrBI5EQaCZWrpr0qOA+WlmSeHl6fHWxEhPE1MTSR5QS7Wc9PT6mRlGRl1ZIi1TAW5FAQEhCqdx2VmQBBy+7wuzO9/fHmTOcOXPOmTOXndt+36/XvnbnzLl8Z/Z7Pufz/VxFKYVGo9FoKpNAsQeg0Wg0mt5DC3mNRqOpYLSQ12g0mgpGC3mNRqOpYLSQ12g0mgpGC3mNRqOpYLSQ15QNIqJE5APFHkdvIiLfEZH7Pd6fISJ/KuSYNOWNFvIa34jIdhHpEpF3ReRNEXlIRI6wvD9VRJ4TkXdEpENEnhWR6bZznB4X1tcV/hOUPkqp/1ZKfRlAREbFv6sqy/uLlVJnF2+EmnJDC3lNpkxTSh0BfBSYDNwIICKfBR4GFgHDgaOB7wLTbMdfCuwHZhZqwBpNX0YLeU1WKKXCwB+AD4uIAHcB31dK3a+UOqCUiimlnlVK/Zd5jIj0Bz4LfB34oIhM9rqGiMwWkTdEpF1ErrC9VyMiPxCRnfFVxU9FpNby/vkiskFE3haR10TkM/HtjSKyVET2i8irImId3zwReVhEfhlfjWwSkQ+JyPUiskdEdonI2Zb9V4nI7SLyUvw6j4vIIMv700Vks4hE4vv+m+W960QkHL/OVhE50zKGX8Z3ey7+OxJfPX1cRC4TkRcs5zlVRNaIyIH471Nt4/u+iLwYv86fRGSI13euqTy0kNdkhYiMAM4B1gNjgBHAI2kOuwB4F0PjX46h1bud/zPAtcCngQ8CZ9l2mQ98CJgAfABowlg5ICInY6woZgP1wCeB7fHjfgO0AY0YD5z/FpEzLOedBvwCGBj/bMsx7pMm4HvA/9rGMRO4AjgG6AZ+FB/Dh4BfA1cDDcCTwDIRqRaRMcBVwElKqQHAVMv4rHwy/rteKXWEUuovtu9oEPBE/JqDMR60T4jIYMtuFwOXA0OBaozvVNOXUErpH/3j6wdDEL0LRIAdwE+AWmAKoIB+aY5/Grg7/vcXgQ4g5LLvA8B8y+sPxa/xAUCA94DjLe9/HNgW//t/gYUO5xwB9AADLNtuBx6K/z0PeMry3rT45w3GXw+Ij6E+/nqVbYwnAIeAIHAT8FvLewEgDJwe/wx7MB5cIdsY5wG/jP89Kn69Ksv7lwEvxP++BHjJdvxfgMss47vR8t7XgD8Wex7pn8L+aE1ekyktSql6pdSxSqmvKaW6gH3x945xOyiu+TcDi+ObHgf6Aee6HNII7LK83mH5uwGoA9bFTSER4I/x7WAI89dczrlfKfWO7bxNltdvWv7uAvYqpXosrwGOsOxjH2MIGBK/VmLMSqlYfN8mpdSrGBr+PGCPiPxGRBodxpuOpGu4fJ7dlr87bWPX9AG0kNfkg60YAuxCj30uwZhvy0RkN/A6hpB3M9m8gSGsTUZa/t6LIXDHxR849Uqpo5ThECY+luMdztkODBKRAbbzhj3GnQ77GKPx8bUDx5pvxP0WI8xrKaV+pZQ6Lb6PAu5wOHe6ErFJ17CMIZfPo6kwtJDX5IxSSgHXADeJyOUicqSIBETkNBG5L77bpcAtGDZ08+dC4BybDdnkt8BlInKCiNQBN1uuFwP+D1goIkMBRKRJRKbGd/kZcLmInBkfR5OIjFVK7QL+DNwuIv1E5CPAfwK/JHu+ZBnj94BH4pr/b4Fz42MIAd8GDgJ/FpExInKGiNQA72M8sGIO5+6Ibz/O5dpPAh8SkYtFpEpELsIwGf0+h8+jqTC0kNfkBaXUI8BFGE7Idgyzx63A4yJyCobG+WOl1G7Lz1LgVQz7vP18fwDuBlbE91lh2+W6+PbVIvI2hr1/TPzYlzCcjQuBA8CzHNZ4v4hh624HfgfcrJR6OoeP/gvgIQyzSD/gm/ExbAW+BPw/DM1+Gkb46SGgBsNxvDd+3FDgeofvoBO4DXgxbpY6xfb+PuA8jAfIPmAOcJ5Sam8On0dTYYihhGk0mkwRkVUYTlLXDFWNpthoTV6j0WgqGC3kNRqNpoLR5hqNRqOpYLQmr9FoNBVMVfpd8s+QIUPUqFGjinFpjUajKVvWrVu3VynVkH7PwxRFyI8aNYq1a9cW49IajUZTtoiIPcM5Ldpco9FoNBWMFvIajUZTwWghr9FoNBVMUWzyGo1GkynRaJS2tjbef//9Yg+l1+nXrx/Dhw8nFArlfC4t5DUaTVnQ1tbGgAEDGDVqFEZRz8pEKcW+fftoa2tj9OjROZ9PC3mNJk+0rg+zYPlW2iNdNNbXMnvqGFomNqU/UOOL999/v+IFPICIMHjwYDo6OvJyPi3kNZo80Lo+zPWPbaIravQXCUe6uP6xTQBa0OeRShfwJvn8nNrxqtHkgQXLtyYEvElXtIcFy7cWaUQajYHW5DWaPNAe6cpouxva5FO67Nu3jzPPPBOA3bt3EwwGaWgwkk9feuklqqurXY9du3YtixYt4kc/+lFBxmpFC3mNJg801tcSdhDojfW1vs+hTT6lzeDBg9mwYQMA8+bN44gjjuDaa69NvN/d3U1VlbNInTx5MpMnTy7EMFPQ5hqNJg/MnjqG2lAwaVttKMjsqWN8n0ObfPJL6/owU+avYPTcJ5gyfwWt6/Pf+vayyy7jK1/5Ch/72MeYM2cOL730Eh//+MeZOHEip556Klu3Gv+7VatWcd555wHGA+KKK67g9NNP57jjjut17V5r8hpNHjA17VxMLfky+WgKuypqa2vjz3/+M8FgkLfffpvnn3+eqqoqnn76ab7zne/w6KOPphyzZcsWVq5cyTvvvMOYMWP46le/mpeYeCe0kNdo8kTLxKacBEg+TD4aA69VUb6F/Oc+9zmCQWMVd+DAAS699FL+9a9/ISJEo1HHY84991xqamqoqalh6NChvPnmmwwfPjyv4zLR5hqNpkTIh8lHY1DIVVH//v0Tf9900000Nzfz8ssvs2zZMtfs3JqamsTfwWCQ7u7uvI/LRAt5jaZEaJnYxO0XjKepvhYBmupruf2C8drpmgVuq5/eXhUdOHCApibj//XQQw/16rX8os01Gk0JkavJR2Mwe+qYJJs8FGZVNGfOHC699FJuvfVWzj333F69ll+K0uN18uTJSjcN0Wg0mfDKK6/wb//2b773L/ecA6fPKyLrlFIZxWJqTV6j0VQkelVkoG3yGo1GU8FoIa/RaDQVTF6EvIjMEpHNIvKyiPxaRPrl47wajUajyY2chbyINAHfBCYrpT4MBIEv5HpejUaj0eROvsw1VUCtiFQBdUB7ns6r0Wg0mhzIObpGKRUWkR8AO4Eu4E9KqT/Z9xORK4ErAUaOHJnrZTUajaag5FJqGIwiZdXV1Zx66qm9PlYrOQt5ERkInA+MBiLAwyLyJaXUL637KaXuA+4DI04+m2uVe9yrRqMpX9KVGk7HqlWrOOKIIwou5PNhrjkL2KaU6lBKRYHHgLx/ihtbN3Hg+St5dsSneX38eTw74tMceP5KbmzdlO9LaTSaSmDbYmgdBb8KGL+3Lc77JdatW8enPvUpJk2axNSpU3njjTcA+NGPfsQJJ5zARz7yEb7whS+wfft2fvrTn7Jw4UImTJjA888/n/exuJGPZKidwCkiUodhrjkTyGs6a+v6MB/cOZeZg5/EbH1YRYyZg59k0U5oXX+f1ug1pcm2xbDxBujcCXUj4cTbYPSMYo+q8tm2GF66Eno6jdedO4zXkLfvXynFN77xDR5//HEaGhpYsmQJN9xwAw888ADz589n27Zt1NTUEIlEqK+v5ytf+UrG2n8+yIdN/q8i8gjwN6AbWE/cLJMvFizfyrMj/oi9t60IzBj8Rz6wZAO3LNvMzdPGaWGvKR0KIGg0Lmy84fD3btLTaWzP03d/8OBBXn75ZT796U8bp+/p4ZhjjgHgIx/5CDNmzKClpYWWlpa8XC9b8lLWQCl1M3BzPs7lRHuki+CImON7QWK8Pn4a7dEh3Ll8Jmt3XMqtLeN7aygajX8KIGg0LnTuzGx7FiilGDduHH/5y19S3nviiSd47rnnWLZsGbfddhubNhXPrFwWGa+N9bX0uAxVBAKiGF7dwfzh9/LOlp/3SpsvjSZjCiBoNC7UuUTwuW3PgpqaGjo6OhJCPhqNsnnzZmKxGLt27aK5uZk77riDAwcO8O677zJgwADeeeedvF3fL2Uh5GdPHcPifZ8hXcHMusBBZg9bxC3LNhdmYBqNFwUQNBoXTrwNgnXJ24J1xvY8EQgEeOSRR7juuus48cQTmTBhAn/+85/p6enhS1/6EuPHj2fixIl885vfpL6+nmnTpvG73/2uLB2vvU7LxCZu3DGfRTsNG3wQw3Rjt9EDNIb28lZnlNb1YVomNiXCLifzJNc3/oKjqzoQ7QDTFIITb0u2yUPeBY3GBfPe7iWn97x58xJ/P/fccynvv/DCCynbPvShD/H3v/89L9fPhLIQ8gC3toxn1NyvcXP71wB4YezlDK/uSNmvPToEINHh/vrHNvHp/k9z+/B7qQscNHbSDjBNIehlQaNJw+gZ+rumjIQ8GO3QzEbHd+6eyXyr4AY6YzXcuXsmYDhr5y3dTFe0hznDFiXtB2gHmKYwVKig0YmJ5UNZ2ORNrI2Ol0aamdt2FW2HGogpoe1QA3PbrmJppBmAo2pDRLqMTumNob3OJ9QOMI0mY1rXh7n+sU2EI10oIBzp4vrHNhUk4KEYneyKQT4/Z1lp8qamsGD5VsKRLpZFmhNCfXr9SuYMW8TdI+6iPTqEFW+fxBnHrqExtJcYQgCHL007wDSajFmwfGtS71SArmgPC5Zv7VVtvl+/fuzbt4/BgwcjTg65CkEpxb59++jXLz8V28tKyENySy9zyTiJJ5NMN8OrO7jEkh0bQKGUzVGrHWAaTVa0x02mfrfni+HDh9PW1kZHR6ovrtLo168fw4cPz8u5yk7IW0kI/NYroDPZ5u6UHdutAgRQBPprB5hGky2NFt+YfXtvEgqFGD16dK9eoxIpK5u8Kz5t6wEUS0/YBS3btYDXaLJk9tQxhILJWlQoKMyeOqZII9J4UdaavElnVSN13emdPu+HGksuAkBHKWjKEruLq2/4Q8uSstfkW9eH+e6Oi+mM1SRtT3FOB+uoO+mOwg3MB8WMUqhIClBaVmM4XqOx5BssGlOJ3BRNaVH2mrwRafMpDvXEmDNsEY2hvbRHh/D6+42cOmATQYkhEuS1AZ9l5pJjaI88kawxF7EUbLGiFCqSCqn4WA4ru0wdr63rw8xbujkR0hwQiCkj76UUP1+lUfaavDmxlkaaOW3Lgxy3aRl37p7J5CO2UCUxBED10Lh/CZN4MqExz1qygW/Ov5auF//TEAgo6NxB9+r/KpgGWKwohYrEq+JjmVAuKzs3B6vT9tb1YWY/vDEh4MEQ8FC6n6/SKHsh7zSxnDJcawMHubnxPl4Yezmvj5/G82Mv5+bG+6i17Velumh/7hpGz32CKfNX9OoEdLtZAiIFuX5FUQEVH71WdqWENSnRpDYUdHS8Opl2rJTi56s08iLkRaReRB4RkS0i8oqIfDwf5/WD04Rzy3AdFHyH4dUdidLEg4LOZT+HVXUURJNyGjtAj1IlrcmVJBVQ8bFcVnYtE5u4/YLxNNXXIhhml9svGO9odvEz9lL7fJVGvmzy9wB/VEp9VkSqgbp0B+QLexasYBQpcype5hQ774RZ5Ax610ZuHXt7pIuACD02j7G20fukAio+Fiv+PBusSYleuH0m+z6a3iNnTV5EjgI+CfwMQCl1SCkVyfW8mdAysYkX555BU30tCqN4mT3axg17FI61yJlJb2oa5ti3zT+XmEu9Cq3p+GD0DDj5Pqg7FhDj98n3lZXTNRMzSLkwe+oYQgH3EgTl/vnKgXxo8qOBDuBBETkRWAd8Syn1Xh7OnRFWJyzAXSMWUiXObQNNJO7pFyAcbeDO3TMTx5sUStMoJ02uJCnzio/2lV1jfS3NYxtYsHwrs5ZsKMlom3TRQObfOrqmeEiu1c5EZDKwGpgSb+p9D/C2Uuom235XAlcCjBw5ctKOHTtyuq4TU+avSBKSr4+fRkD8fb62Qw2ctuXBlO0CLLxoQkEmohldYXW+1YaCrvZOTWVT6vOh1MdXiYjIOqXU5EyOyYfjtQ1oU0r9Nf76EeCj9p2UUvcppSYrpSY3NDTk4bKp2Je7Vtt6OtyctQoKNmEzcWhpKp9Sj7Yp9fFpDHI21yildovILhEZo5TaCpwJ/CP3oWWOfWno1FjEjT09zg+epgKbSvw6tDSVT6lH25T6+DQG+Yqu+QawOB5Z8zpweZ7OmxUHuw07vGlbv2fED10jaQAI1rGr8SZq/xVMWXpqp5CmEDjZtkvdR1Pq49MY5CVOXim1IW6K+YhSqkUp9VY+zpsN9iXk0kgz4aiHeSgehXHS6VcV3lSia61ocM90bR7bUNLRNpUYDVSJlH3tGjtOS0VHs02wLiXErqCmkgqptaLJHTfb9sotHdx+wfiSrWXjFA1USuPTGFSckHdaQi6NNDOorpp5I39VsEJk1uX3UbUhRCDSGT18I+zwqLWihXyfwsu2Xeo+mlIfn6YChfzsqWMcw7omfOLrMPG/XY9rXR9mw/M/5stH3k9jaC/vhxqN0sRZCFx7aJm1OJO5FD//33bi6CYoo1ormvzQG7btcqhmqSkMZV+gzI5TGOKFk5pYsHyra9Gv1vVhXnj6HuYMuitR26auO5y+IqWLTd1p+W2lK9rDm90ufoIyqrWi8aZ1fZgp81ekLTaXb9t2uVSz1BSGnJOhsmHy5Mlq7dq1BbmWW8LGhZOaWLmlI1Ez5tkxlznWu6HuWKNdoB27TR0Sdv7R/1uftlHO+fUruWf0/zge3xvmGq3ZFZZME4Xy+f+xJwWaNNXX8uLcM7I6Z9b47Neg56c/skmGqjhzjR03p9bi1TsTgrhHKddkKFfziUf98sb6B9IWZVrLOXDyxII0LLELHFOzg8IlevU1Mm0Ik0/bdsnEr/sMLmhdH2b2IxuJ9hh3ZDjSxexHNgJ6fuaDihHybpqA28S2a9pulSupG+l8bo/65U5+ASuJpfjoM/wJdbs21HgOtD/p++GgO1AVnmIK2pKJX/dq5GKZr7cs25wQ8CbRHsUtyzbr+ZkHKsIm72WD9DuxnSpXdksta+qvdTx3Z1Wj84nqRib8AvW1ocRmsxBfxvH3pjZk6V7Fq/+T/PqlKz19ByWj2fUhMumelG96w8bvx7eQgs9GLm91Rh13c9uuyYyKEPJemqrThHeKalkaaWZu21W0HWogpoTOqiaqTvk/rl59guO579w907ChW7HVLzczb8GoumfeaBlpJ07akJ00be6KKXD6KsVMFMpnDaScnLgV0MilEqgIIZ8uztg+4WecMtJV0F+0+zcEZsSo+3wbjJ7heu6H2qd41i/PW/EmvyGVHvvpzMTCU+xic9Y+BS/OPSPr6+Y0j0+8La0iBCSteP1s12RGRdjk09kgnZxav1ztLBTtQt3t3AK0Rk6nxSnyxuE86bYDzpEIdSPjppk0eGhHOjOxOFRColC6eewZFWPa3dMEF8ybPo7ZD29M6gUbCgjzpo9zHZeOxvFPRQh5twQoL021yadzavbUMcxasiHFUasgxXFpnXhOrfyczp/ALRJh9KWw7efeJhsfbe4qQeBoCo+XAuUrasulkYtdSF908ohESHM6oa2jxTKjIsw12SyN/ZowWiY2uca8W7Ucu+3SScB7PnjcIhHan0w1C33gq2Xd5k5TJLIoiOd1n7iZcm5ZttnznE52/kfXhZk9dUyKecnJ6avr2GdGRWjykLmmmokJw4/W75blGhQhplT6JaVXJEKZt7XTlABZFsTzuk9mLdngeMxbnVFa14dd57rfkF43jd0tNFlHizlTMUI+G/w+GJrHNiQlT0GqVu42wXqU8uxlaWomS4a5x+lrNLnSueY66nzErJvzMRzpIhg3ObrNXzdTDqSaMq349Ve5PQyCmZpC+zgVYa7pTVrXh3l0XThJwAtw4aTkB4TXBAtHupi1ZAM3tm5KObe5bHWK0/dja9doTNzi2VvXh+kXbXc+yLKCtM5HOGxydAub9PJ5eWnVfkN6vRQnHS3mn7wJeREJish6Efl9vs5ZCtyybHOKNqGAlVs6kmycTx03k88Oftb1PApYvHpn0o1i1VTscfra1q7JBK949gXLt7r3O7asFL0K6znZvFsmNrmGOXopPX79YW7nGFgXoqYqkPRa90J2J5+a/LeAV/J4vqLTuj7smnU3mSeTMlHrusPMH34vlzW+6Ho+MyLHxK6pLI00c9qWBzl+0zKjKJoW8BqfeNm5211Wip2xmqSVYjqbttP786aPy1ir9hso4fQwCAWFd9/vTirf/e773dyybHPmGbl9hLzY5EVkOHAucBtwTT7OWQp4eeuvb/xFSjRMlepi3shf8VTnWa62ynCkiynzV2QXZqnRuOAmoMORLprqaxP9jucMW0RjaC/t0SHc//aXmWdRJLxs7Ob7drLNwfDjD3M691vvHaQzGkvaLxpTCWVMh1Omki/H693AHGCA2w4iciVwJcDIkeXhTPTSbI6ucnCSQqJAmVNsPRj2fLvN04pVC9IJHxq/eCXtNY9t4NF1YZZGmhPC3ix7bMWrsJ6Xdm4V2OacnbVkQ17mrP3cV7tE9FjRxfeSyVnIi8h5wB6l1DoROd1tP6XUfcB9YNSTz/W6hcDtxqmvDSFumajxAmVrd+xPicgRUqtfgnOYpd+ED/0g0IB30p7fXrFWzdlPdI2dfCQpOc1n65j8osMpD5Nz0xARuR24BOgG+gFHAo8ppb7kdkwhm4bkgmfjh/pVrk1DTFu6fcK6TVIBts0/N2mbn8YPmTam0FQ2o+Y+4bjdaX71Brk2K3Gaz6GAgJBSijgdRWmQUgCK0jREKXU9cH18AKcD13oJ+HLC296Yvi6H3e7odhM42Tr9xBLrOvEaK35LdfQWuZa0dprP1no2ftHhlMn06WQoP3g6iDLMRM2kxo6fxg+6TrzGSjY1nDLFyzyYa7OSbOZtbSjIR0cexerX36JHKYIiKTksfZ28JkMppVYppc7L5znLEbeklExq7PiJJdZ14jVWeru8cbra8rmWtM503poC/W87DySCGHqU4tF1YR1GaaHiG3n7Ja0DM4OGxHZtSoAZp4zk1pbxKfvnMiZtk++7FMPh7tdPlO243GzybiYbwX31oG3ylmO0kPchLO3FnTASSe7cfw0TPvF1X3Z3E7+RCpmMPdtoCE15UqyH++i5T7iGBWfi2PV6ENzYuolf/3VXwvRyynEDefG1/Y7nqa8NcaArmpcxlQvZCHlduwYf3W8cygDXBQ7y5SPvT6npkc6umFH7NB+0TGxKLJPT1RrRVAbFKrWbD/Ogl8nHrBNlNb382UXAA7x3qJv6uszLKvQ1tJDHhwPTpQxwY2hv0s3Vuj5MQJwaCyaT7xtS19fuWxTL4Z6PNpJuc/Xbv93I1Us2ONaJciPaozgY7Ulp5amja5LRQh4fGopLuV+z6FN7pCuhoThlsToRjh+TD3SUTd+iWA53P45dt6ADE6/KkpkwvX4lL4y9nJfHnsvzYy9nev1KAM6vX8mfT7iClldG+G6MUunoEEq8Q89a14fZsPNi5gy6i7rAwcT7nbEa7tw9EzBuLq8Kfm7kK4M119A1TXlRiFBJN7xCiv1kvKarj+OH6fUrmT/83sT9OLy6g/nD72VS3St8btAz1Kn4feqzMUqlozV53DUUMATxQ+1TksoAtx1qYG7bVSyNNCduLi+t2c2AYzeptK4PM/uRjUn2ytmPbEyr8edjGa0pH3o7VDJb/JgNneZqpswZtihJ4QLDRzZj8B9Tticao/Rh+rQmn05rnjJ/RVK9d7O4k1lrxhrF4lVb46jaUFJpVCtmVcrZU8dwy7LNKenb0R7FLcs2e97A2VYC1JQvpdiY3Y/Z0D5X3SqxetEY2uu4PUjMcbtra80+Qp8V8n6Wlm6TNqZUSniWVwW/SFfUtTiZ9dpu5h63mvZWSvGm1/Qt/JoN7ZUl3Sq2utEedW6V6Rrz0MdbaPZZc42fpWUmDi7rEtoJhbvZxry2RlPOZGM2bJnYxKnHD3J8L+Bywzi2ynRDt9Dsu0Lez9Iy00nbMrGJF+ee4SrMFbg+BLxwa7Gm0ZQS2foKtu9zWzE772+2ynSz8igFMSXs7h4KJ99Ha+R0z4ifSqfPmmv8LC2ztXWnS7U+/vonHe2QIlAlyWncoYAwb/o4359LoykmTmbDdL6vbEJ9l0aaubnxPgZXvZPyXjjawGlbHkSAheMn5Fzjvtzps5q8Xy3d1M63zT+XF+ee4WtipDu3m6NJKVjwuROTNKEFnzvxcGmFeNNwHf+rySsZzK10cfBO+3sVNYPsQn2n169kQCD14XAwFkwb2tzXEgX7rCbfmxEp6c7tVve7qb7W2YFqr52j4381+SKDuZVN5yc/PQ9GDXa+H4IB4dwjVyT1pb1z90yWRpqZM2wR1YHulGPei9UlhTbPcmkX2JcSBfuskIfejUjxOnfGySwOtXMS8b9ayGtyIYO5lU2TmnS+r9b1Ydf6NJ8f/BzfPfpeai1JT3cMv5dTB2ylKeTcY7k++K6v0Oa+lCjYp4V8LuRSUtXcb97SzYn4+X4hD8uZW5xvH4//1eSBDOZWNuUz0vm+FizfmhQ+Ob1+ZUJzjyFUSXLse23gIBcNfMI1uCHQf2RSieFiZgeXCvlo5D0CWAQcjRFAcp9S6p5cz1tK2AV689gGlqzZlUhcMjNTITNnzsHuwxP4rc4os5ZsYO2O/al15z2ahms0OZHB3MqmfEY6Idse6bIIdkM7N0MnAy7R8+KyvVtquXXnxTw094mkstsXTmpi5ZaOPpsomI9G3scAxyil/iYiA4B1QItS6h9ux5RaPXkvnGp3uzGwLsT6757t67xudecFWHjRhNSGJWmahms0WZHB3Mq2jr1VSTqqNoQIfLL6T1zf+AuGBvcA7jHxflHA7PBsHtn3qZT3KqmRTkk0DRGRx4F7lVJPue1TTkI+XRMQO9t9Nipwa8BgMrAuRKQzeljzqF/lqzOVJjOK0WGp5PDZ9Qxy7/w0++GN/PuRK5IKjGVCTDk/ENqjQzn1lQdcj6uUTlHZCPm82uRFZBQwEfirw3tXAlcCjBxZPmaG3vLCp6vGZ5YySEQwXHA6LS3be2UsfZVsokUqEq+G9JYHQGdVIxt2z6Q9MiUrAW+WL3AqMOZFtwoQQNEeHcIzb59kVJq0HN8ttcx/4xLPc/SlaBo7eYuTF5EjgEeBq5VSb9vfV0rdp5SarJSa3NDQkK/L9jqZeOEzyUydPXWMZ5kDK30trrdQ6BjqNJimnM4dgKKuO8ycQXcxrX5lRt3HzIepuXJ1KzDmRGeshmt2zeK4Tcs4bcuD3Nz+taSKsNQdy60dVyeKB7rRl6Jp7ORFyItICEPAL1ZKPZaPc5YKTolNoYCkLBkzzUxtmdjEjFOcVzRmQ4TXx0/jhXhDhL6sifQWutlKGlzaXs4Ztgg43NEpXWKU/WFqNttxI6aMxEBrSW8rSyPNnLblQT6x62lo2c7P26d4nq+vRdPYyUd0jQA/A15RSt2V+5BKC7fEJqdtmS7xb20ZT8/ri/n6oAcSyR725ajZEGFQXTVQeY2Ji4lutpIGj7aXJva+wmDcM1bbvd33dOfumSk2edM1GI42JBKe0vHewW5a14c9TZ9BkURviCnzV/RJ30s+omtOA54HNkGioPN3lFJPuh1TTo7XXmXbYrpX/xdV6vAEdXMsdVY1Uff5tgIOrvLJNlqkz9A6yjG8su2QURvGCTMRKV1E2i2NP+GSwU8mzfXOWI2j5u6FAB8Y2p9/7XnP9f2FF02omP9zNo7XnM01SqkXlFKilPqIUmpC/MdVwGssbLwhScCDeyhZXXffqpxXCEq1w1LJcOJtRjilBWvbSyfaI12+WmGeeeSalLleFzjIdcN+kVGlVgWuAh50/RrohRBKP2hNPs6vAnj3o7ci8PFfOMYu9/kQQE3vEY+uib23M6l2jBtN9bVJJhprBqt5/POHzmbd8Z8hIKlzP6aEwIwYa1bdS+P273GMrWZNptx90QTPpiTT61dy3bBFHBPay56eBnaNvImTTr8q4+sUiqKHUGoyxC3b0BFF55rrqLMIeR0CqLGT94d+PLzy+DR5HWCYRsKRrkS2qVvD7e++EWB3dwONoT0p59jT08CwbYs56c3roLoz6TggK0HvZrO3j29Y1R6OaruWNasoaUGfKX221HBJ4LAcTnltoV+0PSmCoa8vQzXJ+Cnrmy1+nNHmQ8B0xro13L664SHueOOSlO5OSsHA0CFY9y3PqB6/TK9fySmbPsoLI8/ixbFXML1+ZdL7TuOrDRxkxM7vZ3SdUkcL+WIyeoaRPl53LGDE/B5+nUp7dEiSANchgBorvfnQdwolBiMze2Cdc36IWzx8Y2gvj8e7O+3vHpCIrBGBGhVBHdyX0fmcMLX0YVV7EBRN1Xu4Y/i9nG8R9G7nGxp0rnBZrmghnyGZNk1Iy+gZ0LIdLo4Zv0fPgBNvS9FyTIeXVYBn0oNWU/n05kPfyUl990UTWP/ds4m4NJp3i4c3ty+NNNMZ65fSgNutIXe6+Horblr6PWMfSTh23c63p6d8kjX9oIV8BvTmcjiJ0TO4c/81iaw+a1KIVYBn0zhZU7n09kPfrUua2/l/vP8KulyUlcTYXLRpezyIeVydV0luC65af+cO7j7lH9SGgo4NwbtiNewaeZOva5QLWshnQCFt4BM+8XU+/dqiRDq3tduNiZN2deEko1FCX21a3JfJ9aGf7Sq1eWxDSomOUFB4ZO8nuc5SgsApg9VNm97fM8DxuM5ozHF/k6AIgrc2PvGNObQMXJVoCG5eZ3f3UF4e/oOKcrpCHw+hzDQSwa1y5PT6lcw/7jfUdbfntUJkpuPTyT2abKNrsp07N7ZuYvHqnZ4hk14RMfYIF/CfFOV0rWWRZrbNP9e5hLIFa0JXOd0jJVFq2A+lIOSzmdROZYedJmm31HJrx9X8vD3zan254FYWuVLKrGp6j2zmjrWyJGQmsAWyfjC4XSsWP2Gg/7Fw4m2s2fEWk8PfcLTxx5Rw3KZlrp+zVPNPipLxWq5kY3pxWg47OXiqVBdfPvL+3rXbO6CjbTTZ4jRHptevZMmwLxhJe62jDO3Ygr11n1vIpDX0sam+lu3zz2XhRRMSDlCz4JjVNJkOp2sFJJ4x3rkD9ZdL+OfLfyQcdTbb2M1E1s9fMN9bgeizQj4bgWjawK14hYmZZGq3z9Y26uYAU6Dt8xpP7HPH1JSHV3cAykjae+nKJEFvv1fS3QtW/4DpxL37ogmu0TSe400TTikovjjwCZ55+yTXSLWk81k+f6Xln/RZIZ9tJELLxKak2hrpwsQSrx0eHnZhfmPrJibc8ieuXrIhKy3CLZaZDM+j6XvY545jY4+eTqP8cBz7veJ1LzjVBUrUmfdhMbY/B/yEUwbEqJEzN43z1+6crrQVcZ8V8rlEIliPdQrDSqcpgPOS8JerdxLpSo059qtFWKNtnChnbUTTu9gjtdxDEA+XH7bfK8+8fVLCLp4gWMfwTy1MhFy2rg8z4ZY/MWruE1y9ZIOv3smQWuHJ6b5zojG019McNLAulPLwqbT8kz4r5HOpQGg9dlmkmTv3X0N7dKhvTQGcl4Re+NUizGWwG5n0q9X0Laxx8IH+Li0660Ym7W92Q5tev5LPDXrGVllSYPSliUgzs8erkyKTKWb44+7uoZ77pdP433cIyay0/JM+XaCsZWJT1h7z5GPPpXX91x1raNfXhpg3fVzKdTJd+mWqRZhFopy2azRpOfG21BDEYJ2x3cK86eO4/rFNLn1bFbQfrjq+YPlWoimqfvY89d5ZnDF+lnFvvfQ1ePWnWHX+bqnlB7sv9TyHubq13p/m3/OWbk48kPr5TMIqRfq0kE8ig471Trh1kHJ7iKRr5G2lNhSkeWxDRp1tnAS813aNJglz7qe5JxJZr/9Ib97Jl01bIPUeOPkn0DCFzjXX0S/aTnt0CPe//WU6Bk9DIvs9K2jax9W6Ppwk4AHe6oyWbYXXvAh5EfkMcA8QBO5XSs3Px3kLhj1xwowkAP+CfttiWnbcQMvInTDWvCHOcI239dM9Bwyb4bkfOYZH14UzKinc5PIQyaQhQylRqnHLFU28zLAX5v9l8rAh8UgcGxbzTiaKjRuecfuR05m94T7bamF/2nMeVXu4wJpT/oyJk9ZfDuSjx2sQ+DHwaaANWCMiS5VS/8j13AXDoWFxIpLAY5InJjhPMn/4vdSay9X4Q2LNjre4/pkPeApnq+BqHtvAyi0dKYJsyvwVriFdbhPO6SFSrnZFXTe/NLH+X+4ktW+r3bwze+oYZj+8kX8/ckUi+emtniOokSj9A+8DRjmDW9qvdIyVF1KdvVbmLd2clTnIasFM5ysrR59WPjT5k4FXlVKvA4jIb4DzgfIR8i4Ni123kzzBl4xddFjAm/R0MmLn9+mKPpC02Sqc/foEso3ph9ybjZcCXnHL5fh5KgXr/8UUyqbwDvQfeVjAt46Czp201I1k/GmfonH/ksT9MrjqnaRzDq56hwXD7046JxgCfsYpIz3/39k6dK1VNNMJ8XL0aeVDyDcBuyyv24CP2XcSkSuBKwFGjnTx3BcLtw5Nde7jtE7wTOtSZ+N0dZp8lza+CK1fd7WZ5uJYLiUqLW65EmhdH06Zk0sjzSyNNCPgXD+mcwfHd/4CAt7adk2ghznDFiWE/GWNLxqO3c52aM1fbSgTa1CDW8CCSTn6tArmMlZK3aeUmqyUmtzQUGL1mt06NNkiCaxYBUymdakzjZRxCun67OBnubHh7vjDyTkjsVKotLjlcsdcxbqR+L84mUF99jQ2FafLGl9k3tH3xBvZp85zM6Fw1NwnMv0YQKoJM50Qr68N5befRAHIh5APAyMsr4fHt5UPbh2aPLQFq4BxTMwI1rFr5E15ibd1iun/3rG/okrZNNmeTjrXXJfRucuBSotbLjfsmdnzlm6mK9rD9PqVvDD2cl4fP40Xxl7O9PqVyf8XD3NnOtqjQ6gNBY26Ny7+MmtCoR+a6mv50ikjXXNjbmx1f3ABhALCe4e6y66mTT7MNWuAD4rIaAzh/gXg4jyct7D4iCSwYnVsmsvK64b9gsbqDiRuOjlp9AxuP6p3okJqu9sdt5t9YN2uUY5RKpXkXyg3nJzekFoF0my2fcnwUYQZz5T5K1jiEnGjEMRDoz8YC/KD3Zdy4aQm+r3XnlrTAKBzp6+EQrc8FafPuXi1+0OpvjaEiBFKaaUcfEM5C3mlVLeIXAUsxwihfEAptTnnkZU4dsGzjnNYc+KslH92PuziTjda+7AGmqpTu92bfWCdrlnOUSqV4l8oN9wEqVvFyXF75jMzHlHmFHHTGavhd5Gz+I+hG6jrbmdft1t0zemweidfGesempnOJyPAhpvP9v05vQw1G24+m9EuJqFS9w3lJU5eKfUk8GTaHSuMQgkepxvtjt2XON5A9j6w6c7jqonkmBymqQzc5pJbsEG/aLtrxI21VvxP3jXi3ae71LE3uXO384Pi+699gfq6UIpmnTRGPz6b+Dx/fsRO2oc517I3c0vcAiBK3TekM17LAKcbzesGckt4crthUyZuPpLDNBWBm2Db3d1AY8h5JWnFjLixE450JRqVpGsgMrftKod5/gkCRAkFhWhPqg7uy2djmecBOWxyMsdtYp7HKfdEMNofljJ9tjNUOeHWtad/dZDOQz1Jy0yv7lZu5xFg4UUTDh/TOsolpPRYaNme1WfQFJd8twVcdOarTHxjTpLzvzNWw/VtV/G4j6YfVsFuvp5Wv5I7rEmFpG8FWF8bon9NFeFIVyL8scnh8zl+/h1THOe5tTWgOTbzmLU79ie1OzS/j0K1D9SdoSoUp+iSUFA41B1LuVEunORuQpo9dYyjD0tBcgniLJLDNKVLLp2O3JrFX736BK7Z8bWUOu2Px+PkvbALeOKvv9P4i5SkQntnKTtmApQAw47qx90XTUiUNU73+ZXLfLaboqzHPPH3N1LGXuolvLW5psj40bCcokveO9idkuGngJVbnBOwzPNcvWSD43vm8rk90sVfTmhgWFXqUtwrOSwflGPkTzmQa8aw1fdk1ezDOJtiFMbDwK1ch5sN3i15MF0XKPN8boEEbp//zW7nee6W99IV7XGN5ill56sW8kXEKdrlhafv4ex//Yq67vYkh6fdyZutp98ro8+8Wf67/ZKUZbNTclg+hXI5R/70Gnlyfvv2xfjAT9hiusbxbmbDPT3OQnd3t3+bt9PDy+1z3t5+SXLNKZwb/vihlJ2v2lxTROw3zPT6lXzvmB+5ZvdZybafq5+07KWRZq5LNGRwTg7Ld7PjSuuraSWrnr2mUzAPGc1uc0XiY8uEdEqEH4enW3LbrpE3OWaeN37yLgbWhfCLVai3rg+7mo/Wcg4vD/8B4UOpDX8yKVFT6ol5WsgXEfsN46evpkm2/Vz9lhpeGmnm4/94AC6O0Xrsi0xZckySkMq3UK7U+jRZPwy9KqNmiG9fjA+8NNZ03dXMh92sJRvoFwoYCUaW4046/SrXzPObp41zne9OTPze4RaDTmqNWdHypNOvYs2Jf+Ojr/4x0RpwYF2IGR8b6et6QRFPP1gpoIV8EbHfMH76app49XOdXr+Sp46fyfR/jDAiZSzan9PDwU1paayvdRVSbkvgbIVypdanyfphmEfnd8vEJtdEn0z/X25auJPD04p9Hr3VGeVgd4yFluNa14cNZWL1j5my82laj30xsXq0O4DTVYP0ip8H4wE3a8kGpsxfwdod+znYfbgN4FudUR5dF+bCSU2uSpFZ0uFfHz6Pr0bOYs2qexOfs9Rq22ghX0TsN4xrP0oXh6fZk9M63c108+HVHQQkdZnvFC1xx6RNvDj2CscaJG5Cyu0my1YoV2p9mqxXKG5O7iyd327CKtP/V7a9kdM97PyseKw9aH/4+RMz0uxNrPV2nh97OZN4kl+u3uk4tpVbOpg9dUzKXLffY03Ve/hw27WsWXVvymeYtWRD2po4vY12vBYQJ0fl7ReMT2y7/+0vc2PD3cmFxywOTzdHpzViwcvk0xo5Pen4hRdNoKV+Fbz0fag2TAPDqzu4Z8QPORT4GTX19zIrUu/4WXqUojYUzFtTkkqtT5N1lqTPHqt+yWcTmWwyvdM97DKJAHIzF6bDrd4O4Jqwdf1jm1L8WE73WG3goGP/CAUsXr2TyccOKtpc1slQBcItqSRFC0pEVOwACYLqgbpjWVN/baImiP14IHHu18dPMzR4GwrhhFeeTDl+w4cvpqbnLedBB+uY9+a3eKh9SspbZsJJpQnlfOP7/+5EnktLrFl1LyN2fp+hwQ729DSwa+RNhh28ALhF1JiROKPnPuFqO982/9zEa6/2fOl4YezljnVw7MlP6XC7x2JKOG7TMsdj0kUc+SWbZCityRcI35qKeRPbygp8+N1r+XT/5Mw/83hz8nz7txtpjzoXdHqzuyHl+p/u/zTV3W+5G+V7OpkzbBFLOj7pqAHqomHpyWmFkmFlVE+2LeakN6+DKmNODavaw7A3r4NtAwtSqmL21DG88PQ9XN3wUKI8wd0dl3HaWd8C/K94stHgE+dy8Xmli8O3E+k5gkG2jlYAMYTXx09LKjFiEo50eVaH7U20Tb5A+LHNmk6btmdnpURW1Lpk/pnHt0xsIqaUY237zlgNt7dfknLsnGGL0oaK1XW3Z2WD1RzGakv2ck72KnmM1smGlvpVSXZs01TSUr8K8O+TySXays3n5eoLi2O14//thC8yIPCu435VEkv6bNPrVyad46SNH0X9KpASDNHbaE2+QLhVzDM1FesytHGEf43Dquk01tc6Fi67/+0vs5azgOQbxJcGUzdSa+yVQLFLVWy8IaXJTZXqMh4y8WQ/SF3xAIlM7Mb6WvqFAnRFYymn94NTRctDsSqOCr7NtvHnAfBWzwDmWRqJn1+/ktstxzhp8E6Y5RiWRppTfAF07qDrxf/k5R1vFcRcpoV8AWhdH+bd97tTtoeCkpjI1mWom8nFrnHYNR3TuWat/GfafydAii3zjahzTXqTrlgNL9dfy0m+P2kyukxBCZFFH+O84uMhY1cm3BqWZMvSSDOnHT+Ez4d+iurcSaR7AEcE3mFA8PC9OajqHX4w/B4A/vD2Gcw79lfUqYNup/TEVKLcHLWN279P6/r/6PV7QptrCsCC5VuJxlIdNf2rqxL/YOsy1M3kYk+3rqlK/vd5hbc5vdc+KjXDMKZAxX+6YtUsWbMzq1hfp5C4QoSTlWKcckmQRR/jvJJFSGim9nczTDTgYYK857WToWU7cnGMgUcOJOTQVLw60G2YRgXq1W7f17djKmVuK+ZjQh0FyejOSZMXkQXANOAQ8BpwuVIqkodxVRRudsQDlgJjVseTV614K5GuaEp9Fy/TSup7ZxiOt3gER6TnCOqki+qAodkMqnqH7x3zI+Yuj7Fg+TlJmng6Ld3pBlXAL+Mt1m5tGe/6fWVL6/owsx/emHighiNdzH54Y+Kz2/ftU6sM07nay41gXL/XLEJCM7G/26NX3Bp7J53Tw1TVGNpLtEe5FjGzcihWhUJREzg8361KmdfKvBAZ3blq8k8BH1ZKfQT4J3B97kOqPPxkc9odT0sjzZy25UGO27QskW7thFv2pG+NdvQMo0b8xTHe7emXEPAmpm3RmpzipKVfvWQDE7/3p8R1vCbv4tXZrQ7SMW/p5pQVUzSmmLc0uRtlvuvulA2W/zUt23tFwLt+r6NnuJYscMOr5o4VJwetr+Qvj1WEqYXf3n5JygroUKyKfd0DEvVurm37FrPbrk6UXd7dPTSpBr7XyrwQGd05afJKqT9ZXq4GPpvbcCoTP4koVseTvVtOOuwROvOWbk4qQ+y3omO6EDPrA8VpGf1W5+GVhVdJWbNmSr41Z3vpZbftuZbe1TiT9nvNMCTU7b65cFJTUulip1WYr+SvE2+D1ZeDSp4fh2JVCS18LefAyRMTK6DOqka+u+NiHtn3qZTxPvXeWdx+wXhm2erluK3Mjf17P6M7n47XK4Albm+KyJXAlQAjRxbI2VMi+I2Vttft3vD8j/nykffTGNrL7u4G5r9xiaNG7xShY8ePEHs/1BivgJmM1eGbbnlpXmf21DEpkz3pnB7n6W1TSqUWQys2+f5ec8kx8HWs+cBZ+y1UdB+o1OiazkPdtEZOpyXeEa0OOG19mL/ElbGW+lVcO+znNIb2ciAwjHuevwxFavLg0kgzzx08m/41VYnx3H5BYUyEaYW8iDwNDHN46wal1OPxfW4AugHX4E+l1H3AfWBkvGY12jLGbxiiKeAm82RSrevG0B7uOvYnVAcDSVqEVTtJ56hKl5BRd9IddK/+r5SWblaHr/lA8Yp0CEe6WLB8q+dKxG2Zmktd+YEuYar2MrXl2pC51OmN7zWX8F1fx8ZXF2bZZfsq2Lo6tfq9ADY8/2PmDPp/iciZgeoN5gy6i/2dh1KUsdpQkHnTx5VmMpRS6iyl1IcdfkwBfxlwHjBDFaNGQgVhtWnOHrYopRValerihuG/TBRMspc59RNi5ml7Hj2DqlP+D+qORSGEDw1N6a/51nsHiXQe8ryGpBmLV82UXEoY3zxtHKFgssU2FBRunjYuaVulFkMrNuX+vbZMbKJ/Tarea59/5n365SPvTwmNdGpXWOxyxLlG13wGmAN8SinVmW5/jTdWAedmHz8qtjtRMKlHKR5dF04UP/Lq+mTSFe3h6iUbuHrJBoIiXPPBl7iw6n+S65m0bEeANevDPLdjM3BYs+lMk4ji5EuYXr+S64Yt4pjQ3sPXmPgZx+NzWfJnYhbzs5/GJ/EaOy2dOzl7fCN37p7Jz9unlOX36mf+mfep3zIJ9vu00ORqk78XqAGeEkO7XK2U+krOo+qjWCeS34Qoq63dKuCn169MG4J57lEruLz6cCbesKo9HNV2LWtWwUmnX0XLxCYWLN/q6tA0CYoQU8pxuW7P9ktXMyXXJb/f5b3O4s0TZgereGhkXXeYeUffw7xp42D0uWkOLj38zD/zPvV7j0JxHfs5hVAqpT6glBqhlJoQ/9ECPgesE8lvQhQcnnSmGcde79qplgZ4ZeJ9LxF66ccE1KNUovmDPXQtk25XUP5L/j5HkWvi5Bs/88+8TzO5R+GwT6zQ6IzXEsI6wZZGmpnbdlUi9tbaf9KOOelMTd5JsDrZCt0z8fYm4pz9tro0bf32m8TtGqpzp2Mcf7ZNKTRFotg1cfKMn/lnznH7PfqWHMOd+69xzWmBND6xXkLXrikh7LZiaw0aL5rHGt3sm+JLTb+2Qj/LTYWznd1OV7SHeUs307+mKtE5qkcp9vQ4Zwy2H2pIrBLsETTalFJGFLsmTi+Qbv5Z79NlkWbWcTgbfB4wIcdQ5nyjNfkSwyxLu/CiCWn7WJqs3GIIalPD8FtS1Wm5GVPwzNvJJckUJDSb+trkcEQrka4o4UgX0+tX8uyYy3h9/DQGhg5BoDppv65YDXfsTi59bD4kNGVGsWviFAmv8tHmasCNQudjaCFfgpghWukiZUzCkS6mzF/BrCUb6BcK8MM3nW2Ff+5/DV86ZWTi4fHEgTN4oftcrJUAAgKfG/RMkv3erAuy8KIJjiFmVuz+gBoVMaqdVQ/GTGd3MztFuqKVX1qg0siiXEElY5YTmRWPXnOi0PkYuv1fCeLX4emGANMs0TVvRIfwYv9r+Pzn5qbu3DrKcblttkSzthic/chGoj3KM3LHrcUadcca9VLSfL58tUnTaAqNn9aEvls/uqDb/1UIuS7nFCTs+eak+rzbpHJxkDWG9ib6uLbUr+Kt57/KP0843DDBVFLszZBdG5FYrjN76hiuXrLBcbdca4ZrNMXCLePcGmJcjLwBba4pMH6qQ2aznLO2KHth7OUJc4tTtp71+p1VjY7n29PTkBDwrL6cgVXvIELix4o1csfNH2C9TsvEppRSAyZmerlGU264KWcxpYra+lEL+QLit8StU6yuF+ni4s3J53T97+64mEP0SzpfZ6yG/26/hNkPb6RzzXUpVfqcMDV4v7HDN08b5xieaVaodEI3BNGUMn5KihcDLeQLiFddFqsAW7B8KxdOanJ13ARFEjG8taFA2rh4c5I5Xf+RfZ9i9s6vO8bjR2OKftF2X59NULww9nIAx/j+n7cnV+ZrmdiUUYXKPlsDXlM2lGoin7bJFxC35ZwpsKyVFx9dF3aNrulRymjfFxd4XnHxAolJ5nb9xyPNPG6LdjGdq+Kzqr3IYfv83LarOG3Lg0nvB0UYPfeJJLtkUwYlDHQNeE2pU6o1kbSQLyBudTGCIo4CzK3gmL3Ko1dS06nHHy6K5NXIw0pKd3kHzGG52eftIZLm5whHupj9iNGSz1djB/Oz6BrwmjKgFBP5tLmmgLgt57w0dvv+TtmnXnbwv+08kDBpOF3fySDkWG+Gww2+93cP4Fu7vu04Zji8ggDnpsrRHsUtyzZnVMKgVO2dGk2po4V8AXETam79KK3vm/s7lfGdM2wR/eQg3SqAUiTZ1a3RNU7XP/X4QSnXda03gzB60+/56D9+HX/tTHt0SCJLNuayk9ncwytz0Eqp2js1mlJHm2sKjNtyzs1sYd/fmkhkN6sEiBFT0BjqSDhdl0aaCUe6GDX3CYIifPFjI5KSjabMX5EylnQ1babXr+SuEQsdtfSYIhFJM4knWTLWu9yxX0rV3qnRlDpayJcAmQgwqx3byawScElSAsP888vVRlLSrS1GFquTTfvO3TNTbPKHYlXUBt7n9fHnJV3HiaWR5pQHkH08XjVw3ChFe2dv0tu9bjV9Ay3kS4RMml2A8UBwzS6N4+YE/fVfdyWEvJMz1t5dPtJzBP0DnQyueod0tEcbEse6hXX+4e0zmDd9nNPhmji59LrVaKzkxSYvIt8WESUizumOmrxi2rHfDzlnq1pxehBYHb1uiVdLI82ctuVBjtu0jM5YP2oC7vU4TKxJT15hnQs+d6IWVGnIpdetRmMlZyEvIiOAs4Hy7BJQRtgzPjcPnUu3eEeXuJUZsDbp8Eq8AneBbaVbBZKqS7pdd09PQ9pzaXTIqCZ/5EOTX4jRzLvw5Sz7EE4ZnzOf+QD/2zWHtkMNKEVKJItXKzJTI2xdH/ZMvAJ3gW29zjW7ZiWZhdzCOpdHJnPSxo+ifhUwKmBuW+x57r6KDhnV5IuchLyInA+ElVIbfex7pYisFZG1HR0OpWg1nrgu3/95EqdteZDRm37P1bu+7atdIByuQX/1kg2epVHBWWCbEfa7u4cyt+0qgKQCaZBa3uDh/WfyuUHP0FS9x8ik7dxhNIHWgj4FHTKqyRdp68mLyNPAMIe3bgC+A5ytlDogItuByUqptGt7XU8+c0bPfaJoS6W6UIB/XLHfaM7cudNo7XbibTB6Bq3rw6xYvjAlGqczVpPykHGrNd+tAgRFIZbzanR0jSaVbOrJZ900RETGA88AZqv24UA7cLJSarfXsVrIZ06ujURy5e6LJrgKmPafH01jKLWPq9l4xOT18dMISJr5Fqzr052FNBovCto0RCm1CRhqufh2fGrymsxxqvPip8F2vvj2bzeydsd+Vm7pSGiWzWMbWLmlg+dHOJvf/DYOT6Kn01gxaCGv0eQFHSdfotiX6s1jG+gXCiQJ+UKab6yJVIn2f+/tpX3YECI9RzDIIYbeqXH4HcPvpdaj8Bng2q1Ko9FkTt6EvFJqVL7O1SfZtjhh8+6sauSFHRcTjnwKMJykpoDNJ169Wr2OsWeyHowFORSrojrQndjPKbLn+UNns6xmCFPeu4tjQnuJIVRJLPUidSNz/3AajQbQjbxLg22LjSiTns7EJjMcsj3akFPNFzecygk7OUvtuDlP93UPoCvWL6cHBqBt8hqNB7qRdxnSuj7MKZuuYVhVZ9J2rxo0bgysC/F+NJY2JBK8yw54XcctMWpg8F0mxatT+sW8znXDFtFYvVdH12g0vYAuNVxEWteHmf3wRoYGvZ2R1lZ+btSGgtw8bZxn6WIrXmUHAOprQ47lDtwSo9IlTLlhlk+Qi2PQsl0LeI0mz2hNvojMW7qZaEz5ijppDLm/32SLoW6Z2JQScmm1v78RbfB0lgowb/o4mg78jhE7v8/QYEfCBONUodIrs9YPOotTo+k9tCZfRCJdRuMM54zSZGIe/yqnZhvWjEnT9j28uoOAKJqq93BEoJNDseRnvCmsFdBSv4qT3ryOYVV7CIhKMhs5NerO1megszg1mt5Fa/IlgLW0b1OoI6VvKkAQhygUDDu8iT3s8sJJTazc0uFof68O9BDpGcCeQ6nO0qb6Wtj49SRHMBw2G5225cG8OIKb4qGhC5ZvZdaSDTqrU6PpBbSQLyID60KJNnhLI80sjTS7Rq+Eo6nVG4MB4eZpRl12p/rjj64Lc/sF4xn+irP9/ajgu3z8n7917EilXtnp2P/Vy2yUCQFJTfDKpGa6TvnXaPyhzTVF5OZp4wgFk0Wpk+nGbOn3wtjLmV6/MrG9J6a4eskGpsxfwS3LNrvXH3eJO28/1MCFk5ocG2m/2e1eEtg6hmyJqexrpjtV5Lz+sU2J8skajeYwWpMvItYuT6aT1G66UaQPp/SqaROOdMGJt9H14n8mZZp2xmq4Y/clrNvdkdTz1eT29ktYOOKHKW3+AkLaMEs/iGRfM93r4aC1eY0mGa3JFxmzy5M17NEMKwxHG1KErJ9wSitBERg9w9VZ6iZQ13KOo7kG/DURSYdS7lE1ARFPrVw31NBo/KOFfIngVD88XSy7H8xmIGs5J9HOz+o4dRO0s6eOoT061PG990ONSSYeq/M3E9xaD/Yo5Wl+0Q01NBr/aCFfIrRMbEokMpnC062HayaJR+YKIdMmFC0Tm2gfdRNd9tDOYB11J93Bi3PPYNv8c3lx7hncPG2co7D2YmBdKPGZnVoPetnmdUMNjcY/2iZfQrRMbEq2KW+7I6WmTVc8lt1eZjgUFFAQtfQAtAo+q/3fb0RK+Kj/4Im/t/PlI++nMbSX90ON1J10R0pWqnmOb/92o2cbQetYzaiglolNzFqywfn6LuaXbD6LRtNX0QXKSh1LdUp7Rya7kIP8CT57SCbAZwc/y/eO/RV13e1JYzHx6l5VXxviQFfUcVxuDVEEWGhpVqLDJjV9nYJ2hsoFLeRLD1OAhiNdBEVSNHI/FSPdhHV9bYgNN5/tee1ZSzY4PiCa6mt5ce4Zjg+d2lAwEfKp0fQFdBVKTVbYBaiTycUpa9bexcmpe1VtKMh5Jx7DlPkrUjpKWTVyN1XDjJjRYZMaTXbkLORF5BvA14Ee4Aml1JycR6UpKE4C1I5rRI+li5OTrbx5bAOPrgsnZbVaG6CYiUzW7N+k68YdxzpsUqPJjpyEvIg0A+cDJyqlDoqIc8ydpqTxIyjf6jmCwQ5VK/d1H8H0+SsS9nG783jK/BVpHyBd0R5qqgLUhoKOJRbAEPZOpiAdNqnReJNrCOVXgflKqYMASqk9uQ9JU2j8CMqAS2qU4F1WwK+mfaArmhJCarW367BJjSY7cjXXfAj4hIjcBrwPXKuUWuO0o4hcCVwJMHKk7uFZSjjZ0u0cFUzV4gEGBt/hhbGXc+fumSxYXp1iH3fTwO001temhpBa0GGTGk12pBXyIvI0MMzhrRvixw8CTgFOAn4rIscph5AdpdR9wH1gRNfkMmhNfnGqoWPHrbGJyOGaOte3ASTXwfHzAPGrkXs9BDQajTM5hVCKyB+BO5RSK+OvXwNOUUp51qPVIZSljVO8u2MIpY3d3UMZNvPNlO32+Han6BotvDWa9BQjhLIVaAZWisiHgGog9+pVmqLiZGJZGmlmUF01X6n/GUcH9zg2Njm6yvnZrjVwjaZ45Op4fQA4TkReBn4DXOpkqtGUF25Ozgmf+Dqrx//NtXCZuNSt12g0xSMnTV4pdQj4Up7GoikR0jk51xy4iUFt1ybVpydYZ5Q5KAVcSkFoNH0RXdZAkx2lKki3LYa/XgGxQ4e3BarhYw+Uxvg0mhzQtWs0mkeGwKF9qdurB8NntbtIU95kI+R1PXlNZeEk4L22azQVjhbyGo1GU8FoIa+pLEKDM9uu0VQ4WshrKovJ94DYes5KyNiu0fRBtJDXVBajZ8ApD0LdsYAYv095UEfWaPosummIpvIYPUMLdY0mjtbkNRqNpoLRQl6j0WgqGC3kNRqNpoLRQl6j0WgqGC3kNRqNpoIpSu0aEekAdhT8wu4MoTzq4Otx5hc9zvyix5lfnMZ5rFKqIZOTFEXIlxoisjbToj/FQI8zv+hx5hc9zvySr3Fqc41Go9FUMFrIazQaTQWjhbzBfcUegE/0OPOLHmd+0ePML3kZp7bJazQaTQWjNXmNRqOpYLSQ12g0mgqmTwp5EVkiIhviP9tFZIPLfttFZFN8v4I3pRWReSIStoz1HJf9PiMiW0XkVRGZW4RxLhCRLSLydxH5nYjUu+xXlO8z3fcjIjXxOfGqiPxVREYVamyWMYwQkZUi8g8R2Swi33LY53QROWCZD98t9Djj4/D8P4rBj+Lf599F5KNFGOMYy/e0QUTeFpGrbfsU5fsUkQdEZI+IvGzZNkhEnhKRf8V/D3Q59tL4Pv8SkUt9XVAp1ad/gB8C33V5bzswpIhjmwdcm2afIPAacBxQDWwETijwOM8GquJ/3wHcUSrfp5/vB/ga8NP4318AlhThf30M8NH43wOAfzqM83Tg94UeW6b/R+Ac4A+AAKcAfy3yeIPAboxEoqJ/n8AngY8CL1u23QnMjf891+keAgYBr8d/D4z/PTDd9fqkJm8iIgJ8Hvh1sceSAycDryqlXldKHQJ+A5xfyAEopf6klOqOv1wNDC/k9dPg5/s5H/h5/O9HgDPjc6NgKKXeUEr9Lf73O8ArQFMhx5BHzgcWKYPVQL2IHFPE8ZwJvKaUKokse6XUc8B+22brHPw50OJw6FTgKaXUfqXUW8BTwGfSXa9PC3ngE8CbSql/ubyvgD+JyDoRubKA47JyVXzJ+4DLEq4J2GV53UZxhcMVGFqcE8X4Pv18P4l94g+rA0DRmsLGzUUTgb86vP1xEdkoIn8QkXGFHVmCdP/HUpuTX8BdkSuF7xPgaKXUG/G/dwNHO+yT1fdasZ2hRORpYJjDWzcopR6P//1FvLX405RSYREZCjwlIlviT+GCjBP4H+D7GDfV9zFMS1fk8/p+8fN9isgNQDew2OU0vf59ljsicgTwKHC1Uupt29t/wzA5vBv3z7QCHyzwEKGM/o8iUg1MB653eLtUvs8klFJKRPIW216xQl4pdZbX+yJSBVwATPI4Rzj+e4+I/A5j6Z/XyZxunCYi8n/A7x3eCgMjLK+Hx7flFR/f52XAecCZKm5AdDhHr3+fDvj5fsx92uLz4ihgXy+PKwURCWEI+MVKqcfs71uFvlLqSRH5iYgMUUoVtNiWj/9jQeakT/4d+JtS6k37G6XyfcZ5U0SOUUq9ETdt7XHYJ4zhRzAZDqxKd+K+bK45C9iilGpzelNE+ovIAPNvDOfiy0779hY2O+Z/uFx/DfBBERkd11q+ACwtxPhMROQzwBxgulKq02WfYn2ffr6fpYAZqfBZYIXbg6q3iPsAfga8opS6y2WfYaavQEROxrh/C/ow8vl/XArMjEfZnAIcsJgiCo3rar0Uvk8L1jl4KfC4wz7LgbNFZGDcdHt2fJs3hfYsl8oP8BDwFdu2RuDJ+N/HYURibAQ2Y5glCj3GXwCbgL/HJ8Ex9nHGX5+DEY3xWpHG+SqGrXBD/Oen9nEW8/t0+n6A72E8lAD6AQ/HP8dLwHFF+A5PwzDL/d3yPZ4DfMWcp8BV8e9uI4aD+9QijNPx/2gbpwA/jn/fm4DJhR5nfBz9MYT2UZZtRf8+MR46bwBRDLv6f2L4gJ4B/gU8DQyK7zsZuN9y7BXxefoqcLmf6+myBhqNRlPB9GVzjUaj0VQ8WshrNBpNBaOFvEaj0VQwWshrNBpNBaOFvEaj0VQwWshrNBpNBaOFvEaj0VQw/x/x7cUjcIsiRgAAAABJRU5ErkJggg==\n",

      "text/plain": [

       "<Figure size 432x288 with 1 Axes>"

      ]

     },

     "metadata": {

      "needs_background": "light"

     },

     "output_type": "display_data"

    }

   ],

   "source": [

    "# DATA LOADING AND BASIC INSPECTION\n",

    "# We use wheat data from BGLR (https://cran.r-project.org/web/packages/BGLR/BGLR.pdf)\n",

    "'''\n",

    "Matrix Y contains the average grain yield, column 1: Grain yield for environment 1 and so on.\n",

    "Matrix X contains marker genotypes.\n",

    "'''\n",

    "\n",

    "# load the dataset as a pandas data frame\n",

    "X = pd.read_csv('DATA/wheat.X', header=None, sep='\\s+')\n",

    "Y = pd.read_csv('DATA/wheat.Y', header=None, sep='\\s+')\n",

    "\n",

    "# data partitioning into train and validation\n",

    "itrait=0 # first trait analyzed\n",

    "X_train, X_test, y_train, y_test = train_test_split(X, Y[itrait], test_size=0.2)\n",

    "print(X_train.shape, y_train.shape)\n",

    "print(X_test.shape, y_test.shape)\n",

    "\n",

    "# print basic statistics: max, min, mean, sd\n",

    "print('       min max mean sd')\n",

    "print('Train:', y_train.min(), y_train.max(), y_train.mean(), np.sqrt(y_train.var()))\n",

    "print('Test:', y_test.min(), y_test.max(), y_test.mean(), np.sqrt(y_test.var()))\n",

    "\n",

    "# do basic histograms\n",

    "plt.title('Train / test data')\n",

    "plt.hist(y_train, label='Train')\n",

    "plt.hist(y_test, label='Test')\n",

    "plt.legend(loc='best')\n",

    "plt.show()\n",

    "\n",

    "# marker PCA, use whole X with diff color for train and test\n",

    "X = np.concatenate((X_train, X_test))\n",

    "pca = PCA(n_components=2)\n",

    "p = pca.fit(X).fit_transform(X)\n",

    "Ntrain=X_train.shape[0]\n",

    "plt.title('PCA decomposition')\n",

    "plt.scatter(p[0:Ntrain,0], p[0:Ntrain,1], label='Train')\n",

    "plt.scatter(p[Ntrain:,0], p[Ntrain:,1], label='Test', color='orange')\n",

    "plt.legend(loc='best')\n",

    "plt.show()\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 4,

   "metadata": {},

   "outputs": [

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6X0lEQVR4nO2deZwV5ZX3f6ebbmgQaZBFuNhC1LSKKA0dgyGba0dx6bgRI/MmTmaczDuTxGU6gdE36oREZzrJmJg46kxikokxuGBHQxI0bokGiWCDiNCugNwGaYRm6waa7vP+UVW369at5alb213O9/PpT9+qW7fq1FNPnec85znPeYiZIQiCIJQPFUkLIAiCIMSLKH5BEIQyQxS/IAhCmSGKXxAEocwQxS8IglBmDElaABXGjh3LU6ZMSVoMQRCEomLVqlU7mHmcdX9RKP4pU6Zg5cqVSYshCIJQVBDRJrv94uoRBEEoM0TxC4IglBmRKX4i+ikRbSei10z7WoloAxG9SkSPEVFtVNcXBEEQ7InS4v8ZgM9Y9j0F4BRmPhXAGwAWRnh9QRAEwYbIFD8z/wnATsu+J5n5sL75EoDJUV1fEARBsCfJqJ6/BbDY6UsiuhbAtQBQV1cXl0yR09aeRuuyDnR292JSbQ1amurR3JBKWixBEMqIRAZ3iegmAIcBPOB0DDPfx8yNzNw4blxOGGpR0taexsIla5Hu7gUDSHf3YuGStWhrTyctmiAIZUTsip+IvgjgQgBXc5nlhG5d1oHevv6sfb19/Whd1pGQRIIglCOxunqI6DMAvg7gU8zcE+e1C4HO7l5f+wVBEKIgynDOBwEsB1BPRFuI6EsAfgRgJICniGg1Ed0T1fULkUm1Nb72C4IgREFkFj8zX2Wz+ydRXa8YaGmqR8sja9DXP+jhqqmqREtTfYJSCYJQbsjM3Rhpbkjh0pmDETyp2hrcful0ieoRBCFWRPHHzMy60QCAKxsn48UFZ4nSFwQhdkTxC4IglBmi+AVBEMoMUfyCIAhlhih+QRCEMkMUvyAIQpkhil8QBKHMEMUvCIJQZojiFwRBKDNE8QuCIJQZovgFQRDKDFH8giAIZYYofkEQhDJDFL8gCEKZIYpfEAShzBDFLwiCUGaI4hcEQSgzRPELgiCUGaL4BUEQygxR/IIgCGXGkKQFEIRCpK09jdZlHejs7sWk2hq0NNXL+shCySCKXxAstLWnsXDJWvT29QMA0t29WLhkLQCI8hdKAnH1CIKF1mUdGaVv0NvXj9ZlHQlJJAjhEpniJ6KfEtF2InrNtG8MET1FRG/q/0dHdX1ByJfO7l5f+wWh2IjS4v8ZgM9Y9i0A8DQznwDgaX1bEAqKSbU1vvYLQrERmeJn5j8B2GnZfQmAn+uffw6gOarrC0K+tDTVo6aqMmtfTVUlWprqE5JIEMIl7sHdCcy8Vf+8DcAEpwOJ6FoA1wJAXV1dDKIJgoYxgHvd4tUAgJRE9QglRmKDu8zMANjl+/uYuZGZG8eNGxejZIKQHb3z4oKzROkLJUXciv99IpoIAPr/7TFfXxAEoeyJW/E/DuAL+ucvAPhNzNcvGNixryMIghAtUYZzPghgOYB6ItpCRF8CcAeAc4noTQDn6NuCIAhCjEQ2uMvMVzl8dXZU1ywmiJKWQBCEckVm7iaEuHoEQUgKUfyCIAhlhij+hBBXjyAISSGKPyHE1SMIQlKI4hcEQSgzRPEnhLh6ChuWLplQwojiTwjRK4IgJIUofkGwQRpmoZQRxZ8Q4uoRBCEpRPEnhFiUgiAkhSh+QbBB2mWhlBHFLwiCUGaI4hcEGyScUyhlRPELgiCUGaL4BcEGsfeFUkYUf8wQJI5TEIRkEcUfMyy2ZFEgLn6hlBHFLwiCUGaI4o8ZcfUUB9IzE0oZUfwxIwpFEISkEcUvCDaIj18oZUTxx4y4egRBSBpR/DEjrh5BEJJGFL8gCEKZIYo/ZsTVUxyIj18oZRJR/ER0PRGtI6LXiOhBIhqWhBxJIK4eQRCSJnbFT0QpAF8F0MjMpwCoBPC5uOUQBEEoV4YkeN0aIuoDMBxAZ0JyeNLWnkbrsg50dvdiUm0NWprq0dyQyvt84uopDqRnJpQysVv8zJwG8F0AmwFsBbCbmZ+0HkdE1xLRSiJa2dXVFbeYADSlv3DJWqS7e8EA0t29WLhkLdra03mfUxSKIAhJ46n4SWM+EX1T364jotPzvSARjQZwCYCpACYBGEFE863HMfN9zNzIzI3jxo3L93KBaF3Wgd6+/qx9vX39aF3WkYg8QnzI4K5QyqhY/HcDOAPAVfr2XgA/DnDNcwC8y8xdzNwHYAmAjwU4X2R0dvf62q+CuHoEQUgaFcX/UWb+JwAHAICZdwGoDnDNzQBmE9FwIiIAZwNYH+B8kTGptsbXfhXE1VMcyFMSShkVxd9HRJXQ3wUiGgdgIN8LMvMKAI8AeAXAWl2G+/I9X5S0NNWjqjLbQq+pqkRLU31CEgmCIARHRfH/EMBjAMYT0bcBvADgO0Euysy3MPOJzHwKM/8NMx8Mcr6oaG5I4bOmCJ5UbQ1uv3S6RPWUAbLYulDKeIZzMvMDRLQKmkuGADQzc0G6ZqKgoW40Hlq5BZ/7yDG447JTA59PXD2CICSNp+InojoAPQCeMO9j5s1RClYoiOFXHIQ930Ieu1DKqEzgWgrtPSAAw6CFYXYAmBahXCWLuHrCx5hvYYTeGvMtAARS/oJQqnj6+Jl5OjOfqv8/AcDpAJZHL1ppIq6e8IlivoX09IRSxvfMXWZ+BcBHI5BFEPIiivkWglDKqPj4bzBtVgCYiQLOrVPoiKsnfCbV1iBto+SDzLeQjplQyqhY/CNNf0Oh+fwviVIoQfBDS1M9aqoqs/bJfAuhGGhrT2POHc9g6oKlmHPHM4HygPlBJZzztjgEEYR8MQZwr1u8GgCQqh2GlqYTA0b1iMkvREuSQQmOip+InoBLh5eZL45EIkHIg+aGVEbxv/CNs6BlAxGEwsUtKCExxQ8tdbIgFB3MgKreDzv+XygN4qgXSQYlOCp+Zn4+8qsLQoK4dbU/9eFkUoELyROXCyaSoARFVPLxn0BEjxDR60T0jvEXuWSCkCeq3nlZb0GwI656kWRQgkpUz/0A/gvAYQBnAvgFgF9GKZQgxIFbV1uGdsuXuFwwzQ0p3H7p9Mx2GEkgVVFR/DXM/DQAYuZNzHwrgLnRiiUI+aOaWTOK9RaE4ifOemFW8i8uOCu28SUVxX+QiCoAvElE/0xEnwVwRMRyCSVIXDHLqta623oLkpa5fGlpqsfQIdmqsdTmhago/q8BGA7gqwBmAZgP4AtRCiWUHlEsXB+U5oYUPn96XWY7zq62ULg0N6Rw/TknZLZLsV6oKP5+Zt7HzFuY+RpmvoyZX4pcsgJDwsKDEedAqh9j/SNTxwAALph+dFZXW+z98ubskyYAAI4ff0SsLpi4UFH83yOi9UT0LSI6JXKJChTp+Qej0BOpSQ4loZxQSct8JrRoni4A9xLRWiK6OXLJhJIizgGzfNItWH8jDb0AlO4SnEppmZl5GzP/EMCXAawG8M0ohRJKD0mkJhQTpe7aVUnLfBKAeQAuA/ABgMUAboxYLqHEyE2kFl16hHyMNKurR5K0CUDpjvWoLL34UwC/BtDEzJKHX8gbcyK1FxeclawwOiXakxcEV1TSMp9hfCaia5n5vmhFEoQCQBoEoYTxu/TilyORQhASotR9uUK+lHbF8Kv4QykNIqrVE79t0ENFz/D+lSCo4cd943SsGPwCgJKtCK6uHiIaB+BYAG8xczeAi0K67g8A/IGZLyeiamgzgwUhFPIamC1tA0/wSan3BB0tfiL6OwDrANwFYAMRXczMW4JekIhGAfgkgJ8AADMf0hsVQUiOErXshGCUarVwc/VcB2CaPrj7MQALQ7rmVGiTwe4nonYi+h8iGmE9iIiuJaKVRLSyq6srpEsL5UAYkToS7SOUMm6K/xAzdwEAM78DYGhI1xwCYCaA/2LmBgD7ASywHsTM9zFzIzM3jhsnqyEJEVPiXXtBMOPm459MRD902mbmr+Z5zS0AtjDzCn37EdgofkHIFz/GutOxMoGrvCl1O8BN8bdYtleFcUFm3kZE7xFRPTN3ADgbwOthnFsoPeJeDN18PUEoVdwWW/95hNf9CoAH9IiedwBcE+G1hCIl30Wv/STWMlt21uuZ95daWl5BjbJO0hY2zLxa99+fyszNzLwrCTmEwiaOHP7m19ruesZ+QSglElH8xUipx/UWIvnm8M/HRqMA1xOEYkMlSZsgJMKk2hqkbZRuNDn8472eUNiQgqUX9/hTmLhN4BpCRP9ARH8golf1v98T0ZeJqCpOIYXyxG0xdDfydcvarRlg7E+CuBanF5xxqkqFuIa0H9xcPf8LYAaAWwFcoP/dBuA0AL+MWjBBaG5IYe70iZlt5UWv88zY0NyQwu2XTreVI26KXbEUO172fpxrSEeBm+Kfxcz/yMwv6Qutb9E//yOAhrgEFMqbU1KjAABf+vjUSBa9tkZtFEpXvdgVS6ng1Hss9vEgN8W/k4iuIKLMMURUQUTzAEgUjhArvjJulsDkq2JXLKVOnGtIR4Gb4v8cgMsBvE9EbxDRGwC2AbhU/04Qih6VQbwkKHbFUurYjQdVVRB6Dh0uijEZR8XPzBuZeR4zjwNwBoAzmHm8vu/d+EQUBH/htP7y8Rdm70AWp08Wr1phjAcZz2h4VQVAwK6evqIYk1GK42fmD5j5A2ObiM6NTiRBiJ+oLX+/ETqGYjGimsYeMVRtYFsIFTe3YXNDCudNmwAAGFpVib7+7GMLeUwm3zj+nwCoC1MQQQiLfGx4L8s/SMx2vqknmhtS+MXyjXhlczfumT8TjVPGKN6NEBS/PcFdPX22+wt1TMZR8RPR405fATgqGnEEofDIV3EbuEXoeP3e6IkUgkOqmCcs+cVveY8eXmWr/At1TMbN4v8EgPkA9ln2E4DTI5NIEAKSj9/ezdUTRHEDwSJ0KnSxkh6KCNr4lSrGc5l76kQ8uiqdVU8KeUzGzcf/EoAeZn7e8vccgMJ0XAklRz4KL2wdGTS0MkiEDulTiQYS1Pxt7Wnc+NCasppX4Le4G48dkzX5L1U7rKDHZNyies5n5mcdvvtkdCIJ+SJT/KMhaGhloAidhC1+w9LvdxCgUH3YYeFV7uaOolnJv/CN8Ccbholk5ywRSnWKfz7BNmEryXxzBhlYU0Eop56AydWTkJffKVW1QaH6sIMTrLyTds154an4iWgvEe2x/L1HRI8R0YfiEFLwplSn+BfCC9TckMLnTx8MYvOjuM3nMPCTeoIyJr/ypULFzaIvZB92WHjVP6fvC6DauqISznkntHVyfwWt4/k5AMcBeAXATwF8OiLZSppC80MXOn4M/3ysY6/znz71KPx8+SZcMP1o3H31LN/nD0pSisQpVXUlUUH7sIMS1ODQAgwKc1Y4oKb4L2bm00zb9xHRamb+BhH9a1SCFQrFkvclrFzyN7etxYMr3kM/MyqJcNVHj8Gi5tyMlXET9VNQPX/cPRBK2Mff0lRvuxzl9648rWSVvh+cXJGFrjVUFH8PEV0J4BF9+3IAB/TPhX5/kRE0pjlsW8DuBfXbFb+5bS1++dLmzHY/c2a78dgxxRPDHUGtTCqlT0Umjj++V81aty+blcqqF0Dph3AmaQjEMV9CRfFfDeAHAO7Wt5cDmE9ENQD+OVRpioQwYprDri/Gdb/+yKs41D+AsUdU4+a5J/uqMA+ueM92/wMvbc6KUU4ihjtqvVuonXKjwRlQrDBBlYZd3X50VXEHCCRBvg1CXPMlPAd3mfkdZr6Imcfqfxcx81vM3MvML4QmSRFRqAOpzQ0pzDimFgDwo8/P9F1RnEL2GCjI+3XCzzsXlwslaDI4ld+HEdnlVLdLCZWwZ9XH5ezqye95x6VbPC1+IpoM4C4Ac/RdfwbwNWbeEqokRUQYA6mRWZgBTlxJ5Kj87Qg6cFxWKQB86AFzuVQP0WwzlZ8HnWEMlE4wgBN+LWqvBtcxqifPdj6uIA2VOP77ATwOYJL+94S+r2wJI1d6IQ6OXPXRY2z3j6jOXYcWCBbDHeW8Az8vXVy+e1WRrOVy8PAAAOAvb+3w/G0YSqN04/I1VC1qvxZ7WPUornUYVBT/OGa+n5kP638/AzAuVCmKjFLNlb6oeTrO+NBgBshKIsyfXYdvf3Z66PcbZZfWz0ublKvHyd3gNGFqySveDWIYSqOlqR5Dh3irBUPeYpstrto4Go9LNV13WPUoLt2iMrj7ARHNB/Cgvn0VgA9cji95jC7hTY+txf5D/RhVU4XbLp6WaFRPhoAV8LJZx2D5OzsBAG/ffsHgaZlx/UNrAGgTmIK6ZZRfwALrG+X7gpt/5uZucCqXD/Yf8rxGS1M9/uXhNThsGgn2qzSaG1LYvHM/vv/UmwCc3X9GA11sidv8hj17uXocffx51hOj3K5bvBpAOO+aHSoW/98CuBLasotboYVzXhP0wkRUSUTtRPTboOdKguaGFC6dORkAcON5H855MF6WUNTqLN+GxamiXzIjv5mnTkTZpY3Cig/aUJtlcuvtON3/mBHVntdobkjhnJMmZLb9zjA26ux/6kp/cu0wx+Rwnd29vntthdA7ULWoVeuQ88xdtR6eHfnO8vaDp8XPzJsAXBz6lYGvAVgP4MgIzp0obhZdXISp+9ra0/iPZRuytv1WRmPA0uDME8cppbGlgg20zB+33s5/zpthO2HqszMmKZ375ElH4g/rtuErZx2PG89Tt/StdRYA0rsPoNYlz7xqr62tPY1bH1+H7t7B8yTVO/BrUau+R1bL39wg+BlQtr4n+bxrKjha/ER0FxH90OkvyEX1SKG5AP4nyHmSxskNoWIJFWJUjx1Gpe3sPpDZ53cQ1jxgafDoqjQumzVYoZ2sUy9Xj9mSGvyNf1QH5zq7e/KyWs334dbbsSZ0G6b720//ULRrH9nVWWbgQF+/bZVqaapX6rUZz96s9A2SCglWsaj9uhitlr95U7VnZPeeRJVo0c3VsxLAKpe/INwJ4OsABpwOIKJriWglEa3s6uoKeLlosb4YKpZQYXmunQljENbpHM9uGHyuXl1aO8VsjYAxeHLdNmXZ/PJa5568IpHMiqGlqR7VLtk+zeXwyQ+P038fbY1xqrO9fQO2dbW5IaXkNvHK7lmq4aPm56XaM4pzfpCjq4eZf27dR0RHM3Ogt4qILgSwnZlXEdGnXa5/H4D7AKCxsbFY9CSA8PLm5EXQ5FKW7TBCBMM4h53ec1Iq9zz/Nq6ZM1X53H6wzqD1GycPaEpzzZZu3P/iRgDu7oa4cvU4uXTcMOS98eE16B9gTDhyKBaef1LWfXg940INH3Urb/M8i2F6w+fWY1TVB3EmWvSbj/93IVxzDoCLiWgjgF8DOIuIfhnCeQsGFUuoYFMQWCp8GIOwUQ3kOr0Q2/ccDHTesOQwY1UkH52quW7OO3mCa2/HGOOIUu+3taex78BhX78xXF2A1nABwEP/cEbOfbg942IIgbY+N2sv0zA8Xt64M/t3ps9O+uDME8dlyhDQGl87omgc/Sr+wPqKmRcy82RmngItxfMzzDw/6HkLCZWFNwol26QXQRchMc4RNDb51y9vzvGrO70Q448cqnzeMMJFVV7MfK9Tob+hfi1+leON8ZHrFq9Gn2oyIOP8GHR17Tvo3FOwe/aAtjh5MaZ1duplLn11a9a2ufwNfVBdqT3MsUdU47JZKTy6Kp3VE9h34HDgd00Vv4r/v0OXoGjw1+bFEZJlYB7gXP1ed6BzWRVUc0MKV8yanLVP9YU15Lp+8eqcSUGq53hty24AwL6D/Tl+dSelcu0nw18fyKkrr/pi5uuq8bvmrmottRtItMOufM309vVjd69zb8FqBBm0f/O8olD61vfBqXeX4yazPK7mhhRm1NUC0PJoPbuhK6cB6RtgjKjO9r5H1Tj6UvzMfLf3Ub7O9xwzXxjmOeMkiN81LFePtet5qF8bL3/hzfAGxGcdOyZrW1Xpm+WyRnWoVubn3si9D7Nf3U6pAMCM257ElAVLMWXBUjT825OBIyNWvJs7Z5EAXDYr5Xkvbe1pnPP957O2Df70Rpd7hFBEC3B5DboCwJAKcixfM/0evYViUPBWnN5tp97daIubxq6HZwz4EpwbkN15vid+kTV3wyCPRB1hvchOL/DP/rIx0ckyKopFhT0Ovud0d69jjPN3lq7Pamh29fSh5ZE1rmXg9QR/s7ozZx8DWZFJdhgN4Nbd2eGwv1qxCQBw4PCAa4RQ5y5NQXz1wfas5xh0MpTKuMT4kUOVFE9lRdQjVvHj5Jpz6mXOPXWi9zlNaSDiysnjhCj+IsfpBd53sD/RhdfDikQ4cpjzHEOne+q3eWf7+jlQWNzO/fZ+bK/7dArRe8Em6Zpd6N6ruqsLGHyON7etDZzgTkXBjBruPVu4pqoSo2pUMr+UBkYv02jrDBdmo6VXbNdj2LFPCzqYd+9y7NyfG4AQ52C3KP6ECMtGUrUQVOOBwwgbbGtPZ1aOCsqnP+ycD9BvjHOQxmhEtf2rMqrGPhLD65pO3hHr8dY8Ob19/XhwxXu2jcl1i1fje0+9AQB4ruN9V7mcLFc/TBw1DLdfOh1HDM0tA2uPpNhwew+aG1KZ5376VE3h58zctfymrT2NTTt7Mt/19uVOYZpZNyo2t5gofiSTQyTMqBvVF1gp7DCgPIZrw09efzempUa5fh9bymGHhuzQ4X7XuuN0TSfviIqMKmX7Wude3NzmnCbEOj5Sa9OAeU0a++MNn8qOVNMPt0u5HTdhvdNORWBk7XTOx5/9ReuyDk+j6i9v74ytV172ij/KvPBxYH2BjZAxO/JVfH5UuJNvvzKixPd292SnVKsqKVA3ev9B+/GKnr4B17rjFMr68ePH5pwr7K6+01KaBmal/Y3zT/R9fqdIo7DGd/IljHc6qNmiOhHS+pu4UliUveJPahnFMNWg+QU2ll60Kr8www7drCln14bXSkb+XzWne7KGn46qqULr5ac5rLCkdi2V7JhAbt2xizy6/dLpuHr2sVn7/GbSVCGsXpcT1rMbbXvSaRgKcWlUVaMrrrIre8Xv1A1V6Z4Gea2s+dnDcjUZ0Qj1E0Zm9oWtVKzW1PWLV2fcCvlGK/jVUW73NKNudNa21SWRD82K2TGB3JfXem3r9lfOOj6SuR5R9bIM2CHTVtDIlKgilvylCFGrkE7RP9aftzTVKwX/SVRPDLhVKD8vTZDXK2xXkzFoOP7IYQCAL3/quNCVSk4WRwAPvLTZdVKVV2/D6TVz2u92T1EYuh+1ZMck5MZuG1izU9oNbpplHFIRzWvotJSmKl7l6DfkUYUw3ocwQyWdisD6zudm58ydCJkaNSyzXVOV+8xrqiokqidqjArmhFM32WyNBLHMjYoTdrfUaqmEYfSt2rTT8xjDP+k0qSrOSTzWl06lDFSX2DMYWlWBWy6a5jrF/ua2tbh+8eqc3qO13gypDP6QrIrk5EkjsajZe/JVEJwik+xSlqgSxvsQRooQVdvBsXG02W8Oj13/rfNxzknjs77/TnN8KSzKVvF7DUDZVVarNdJzSPt9PmkSjHoRdUY+3zlebGrs719TS8hqyJxP5XXqWhfq1CBm7T4vMk3cMbuf2trTeOClzbYKxKrENmzdo+zaqK2pslVqt196ata+T304W6lEgdu4jTVliSphvA8qubJU8Xp/jDLwCue0O5fVqxBnCuLymXlhwa0iOVkHzgmaOrH87Q/Q2d2LST7XyAw7hXPgaASbE3QrpusN4p/06+pJGqOcTp1ciyXtnfjCGcfitktOyXzfuqzDUfbO7t6sFBBPmBJ8ua3WVlNViVsvngYAuOGh1RmL21BqxqpScRGFSy2s98FcHn4aHgOve8snXba5oWxrT+PpDduzvr+pbS0qYpoFXbYWv1NFqiRytA7cFqtQ8Um2tafxnd+tBwD8bu3WjE88zIx81u73m+/vDTxw7JQu1kzQUMQwlUiOvzUCBWW8xG7K3YlRNVV44KVNjt879USNetnckML4kYP+4qjcA2937cta2cyKtZcWRjmH4aaJE+c4fvfftS7rwGHLy9rbNyDhnFHjNAD1vSvtw/6AYLNkrUvQ9Rzqz1h2l80cDD8MHIFjqXHPv9EVeOD4M6cc7fp9GFFDYaRHHjxX+Oe22mFWN4d1jMCtrnT39uGQXV4JD7LTekffFzo84H4Vq5ERhkSGm8YoTWN2cPyJ3gJG9djsN9eZOBddsaNsFb9RwYZXV+bsdyLILFm3QauZevjhFbMmB47AsVa3XKvCfaDMrhrPtIRHTqodtDYvnZkKJWooH2vRKWLGdQHUEK4LOA9sGrQ01eekog6TqFfkUqH5xy9g6oKlGddMWMtDNjekMGKo5oV+8vpPJpLd0/tWvGbuuu+TJG1FhHXQyNpomElqWTWjcrnV23R3b6D5As+3nJnX79x4Yk2nL5fUlAVLbSNm7HAqi7b2NBYtfR0A8IfXtoU6W7u5IYV/OS8/94SKceGll+JoGLbtOQjGYFrmp9a75wfKh+TbN3cJnAa47faaj7Vz8Q6TcM7oMVwvRmSOeb8Z63M1Wx+XzNA+q/jo42rhd+lZ//5kk8fejKPbR0FjmCswhRR30/LIqzmTwqYsWIq7n33L8TeqA8J2t2Q8f2MBjd6+/tBTdZx78gTfvxkzolopB37Ui6/nw4+eftN2v9s4gRdJ3aZyOKficW3taWz6oCdrn3WG+aJLTpFwzqhxitDJZ3CluSGVSZFw9JH2PkmVQasw6njalPfdC7Pbx3Cb/L/frPP8nfllDDJPwE3JGpdwysfvhtdkGqCwpvXX1lRhsu4+27n/UKKpBYKw92B/1noBBl712m6mrlGt7Bq4JJIqWvGK6jHLbRgZZrfrwiVrc8rlIh+zw4NStoo/TNfLjGNqM2lan/jKx21bbTs3URSDVn4tJGNBE7tl+JxeKNVlAL2ISsGpRJv4ef5z7ngGK97JXYHLUw7F486dNgGdpgZbyX1VeAY/gMFnqvps7WbqXr94NfYe1Bp861q2USZVNDco19z/VwDAjn2HXBsXx16n6Yvbnlhna2RY58fE+UzLVvGH7XoxGnPDurSzVMxK/oLpEwtiSTqCfcUEgBsfWoOpC5ZmQlANfK7L7UhUEQy5UT25+Hn+6e5e/K9L+GU+DDMN/D617n3lMjWU0wf7DwW6vnVgvH3zrkDnMzCeqeqztet5mYviW0tfz1K6UfXUrA3KPlM2VrvGxa1HYj1vznq8OqrzY6KgbCdwtTTV418eXpMT9eJvcGXwtwMDhsJH1n8lYpqeSrBXik4V00hbYa2g2T5+rXK3LuvImsCmgtNknaDkxvHnPoyWpnosXLI2S4kYrje7iVBO4ZfmU1vL4Ysfm+Io48dPGIs/rtcm8FjXI3bCUE5WxWe3BKWde8uQL93dm1MXHgvJXWI0nKrP1quBONA3gFsfX5cp13zmTajgNZO/t68fNz60BoAlrNYjqsetQaqpqrBdkCUOytbib25I4ZyTcgff8rHCCZRRkipRNdpvNMyTun6zOo2Gf3tSyXepusKReSJgWD1Jc1bGzTv323a9VYgrgsHu5bS63rS0B/m73t7u2pdTDt+29JTMmOP+7RZBscPXuJTlntva02h5eM1g6KXl8L485hVYGTpkMCrFT+PvRXdvX6Zcg5zHDZWGo585Y/lnfPweb5Xbea1zP+zqaVTjGGWr+AHgpIlHhnauwZmczq4e88NbunZrZu1Uw6Lu62fs6ulTmgGsusKRWQynZFl2+V/cMCewW5veY9v1VsFJyd45b0YmBj6fWHhVFWa+/vmnHO1b6be1p/GDp7WlDl94c4evxUfMr/y50yZAJU9bkBTitz6+Dn1h+egc+No5J2TKULUsw1gC0s/MXquLy3i/8p2cuWHrXgBa+WYPaLPnea0RhXYY73fLw2tCVf5lrfjDTFduvFMDDha/NRtoz6F+PPDSZs/upZ0152eFI7McLU31qKzIDT299eJpuP3S6a6rd5kxu3pUKq9fmhtS+MgUbS3TximjPY7OJYpUAtU2mnnhkrXY3asNQga5xMMrtzguxWjGKVW4SgpxVXdSEM6sV08MZ+6pqoSvOjH2iGrlnppdEINhXNnF1TvR2d2L3kPaczdcxbt6+rJchEada2mqd/TkWnt6br2HvgHGrY97R9ypUraKv609jf/+8zu+f+c0mDPo4892+Rh4DWI5YddVzNef2dyQwpzjsvPKD9PT+TY3pDKrd3lhVvxuk9jMPPbKFjUhI0A1vYHTTOBUbQ2unl2Xsz/I8oJbd2c/QxVXr1OqcK+VtuIKd/TTwDq5Bf/+E1OVfm8YKXdfPUu5d+E2MNzckMKVjZMdfpnNpNqarMFfO556XZvM1tyQsq07NVWVOHeav3keYTbesSt+IjqGiJ4loteJaB0RfS1uGYyWf28IMeIGAzk+/uwD81XWdl3FMCd97erpy1g9O/YdVPqNuQxOTY1Scse0PplcbLqTh8OsEJ9Y04mWR9bYuk1eXHAWTp9yVM7+ILy5fV9o5/LKdx/XvIB88gdZe7WzP+RdzjVVlUiNHuZ5nBWvEN5Zx47JnN/t2i1N9Z5RWGajclHzdDQeW5v1fW9fP55clz3TudTDOQ8DuJGZTwYwG8A/EdHJcQrg5irxso7M1q7xsf29XZmKcPk9f0FbezrnIeajrAnAmSeOy9nv5he1diuHWHwIq2zC9nr7+nHr4+uw8YP9SnKZy2DK2BG44dwPe/6mszt7Ylmck27m3bs853pW11vfALsOcP71Xf9x/G4cCDGao6WpPte9ZfqsYnSoDjBHgVk+L4VqJAQ8asRQ39dRDeE9eeJIVNsYM7U1VRm3kpdrrmtvthFVN2ZEzjG7LRa8l953WvEtH2JX/My8lZlf0T/vBbAeQKwB7W4vgtdkELuKueSVwePf33MQLY+swRNrOrOOcZq5+4njBy0c68NgAI+uSmfJY4TkmRsuN4vvzPrshmO/Qxe1u7dPOZa8+ccvZj4TAefoqQmmjh3hWDmt70nY6RHMfHfZhqzt7XsP5lzPzzgJAPxmdaf3QT6wW3ovX5obUq7WoorRcZ5Pt4Md+VqsFaYxCq/JgUZCQOMoP+N0Tu/gmSeOw5w7nsENergmAPzdx3NdTs9//cyMW2nkMPdI+HEj/TdMblRVEm65aFpo50vUx09EUwA0AFhh8921RLSSiFZ2dbnnnfGL24vgNRnErmL2WzRmXz/j64++mtk24qztZu5e3ji4LupQGyvemlbBboZtz6FBl5XXuzdiaLAICkBr3Aze7dqfUep7ew9hn4P7zCqXOS46bA4czi0Fo1djhMD6nT+w02OyzYihlWg6WX1ws37CSF/X98LtudspPOuAcEOd/0H0sOi36UV7YfRw/MRnGO+geQyXwFj88ntZ9WH1e7sdLjr4cXi1u+L/kk3D4YXT+OHRo4ah9XLndPH5kJjiJ6IjADwK4Dpm3mP9npnvY+ZGZm4cNy7X3ZEvbe1p7D/o7tsPY0ap+Rka1qbdzF1zQ+JkgRryOFmpThOwAOCld7LXy51l84LXVFXm3Y1c8e5OPK1PRNrVe9hXyKDXoGTYqMSD2zF1wVLPrv0Vs47B33/yuKx9412svo739/qUwh23GaTNDSlcNitbadSN8e8j95Yh+Dn8pgPJJzLPfIWevoEcF18/Mx5a+Z6rbGZjy46zTwpv6cvH/2lO6LP8E1H8RFQFTek/wMxL4rqudTEUJ8w9AnPFamtP48zvPpfZ3rhDzSfu1os4bKp0Tn57Q558GqS9lkbuuPFHZG0b4XC3XDQt75forme0WHZrz6dUYNi7+Lziz3/xpdMdv7PO2Ayy4t6cO57BCTf93vH7tvY0Hl2V7VbbvFM9mZ8qYSwOY6f4J44Kr5FqXdah5NLcsS83JYZZNqeerUHObSg8X4b92FcUBlISUT0E4CcA1jPz9+O8topft6qCMj6/qQuW4j/+oCnsjTu0mZlbTcm0Xt6UO1DqhJ3SbmtP49tLB2d3Hn1kdc5grHlySj4DxCOHundJjXC45oYUjh0z3PaYUR4Df3sOaGVqnSNgZliePu2VG9XLOE4MV50Za8M54GP8VmV5S8BeCdr1Yt409Sjs6r1VmajOtnYjHIs/d98TX/l47rUcfm+e0W6HqvF01IjqnH1z73oho5hDmOiMUTXZ7+ZvX+20fQ6/X7s1Z19QkrD45wD4GwBnEdFq/e+COC6sOi3b8PkxBkfel7+zM/fl8WHhWpX25g/25/Q+Nu88gFNTg7OJU7U1uGxWCq3LOjB1wVLsP3hYeZKJwRnHjVE6rq09jfd29th+d8L43IgEO0YOHeKo4L95ob/Are17tAb24OHg0S/jRw4NNSICAOYqJNnzYwHv3K8Wo/3Vs49XOu7Pb+7IKKm4lvMz322+A/d2LqtDLnXgT2/syCj6Gbc9mQnJdZJFxXiqJMKVprE3g227D2DhkrW4uc27kVR58l87Ozsa7q6n37I1TO/9k//5Rl4kEdXzAjMTM5/KzDP0v99FdT2zBVCh4MsY4HDylpgh5OYueXnjLlsr7K0uzX00Ymgl0t29eOClzZmK3N3b51u25W8P+vjtwkyJBvO4OJ16zRaHwS6dI/UIhxFDh2DRJafYHtM0zX3dXiubHBqhfOjaexDMuQvmBCGnJ29z6iiGMG567DWl4w4PcMa9GNdyfmalne/cAbsy6+vPVfzGcfc8/7bn+3Hd4tVo+LcncXPbWtvxPWu1mHHMKHzsePv5BL192ox71fsw9I858s/A6tba5rCWxvY9avNr/FDSM3etOW2iGkz0ioG+enadTfZEe4yFR4ywSxWJ3cI5zT7+hUtexdtd2ROHrrxnOW58aI3roKxXY2O8mNt29zr6TxsX/dH1HFbCsPQNDKUQ5jp+mxTGd6Kobn6GUQxL3z6qJ0ypNPzOHbDDrhdtHlez9iRU68munj788qXNOeN7o4dX4evnn5i17xgHl6eByiNgsGMUnhNOdmnYoaFAiSt+v7HafjF82k0eMdDPbuiKdMKS6kvW2zeAVywTuMJoEI2Byn4Gvvm4vUXq9wpRLFQeZpKylZt2ZaWhsF1cO+EVYw1L3xpKDACnpMJLUGjwfMdg2LVqKg8ge/znxodzQ3zNjy3s+R811ZU41yZLb9BGm9lb/1iv4VQ9r5kzJZgwNpS04o/at/mRY7XQSK8YaCMnSVT5avzUUacJXGER1oxUp4HmQoEB/Otj7r7eS370ouv3UVJpCVKwul4m1YZfvg+s0Fwgbe1p7PeRvM/PWJkRIRdWk7q1O9e94nZu1Y7SMxu2e1r6qobBp30kv1OlZBV/nCkBFihERPT29Sear8YgjAlcUTN1wdJQffx2hOHpMIdk2vr4Q7hGvpwwfgQeXZV2TN0dhRvKyPXk17/vd9zKbNCpZpR1YuKoYbbPyW7fhCOH4urZdUpjRT90WHg+6xqKtx1FmHRJKn7DtxZ2cVl9+S+9u9PhSHus+WqSYHJMA31BYITr4zcwp0mwy5gYhLe374Of5qQqSOC+Ahu27XN1M4S1brKVtvZ05D3tSbU1Ga35N6bnmE8ajBvO/XAms64XD//Dx7CoeTquOt277qjUX9UnEMWzKsmlF6Py7UehjLywWy4xCB3v22eFrKqk0KOZCokKAr7dfApueFhLpfHshnDTgLz49gdoqNuu/gOL3h9SQTnLgEZJVFdqXdYR2ZKaBt09hzLnN5dZPssYOk3OswsrveLev2D7noPKcy48ry0Wf7iEaXGYu3VRDhQDue6CClKf2BOUER65R4odZmDuaZMy22Erpv4Bzvi4VbA2slPHqs2VCCsax8gXHzad3b2+FjXJB/P4wa/+ql7mdvy/37yGP67PLQs7Vfv+noOua1T7RdXHH4XFX5KKP6y45aFDKnCJSVlEjuX51k8YGVol88KaIrbUOHrUsFBXMLLDbpq/KqpWXaFnxZhUW4PmhpSSOyQMgvZSD/QN4Gd/2RiOMD4xsgJ4YTONITAlqfiDruNpWCsfO+4onKa4KpUbXgtlGFir8O6eQ6EMQqrglZZBlQlHhh9zHAY9h/rx4F9zE28VCn2K+R0KWe9XmxZbj9DgDx1r7nwABVXQV967PPQF10tS8Rtxy8PyiAWvqiTM+4g2XXvb7gP47pNvZH2XD17ZQJ3o1LuWcRDHmqxJUug9mvd2xpNWIUqqKytw/eLVmHHbk/iFwuzWQsUpdUmSGCHhYSn/klT8gKb8T5hwhPeBFvoHODPKv37b3iyFkW+3stSVqpn3I5heLhQH+w4ezsySLuZMra9s7sY1P3s5aTFy8ForxA8lq/jb2tNYm85J8+/JAAO/KmCXgFAYhLmCliCoElbgSsnW3qgH8oTyxi5xmCBETViBKyWr+MvJvaKKNde/kD8JTOkAoC2cowIBOH1KcsspCtFgzfKbLyWr+IVcrAuvC8XHnfMalI5jaKm/hdKhpqoitCUYS1bxh73wRinw1HofM0uFguSwj2W9ind4VbBjWIAQdSslq/hvuWha0iIIQuh88f7CizYR4qE7xMmcJav4w16VXhAEIUnCXEmtZBU/EO5Se4IgCElRVUmhDewCJaz429rTJZ1tUhCE8qH18tNC9WKUrOIPa4abIAhCkkThtyhZxR9lPnBBEOJhzPDSTheuAiN8Q7ZkFX+l05L1giAUDTt78ktwWGqEvapZySr+/oiWlhMEQYibMCN6gBJV/G3t6djy2JcjqusLCIIQDmFG9AAJKX4i+gwRdRDRW0S0IOzzty7rkFmLERJ2JRQEwZkwUzUYxK74iagSwI8BnA/gZABXEdHJYV4jbH+YMIhRCSUlhiDEw+EBDnX1LSAZi/90AG8x8zvMfAjArwFcEuYFwvaHCRoVBNx+6akAtJQYMkFOEKKnr59LIqonBcC80skWfV8WRHQtEa0kopVdXV2+LhB0zV0hl9HDq/D9K2dkupzNDSm0Xn5awlIJQnlQNlE9zHwfMzcyc+O4cf7SCRtr7roNQlYQMH92XUE2EH4Wdwo7bLUCWtkY554/uw4b75iL9m+el+NnbG5IxT7QW1VBSj2NSiKcMH5EDBKFS1UFYbiPCkAARlTHU4dHD6/C/Nl1qJJ1HWKnFKJ60gCOMW1P1veFSnNDCi8uOAsb75iLO+fNQG3NoE/asF4XNU/3bCBqqiox57gxOVFCVRWE0cOrQNCiXObPrnM8DwGYc9wYpGprQABqa6oclVdNVSVar5iB+bPrPCOTaqoq8b0rT8PGO+Zi/uw6X43AiOpK3DlvBu6cNyMjV6q2Bt+fNwPv3D4XG++Yi7dvvwCLmqe7nqelqd5RERiKwlz2ADC8qiIzRmDITDbfG2VlLufWK05D6+Wn5ZzTjFEuT93wacw5boyr/DVVlbYyuh1vrQ/VleQo753zZmDjHXMzf3bP1dg27u/1b52feS7AYBmZ65lx/v+cNwPf/uz0UAyYalOdND8D4z7av3keFjVPR+sV7uVvvS83+Z0a56pK0g0zdxVlbShHD6/KqtdmOaxyhdV8GdU/VVtjqyvM36vck5WaqsrQAyqIY453J6IhAN4AcDY0hf8ygM8zs+NaiY2Njbxy5cpY5GtrT6N1WQc6u3sxqbYGLU31aG5IOe73cw67Y9LdvagkQj8zUiEcaxx/6+PrMquQjR5ehVsumhZZxtK4r2e+bphlaPfcACg/d79yR3XOUTVVINLS+E6qrcGZJ47Dsxu6Cva+VOqP33dARS5rmR063I+ePm29gwrS1t82X8uuHPMtM7v7qdWf266ePt/36AQRrWLmxpz9cSt+XZgLANwJoBLAT5n5227Hx6n4BUEQSgUnxZ9IIgxm/h2A3yVxbUEQhHKnYAd3BUEQhGgQxS8IglBmiOIXBEEoM0TxC4IglBmJRPX4hYi6AGzK8+djAewIUZy4KWb5i1l2QORPkmKWHSgc+Y9l5pwZsEWh+INARCvtwpmKhWKWv5hlB0T+JClm2YHCl19cPYIgCGWGKH5BEIQyoxwU/31JCxCQYpa/mGUHRP4kKWbZgQKXv+R9/IIgCEI25WDxC4IgCCZE8QuCIJQZJa34o17UPShEdAwRPUtErxPROiL6mr5/DBE9RURv6v9H6/uJiH6o38+rRDQz2TvQ1lAmonYi+q2+PZWIVugyLiaian3/UH37Lf37KYkKrslUS0SPENEGIlpPRGcUWdlfr9eb14joQSIaVsjlT0Q/JaLtRPSaaZ/v8iaiL+jHv0lEX0hQ9la97rxKRI8RUa3pu4W67B1E1GTaXxg6iZlL8g9ayue3AXwIQDWANQBOTloui4wTAczUP4+Etk7ByQD+A8ACff8CAP+uf74AwO+hrSExG8CKAriHGwD8CsBv9e2HAHxO/3wPgH/UP/9fAPfonz8HYHEByP5zAH+nf64GUFssZQ9tudJ3AdSYyv2LhVz+AD4JYCaA10z7fJU3gDEA3tH/j9Y/j05I9vMADNE//7tJ9pN1fTMUwFRdD1UWkk5KrOLG8KDOALDMtL0QwMKk5fKQ+TcAzgXQAWCivm8igA79870ArjIdnzkuIXknA3gawFkAfqu/pDtML0PmGQBYBuAM/fMQ/ThKUPZRuuIky/5iKXtj7eoxenn+FkBToZc/gCkW5emrvAFcBeBe0/6s4+KU3fLdZwE8oH/O0jVG2ReSTiplV4/Sou6Fgt71bgCwAsAEZt6qf7UNwAT9c6Hd050Avg5gQN8+CkA3Mx/Wt83yZWTXv9+tH58UUwF0Abhfd1X9DxGNQJGUPTOnAXwXwGYAW6GV5yoUT/kb+C3vgnoOJv4WWg8FKALZS1nxFw1EdASARwFcx8x7zN+xZhoUXMwtEV0IYDszr0paljwZAq3r/l/M3ABgPzRXQ4ZCLXsA0H3hl0BrwCYBGAHgM4kKFZBCLm83iOgmAIcBPJC0LKqUsuKPZVH3oBBRFTSl/wAzL9F3v09EE/XvJwLYru8vpHuaA+BiItoI4NfQ3D0/AFBL2rrKQLZ8Gdn170cB+CBOgS1sAbCFmVfo249AawiKoewB4BwA7zJzFzP3AVgC7ZkUS/kb+C3vgnoORPRFABcCuFpvuIAikL2UFf/LAE7QoxyqoQ1oPZ6wTFkQEQH4CYD1zPx901ePAzCiFb4Azfdv7P8/esTDbAC7Td3kWGHmhcw8mZmnQCvbZ5j5agDPArhcP8wqu3FPl+vHJ2bdMfM2AO8RUb2+62wAr6MIyl5nM4DZRDRcr0eG/EVR/ib8lvcyAOcR0Wi913Oevi92iOgz0FydFzNzj+mrxwF8To+kmgrgBAB/RSHppCQGFuL6gxYZ8Aa0kfSbkpbHRr6PQ+vavgpgtf53ATTf69MA3gTwRwBj9OMJwI/1+1kLoDHpe9Dl+jQGo3o+BK2SvwXgYQBD9f3D9O239O8/VAByzwCwUi//NmhRIkVT9gBuA7ABwGsA/hdaFEnBlj+AB6GNR/RB63F9KZ/yhuZPf0v/uyZB2d+C5rM33t17TMffpMveAeB80/6C0EmSskEQBKHMKGVXjyAIgmCDKH5BEIQyQxS/IAhCmSGKXxAEocwQxS8IglBmiOIXBBeI6CY9A+arRLSaiD5KRM8R0UrTMY1E9Jz++dNEtFs/dj0R3ZKY8ILgwBDvQwShPCGiM6DNypzJzAeJaCy0rIoAMJ6Izmfm39v89M/MfKGe+2c1ET3BzK/EJbcgeCEWvyA4MxHADmY+CADMvIOZO/XvWqFN0nGEmfdDS5x2fKRSCoJPRPELgjNPAjiGiN4goruJ6FOm75YDOEREZzr9mIiOgpZLfl3EcgqCL0TxC4IDzLwPwCwA10JL4bxYT8plsAjAzTY//QQRtUNrOO5gZlH8QkEhPn5BcIGZ+wE8B+A5IlqLwYRiYOZniGgRNKvezJ+Z+cL4pBQEf4jFLwgOEFE9EZ1g2jUDwCbLYYugZWgUhKJBLH5BcOYIAHfpi2gfhpaN8VpoufsBAMz8OyLqSkY8QcgPyc4pCIJQZoirRxAEocwQxS8IglBmiOIXBEEoM0TxC4IglBmi+AVBEMoMUfyCIAhlhih+QRCEMuP/AxZlPwZ6LO9kAAAAAElFTkSuQmCC\n",

      "text/plain": [

       "<Figure size 432x288 with 1 Axes>"

      ]

     },

     "metadata": {

      "needs_background": "light"

     },

     "output_type": "display_data"

    }

   ],

   "source": [

    "# OPTIONAL: SNP preselection according to a simple GWAS\n",

    "pvals = []\n",

    "for i in range(X_train.shape[1]):\n",

    "    b, intercept, r_value, p_value, std_err = stats.linregress(X_train[i], y_train)\n",

    "    pvals.append(-np.log10(p_value))\n",

    "pvals = np.array(pvals)\n",

    "\n",

    "# plot GWAS\n",

    "plt.ylabel('-log10 P-value')\n",

    "plt.xlabel('SNP')\n",

    "plt.plot(pvals, marker='o')\n",

    "plt.show()\n",

    "\n",

    "# select N_best most associated SNPs\n",

    "#N_best = X_train.shape[1] #all SNPs\n",

    "N_best = 100\n",

    "snp_list = pvals.argsort()[-N_best:]\n",

    "\n",

    "# or select by min_P_value\n",

    "min_P_value = 2 # P = 0.01\n",

    "snp_list = np.nonzero(pvals>min_P_value)\n",

    "\n",

    "# finally slice X\n",

    "X_train = X_train[X_train.columns[snp_list]] \n",

    "X_test = X_test[X_test.columns[snp_list]] \n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 5,

   "metadata": {

    "scrolled": true

   },

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "\n",

      "MSE in prediction = 0.7821123801773443\n",

      "\n",

      "Corr obs vs pred = 0.44122221745754225\n"

     ]

    },

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAApZUlEQVR4nO3de7wcdX3/8deb5AABlIBEJYdLqCKIUAikIMLPyq2AWokgBURb2lpqlRatP34G6U+oj1pj+VlrRWsp3hAr2AqRChrRoCgV4URACBehSCAHlECIXBI0hM/vj5lD9uzZnb2c2Z2Z3ffz8TiPszs7Z+a7s2fnM9/P9zKKCMzMzJrZrOgCmJlZuTlQmJlZJgcKMzPL5EBhZmaZHCjMzCyTA4WZmWVyoLCBIOk8SZcUXY5OSDpN0g+LLke3ao+5pF0kPSVpRh/2e7+kI3u9H9vEgWKIVekLl55Ub5O0TtIvJP2LpNlFl6vsJH1P0jPpSfxRSZdL2jHv/UTEAxGxTURsbFGe10lalff+022/Q9KdkraoWfYiSY9IOqYX+xwWDhRWepLeB3wUOAvYFng1sCtwjaTN+1iOmf3aV87OiIhtgFcAs4GP169Q4ff2vIi4CBgHPliz+J+AqyPiW4UUakA4UNgUkraT9A1JqyU9nj7eqeb10yTdJ+lJST+XdGq6/OWSvi/pV+nV62U1f/MaSTelr90k6TVtluWFwN8CfxkR34qIDRFxP/AHwDzgbTWrbynpsrRcP5G0b8123i9pPH3tbklHpMs3k7RI0v9IekzSVyVtn742T1JI+lNJDwDLJH1T0hl1ZbxV0vHp4z0lXSNpTbqfP6hZ70WSrpT0hKQbgZdlvO+m+1Hi4+mV8hNpTWvvVscyItYAXwP2Trd3f3pcfgo8LWmmpFdL+m9Ja9P9va5m/7uln++Tkq4Bdqh5beJYzUyfby/p85IeSv+HlkjaGvgmMDet4TwlaW7WZ5Bu6+2SVqavndPibb4DeJek/SQdDRwBvLfVsbEWIsI/Q/oD3A8c2WD5i4ATgK2AFwD/ASxJX9saeALYI32+I/Cq9PFXgHNILkC2BA5Nl28PPA68HZgJnJI+f1H6+iLgG03KeAzwLDCzwWtfBL6SPj4P2AC8BRgB/jfw8/TxHsCDwNx03XnAy9LHZwI3ADsBWwD/WrPNeUAAF6fvexbwh8D1NWXYC1ib/u3W6X7+OH2f84FHgb3SdS8FvpqutzfJ1e8Pm7zvrP0cDSwnqR0IeCWwY5PtfA94R/p4B2AZ8KWaz/8WYOf0vY0CjwGvTz/Do9Lnc9L1fwT8Y1qG1wJPApfUHauZ6fOrgMuA7dLP4HfT5a8DVtWVMesz2At4Kt3fFun+n6XB/23N9v4S+En6+S8s+ns2CD+FF8A/BX74TQJFg/X2Ax5PH2+dnrBOAGbVrXcxcCGwU93ytwM31i37EXBaG/t+G/CLJq8tBq5JH58H3FDz2mbAw8D/Al4OPAIcCYzUbeNO4Iia5zuSBJyZNSe/36p5/QXA08Cu6fMPA59LH58E/KBu+/8KnAvMSLe7Z81rf0/zQJG1n8OBn5Gk4DZrcfy+B6xLP7Nx4MtsOvHfD/xJzbrvJw0iNcuWAn8E7JKeoLeuee3faRAo0mP4HLBdg/K8jqmBIusz+CBwac1rWwO/ITtQCPgxcEXR37FB+XHqyaaQtJWkf02r+08A1wGzJc2IiKdJTojvBB6WdJWkPdM//T8kX9IbJa2Q9Cfp8rnAyrrdrCS5gm3lUWCHJjn0HdPXJzw48SAingNWkdQi7gXeQxJMHpF0qaS56aq7AlekqZa1JCetjcBLmmz3SZKr5ZPTRaeQnHwntnXQxLbS7Z0KvBSYQ3Lie35bTD0mz8vaT0QsAy4APpW+nwvTFF0zfxURsyNiNCJOjYjVjd5bWv4T68p/KMlxnktysfB0G+XfGVgTEY9nlKlW1mcwl8nH/2mSWk5TkUSLO4EVbe7fWnCgsEbeR5KuOSgiXkhS7YckCBARSyPiKJITyF3Av6XLfxERfxYRc4E/Bz4t6eXAQyQng1q7kFzhtvIj4NfA8bULJW0DHAt8t2bxzjWvb0aSyngoLdu/R8ShaTmCpHEckpPQsemJdOJny4ioLVv9FMtfAU6RdDBJiu3amm19v25b20TEXwCrSa7Id67Zzi4t3nuz/RAR/xwRB5CkZl5B0tDfjdr39iBJjaK2/FtHxGKS2tl2aTtDq/I/CGyvxr3SGk1XnfUZPMzkz3UrktSo9ZEDhY1I2rLmZyZJ2mM9sDZtVDx3YmVJL5F0XHrC+DVJ/vi59LUTtanR+3GSk8JzwNXAKyS9NW0wPYnkBPeNVoWLiF+RNGZ/UtIxkkYkzSPJ9a8CvlSz+gFpY+9MkhrEr4EbJO0h6XAl3SafSd/bc+nffAb4sKRd0/cwR9JxLYp1NUnA+RBwWVp7IX0/r0gbX0fSn9+R9MpIuo1eDpyX1tj2IknpdLyfdJsHSRohSU89U/N+puMS4PclHS1pRvr/8DpJO0XESmAM+FtJm0s6FPj9RhuJiIdJGq0/raRjxIikiYuNXwIvkrRtzZ9kfQb/CbxR0qFKerh9CJ+3+s4H3K4mOXFO/JxH0qVwFkla5wagtmvhZsBfk1yprwF+F/iL9LXfAX4s6SngSuDMiLgvIh4D3khSU3mMJEX1xoh4FEDSByR9s1kBI+IfgA8A/4+kIf3HJFehR0TEr2tW/TpJWmyi4fz4iNhA0gi6OH0/vwBeDJyd/s0n0rJ+W9KT6fs9KOuApfu8nKTN499rlj8J/B5JuuihdF8fTfcPcAawTbr8C8Dnu9kP8EKSWtzjJOmfx4Dzs7bVjoh4EDiO5FivJjnGZ7HpPPFWkmOzhuTi4eKMzb2dpJ3hLpL2ofek+7iLpKZ0X5pqmkvGZxARK4B3k7z/h9P33JNxGNacknSemZlZY65RmJlZJgcKMzPL5EBhZmaZHCjMzCxT5ScCa2SHHXaIefPmFV0MM7PKWL58+aMRMafRa4UGCkmfI+k2+UhETJnULJ2Q7Oskc7YAXB4RH2q13Xnz5jE2NpZjSc3MBpukpjMFFF2j+ALJVARZ/bF/EBFv7E9xzMysXqFtFBFxHcngHTMzK6kqNGYfrGRe/G9KelWzlSSdLmlM0tjq1aubrWZmZh0qe6D4Cck0y/sCnwSWNFsxIi6MiAURsWDOnIbtMWZm1oVSB4qIeCIinkofX00ygd0OLf7MzMxyVHRjdiZJLwV+GREh6UCSwJY5F721b8nN45y/9G4eWrueubNncdbRe7Bwfju3iDCzYVJ099ivkNzxagdJq0hmpBwBiIjPkNzW8i8kPUsys+nJ4VkMc7Hk5nHOvvw21m/YCMD42vWcffltAA4WZjZJoYEiIk5p8foFJN1nLWfnL737+SAxYf2GjZy/9G4HCjObpNRtFNY7D61d39FyMxteDhRDau7sWR0tN7Ph5UAxpM46eg9mjcyYtGzWyAzOOnqPgkpkZmVV6l5P1jsT7RDu9WRmrThQDLGF80cdGMysJQcKM+uKx+EMDwcKM+uYx+EMFzdmm1nHssbh2OBxoDCzjnkcznBxoDCzjnkcznBxoDCzjnkcznBxY7aZdczjcIaLA4WZdcXjcIaHU09mZpbJgcLMzDI5UJiZWSYHCjMzy+RAYWZmmRwozMwskwOFmZllcqAwM7NMDhRmZpbJgcLMzDJ5Cg+znPiObzaoHCjMcuA7vtkgKzT1JOlzkh6RdHuT1yXpnyXdK+mnkvbvdxnN2uE7vtkgK7pG8QXgAuDiJq8fC+ye/hwE/Ev626xUynbHN6fBLE+F1igi4jpgTcYqxwEXR+IGYLakHftTOrP2lemObxNpsPG16wk2pcGW3Dze97LYYCh7r6dR4MGa56vSZVNIOl3SmKSx1atX96VwZhPKdMc3p8Esb0WnnnITERcCFwIsWLAgCi6ODZky3fGtiDSYU12DreyBYhzYueb5Tukys6716qRWlju+zZ09i/EGQaFXaTD3+Bp8ZU89XQn8Ydr76dXAryLi4aILZdU1DPn7fqfBnOoafIXWKCR9BXgdsIOkVcC5wAhARHwGuBp4PXAvsA7442JKaoMi66TW6dVvWdMt/U6Dla3Hl+Wv0EAREae0eD2Ad/epODYE8jqplT3d0s80WL9TXZ0qa0CvkrK3UZjlKq+TWp41k6o76+g9JgVNKK7HF0wODNvOGuHp3zzLho1J/5ayBfSqKHsbhVmu8srfO92yycL5o3zk+H0YnT0LAaOzZ/GR4/cp5ERc3wa1dv2G54PEhHbaT5bcPM4hi5ex26KrOGTxsoFqw+qGaxSWadCq7Xnl78uebum3svT4alTTayQroJc9rVgEBwpralC/MHmc1MqWbrFEuzW6rIDutOJUTj1ZU+722FyZ0i22STs1ulYB3WnFqVyjsKb8hclWlnSLbdKopjeymdhmy5msXbehrVSj04pTOVBYU/7CWKeKbtPKow3KacWpHCisKX9hulf0CbMIZWnTmm5Nr0zzdpWFA4U15S9Md8pywuy3QWoEdlpxMgcKy+QvTOcG6YTZCbdpDS73ejLLWbMT4/ja9QM9cKtXN2/y4LfiOVCY5SzrxDhoM9XW6sWstcMw228VOFCY5azRCXPCII9D6cXYEo/lKQe3UZjlbOLE+J7Lbmn4ers5+yr2nMq7TasM7R5V/Bzy5hqFWQ8snD/K6DRy9k65JHrV7tEufw4JBwqzHplOzt4pl0S/79ZXz59Dwqknsx6ZzjiUMqRc+iUrtVP0WJ5h+hyyOFDYQClbPrnbnP2wTJ/SzuDEIsfyDMvn0IpTTzYwBimfXHTKpV/KntqpyufQ67EmDhQ2MMp+0unEsExjXvbUThU+h35cIDn1ZAOj7CedTg3D9Cmztxrh8XUbGi4vi7J/Dv2YMsY1ChsYRXeltM5FdLbcpurHBZIDhQ2MquSTbZNfrZ9am8hablP14wLJgcIGRm0+GWCG9HwVvIwN2p7szrXAPPTjAsmBwgbKwvmjz39xNqb5izL2fuqmAXIQA4trgdPXjwZ3RYHJQEnHAJ8AZgAXRcTiutdPA84HJr4RF0TERa22u2DBghgbG8u5tFYVhyxe1rDv++jsWVy/6PACSjRVp2WsH28AyQm1bD1wulG2sS/DStLyiFjQ6LXCej1JmgF8CjgKWAXcJOnKiLijbtXLIuKMvhfQKqsKvZ86LWOZboaU94m97L2KqqDXwbbI1NOBwL0RcV9E/Aa4FDiuwPLYgKhC3rvTMpYl+A3SoMZB0Y/PpMhAMQo8WPN8Vbqs3gmSfirpPyXt3J+iWZVVIe/daRn7HfyatYcM0qDGQdGPz6Tsjdn/BcyLiN8GrgG+2GxFSadLGpM0tnr16r4V0MqnCqNpOy1jP4Nf1hVqWWo2tkk/PpMiR2aPA7U1hJ3Y1GgNQEQ8VvP0IuAfmm0sIi4ELoSkMTu/YloV5Z337kUOuJMyTmcW1U7LnnWF6knyyqcfn0mRgeImYHdJu5EEiJOBt9auIGnHiHg4ffom4M7+FtGsvRlO+6Gb4NdN2bOuUD9+0n4Ne1+VKa03bM46eo+efyaFpZ4i4lngDGApSQD4akSskPQhSW9KV/srSSsk3Qr8FXBaMaW1YTadHHDRYx+6KXtWe0gV0nrDph+fSaGTAkbE1cDVdcs+WPP4bODsfpfLrFa3OeAy1ES6KXujK9SRzcS63zzLbouu8liHEup1F+OyN2abFa7bHkdF9xBacvM4m0kNX8sqe/0V6uxZIyB4fN0Gd4kdUg4UNnQ6TQd12+OoyB5CE7WZjQ1mXmin7Avnj3L9osP5+eI3sPUWM9mwcfJ23CV2uPh+FDZUukkHddvjqMgeQo1qM5BMlNhp/tpdYs2BwoZKt1NhdJMD7kdvlGaancSfi+j4fZSxS6znh+ovp55sqPTz6rjIHkJ5juQu20h3TyPSf65RWKVM90qy31fHRU14l2dtZjqD/XqhTBMkDgsHCquMPLqbFpkOytKLGVkhv5N7mWZ4dZtJ/zlQWGXkcSVZtqtjaB4Ax1au4dq7Vnddzn6d3PvdXlDGNpNB50BREm6cay2vK8l+Xx23+mybBcAv3/AAE51Si5o2pJUiBhV2Uiv09yofbswuATfOtacK95mo185n2yzQ1Y+AqB+7UPT0IFDMoMJ2Own4e5WfzBqFpO2zXo+INfkWZzi5ca49ZW1fyNLOZ9ssldLIRFApw/QgteVpd3le2qkV+nuVn1Y1iuXAWPp7NfAz4J708fLeFm14uHGuPVWckK6dz7ZR99NmJmpPvbiS76aGUuZanr9X+cmsUUTEbgCS/g24Ip3ED0nHAgt7Xroh4ca59pWp90072vlsJ97P3/7XCh5ft6HptmprT52cBNvJ03dbQ2m3lldEW4G/V/lpt43i1RNBAiAivgm8pjdFGj5lG9Bk+Wn3s104f5StNm9+3VZfe2r3Sr7dPH23NZRGtbwTDhjl/KV3P18z+Zslt3XcVpBH+4u/V/lpt9fTQ5L+BrgkfX4q8FBvijR8ythlsyrK3qulk8+2WS1BwPWLDp+07LA953DJDQ9MWfewPedMet5unn46aZraWl6jmklt762sMkzIq/3F36v8tBsoTgHOBa4g6YxxXbrMclK1lEoZlKVBt5WJz3YiqL33sls4f+ndU05anaRKrr2r8X3h65e3GwDyStM0CkzN7kvcrGx5NkL7e5WPtlJPEbEmIs4EDo2I/SPiPe7xZEUr+n4PnWiUAnrvZbcwrya10kmqpJMA0Ej98rzSNJ00FDcrmxuhy6etQCHpNZLuIL1ntaR9JX26pyUza6FKJ5SsK+3amlC7vbryDgB59ShrVq762ydlBaEy96QaVu2mnj4OHA1cCRARt0p6bc9KZdaGfvVqyaMdpFXwmqgJXb/o8La23W5vo07y9HmkaZqV64QDRtuejqSK42UGXdtTeETEg5p8W8Wpd0Uxy1nWSbofJ5S82kHaGVTXSU2o3wGgF+Xq5TYsX+0GigclvQYISSPAmaRpKLNWur0ib3WS7scJJa+G1UZBrd6gpFbyCExuhC6XdgPFO4FPAKPAOPBt4F29KpQNjulckbdzkm51Qplu2ijPiQgheU/ja9cjJvcG6rQmVJUeXzYY2h1wt0dEnBoRL4mIF0fE24BX9rJgNhim0zNpuifpPCaFy7NhdeH8Ua5fdDj3L34DHz9pv2k1HFepx5dVX7s1ik8C+7exzGyS6Zzsp9tYnUfaqFftINNNrVSpx5dVX6vZYw8mmapjjqS/rnnphUB7s5jZUJvOyX66J+k8TqbttoP45j02yFrVKDYHtknXe0HN8ieAt0x355KOIWn7mAFcFBGL617fArgYOAB4DDgpIu6f7n6tf6Zzsp9uY3VeJ9N22kE6aS/II6j0oqZT9ulQrDiKaDbAvmYladeIWJnrjqUZJNOWHwWsAm4CTomIO2rWeRfw2xHxTkknA2+OiJNabXvBggUxNjaWZ3FtGoo6AdWfwCE5meY9Nfkhi5c1DEijs2dNmaMpzzLleVz7daysvCQtj4gFjV5rt43iIkknRsTadIPbAZdGxNHTKNeBwL0RcV+6zUuB44A7atY5DjgvffyfwAWSFO1EN+ubViesoro6NquRQHJyzytwdZLiKus8Rr7Jj2VpN1DsMBEkACLicUkvnua+R4EHa56vAg5qtk5EPCvpV8CLgEenuW/LSdW6aY6tXMPXlo93Xd5GQbGTFFdZG6HLWi4rh3a7xz4naZeJJ5J2pfmkkIWQdLqkMUljq1c3nlmzKGW4t3GvlLmbZqPusV++4YGuy9usu+1he85pe0K9ss5jVNZyWTm0GyjOAX4o6UuSLiGZZvzsae57HNi55vlO6bKG60iaCWxL0qg9RURcGBELImLBnDlzGq1SiEG/wXuZr0TzmPJ6wpKbx3nfV29tGGSuvWt12xPqlfVmOmUtl5VDW6mniPiWpP2BV6eL3hMR003/3ATsLmk3koBwMvDWunWuBP4I+BFJL6tlVWufGPTcb5m7aeYx5TVsCvYbm/zrPbR2fdvtBWWZx6hRCu0jx+9TeLmsnFqNo9gzIu5KgwRsuqvdLpJ2iYifdLvjtM3hDGApSffYz0XECkkfAsYi4krgs8CXJN0LrCEJJpVS5ivuPJR5ps9mQazT6TMaBfv6/XSi6HmMmrUrfeT4fab00jKD1jWK9wF/BnyswWsBTOu/Kr0P99V1yz5Y8/gZ4MTp7KNoZb7izkOZr5DzmPIasoN6WYJiJwa9lmv5a2scRdWUaRzFMPVP76Rff9a6nY4PyDrGMP0g1mycxAyJj/3BvoV8jtMZQ7HboqsattUI+PniN+RaTquOrsdRSDo+6/WIuHw6BRsGZbni7rVOuslmrQt03N026wq53RsBZWlWM+lXsK8PCoftOWdaXXwHvZbbikegdy6zRiHp8+nDF5PM+bQsfX4Y8N8R8cbeFq87ZapRDItORidnrQu0vZ0J/bhCLtPo8vo2lglZx6jVNge1lltvmN97K13XKCLij9MNfBvYKyIeTp/vCHwh53JahXXSaN9NA3/Wa+1cIU831VVUA3SeXXwnDEsttxG3z3Sn3ZHZO08EidQvgV2arWzDp5N0Rqt1O02LtOp5lXeqq5/y6uJbr+ieV0UZ9F6IvdLugLvvSloq6TRJpwFXAd/pXbGsajoZsJW1bjcDvxbOH80c8JZ1FVnmkeXQ/OSvuudV7H1VBI9A7067A+7OkPRm4LXpogsj4oreFcuqppN0RjvrdpoWybpCzjvVlSXvtoy8uvj2soxVUuZxP2XWdvfYdH6n3SPiO5K2AmZExJM9LV2X3JhttfJuPG+mVw2lnk48X2ULlGUpT1Zjdrv3o/gz4HRg+4h4maTdgc9ExBH5FjUfVQoUZfknGWStxlnkdeLspOdXUapQxmFSpsCdx/0o3k1y/4gfA0TEPTlMMz6UagPD7K1GeOqZZ9nwXBKsy9aQOih6kepqpAoNpVUo4zCpSi+sdgPFryPiN1LShJbO5Dp4Q7p7rP7q4fF1G6as069/kmGryWS1YeTVA6gKA9m2nTXC2vVT/+/KVMZhUpXA3W6g+L6kDwCzJB0FvAv4r94VazC1mlxuQq//Sfp1s6FhC0bNGkoP23NOrnfU69aSm8d5+jfPTlk+spncmFuQKlxcQPvdY98PrAZuA/6cZCK/v+lVoQZVuwGg1/8k/egSOuj34WikUTfdEw4Y5WvLx0txHM5fejcbNk5NBGyz5cyBDuBlVpX7gLSsUUiaAayIiD2Bf+t9kQZXs6uHWv34J+lHdbcquVfIt+ZTn8Y6ZPGyUhyHJTePN/3fa5QCrf/bYaoZ9lNVRsm3DBQRsVHS3en9Jx7oR6EGVaPUxMgMsfXmM/nV+g19+yfpR3W3yNxrp7PY9jINV4Yc9MR7bGaG6ofvTf3bso5cHwRVGCXfbhvFdsAKSTcCT08sjIg39aRUA6osVw/9GHRUVO610xNbr2s+ZchBt2oba3bnvmZ/W9aaofVOu4Hi//a0FEOkDFcP/QhYRY2A7fTE1usr/jKMBG71XkYzglYZakRWvFb3o9gSeCfwcpKG7M9GxNRuE1Y5vQ5YRdWeOj2x9fqKvwy1yKy2sVZBqww1IiteqxrFF4ENwA+AY4G9gDN7XSgbDL0IRs3aHyaWN0uiNDux9eOKv+haZKP3CDB71gjnvelVmWUrQ43IitcqUOwVEfsASPoscGPvi2TWWLP2h7GVaybd8a1e1omtDFf8vTad9zgMx8daa3WHu59ExP7NnpdVleZ6queuiM1l3bu6WYPsqI+hWVumM9fTvpKemNgOycjsJ9LHEREvzLGcQ89dEbM1a2doFiQEnujOLAetboU6I+t1y5e7ImZr1rDarEZR1QZX1yqtbNqdwsP6YBC6Ii65eZxDFi9jt0VXccjiZblOVdFsuoNTDtq5EtMgtGMYpz6x8mt3HIX1QdW7IvY6dZbVsLpg1+0H4iq8zLVK13SGlwNFiVS9K2I/TnLNupoW3QU1L2WtVbr9bLgVknqStL2kayTdk/7ersl6GyXdkv5c2e9y9luj2UerdIvKsp7kqqRZ7bHoWmU/Zhy28iqqRrEI+G5ELJa0KH3+/gbrrY+I/fpasoJV+cq46qmzMihrrdIXAcOtqMbs40hGfZP+XlhQOSxHVZlbv8zKWqssa03H+iNzwF3PdiqtjYjZ6WMBj088r1vvWeAW4FlgcUQsydjm6cDpALvssssBK1euzL3c1pobPPNXhmNa30YxYbutRjj397OnAbFqyBpw17NAIek7wEsbvHQO8MXawCDp8YiY0k4haTQixiX9FrAMOCIi/qfVvqs8MtusVqMT9KyRGYXUMpbcPM55V66Ycs/tosrTa2UI0P2UFSh6lnqKiCMjYu8GP18Hfilpx7RwOwKPNNnGePr7PuB7wPxeldesjMrUiLxw/ihbbzG1WXMQG7U9nmWyohqzrwT+CFic/v56/QppT6h1EfFrSTsAhwD/0NdSWil1e6VXxSvEvBuRp3sMhqVRu8zjWYpQVGP2YuAoSfcAR6bPkbRA0kXpOq8ExiTdClxL0kZxRyGltdLo9kqvqleIeTYi53EMhqVRe1gCYrsKCRQR8VhEHBERu6cpqjXp8rGIeEf6+L8jYp+I2Df9/dkiymrl0m0qpkwpnE7k2ZOsnWPQagqWYenZNiwBsV2e68kqpdsrvapeIebZXbbVMWinxlHW7rt5G5aA2C5P4WGV0u2gvioPBsxrEGarY9BuXr7Kg0Lb5Rs2TeYahVVKt1d6vkJsfQyqWuvqlYXzR7l+0eH8fPEbuH7R4UMbJMA1CquYbq/0hvkKsban07azRthyZDPWrtsw5RhUudZlvVXIyOxe84C7/qlil9Nh0smAvTIN7rP+K2TAnQ2+qnY5HSad9PYaloZq65xTT9Y1D0oqv07bHYahodo65xqFdc2Nn+Xn8QCWBwcK65pPQuXn3l6WBwcK65pPQuXndgfLg9soSq7MvYqGuctplbjdwabLgaLEqnBDe5+EzAafU08lVtWJ7MxssDhQlJh7FZlZGTj1VGLDOqVCmdtlzIaRaxQlNoy9ijza26x8HChKbBi7Nrpdxqx8nHoquWHrVeR2GbPycY3CSsWjvc3Kx4HCSmUY22XMys6pJysVj/Y2Kx8HCiudYWuXMSs7p57MzCyTA4WZmWVy6skqxyO3zfqrkBqFpBMlrZD0nKSGN/NO1ztG0t2S7pW0qJ9ltHLyyG2z/isq9XQ7cDxwXbMVJM0APgUcC+wFnCJpr14VaMnN4xyyeBm7LbqKQxYv84mnpDxy26z/Ckk9RcSdAJKyVjsQuDci7kvXvRQ4Drgj7/JU4b4PlhiWkdtOr1mZlLkxexR4sOb5qnRZQ5JOlzQmaWz16tUd7chXqdUxDCO3nV6zsulZoJD0HUm3N/g5rhf7i4gLI2JBRCyYM2dOR3/b66tUp7XyMwwjt33hYmXTs9RTRBw5zU2MAzvXPN8pXZa7Xt73wWmtfA3DyO1hSa9ZdZS5e+xNwO6SdiMJECcDb+3Fjs46eo9JJ3PI7yo16+pwkE5u/TToI7eH9YZVVl5FdY99s6RVwMHAVZKWpsvnSroaICKeBc4AlgJ3Al+NiBW9KE+v7vuw5Obxhl948NXhIMorxTgM6TWrFkVE0WXI3YIFC2JsbKzQMtSnnOqNzp7F9YsO73OprFcafd6zRmZ0fcHhXk/Wb5KWR0TDcW1lTj1VWqOU0wRfHQ6evFOMg55es2pxoOiRrNTSoN/OtEryunJ3A7QNsjKPo6i0Zg2Po7NnOUiURJ7jFYZhfIcNLweKHnGDZP7yHo+S53gFf942yJx66pFh6O/fT70Yj5Jnusiftw0yB4oecoNkfnoxHiXv8Qr+vG1QOVAUrFljqrtHTtaLxuJeDrQ0GyQOFAVqlk4ZW7mGry0f97QfNXoxWtnpIrP2eMBdgQ5ZvKzhyW+GxMYGn8swD9LLe0CbmU3mAXcl1Sxt0ihIZK0/DHz1b1YcB4oCNUunNKtRDHuffDcWmxXD4ygK1Kzv/SkH7ew++WZWGq5RFCgrnbJg1+2dZqkw91qzQeLGbLOcueHdqiirMdupJ7Oc+VamNmgcKMxy5plkbdA4UJjlzDPJ2qBxoDDLmWeStUHjXk9mOfPgQBs0DhRmPeDBgTZIHCgsk8cDmJkDhTXVi5sFmVn1uDHbmvJ4ADMDBwrL4PEAZgZOPbVtGHP1vbhZkJlVTyE1CkknSloh6TlJDecWSde7X9Jtkm6RVNjkTRO5+vG16wk25eqX3DxeVJH6wuMBzAyKSz3dDhwPXNfGuodFxH7NJqvqh2HN1S+cP8pHjt+H0dmzEMkd9jyxndnwKST1FBF3AkgqYvcdG+ZcvccDmFnZG7MD+Lak5ZJOz1pR0umSxiSNrV69OtdCeO4eMxtmPQsUkr4j6fYGP8d1sJlDI2J/4Fjg3ZJe22zFiLgwIhZExII5c+ZMu/y1nKs3s2HWs9RTRByZwzbG09+PSLoCOJD22jVy5bl7zGyYlbZ7rKStgc0i4sn08e8BHyqqPM7Vm9mwKqp77JslrQIOBq6StDRdPlfS1elqLwF+KOlW4Ebgqoj4VhHlNTMbZkX1eroCuKLB8oeA16eP7wP27XPRzMysTtl7PZmZWcEcKMzMLJMDhZmZZVJEFF2G3ElaDazs4S52AB7t4farxsdjMh+PyXw8Jivr8dg1IhoOQhvIQNFrksaKnHuqbHw8JvPxmMzHY7IqHg+nnszMLJMDhZmZZXKg6M6FRRegZHw8JvPxmMzHY7LKHQ+3UZiZWSbXKMzMLJMDhZmZZXKg6JKk8yXdJemnkq6QNLvoMhWp3fugDzJJx0i6W9K9khYVXZ6iSfqcpEck3V50WcpA0s6SrpV0R/pdObPoMrXLgaJ71wB7R8RvAz8Dzi64PEXr5D7oA0fSDOBTJDfZ2gs4RdJexZaqcF8Ajim6ECXyLPC+iNgLeDXJzdgq8T/iQNGliPh2RDybPr0B2KnI8hQtIu6MiLuLLkeBDgTujYj7IuI3wKVAJ3dzHDgRcR2wpuhylEVEPBwRP0kfPwncCVTiJjcOFPn4E+CbRRfCCjUKPFjzfBUVOQlY/0maB8wHflxwUdpS2jvclYGk7wAvbfDSORHx9XSdc0iqlF/uZ9mK0M7xMLNskrYBvga8JyKeKLo87XCgyNDqvt+STgPeCBwRQzAgJY/7oA+wcWDnmuc7pcvMnidphCRIfDkiLi+6PO1y6qlLko4B/g/wpohYV3R5rHA3AbtL2k3S5sDJwJUFl8lKRJKAzwJ3RsQ/Fl2eTjhQdO8C4AXANZJukfSZogtUpGb3QR8WaceGM4ClJI2UX42IFcWWqliSvgL8CNhD0ipJf1p0mQp2CPB24PD0nHGLpNcXXah2eAoPMzPL5BqFmZllcqAwM7NMDhRmZpbJgcLMzDI5UJiZWSYHCrMGJO0k6euS7pH0P5I+IWlzSadJuqDo8tWT9FTRZbDB5UBhVicdGHU5sCQidgdeAWwDfLhH+/MMCVZqDhRmUx0OPBMRnweIiI3Ae0kmf9wK2FnS99LaxrkAkraWdJWkWyXdLumkdPkBkr4vabmkpZJ2TJd/T9I/SRoDzpG0UtJmNdt6UNKIpJdJ+lb69z+QtGe6zm6SfiTpNkl/1+8DZMPFVzJmU70KWF67ICKekPQAyXfmQGBvYB1wk6SrgF2BhyLiDQCStk3n9fkkcFxErE6Dx4dJAg7A5hGxIF1/f+B3gWtJ5g9bGhEbJF0IvDMi7pF0EPBpkkD2CeBfIuJiSe/u3aEwc6Aw68Y1EfEYgKTLgUOBq4GPSfoo8I2I+IGkvUkCyjVJNosZwMM127ms7vFJJIHiZODT6SyjrwH+I/17gC3S34cAJ6SPvwR8NNd3aFbDgcJsqjuAt9QukPRCYBeSKeXr572JiPhZWit4PfB3kr4LXAGsiIiDm+zn6ZrHVwJ/L2l74ABgGbA1sDYi9mvy955/x/rCbRRmU30X2ErSH8Lztzn9GMmtPdcBR0naXtIsYCFwvaS5wLqIuAQ4H9gfuBuYI+ngdDsjkl7VaIcR8RTJDLSfIKmRbEzvVfBzSSemfy9J+6Z/cj1JzQPg1FzfvVkdBwqzOum9Rd4MnCjpHpJ7oj8DfCBd5UaSewr8FPhaRIwB+wA3SroFOBf4u/SWqG8BPirpVuAWklRSM5cBb2NySupU4E/Tv1/Bpturnklyz+Xb8J30rMc8e6yZmWVyjcLMzDI5UJiZWSYHCjMzy+RAYWZmmRwozMwskwOFmZllcqAwM7NM/x80shS/i82zwwAAAABJRU5ErkJggg==\n",

      "text/plain": [

       "<Figure size 432x288 with 1 Axes>"

      ]

     },

     "metadata": {

      "needs_background": "light"

     },

     "output_type": "display_data"

    }

   ],

   "source": [

    "# Standard penalized methods (lasso using scikit-learn)\n",

    "from sklearn.metrics import mean_squared_error\n",

    "\n",

    "# alpha is the regularization parameter\n",

    "lasso = linear_model.Lasso(alpha=0.01)\n",

    "lasso.fit(X_train, y_train)\n",

    "y_hat = lasso.predict(X_test)\n",

    "\n",

    "# mean squared error\n",

    "mse = mean_squared_error(y_test, y_hat)\n",

    "print('\\nMSE in prediction =',mse)\n",

    "\n",

    "# correlation btw predicted and observed\n",

    "corr = np.corrcoef(y_test,y_hat)[0,1]\n",

    "print('\\nCorr obs vs pred =',corr)\n",

    "\n",

    "# plot observed vs. predicted targets\n",

    "plt.title('Lasso: Observed vs Predicted Y')\n",

    "plt.ylabel('Predicted')\n",

    "plt.xlabel('Observed')\n",

    "plt.scatter(y_test, y_hat, marker='o')\n",

    "plt.show()\n",

    "\n",

    "# Exercises\n",

    "# - Implement an internal crossvalidation to optimize alpha\n",

    "# - try different sizes of most associated SNPs\n",

    "# - implement ridge regression instead of lasso"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 5,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "\n",

      "The \"oldie\" way:\n",

      "Model: \"sequential\"\n",

      "_________________________________________________________________\n",

      "Layer (type)                 Output Shape              Param #   \n",

      "=================================================================\n",

      "dense (Dense)                (None, 64)                10368     \n",

      "_________________________________________________________________\n",

      "activation (Activation)      (None, 64)                0         \n",

      "_________________________________________________________________\n",

      "dense_1 (Dense)              (None, 32)                2080      \n",

      "_________________________________________________________________\n",

      "activation_1 (Activation)    (None, 32)                0         \n",

      "_________________________________________________________________\n",

      "dense_2 (Dense)              (None, 1)                 33        \n",

      "=================================================================\n",

      "Total params: 12,481\n",

      "Trainable params: 12,481\n",

      "Non-trainable params: 0\n",

      "_________________________________________________________________\n"

     ]

    }

   ],

   "source": [

    "# Implements a standard fully connected network (MLP) for a quantitative target\n",

    "# Currently, there are two ways of specifying a model in keras\n",

    "# 1. 'old' ways\n",

    "from tensorflow.keras.models import Sequential\n",

    "from tensorflow.keras.layers import Dense, Activation\n",

    "\n",

    "# no. of SNPs (inut features) in data\n",

    "nSNP = X_train.shape[1] \n",

    "\n",

    "print('\\nThe \"oldie\" way:')\n",

    "# Instantiate\n",

    "model = Sequential()\n",

    "# Add first layer that contains 64 neurons\n",

    "model.add(Dense(64, input_dim=nSNP))\n",

    "model.add(Activation('relu'))\n",

    "# Add second layer (32 neurons)\n",

    "model.add(Dense(32))\n",

    "model.add(Activation('softplus'))\n",

    "# Last, output layer, it is a single neuron since variable to predict is a float scalar\n",

    "model.add(Dense(1))\n",

    "# list some properties\n",

    "model.summary()"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 6,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "\n",

      "The \"new\" sequential way:\n",

      "Model: \"model\"\n",

      "_________________________________________________________________\n",

      "Layer (type)                 Output Shape              Param #   \n",

      "=================================================================\n",

      "input_1 (InputLayer)         [(None, 161)]             0         \n",

      "_________________________________________________________________\n",

      "dense_3 (Dense)              (None, 64)                10368     \n",

      "_________________________________________________________________\n",

      "dense_4 (Dense)              (None, 32)                2080      \n",

      "_________________________________________________________________\n",

      "dense_5 (Dense)              (None, 1)                 33        \n",

      "=================================================================\n",

      "Total params: 12,481\n",

      "Trainable params: 12,481\n",

      "Non-trainable params: 0\n",

      "_________________________________________________________________\n"

     ]

    }

   ],

   "source": [

    "# 2. 'New' way: adding layers is done recursively\n",

    "# https://keras.io/getting_started/intro_to_keras_for_engineers/\n",

    "from tensorflow.keras import layers\n",

    "from tensorflow.keras import Model\n",

    "\n",

    "# define input, 'None' if input dimension can vary\n",

    "input = keras.Input(shape=(nSNP))\n",

    "x = layers.Dense(64, activation='relu')(input)\n",

    "x = layers.Dense(32, activation='softplus')(x)\n",

    "output = layers.Dense(1, activation='linear')(x)\n",

    "\n",

    "# instantiate\n",

    "model = Model(inputs=input, outputs=output)\n",

    "\n",

    "# summary\n",

    "print('\\nThe \"new\" sequential way:')\n",

    "model.summary()"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 7,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Epoch 1/20\n",

      "15/15 [==============================] - 0s 5ms/step - loss: 0.8945\n",

      "Epoch 2/20\n",

      "15/15 [==============================] - 0s 5ms/step - loss: 0.8381\n",

      "Epoch 3/20\n",

      "15/15 [==============================] - 0s 1ms/step - loss: 0.8624\n",

      "Epoch 4/20\n",

      "15/15 [==============================] - 0s 3ms/step - loss: 0.8408\n",

      "Epoch 5/20\n",

      "15/15 [==============================] - 0s 3ms/step - loss: 0.7929\n",

      "Epoch 6/20\n",

      "15/15 [==============================] - 0s 3ms/step - loss: 0.7480\n",

      "Epoch 7/20\n",

      "15/15 [==============================] - 0s 3ms/step - loss: 0.6971\n",

      "Epoch 8/20\n",

      "15/15 [==============================] - 0s 3ms/step - loss: 0.7321\n",

      "Epoch 9/20\n",

      "15/15 [==============================] - 0s 6ms/step - loss: 0.6327\n",

      "Epoch 10/20\n",

      "15/15 [==============================] - 0s 5ms/step - loss: 0.6795\n",

      "Epoch 11/20\n",

      "15/15 [==============================] - 0s 2ms/step - loss: 0.7086\n",

      "Epoch 12/20\n",

      "15/15 [==============================] - 0s 3ms/step - loss: 0.6861\n",

      "Epoch 13/20\n",

      "15/15 [==============================] - 0s 4ms/step - loss: 0.6824\n",

      "Epoch 14/20\n",

      "15/15 [==============================] - 0s 3ms/step - loss: 0.6399\n",

      "Epoch 15/20\n",

      "15/15 [==============================] - 0s 2ms/step - loss: 0.6339\n",

      "Epoch 16/20\n",

      "15/15 [==============================] - 0s 4ms/step - loss: 0.5835\n",

      "Epoch 17/20\n",

      "15/15 [==============================] - 0s 4ms/step - loss: 0.5787\n",

      "Epoch 18/20\n",

      "15/15 [==============================] - 0s 5ms/step - loss: 0.5875\n",

      "Epoch 19/20\n",

      "15/15 [==============================] - 0s 5ms/step - loss: 0.6990\n",

      "Epoch 20/20\n",

      "15/15 [==============================] - 0s 5ms/step - loss: 0.5846\n"

     ]

    },

    {

     "data": {

      "text/plain": [

       "<tensorflow.python.keras.callbacks.History at 0x7f0dfdb885b0>"

      ]

     },

     "execution_count": 7,

     "metadata": {},

     "output_type": "execute_result"

    }

   ],

   "source": [

    "# Once defined, a model needs to be compiled, trained and validated\n",

    "# Model Compiling (https://keras.io/models/sequential/) \n",

    "# Stochastic Gradient Descent (‘sgd’) as optimization algorithm\n",

    "# Mean Squared Error ('mse') as loss, ie, quantitative variable, regression\n",

    "model.compile(optimizer='SGD', loss='mse') \n",

    "\n",

    "# training: this can take a while!\n",

    "model.fit(X_train, y_train, epochs=20)"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 8,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "1/1 [==============================] - 0s 213ms/step - loss: 1.3725\n",

      "\n",

      "MSE in prediction = 1.372475028038025\n",

      "\n",

      "Corr obs vs pred = 0.32896922412031954\n"

     ]

    },

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAozklEQVR4nO3df7wcdX3v8debw0EOYAmUtMKBQBQEQYTAufww1grWBlEkQhERb6VaqLfSWq+Xayhef1S8xOa21FYt4o9SRQVaMWLRRjCIFqWSGCgGiKD8SA5U+ZWCJEoSPvePmRP2nOzO7jk7Oz9238/H4zzO7uzszHdmd+cz38/3+51RRGBmZtbKdmUXwMzMqs2BwszMMjlQmJlZJgcKMzPL5EBhZmaZHCjMzCyTA4XVjqTLJF1YdjmmQ9IHJF1edjlmqnGfS/otSWsKWm9I2r+IdVlrDhQ2iaT7JD0taY8p01elP9r90uctD9bpfE9J+oWkcUl/LWmow/VL0nmS7pa0UdIDki6S9JyuN67PpZ/dxnS//yz9jHbJez0R8d2IOLCD8pwl6d/yXn+67AslfWvKtBdKekLSob1Y5yBzoLBm7gXOmHiS/vB2muYyDouIXYBXAm8Czu7wfX8LnAP8PvBc4NXpMq6a5vq7Imn7IteXo5PS/X4EMAa8d+oMNd62Rh8CnifpbEhOMIBPAX8dEbeXWrI+5EBhzXye5EA94S3A52ayoIi4C/gu8OJ280o6APhj4MyI+H5EbI6I1cCpwAmSjm+YfQ9J10l6UtKNkvZNlyFJF0v6eXp2ebukF6evPUfS/0trKT+TdImkkfS1V0haJ+k9kv4T+AdJd0p6bUP5tpf0sKQj0ufHSPqepPWSbpP0ioZ556blelLSdcCkGtqU7W65Hkk7Srpc0qPpem6R9Jsd7Pdx4BsT+z2t5b1D0t3A3em010q6NV3u9yS9pKEM8yT9MC3/lcCODa+9QtK6huf7SLo6LfOjkj4m6UXAJcCxaQ1nfbvPIH39PEkPSXpQ0lsztu9XwFuBxZL2Ijm52A34cLt9Y9PnQGHN3Az8mqQXpSmjNwIzyq9LOhj4LWBV+vwTkj7RYvZXAusi4geNEyNibVqmVzVMPpPkrHIP4FbgC+n03wVeDrwQ2BV4A/Bo+tridPrhwP7AKPC+hmU+D9gd2JfkwPMlGmpWwALgkYj4oaRR4FrgwvQ9/wv4sqTZ6bxfBFam5fsQSbBtpeV60vftCuwD/DrwdmBjxrKA5OANnEi631MLgaOBgyXNAz4L/FG63E8C16QH8h2ApSQnDLsD/0QSrJutZwj4F+B+YD+SfXpFRNyZlvX7EbFLRMxK39LyM5B0Asl+fBVwAPA7WdsYEf8OXJaW88PAWyNiU/aesRmJCP/5b+sfcB/JD/S9wEXACcB1wPZAAPul810GXNhiGQE8ATwO/ITkYLpdB+t+L3Bzi9euAD7VsO4rGl7bBdhCcjA9HvgxcEzjOgEBTwEvaJh2LHBv+vgVwNPAjg2v7w88CeyUPv8C8L708XuAz08p4zKSA/scYDOwc8NrXwQub7FtWet5K/A94CUdfna/ANaTHLg/AYw0fCbHN8z798CHprx/DfDbJIH2QUANr31v4vNO99W6hn34MLB9k/KcBfzbND6DzwKLG157YVru/TO2eSTd1ovL/u30818/5CqtNz4PfAeYy8zSTkdExD3TfM8jwJ4tXtuTpO1kwtqJBxHxC0mPAXtFxHJJHwM+Duwr6WqSs9QdSdpZVibpbCA5cDU2sj8cEb9sWO49ku4ETpL0NeB1wLz05X2B0ySd1PD+YeAGYC/g8Yh4quG1+0kC2TbarOfz6fuukDSLpGZ3QbQ+c14YEde3eG1tw+N9gbdI+pOGaTukZQ9gPNIjcUP5m9kHuD8iNrd4vdFssj+DvUhqYe3WuVVEbJR0L7C6g/XbDDn1ZE1FxP0kB+YTgasLWu1yYB9JRzVOTNMoxwCNvVz2aXh9F5IUyYMAEfG3EXEkcDDJWel5JEFoI3BIRMxK/3aNpOF3QrNLKU+khU4G7mgIfmtJahSzGv52jojFwEPAbpJ2bljOnDbb3nQ9EbEpIj4YEQcDLwVey+T2o+lo3L61wIenlH+niPhSWv5RNRzNM8q/FpjTooF86v5s9xk8xORg2m6fWUEcKCzL20jSFU+1eH0obWyd+Nuhm5VFxI9JGkC/kDYUD0k6BPgycP2UM+UTJb0sXeeHSFJWayX9N0lHSxomSXP8EngmIp4h6RVzsaTfAJA0KmlBm2JdQdLu8T9I0kcTLiepASxIy7lj2si7dxpkVwAflLSDpJcBJ2276PbrkXScpEPTtoAngE3AM22W1YlPAW9P95Uk7SzpNZKeC3yfJHX2p5KGJZ0CHNViOT8gOcAvTpexo6T56Ws/A/ae+F508BlcBZwl6WBJOwHvz2E7LQcOFNZSRPwkIlZkzLKI5Axx4m95u2WmvVwuyZjlXODTJAfiXwD/CnybbRtTv0hyIHkMOBJ4czr910gORo+TpC4eBZakr70HuAe4WdITwPVA5niAiHiI5MD5UuDKhulrSc7+/5wkR7+WpOYy8Zt6E0nD8WNpOTPTd63WQ9LA/s8kQeJO4EaSdFRX0s/1bOBjJPvqHpI2BSLiaeCU9PljwOm0qFVGxBaSILg/8ACwLp0fku/DauA/JT2STmv5GUTEN4C/Sd93Dx18n6wYmpyGNDMzm8w1CjMzy+RAYWZmmRwozMwskwOFmZll6ssBd3vssUfst99+ZRfDzKw2Vq5c+UhEzG72Wl8Giv32248VK7J6dZqZWSNJLUfCO/VkZmaZHCjMzCyTA4WZmWVyoDAzs0wOFGZmlqkvez2ZWXNLV42zZNkaHly/kb1mjXDeggNZOG+07GJZxTlQmA2IpavGOf/q29m4aQsA4+s3cv7VtwM4WFgmp57MBsSSZWu2BokJGzdtYcmyNSWVyOrCgcJsQDy4fuO0pptNcKAwGxB7zRqZ1nSzCQ4UZgPivAUHMjI8NGnayPAQ5y3IvMmfmRuzzdrpl55CE2Xuh22xYjlQmGXot55CC+eN1rLcVi6nnswyuKeQmQOFWSb3FDJzoDDL5J5CZg4UZpncU8jMjdlmmdxTyMyBwqwt9xQaLP3SHTpPDhRmZql+6w6dFweKPtfp2ZHPosyyu0MP8u/BgaKPdXp25LOo6XNg7U/uDt2cez31sU4Hi3lQ2fRMBNbx9RsJng2sS1eNl10065K7QzfnQNHHOj078lnU9DiwbmvpqnHmL17O3EXXMn/x8toGTXeHbs6Boo91enbks6jpcWCdrJ9qWAvnjXLRKYcyOmsEAaOzRrjolEMHPq1YahuFpM8CrwV+HhEvbvK6gI8CJwIbgLMi4ofFlrK+zltw4KS2B2h+dtTpfJbYa9YI402CQt6BtS7tIP3WAOzu0Nsqu0ZxGXBCxuuvBg5I/84B/r6AMvWNTs+OfBY1PUWkJ+p0lu4aVv8rtUYREd+RtF/GLCcDn4uIAG6WNEvSnhHxUDElrL9Oz458FtW5IkZr1+ksvagalpWn6t1jR4G1Dc/XpdMcKKxUvQ6sdTpL7yR1OZ00Wrt565KS60RdtqXqgaJjks4hSU8xZ86ckktj1p06naW3q2FNZ5xOu3n7acxPnbal6oFiHNin4fne6bRtRMSlwKUAY2Nj0fui2aAo46wvjw4GU8t93EGzueGuh3uyHc1qWBPrbxbwWqXR2qXc6pSSa6dO21L1QHENcK6kK4Cjgf9y+4QVqayzvm7bQZqV+/KbH9j6eruz+m4D49T1N9MsjdYu5VanlFw7ddqWsrvHfgl4BbCHpHXA+4FhgIi4BPg6SdfYe0i6x/5BOSW1QVXmWV837SDNyj1Vs+3IKzB2sv5mabR2Kbc6peTaqdO2lNo9NiLOiIg9I2I4IvaOiM9ExCVpkCAS74iIF0TEoRGxoszy2uCp01lfo07LN3W+vEadt1t/qzRau67H/TRyuk7bUvXUk5WkLr0xeq3VWd92EnMXXVvZfdOq3M3ma5RXYMxa/2jGPmuXcmv2+nEHzWbJsjW868pbK/t5NFOnm2IpGaLQX8bGxmLFinpVPqp0YG6WXx4ZHhrIQXid5NqruG9mWu75i5c3PcCPzhrhpkXHd7X+Xuwnf1fzI2llRIw1e63skdlG9Ubh1v2id3leoG7qqPUhaZt5qrhvJsqdpdnBNK90SFGj/ev+Xa0Lp54qoGrd5HqVly+i1tSLXkqNjcpzF13bdJ4qtllMdCdtVUOYSepnuuvv9fe31X4fX7+R+YuXVzaVUzcOFBVQtQbTXvTGKKqbaa+Dbqt9s+vIcNfL7lazQDyT8Rh1upxLVltIlQew1Y1TTxVQtct896I3RlEpgl4H3fMWHMjwdtumn556enOpF+xrlb4E+vqCj82+q42chsqHaxQV0Oqs77iDZjN/8fLCG7h70RujqFpTr2pDjftieEhsemZyJ5BNW6LUEbVZgfimRcf3TWCYqvG72qpmUcW0YN04UFRAqy5/X145Xtp1YPJOPxQ1uCjve2s0S5m10uyAVFRvtqqlL4s08V1t1WOrigPY6sapp4pYOG+UmxYdz72LX8NNi47nhrse7qveHEUNLsq7t00nI4wnTD0gFdmbrZP0Zb/crrSVOg1gqxvXKCqq384QixxclGdtqNP93eyAVGRvtnY1qU6uylqVcTwzVacBbHXjQFFRdboOTKeyDuB5HqjyXFarz2G3nYbZaYftM9dRZLCfWPcHrlnN+o2bANhx+NmEQbvOBHW53HU7deqxVScOFBU1SPexzrPrbN7dcFt9Du8/6ZC2yysj2P9q8zNbHz++YdPWbc8KWlUbx2PV4zaKihqk+1jn2XU272643XwOzXLm4tnBYHm3EWRte1YbRp3TnP3a7lK17XKNImd5pj0GpRrdzYFq6v7uRRfJmX4OU7tuCpjoVNuL9E7Wfrz49MNb1lBbdS2tepqzTneIm44qbpdrFNPQLspX7ZpNdTHTAYfN9ve2Q+E6W1avTPRmG501wtTLb+bdiy1rP2bVjOraW6hfr/NUxe1yjaJDnUR553pnZqbtMc32d8CkM/dOl9VrRaR32u3HVjWjuvYWqnPKLEsVt8uBokOdBIEqfsB1MNN7DLTar0Fyxlylg14RDdvdHPDrmObsx56BUM3tcqDoUCdBoIofcF00Hqg6zdG22t/TvXdCEYrqxdbNAb9uYyn6tWdgFbfLbRQd6iSPXtdcb9V0mqOt0/6uei+2OravVX2fzlQVt8s1ig51EuXrmuutmk5qbxNnvxs3bWFIYktE5i02q6DK6Z1WwfndV90GVLcXUZX3aTeqtl0OFB3qNAhU7QOuumbpjnYpvKmpqS0RW4O29/3MtArOWyJK75pp5XOgmAYHgXy1aos49cjRSVfOhcm1N/cuy1/WGJRe79u6tY0MIgcK67lWB4JWB/wb7nqYi045tOl7lq4ar+x9B3p1wCviQNostdqoV/u2ioPLbFsOFBXTb2dXWQeCrLaIZrW3iWW1Umbvsl4d8Io6kE4s691X3caWmDo0sHf71rXDenCvpwqpY8+TdmZ6/aFOlzWh7N5OeYymbTbyv8hRugvnjfJXbzis0J5kHntUDw4UFVLFofvdyjoQTLd7a9bBo+zug90e8FqdJBSdZiu6a2bV7hdvzTn1VCH9eHaV1YNput2JswbYlZ2m6HawZauThImuvzNd7kwU2WmjioPLbFsOFBXSjyO7Z3r9oZksq0ydlq1VG1RW99SR4aFtgsiGpzezdNV46QFyJt679Ha+9O9r2RLBkMQxz9+N+x7d2Dftcv2o1NSTpBMkrZF0j6RFTV4/S9LDkm5N//6wjHIWpU4jjTuVZyqjiiNWp1O2rDaoVicDE8uZNTI8afrETYnq1n713qW3c/nND2ytJW2J4KafPMZxB83eer/4KnyeNpmiSbW2kBVLQ8CPgVcB64BbgDMi4o6Gec4CxiLi3Okse2xsLFasWJFjaYvTb72e7FnzFy9vmTprVSOZCDaHf/CbW29xOvW9072uVZnfsRec//WmqbQhiZ9cdGIhZbDmJK2MiLFmr5WZejoKuCcifgog6QrgZOCOzHf1ubzyww441dMqvTS+fiPvuvJWdh0ZZsfh7Vi/YdM2Y0eaBYmsZULz7wCUe3/sZkEia7pVQ5mBYhRY2/B8HXB0k/lOlfRyktrHuyJibZN5rEEVBzE5cGWPfg5g/cZNjAwPcfHph0/aN1m93nYdGWb+4uVNByY2+w48Z/vtSh230KpxfkitbjllVVD17rFfA/aLiJcA1wH/2GpGSedIWiFpxcMPP1xYAauoat1s+3F8yEw0a4OaqtnnlFVreOrpzU33a6vvwExqJnk64+h9pjXdqqHMQDEONH479k6nbRURj0bEr9KnnwaObLWwiLg0IsYiYmz27Nm5F7ZOqtbNttVB64NfW11KecoytcG7lamfU6uG7u0Em7ZMPjufCDTT/ayL6ll34cJDefMxc7bWIIYk3nzMHC5ceGgh6+9X7W7T3K0yU0+3AAdImksSIN4IvKlxBkl7RsRD6dPXAXcWW8R6qlo321YHrcc3bOppF88y0l3t1tnYBtWqcXvq59SqoTvrukyzdhrm8Q3b1h523mGIZ4JSuxhfuPBQB4YcFZFqLq1GERGbgXOBZSQB4KqIWC3pLyS9Lp3tTyWtlnQb8KfAWeWUtl667Wab99lJVoDqVTqsjHTXdNfZ6efUquvtaMao5lZtw8ND21W2i7HNTBGp5tK6x/ZSnbvH5mWmZ9NTz05gcjfNmZblz668telrAu5d/JoZLTdLVlfUXt0mdSbr7KbWk/VZvevKW2n2y+7V/rbyzF10bS6fdVW7x1oPzbSbbV5X85x6ABwZ3o6Nm57ZZr526bCZHkjLaKeZyTq76Q6ddQmUJcvWVCr9aL1TRKrZgcImyeMA2yxnOjwkhrcTm5559tynXTqsm9xr3j+eTgJWq3UGSW2jF20krQJNlS93Yvkq4rOuevdYK1geV/NsVivZtCXYZcftp5Ub7yb3muflUDpte8jq/lp0l+AqX+7E8lXEZ+0ahU2Sx9lJq9rH+g2bWPW+3+16OZ3UbqZ7ZdosnabjGtfZrGZR9A15fOvewdHrz9qBokB1GJ2cxwE2r7RPt8vJ68cznYA1sc5WDYx1vmS8DS4HioJU8bIarXR7gM0rZ9rtcvIKzDMJWFUby2LWDbdRFKRql9Xopbxypt0sJ89xFDNp7zhvwYEMD00efz08JDcmd6HXo4+tNdcoUr1OC+XVm6jqqasJeaV9yu7mO1GGiWVOa99PzT3135ClwtSpRt6PHCgo5kvYbSrCP5TpyXscxXQD1pJlayZ1BQbY9EwU2pjdT/IM/DZ9Tj1RTFqo2+6aeZZxEKrweXTz7UbVLsxYd96f5XKgoJgvYbd5+7zKOCiX/C77trJlB6p+4/1ZLgcKivsSLpw3yk2Ljp/RvYHzKuOgNKpPDcyz0rvHvevKWwupRZUdqPqN92e53EZBPS53kFcZO62Z1KnhvJWJdoVO2ne6vThfq/fWfR9WhfdnuRwoqMeXMK8ydtKoXseG83YH66yG0Hbbm7Xsdu+t6v6qI+/P8jhQpOrwJcyjjJ3UTOrWw6TdwbpdLapdOi5r2XXbV2YzkdlGIWn3rL+iCmn56aRRvW49TNod6Nu172Rtb7tl121fmc1EuxrFSpJhQgLmAI+nj2cBDwBze1k46412NZNOx3xUpR2j3cG6XS0qa3vbLduX6rBBkFmjiIi5EfF84HrgpIjYIyJ+HXgt8M0iCmjF66SHSZW62barMbSrRWVtb7tluzeODYJO2yiOiYizJ55ExDck/WWPymQl66ThPI/cfF41kk7aXbJqUe22N2vZdegIYdatju6ZLWkZ8F3g8nTSmcDLI2JBD8s2Y75ndu91e5/evO/N3cs0WFVSbNYb/nwTedwz+wzg/cBXSNosvpNOswHVbTtG3r2FetlrrQ494mxm6tgVvAwdBYqIeAx4p6SdI+KpHpfJaqCTdE/Wj9C9haydIs703b25Mx1dwkPSSyXdAdyZPj9M0id6WjKrtE662Wb9CH3tHstSVGcJn7B0ptPU08XAAuAagIi4TdLLe1Yqq4V2KZmsH+HFpx9e+cumWHmKOtN39+bOdHxRwIhYO2XSlqYzmqWyag153QVvUAzCpeEbFXWm7+7Nnem0RrFW0kuBkDQMvJM0DWXWSrt2jH5tJM47tz6IDa5Fnem7e3NnOg0Ubwc+CowC4ySD7f64V4Wy3ii6G+Ag/gh7cVBvlYb54NdW9+2+LPKKzv16wpKnTgPFgRFxZuMESfOBm/IvkuVpIjiMr9+IePa2zd0ewDoNOoP2I+xFbr1VuuXxDZtYumq8L/fvIJ5kVFmngeLvgCM6mDYtkk4gqakMAZ+OiMVTXn8O8DngSOBR4PSIuK+bdfZCVQfsTD27nTpAbqYHsEFMhXSqF7n1VmkYoOXnV9Xv5HQM2klGlbW7euyxkt4NzJb0Pxv+PkBycJ8xSUPAx4FXAwcDZ0g6eMpsbwMej4j9SXpefaSbdfZCla55NFWzs9upZnIAG5S75M1EL7r9ZqVbmn1+Vf5OWj216/W0A7ALSc3juQ1/TwC/1+W6jwLuiYifRsTTwBXAyVPmORn4x/TxPwOvlKQu15urKh80OwkCMzmAue95a73oRbNw3iizRoabvtbs86vyd9LqKTP1FBE3AjdKuiwi7s953aNAY5fbdcDRreaJiM2S/gv4deCRnMsyY1U+aGalLGDmB7B+7HueV6qmV7n1D7zukG0ad4e3Exue3szcRddOWk+Vv5NWT522UXxa0mkRsR5A0m7AFVW6KKCkc4BzAObMmVPYeqt80GzWc2SiQXu0x1drrZO821x6kVufGoB2HRnmqac38/iGTcDkMlf5O2n11OmAuz0mggRARDwO/EaX6x4H9ml4vnc6rek8krYHdiVp1N5GRFwaEWMRMTZ79uwui9a54w6azdRcWFUOms0GtV18+uHct/g13LTo+BkfzPptsFxdUjUL541y06LjuXfxa9j5Oduzacvk7gkTZfYgMstbpzWKZyTNiYgHACTty7adaKbrFuAASXNJAsIbgTdNmeca4C3A90naRJZHJ9dFL8jSVeN8eeX4pB0h4NQjq9Nbo1c9R/qpR0odUzVZZXbXUstbp4HiAuDfJN1Iciz8LdI0z0ylbQ7nAstIelB9NiJWS/oLYEVEXAN8Bvi8pHuAx0iCSWU0OxMN4Ia7Hi6nQDYjdUzVtCtzPwVyK1+nlxn/V0lHAMekk/4sIrpuUI6IrwNfnzLtfQ2Pfwmc1u16eqWOZ6KDpNMG6jq2udSxzFZfmYFC0kERcVcaJAAeTP/PSVNRP+xt8aqtjmeig2I6DdR1TNXUscxWX5m3QpX0qYg4W9INTV6OiDi+d0WbuaJuhZr37TwtP/MXL28axEdnjXDTokp+bc1KNeNboUbE2en/43pRsLrzWV3+8hrP4LSgWX7apZ5OyXo9Iq7Otzj140bD/OQ5nsFpwerph+tPDap24yhOSv/eRtID6cz079PAW3tbNBs0eY5n8FiCavH1p+otM1BExB9ExB8Aw8DBEXFqRJwKHJJOM8tNnumifhsUWHd1GdRozXU6jmKfiHio4fnPgOKuk2EDIe90kdOC1eE2o3rr9BIe35K0TNJZks4CrgWu712xbBA5XdS/enH5dStOR4EiIs4FLgEOS/8ujYg/6WXBbPA4XdS/fBJQb52mngB+CDwZEddL2knScyPiyV4VzAZTVdNF7rHTHXclr7eOAoWks0mu7bQ78AKS+0RcAryyd0Uzqwbf+jUfVT0JsPY6baN4BzCf5M52RMTddH+ZcbNacI8dG3SdBopfpbcrBbbeG6Iyl/s26yX32LFB12mguFHSnwMjkl4F/BPwtd4Vy6w63GPHBl2ngeI9wMPA7cAfkVwa/L29KpRZlbjHjg26to3ZkoaA1RFxEPCp3hfJrHNF9EZyjx0bdG0DRURskbSm8VaoZlVQZG8k99ixQdbpOIrdgNWSfgA8NTExIl7Xk1JVgPvNV19WbyR/Vmb56TRQ/J+elqJi3G++HppdFwrcG8ksb+3uR7Ej8HZgf5KG7M9ExOYiClYmn6lW39JV44jmfbTdG8ksX+1qFP8IbAK+C7waOBh4Z68LVTb3m6++JcvWNA0SAvdGwqlTy1e7QHFwRBwKIOkzwA96X6Ty+e5o7ZV9IGoVtAOnB506tby1G0exaeLBIKScJrjffLYq3K2sVdAedTD3JUcsd+0CxWGSnkj/ngReMvFY0hNFFLAMvtx1tiociBzMW3Pq1PKWmXqKiKGs1/uZ+823VoUDkQfBtebUqeVtOvejMAOqcyByMG/uvAUHTmqjANe2piq7ja1uOr3Wk9lWnaR9lq4aZ/7i5cxddC3zFy8vtP1i0Dl1mq0KbWx14xqFTVu7tE8Vet0M+hmja1uteZzU9JUSKCTtDlwJ7AfcB7whIh5vMt8WkoF+AA/08yVD6ibrQFT2D7EKgcqqqwptbHVTVuppEfCtiDgA+Fb6vJmNEXF4+ucgURNl/xCr0CvLqsv3F5m+sgLFySSjvkn/LyypHNYDZf8Qyw5UVm3uWj19ZQWK34yIh9LH/wn8Zov5dpS0QtLNkhYWUzTrVtk/xLIDlVWbG/unr2dtFJKuB57X5KULGp9EREhqdf/tfSNiXNLzgeWSbo+In7RY3znAOQBz5szpouTWrbLHOLh7qLXjxv7pUUSrY3QPVyqtAV4REQ9J2hP4dkRk/oolXQb8S0T8c7vlj42NxYoVK/IprNXSoPd6MpsuSSsjYqzZa2V1j70GeAuwOP3/1akzSNoN2BARv5K0BzAf+MtCS2m15TNGs/yUFSgWA1dJehtwP/AGAEljwNsj4g+BFwGflPQMSVvK4oi4o6TyWsFcIzCrjlICRUQ8CryyyfQVwB+mj78HHFpw0awCPA7CrFp8CQ+rHI+DMKsWBwqrHI+DMKsWBwqrHI+DMKsWBwqrnLIH7JnZZL56rFVO2QP2zGwyBwqrJI+DMKsOp57MzCyTA4WZmWVyoDAzs0wOFGZmlsmBwszMMjlQmJlZJgcKMzPL5EBhZmaZHCjMzCyTA4WZmWVyoDAzs0wOFGZmlsmBwszMMvnqsdb3lq4a9yXLzbrgQGF9bemqcc6/+vat9+AeX7+R86++HcDBwqxDTj1ZX1uybM3WIDFh46YtLFm2pqQSmdWPA4X1tQfXb5zWdDPblgOF9bW9Zo1Ma7qZbcuBwvraeQsOZGR4aNK0keEhzltwYEklMqsfN2ZbX5tosHavJ7OZc6Cwvrdw3qgDg1kXnHoyM7NMpQQKSadJWi3pGUljGfOdIGmNpHskLSqyjGZmligr9fQj4BTgk61mkDQEfBx4FbAOuEXSNRFxRzFFrC6PNDazIpUSKCLiTgBJWbMdBdwTET9N570COBkY6EDhkcZmVrQqt1GMAmsbnq9Lpw00jzQ2s6L1rEYh6XrgeU1euiAivtqD9Z0DnAMwZ86cvBdfGR5pbGZF61mgiIjf6XIR48A+Dc/3Tqe1Wt+lwKUAY2Nj0eW6K2uvWSOMNwkKHmlsZr1S5dTTLcABkuZK2gF4I3BNyWUqnUcam1nRyuoe+3pJ64BjgWslLUun7yXp6wARsRk4F1gG3AlcFRGryyhvlSycN8pFpxzK6KwRBIzOGuGiUw51Q7aZ9Ywi+i9LMzY2FitWrCi7GGZmtSFpZUQ0HdfmS3jkwOMazKyfOVB0yeMazKzfVbkxuxY8rsHM+p0DRZc8rsHM+p0DRZd8BzUz63cOFF3yuAYz63duzO6S76BmZv3OgSIHvoOamfUzp57MzCyTA4WZmWVyoDAzs0wOFGZmlsmBwszMMjlQmJlZJgcKMzPL5EBhZmaZHCjMzCyTA4WZmWVyoDAzs0wOFGZmlsmBwszMMvnqsSVZumrclyY3s1pwoCjB0lXjnH/17VvvtT2+fiPnX307gIOFmVWOU08lWLJszdYgMWHjpi0sWbampBKZmbXmQFGCB9dvnNZ0M7MyOVCUYK9ZI9OabmZWJgeKEpy34EBGhocmTRsZHuK8BQeWVCIzs9bcmF2CiQZr93oyszooJVBIOg34APAi4KiIWNFivvuAJ4EtwOaIGCuqjL22cN6oA4OZ1UJZNYofAacAn+xg3uMi4pEel8fMzFooJVBExJ0AkspYvZmZTUPVG7MD+KaklZLOKbswZmaDqGc1CknXA89r8tIFEfHVDhfzsogYl/QbwHWS7oqI77RY3znAOQBz5syZUZnNzGxbPQsUEfE7OSxjPP3/c0lfAY4CmgaKiLgUuBRgbGwsul23mZklKts9VtLOwHYR8WT6+HeBv+jkvStXrnxE0v09LWB7ewCD2gjvbR9cg7z9dd/2fVu9oIjiT74lvR74O2A2sB64NSIWSNoL+HREnCjp+cBX0rdsD3wxIj5ceGFnSNKKfurOOx3e9sHcdhjs7e/nbS+r19NXeDYINE5/EDgxffxT4LCCi2ZmZlNUvdeTmZmVzIGidy4tuwAl8rYPrkHe/r7d9lLaKMzMrD5cozAzs0wOFGZmlsmBoockLZF0l6T/kPQVSbPKLlNRJJ0mabWkZyT1ZZfBqSSdIGmNpHskLSq7PEWS9FlJP5f0o7LLUiRJ+0i6QdId6ff9nWWXqRccKHrrOuDFEfES4MfA+SWXp0gTVwhuOpK+30gaAj4OvBo4GDhD0sHllqpQlwEnlF2IEmwG3h0RBwPHAO/ox8/dgaKHIuKbEbE5fXozsHeZ5SlSRNwZEWvKLkeBjgLuiYifRsTTwBXAySWXqTDpNdgeK7scRYuIhyLih+njJ4E7gb670YwDRXHeCnyj7EJYz4wCaxuer6MPDxjWmqT9gHnAv5dclNxV9lpPddHJVXIlXUBSRf1CkWXrtZyuEGxWe5J2Ab4M/FlEPFF2efLmQNGldlfJlXQW8FrgldFng1byuEJwHxkH9ml4vnc6zfqcpGGSIPGFiLi67PL0glNPPSTpBOB/A6+LiA1ll8d66hbgAElzJe0AvBG4puQyWY8puU3nZ4A7I+Kvyy5PrzhQ9NbHgOeS3HTpVkmXlF2gokh6vaR1wLHAtZKWlV2mXko7LZwLLCNp0LwqIlaXW6riSPoS8H3gQEnrJL2t7DIVZD7w34Hj09/4rZJOLLtQefMlPMzMLJNrFGZmlsmBwszMMjlQmJlZJgcKMzPL5EBhZmaZHCjMmpC0t6SvSrpb0k8kfVTSDpLOkvSxsss3laRflF0G618OFGZTpIOorgaWRsQBwAuBXYAP92h9vkKCVZoDhdm2jgd+GRH/ABARW4B3kVzYcSdgH0nfTmsb7weQtLOkayXdJulHkk5Ppx8p6UZJKyUtk7RnOv3bkv5G0grgAkn3S9quYVlrJQ1LeoGkf03f/11JB6XzzJX0fUm3S7qw6B1kg8VnMmbbOgRY2TghIp6Q9ADJb+Yo4MXABuAWSdcC+wIPRsRrACTtml4D6O+AkyPi4TR4fJgk4ADsEBFj6fxHAL8N3EBybbBlEbFJ0qXA2yPibklHA58gCWQfBf4+Ij4n6R292xVmDhRmM3FdRDwKIOlq4GXA14G/kvQR4F8i4ruSXkwSUK5LslkMAQ81LOfKKY9PJwkUbwQ+kV6R9KXAP6XvB3hO+n8+cGr6+PPAR3LdQrMGDhRm27oD+L3GCZJ+DZhDcrn4qde9iYj4cVorOBG4UNK3gK8AqyPi2Bbrearh8TXA/5W0O3AksBzYGVgfEYe3eL+vv2OFcBuF2ba+Bewk6fdh621O/4rkdp8bgFdJ2l3SCLAQuEnSXsCGiLgcWAIcAawBZks6Nl3OsKRDmq0wIn5BcgXaj5LUSLak9zW4V9Jp6fsl6bD0LTeR1DwAzsx1682mcKAwmyK9b8jrgdMk3U1yv/NfAn+ezvIDkvsP/Afw5YhYARwK/EDSrcD7gQvTW6L+HvARSbcBt5Kkklq5Engzk1NSZwJvS9+/mmdvr/pOkvsz347vpGc95qvHmplZJtcozMwskwOFmZllcqAwM7NMDhRmZpbJgcLMzDI5UJiZWSYHCjMzy/T/AYRks80YcZT2AAAAAElFTkSuQmCC\n",

      "text/plain": [

       "<Figure size 432x288 with 1 Axes>"

      ]

     },

     "metadata": {

      "needs_background": "light"

     },

     "output_type": "display_data"

    }

   ],

   "source": [

    "# cross-validation: get predicted target values\n",

    "y_hat = model.predict(X_test, batch_size=128)\n",

    "\n",

    "mse_prediction = model.evaluate(X_test, y_test, batch_size=128)\n",

    "print('\\nMSE in prediction =',mse_prediction)\n",

    "\n",

    "# correlation btw predicted and observed\n",

    "corr = np.corrcoef(y_test,y_hat[:,0])[0,1]\n",

    "print('\\nCorr obs vs pred =',corr)\n",

    "\n",

    "# plot observed vs. predicted targets\n",

    "plt.title('MLP: Observed vs Predicted Y')\n",

    "plt.ylabel('Predicted')\n",

    "plt.xlabel('Observed')\n",

    "plt.scatter(y_test, y_hat, marker='o')\n",

    "plt.show()\n",

    "\n",

    "# Exercises\n",

    "# - Check predictions across environments (Y[0] is first environment, etc)\n",

    "# - Try to improve model with other activation functions and|or no. of neurons∫ "

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 10,

   "metadata": {},

   "outputs": [],

   "source": [

    "# Controlling overfit: regularization, dropout and early stopping\n",

    "\n",

    "# deletes current model if exists \n",

    "if 'model' in locals(): del model\n",

    "\n",

    "# define input\n",

    "input = keras.Input(shape=(nSNP))\n",

    "x = layers.Dense(64, activation='relu',\n",

    "                     kernel_regularizer=keras.regularizers.l2(0.01),\n",

    "                     activity_regularizer=keras.regularizers.l1(0.01))(input)\n",

    "x = layers.Dense(32, activation='softplus')(x)\n",

    "x = layers.Dropout(rate=0.2)(x)\n",

    "output = Dense(1, activation='linear')(x)\n",

    "\n",

    "# instantiate\n",

    "model = Model(inputs=input, outputs=output)\n",

    "\n",

    "# Model Compiling (https://keras.io/models/sequential/) \n",

    "model.compile(loss='mean_squared_error', optimizer='sgd')"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 11,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Epoch 1/20\n",

      "14/14 [==============================] - 1s 55ms/step - loss: 3.4939 - val_loss: 2.2004\n",

      "Epoch 2/20\n",

      "14/14 [==============================] - 0s 17ms/step - loss: 2.0722 - val_loss: 2.1860\n",

      "Epoch 3/20\n",

      "14/14 [==============================] - 0s 16ms/step - loss: 2.0977 - val_loss: 2.0887\n",

      "Epoch 4/20\n",

      "14/14 [==============================] - 0s 15ms/step - loss: 2.0177 - val_loss: 2.0601\n",

      "Epoch 5/20\n",

      "14/14 [==============================] - 0s 11ms/step - loss: 1.9894 - val_loss: 2.1574\n",

      "Epoch 6/20\n",

      "14/14 [==============================] - 0s 13ms/step - loss: 2.0298 - val_loss: 2.0182\n",

      "Epoch 7/20\n",

      "14/14 [==============================] - 0s 11ms/step - loss: 1.8288 - val_loss: 2.0265\n",

      "Epoch 8/20\n",

      "14/14 [==============================] - 0s 14ms/step - loss: 1.7566 - val_loss: 2.0746\n",

      "Epoch 9/20\n",

      "14/14 [==============================] - 0s 12ms/step - loss: 1.8612 - val_loss: 1.9724\n",

      "Epoch 10/20\n",

      "14/14 [==============================] - 0s 17ms/step - loss: 1.7398 - val_loss: 1.9568\n",

      "Epoch 11/20\n",

      "14/14 [==============================] - 0s 13ms/step - loss: 1.7656 - val_loss: 1.9442\n",

      "Epoch 12/20\n",

      "14/14 [==============================] - 0s 10ms/step - loss: 1.7252 - val_loss: 2.0157\n",

      "Epoch 13/20\n",

      "14/14 [==============================] - 0s 16ms/step - loss: 1.7534 - val_loss: 1.9209\n",

      "Epoch 14/20\n",

      "14/14 [==============================] - 0s 16ms/step - loss: 1.6702 - val_loss: 1.8983\n",

      "Epoch 15/20\n",

      "14/14 [==============================] - 0s 22ms/step - loss: 1.7617 - val_loss: 1.9076\n",

      "Epoch 16/20\n",

      "14/14 [==============================] - 0s 15ms/step - loss: 1.6871 - val_loss: 1.8824\n",

      "Epoch 17/20\n",

      "14/14 [==============================] - 0s 29ms/step - loss: 1.6876 - val_loss: 1.9560\n",

      "Epoch 18/20\n",

      "14/14 [==============================] - 0s 14ms/step - loss: 1.7251 - val_loss: 1.8641\n",

      "Epoch 19/20\n",

      "14/14 [==============================] - 0s 17ms/step - loss: 1.6578 - val_loss: 1.8538\n",

      "Epoch 20/20\n",

      "14/14 [==============================] - 0s 12ms/step - loss: 1.6470 - val_loss: 1.9490\n",

      "1/1 [==============================] - 0s 79ms/step - loss: 2.3078\n",

      "\n",

      "MSE in prediction = 2.307785749435425\n"

     ]

    }

   ],

   "source": [

    "# Controlling overfit: early stopping\n",

    "# Split the train set into proper train & validation\n",

    "X_train0, X_val, y_train0, y_val = train_test_split(X_train, y_train, test_size=0.1)\n",

    "\n",

    "# This is the number of times all data are considered for parameter estimation\n",

    "nEpochs=20\n",

    "\n",

    "# Early stopping means not enough iteration to achieve convergence are allowed\n",

    "early_stopper = keras.callbacks.EarlyStopping(monitor='val_loss', patience=10, min_delta=0.01)\n",

    "model.fit(X_train0, y_train0, epochs=nEpochs, verbose=1, validation_data=(X_val, y_val), callbacks=[early_stopper])\n",

    "\n",

    "# cross-validation\n",

    "mse_prediction = model.evaluate(X_test, y_test, batch_size=128)\n",

    "print('\\nMSE in prediction =',mse_prediction)\n",

    "\n",

    "## In this case neither l1 nor l2 regularization helps"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 12,

   "metadata": {},

   "outputs": [],

   "source": [

    "## We will use keras tuner to optimize hyperparameter search\n",

    "# https://keras.io/keras_tuner/\n",

    "# see also Chollet book 2nd ed. chapter 13\n",

    "'''\n",

    "To put the whole hyperparameter search space together and perform hyperparameter tuning, \n",

    "Keras Tuners uses `HyperModel` instances, i.e., a reusable class object \n",

    "I found defining a 'class' is much better than a function for using keras tuner,\n",

    "and much more elegant!\n",

    "'''\n",

    "from kerastuner import HyperModel\n",

    "\n",

    "class build_model(HyperModel):\n",

    "    def __init__(self, input_dim):\n",

    "        # input and output sizes, if needed, are defined\n",

    "        self.input_dim = input_dim\n",

    "        \n",

    "    def build(self, hp):\n",

    "        model = keras.Sequential()\n",

    "        \n",

    "        # first layer\n",

    "        model.add(\n",

    "            layers.Dense(\n",

    "                input_dim=self.input_dim,\n",

    "                # Define the hyperparameter search for # neurons, units is integer\n",

    "                units=hp.Int(\"units\", min_value=32, max_value=128, step=32),\n",

    "                # choose activation function to use\n",

    "                activation=hp.Choice(\"activation\", [\"relu\", \"tanh\"])\n",

    "            )\n",

    "        )\n",

    "\n",

    "        # second layer, dropout rate is float\n",

    "        model.add(layers.Dense(32, activation=\"softplus\"))\n",

    "        model.add(\n",

    "            layers.Dropout(\n",

    "                rate=hp.Float(\"rate\", min_value=0.0, max_value=0.5, step=0.10)\n",

    "            )\n",

    "        )\n",

    "        \n",

    "        # output layer\n",

    "        model.add(layers.Dense(1))\n",

    "       \n",

    "        # options in compiling\n",

    "        model.compile(loss='mse', optimizer='sgd', metrics=[\"mse\"],)\n",

    "        \n",

    "        return model\n"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 13,

   "metadata": {},

   "outputs": [],

   "source": [

    "hypermodel = build_model(input_dim=nSNP)"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 14,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Search space summary\n",

      "Default search space size: 3\n",

      "units (Int)\n",

      "{'default': None, 'conditions': [], 'min_value': 32, 'max_value': 128, 'step': 32, 'sampling': None}\n",

      "activation (Choice)\n",

      "{'default': 'relu', 'conditions': [], 'values': ['relu', 'tanh'], 'ordered': False}\n",

      "rate (Float)\n",

      "{'default': 0.0, 'conditions': [], 'min_value': 0.0, 'max_value': 0.5, 'step': 0.1, 'sampling': None}\n"

     ]

    }

   ],

   "source": [

    "'''\n",

    "After defining the search space, we need to select a tuner class to run the search. \n",

    "You may choose from RandomSearch, BayesianOptimization and Hyperband.\n",

    "\n",

    "To initialize the tuner, we need to specify several arguments in the initializer.\n",

    "\n",

    "hypermodel: The model-building function, which is build_model in our case.\n",

    "objective: The name of the objective to optimize (whether to minimize or maximize is automatically inferred for built-in metrics). We will introduce how to use custom metrics later in this tutorial.\n",

    "executions_per_trial: The number of models that should be built and fit for each trial. \n",

    "                      Different trials have different hyperparameter values. \n",

    "                      The executions within the same trial have the same hyperparameter values. The purpose of having multiple executions per trial is to reduce results variance and therefore be able to more accurately assess the performance of a model. If you want to get results faster, you could set executions_per_trial=1 (single round of training for each model configuration).\n",

    "overwrite: Control whether to overwrite the previous results in the same directory \n",

    "           or resume the previous search instead. \n",

    "directory: A path to a directory for storing the search results.\n",

    "project_name: The name of the sub-directory in the directory.\n",

    "\n",

    "'''\n",

    "\n",

    "tuner = kt.Hyperband(hypermodel,\n",

    "                     objective=\"val_mse\",\n",

    "                     max_epochs=10,\n",

    "                     overwrite=True)\n",

    "\n",

    "# search summary\n",

    "tuner.search_space_summary()"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 15,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Trial 30 Complete [00h 00m 03s]\n",

      "val_mse: 0.805836021900177\n",

      "\n",

      "Best val_mse So Far: 0.805836021900177\n",

      "Total elapsed time: 00h 00m 58s\n",

      "INFO:tensorflow:Oracle triggered exit\n"

     ]

    }

   ],

   "source": [

    "# start search\n",

    "'''\n",

    "Then, start the search for the best hyperparameter configuration. \n",

    "All the arguments passed to search is passed to model.fit() in each execution. \n",

    "Remember to pass validation_data to evaluate the model.\n",

    "'''\n",

    "\n",

    "tuner.search(X_train0, y_train0, epochs=5, validation_data=(X_val, y_val),\n",

    "            # Use the TensorBoard callback.\n",

    "            # The logs will be write to \"/tmp/tb_logs\".\n",

    "            callbacks=[keras.callbacks.TensorBoard(\"/tmp/tb_logs\")],\n",

    ")"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 16,

   "metadata": {},

   "outputs": [

    {

     "data": {

      "text/html": [

       "\n",

       "      <iframe id=\"tensorboard-frame-ba82ea3297b11799\" width=\"100%\" height=\"800\" frameborder=\"0\">\n",

       "      </iframe>\n",

       "      <script>\n",

       "        (function() {\n",

       "          const frame = document.getElementById(\"tensorboard-frame-ba82ea3297b11799\");\n",

       "          const url = new URL(\"/\", window.location);\n",

       "          const port = 6006;\n",

       "          if (port) {\n",

       "            url.port = port;\n",

       "          }\n",

       "          frame.src = url;\n",

       "        })();\n",

       "      </script>\n",

       "    "

      ],

      "text/plain": [

       "<IPython.core.display.HTML object>"

      ]

     },

     "metadata": {},

     "output_type": "display_data"

    }

   ],

   "source": [

    "# click on HPARAMS, click on TABLE VIEW and in sorting by ascending validation.epoch_mse\n",

    "%load_ext tensorboard\n",

    "\n",

    "%tensorboard --logdir /tmp/tb_logs"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 17,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Model: \"sequential\"\n",

      "_________________________________________________________________\n",

      "Layer (type)                 Output Shape              Param #   \n",

      "=================================================================\n",

      "dense (Dense)                (None, 128)               20736     \n",

      "_________________________________________________________________\n",

      "dense_1 (Dense)              (None, 32)                4128      \n",

      "_________________________________________________________________\n",

      "dropout (Dropout)            (None, 32)                0         \n",

      "_________________________________________________________________\n",

      "dense_2 (Dense)              (None, 1)                 33        \n",

      "=================================================================\n",

      "Total params: 24,897\n",

      "Trainable params: 24,897\n",

      "Non-trainable params: 0\n",

      "_________________________________________________________________\n"

     ]

    }

   ],

   "source": [

    "'''\n",

    "Query the results\n",

    "When search is over, you can retrieve the best model(s). \n",

    "The model is saved at its best performing epoch evaluated on the validation_data.\n",

    "'''\n",

    "# Get the top 2 models.\n",

    "models = tuner.get_best_models(num_models=2)\n",

    "best_model = models[0]\n",

    "\n",

    "# Build the model.\n",

    "best_model.build()\n",

    "best_model.summary()"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 18,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Results summary\n",

      "Results in ./untitled_project\n",

      "Showing 10 best trials\n",

      "Objective(name='val_mse', direction='min')\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 128\n",

      "activation: relu\n",

      "rate: 0.1\n",

      "tuner/epochs: 10\n",

      "tuner/initial_epoch: 0\n",

      "tuner/bracket: 0\n",

      "tuner/round: 0\n",

      "Score: 0.805836021900177\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 32\n",

      "activation: tanh\n",

      "rate: 0.0\n",

      "tuner/epochs: 4\n",

      "tuner/initial_epoch: 0\n",

      "tuner/bracket: 1\n",

      "tuner/round: 0\n",

      "Score: 0.8587983250617981\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 96\n",

      "activation: tanh\n",

      "rate: 0.2\n",

      "tuner/epochs: 4\n",

      "tuner/initial_epoch: 0\n",

      "tuner/bracket: 1\n",

      "tuner/round: 0\n",

      "Score: 0.8632120490074158\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 96\n",

      "activation: relu\n",

      "rate: 0.2\n",

      "tuner/epochs: 4\n",

      "tuner/initial_epoch: 0\n",

      "tuner/bracket: 1\n",

      "tuner/round: 0\n",

      "Score: 0.8777315616607666\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 96\n",

      "activation: relu\n",

      "rate: 0.4\n",

      "tuner/epochs: 10\n",

      "tuner/initial_epoch: 0\n",

      "tuner/bracket: 0\n",

      "tuner/round: 0\n",

      "Score: 0.8842925429344177\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 64\n",

      "activation: tanh\n",

      "rate: 0.1\n",

      "tuner/epochs: 10\n",

      "tuner/initial_epoch: 0\n",

      "tuner/bracket: 0\n",

      "tuner/round: 0\n",

      "Score: 0.8985827565193176\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 64\n",

      "activation: tanh\n",

      "rate: 0.30000000000000004\n",

      "tuner/epochs: 4\n",

      "tuner/initial_epoch: 2\n",

      "tuner/bracket: 2\n",

      "tuner/round: 1\n",

      "tuner/trial_id: 70047aa37f3d81e0f41b100beca637b1\n",

      "Score: 0.9127152562141418\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 64\n",

      "activation: tanh\n",

      "rate: 0.30000000000000004\n",

      "tuner/epochs: 10\n",

      "tuner/initial_epoch: 4\n",

      "tuner/bracket: 2\n",

      "tuner/round: 2\n",

      "tuner/trial_id: 5f77a16951a54fe8483bc612bd85bc47\n",

      "Score: 0.9379919171333313\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 64\n",

      "activation: tanh\n",

      "rate: 0.30000000000000004\n",

      "tuner/epochs: 2\n",

      "tuner/initial_epoch: 0\n",

      "tuner/bracket: 2\n",

      "tuner/round: 0\n",

      "Score: 0.9548115730285645\n",

      "Trial summary\n",

      "Hyperparameters:\n",

      "units: 96\n",

      "activation: tanh\n",

      "rate: 0.2\n",

      "tuner/epochs: 10\n",

      "tuner/initial_epoch: 4\n",

      "tuner/bracket: 1\n",

      "tuner/round: 1\n",

      "tuner/trial_id: a5be4becbb3a1fa352da9966dbdcfce3\n",

      "Score: 0.9745121002197266\n"

     ]

    }

   ],

   "source": [

    "# print summary\n",

    "tuner.results_summary()"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 19,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "N / class\n",

      "Test  [86, 331, 61]\n",

      "Train  [35, 71, 14]\n",

      "All  [121, 402, 75]\n",

      "1/1 [==============================] - 0s 188ms/step - loss: 2.4009\n",

      "\n",

      "MSE in prediction = 2.400940179824829\n",

      "\n",

      "Probabilities matrix\n",

      " [[0.23670393 0.62736833 0.13386099 0.00206673]\n",

      " [0.08689652 0.78197324 0.12901944 0.00211086]\n",

      " [0.03321762 0.58404362 0.35806727 0.02467138]\n",

      " [       nan        nan        nan        nan]]\n"

     ]

    },

    {

     "name": "stderr",

     "output_type": "stream",

     "text": [

      "/home/miguel/anaconda3/envs/tensor/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3474: RuntimeWarning: Mean of empty slice.\n",

      "  return _methods._mean(a, axis=axis, dtype=dtype,\n",

      "/home/miguel/anaconda3/envs/tensor/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in true_divide\n",

      "  ret = ret.dtype.type(ret / rcount)\n"

     ]

    },

    {

     "data": {

      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWsAAAEGCAYAAACjLLT8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAY8klEQVR4nO3dfbxdVX3n8c/3XkAQEFCgQ5OAEYMMD5Un4wPjEzPUFC20U61BrMVSIq8hFaV1xIdipTCvUac4jKa+SIEiUyEIoqYQjVQhVHwgwfKUMJGYtiTRFhEQEQWC3/nj7JucXO49Zx/uOWfvffN989qvnL3P3mv9znmF31lZe+21ZJuIiKi3kaoDiIiI7pKsIyIaIMk6IqIBkqwjIhogyToiogF2qDqAySxb/UCGqRTe+ckVVYdQGzefd0LVIdTG7H13rTqE2th5BzTVMnY5cmHpnPOLf/r0lOvrVW2TdUTEUKneHQ1J1hERABp6Y7knSdYREZCWdUREI6RlHRHRACOjVUfQUZJ1RASkGyQiohHSDRIR0QBpWUdENEBa1hERDZCWdUREA2Q0SEREA6RlHRHRACPps46IqL+0rCMiGiCjQSIiGiA3GCMiGiDdIBERDZBukIiIBqh5y7re0UVEDItUfutalOZJWitpnaRzJnj/k5LuKLbvS3qkW5lpWUdEQN9a1pJGgUXA8cBGYKWkpbbXjJ1j+71t5/8JcGS3ctOyjoiA1miQsltnc4F1ttfbfhJYApzU4fyTgau6hlf6g0RETGcaKb1JWiBpVdu2oK2kGcCGtv2NxbFnVikdAMwGvtEtvHSDRERAT6NBbC8GFveh1vnAtbaf7nZiknVEBPRzNMgmYFbb/szi2ETmA2eWKTTdIBER0M/RICuBOZJmS9qJVkJe+szqdDCwF/DtMuENrGVdBHISW/tqNgFLbd87qDojIp61PrWsbW+WtBBYDowCl9leLek8YJXtscQ9H1hi22XKHUjLWtL7ad0BFXBbsQm4aqIxh23Xbem0/8o1VwwitIiICWlkpPTWje1ltg+yfaDtC4pj57Ylamz/he1J8+F4g2pZnwYcavup9oOSLgRWA/9zoovaO+2XrX6g1K9NREQ/qOaPmw+qz/pXwK9PcHy/4r2IiHpRD1sFBtWyfg/wdUn3sXW84f7Ai4GFA6ozIuJZq3vLeiDJ2vZXJR1E60me9huMK8uMJ4yIGLbtMlkD2P4V8J1BlR8R0U8jJW4cVikPxUREQGV90WUlWUdEsB13g0RENEmSdUREAyRZR0Q0QJJ1REQDaCTJOiKi9tKyjohogCTriIgmqHeuTrKOiIC0rCMiGiHJOiKiATI3SEREE9S7YZ1kHREB9e8GqXe7PyJiSCSV3kqUNU/SWknrJlt3VtLvS1ojabWkK7uVmZZ1RAT9a1lLGgUWAccDG4GVkpbaXtN2zhzgA8Cxth+WtG+3cpOsIyLo6+Pmc4F1ttcDSFoCnASsaTvndGCR7YcBbD/QrdB0g0RE0Fs3iKQFkla1bQvaiprB1rVnodW6nrFtbRwEHCTpVknfkTSvW3xpWUdE0Fs3iO3FwOIpVLcDMAd4HTATuEXS4bYfmeyCtKwjIujrDcZNwKy2/ZnFsXYbgaW2n7L9z8D3aSXvSSVZR0RAa5x12a2zlcAcSbMl7QTMB5aOO+dLtFrVSNqbVrfI+k6F1rYb5PD99qg6hNp49PYVVYdQG9Ibqw4hpql+jQaxvVnSQmA5MApcZnu1pPOAVbaXFu/9pqQ1wNPA+2z/pFO5tU3WERHDNNLHxQdsLwOWjTt2bttrA2cXWylJ1hER1P8JxiTriAig5rk6yToiAtKyjohohJrn6iTriAjo7w3GQUiyjoggyToiohHSDRIR0QC5wRgR0QBJ1hERDVDzXJ1kHREBucEYEdEI6QaJiGiAmufqJOuICEjLOiKiEWqeq5OsIyKg/i3rrst6STpW0q7F67dLulDSAYMPLSJieEZGVHqrJL4S53wGeFzSS4E/BX4AXDHQqCIihkwqv1WhTLLeXCxBcxLwaduLgN0HG1ZExHD1cXXzgSiTrH8m6QPA24EbJI0AOw42rIiI4epny1rSPElrJa2TdM4E758q6ceS7ii2P+5WZpkbjG8F3gacZvvfJO0PfKLEdRERjdGvFrOkUWARcDywEVgpaantNeNOvdr2wrLllknWPwMusv20pIOAg4GrylYQEdEEfezemAuss72+KHcJrW7k8cm6J2W6QW4BniNpBvA14A+Ay6dSaURE3fQyGkTSAkmr2rYFbUXNADa07W8sjo33e5LuknStpFnd4ivTspbtxyWdBvy17Y9LurPEdRERjdFLw9r2YmDxFKr7e+Aq209IehfwWeC4TheUaVlL0iuBU4AberguIqIx+jgaZBPQ3lKeWRzbwvZPbD9R7F4CHN2t0DJJ9yzgA8AXba+W9CLgphLXTUjSOzu8t+WfFp/77CXPtoqIiJ71cTTISmCOpNmSdgLmA0u3rUv7te2eCNzbrdCu3SC2b6HVbz22vx54d9dwJ/dR4G8nqWvLPy02PPSEp1BHRERPRvp0g9H2ZkkLgeXAKHBZ0dA9D1hleynwbkknApuBh4BTu5XbNVlL2gf478ChwM5tAU3avyLprsneAn6tW50REcPWz8fIbS8Dlo07dm7b6w/Q6rEorcwNxs8BVwNvAs4A/hD4cZdrfg14A/DwuOMCvtVLgBERw1DzhWJKJesX2L5U0lm2VwArJK3scs31wG627xj/hqSbew8zImKw6j7rXplk/VTx548kvRH4IfD8ThfYPq3De28rH15ExHDUPFeXStbnS9qD1ox7nwKeB7x3oFFFRAyZqHe2LjMa5Pri5U+B1w82nIiIajS2z1rSp4BJh8/ZnsrwvYiIWqlqUYGyOrWsVw0tioiIivVrnPWgTJqsbX92mIFERFSp5rm61BqMN0ras21/L0nLBxpVRMSQ1X2lmDKjQfax/cjYju2HJe07uJAiIoav7i3rMsn6aUn7274foFjZPPN2RMS0MlrzbF0mWX8I+KakFbQeF381sKDzJRERzdL4Jxhtf1XSUcArikPvsf3gYMOKiBiumo/cK9WypkjO13c9MSKioRrfso6I2B7UPFcnWUdEQINb1pK6zaz3UP/DiYioxmjNO607taxvpzVET8D+tBYSELAncD8we9DBRUQMS71TdYcnGG3Ptv0i4B+A37a9t+0X0Fox5mvDCjAiYhhGpNJbJfGVOOcVxXpiANj+CvCqwYUUETF8fVzdHEnzJK2VtE7SOR3O+z1JlnRMtzLL3GD8oaQPA39X7J9Ca7WYiIhpo183GCWNAouA44GNwEpJS22vGXfe7sBZwHfLlFumZX0ysA/wReC64vXJ5UOPiKi/Pras5wLrbK+3/SSwBDhpgvP+EvgY8Msy8ZV5gvEh4CxJu9r+eZlCIyKappfRIJIWsO20G4ttLy5ezwA2tL23EXj5uOuPAmbZvkHS+8rU2TVZS3oVcAmwG7C/pJcC77L938pUEBHRBL10gxSJeXHXEyeuZwS4EDi1l+vK9Fl/EngDsBTA9p2SXtNrgL3afZc8rzNmv9efUHUItXHlXZuqDqE23vfaA6sOoTZ23qFMj25nUy9hi03ArLb9mcWxMbsDhwE3Fz8Q/wFYKulE25Ou0FV2bpAN4351ni4ZdEREI/TxCcaVwBxJs2kl6fnA28betP1TYO+2em8G/qxTooZyyXpD0RViSTvSunt5b8/hR0TUWL8eYLS9WdJCYDkwClxme7Wk84BVtpc+m3LLJOszgItodZpvovVATPqrI2Ja6efj5sWzKcvGHTt3knNfV6bMMsn6JbZPaT8g6Vjg1jIVREQ0Qc2nBinVp/6pksciIhqrn08wDkKnWfdeSeux8n0knd321vNo9cNEREwbVc35UVanbpCdaI2t3oHWUJMxjwJvHmRQERHD1sehewMxabK2vQJYIely2/86xJgiIoau5g3rUj8ml0jac2xH0l6Slg8upIiI4RsdUemtCmVGg+xt+5GxHdsPS9p3cCFFRAxf3UeDlEnWv5K0v+37ASQdQGsFmYiIaaPJNxjHfAj4pqQVtFa+eTXbzjYVEdF4Nc/VpaZI/Woxnd8rikPvsf3gYMOKiBiuuneDTHqDUdLBxZ9H0Vow94fFtn9xLCJi2lAP/1WhU8v6T4HTgb+a4D0Dxw0kooiICvRhltWB6jTO+vTiz9cPL5yIiGr0cYrUgej0uPl/7XSh7ev6H05ERDXq3mfdqRvkt4s/96U1R8g3iv3XA9+itXhuRMS0UPOGdcdukHcCSPoacIjtHxX7+wGXDyW6iIghmQ7jrGeNJerCv9MaHRIRMW2MNvUGY5uvF3OBXFXsvxX4h8GFFBExfCMVDckrq8xDMQsl/S4wtqL5YttfHGxYERHDVfNekNJTuH4PuMH2e4HlknbvdkFERJOMqPzWjaR5ktZKWifpnAneP0PS3ZLukPRNSYd0ja9EpacD1wIXF4dmAF/qHm5ERHOMSKW3TiSNAouA3wIOAU6eIBlfaftw20cAHwcu7Bpfic9wJnAsrRVisH0freF8ERHTRh/XYJwLrLO93vaTwBLgpPYTbD/atrsrJWYyLZOsnygqLD6QdihTsKSDJf1nSbuNOz6vRJ0REUPVy+IDkhZIWtW2tc9EOgPY0La/sTi2DUlnSvoBrZb1u7vFVyZZr5D0QWAXSccD1wB/3+kCSe8Gvgz8CXCPpPZflf/R4botX8Dll/5NidAiIvpjpIfN9mLbx7Rti3utz/Yi2wcC7wc+3O38MkP33g/8MXA38C5gGXBJl2tOB462/ZikFwLXSnqh7Ytg8vExxQdeDPDIL57OAgcRMTR9nBtkEzCrbX9mcWwyS4DPdCu0Y7IuOspX2z4Y6KWpO2L7MQDb/yLpdbQS9gF0SNYREVXpY2JaCcyRNJtWkp4PvG2buqQ5xf0/gDcC99FFx24Q208DayX1+sTiv0s6oq2cx4A3AXsDh/dYVkTEwPVrNIjtzcBCYDlwL/B526slnSfpxOK0hZJWS7oDOBv4w27xlekG2QtYLek24OdtAZ04+SW8A9g8wQd4h6SLJ74kIqI6/fwnv+1ltLqM24+d2/b6rF7LLJOs/7zXQm1v7PDerb2WFxExaCM1nyO103zWOwNnAC+mdXPx0qJ1HBEx7dR8HqeOLevPAk8B/8jWJ3F6brpHRDRBY1eKoTWH9eEAki4FbhtOSBERw1fvVN05WT819sL25rr/6kRETEXdc1ynZP1SSWPPr4vWE4yPFq9t+3kDjy4iYkhGm5qsbY8OM5CIiCrVO1WXG7oXETHt1bxhnWQdEQHTYFmviIjtQVrWERENoLSsIyLqr7GjQSIitic1z9VJ1hERkGQdEdEI6bOOiGiAms+QmmQdEQF0XQGmaknWERHUvxuk7vNtR0QMxYjKb91ImidpraR1ks6Z4P2zJa2RdJekrxeLiXeO79l9rIiI6UU9/NexHGkUWMTWRVtOlnTIuNP+CTjG9m8A1wIf7xZfknVEBK2he2W3LuYC62yvt/0ksAQ4qf0E2zfZfrzY/Q4ws1uhSdYREbSmSC27dTED2NC2v7E4NpnTgK90K7S2Nxj33CXTaY9Zf+EJVYcQMe318ri5pAXAgrZDi20v7rVOSW8HjgFe2+3c2ibriIih6mEwSJGYJ0vOm4BZbfszi2PbVif9F+BDwGttP9GtznSDRETQvxuMwEpgjqTZknYC5gNLt6lLOhK4GDjR9gNl4kvLOiKC/s0NUiwwvhBYDowCl9leLek8YJXtpcAngN2Aa4qFeu+3fWLH+Gz3J8L+q21gEVE7U061K9f/tHTOedmL9hj6EzRpWUdEQO1XzE2yjoggc4NERDRCvVN1knVEREvNs3WSdUQE9Z91L8k6IoIs6xUR0QhJ1hERDZBukIiIBkjLOiKiAWqeq5OsIyKA2mfrJOuICNJnHRHRCGUWwq1SknVEBKQbJCKiCdINEhHRABm6FxHRADXP1UnWERFA7bN1knVEBPVffCCrm0dE0GpYl926liXNk7RW0jpJ50zw/mskfU/SZklvLhNfknVEBPQtW0saBRYBvwUcApws6ZBxp90PnApcWTa8dINERNDXoXtzgXW21wNIWgKcBKwZO8H2vxTv/apsoQNrWUuaK+llxetDJJ0t6YRB1RcRMRVSL5sWSFrVti1oK2oGsKFtf2NxbEoG0rKW9BFa/wTYQdKNwMuBm4BzJB1p+4JJrlsALAC4+OKLWbBgwUSnRUT0XS/3F20vBhYPLJgJDKob5M3AEcBzgH8DZtp+VNL/Ar4LTJisx30BHlBsERHP0MdukE3ArLb9mcWxKRlUN8hm20/bfhz4ge1HAWz/AijdRxMRMSy9dIN0sRKYI2m2pJ2A+cDSqcY3qGT9pKTnFq+PHjsoaQ+SrCOihvo1dM/2ZmAhsBy4F/i87dWSzpN0IoCkl0naCLwFuFjS6q7x2f3vbZD0HNtPTHB8b2A/23eXKCbdIBFR1pT7MDY+/ETpnDNzr+cM/QmagfRZT5Soi+MPAg8Oos6IiKmp9xOMGWcdEUEWH4iIaISaTw2SZB0RAVl8ICKiGeqdq5OsIyKg9rk6yToiAtJnHRHRCKp5tk6yjogg3SAREY1Q84Z1knVEBGToXkREI6RlHRHRAEnWERENkG6QiIgGSMs6IqIBap6rk6wjIoDaZ+sk64gI0mcdEdEIWXwgIqIJkqwjIuov3SAREQ1Q96F7I1UH0IHqsEl6V9Ux1GXLd5HvosbfxZTtvAMqu/Wjvl7VOVnXxYKqA6iRfBdb5bvYKt/FECRZR0Q0QJJ1REQDJFl3t7jqAGok38VW+S62yncxBLJddQwREdFFWtYREQ2QZB0R0QBJ1pOQNE/SWknrJJ1TdTxVknSZpAck3VN1LFWSNEvSTZLWSFot6ayqY6qKpJ0l3SbpzuK7+GjVMU136bOegKRR4PvA8cBGYCVwsu01lQZWEUmvAR4DrrB9WNXxVEXSfsB+tr8naXfgduB3tse/F5IE7Gr7MUk7At8EzrL9nYpDm7bSsp7YXGCd7fW2nwSWACdVHFNlbN8CPFR1HFWz/SPb3yte/wy4F5hRbVTVcMtjxe6OxZaW3wAlWU9sBrChbX8j2+n/lDExSS8EjgS+W3EolZE0KukO4AHgRtvb7XcxDEnWET2StBvwBeA9th+tOp6q2H7a9hHATGCupO22i2wYkqwntgmY1bY/szgW27mif/YLwOdsX1d1PHVg+xHgJmBexaFMa0nWE1sJzJE0W9JOwHxgacUxRcWKm2qXAvfavrDqeKokaR9Jexavd6F1M/7/VRrUNJdkPQHbm4GFwHJaN5E+b3t1tVFVR9JVwLeBl0jaKOm0qmOqyLHAHwDHSbqj2E6oOqiK7AfcJOkuWo2bG21fX3FM01qG7kVENEBa1hERDZBkHRHRAEnWERENkGQdEdEASdYREQ2QZB2TkjRT0pcl3SfpB5IuKsadI+lUSZ+uOsbxJD3W/awt5/6FpD8bZDwR/ZJkHRMqHgC5DviS7TnAQcBuwAUDrHOHQZUd0XRJ1jGZ44Bf2v5baM0DAbwX+CNJzy3OmSXp5qLl/REASbtKuqGY5/geSW8tjh8taYWk2yUtL6Ybpbj+f0taBXxI0r9KGmkra4OkHSUdKOmrxfX/KOng4pzZkr4t6W5J50/2YSS9Q9JdRVz/d4L3T5e0snj/C2OfUdJbis9xp6RbimOHFnM531GUOacv33hEB2nJxGQOpTVf8xa2H5V0P/Di4tBc4DDgcWClpBuAA4Af2n4jgKQ9ivk0PgWcZPvHRQK/APijopydbB9TnH8U8Fpac028CVhu+ylJi4EzbN8n6eXAX9P6QbkI+IztKySdOdEHkXQo8GHgVbYflPT8CU67zvbfFOefD5xWxHwu8Abbm8YerwbOAC6y/bmiW2i0xPcZMSVpWcdU3Gj7J7Z/QavL5D8BdwPHS/qYpFfb/inwElpJ/cZiSs0P05oca8zV416/tXg9H7i6mOXuVcA1xfUX03rcGVqPgF9VvH5Gi7lwHHCN7QcBbE80N/dhRYv9buAUWj9WALcCl0s6na1J+dvAByW9Hzig+PwRA5VkHZNZAxzdfkDS84D9gXXFofFzFdj294GjaCXt8yWdCwhYbfuIYjvc9m+2XffzttdLgXlF6/do4Bu0/p4+0nb9Ebb/Y3u9U/uoAFwOLLR9OPBRYOfiA51B68dlFnC7pBfYvhI4EfgFsEzScX2oP6KjJOuYzNeB50p6B2xZ6uyvgMttP16cc7yk5xezrv0OcKukXwcet/13wCdoJe61wD6SXlmUtWPRNfEMxeojK2l1b1xfzJn8KPDPkt5SXC9JLy0uuZVWCxxaLeKJfAN4i6QXFNdP1A2yO/CjostmSzmSDrT9XdvnAj+m1U//ImC97f8DfBn4jUnqjeibJOuYkFszfP0urSR3H601KX8JfLDttNtoze18F/AF26uAw4Hbiu6KjwDnF0ujvRn4mKQ7gTtodWtM5mrg7WzbPXIKcFpx/Wq2LrN2FnBm0X0x4Wo+xYyJFwAriusnmt70z2mt+nIr2071+Yni5uU9wLeAO4HfB+4pPuNhwBUdPktEX2TWvYiIBkjLOiKiAZKsIyIaIMk6IqIBkqwjIhogyToiogGSrCMiGiDJOiKiAf4/cU5yBLhcTMQAAAAASUVORK5CYII=\n",

      "text/plain": [

       "<Figure size 432x288 with 2 Axes>"

      ]

     },

     "metadata": {

      "needs_background": "light"

     },

     "output_type": "display_data"

    }

   ],

   "source": [

    "# MLP example with multiclass target\n",

    "\n",

    "# Let us round class vector and make a few classes, all positive numbers\n",

    "# just a trick to convert to classes\n",

    "yi_train=[int(round(x-np.min(y_train))/2) for x in y_train]\n",

    "yi_test=[int(round(x-np.min(y_test))/2) for x in y_test]\n",

    "\n",

    "# N obs / clas; this may result in some very rare classes so consider merging\n",

    "print('N / class')\n",

    "print('Test ',[yi_train.count(i) for i in range(max(yi_train+yi_test))])\n",

    "print('Train ',[yi_test.count(i) for i in range(max(yi_train+yi_test))])\n",

    "print('All ',[(yi_test+yi_train).count(i) for i in range(max(yi_train+yi_test))])\n",

    "# WARNING: make sure all clases in test are in train!!\n",

    "\n",

    "# convert to classes, i need to make all classes equivalent\n",

    "n_train=len(yi_train)\n",

    "itemp = keras.utils.to_categorical(yi_train+yi_test)\n",

    "i_train = itemp[:n_train,:]\n",

    "i_test = itemp[n_train:,:]\n",

    "\n",

    "# no. of SNPs in data\n",

    "nSNP=X_train.shape[1]\n",

    "nClasses=i_train.shape[1]\n",

    "\n",

    "# Instantiate\n",

    "model = Sequential()\n",

    "\n",

    "# Add first layer\n",

    "model.add(Dense(64, input_dim=nSNP))\n",

    "model.add(Activation('relu'))\n",

    "# Add second layer\n",

    "model.add(Dense(32))\n",

    "model.add(Activation('softplus'))\n",

    "# Last, output layer\n",

    "model.add(Dense(nClasses, activation='softmax'))\n",

    "\n",

    "# Model Compiling \n",

    "model.compile(loss='categorical_crossentropy', optimizer='adam')\n",

    "\n",

    "# training\n",

    "model.fit(X_train, i_train, epochs=100, verbose=0)\n",

    "\n",

    "# cross-validation: get predicted target values\n",

    "i_hat = model.predict(X_test, batch_size=128)\n",

    "\n",

    "mse_prediction = model.evaluate(X_test, i_test, batch_size=128)\n",

    "print('\\nMSE in prediction =',mse_prediction)\n",

    "\n",

    "# do a heatplot, obs vs expected class distribution\n",

    "# collect all results by class\n",

    "# compute average prediction by class\n",

    "heat = np.zeros([nClasses,nClasses])\n",

    "for i in range(nClasses):\n",

    "    iclass = np.nonzero(i_test[:,i]>0) # samples of i-th class\n",

    "    for j in range(nClasses):\n",

    "        heat[i,j] = np.mean(i_hat[iclass,j])\n",

    "\n",

    "# plot observed vs. predicted targets\n",

    "print('\\nProbabilities matrix\\n',heat)\n",

    "plot = sns.heatmap(heat, cmap=\"Blues\")\n",

    "plot.set(xlabel='Observed class', ylabel='Predicted class')\n",

    "plt.show()"

   ]

  },

  {

   "cell_type": "code",

   "execution_count": 21,

   "metadata": {},

   "outputs": [

    {

     "name": "stdout",

     "output_type": "stream",

     "text": [

      "Model: \"sequential_2\"\n",

      "_________________________________________________________________\n",

      "Layer (type)                 Output Shape              Param #   \n",

      "=================================================================\n",

      "conv1d (Conv1D)              (None, 53, 32)            128       \n",

      "_________________________________________________________________\n",

      "max_pooling1d (MaxPooling1D) (None, 26, 32)            0         \n",

      "_________________________________________________________________\n",

      "flatten (Flatten)            (None, 832)               0         \n",

      "_________________________________________________________________\n",

      "dense_6 (Dense)              (None, 64)                53312     \n",

      "_________________________________________________________________\n",

      "activation_2 (Activation)    (None, 64)                0         \n",

      "_________________________________________________________________\n",

      "dense_7 (Dense)              (None, 32)                2080      \n",

      "_________________________________________________________________\n",

      "activation_3 (Activation)    (None, 32)                0         \n",

      "_________________________________________________________________\n",

      "dense_8 (Dense)              (None, 1)                 33        \n",

      "=================================================================\n",

      "Total params: 55,553\n",

      "Trainable params: 55,553\n",

      "Non-trainable params: 0\n",

      "_________________________________________________________________\n"

     ]

    }

   ],

   "source": [

    "#--> CNN example\n",

    "nSNP = X_train.shape[1] \n",

    "nStride = 3  # stride between convolutions\n",

    "nFilter = 32 # no. of convolutions\n",

    "\n",

    "# Instantiate\n",

    "model_cnn = Sequential()\n",

    "\n",

    "#WARNING!!! I need this to match dimensions \n",

    "#https://stackoverflow.com/questions/43396572/dimension-of-shape-in-conv1d\n",

    "X2_train = np.expand_dims(X_train, axis=2) \n",

    "X2_test = np.expand_dims(X_test, axis=2) \n",

    "\n",

    "# add convolutional layer\n",

    "model_cnn.add(layers.Conv1D(nFilter, kernel_size=3, strides=nStride, input_shape=(nSNP,1)))\n",

    "# add pooling layer: takes maximum of two consecutive values\n",

    "model_cnn.add(layers.MaxPooling1D(pool_size=2))\n",

    "# Solutions above are linearized to accommodate a standard layer\n",

    "model_cnn.add(layers.Flatten())\n",

    "model_cnn.add(layers.Dense(64))\n",

    "model_cnn.add(layers.Activation('relu'))\n",

    "model_cnn.add(layers.Dense(32))\n",

    "model_cnn.add(layers.Activation('softplus'))\n",

    "model_cnn.add(layers.Dense(1))\n",

    "\n",

    "# Model Compiling (https://keras.io/models/sequential/) \n",