gp_benchmarks_meta_features.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports and Configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-18T20:48:06.336532Z",
     "start_time": "2018-04-18T20:48:04.980871Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import requests\n",
    "import zipfile\n",
    "import io\n",
    "import re\n",
    "import timeit\n",
    "\n",
    "from math import sqrt, sin, cos, log, pi, e\n",
    "from gplearn.genetic import SymbolicRegressor\n",
    "from gplearn.functions import make_function\n",
    "from gplearn.fitness import make_fitness\n",
    "\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from multiprocessing import Pool\n",
    "\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "\n",
    "import os\n",
    "\n",
    "# Change the current folder to the folder where we can find the results\n",
    "current_folder = !pwd\n",
    "os.chdir(current_folder[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-24T19:21:36.997920Z",
     "start_time": "2018-03-24T19:21:36.985871Z"
    }
   },
   "source": [
    "# Retrieving the data\n",
    "We selected the datasets used in experiments in GP in the last five years of GECCO.\n",
    "\n",
    "Here are the auxiliary functions used to recover the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-31T19:54:45.059288Z",
     "start_time": "2018-03-31T19:54:44.488808Z"
    }
   },
   "outputs": [],
   "source": [
    "def synthetic_gen(name, rnd, function, training_gen, test_gen=None):\n",
    "    \"\"\"Generates training and test sets for synthetic datasets.\n",
    "    \"\"\"\n",
    "    print(\"Generating the \\\"\" + name + \"\\\" dataset\")\n",
    "    training_set = []\n",
    "    test_set = []\n",
    "    for i in range(training_gen.n):\n",
    "        inst = training_gen.generate(rnd)\n",
    "        training_set.append(inst + [function(*inst)])\n",
    "    if test_gen:\n",
    "        for i in range(test_gen.n):\n",
    "            inst = test_gen.generate(rnd)\n",
    "            test_set.append(inst + [function(*inst)])\n",
    "    else:\n",
    "        test_set = training_set\n",
    "    return {\"training\": pd.DataFrame(training_set),\n",
    "            \"test\": pd.DataFrame(test_set)}\n",
    "\n",
    "class U:\n",
    "    def __init__(self, ini, end, n):\n",
    "        self.ini = ini\n",
    "        self.end = end\n",
    "        self.n = n\n",
    "    \n",
    "    def generate(self, rnd):\n",
    "        return [rnd.uniform(ini, end) for ini, end in zip(self.ini, self.end)]\n",
    "\n",
    "class E:\n",
    "    def __init__(self, ini, end, step):\n",
    "        self.ini = ini\n",
    "        self.end = end\n",
    "        self.step = step\n",
    "        \n",
    "        mesh = np.meshgrid(*[np.arange(ini, end+step, step) \n",
    "                           for ini, end, step in zip(self.ini, self.end, self.step)])\n",
    "        self.points = [dim.reshape(1,-1)[0] for dim in mesh]\n",
    "        self.index = 0\n",
    "        self.n = len(self.points[0])\n",
    "    \n",
    "    def generate(self, rnd):\n",
    "        inst = [self.points[i][self.index] for i in range(len(self.points))]\n",
    "        self.index += 1\n",
    "        return inst\n",
    "\n",
    "def get_data(name, url, rnd=None, pd_sep=',', pd_header=None, pd_skiprows=None, dataset=None):\n",
    "    print(\"Reading the \\\"\" + name + \"\\\" dataset\")\n",
    "    if dataset == \"BOH\":\n",
    "        from sklearn.datasets import load_boston\n",
    "        boston = load_boston()\n",
    "        df = pd.DataFrame(boston['data'])\n",
    "        df = pd.concat([df, pd.Series(boston['target'])], axis=1)\n",
    "    elif dataset == \"CCP\":\n",
    "        # Get the file object from an url\n",
    "        r = requests.get(url)\n",
    "        # Create a ZipFile object from it\n",
    "        z = zipfile.ZipFile(io.BytesIO(r.content))\n",
    "        # Read from a xlsx file inside the zip file\n",
    "        df = pd.read_excel(z.open('CCPP/Folds5x2_pp.xlsx'))\n",
    "    elif dataset == \"CST\":\n",
    "        df = pd.read_excel(url)\n",
    "    elif dataset == \"ENC\":\n",
    "        df = pd.read_excel(url)\n",
    "        # Drop Y1\n",
    "        df.drop(\"Y1\", axis=1, inplace=True)\n",
    "    elif dataset == \"ENH\":\n",
    "        df = pd.read_excel(url)\n",
    "        # Drop Y2\n",
    "        df.drop(\"Y2\", axis=1, inplace=True)\n",
    "    elif dataset == \"YAC\":\n",
    "        # Get the data as text\n",
    "        text = requests.get(url).text\n",
    "        # Split in rows (remove the last line)\n",
    "        e_re = re.compile(\"\\s*\\n\\s*\")\n",
    "        rows = e_re.split(text)[:-1]\n",
    "        e_re = re.compile(\" +\")\n",
    "        # Split cells per row\n",
    "        df = pd.DataFrame([e_re.split(row) for row in rows])\n",
    "    else:\n",
    "        df = pd.read_csv(url, header=pd_header, sep=pd_sep, skiprows=pd_skiprows)\n",
    "        if dataset == \"ABA\":\n",
    "            # Get dummy variables for the first column\n",
    "            df_dummies = pd.get_dummies(df.iloc[:,0])\n",
    "            # Drop the first column\n",
    "            df.drop(df.columns[0], axis=1, inplace=True)\n",
    "            # Concatenate the dummy variables with the data\n",
    "            df = pd.concat([df_dummies, df], axis=1)\n",
    "            df = df.sample(500, random_state=rnd, axis=0)\n",
    "        elif dataset == \"CPU\":\n",
    "            # Drop the first two columns\n",
    "            df.drop(df.columns[[0,1]], axis=1, inplace=True)\n",
    "        elif dataset == \"FFR\":\n",
    "            df.drop([\"month\", \"day\"], axis=1, inplace=True)\n",
    "        elif dataset == \"OZO\":\n",
    "            # Imputation (replance NaN's by the mean of the column)\n",
    "            df.fillna(df.mean(), inplace=True)\n",
    "    return df.apply(np.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-24T19:35:03.633506Z",
     "start_time": "2018-03-24T19:35:03.624408Z"
    }
   },
   "source": [
    "\n",
    "The synthetic data is generated on the run and the real-word datasets are loaded from online resources, with exception of the Ozone dataset, which is loaded from file, and the Boston dataset, which is loaded using [Scikit](http://scikit-learn.org/stable/datasets/index.html). \n",
    "\n",
    "The Ozone dataset was loaded inside an R evnironment from the [mlbench package](https://cran.r-project.org/web/packages/mlbench/index.html) and then saved as file on the folder `./data`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-31T19:55:55.461892Z",
     "start_time": "2018-03-31T19:54:48.446814Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating the \"Meier-3\" dataset\n",
      "Generating the \"Meier-4\" dataset\n",
      "Generating the \"Nonic\" dataset\n",
      "Generating the \"Sine\" dataset\n",
      "Generating the \"Burks\" dataset\n",
      "Generating the \"R1\" dataset\n",
      "Generating the \"R2\" dataset\n",
      "Generating the \"R3\" dataset\n",
      "Generating the \"Poly-10\" dataset\n",
      "Generating the \"Koza-2\" dataset\n",
      "Generating the \"Koza-3\" dataset\n",
      "Generating the \"Korns-1\" dataset\n",
      "Generating the \"Korns-4\" dataset\n",
      "Generating the \"Korns-7\" dataset\n",
      "Generating the \"Korns-11\" dataset\n",
      "Generating the \"Korns-12\" dataset\n",
      "Generating the \"Vladislavleva-1\" dataset\n",
      "Generating the \"Vladislavleva-2\" dataset\n",
      "Generating the \"Vladislavleva-3\" dataset\n",
      "Generating the \"Vladislavleva-4\" dataset\n",
      "Generating the \"Vladislavleva-5\" dataset\n",
      "Generating the \"Vladislavleva-6\" dataset\n",
      "Generating the \"Vladislavleva-7\" dataset\n",
      "Generating the \"Vladislavleva-8\" dataset\n",
      "Generating the \"Pagie-1\" dataset\n",
      "Generating the \"Keijzer-1\" dataset\n",
      "Generating the \"Keijzer-2\" dataset\n",
      "Generating the \"Keijzer-3\" dataset\n",
      "Generating the \"Keijzer-4\" dataset\n",
      "Generating the \"Keijzer-5\" dataset\n",
      "Generating the \"Keijzer-6\" dataset\n",
      "Generating the \"Keijzer-7\" dataset\n",
      "Generating the \"Keijzer-8\" dataset\n",
      "Generating the \"Keijzer-9\" dataset\n",
      "Generating the \"Keijzer-10\" dataset\n",
      "Generating the \"Keijzer-11\" dataset\n",
      "Generating the \"Keijzer-12\" dataset\n",
      "Generating the \"Keijzer-13\" dataset\n",
      "Generating the \"Keijzer-14\" dataset\n",
      "Generating the \"Keijzer-15\" dataset\n",
      "Generating the \"Nguyen-1\" dataset\n",
      "Generating the \"Nguyen-2\" dataset\n",
      "Generating the \"Nguyen-3\" dataset\n",
      "Generating the \"Nguyen-4\" dataset\n",
      "Generating the \"Nguyen-5\" dataset\n",
      "Generating the \"Nguyen-6\" dataset\n",
      "Generating the \"Nguyen-7\" dataset\n",
      "Generating the \"Nguyen-8\" dataset\n",
      "Generating the \"Nguyen-9\" dataset\n",
      "Generating the \"Nguyen-10\" dataset\n",
      "Reading the \"Abalone\" dataset\n",
      "Reading the \"Airfoil\" dataset\n",
      "Reading the \"Boston\" dataset\n",
      "Reading the \"Combined-cycle\" dataset\n",
      "Reading the \"Computer-hardware\" dataset\n",
      "Reading the \"Concrete-strength\" dataset\n",
      "Reading the \"Energy-cooling\" dataset\n",
      "Reading the \"Energy-heating\" dataset\n",
      "Reading the \"Forest-fire\" dataset\n",
      "Reading the \"Ozone\" dataset\n",
      "Reading the \"Wine-quality-red\" dataset\n",
      "Reading the \"Wine-quality-white\" dataset\n",
      "Reading the \"Yacht\" dataset\n"
     ]
    }
   ],
   "source": [
    "seed = 1234\n",
    "rnd = np.random.RandomState(seed)\n",
    "# Loading synthetic datasets\n",
    "data_synt = {\"meier-3\": synthetic_gen(\"Meier-3\", rnd, \n",
    "                                      lambda x_1,x_2: (x_1**2*x_2**2)/(x_1+x_2), \n",
    "                                      U([-1, -1], [1, 1], 50), U([-1, -1], [1, 1], 50)),\n",
    "            \"meier-4\": synthetic_gen(\"Meier-4\", rnd, \n",
    "                                     lambda x_1,x_2: x_1**5/x_2**3, \n",
    "                                     U([-1, -1], [1, 1], 50), U([-1, -1], [1, 1], 50)),\n",
    "            \"nonic\": synthetic_gen(\"Nonic\", rnd,\n",
    "                                   lambda x_1: sum([x_1**i for i in range(1,10)]), \n",
    "                                   E([-1], [1], [2/19]), U([-1], [1], 20)),\n",
    "            \"sine\": synthetic_gen(\"Sine\", rnd,\n",
    "                                  lambda x_1: sin(x_1), \n",
    "                                  E([0], [6.2], [0.1])),\n",
    "            \"burks\": synthetic_gen(\"Burks\", rnd,\n",
    "                                   lambda x_1: 4*x_1**4 + 3*x_1**3 + 2*x_1**2 + x_1, \n",
    "                                   U([-1], [1], 20)),\n",
    "            \"r1\": synthetic_gen(\"R1\", rnd,\n",
    "                                lambda x_1: (x_1+1)**3/(x_1**2-x_1+1), \n",
    "                                E([-1], [1], [2/19]), U([-1], [1], 20)),\n",
    "            \"r2\": synthetic_gen(\"R2\", rnd,\n",
    "                                lambda x_1: (x_1**5-3*x_1**3+1)/(x_1**2+1), \n",
    "                                E([-1], [1], [2/19]), U([-1], [1], 20)),\n",
    "            \"r3\": synthetic_gen(\"R3\", rnd,\n",
    "                                lambda x_1: (x_1**6+x_1**5)/(x_1**4+x_1**3+x_1**2+x_1+1), \n",
    "                                E([-1], [1], [2/19]), U([-1], [1], 20)),\n",
    "            \"poly-10\": synthetic_gen(\"Poly-10\", rnd,\n",
    "                                     lambda x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_10: \n",
    "                                         x_1*x_2+x_3*x_4+x_5*x_6+x_1*x_7*x_9+x_3*x_6*x_10,\n",
    "                                     U([0]*10, [1]*10, 330), U([0]*10, [1]*10, 170)),\n",
    "            \"koza-2\": synthetic_gen(\"Koza-2\", rnd,\n",
    "                                    lambda x_1: x_1**5-2*x_1**3+x_1, \n",
    "                                    U([-1], [1], 20), U([-1], [1], 20)),\n",
    "            \"koza-3\": synthetic_gen(\"Koza-3\", rnd,\n",
    "                                    lambda x_1: x_1**6-2*x_1**4+x_1**2,\n",
    "                                    U([-1], [1], 20), U([-1], [1], 20)),\n",
    "            \"korns-1\": synthetic_gen(\"Korns-1\", rnd,\n",
    "                                     lambda x_1, x_2, x_3, x_4, x_5: 1.57+24.3*x_4,\n",
    "                                     U([-50]*5, [50]*5, 10000), U([-50]*5, [50]*5, 10000)),\n",
    "            \"korns-4\": synthetic_gen(\"Korns-4\", rnd,\n",
    "                                     lambda x_1, x_2, x_3, x_4, x_5: \n",
    "                                         -2.3+0.13*sin(x_3),\n",
    "                                     U([-50]*5, [50]*5, 10000), U([-50]*5, [50]*5, 10000)),\n",
    "            \"korns-7\": synthetic_gen(\"Korns-7\", rnd,\n",
    "                                     lambda x_1, x_2, x_3, x_4, x_5: \n",
    "                                         213.80940889*(1-e**(-0.54723748542*x_1)),\n",
    "                                     U([-50]*5, [50]*5, 10000), U([-50]*5, [50]*5, 10000)),\n",
    "            \"korns-11\": synthetic_gen(\"Korns-11\", rnd,\n",
    "                                     lambda x_1, x_2, x_3, x_4, x_5: \n",
    "                                         6.87+11*cos(7.23*x_1**3),\n",
    "                                     U([-50]*5, [50]*5, 10000), U([-50]*5, [50]*5, 10000)),\n",
    "            \"korns-12\": synthetic_gen(\"Korns-12\", rnd,\n",
    "                                     lambda x_1, x_2, x_3, x_4, x_5: \n",
    "                                         2-2.1*cos(9.8*x_1)*sin(1.3*x_5),\n",
    "                                     U([-50]*5, [50]*5, 10000), U([-50]*5, [50]*5, 10000)),\n",
    "            \"vladislavleva-1\": synthetic_gen(\"Vladislavleva-1\", rnd,\n",
    "                                             lambda x_1, x_2: \n",
    "                                                 e**(-(x_1-1)**2)/(1.2+(x_2-2.5)**2),\n",
    "                                             U([0.3]*2, [4]*2, 100), E([-0.2]*2, [4.2]*2, [0.1]*2)),\n",
    "            \"vladislavleva-2\": synthetic_gen(\"Vladislavleva-2\", rnd,\n",
    "                                             lambda x_1: \n",
    "                                                 e**(-x_1)*x_1**3*(cos(x_1)*sin(x_1))*(cos(x_1)*sin(x_1)**2-1),\n",
    "                                             E([0.05], [10], [0.1]), E([-0.5], [10.5], [0.05])),\n",
    "            \"vladislavleva-3\": synthetic_gen(\"Vladislavleva-3\", rnd,\n",
    "                                             lambda x_1, x_2: \n",
    "                                                 e**(-x_1)*x_1**3*(cos(x_1)*sin(x_1))*(cos(x_1)*sin(x_1)**2-1)*(x_2-5),\n",
    "                                             E([0.05]*2, [10, 10.05], [0.1, 2]), E([-0.5]*2, [10.5]*2, [0.05, 0.5])),\n",
    "            \"vladislavleva-4\": synthetic_gen(\"Vladislavleva-4\", rnd,\n",
    "                                             lambda x_1, x_2, x_3, x_4, x_5: \n",
    "                                                 10/(5+(x_1-3)**2+(x_2-3)**2+(x_3-3)**2+(x_4-3)**2+(x_5-3)**2),\n",
    "                                             U([0.05]*5, [6.05]*5, 1024), U([-0.25]*5, [6.35]*5, 5000)),\n",
    "            \"vladislavleva-5\": synthetic_gen(\"Vladislavleva-5\", rnd,\n",
    "                                             lambda x_1, x_2, x_3: \n",
    "                                                 30*(x_1-1)*(x_3-1)/((x_1-10)*x_2**2),\n",
    "                                             U([0.05, 1, 0.05], [2]*3, 300), \n",
    "                                             E([-0.05, 0.95, -0.05], [2.1, 2.05, 2.1], [0.15, 0.1, 0.15])),\n",
    "            \"vladislavleva-6\": synthetic_gen(\"Vladislavleva-6\", rnd,\n",
    "                                             lambda x_1, x_2: 6*sin(x_1)*cos(x_2),\n",
    "                                             U([0.1]*2, [5.9]*5, 30), \n",
    "                                             E([-0.05]*2, [6.05]*2, [0.02]*2)),\n",
    "            \"vladislavleva-7\": synthetic_gen(\"Vladislavleva-7\", rnd,\n",
    "                                             lambda x_1, x_2: \n",
    "                                                 (x_1-3)*(x_2-3)+2*sin((x_1-4)*(x_2-4)),\n",
    "                                             U([0.05]*2, [6.05]*2, 300), U([-0.25]*2, [6.35]*2, 1000)),\n",
    "            \"vladislavleva-8\": synthetic_gen(\"Vladislavleva-8\", rnd,\n",
    "                                             lambda x_1, x_2: \n",
    "                                                 ((x_1-3)**4+(x_2-3)**3-(x_2-3))/((x_2-2)**4+10),\n",
    "                                             U([0.05]*2, [6.05]*2, 50), E([-0.25]*2, [6.35]*2, [0.2]*2)),\n",
    "            \"pagie-1\": synthetic_gen(\"Pagie-1\", rnd,\n",
    "                                     lambda x_1, x_2: 1/(1+x_1**(-4))+1/(1+x_2**(-4)),\n",
    "                                     E([-5]*2, [5]*2, [0.4]*2)),\n",
    "            \"keijzer-1\": synthetic_gen(\"Keijzer-1\", rnd,\n",
    "                                       lambda x_1: \n",
    "                                           0.3*x_1*sin(2*pi*x_1),\n",
    "                                       E([-1], [1], [0.1]),\n",
    "                                       E([-1], [1], [0.001])),\n",
    "            \"keijzer-2\": synthetic_gen(\"Keijzer-2\", rnd,\n",
    "                                       lambda x_1: \n",
    "                                           0.3*x_1*sin(2*pi*x_1),\n",
    "                                       E([-2], [2], [0.1]),\n",
    "                                       E([-2], [2], [0.001])),\n",
    "            \"keijzer-3\": synthetic_gen(\"Keijzer-3\", rnd,\n",
    "                                       lambda x_1: \n",
    "                                           0.3*x_1*sin(2*pi*x_1),\n",
    "                                       E([-3], [3], [0.1]),\n",
    "                                       E([-3], [3], [0.001])),\n",
    "            \"keijzer-4\": synthetic_gen(\"Keijzer-4\", rnd,\n",
    "                                       lambda x_1: \n",
    "                                           x_1**3*e**(-x_1)*cos(x_1)*sin(x_1)*(sin(x_1)**2*cos(x_1)-1),\n",
    "                                       E([0], [10], [0.05]),\n",
    "                                       E([0.05], [10.05], [0.05])),\n",
    "            \"keijzer-5\": synthetic_gen(\"Keijzer-5\", rnd,\n",
    "                                       lambda x_1, x_2, x_3: 30*x_1*x_3/((x_1-10)*x_2**2),\n",
    "                                       U([-1, 1, -1], [1,2,1], 1000),\n",
    "                                       U([-1, 1, -1], [1,2,1], 10000)),\n",
    "            \"keijzer-6\": synthetic_gen(\"Keijzer-6\", rnd,\n",
    "                                       lambda x_1: sum([1/i for i in range(1, x_1+1)]),\n",
    "                                       E([1], [50], [1]),\n",
    "                                       E([1], [120], [1])),\n",
    "            \"keijzer-7\": synthetic_gen(\"Keijzer-7\", rnd,\n",
    "                                       lambda x_1: log(x_1),\n",
    "                                       E([1], [100], [1]),\n",
    "                                       E([1], [100], [0.1])),\n",
    "            \"keijzer-8\": synthetic_gen(\"Keijzer-8\", rnd,\n",
    "                                       lambda x_1: sqrt(x_1),\n",
    "                                       E([0], [100], [1]),\n",
    "                                       E([0], [100], [0.1])),\n",
    "            \"keijzer-9\": synthetic_gen(\"Keijzer-9\", rnd,\n",
    "                                       lambda x_1: log(x_1+sqrt(x_1**2+1)),\n",
    "                                       E([0], [100], [1]),\n",
    "                                       E([0], [100], [0.1])),\n",
    "            \"keijzer-10\": synthetic_gen(\"Keijzer-10\", rnd,\n",
    "                                       lambda x_1, x_2: x_1**x_2,\n",
    "                                       U([0]*2, [1]*2, 100),\n",
    "                                       E([0]*2, [1]*2, [0.01]*2)),\n",
    "            \"keijzer-11\": synthetic_gen(\"Keijzer-11\", rnd,\n",
    "                                       lambda x_1, x_2: x_1*x_2+sin((x_1-1)*(x_2-1)),\n",
    "                                       U([-3]*2, [3]*2, 20),\n",
    "                                       E([-3]*2, [3]*2, [0.01]*2)),\n",
    "            \"keijzer-12\": synthetic_gen(\"Keijzer-12\", rnd,\n",
    "                                       lambda x_1, x_2: x_1**4-x_1**3+(x_2**2/2)-x_2,\n",
    "                                       U([-3]*2, [3]*2, 20),\n",
    "                                       E([-3]*2, [3]*2, [0.01]*2)),\n",
    "            \"keijzer-13\": synthetic_gen(\"Keijzer-13\", rnd,\n",
    "                                       lambda x_1, x_2: 6*sin(x_1)*cos(x_2),\n",
    "                                       U([-3]*2, [3]*2, 20),\n",
    "                                       E([-3]*2, [3]*2, [0.01]*2)),\n",
    "            \"keijzer-14\": synthetic_gen(\"Keijzer-14\", rnd,\n",
    "                                       lambda x_1, x_2: 8/(2+x_1**2+x_2**2),\n",
    "                                       U([-3]*2, [3]*2, 20),\n",
    "                                       E([-3]*2, [3]*2, [0.01]*2)),\n",
    "            \"keijzer-15\": synthetic_gen(\"Keijzer-15\", rnd,\n",
    "                                       lambda x_1, x_2: (x_1**3/5)+(x_2**3/2)-x_2-x_1,\n",
    "                                       U([-3]*2, [3]*2, 20),\n",
    "                                       E([-3]*2, [3]*2, [0.01]*2)),\n",
    "            \"nguyen-1\": synthetic_gen(\"Nguyen-1\", rnd,\n",
    "                                       lambda x_1: x_1**3+x_1**2+x_1,\n",
    "                                       U([-1], [1], 20),\n",
    "                                       U([-1], [1], 20)),\n",
    "            \"nguyen-2\": synthetic_gen(\"Nguyen-2\", rnd,\n",
    "                                       lambda x_1: x_1**4+x_1**3+x_1**2+x_1,\n",
    "                                       U([-1], [1], 20),\n",
    "                                       U([-1], [1], 20)),\n",
    "            \"nguyen-3\": synthetic_gen(\"Nguyen-3\", rnd,\n",
    "                                       lambda x_1: x_1**5+x_1**4+x_1**3+x_1**2+x_1,\n",
    "                                       U([-1], [1], 20),\n",
    "                                       U([-1], [1], 20)),\n",
    "            \"nguyen-4\": synthetic_gen(\"Nguyen-4\", rnd,\n",
    "                                       lambda x_1: x_1**6+x_1**5+x_1**4+x_1**3+x_1**2+x_1,\n",
    "                                       U([-1], [1], 20),\n",
    "                                       U([-1], [1], 20)),\n",
    "            \"nguyen-5\": synthetic_gen(\"Nguyen-5\", rnd,\n",
    "                                       lambda x_1: sin(x_1**2)*cos(x_1)-1,\n",
    "                                       U([-1], [1], 20),\n",
    "                                       U([-1], [1], 20)),\n",
    "            \"nguyen-6\": synthetic_gen(\"Nguyen-6\", rnd,\n",
    "                                       lambda x_1: sin(x_1)+sin(x_1+x_1**2),\n",
    "                                       U([-1], [1], 20),\n",
    "                                       U([-1], [1], 20)),\n",
    "            \"nguyen-7\": synthetic_gen(\"Nguyen-7\", rnd,\n",
    "                                       lambda x_1: log(x_1+1)+log(x_1**2+1),\n",
    "                                       U([0], [2], 20),\n",
    "                                       U([0], [2], 20)),\n",
    "            \"nguyen-8\": synthetic_gen(\"Nguyen-8\", rnd,\n",
    "                                       lambda x_1: sqrt(x_1),\n",
    "                                       U([0], [4], 20),\n",
    "                                       U([0], [4], 20)),\n",
    "            \"nguyen-9\": synthetic_gen(\"Nguyen-9\", rnd,\n",
    "                                       lambda x_1, x_2: sin(x_1)+sin(x_2**2),\n",
    "                                       U([-1]*2, [1]*2, 100),\n",
    "                                       U([-1]*2, [1]*2, 100)),\n",
    "            \"nguyen-10\": synthetic_gen(\"Nguyen-10\", rnd,\n",
    "                                       lambda x_1, x_2: 2*sin(x_1)*cos(x_2),\n",
    "                                       U([-1]*2, [1]*2, 100),\n",
    "                                       U([-1]*2, [1]*2, 100))\n",
    "            }\n",
    "\n",
    "seed = 4567\n",
    "rnd = np.random.RandomState(seed)\n",
    "# Loading real-world datasets\n",
    "data_real = {\"abalone\": get_data('Abalone', \n",
    "                                 'https://archive.ics.uci.edu/ml/'\n",
    "                                 'machine-learning-databases/abalone/abalone.data',\n",
    "                                 rnd=rnd, dataset='ABA'),\n",
    "            \"airfoil\": get_data('Airfoil',\n",
    "                                'https://archive.ics.uci.edu/ml/'\n",
    "                                'machine-learning-databases/00291/airfoil_self_noise.dat',\n",
    "                                pd_sep=\"\\t\"),\n",
    "            \"boston\": get_data('Boston', \"\", dataset=\"BOH\"),\n",
    "            \"combined-cycle\": get_data('Combined-cycle', \n",
    "                                       'https://archive.ics.uci.edu/ml/'\n",
    "                                       'machine-learning-databases/00294/CCPP.zip',\n",
    "                                       dataset=\"CCP\"),\n",
    "            \"computer-hardware\": get_data('Computer-hardware', \n",
    "                                          'https://archive.ics.uci.edu/ml/'\n",
    "                                          'machine-learning-databases/cpu-performance/machine.data',\n",
    "                                          dataset=\"CPU\"),\n",
    "            \"concrete-strength\": get_data('Concrete-strength', \n",
    "                                          'https://archive.ics.uci.edu/ml/'\n",
    "                                          'machine-learning-databases/concrete/compressive/Concrete_Data.xls',\n",
    "                                          dataset=\"CST\"),\n",
    "            \"energy-cooling\": get_data('Energy-cooling', \n",
    "                                       'http://archive.ics.uci.edu/ml/'\n",
    "                                       'machine-learning-databases/00242/ENB2012_data.xlsx',\n",
    "                                       dataset=\"ENC\"),\n",
    "            \"energy-heating\": get_data('Energy-heating', \n",
    "                                       'http://archive.ics.uci.edu/ml/'\n",
    "                                       'machine-learning-databases/00242/ENB2012_data.xlsx',\n",
    "                                       dataset=\"ENH\"),\n",
    "            \"forest-fire\": get_data('Forest-fire', \n",
    "                                    'https://archive.ics.uci.edu/ml/'\n",
    "                                    'machine-learning-databases/forest-fires/forestfires.csv',\n",
    "                                    pd_header='infer', dataset=\"FFR\"),\n",
    "            \"ozone\": get_data('Ozone', './data/ozone.data', dataset=\"OZO\"),\n",
    "            \"wine-quality-red\": get_data('Wine-quality-red', \n",
    "                                         'https://archive.ics.uci.edu/ml/'\n",
    "                                         'machine-learning-databases/wine-quality/winequality-red.csv',\n",
    "                                         pd_header='infer', pd_sep=';'),\n",
    "            \"wine-quality-white\": get_data('Wine-quality-white', \n",
    "                                           'https://archive.ics.uci.edu/ml/'\n",
    "                                           'machine-learning-databases/wine-quality/winequality-white.csv', \n",
    "                                           pd_header='infer', pd_sep=';'),\n",
    "            \"yacht\": get_data('Yacht', 'https://archive.ics.uci.edu/ml/'\n",
    "                              'machine-learning-databases/00243/yacht_hydrodynamics.data', dataset=\"YAC\")\n",
    "            }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Analyzing  GP performance on the datasets\n",
    "## Pre-running analysis : Comparing function times (used by GP)\n",
    "Just to have an idea of the time spent by functions (from a possible function set) inside GP. Notice that we analyzed the functions as they are implemented by `gplearn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-24T19:23:33.576526Z",
     "start_time": "2018-03-24T19:23:31.589629Z"
    },
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sin time:        0.1362\n",
      "Cos time:        0.1267\n",
      "Log time:        0.3099\n",
      "Sqrt time:       0.0646\n",
      "Prot. div time:  0.2343\n",
      "AQ time:         0.1254\n",
      "Mul time:        0.0354\n",
      "Sub time:        0.0366\n",
      "Add time:        0.0352\n"
     ]
    }
   ],
   "source": [
    "fmt = \"{:<17}{:.4f}\"\n",
    "print(fmt.format(\"Sin time:\",\n",
    "                 timeit.timeit(stmt=\"[np.sin(x) for x in X]\", \n",
    "                                      setup=\"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000)\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"Cos time:\", \n",
    "                 timeit.timeit(stmt=\"[np.cos(x) for x in X]\", \n",
    "                                      setup=\"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000)\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"Log time:\", \n",
    "                 timeit.timeit(stmt=\"[_protected_log(x) for x in X]\", \n",
    "                                      setup=\"from gplearn.functions import _protected_log;\"\n",
    "                                      \"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000)\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"Sqrt time:\", \n",
    "                 timeit.timeit(stmt=\"[_protected_sqrt(x) for x in X]\", \n",
    "                                      setup=\"from gplearn.functions import _protected_sqrt;\"\n",
    "                                      \"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000)\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"Prot. div time:\", \n",
    "                 timeit.timeit(stmt=\"[_protected_division(x,y) for x,y in zip(X,Y)]\", \n",
    "                                      setup=\"from gplearn.functions import _protected_division;\"\n",
    "                                      \"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000);\"\n",
    "                                      \"Y = np.random.rand(1,1000);\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"AQ time:\", \n",
    "                 timeit.timeit(stmt=\"[np.divide(x, np.sqrt(1+y**2)) for x,y in zip(X,Y)]\", \n",
    "                                      setup=\"from gplearn.functions import _protected_division;\"\n",
    "                                      \"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000);\"\n",
    "                                      \"Y = np.random.rand(1,1000);\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"Mul time:\", \n",
    "                 timeit.timeit(stmt=\"[np.multiply(x,y) for x,y in zip(X,Y)]\", \n",
    "                                      setup=\"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000);\"\n",
    "                                      \"Y = np.random.rand(1,1000);\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"Sub time:\", \n",
    "                 timeit.timeit(stmt=\"[np.subtract(x,y) for x,y in zip(X,Y)]\", \n",
    "                                      setup=\"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000);\"\n",
    "                                      \"Y = np.random.rand(1,1000);\", \n",
    "                                      number = 10000)))\n",
    "print(fmt.format(\"Add time:\",\n",
    "                 timeit.timeit(stmt=\"[np.add(x,y) for x,y in zip(X,Y)]\", \n",
    "                                      setup=\"import numpy as np;\"\n",
    "                                      \"X = np.random.rand(1,1000);\"\n",
    "                                      \"Y = np.random.rand(1,1000);\", \n",
    "                                      number = 10000)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Given that Log and Protected Division are the ones with largest execution time, we selected Addition, Subtraction, Multiplication, Analytic Quotient, Square Root and Sine as our function set.\n",
    "\n",
    "## gplearn: Configuration\n",
    "`gplearn` is Python library for GP, compatible with Scikit learn. We had to specify the Analytic Quotient function manually before the evolution and the parameters used by GP are hard-coded in the function `run_gp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def _aq(x1, x2):\n",
    "    return np.divide(x1, np.sqrt(1 + x2**2))   \n",
    "\n",
    "aq = make_function(function=_aq,\n",
    "                   name='aq',\n",
    "                   arity=2)\n",
    "\n",
    "def run_gp(training, test, rnd=None):\n",
    "    # Normalization constant used to compute the normalized RMSE\n",
    "    norm_const = np.sqrt(training.shape[0]/(training.shape[0]-1))/training.iloc[:,-1].std()\n",
    "    \n",
    "    # Genetic programming instance used to perform symbolic regression in the datasets\n",
    "    function_set = ['add', 'sub', 'mul', aq, 'sqrt', 'sin']\n",
    "    gp_param = {'population_size': 1000,\n",
    "                  'generations': 50, \n",
    "                  'stopping_criteria': 0.00,\n",
    "                  'const_range': (-1,1),\n",
    "                  'init_depth': (2,6),\n",
    "                  'init_method': 'half and half',\n",
    "                  'metric': make_fitness(function=lambda y, y_pred, w:\n",
    "                      np.sqrt(np.average(((y_pred - y) ** 2), weights=w)) * norm_const,\n",
    "                                           greater_is_better=False),\n",
    "                  'tournament_size': 10,\n",
    "                  'p_crossover': 0.85, \n",
    "                  'function_set': function_set,\n",
    "                  'p_subtree_mutation': 0.05,\n",
    "                  'p_hoist_mutation': 0.05, \n",
    "                  'p_point_mutation': 0.05,\n",
    "                  'max_samples': 1, \n",
    "                  'verbose': 0,\n",
    "                  'parsimony_coefficient': 0.001, \n",
    "                  'random_state': rnd,\n",
    "                  'n_jobs': 1}\n",
    "    est_gp = SymbolicRegressor(**gp_param)\n",
    "    est_gp.fit(training.iloc[:,:-1], training.iloc[:,-1])\n",
    "    \n",
    "    #norm_const = np.sqrt(training.shape[0]/(training.shape[0]-1))/training.iloc[:,-1].std()\n",
    "    stats = {'size': est_gp._program.length_,\n",
    "             'tr_rmse': est_gp._program.raw_fitness_ / norm_const,\n",
    "             'tr_nrmse': est_gp._program.raw_fitness_}\n",
    "    y_est = est_gp.predict(test.iloc[:,:-1])\n",
    "    rmse = np.sqrt(np.average((y_est - test.iloc[:,-1]) ** 2))\n",
    "    norm_const = np.sqrt(test.shape[0]/(test.shape[0]-1))/test.iloc[:,-1].std()\n",
    "    \n",
    "    stats['ts_rmse'] = rmse\n",
    "    stats['ts_nrmse'] = rmse * norm_const\n",
    "    \n",
    "    return stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "We implemented the experiments so that they run in parallel, taking advantage of multi-core processors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def evaluate(dataset_lst, n_jobs=None, n_rep=30, cv=None):\n",
    "    stats = {}\n",
    "    experiments = []\n",
    "    seed = 1234\n",
    "    rnd = np.random.RandomState(seed)\n",
    "    kf = KFold(n_splits=cv, shuffle=True, random_state=rnd)\n",
    "    \n",
    "    for data_name, data_points in dataset_lst.items():\n",
    "        stats[data_name] = {}\n",
    "        # Reset the pseudo-random number generator for each dataset (used by GP)\n",
    "        rnd = np.random.RandomState(seed)\n",
    "        \n",
    "        if cv is not None: \n",
    "            partitions = list(kf.split(data_points))\n",
    "        \n",
    "        for i in range(n_rep):\n",
    "            if cv is not None: \n",
    "                training = data_points.iloc[partitions[i % cv][0],:]\n",
    "                test = data_points.iloc[partitions[i % cv][1],:]\n",
    "            else:\n",
    "                training = data_points['training']\n",
    "                test = data_points['test']\n",
    "            experiments.append([data_name, i+1, training, test, rnd.randint(10**6)])\n",
    "    if not n_jobs or n_jobs > 1:    \n",
    "        with Pool(processes=n_jobs) as pool:\n",
    "            for exec_stats in pool.imap(worker, experiments):\n",
    "                data_name = exec_stats[0]\n",
    "                for stats_name, stats_value in exec_stats[1].items():\n",
    "                    if stats_name not in stats[data_name]:\n",
    "                        stats[data_name][stats_name] = []\n",
    "                    stats[data_name][stats_name].append(stats_value)\n",
    "    else:\n",
    "        for experiment in experiments:\n",
    "            exec_stats = worker(experiment)\n",
    "            data_name = exec_stats[0]\n",
    "            for stats_name, stats_value in exec_stats[1].items():\n",
    "                if stats_name not in stats[data_name]:\n",
    "                    stats[data_name][stats_name] = []\n",
    "                stats[data_name][stats_name].append(stats_value)\n",
    "    return stats\n",
    "\n",
    "def worker(param):\n",
    "    print(\"Dataset \\\"\" + param[0] + \"\\\", exec.\", param[1])\n",
    "    return [param[0], run_gp(training=param[2], test=param[3], rnd=param[4])]\n",
    "\n",
    "def write_stats(stats_dict, path='./performace_metrics/'):\n",
    "    \"\"\"Write statistics returned by worker/run_gp methods.\n",
    "\n",
    "    Extended description of function.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    stats_dict : dict\n",
    "        A dictionary containing results from different datasets. The format is\n",
    "        in the form stats_dict['data_set_name']['statistics_name']\n",
    "    arg2 : str\n",
    "        Description of arg2\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    bool\n",
    "        Description of return value\n",
    "    \"\"\"\n",
    "    for key, data in stats_dict.items():\n",
    "        df = pd.DataFrame(data)\n",
    "        df.to_csv(path + key + '.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "## Running GP for each dataset\n",
    "This stage should take a few hours. That is why the results are saved in a local folder for future accesses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "%%script false  # Remove this cell magic to run GP again\n",
    "\n",
    "stats_synt = evaluate(data_synt, n_rep = 30)\n",
    "write_stats(stats_synt)\n",
    "\n",
    "# The real-word dataset is evaluated using cross-validation\n",
    "stats_real = evaluate(data_real, n_rep = 30, cv = 5)\n",
    "write_stats(stats_real)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extracting Meta-Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-18T20:48:26.921552Z",
     "start_time": "2018-04-18T20:48:24.809604Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading the vladislavleva-8 dataset\n",
      "Loading the vladislavleva-1 dataset\n",
      "Loading the koza-3 dataset\n",
      "Loading the keijzer-6 dataset\n",
      "Loading the koza-2 dataset\n",
      "Loading the nguyen-4 dataset\n",
      "Loading the vladislavleva-2 dataset\n",
      "Loading the keijzer-1 dataset\n",
      "Loading the r2 dataset\n",
      "Loading the nonic dataset\n",
      "Loading the korns-1 dataset\n",
      "Loading the meier-3 dataset\n",
      "Loading the nguyen-9 dataset\n",
      "Loading the keijzer-9 dataset\n",
      "Loading the korns-7 dataset\n",
      "Loading the nguyen-5 dataset\n",
      "Loading the nguyen-6 dataset\n",
      "Loading the nguyen-2 dataset\n",
      "Loading the sine dataset\n",
      "Loading the keijzer-15 dataset\n",
      "Loading the nguyen-1 dataset\n",
      "Loading the vladislavleva-4 dataset\n",
      "Loading the vladislavleva-6 dataset\n",
      "Loading the keijzer-2 dataset\n",
      "Loading the nguyen-7 dataset\n",
      "Loading the poly-10 dataset\n",
      "Loading the nguyen-10 dataset\n",
      "Loading the keijzer-10 dataset\n",
      "Loading the keijzer-8 dataset\n",
      "Loading the meier-4 dataset\n",
      "Loading the keijzer-12 dataset\n",
      "Loading the nguyen-8 dataset\n",
      "Loading the keijzer-13 dataset\n",
      "Loading the keijzer-14 dataset\n",
      "Loading the vladislavleva-5 dataset\n",
      "Loading the burks dataset\n",
      "Loading the keijzer-7 dataset\n",
      "Loading the r3 dataset\n",
      "Loading the keijzer-4 dataset\n",
      "Loading the keijzer-5 dataset\n",
      "Loading the korns-11 dataset\n",
      "Loading the r1 dataset\n",
      "Loading the vladislavleva-3 dataset\n",
      "Loading the korns-4 dataset\n",
      "Loading the korns-12 dataset\n",
      "Loading the keijzer-11 dataset\n",
      "Loading the vladislavleva-7 dataset\n",
      "Loading the keijzer-3 dataset\n",
      "Loading the nguyen-3 dataset\n",
      "Loading the pagie-1 dataset\n",
      "Loading the boston dataset\n",
      "Loading the concrete-strength dataset\n",
      "Loading the energy-cooling dataset\n",
      "Loading the forest-fire dataset\n",
      "Loading the combined-cycle dataset\n",
      "Loading the ozone dataset\n",
      "Loading the energy-heating dataset\n",
      "Loading the abalone dataset\n",
      "Loading the airfoil dataset\n",
      "Loading the yacht dataset\n",
      "Loading the computer-hardware dataset\n",
      "Loading the wine-quality-red dataset\n",
      "Loading the wine-quality-white dataset\n"
     ]
    }
   ],
   "source": [
    "# Geta data from file\n",
    "names_synt = ['vladislavleva-8', 'vladislavleva-1', 'koza-3', 'keijzer-6', 'koza-2', \n",
    "              'nguyen-4', 'vladislavleva-2', 'keijzer-1', 'r2', 'nonic', 'korns-1', \n",
    "              'meier-3', 'nguyen-9', 'keijzer-9', 'korns-7', 'nguyen-5', 'nguyen-6', \n",
    "              'nguyen-2', 'sine', 'keijzer-15', 'nguyen-1', 'vladislavleva-4', \n",
    "              'vladislavleva-6', 'keijzer-2', 'nguyen-7', 'poly-10', 'nguyen-10', \n",
    "              'keijzer-10', 'keijzer-8', 'meier-4', 'keijzer-12', 'nguyen-8', \n",
    "              'keijzer-13', 'keijzer-14', 'vladislavleva-5', 'burks', 'keijzer-7', \n",
    "              'r3', 'keijzer-4', 'keijzer-5', 'korns-11', 'r1', 'vladislavleva-3', \n",
    "              'korns-4', 'korns-12', 'keijzer-11', 'vladislavleva-7', 'keijzer-3', \n",
    "              'nguyen-3', 'pagie-1']\n",
    "\n",
    "names_real = ['boston', 'concrete-strength', 'energy-cooling', 'forest-fire', \n",
    "              'combined-cycle', 'ozone', 'energy-heating', 'abalone', 'airfoil', \n",
    "              'yacht', 'computer-hardware', 'wine-quality-red', 'wine-quality-white']\n",
    "\n",
    "data_synt = {}\n",
    "data_real = {}\n",
    "for k in names_synt:\n",
    "    print('Loading the ' + k + ' dataset')\n",
    "    data_synt[k] = {'training': pd.read_csv('data/' + k + '_tr.csv'),\n",
    "                    'test': pd.read_csv('data/' + k + '_ts.csv')\n",
    "                    }\n",
    "\n",
    "for k in names_real:\n",
    "    print('Loading the ' + k + ' dataset')\n",
    "    data_real[k] = pd.read_csv('data/' + k + '.csv')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-18T20:48:31.984863Z",
     "start_time": "2018-04-18T20:48:29.852958Z"
    }
   },
   "outputs": [],
   "source": [
    "def get_linear_r2(df):\n",
    "    \"\"\"Computes the coefficient of determination of a linear model\n",
    "    \n",
    "    The model is fitted using the `n_columns-1` first columns as input for the model and the \n",
    "    last column as output\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    df: pandas.DataFrame\n",
    "        Input dataset, where the `n_columns-1` first columns correspond to the features and \n",
    "        the last column corresponds to the response variable.\n",
    "    \"\"\"\n",
    "    lr = LinearRegression()\n",
    "    lr.fit(df.iloc[:, :-1], df.iloc[:, -1])\n",
    "    return lr.score(df.iloc[:, :-1], df.iloc[:, -1])\n",
    "\n",
    "def get_meta_data(data_set_name, training, test):\n",
    "    if not isinstance(training, (list, tuple)):\n",
    "        training = [training]\n",
    "    if not isinstance(test, (list, tuple)):\n",
    "        test = [test]\n",
    "        \n",
    "    meta_data = {}\n",
    "    meta_data['N. instances (Training)'] = np.mean([fold.shape[0] for fold in training])\n",
    "    meta_data['N. features'] = np.mean([fold.shape[1] - 1 for fold in training])\n",
    "    meta_data['Target std (Training)'] =  np.mean([fold.iloc[:,-1].std() for fold in training])\n",
    "    meta_data['Target skewness (Training)'] = np.mean([fold.iloc[:,-1].skew() for fold in training])\n",
    "    meta_data['Linearity measure (Training)'] = np.mean([get_linear_r2(fold) for fold in training])\n",
    "    training_corr = [fold.corr() for fold in training]\n",
    "    meta_data['Mean absolute correlation attribute-target (Training)'] = np.mean([fold_cor.iloc[-1, :-1].abs() for fold_cor in training_corr])\n",
    "    \n",
    "    meta_data['N. instances (Test)'] = np.mean([fold.shape[0] for fold in test])\n",
    "    meta_data['Target std (Test)'] = np.mean([fold.iloc[:,-1].std() for fold in test])\n",
    "    meta_data['Target skewness (Test)'] = np.mean([fold.iloc[:,-1].skew() for fold in test])\n",
    "    meta_data['Linearity measure (Test)'] = np.mean([get_linear_r2(fold) for fold in test])\n",
    "    test_corr = [fold.corr() for fold in test]\n",
    "    meta_data['Mean absolute correlation attribute-target (Test)'] = np.mean([fold_cor.iloc[-1, :-1].abs() for fold_cor in test_corr])\n",
    "    \n",
    "    gp_output = pd.read_csv('performace_metrics/' + data_set_name + '.csv', index_col=0)\n",
    "    meta_data['Median NRMSE  (Training)'] = gp_output['tr_nrmse'].median()\n",
    "    meta_data['Median NRMSE (Test)'] = gp_output['ts_nrmse'].median()\n",
    "    meta_data['Function size'] = gp_output['size'].median()\n",
    "    return meta_data\n",
    "\n",
    "def get_cross_validation_sets(data_points, kfold, cv=5):\n",
    "    training = []\n",
    "    test = []\n",
    "    partitions = list(kf.split(data_points))\n",
    "    for i in range(cv):\n",
    "        training.append(data_points.iloc[partitions[i][0],:])\n",
    "        test.append(data_points.iloc[partitions[i][1],:])\n",
    "    return training, test\n",
    "\n",
    "seed = 1234\n",
    "rnd = np.random.RandomState(seed)\n",
    "kf = KFold(n_splits=5, shuffle=True, random_state=rnd)\n",
    "\n",
    "meta_info = {}\n",
    "for key, data in data_real.items():\n",
    "    training, test = get_cross_validation_sets(data, kf)\n",
    "    meta_info[key] = get_meta_data(key, training, test)\n",
    "    meta_info[key]['Nature'] = 'Real'\n",
    "\n",
    "for key, data in data_synt.items():\n",
    "    meta_info[key] = get_meta_data(key, data['training'], data['test'])\n",
    "    meta_info[key]['Nature'] = 'Synthetic'\n",
    "\n",
    "meta_info = pd.DataFrame(meta_info).transpose()\n",
    "# Impute NaN's with 0 (NaN's appear when a feature has no variation)\n",
    "meta_info['Mean absolute correlation attribute-target (Training)'].fillna(value = 0, inplace = True)\n",
    "meta_info['Mean absolute correlation attribute-target (Test)'].fillna(value = 0, inplace = True)\n",
    "# Cast the columns to numeric types\n",
    "for k in meta_info.keys():\n",
    "    if k != 'Nature':\n",
    "        meta_info[k] = pd.to_numeric(meta_info[k])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-18T21:15:19.625077Z",
     "start_time": "2018-04-18T21:15:06.145806Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 5 folds for each of 10 candidates, totalling 50 fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  50 out of  50 | elapsed:   13.1s finished\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Random forest best hyperparameters: {'n_estimators': 120}\n",
      "\n",
      "Feature importance:\n",
      "Target skewness (Training)                               0.290170\n",
      "N. instances (Training)                                  0.184041\n",
      "N. features                                              0.122526\n",
      "Linearity measure (Training)                             0.075096\n",
      "Target skewness (Test)                                   0.056585\n",
      "Target std (Training)                                    0.054541\n",
      "Linearity measure (Test)                                 0.052174\n",
      "Target std (Test)                                        0.051961\n",
      "Mean absolute correlation attribute-target (Test)        0.044588\n",
      "Mean absolute correlation attribute-target (Training)    0.037164\n",
      "N. instances (Test)                                      0.031154\n",
      "dtype: float64\n"
     ]
    }
   ],
   "source": [
    "X = meta_info.drop(['Median NRMSE  (Training)', 'Function size', \n",
    "                    'Nature', 'Median NRMSE (Test)'], axis = 1)\n",
    "y = meta_info['Median NRMSE (Test)']\n",
    "\n",
    "# Tuning the number of trees in the forest\n",
    "n_estimators = [int(x) for x in np.linspace(start = 10, stop = 1000, num = 10)]\n",
    "\n",
    "random_grid = {'n_estimators': n_estimators}\n",
    "\n",
    "seed = 789\n",
    "rnd = np.random.RandomState(seed)\n",
    "\n",
    "rf = RandomForestRegressor(random_state=rnd)\n",
    "rf_random = RandomizedSearchCV(estimator = rf, \n",
    "                               param_distributions = {'n_estimators': n_estimators}, \n",
    "                               n_iter = 10, cv = 5, verbose=1, \n",
    "                               random_state=rnd, n_jobs = -1)\n",
    "\n",
    "rf_random.fit(X, y)\n",
    "print('Random forest best hyperparameters:', rf_random.best_params_)\n",
    "\n",
    "rf = rf_random.best_estimator_\n",
    "\n",
    "feat_importance = pd.Series(rf.feature_importances_, X.columns.values)\n",
    "feat_importance.sort_values(inplace=True, ascending=False)\n",
    "print(\"\\nFeature importance:\") \n",
    "print(feat_importance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-02T11:27:40.834127Z",
     "start_time": "2018-04-02T11:26:58.728475Z"
    }
   },
   "outputs": [],
   "source": [
    "# Prevent Type 3 font errors in GECCO submission\n",
    "matplotlib.rcParams['pdf.fonttype'] = 42\n",
    "matplotlib.rcParams['ps.fonttype'] = 42\n",
    "\n",
    "# Use default seaborn pallete\n",
    "sns.set(color_codes=True)\n",
    "\n",
    "plot_keys = list(meta_info.keys())\n",
    "[plot_keys.remove(k) for k in ('Median NRMSE  (Training)', 'Nature', 'Function size')];\n",
    "\n",
    "sm = cm.ScalarMappable(cmap=cm.get_cmap('RdYlBu'))\n",
    "sm.set_array(meta_info['Median NRMSE (Test)'])\n",
    "\n",
    "colors = sm.to_rgba(meta_info['Median NRMSE (Test)'])\n",
    "markers = meta_info['Nature'].replace(['Real', 'Synthetic'], ['o', '^'])\n",
    "\n",
    "for i in range(len(plot_keys)):   \n",
    "    for j in range(i+1, len(plot_keys)):\n",
    "        markers = meta_info['Nature'].replace(['Real', 'Synthetic'], ['o', '^'])\n",
    "        \n",
    "        for k in range(meta_info.shape[0]):\n",
    "            plt.scatter(meta_info[plot_keys[i]][k],\n",
    "                        meta_info[plot_keys[j]][k], \n",
    "                        c = colors[k],\n",
    "                        marker = markers[k],\n",
    "                        linewidth=0.2, \n",
    "                        edgecolor='black')\n",
    "        \n",
    "        plt.xlabel(plot_keys[i])\n",
    "        plt.ylabel(plot_keys[j])\n",
    "        cbar = plt.colorbar(sm)\n",
    "        cbar.ax.set_ylabel('NRMSE (Test)')\n",
    "        \n",
    "        plt.savefig('plots/' + plot_keys[i] + ' VS ' + plot_keys[j] + '.pdf', \n",
    "                    bbox_inches=\"tight\")\n",
    "        plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-02T11:34:26.800115Z",
     "start_time": "2018-04-02T11:34:22.019251Z"
    }
   },
   "outputs": [],
   "source": [
    "# Ploting the distribution of the meta features\n",
    "for k in plot_keys:\n",
    "    #meta_info[k].plot(kind='hist', bins=25, alpha=0.5, color='blue')\n",
    "    df_aux = pd.DataFrame({'Real': meta_info[meta_info['Nature'] == 'Real'][k], 'Synthetic': meta_info[meta_info['Nature'] == 'Synthetic'][k]}).plot(kind='hist', bins=25, stacked=True)\n",
    "    plt.title(k)\n",
    "    plt.savefig('plots/Hist. ' + k + '.pdf', \n",
    "                bbox_inches=\"tight\")\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-01T22:37:09.207209Z",
     "start_time": "2018-04-01T22:37:08.304332Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.22900379572563082\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAADtCAYAAADZRzznAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXl4nHd19/25t9k3rZZt2ZZtWV5iO7GzxyYOMSSBLAWyEWhYn6YN21u4KA/Nm0IpbSkNdGMrhKdP3gItKQFKQkJwQkK2QuLETrzLqyxrtaTRSLPP3Mv7xz0zGo1G0mySxs58r8tXopn7/t2/mbnnzDnne873CIZhGNRQQw011DAtxIXeQA011FBDtaNmKGuooYYaZkHNUNZQQw01zIKaoayhhhpqmAU1Q1lDDTXUMAtqhrKGGmqoYRbIMz05NBScr33UUEMN5ziamtxlr7H18m8XfOyelz9W9vUKxYyGsoYaaqhhPiGIwkJvIS9qhrKGGmqoGghSzVDWUEMNNcwMsTppk5qhrKGGGqoGNY+yhhpqqGEWiErNo6yhhhpqmBm10LuGGmqoYWbUQu8aaqihhllQKw+qoYYaapgFNY+yhhpqqGEW1DzKGmqooYZZICjSQm8hL2qGsoaSYRgGhmEgVilTWcO5h1roXcN5AcMwEAQADV3XEAQBSbJgGGT+1VBDqaiF3jWc09B1HdAxDA2Xy4FhSITDUQRBQBBIGU8mGcya0ayhWNQ8yhrOOUwYRx0wUt7kzDdyPqMJArpes5o1FIAqTePUDGUNkzCTccw2krMZTPMYEEUBj8fB6Gg4tf5c7LqG8wW1FsYaqhYmKaPNahzLQdpREISJkLxmNGvIRS1HWUNVwTSOE95jqcbRYlEQBDHlic6O7NC8ZjRrmAKp5lHWsMAwjaGBrmtlGUdJEpFlGYfDhqpqyLKErhskEkkSCRWjQBYnn9GskUBvbtQ8yhoWDKa3ZxrHpqY6RkbGMIzijKOiyNhsVmw2K2Cgqhqjo2OAgKIoyLKExSLj8TjRNI1EQiWZTBa8fi4JBAKCIKCqNVfzzYQa613DvCK7nKdUSJKE3T5hHKPROH5/AJvNiiCYTLaY8gBUVUNVNSCOoshYLDIOh5VksvjrCwIoiml4Q6FYzdN8E6HmUdYw5yiknMf0JKdfQxRFbDYrdrsVURSIxRIEAuMpI1gYkkmVZFJFEMwcpqJY8XpdJJMqiUSyqLXylRvV8pnnL8RaC2MNc4FCy3lmgiAI2GwWbDYriiITiyUIBsMkEoWHzvlgGKbRtNksjI9HsFjkjAFOJEyjqWmFW7200ayRQOcvxJpHWUOlkGaszbC69HIeq9WC3W7FYlFIJJJEIjHi8cQc7NfcczyeJB5PIooiFouM02kHjIzRzC5KzzaGuagx5+cvaoayhrJQqXIeQQCXy4nVqpBMqsRiccbGQgUz1ZWAruvEYglisQSSJGKxKLjdjknMeaHINZqCICCKIolE6bnZGhYONTKnhqKRbRxdLjvxeIJEwnSbijGOsixht9uw2SyIooimaQwPhwuufZxLaJpONBonGo2nmHMFj8eKrutFG+9sEkjTtBoJdA6i5lHWUDCyy3nSEEWhKDkzSZogZUAgFovj94/j87mJxeJVYSRzMcGcg91u5kt9vgkSqFAGPU1Y1YQ6zj3UDGUNMyKXlMnFbGw1mDdZutZRkqRMWJ1MFh7KVgs0TUcQVKLROIqiYLVacDjElNFUp2XO871HNaN57kCWa6x3DTkoTp3HAKY+JwhChpRRFJl4PEEoFCmbsV5oZBs2M2+ZRBTN4vYJ5tzMZ05mzgXy/dBkr5tNANXKjaoLYi1HWQOUXs6T61FarWY5j9VqMtbRaIzR0UIY6/wGtxqR6/HpukE8niAeT0zLnBeKyQSQeS2LRSYWO/e87/MJtdD7TYxKlfNIkojH48Rms6KqGtFonPHx4hjrQkL4cwGTmXMJq1XG7XakPEQdQRCK6jkXRQGn00oyqdbKjRYQtc6cNxnSxtFuV1LkiVaScZRlKUPKGAZEIlGGhwNVScbkx9zf+JqmEYmY7ZNmCkLC43GiqlqKBJrdS0x7lbUazYVFzaN8EyC31hHAZnOSTCaLMmyiKGZ6rEVRIBqNE4vF0TSdcDg6R7ufS8xvjWYyCdGo2QlktSo4HLYC2ien5jZr6kbzD7km3Hv+Il85TxqFjE+AdBuh6TnKsslYj49PMNYul6PS265ylOZZZIfcZt5SRRCEWdsnZ+oESj9fY87nHjWP8jzDbOU8acyWEzR7rG1YLDLxeJJwOJq3jbBQg1sIKrXOXMI0XJWxQpPbJwUsFgWn0wYIRXcCpfeWbTQlSUTTDDStZjXLRbWOPq4ZyiJQyrCtfCyzxWKWuFitFpJJs1ZwbCw4L22ElbyGKAqZgvZiFIEWErpu5GmftGc8w2JIIPN488cuXd9ZKzcqDzUy5xxFueo8aY8yW/hW0zRisTjBYLjg6YSGUR1hiSCYBJOimLm/RCKJLMuIokA8bnpn2fnYanZes9snrVYFm80yRXi4EJuZW2pUI4FKRzXc4/lQM5R5ULlyHinVv+zCMPSM8G0x0mJZuyrhnMoh2wvWdR1VVQkEgoiiiKIoiKKI1Wp6Z7niFvOZxzMFhYt/fzVNR9N0QqFojvCwmjKa04fnuV5ojTkvHbWC8ypHNmNtsUjY7XYCgSBQnHHMFb7VdYNIJEY4HClzf/PvnaUVzu12K7puEI3GCAbD2O22Ke+JrucXt9A07RzKiZr/ny08rCjFMOf5160ZzcIhy5XJUSaTSe677z56e3tJJBLce++97Ny5M/P8M888w7e+9S1kWebWW2/ljjvumHlfFdnVOYp85Tzm41KR0mXTC9+63c68bHgJu6USNYmzzcpJ94vb7bZMadLoaHEK59niFjabBatVSimcJ2fs015YTC0PMtsnZ2fOTY+ygCvkMZqiKNbmAmWhUqH3o48+is/n44EHHiAQCPCud70rYyiTySRf+cpXeOSRR7Db7dx1111ce+21NDY2Trvem9JQzlTOA4UzzIUL31aTR5X/Gz2Zfa+MwjmY3pmiyIRC0Uw4O8E2J6fN0aZn8hSL2cp8SsVMwsPmILTi9ykI4PU6GBsLoeu1ciMAoUKs9w033MD1118PmJ+dJE2IbZw4cYLly5fj9XoBuPjii9m9ezfveMc7pl3vTWMoCy3ngZkNpcWipEgZS0HCt+ZaZW6euQm9FUXO6FROsO/jc25ophPrnfoezp/VKMbA5goPu1x2nE47uq6XNLIXIG0f3uw1mpXKUTqdTgBCoRCf+tSn+NM//dPMc6FQCLfbPenYUCg043rntaEsrZxnaniaLXybztUND0cKIg0qV/9YmdA7Xdju8ThTryU+7y2RuWK9JuNsRVVnJ05mx8zqQZU+T9N0DANCoXDK05zIzZZTo/lmVTeqJOvd39/Pxz/+cd73vvdx8803Zx53uVyEw+HM3+FweJLhzIfzzlBWYthW2gt0Ou1ThG81bWHya+V4lIIgZFoi01JsxeYd5wrZ+czslkOTWS/NpSrFEysnZE+fm/1aTOZ8MglUzOjefOpGcP4bzUoZyuHhYT7ykY/whS98gSuvvHLSc6tXr+b06dMEAgEcDgevvvoqH/3oR2dc77wwlKZhM1AUkVgsWXI5zwSRYUUURURRLFv41jCMBes2yM6hpnUq7XYb8Xi8KoxkLrKJE6fTlhp1q8yaz8zGwhDsUy86mTnPLzxcaHH7m4k5l6XKfFf+9V//lfHxcb797W/z7W9/G4Dbb7+daDTKnXfeyec//3k++tGPYhgGt956K4sWLZp5XxXZ1QIgl7GWJAmn00ksFi+asc4nfOvzeQgGw7MvUNA1yl+j0BDezDua3uNE3nEih5r2kKsZhmGk8n0qmqZNymfG48mCC8GLQ6kh+8ztljMJD5fyA5xtNOvqnIyOmvfo+WI0K+VR3n///dx///3TPn/ttddy7bXXFrzeOWUopyvnSaPYUa2lCd8WB/MLNLeGSZLEVA7VimEYxGLzn3ecK2TnM83BYaahmT6fWZrBmyu2PBu5wsNWq4IgCHg8jrwje2eDOXEy/f/nBwkk1TpzSsds5TxQGLuc7W2VKnxbLOaqUDxbbSg9HycQKCzvON1+vvmtH/GJj7+/wjutHJJJjWRSywpnFRwOa05N47lhKNIMuSxLRCLxjPCwruvE42pKQLi4nvPzQd2o1sJYJIop54HpQ9Ns4VtTEKEwbyu9XvlGtDKsd3o/k/OOScLhCPF44fWO072cF1/awxNPhtly0V62bdsy6zpznwuc3jPMDWdNNSB7+tmS3u9SWx/Lg/kas4WH016z2T45vfDwTPdmPqN5rgQXYpV2cVWVoSy1nAcmG8p8wrfFsrzpEqFyDWUlftEVxSxuTv83Go3NWLtZCv77F/tB9PHzX+yb1VAWd9m5vfFz1YCcTjsOh61oYYtSUc49ks/7TXvNwIzCw4V6ztlCHVD9RrNSLYyVRlUZSsMonbFOH1tf780rfFv8XipVKF6ahzMxl9sGmCSGqmr4/WPl7ohc4/XiS3s4eMTsLjl4BF56qTCvstpg1jTqhMOJlKeZLWwxc3nOQoTssxnZmdonNU0rOjRP/1cQyHQCVRuqVI6yugwlFC8qm916B0wrfFssKimUWyjSPeN2uw1ZljI6lcmkmhGamAs8/8I+li3VADOEf+6FNypmKM23cD7JFTOcNT2zifIcm82Kw5Ff3bwczIeBzSc8nI6WbDbLFGm72VCl0S0AUpVursoMZWGM5XTCt42NvjK7OrJRKRGK2Q1uLgOfz9hXihTKt859n/9Q+QvPcs2FwvT5THPEbTyezKpMmL+OHigtbE+nGlRVw2azIAjCFGm7+RCAnivUyJwyUYjwbSW9wMoapnwk0wQDr2ka0WhsFgZ+7suMioFZFG7PlCSlUwPVjNx8ptWqZIR6BQGSJWiAVKKjp7RzTSM7WdpOLlh4uFptac1QloC0HqLNZgWMWYVvK2soK9mjbSJNMqWLvqPRWFFCvtUQlZghrA1JEolGY4RCkQwTP53K+XyiUOOjaTqRSByTaTZzmU6nSZoUZ/TL+VBKPzf3dU60T8ZnFR6uViMJNda7IKR/JR0Oe+aLF4slCq4PNCo4LqGSRlcURerqPCmtynhJbZGVvbmLe1253T6JhBmuBoPhVHumkJEfy1Y5j8eLG9M7dY/zEwonkyqapmQKwWea1jjlajN05cy607LOnf51ziY8XEzP+XyjxnrPgGAwyG9/+xuefvpX+P1+fvrTn5ekh6jrlWGqofzQ28yj2jLdF9NrVRa8IyqXM539uMneL5NUhhwOO1KentxclXOrVUFRzNBclqWilcHn2/MxDPJMazTzmfH4xA9E7j5LRamam9n7ne35fMy5wwF+f3mK+3OFc6IzJ1c+/f3v/xDbt+/IPP/wwz/iscd+gc/nA+Bzn7uP5cvbyt7EQw89yMBAP+9+921cdtnlJfdYL3ToPSHHls47xgmFwtTX+8pm4ueqwycXada9HO8XJkJBWZZwOKw5oXlxrXrFoFIGNjufOTX/ly2fVm6t7fwQQdnMeTXlunNxToTeufLpN998yyRD2dl5hPvv/xLr1q2v6CY++cnPAKDrKuZAr9JQ2ZKewtbKLm4XhLQc20TeURCKV7+eb0wV0pit7714pjYUimaF5hMCF5VQUa8EZjI8U/N/ZrSQTJq52IUic0pNa1QzK35OkDlT5dMnR+adnYf54Q//LyMjI1x11XbuvvvD87fTAlBp1nu6X15BAKvVmlEcisUSMxS3V7aFsVKYKbSeK+QLzdMCF+cCaw7Z+T8h1TljQRQFDMNadE620l09haKK7eS5YShz5dP/6I/unXTwzp3X8Z733IHT6eS++z7LSy+9wLZtb5m/3c6Cyofekx/Lrt+ceUZO9joV2U7FIMsSsmy2QpYTWpeLtJcmCBPva3ZoXjrmRz0oHcqC2UVlGMY81zOWXr9ZzTgnDCVMlk/fseP6zOOGYXDHHe/D5XIBcOWV2zl2rLPChlIoy9jNRY7SzDuaYWnaIxofDy9I+FLqa8sOrbNnclcDJhMoE6E5CChKcQQQlN/RUwoMg5x85sQ4CFM/M/8P0UJ5lNWMSgn3VhqTdpWWT/+zP/szbrvttkkHhsNhPvCBO4lEIhiGwZ49u1m7dt28bnY2VMpQmgKr5uCtujoPhgF+/zgjI2NEIrEFy/EUc11RNAUiGht9eL1uNE1neDhANFqd6uYwEZqPjYVJizF7vS4cDrOlsxqRa+xUVSMSMUVL4vEkVqsy7WuoRMH5+QZJFAr+N5+Y5FHmyqcnEiq33PJuotEof/AH7+Geez7Gpz71JyiKwiWXXMaVV26f183OhnIMZW7eUVU1kslkRTyvykm2zY7JrHViwULrcmEYZt5U1/W8oXmlWfO58NBy85m5Rfnz3f54LqBaWW/BmOHdHhqa3/BM13V0PVGysTP7pS2MjRW+76l94zFisUSm93psbOYxloWgqamekZHRsr/czc31DA2NTvmC5Bv/EIvF867hdNoRBIFQqPQ6OrMzR8oUnCvK9GIdph6ohVAoWtQ13G4HkUhsUrF3esytxSJnBG5z85k+n4tAoPjPzOt1pVpIizvP1DnVC9YETacXLBY5VSWRyOo3Lxz53p9CoetzI7XW1DTzJMNC8JPOswUfe/va5rKvVyiqouA8DbOUpvRflEI9ylwx32g0NqVvvFJsdXqtSiff87VDFiZIPL9SVpV0EKYbC5Gr1VjqPufDQctm/r1eJ5Ik4vE4UdXpRXrz4bz1KKszRVldhrJczNR1MjFh0VaQmO9M5UHF76tyAhumrJy1jNDaICc1PeeYi+/zxFiIyWEtmJ/1XBW056JcAxuJxDGMWI5IbzIzqXGurlutkM8V1vtcRj6PMluvMh5PFNwaWSnh3un2VQzSobUoCtjtNiKRGKOj45XZXJWiUEOQ3XEiSSJut2PS7JlCS41K/3wqk2fMbTV0OMxIYbpRveetR1mlOcqqMpTlhrppg5RmrG22bL3K8aLr5OZbuDcb+UJrXdcZGwtWTHR2fjB/76HZJWMwNhYuOjSvljKdXMOfPar3fNCbnA01QznHMCXZzBIMj8dVdqfJQvWNm5JlpgecG1qbYyHORZSqAlT6OdOF5pVnzeeu8Ds7J5tdn6mq5VUxVLOdrdLI+9w2lIIgZNje9MhWwzAYGQmUvXalRShmMpRpEd+09zO9B1y+glAlX5ckSTgcNgyDOSpBKv8bPZ2Hpml6RkZtIcR3iz1vot+clCqTjNfrIplMEo+bM3TOB9RGQRSM2X+hJ49sTRAKRTK5KIejUl5XZT3KXJjkkg273RTTKIS1royRK8/YplMbFosFm81CPJ7M5NQqPY+m+L3NbHwms+ZyxtM0PbT5drNK90STSRWr1UIwGElJwc2czzyXUPMoy0S+WsFKj2zNRuXHSphr5YbW5UyKnE9k5301Tc8U46frKNP6jS5Xdr9zdSgD5cPkYnDTQ/N4nEUbm1JJlUp05ZjecoJ4PDEln2m2TuYfBVHNobdcpfVBVW0oJUnM6DsahkEsNnvesVJdMJUsDxIEM79ohkqlkUvmnuaXYEqz7GmFoUgkzvDwaCbVkY1c/UarVcnocpaC0gxJ8V6aYRgkk2arYTgczQnNs3UnK4tKs9a5+Uyz1Cj/KIhqRtV6lKqqIsvVYy9DoSDd3V1s374tk3csdBQETHhv5RvK8sqDskNrSRKJx5NzLmNWCAoJ37NTG6XUamYrA5ktoUrGWyulC6VQlCuIkT80t83Imi9Emc5s18zOZ1osSsZopqdOQvVWTVQt6y1Wiav7xhuv8/DDP+K113Zz/fXXs2XLloLbwrJRyfrHUpAvtDa/bOXMj5nAXHmU6UFudrstMxWy3NSGSfJomVa97KmHM6nqzDfyGdh8fdqCMDUPWDqZU45yUOHn5o7qNQ1mce2k84mq9SgjkUhGOg0mj4OIRKJ88IMfnaRy/uKLz/PQQ99HkiRuvPEWbrnl3RXZyJkzXVx11TY+//n7cDodJRlJmJuyntluytlYa7u9cvNuKglBINOtZHrvMfz+sTlhULOnHk50oUx4OQvtaU+HmVnz0nOw8820p1Mj1ZyfhCr2KOPx+CRDmT0O4sSJHj70ofdlDKWqqnzjG//Agw/+O3a7nXvv/Sjbt19NfX1D2Ru56aZ3AdU1DiIdpua7ufKx1iMj+UfPVrrUqDwYCIKIx+PCZrOQSKiEw9GyZ/oUg3QXStrLyZ7aOGF8SmGFS1WOKsxDyxeap0vUiu81L6f+8vwU7YUqNpT19fWTHphpHERX1ymWLl2Gx+MBYPPmC3n99b1ce+3b5m/Hs2Bu5uZM3JSlsNaV2lM2e14s0obd6bQhCALhcHTBc6b5CKC0V17KyzR/1ObHgKRDc5/PhaZp04bmM++1tGufr33eUMW93rlfvOxxEPff/78njYMIh8OTvE+Hw0k4XL4MWSVR6dAbhMx0xdkLwucaxV/QzEvZMjWnkUgMRZEJh6srTzVBAJk5QZBxu+eeAKoEZgrNp2PN5ytHmYsqfhuBKs5R6rpOLqGTHgdx883v4brrbsg87nQ6iUQmRslGIpMNZ6VQDeMgRNGUfKurM73nmULr+dpToSF8uqzKbremwsUJYsZqtaAo81flUOosGrvdSjgczRBAxcqQFbfHyjDXpbDmpaCcCYzVjqoNvXMfSI+D+MIXvkB7+8ZJz7W1raSn5wzj42PY7Q5ef30vd911d4W3VBlhjFKRXRpjGKYXHY3mF8EtBpX5/GcmhdIam7JszuSeXkau/M3Mx/1cLAE03yHpTNebjTUvd17O+QqpSl+bnPthZY+DSIcN2eMgPvGJT/OZz3wSXde58cZbaGqaP5XhQlCKocwfWgfx+TwVaccz3+PKlGHlvjRZNtsH0x1Ls0+GrET/dNlLFI1sAsgcQFYpAmhuB5JNx5obBmVUF5yfEmtQeY/yjTfe4Gtf+xo/+MEPJj3+0EMP8ZOf/IQ0R/OlL32JVatWTbuOnNthcf/993P//fcD+UdBbN9+Ndu3X13u/mdA+aIPhYy8zBbyNev84lNC63LIk9w9VUq4FybEQNJ7L1TdvJJ7mWtM976bivQJotHEFCm1+e4xL9bAZofmTqcNURQz3VrFhOY1MqcwPPjggzz66KPY7fYpzx04cICvfvWrbNy4Mc+ZefZVsV1VCOWOg9B1A2kG/z1XUGMmId/KFa9XLm9qt9twu51FiRCfr5gspaZgsykIgulxFqfbWGq5TellOmlPM5lUCxLqnXTV81S0FyrrUS5fvpxvfOMbfO5zn5vy3MGDB/ne977H0NAQ11xzDX/8x38841p5yZxzGaZxm/x6skNrVS2m66RyJEypyCZmzC9SgtHR8fP2i5JGMa8vLQ5hGDqKoiBJUka30SSAZvbU5ksqLedswJg2NJ+JNTcNZanXrW5UkvW+/vrr6enpyfvcjTfeyPve9z5cLhef+MQnePbZZ3nrW9867VpV51GWi3SOMt+MnGJZ63R5UKX2VAyyx86m5/uk2wvPdyNZOgQMI00AkfIyrTgcQqbMqFokyPIZ2amsuTIDa14dr6PSmA/W2zAMPvjBD+J2m1Mjd+zYwaFDh84tQ1muB5dWgm5srCs7PK1UjtLE7OtMJZVijI5OEDOVKDOqZJ1p2tsFYYaynYVLiE70OYspAii/p1a6VNrc1ULmsubZofn5HHoL8yDYEQqFuOmmm3jiiSdwOBy8/PLL3HrrrTOeM4X1PheRNjA2mzU1G9zA7w8suIJQIetkEzOler7zDUUxS3WsVivxeAJFMb/I049ZmJ97bLpQOHtE7GRPLVmypsB8IV9oLooCDoetJBm4av+6C8Lc3fePPfYYkUiEO++8k09/+tN84AMfwGKxcOWVV7Jjx44ZzxWMGaxJPtZ7PqBps9ct5guto9EYoijidjvx+8fK3ofTaUcQBEKhSFnrmIK2jkl7SrO1VqupEh6Nxqb1fCPhEMf2P47XkSCSdNC27nocztIK/RVFLun9Mcc+WLHZbBiGjqbpjI6OZwn3ml6bxSKjaXoWUaEgy2ImHC4UPp+LQKC4ri+r1dxHIXWvadLHYlEQBDPsDYdjRf24povKw+FYUfsEcLsdRCKxkn4QfT4X4XAs9d5KRbHmmjZ3xrKpyV32Gn2Rwse4LHH4yr5eoai60NvE9GzibKx1uax5NgotNZp9HfO1TExWtGEYBtFojPHx8KxfzmP7H+e6yyTADhjseuVxLrzizrL3VQjMHJ+pLhSNxvD7A1gsFhQlV7h3steWLg7XNB3DKM4YzEf5iykEbfaZm69PTHUAmcXsleqimQ7lvsaZQvNzeRzEfITepaBKDeVkZIfWs2klVl5mrfx1FEVGliUaG31FCxEDuKwxwJnzd2kopI4yN1eaVheKREIcP/o8HrdBNKawYuVVuFyeKeenv8TpciZFUXC5xCJCxVILx0vL3RmGkSF70p5+IRMbK8F6l4tiWfM3c+hdDmRN06bI+lcD8uXu/P7Zc3eVJWBKN7q5xt0wDIaGRktaKxS3zfh3ccjP5E/WphTz5kqPH32e697Wkvl719PPc9HWm6a9kq4bJJMqhqGTTGoZdaBq06DMNniFEkCpMynV2M2F4O/0rLm597n2kisBkerc45TOnIVEMpnkd797iV//+peEwyG+//3/UzRrXUmF82LLgwRByPRbp42N32/mXOrqvCXvY82mG9n1yuN47XEiqpM1m24sea1cyLI5ctZmMzUVZ9KmdLq0Gf+eCdle5mQNysSczaUpHFMNXr5Ugt1uyxhSTdPLlEqb26Fk+UNz8PvLy7fPNarWo/ziF7/In//5n2OzTfVSDh48wHe+8y9885vfm/T4ww//iMce+wU+n5lM/dzn7mP58rayN/PAA39Lb28P73jHO9mx4xrGxoqXcJsL4d7ZYLK+NqxWC4lEctL4XDBzk+VsyeF0ceEVd6bCWInx8fDsJ80Cu93MPYqiSCRSWAtkOCTN+Hc+5L7ubA1KswXRMsUAzTdm+2yyjY7VquB02jEMA1030PX59YBy9VFnQ3ZoPp9BgWF+AAAgAElEQVSD6UqFKCz0j2Z+yC+88AKf+cxnphjKBx98kJ/97OfYbFP7JDs7j3D//V9i3br1Fd3Mffd9EQBdTxZNAMwFZjK6xREzlSlcL3cdRZFxOu0pkVxzLnQx3np7x9XsenoiR9neUVjP//QKOxrJZHSKASq17rUc+bFCvLRsAshMrVhRFAuiKM4LATSxj9LOOxcIHqFaQ+9//ud/ziiWZ2P58uX8zd88wJe//IUpz3V2HuaHP/y/jIyMcNVV27n77g9XeFuVKaquTB3l5L3kKpyPjQULUDivnBBFsevkCmjEYnFkWSMQKL70y+Fwsfmid2YmK46PV0a0OdcA2WwWBEHI1GbOvZdZfK7RFBpWSSQAjKIVzkvF+VxsDiBWa+i9adOmvE9cf/317NvXmfe5nTuv4z3vuQOn08l9932Wl156gW3b3jKX+ywKhmEgigKaVq6hNG/M7BpCTdOIRGJFKZxXUri3UGSnA+LxZCbXK0mmIapWqKr5/rpcDjTNmORlJhLJOWFty5NZ06dlnaebNLlQ6ubnAs6b8iDDMLjjjvdllM2vvHI7x451VpWh1PXyDVM6HJRlifp6b6aGsFTvpjL5oZlD77T36HCYaZRIJM74+GjOF6tSaYC5RzyeIB6fmKeT1tw0vcx8M7YXVrg3l3XOFhrOzr8u1Lycc8G+ViuZIz700EN0d3cXfEI4HOYDH7iTSCSCYRjs2bObtWvXVXxj5c2TLt1QKoqMx+OiqakOi0VB13WGhvyEQpGyQsBKeAHThfAWi4LXa+5ZUWTGxkIMDweIRKLnqPcxORRWVY1w2MwBa5qG02nF43FgtSoVSWnMhXBvMqkSCkUJBs3victlx+0291ybwDg9BLSC/80n5AMHDrB161aWL18+44G7dj1JNBrhD/7gPdxzz8f41Kf+BEVRuOSSy7jyyu0V3tb8Cj+kC6NNKTPTExseNmseGxrqytrLxJ4q6/FM7Lm4Lh84N4R78yGbwU2Ln6Tl1Mye7fnVlSzk85xu0qRhGMiyVDQBdD6L9gJIVcp6C6qqGtPVUi5Ur7eua+h66eUMPp+HaHS2kQiT2yFjsQTRaGxSTkkQBJqa6jh71l/SPrLR3FzP0FBuGFwc0j3juq5Pu+fZIIoiDQ3eGYvfI+Egx3b/FJcYJKS7WXPpbZN6y7PJnHSv93Qopv86DUkScThsBIOz1/wJgvm+mN6liKqqRCLxot7nUvuunU4biYRa9MAzs3bVmsmBFzNpMp1fjsWKn8Ou6+a/uUIler2j6hsFH2uXLyz7eoVCAv5SUZS8HmUkUvyHUQkYBhiGVrKhtFotaJqW99dakkScTjteryszhGt8PEw8nshbWuJyOSoy2tXhsBONFvcFTkMUTcPhdjuRJDGVewxNu+eZYLLJNiKR6dsgD/3Pf3Ld2hDLGgRWNybYve84LW0XZZ5XFBlJkojHExmyazrIsoQgCEV5TqIooChywWVCaeJEUaRMnlaWxUyt42yw2YpVRDdhtSqoqlY0wy2KIpIkEQpFUVUNRTGL/hVFwjCY8TNVFBmjxHk75veq6NMKhtNpLXsNzehHEIyC/sni4grsujDIhw4dQlEUtm+vdPhcOsoVtsgNvSe350nEYjH8/vEyhjuVtKuiQ950N4jpPcYJhcLY7Tai0dJ7vQshc1xicNIx5t/VD5MdV0kmo5me7cI8tvnN+2Uz18VOmjRx/sbeVVtH+Yd/+Ifs2bNnofdRUaQNpaLIqX5rC4mEOmN73vRrVaYkwzy9kKFnpvc4eSZ3EMOg4vO4dz3+I6678f1THg/pbiCU8/dUpOsz09MQF753e8LgpcuJJEnMmg2eXxloLvquZz4v//WyJ01ObvOcmDR5vpcHVW0dpdVqJRYrx0OpLgiCkCpatmZIjkInFOaHUZEE+mw96FarJRV+TYx+mItOjzSpdOzIHt6yIcTxzr20r90CpHKT//Mwdt3Pf+4apnlxM5qlkTWX3pY536zDtKYEZEUzRxqPc/CJ72ENDxC2NrH2+o/idJWfryoW041XmPDYJisDFZoXnIOdMpNXOB0BlO6VP4/tZPXWUS5atIjBwcGF3scklBJ2WywKDocZpmqaRjyeqEjnyIQwRrl359SQN3twmPmFjjE6Oj7jXirVr6uO7adtq499ux5H6HuOkOElFolwc0cqzF5Zz1OnfVx4zYcAk0RwOOxIkkQymSQWizM2ZpI5b/z827zVdhRcAH5e+s1DXP7ez6a89+L3W2w/c6HI9jJNxtyZSr/MHetd7nlmB1B60qSMxWIKrsTj4gIa+rmDKFSn4rz405/+NG8L47kASRJxuRw0NdXhcjmIx5MMDY2WTJrkQyW7atLL2GwW6uo8NDT4EAQBv38cv3+MWKw4FfDZ8J1//N95Hz+472Uuajf//+rLluCTRrhudYCE/9ik41zCWOr9rcdutxEORxka8qcmHk4c54gPTTpPHu8HyNQNliJ+PJff/3Rh+NhYKDWDBjweJzabpci9zp9uZro0ytRjTSCKAh6PE6fTXnBK5lywqYKgF/xvPiEnEgnuvffeeb1ouUhLmeWGqfF4nNdf34soCjQ1NREOm+UlnZ2dGIbB2rVrEUWToe3pOU0kEqWjY22GtR0YGCQQ8LN27brMcYcOxTl9upv29o6MwQwERhkZGZ702NhYgNHRAG1tKzKPjY4GSCbjNDe3EAi40XUIBs1Rs3a7I5N3Ghoaoq6uDlmWSXtggUCAujqzhjO9XjKZpKWlEb/fHMFgMsoiVqsVSZKQZRlZlhEEgd8+9V988JYVPPf0T9jxttvZ9fN/47p3fwQwOHn0NYIulXjwOG4r9A2G2LCiHq9rsjCKZluEIDBrR1LE2gSMTvo73aHidNqQZQmv1znn4W4pHl4ioeJwQDAYnaQ/OV37YbnXKx8Cmqal7p1iCKBzA0IxYjjzWA8s6LpupH/hcj2nhaqjBNC0BNm/1rIsp/qtzVxNJDJ9naQkCQwODjI4OER9fX2K+UzQ2XkEgPXrNyDLMqqqcuTIIXTdYN269VgsFpLJJJ2dhzOPNTc3MDIyyuHDhzCMiXMjkTAnT54ARDZs2IAoioyOjtLTcwabzcKaNWa30ujoMGNjY9TX17F06TIikQgnThwnFIrQ1NTEokUt6LrOkSOH0TSdlStX4nK5CIWCnDp1CkGYWH9goJ9kMoauCyxbthxN0zhy5DCxWIzFixfjcrmIRiMcO3aMUwce4+53b+GF3d3otk1o3bswlryVi6/YQWNjHXv37uPgy49z1xUOFFlC03Ue/n0Yn9tDaz1EDBfu5Zcz4h9l/fr1SJLptZw4cYLGxnqWL19OKBRBFAUEYPyNX+FMDBOxNrH+nf8rk6NM1/0lEmpmRo35o5aYNgebFsYIhYoryyq1HjJ3Pk+6/VCSxJTxSeQtAfJ6XYyPh4o2lna7OQSvlMFmHo+TUCgyZT+iODEDKJcASmMu5+VAZeoo0X9T+LHizvKvVyCE73znO0YoFOKzn/3slCcX2lAKApOUb6LRGNFofNZfTDP5bcfvD/Daa6/S29uD2+1CkiQMw+DUqS6SyQQrVqxMTW6E7u4uEokkra2tOJ1OdF3n1KlTWK0KjY3NWK0WDMPg5MnjgMjixUux2+3ousrJkycRRYnW1taMse3r68HhsLN8+QpUVSORiNPV1Y0sSyxbthwwGB4eYWRkBK/XR0vLIgDOnDlDLBaluXkRXq8XXdc4efIkAG1tbdTX++jrG2RkZASXy8nateuQZYWBgT7i8QTNzc10HX+N7WtDdPePEgzFOHxkhHuuXcJj+5K423eydOli6uqaCIeCPPfL/4NViNDY0saWa+5ixO+nt7cPm81Ce/vajBHXdZ3W1mX4fD5GR4cZGfEDAh0dHQiCyEB/Lz2/+RnLhQhqw1IueO/HcbrceQukTS/IMm3pzkIbyjSyjU8+lfNSBqABOBxWVFUvSU7O63UyPh6Z0Ss3dT4VZFlOTZpUU3XFRV+uKFTEUGpPFX6s9Pbyr1cghLVr1xoAR44cmfLkQhhKXdd59dVX+NWvHgUMHnjga0Sj8aJuqnxTD3t7ewkGx3E4nFgs5pf3+PGjRKMxWltbqaurB+DEieNEo1GWLFlMfX0jHo+LAwcOMjY2xuLFi2lsbMqcG4vFM8cBnDp1HFEUWbmyDbvdSSQS4cCBg+i6xubNm7FYLIyNBTl+/FjKO12PLMsZT9RqVejoWJ/Z7+ioaURNw2pe02KRsdvduFxOxseDHDt2lHg8zooVKwCDoaERfvXf32f1cjeLGp0gCOx+7QTbOhpw26yM1b+FtjXrOHnyFBaLQnt7O4lEktOnuwkGg/h8dTQ3m6HzwEAvqqrR1rYSWZbo6+slEAgQPbGXxUoSv+5AXrUZRbFinHidzaHjBOIq7fVOzq6+ims+8cUZO0nSpTuKokwq3akWQ5mNqV5mEo/HUaKhtKUk2oq3XMUYZ5MAUlL95Qajo5Hq9yiTTxZ+rHJD+dcrEMJ3v/tdY926dVx99VQR1qGh4LQq5y+++DwPPfR9JEnixhtv4ZZb3l2RDX3hC39Ob28PN954Ezt37sRVQpmJLMt4vU5GRqaOZD17dpC9e/ciCGC1mp0Evb19xGJR3G4vTU2NgMHAwCDhcIhFi5qpr68nHlcZHBwgHA7hdLpYtKgFQRAYHh4gHo/T2NhIfX0j8Xicrq5TxONxfL46GhoaMAyDwcEBYrEoPl8jXq8XwzA9RUGApUtbsdnsRKMxuru7kSSRtrY2JEkkEAgwNDSMy+Wkvb2dhoZ6OjuPEY1GaGpqZvnyFSiKQmfnETRNY/nyFXi9Xva98jzxV/6DN+I+dl68FFEU6R8ao3tM4ab3/jGLFrUQDIY5cOAAAKtWrcbhcDA4OMjZs4PY7Xba29eg63pqbZVly1Zw5umHWTW0l8FIApssMtJ6GRfe/if88sufZI3up73Oidsi85JaR8eH/5zFi1tobm7OGMrQ+DivfuMB5J5u1NblXPKpz+Fyu7PaEAVUVUUQxKI7ojweB6FQrKgcnSAIeDwOxsYKU41P1zim9xqJxIo2eKW2PkLpXqwoiiQSc5u7rIihjD9R+LHWd5Z/vQIhxONxI+1h5eKf/umb/PrXT2Cz2fne9x7KPK6qKu9//208+OC/Y7fbuffej/L3f/+P1Nc3lL0hVVWRZbkslXNJkqirczM8PP2M4JGREfz+EURRxJnqY+7v72d4eAifb8KLi0SCnDnTg9Vqpa1tFWCSL4HAMPX19XR0dBCNmoTP4OAgLpczc1xvby9+/zA+Xx0bNqwjmVTp7DzK+PgY9fUNLF3aCsDx48eIRqMZj1XX9VROVKetbRUOh4OxsTGOHu3E4bDS2NiCrmsMDZ2lr68PRVFobV0GQF9fH6FQiOjR59nekuCNcR/Kmu2IIqxcuRJBkEgkYhw7dgJFsdDW1pbaaw/hcISmpkbq6uoIh8MMDPSjaTqrV7cjSRK9vT10Pf0wF9rjNLstRBM6L43KNG+/Bf3k61wUOkFXIEJC0xlsXMPlt38Ew9DQdZ316y+gtXUZz/3tF9j00tOZz2H/trdxzf/75azPTky1IEoZz63QDqrp8nczoVhDmQ2fz0UyqaX2WvgoC5fLTiw2fY52tmuWYijNtseiTysKFTGUsV8Wfqxt+sF2lYZssVhIJBLkM5ZLl7bmVTnv6jrF0qXLMmVFmzdfyOuv7+Xaa99W/obk8rtPCinpaWhooKGhgfHxcbq7T3Pq1CmcTieybBqEvXv3YLfb6ehoR5ZlBgcH2b9/P16vhzVr2rFarXR3d/P88y9hs9lpa1uOJMmcPTvEG2/sw+FwsGxZK4qi0Nvbw5EjB7FYLDQ1LUIQRE6ePMnu3S/j8XhpaWnJeIV79uyhqamBJUuWIcsyJ04cRVU1li5dymWXXUZDg48XX/wduq7T0bGOiy7amiJ+uhAEgbe+9Vq6Du/BHdGJJWGTO4i1YxlL2zdz6NAhEokomzZtYODAa1y28zZOnjyJYcCOHW9FlmV6es4QCARobl7ExRdfmiK8DqOqKldeuY1FgZOsGnqdwXACr13m8o7L2HjzLby+dwmvP/84K3wqSnMr6y6/gXAkjM/no729nUgkxqFDBxk6fIikbqCkynDknskSf2nG2TCMlKSaDZhQDao0yu3KCYenjrLIR6RUZq/nd1cOAFUwAiYf5NOnT6fyW1NxzTU76e/vm/J4OBzOCPcCOBxOwuHKjAWoBIqpffR4PGzcuImVK1cxMDBAMplg/foLAAiFgpw9O4CiyFx33duRJAm/f5SDBw9iGAIXXngxW7ZcSjQayeQdt269GFmWiUYjnDhxAsMwuPTSy/F6XQSDQQ4ePAwIbNu2HVEUCQRGOXPmDEDKMCUZGBjg6NEjWCwKS5a0omkqJ092sXv3blpbl2CzuZAkkd/97kXC4ShNTU34fD50Xefpp3dx+PdPsb5BRpYshCMJjj30fdZeupNVq9qQZYlHf/KfRF94lEfDSS66fBu6rrFr169R1SRLl7bicNjp7e1h3743EASBlStXYhjw4osvEKKefeP1LJLjhAQPdudSen/9BCtWrGDdDbdz5kwPqqri7e2lsbEJv3+En/3sDTRNZ9Wq1ZyxOvj1SIDFFoUOpw21dXp5v2xJtQnh3uS04yHms1wn+1q5oywmOmny73WuJzCe06hWQ/nAAw+gqir3338/ra2tBZ3kdDqJRCZClUhksuGsDConilEInE4nq1evJh6P09PTw/DwWcDA4bATi8X5j//4DwxDYMWKFciyhK7rPPHELzEMk41WFJMZffLJxxEEkVWrVuF2uxAEiRdffB6bzUpHx1paW5chiiKvv74XSZJYu3YdmzZtJhaLcfz40UwJUnv7GoLBcU6dOoUkyVxxxRWIoohhJDlw4DB2u42dO68DoL+/j2g0wqJFTVxxxeXEbryJgwcPkkgkWLJkMTencqfHjx9F11U64v1csqWVXRHTqNlsFt75TnME7okTx4hEIqxcuYorrriKWCzK8ePHAdi581qsViu9vW/B7x/B7XazatUqZFni2LGj/PanD3Pnx/4Um83GwMAgw8NDNDQ0sHHDBbz83X/m8FO/QPPWo11xNdFggBcbmthww7s4e/Yszc3N0342k7tTKu25VbYLKHuvuQPTyh9lcX6L9gJzqwNXBuTe3l66urrw+/0FG8q2tpX09JxhfHwMu93B66/v5a677p7jrc49BEGgrs7LkiWL0HWNw4ePEo2GWbmygba21ei6zqFDB1FVldWr13DBBea8oc7OIyQSSVasWMEFF2xMHXeAs2fPsmhRC5s2bcZiUejs7CQQCNDY2ITX60XTNJ5//rfEYlFaWpbgdnswDJ2nntpFMhln6dLlOBx2NE3jiSceR1WTXHzxFmRZJhgcZ9euJ5BlmY6ODgAOHjzMb3/7HIpiYfnyFciyzP79+wkGg3g8XpqbG+k7fgh97266RQO3bYiBk0doWLaaX/7yUTRNp729HUmSOXLkCLt378Zut7FsmZn/fOaZ36CqKs3Ni3C73fT3D3DgwD503UCJRVm2/3f87Id1tKxajcfjYdGiRfj9fh7+uy/T2nOCVS4rYmCQ/2lcgbbznXi9HsbGAwy+ZhJimzZdyKpVq6btiTeN4+TxEHa7NVNiVNpnPjeCGPm8zHQN8EJ4lOeMJzrP438LhfyP//iPxGKxgoxktsr5Jz7xaT7zmU+i6zo33ngLTU3TewQLgWImMabFEsxBXGaPeDKp0tKyGLvdRm/vGfx+P263m40bTeN49Ggn8XiM1tZlmVEYp06d4PTpLlpaWti4cTMAp0+fpLe3hxUrlnP11VcTDkc5fbqL8fFxmpqauO66GzLHhUIRmpub2LDBDP3NUqUILS2LM9cdHOxFUeJs3HgxTqeLUCjMvn37MQyDTZs2Y7c7CIWCdHV1AQJvecsORFFkaOgsAwP97Hv2cW7bugktqdI7Os5A5wEuuHQbW7ZsTTHcZs3k5ZdfhsvlYnR0lDNnelAUmeuvvwFRFDl58hTB4DjLli3jyiuvIJFI8osvfIbFCngCfbztPX9Gd3c3IyMjNDY28tYmD+2qj85wFAOdi+wyb7npZgYG+ug+3U3vU4+yPBHjjVdewP+u97J4SQttbStn/MxyPTe3255SjJJKKruZS0zt15ZxuewZj7hQI1bLUS4chEceecR4+9vfnrffeyELzstVOW9qqmNkZGzaUpF8oxTy9Yg7HKaGZTAYZvfuV+jv70dRLNhsZmnRwMAAY2PjNDQ00NjYAAj4/cOMjIykahKbMQyDYDDI6OgIsmxl8eKWVI+3n5GREZxOF0uXmiKkw8PD+P2juN0uWlpaEASztTIWi9HQUMeyZcuIx+OcOHGSYDBEfb15XcMw6O4+QyIRo7m5Ba/Xg6bpdHWdQtM0VqxoY+j0CWwv/4zejh046ltob1+NIAgMDY0wPh7AbrezbJn5g9ndfQZNU6mrq6Ourp5kMkFX12kMw0iF2zLDw0OMjY0x2tfDmpefxCEJHArECG2/gY2XXYHX6yUWi/LSD/+Npac7We6wIAoiv4oIrL3zAzQvWkz/M0+w+Mge+mIJFEFA3Xw5W//ww0QiYaxWO5dccmnBWgQ+nwtV1WbtqMlGqTWbpZ4HZkdPOGxqZiqKTDJpDh+bjQVXFNPIhsPFq33Ntbo5VIb1NgL/WfCxgu+usq9XKOTnnnuO48ePc8stt7B+/fp5u/BcI60jmYtslaFC5nJn5zsvvfQyAIaHhxgdHUWSJNav3wCYpUAjI0P4fHWsX29OpBwaGmJwcACn08HWrVuQJJHu7h56e3uw221s22aKJZteWzd2u43t2816Vr9/hP7+frxeN7fdto1kUk0VhY+hqgbXXvv2zHX9/mHcbg833PAOAE6f7iIYHKe+viHjif7wX/4eb+QsVy9x02ML03bjTRw+fBjQuPzyS3A6XYyOBujuPo0sy1x33XWIokhPTy+jo358vjpuueXClNd5FNC49NJLaWlZxI//5ouciBnY1CjL6n3ELdDSspjx8TEaG5v48Je+ytOf+xgnentIaDrvdTo523eaFdu2M/hfZ2m0KFxS50bVdZ4cOYuqqmzcuAlZtjA4OEh/fx91dfUz5jHTCIWimY6aYvq25xOCkO1lpu9JKyAQjyemVVs/V2cdFYUqzVEKqqoaX/nKV2hubuaee+6Z9OTCepQ6up4o2aNsaPAxNhZEVTVEUcy0QhqGQSQSIxYrTGHIZjND8rGxqe9FIBDg97//HZqm4nA4Mo+NjPix2SwsXdqKYUA0GsHv9yPLEosXLwVMD7avrw/D0FN5OZFYLEZvbw+SZPZ3WywWAoEAp06dQlU12tvbUy2MA/T29iOKBitXrkqx52MMDw8jSQIrV5p1nMPDw4yOBgj0n6Hj5IvsGRhly9I6hiIJjMv/gI2XXYmqanR1dWEYOg0NjXi9PlQ1walTXYiiyLJly5FlmbGxAIODgzgcdtavT6caTvHiY4+xoqOD5c8+ih04FI6T2HkTl1x9DVarjdFRP8PDQ5x55MdcGjKFM87Ekgw0LWLbp/+ME4//lEUH9zAQTSIJkLzwMkI2N+1XXkEwGKGxsRGfr454PEY0GmXlylVs2bIl732Rr8Ywu6MmHepme5mKImG1Fu8Zmm2ClpLGhExXC2l2KllQFDmvyLDZliim9DWLwznjUY78oOBjhYb540VkSZJYuXJlpkSlWlCJcRDpVsb0KIVAIIhaZMPrTIK7Pp+PG254B6FQkMHBQVRVzZQWjY+Pcfq06Z1t2XJxykhr7NnzOoYBmzZtZsuWrRlhDtC57LJL8XiuJR6P89prewgGQ7S3d7B+/UZUVeXw4YMkEkk2btzI1q2XkEgkOHr0MKqqsXXrxVgslkllSVdccRWyLLP7239Fg01micOK6mxg48oGOkNDyLKEz+dly5YtxOMJuru7CQRGaWhoZNMmM8d65IjZ8bNly0UsXtzC+HiQ/fsPoGkajbLCu3tP8KKW5KyvCYsis2mVixHJwOl0Eo8n6OjoYNu2bTxz9DDaS88CsMPn5PjFF7NkyWIs77qTM3YHHeFxoo2LOOttgod+QMPVO9hx09s4e7afoaFhfD4ffbt3c8UVV/L1z36O93z83kxv/UxIJtWU4K2Y42UmSCZNLcrS8n6VZ6DTmqQwkTfP7od/c5QHVecLlAH8fv85Lc2UjfQoBUWREUWBcDiaGaVQCgopNXK53LhcbmKxGKdPn6az8wgejxdFsZBMJnj88cdQFJm1a9elhDl0fvWrxxFFgba2NrxeD5Ik8uSTT5JIJGlubk4Vv1v47W+fRdc1mpqacbvdKIrMs88+i2Fo1Nc34fF4MAyDp57ahShCU1MLLpcLTVN54olfMtzTxUW9h5GsCm1eO092DbPxwhtoa2sjGk0wNtbHa6+9hsPhoLV1GVarja6uU7z66it4vV5WrVqFJIm8+upu/P4AXq8307656/vfpcU/hq0uzoZPfp5wOEh3by+yrJBIJBBFkQMHDjA2FsC2cQvHT3UhHTnCM+4G1i9t4+VvfIfr3vdeGq5+GwN+P06nk/CP/4smw+DX3/lXTgwNcfx3u3nnB+/mob/+W1ydh3mjs5O6Z/fwmyUtHN63jx3veTdXXbWNpqamGT8jXdcz0m8Wi4zNZsXhEFLGsnjM9fiIfKMsDENHVc+P7+m00KsnRZINqaOj4y9feeUV7rnnHurr6yc9uVBTGNMwjMJvYqvVgtvtxO12oKpapmuilLGe2TA9EUtB41ZlWaaxsZHVq9tTDKeFFSvaWLt2HevWrWdszE8sFqO9fTWXXHIpl1yylWAwSDgcwedrYMOGC1izZg2xWIRYLMaiRYvYvPlC1qzpQFXjRKNRli9fxvr1G1i9ugNVTRCLmX3mW7ZsZc2aDmKxKLFYjJaWRVx88SVEjuzB7/cz5OW/Cz4AACAASURBVA/glEWaZZ32q6/DsFjRNI1Vq1axceMmli1bzsiIH4tF4tJLL+GSSy6hpWUxfX39RCJR1nWsY3TXLgLPPc+JI4dJCgLXHHiV9Q4r6ugIg00trLvwIi655FJWrVrF8PBQionfyKZNm6mrr+PQE8/QcGacFesuYPVllzLyzX9HWL+SrVdeyUUXbWHkxElCjzyCiMDmWAJtzRq0H/+aY2oE+0svc7PHS2h4FMt4gq6+HiK/P8BF73wbNp+XoaEhFEVClmf2MIGMEpCqalgscmqqZOFTG8FskxXF4g1tule80HImwzBIJk1ZOkWRUqM4TOWlYsQ/5noCI1RmCiPB1yY2O8s/wbW1/OsVCGH37t2G2T1y6ZQnFzJHCaBpMxunfKMU0irhXq+LRCJZ1DzpfJhJYGM2GIbBmTPdDA6eJRIJ4XSaw6K6u08Ti8Xw+RpwOp0A9PX1EIlEqa9vyPxg9ff3EwqlFX1MhaJoNEJfXz9utzfz2NDQMKOjo6l2SFNs13zMJGHOvPwsjcGzjISjDI+MYlm5lre+98MAnDnTQywWprm5hZaWZmRZ5ujRY+i6xpIlS7HZrESjMR7/xjdpPHiYxRYLiiDwO68dlyKjCALNDgfhJa2sfvv1+P2jeDweWlpMFv/06S4SiTixYT/6936MqBl0E2fY62bLaJyzqxaz5f/5I/x+P0ef/g2Og4cRowkGZRiVLFzkNzjlEtioqEQMUHUDIyLQqcapQ2ZsaSPv/7sv0dDQgKKIjIwE8Ho9XHbZFQWlbtJ5P1XVM0IXE2U701uWUuaVg/nD63LZGB+ffWZ5LtI6lqaBT8u/aQWRVedMjrL/e7MflIKw+J7ZD6oQ5AsvvBDDMBgdHc0oalc7cgdx5Rs9q+uVmi8z+3jX6SAIAitWtLFu3TpGRs4SCoXQdYEVK8yWwNOnTzI+bioUpdnzvr6+jLxa+rGhobMp9tzFVVddydhYiKGhIXp7e7HbHezYsQMQCARG6enpxWKxcM011wCm+EdDQyNWq5X28VHkf/sHht7yNlRVw+VycMstN2G1Wujp6aGvrx9RlNm50+zZT2tjrlzZxlUOJ0vdPk4l4iQNnYubWrnxu99GVVUOHTqMIMDq1aYC0ehogDNnupFlmRtuMGsv/7+/+EsSK5bgUGRa4km2njqLXxDxnurHFk3yrne9m5cdHp5+bh/LsaGg0i4K2BAJx8Y5aygsVaxEJIG91gRLVSseJBp6g/R3Hie5KkFTUyMrVrSlCv4PYbfbZ81jpkPoyaGuZcapjakzS7wnSvfs0mF7epTFRBrBZM3zkVVpVGnqbyoqXEf5xhtv8LWvfY0f/GAySfTMM8/wrW99C1mWufXWW7njjjtmXEf+i7/4C0ZGRvB4PHz961+v6CYrCUkyc482m6nllx4BMR0qN+umtHUURcbhsGG1WojFEtTXN7J0aStHjhwlEAigKAorVkyoDB04sJ+6Oh9Lly5jyZIljIwMs3+/Ka6xenU7TU3NBAKj7N37OvF4ko6OtTQ2NhAKhThw4ACSJLJ+/Qbq6+sJhyPs378fQRDYsGEDDQ0NhEJhvvGXn+M6RUN47QUu/9LXiMdjHDhwAE3T6ejoYN06H6Ojo+zfvx9FUVi3zmS3z54doMtqwUjEWWkxw6tXmxrYt28f9fUNbN5sliAdOXIkI368efNmkskkx44d5bn/+m8Wt7Zy61f/lv7+fr77mf+NX4/SIlqwCiI//od/wb5kES88/FO24sKCwDF0XrWotMo2rk34OBGJ0COoxJvrWCUIbAyrvMwYXmSOPP0sn//3B4nHoxw4cARRFNiw4QJEUeTEiRPIssSiRS0F1WOmI5PJZTtMKQ43Dd78zMuZOHeqwZsYcZtLVlVXSVTBqKBFf/DBB3n00Uex2+2THk8mk3zlK1/hkUcewW63c9ddd3HttdfS2Ng47VrSX//1X/9lc3Mzu3bt4u67J9PtC52jjMWiRCIhFi9ehNPpIJlUCQZDRCKxgopz04xhORAEAYfDTiQyexmIeawNj8eFzWZ2+YyNhYjHzTInc3iVTENDIyMjI+zb9wYDAwMYhoEoSgwPD/Haa69x+nQXsiwjiiKhUIjdu3dz7NhRJEnGZrMgSQKvvPJ7jhzpRBRFFMWU///973/H4cOHADNlYBgGu3e/wpEjR3j1icdY23ec0aTKwJluftc/jGi1Iwgiqqry8ssvc+jQwczsnXg8ziuvvExn51EsFiuN69ZzYKif1yIROltbeMv/+gh2h4O+vn5ee+1Vuru7cbvdiKJId/dpDh8+zPj4GAPHTnDyX3/I4GtvsHf0LFa3k5GuMwg6vBwbY8ihsPJsFGHNCs72D3AyGaHPY8HWVE/9RRuIojE2NEzQ0BEFgfqYyvFwgFGSeJEZXVwPi+oZE1TOnOnD43GjafDyy79PjfnQCIfDnD17lkQijqbpmXRH+n0yNTCn3k9pg6Npekow2YYkiei6jiSJAEVLpaVnHZXSPWS1mmM08nuMRmrEhmnM0+2doiimSu2KvlzRqEiOMvB706ss4J/gvXzGpYaHh7n77rt56qmnuP322zOPHzt2jIMHD3LHHXcgSRInTpxA13XWrFkz7VpyOkT7+c9/Xv6LrBC6u0/z3//9CLt2Pcktt9zCxz72iWnn40wHw6Ck6X9T15l5HjeYRtlut6WMYzLTApm7Tna41t7eTnt7O8FgkMFB01imPbh0aZGiKGzdejFbt16Mqibp7DxCMBhm8+bNXHjhRYiiwJ49e4nF4ilS5oLU6IYjaJrOihXL2bRpEzabhW/ddhMng3GaHFbabDb6h/qxWs2BVK2ty9i0aVPq3E5UVWXlypVcdNFFmf52SZZ435e/jMPhJJFI0tnZCUKSnTuvQRQlzp4d4vTpLmw2K+94xztIJpOcOXOGV/7537DqsM3q5eBLe1l0+23s+ND7GTl5mmV/9T2GVyzC5u/izJPPcu+3/h6Px8P4eJCurlN0vXEA1/Nv4LJ4OZ2M0aXH6LPCEm8jK1uWMRgL4Wlv5U+++lfU1dXR29tHV1cXkiRx0003IwgCR48eJZmM09jYhNvtIRwOc/jwIZxOJ8uWLS/IM8wnypEuGi8W89GGOFESJaRITgd+f4RzQlCjgqH39ddfT09Pz5THQ6EQbvdEPtXpdBIKzax+JgMkEgmSyfyel67rfP3rf8fx48dQFIXPf/4vMiKxAP/0T19j377XMwXXf/d3/1C2ktDDD/8Ij8fLd77zPZYsWVK0kYTCDFzh60xdSBBIlZjYM/N8hodHp2VNs8fVZsPtduN2uwmHw+zfv4/BwUHcbk9q+FmSJ574Jbpu8P+zd+bhUZXn+//MvmXfSUISsgCBhMUVFNxwQbQCQhFXVFqXIu5WxbpVrZZWXKu1VfHnVrCuuHxri1JbrSiKEEBISMhK9n0y+8w5vz/OnJOZkGWSTCAg93XlEuecnPfMzJvnPMv93E9OTg56vQFRFPjkk49xuTyMGzeO6OhIYmOj2bLlW6xWq19wI4ZvPvyA9tkn4XK5sNbXEW2zYtZrqbA5qLb7MEU04vUKqNVqvvvuW+x2O3FxCcTFxaHRaPjmm2+w2+3Ex8cTFxePRuOXWetykJSURGxsLF6vl/feex9RFMjIyCAqKgqXy8W6devxeLx4OzpwbN9DlErN165OnOVtbHxjHbMWLWDL395FdLfT+mMzhVgwbC/j76+8Rmz6GOLjE0hISKDyhyL2a50IXg+xeh02o4mEY6dQeP5cBMHH7LQMVCoVmzZ9gcNhIzk5mZSUZBwOB//3fx/hcrkZNy4HvV7Hzp27sFq/xmw2KYLJlZVVZGdnKaIfoewFWZTDYjGiVquJjo7oV/otnBiskRUEUcllHjY5yoPg+kZERGCzdauf2Wy2IMPZG7Rer9ffYWBT1MUD8d///hu3280LL6xl584dPPvsEzz22BrleHHxbtaseZaYmJiwvZE77lgFHDiJcTAIV46yJ7RajT9XKqnWWK22EMP7/otCFouFE0+c4R9Mth+73UFiYrzi1UmyaS7y8vKYNEkq8hQXF2O12khPT+Pss8/BaNRTXl7O1n9/Qepn/yX17LPJmX0q63/3IFZBTb5BixoVFYKbwvpayrdvo/DkkznrLEmu7fU/Pk702WcTGRnJ2WdLr9XWSqrv0dHRioBHVVU1VmsnqaljmDXrJFwuN7t376Gjw0p6ehqXXnos6x9/irLdpbRE6YjwqonyCiS4RXZ99h+SMjOYuKuBKKJpws2eFDOx5giSO2xceOOF1NXV09TUyHk3XEv2mt8r78vhkIx5TEwMTqeLvXtLEEWR2bNnodfraWiop6lJutelS5ciCALbtxdht9spLCwgPj6etrZ2f/eTluzsbMxmAyUl0jzzpKTkkPvKRRGFtG4w6IiIMPmnH7r7DauHl6Mc2u8eNkYSDoqhzMnJobKykvb2dsxmM9999x3Lly/v93e0MtH8Zz/7Wa+k86KibZx44kwACgoK2bNnt3JMEARqaqpZvfoR2tpaOO+8+Zx//vywvaHhbKpwG0q5BVKjUWO3O2lubh8USb8vj1KGSqVCrVah1RrJzc1FEAQqKyspKytFpVKh0Ugcui++2ITPJ5KYmEhkZCRarYYfftjKli1eEhIk769+4yY0bgdvP/Ekx119JVHj8xnz/dfsd3modbho6vTwbUIkkT/uJH7cOLZv305D+T4sH31KTXo6cVmZFBUVodVKXVtarZba2lq2b9+O0WggP38ikZEWysrK+N//vkKr1ZGRkYlGo2H79u2se2kttvUbmL7qZs6/8Ro2ffgJNX/+G9FoyGn2UvSvTQhpRkQBIrRmYiaO47Tll1NWVsaGDRsYO3YsZrMZm83Kv//9GS6Xi7Fjs1CpRLZt20ZXlxWTyUxaWhqCIPDZZxvxeDyKoauvr2f79m2Aivz8fGJioti1awdWaxdGo5GUlFQ8HjcfffQhKpXon0+kZ/fuH0lLG0taWhrJyckhfaeBcmpyO6Qk/db7jO3hVb0PM6M3FIygofzwww+x2+1cdNFF3HXXXSxfvhxRFFm0aNGA37dWpk4sW7asVxqFzWZTZsqAlIyWPU+n08GiRUtYuvQyBMHHypXXMXGiJDp7qBEOQ6nRSN4jSGG2zeYYUhqgt/uR/1/Ko6oOyKdK0xzHkZWVRU1NDV1dXUREmBVvsry8AqfTRmZmJjNnzsDlkub2fLNxI3F79pJniiS/uQONUc/3O7ahVulI0+pJM6nR+TQ0iipmXnIZJZ//hwU3XMfHt9yBKAjUf/Zv5r7wJ9RqNXa7w19E0nDqqadgNpvo6Ohk164f8fl8HHvscajVaqzWLioqyhFFgVmzZuN471OMGCh+//8YO3kycV1u4qZOo93poN3lID45iV+9+ldA8k47OtqJioriggsuQKfTsW9fKV6vl7S0MRQWTsXtdrNnzx5UKjjhhBP8FKRu+bdzz5WGTFVUVNLZ2UFKSirHH3+Cn7r0IzYbzJw5k4SEBOrr6/15TD0///kiDAY93333Aw6Hk4yMLJKSknA4HD3ymL2nXnrC4/Hh8TgUUnlkpKmXCvTQWx9/EjJrg5h3FArS09N56623AMkZlHHGGWdwxhlnhHwdldjPJ9/UZOWZZ9YwaVIhc+ZIajULF87jvfekSWkS2dWJ2SxVEZ977imys3OZO/e8wb+jXjCcAWM6nZbISEvQyNpQYTTqMZmM/pEOTkwmQ7+SbaFApVKRmBhHY2MLkmGU/thUKnXI16irq8VqteJ02tFoNLhcbiorK2lqavJPfEzgq9dew7tvH+1uNwbUROaP54zlVyOKAv/9x6e4179Njl6PHRUNZ5xK0wefojrlRPK2bMOCii6g+YJzSB8/nqysLCIjI1CrVezatRtRFElNTcVgMCIIAmVle/3zzNPQ6434fF6+/Men2Nd/SIZKj0+jxnTDMoiURgTn5eViNptoa2unrGwfarUk4KFSqbHbu2htbcVsNpGcPAav10djYwNWq5WIiEjGjElBmsm+D0EQSUiIJzo6BrfbQ0VFGSqVlqysLLRarV/mrhW93kBmZqZfgq7CLzqcSGpqKm63m717S3G73WRljUOj0dHYKE3eNBpNpKWlAdDR0UFhYSGZmVnK1E4IfUBY8JhbKUUjir2P7x0Io3mwGISJcL7rdyGfq5q8atjrhQqtzyep6/TlfRUWTuWrr/7LnDlnsXPnDrKzc5Vj1dVV3H//3bz88huIokhR0Xbmzg3nZLSDNw6iu8tH5mk6lc1sNBqGLdDh/xfR0RE4na5BJ/51Oi3jx+eh02kpLS2jpmY/BoOJY489HrVaRUtLC/X1Dcz55S/IyckBoKvLxr59ZbhcLvLzJ1Hm/hhHdDQVHZ14TAaEb39griaGzSXl2NLG0IVIsslClErDRZdfys6dO7HZ7CQnj2HePMlrKykpweVyk5qaSkHB5KDX0tPHkuZTYZxSyH6bFa8gIJaWc8mDv1Gq8W1tHUyYMIGTTjoJr9fLzp07UKvVFBZOBtR0dHSwb18ZarWaU045xS863ER9fQN6fffICtkTjYuLY+rUKX75t2I8Hi+TJhUQFxeLzWZn374yAM4773wiIswKsd5oNLNw4UJ0Oi1FRUXY7U6mTZtKTEwc7e2t7N9fi1qtZcaMmX6x4n3odFoSE5OIjo4O+XvrWYE2GHQIgoDPJxyePMeRxijVnFA8yt4MS1OTVal6l5WVIooiq1bdz9dff0l6+lhmzTqVN998lc8/3+jvwpjHggWLw3ZzguAdVL93IDQaNbGx0TQ3t/V7Xs8uH7vdeUCXT6Bk21Agh9cajRqjUY/RaEQQfDgcrn4LQZJ6tywSLOJ0uoI8EWnmdxOiKChFiI6OTiorK9Dp9Iocmtfr5c8P/Y708Tkkv/gm0ajpEnx8Y7eTgpYojY6sJx9g2umnUVlZTkdHJwkJiUqn1r595dhsXSQmJpKSkgJIhkoeu5uaOkZ5raOjndjYONLTpfxhaWmZX5EonsTERHw+H1VVVbjdTlJSUhRRD2lkRRfJyWMwGo2IokBpaSkqFaSnZ2A0GnE6HVRXV6PVasnMzPQPe2ulqakJi8WieIG1tfvp6uoiPj6BtLRUtFoNxcUl+Hw+EhOTMZlMfsWnenQ6PXl5ORgMRqqrq2hv78BisZCQkITX66WiohyNRkN6+liMRiMOhx2fz8ecOWcQHR0X8ihdGZIikLQ/1Wq1MkBtoJB6OGN1D0b7IoTJo9z+25DPVU29b+CTwgTVo48+Ks6ePZuTTz75gIOHutd7OIZSrVYRHx9LU1NrL8ckfUqz2XhAj3hviIuLxmq1heQByA8c6advuThJRsuIWq3C4XDhcnVTODQaDSaTAb1eh9vtxensn2BvtVr57rsttLe3ERUlsQ9kHuOWf27kmBkn0vjUWtoKMrHvK8cQHY3B2kWXCGmJKbS7XDjHxHPetVdjMJj9veJNtLa2ERMT7RfMlVokm5qaiIiwMGZMKiqV9FpzczMGg5H0dEl/s6OjjcbGRiyWCMWgtbd3YLV2EBcXR3b2OLxegf37a+nokF7Ly8vBYjFTXl5Ba2sb0dExpKdLhk/2WFNSkklMTFRaFAVB0vKU6B7d3mNh4RTMZiMtLc1UV+9Hq5Uq3ABlZWV+RkECKSkpuN0eJf85fvx4YmNjaGtrpaamBrVa659AqWLvXslbTU5OISkpCa1WRUtLG0ajiczMzJAjDrNZilgCu2n0et0A7ZLSfo6IMNPZeYQbyh8eCPlc1fTQzx0uVE8//bSoUqno7Oxk4cKFQSrnh95QDn0chEoFiYnx/pyghEB1c4lfNnCHD0BsbBQ2m6Nf72+g4kxf0Go1mExGRaxVSoOo/d5jaOLCMmThX5fLRWxsLLu+/IaSGx+mIiua5IpGfPmZRBU3kvf0PURljMHr9ZCamkpERBSiCE1NDbS3txEREUliYjIg0tjYSENDIxaLhexsaY5NW1sbVVXVGAx6hSTf3d+tY9IkaQ+1t3dQU1OD2Wxk2rSpeL0+6usbqayUSOETJ05ErVbT0tJCbW0ter2ewsJCDAY9jY0NVFVVo9PpFSNXXV1De7tkRDMyJO5jSUkJDofknUq5Ui27dv2Iy+UkLS2DiAgLVmsX5eXS/PLx4/MQRclINzY2YDJZSE0d42+33IvTKfFEY2Ji8Pl8VFZW4PMJjB07Fo1Gmu/e1WUlOTmZ6OhYfD4f7e1SMerUU08LymP2BovFiNvtPeChq9fr+hXlGI6YxmFlKL+/P+RzVcc+OOz1Ql5LFEWxrq6OiooK8vLygvodD2dDCZCcHE9TU2vQbBy73eE3QKFfJzY2Crvd0ac01kDeY39Qq9UYjQaMRr1ibD0eqZd9sCLDMrxeL9XVVTxwwRKS663EatWY0VAhOOkQfKRMLeD8396J1+tl//792Gx2oqIiSExMRq1W0dVlpa2tDa1WR0rKGARBoKuri4aGBlQqlWIwnU4XtbU1fvX1PNRqDW63m6qqSkBk8uRJ6PUG7HYHpaVluN0ucnNz0Wp1eDwuqqqq8Xi85ObmoNFIbZO1tTV4PD7Gj8/DbDYq1Xy56KJWq2lra6e5uRmdTvL49HodHR1SGsJgMJGUJCko7d+/30+ajyM2VlJkqqjYB2jIzs4mKioKp9NJbW0NWq1eUZWXQndJtSk7OxudTktJSTFOp4tx48aRkBBPfX0jZWX70Gi0isGXevj7z2NaLCZ/g0fvD2hZf1KnC/YyZa0Dq/UIN5Rb7g35XNXxDw17vZDXcrlcol6v5z//+Q92u525c+cqBw+9oRz6OAi9XkdsbJTCc5P6w4dmeGJiIv3hsZQfHKr3GAip7dHg76t243C4lKq6nMeURka4hkRJeu/JZ9A9s5aN1k7iVXqMYxJQ17fSIfiIELVE3XElM885k+RkaQ5NQ0M9DQ2NRERE+g2hCqfTTm3tftRqDePHT8DtlvQvi4v3IIpQUDDZ30Ek0XBUKhWTJ0/yd+142Lp1G263h3HjsoiIiFAKOl6vl6ysLKKiopQZPILgY+zYscpru3f/iM8nkJOTS3x8LBqNhh07duBwOMjJycVkMuPzeaiqqsTj8ZCbO14Z2FZTUxOUo92/v5aWlmaioqLIysoCJE/U6XSRnJxMcnISHo9X4QjLKkiyp2s0mpg+fQo6nY6Skr1YrZ1ERUWTkJCI2+2hpESax56bm4fFYvYzE5xER0eTmZkZ9L2EWi0HqV/bYNABEvVIo1EPaZjZ4WQohc13h3yuesajw14vVCjFnFdeeYW6ujruvrv7Rg+1oZSEVEM3EiqVSsk9iqL0dG5paR92a1l0dAQul8efxxwatUe+P8kIGpAmP/ZvBOUecq1Wg8MxuFD8zyefyQktXWzRqYlNz2Fzey1pGjM+fJKsWG46BQvm0twsDSaTFcI7OztoaKjHYDAp4a0oitTUVKPTaRk7NgOvV0AUBYUMP3ZsBlFRkej1OvbsKcblchMfH09kZBQqlYrq6ko8Hi+xsXHExcX6dTprcLulFEFcXLyyhsfjJjo6lvj4OFQqFZWVFXg8HuLipGKQwaBn3z5p1EVcXDwGgwGfT/S/JpCePhaTSfLaqqur8Pl8jBuX45/700FrawsqlQp5HG5jY6M/3RBBaqqUE92/vxq73UFq6hiSk5Pw+XyUlpZhsznIzc3BYDDQ1dVJU1MTer2B1NRUv2xelSK4nJiYhE6nJTY2FpPJTEaGxMeMjDRjt7sGVQSSctZ6RUzD5XIPak8fVobyf3eGfK76pN8Pe71QoSouLhYjIyNZt24dHR0dPPDAA8rBw8VQBkqauVyS9+jxeElIiKWtrWNYhlIURaKiItBqNXR12YbUGSGNNpWLM5LBHUwFXaYuyb/vcPT9h6bTafnxyy9puOoWzIKITQ0xz6+m8FRpumNJSQlNTY2YzRaioqSe8vb2NurrG4iIsCjUIpvNTmnpXlQqFQUFBajV0kOhrKwUjUZNRkYmZrMZnU5LcXGJf055MvHxUohbXi5NgkxKSlIq5XKOMSYmVhmLW1Ozn9bWFmJiYsjIyAAkD7C1tYXo6GjltYaGOmw2O7GxMUohaf/+GqXwIxv6vXv3+nOWySQlJSlCH9LI3gx/uO2ipKQEkIRI9Hp9EIF9+vRpGAwGKioqlM9l7NgMxfuVyPBpJCTE4/FIBtnr9TF+/AQEQaC6upq2tna/RzkOlUrK65pMJiZOzMPnE4dED9PrtXi9PgwGvb/v3BNS++xhZSi/vCPkc9Wz/jDs9UKFatWqVaJer6eoqIiMjAyeeOIJ5eBoNpSSdyZ5jyqVCrvdicPhDPK4hkvrAanaqFKpsVik/u7B5A9lao9arVKoPcPprJDfs9EojXFwOiV6kUoFBoNBoRG9dNPt6IuKld9zTMpjySPBiW9ppngzoCIyUuq8am/voKpKUi2Si3per5ddu3YhijBx4gSMRhMmk559+/bhcrlIT5faDSXidyUdHe0kJiYplKG6unqam5uIjIwkIyMDt9tNXV0d+/2zddLSJPJ3c3Mz9fW1mEwRSrW7o6Od9vY2IiMjSU1Nw+Px0tTUTFOTVFySRHl10mzxtjZ0Oh0JCXIqoQGrtROLJYLU1DF+z7YGl8tJbGwM8fGJiKLIvn37UKkgOVky8qIoUFJSgtvt9RPYNXR0dNDS0oxarVE80bq6OhwOOxERkSQnp6DVati/v9pPak9Gq9X7GQFSg0FOTi5arTRGNzk5ddB8zJ4TGLVaDQaDDq1WO6Aox2FlKP9ze8jnqk/547DXCxWqjRs3iomJifzwww98//33PP3008rBQ20o4cBxED1FKex2Z59P1cHQeiA0ao/ctSMpsxy4dmBxxuv14XA4R4RYLBthjUat6G4OpQBktXaybdt22tpaiYiQNrrH46K6ukaZqaNWS33mNTVV+HxeUlNT0Wolo9zS0ozdblNUi9xul5/bASOPtAAAIABJREFU2OwfWCZ5jjabjbq6OsxmE+PHT0Sn0+FyOamrk+TeCgoK/OIsdioryzEaDUybNg2Px0trazt795b4c6CSIK/NZqesrBRQUVAwGb1ej8vlpKam2n/fUmOEnGcMFCKur5eMt9kcwfjxeRgMeqqqKmlubiEqKpoxY8YEeY/p6enExsb4i1KSgIbsiTY1NVFXV4/BYGDChPFoNBqqq6v8aYU44uLisdnslJTsRadTM336dNxuH11dXfh8PmJiYhSPu//vu/fRExLXVqIYyepGPUU5DitD+e9bQz5XfdqaAc8JF5Qc5e7du/ntb3/L3/72N+XgaDKUPUUpAosffSEUWg8MrTgjjxNVqzU4nU58PgGj0YBWq/GLJAx8f0OFXq9T1pI8ShU6nXZY69rtdmpra3G73cTFSURzQRAoKSlGp9OSnz8JtVqN0+mmrKwUu92h8BoBGhsb/ET1ONLS0nG5JC+xuroavb6bRiQbm0CjJ3MZtVrJkBiNBrq67Gzfvg1BEMnNzcVsNitFHlGE7OxxmM0WfD4fu3f/6D8vB4tFarksLd2rdAv5fD46OzspKSnB6/UwduxY/6RIaba6IEi5TZVKoiq1trZisZgVOTapem4jKiqapKRkRFGkoqICEP3TMSNwOBxUV1ehVncLiXR0tPsFPEyMHZsBiOzevcefFpLmlUuUMA15eXn98jGNRkmHob/WR8nL1Cv7Qhbl8PkOjphGWAzlxhtDPld95tMDnxQmKIbS5XLxz3/+M6hx/FAbyvLyfXz88XtkZIzl0ksvx+Hom6LTG3pWq3vDcKg9sqK5FPJKJO+uLseIGMjAQlBv3mygJyulB0LjiPaEx+Ohpqaauro6tFpptK7bLXEMBUEkKSmZiAgLKpWKhoZ6rNYuoqIiqd63g+Nmnk17eysdHZJ8VVbWOJxOFw6HncrKKkRRJDt7HGq1BlGUvl+v10deXg5RUVF+LmMpNpudtLQ0Rde0uroKh8NBSkqKEq7W1NTgdDpJSpI8Wem8ahwOO8nJSURFxWA0GmloqAMgNTWNmJg4zGYj+/fXIAgCqanpGI1GxVCDyPjx4zEajUqHU+BYicbGBhoaGjEaTeTlSR6r3LUUFxdPWlqqYswFASZMyCMxMZH29nbKysrQaNRMnJiPx+OloqJSmVOVkZGl5DGNRiNpaWl+1kM3TCZ9yD3ispdpMOjw+QSsVic+38hbyrAYyn/eEPK56rOfHfZ6oULV2toqfvDBByQnJ3PuuecGHZQN5UDivRs2vMcHH7yLRqNh2bLlnHzy7GHf2G9+cyc7dmzjZz+7gPnzFxIXFz/oawRXqyWEg9rTW3FGEASMRuOwDdVAa/VXyAGUfKU0sU9qewyVXhTYEeR0utizpwSrtZOYmBg0Gg0gGwZbkDe5ZfO/8bV/jctwDKeefg4gFYOqqirQ6/X+ENqDw+H05ztF8vLyiIuLRa/XsWPHDjo7uxgzJlXxZsvKyvx95skKhamiohKrtZP4+ISgtsnOzg5iYmKV3GbPVkqdTkd9fR1WaydxcfHExyfg8XgoLS3D6XSQmJikrLFnj+TxpaWlERcneXx79uxGFCEnJxuzWSKw79tXis8nkJWVhdfrpa6ujpqaavR6A9nZErezvFxqB42Pjyc2Ng5B8FFRUY7RaGTMmFR0Oh1NTc00NDRgNBrJzJTyn11dVqKjYzj22GOUh4A8gXEwjgLgjzS8h49H+Y9fhXyueu5zw14vVKhWrVoldnZ2Ehsby4QJE7j00kuVg7Kh/OKLz/nyy/9wzz0PsHPnDl5/fa0i3tvS0swtt6zgxRdfw+1286tfLefFF1/rd/JdKKiqqiA1NR2VSmCoslRRURGK0ZIwdGoPoGgNSormkgHqWZwJNlTCgP3c/a013EKQHKJrNJp+O32C0wjB50m0nRq6uqxYLBble5Ur07GxcXTVfcmJk7R88lUrqghJKGPs2LF+I+2grKwMvV7H+PHjlXxnefk+bDYHCQmJRERYEEVpAqXTaScmJlZ5MDY3N9Le3kFERIRS7W5ubqK9vR2j0ajkQFtb22hra8Vg0Cs9352dVqzWTkwmI2PGjMHjEejsbKOrqwuDwUBycgput5empiba29uwWOTWTCmP2dnZSXR0DMnJSf5iUDUOh4OkpERiY+PQ6SRjKOc/JY/Sx549e3C5PEpqob6+nsbGJiwWs8IqKC8vRxSlkcCxsVLRcc+ePQiClP80Gk3+vnIv0dExZGdnKa2Pg8UQ6cODRlgM5SfXhXyuet6fh71eqNC+//77fPjhh3i9Xu66664gQymjP/He3bt3UVg4Fb1ej16vJy1tLGVle8nPnzysG8vIyAJklfOhQR4CJXuNMsVlMJD7wg0GyVO02Rz9Fmek8EgyNnI/t8ViUniQ/a/VXdWWC0HDGY7WPYJV8hRjY6P8XqkTQRCVtfoT6JB4klL0sHv3brZu3YrRaMBgMKLRaNmy+Qs0jiLamiykJkZT73KRMW4CVVUV6HQGpk2bSkFBIXq9nrKyvXi9XjIzM5g8ebKfprOXri4riYmJTJ4saW0GVsonTz4NkIsyUuFn9mwpYpHbJnU6HSedNBOpFddKbe1+DAYdJ5xwAm63h5aWVsrL96FWqznuuBP8RREHtbX7iYgwMX36HNxur0JW12p1ilZhY2Mj9fUNmEwmRfVdpjmZzWbOPfdcjEYDpaV72blzF3FxcYwfP/GAfvSUlBQ6OjopKioKCufr6+uprq4mLi6W4447FpvNTlHRDlwuJ2PHpmM0WujsrKa6uoKEhETGjcsZ8n44LBBmPcpwQZufn6/00no8Hnw+nxJmyehPvLfnMbPZPOCgnoMFj8dDZGQEBoN+0EYnsGDidLpob7cOOvcoGyq5n9tsNipGNHC2jk6nxWg0oNNp/ZMbrWGdv+Lz+ejqsis51ZgYSWnI4/HS2dkVMvk5Pz+f/Px8mpubaWlpRqPR0lq3ndTcQrq6HJRUtdJmK2LmrNNJS5NI2GVle9Fo1GRn55CVla2Qt91uF1lZmUybNkVpU9yxYwdRUVFkZmYyZkwKbW1tFBUVKcWg+Ph4vyHZjiCI5OXlkpmZicPh5D//+Q8Wi4UpUwrIzMzAbnfwySefYLFEMG7cONLT0/H5BL766ivi4xMYM2YM8fGJCIKXL7/8ivHj84iKiiQnJxu73cnmzV8zcWK+f1xwNl1dXRQVFZGTk0NcXCypqWPo6uqkoaGeqKhooqJiiY6Oo66ujq6uLvR6A+PHT0ClUtHY2IjBYMRsNjNp0mT/92HDYrGQlJREcnIyXq8PrVZLenoq2dlZeDw+fD4BjUaDXq8nIsKE1Tr4/Pdhp/M7SmXWtNOnTwdQKBRut/uAObgWiwW7vbvHVBRFZbZOz2N2u33AQT2DQaiqzr1Re3w+gfb2zh6enbPPhHggT1EQBJxOJ52dwxt3C9K0PqvVpninMTGSZ+fzCRgMcjXTOWRCeygIDK/tdieiKFXpIyMtOJ19fya9ISFBGv7V0dHB1BMk/dEkFWQjedPbt29n6tQp6PVacnNzcTpdFBXtICcnB5VKrXTA1NXVolJpiYgwM25cJunpadTW1mO3O9Dr9URFRTN5cjQulzPogTNlyjQ0Go1f39FAVFQkp59+ipJmkR8ykycXHHDvhYWFB7w2aVJBDy4qHHfcMf7USvd5gR5/X3louUUyEKEOL5MhiioiIizo9ZJIhrx/3G4fUlCkVu7niMPByhMMEtozzjiD8vJyHnroIZYuXdqr+kl/4r35+ZP5y1+ew+Vy4fF4qKwsP6jhgWwgpTnLvRdnuj07rb/FsTsUlo2+ySR7dJ5BeVmDgZyzBJVC9/B6fdjt/YfzQ0Ww4T8wvHY63UGfiUQvcvY5SbInoqOjg0jTcg5X6hMPNr4yUTsQ48Z1vyZX9eV5QQNFAD3bO9vbrcMi8wemTGQP32w2+QeIefypJR0ul3tI0UWokMZISGt5PF4EwYder8dsNqFWy4U5H6BCEIZekBytEA9CdX4oUHk8HrGtrY2tW7dy2mnBMlE9q959ifdu2PAeGza8hyAIXHHFVZx22pyw3WB/4yCGSu0JbAkEAvquR26sp06n84tgaILC70Dj4nA4hzyTJxDBepYDV8rhwFxsqFV7ub8+MK86XKPfX3FJ+ryGJxgSKgwGnd9Aqf26pQNzcoeKwD3ZUyQF+hZREcX+GRwHawwEhKeY41u/LORzNRf9v2GvFyqCZubYbDZ0Op1S2TzUPEoIFu8NB7VHUhnvNgiiKCpP73BRemT05D72R9UZjgCGjP4MTOj3HEwv6suzCzTGLpd7SOMtBkLgGoIgKLlxh8M1omMUent4hbsdVYZWq8VsNvjnMw38nUl8Wb3yYJLbWCXIDkP338ZhZyjfvCLkczWXvDrs9UKFFqS84kcffURxcTHLly8nNTX1oN1A6FCh0agG7T3KCKTKuFwu2ts7lRBTNg6RkREIgk8R1Rgqenp0VqttQAMszVbpUjyLwAp1fwZooPB6sDiwam/AYjErobRs0GWD0dbmOCiTAaXRsENrDAgFPfvoe7IbXC43Lpdb4bXK389gBU5kBCvcO0NWLhcESZHfbndiMMi5d7P/IdztZUrzpuSfwwfiaK1633fffX7KRibnnntukHDvaIE8/Gywfx9y8t1g6N+IBBoHg0GPxWJiKKFw9+ZX+41I56CNiM8nKBVqo9FAdHRkryFtT2M8EnnVwKq9xWLGbDb5vUzHkKYIhgrZGMv8z8DPsTfjPRxDHZhyCOVz9Hp9Ad+PnshIy6CI/QNpBQwG8rwdjUYy3jExUco44QkTJgBqP1948LS4Qwb3QXJ/BwmtwWBg+vTpzJ49e5QaSS0S4Vz0/0GIAxrMnnSbzs7Q6Tay5yD/sUpFDmefquiBxtjnGz73UYaUN5UUkSTjbQZEPB4vWq1GCa+HYoxDRWCawu32YLM50Ou1mM0m9Hp92AU/QhEs7psbOjiNx8ARHBL9qzPkIhZ057UdDpcSrVgspl7pX/JDz2SSquWhRBiDgeQB27HbobW1mUcffQSPx8v8+ReyYMHiYTd/HEyMVo9ywLneow9yp06wwexJ75BDkeHakOC2vu6qcOAf2kjl6GQE/qGB7AEP35vqC4H5UqlX+8C8WbiKUAdSbgangDTY1EPPoW4DNQEMBj3z3/IDV37QDJRGGR5Ef8itRhRh+/ZtfPLJh1x88eWMG5c9QmsGIxw5Ss8LS0M+V3ftumGvFypUgv/R11veZ3QaShmSwSwtLeHTTz9h4cKFTJyYP6xZM/1BrVZhMhkV3uPBqJT3Vb3uWUgJRUkpFAxlBEWg996bNxXKewvXewhOfQQ/SOQHTThC3oEgpyokURHRz90dmX0iXXN05CTDYiifvyjkc3XXrx/2eqFCOxKJ8YOBjo5OVq26g9ra/SxatIjY2ASs1i5GYqPIyjwGg6Qx6fMJ6PWSBqDP5wt7BbZn9bpneC132sjeWExM5JCr9tI1usU8bDb7oN6PVITyHkCm76vI0VOeLpzFoN7CcrlaPhIhb08E5lalAo0LrVbrL/4Ywxp5jCYDGU6MVh6l9lDfwFBhMpm58spfMH36sf4uITmHKaBSiYRj4wRy15zOYKKxzUZQ7rA/AeFQ0E0lMoZcvRYEUamAyl02oXpMPVMHwyVRC4KAzebw34v+gHs5mB5doIhKYLVc6vUPv6Hsj9vZ/SAJLs4F03pCR7CBPIyKNKFilLYwHoY5yoEg+n+kD3ywDrOk5Sfl30IZAAaBXpJ60HmvoZDDB3MvPVMDgXzAkU4dGAxyR4mqV8pNuNEfFzUwNxmesRz4c6NGpYgX6nszGHQYjfK9uEPgu8r5x9FtIMMRersevzDkcw23vTvs9ULFYetR9o3ADRV6pVxKxBsxGELnPsroS/yityKIjIHC66Gi93txA6JSmZd4oiObozMajej12qAWwMhIS1hzqjJ6FpY6O/uulvfkQQ62wBKYqnC7BycqIiOQ1hNYuT8wXdFdoBnNBjKsOBp6H2x0G8v+DGY4DVZP8Yuef4xDCa+HCqmH3InZbFSq5W63F7vdMSJ97HCg19bW1i3y4XZ7g3Kqw+2wCfwsB+OtBvMgJUGNwEFtfaFni2E4+r0DVZ3kdMWrr76G0Whg9uzTMRhkcZqfgIGUMUoN5REYevcHqVLe2trMP//5D046aSZTpkwZMYMV2ActtV+qcbvdww6vB0Jvgr0gKqFiOI10qC2PPTHUnu2eBPFwfJb9heWDbTEcDkQRdu3aydq1L1FcXMyyZVezeHHodJlDjXCE3s5H5od8rvGeD4a9Xqg4gj3KA+FyuVm9+nd89dV/Oeuss4iIkJSlR6piqNNp0em0iCJ+vUEpJ6nRqMNuKEMxWIEE6VBk5/pDTxGNwVaUexL7ZbJ2X+mKYIK4e9AE8f7QW1ju8fjQaKSoZDAthkNBYIFm8uSp/PGPz1BVVUl9fd2IrTlqMWI80+HhJ2Uo1WoNxx13AjfddDtRUVGAEBCSh6dSPlB4LRspSept+GpBQzFYvcnODZRTlXGgwRpeCBrY4x6cO5Q8xZ4E8ZHU7PR6fXi9PnQ6LVqtlBf0+Xwj6kH2RfHJyMgkIyNzRNYdzTjamTOq0V0pH6rBHGz1ergh3VCJ3r3hwHD2wAJHT4M1UtVy+UFjMhlRqVQK7Wgk6USBKZKefNTgCnV4VIOOVIpPOEJvxz3zQj7X9Mgnw14vVPykPMq+0f1El7mYEBq1aKjFIK/XS2enN0gtKJSKcHduDxwOJ1br8EPCYA5kd4HD4XAq4WgoYrrDRU+D5fV60eslkRK1Wh1249zzAdHbCA65Qh0YlkuTPQfbjnh4UHwONcJFOBcEgQceeIDi4mL0ej0PP/wwmZndHvrDDz/M1q1bsVgsADz33HP9TmY4aiiDEBq1KJzVa1ktSKadyF02drtT8UiDh44NvnsmVMjtdi6Xm4gIM1FREQHiD86BLzBE9FdRdjhcQVSncNCLegpihPJwC6yWm0yDIY7/BCk+w0GYQu+NGzfidrtZv34927Zt47HHHuP5559Xju/atYsXX3yRuLi4kK531FD2it6pReXl0hzoWbNODru0mdRl48DhcPi9ugi/tqAkwDHSIwjgwPxjV5fNn1IwYjJFh73q25Mg3pfB6qZdBT9MBtvXr9NpMZuNChl/KPlOURQD9CD1QVM2g8cXywZSw5HUYjjiCFMx5/vvv1emdU6bNo2dO3cqxwRBoLKykvvuu4/m5mYWL17M4sWL+73eUUPZL6QN/u23m3n11Zeord3PDTfcOKLSZpIitaCM2hVFEUEQ8Hp9I2Yke+YfA8N5QThQUHhooWc3QiGI9wb5YWK3OzAa9URESO2jA9GLgqlITlyu8KQPeor5vvfehxQXFzN//oVkZeVw1EAOHuEq5nR1dRER0T0dVqPRKJNj7XY7l112GVdddRU+n48rrriCgoICJk6c2Of1jhrKELB79y4WLlzCqaeejlarDlBGD0+lXLpWz3C+Ox8YqjZmuNbrDYGCwnLoORivbqgE8b4gtf4F04sCPd6eLYYjla4A2eO1M2PGLOrqGrjxxhVMmTKdhx/+/YisdyRD9ITHGYiIiMBmC3zgC8rkWJPJxBVXXKFMm50xYwZ79uw5aiiHi2XLlvfyanioRT1bJ3sL57spNBrMZiOxscZBT0yUMVhF754IDD2DvbrePbXhrjcQehuhISsGjZTyeyACK9gxMfFcddW1XHrpVezbVzZiax7JCFcx55hjjmHTpk3MmzePbdu2MX78eOVYRUUFN998M++//z6CILB161YWLlzY7/WOGsohQ07KD021KNT8XCB8vt7ng4fSr9yTThQOwna3V6fz8zG7xzPIhmskCOIDQVYMUqvVI0Lul9foq4Kt1+uZODE/7GseDFRWVnDNNcvYsOGfvY6uHmmEK/Q+66yz+Oqrr1i6dCmiKPK73/2OtWvXkpGRwZw5c5g/fz5LlixBp9Mxf/588vLy+r3eUR5l2CAG/fRFLZKr1+EYtxqofN7XqNj+Zm2HGxqNBovFpHQjSV0/zhEjiEP/ikHhV0o6sik+NlsXDzxwD7t3/8g773w0aEMZDh5l+y9OC/ncmBf/Pez1QsVRjzJs6FuEw2q1kpiYgNlsxOPx0tVlD4sKe29zdURRxOl0olZr+pwoOBIIHJpls9lRqzXKHPNwjwGG7gKUStV3QSiw2BKopCSR84fCgTxyKT6iKLJ69SNcc80K7r77tkN3H6O0M+eQG0pRFFm4cB7p6WMBKCiYwnXX3RB0zssv/4Wvv/4SjUbLjTfeyqRJBYfiVkNEN3G9qGgrf//7On74YStvvrmOuLj4Eatcu1xuvF4vFouZiAiL32C6sNtHjv/YkyButQY/AKTq9OAEhQfCUASAA1WdBqcIf2RSfD766H3Wr38z6LWUlDHMmXM2eXnj+/itgwPxqHBv76ipqebpp9ewevUTvR4vLt7Dn/70JE899TwNDQ385je/5sUXD97g86GiqGgbq1f/joULF3PuufMwm02Es6c8EL0NAwvOEYbWxx0qggniUo50oAfAcIZ69ddiOBT0FBDpOWq2+2OSPcgjw0D2h4suWkBiYhIAP/64k/z8yfzpT38d1DXCEXo3X3pyyOcmvPHVsNcLFYfcoywu3k1zcyMrV16LwWDgxhtvJSMjSzleVLSN44+fgUqlIiUlBZ/PS1tbG7GxsYfupkPAlCnTeP31t3q8Gl4Rjv74iMHamN38x1CMWl8YSgFKRl+Cwv1V7uV7Nxgkg9xbi+FQEDjHXZakW7PmcaKjoznvvPnExMTyUzGQMtavf1/59+LFP2PNmmcPyX0cnZlD7y7/rbfeyWWXXcUZZ5zJ9u3b+O1v7wvyGG22LqKjY5T/N5st2Gxdo95Q9g45tyX4vZbBi3D0LOAMlH+U+rjt2O0qf6V88IPIhkoQ7w09xY27K/fdIiJDaTEcKtxuqZf7/PMv4G9/e5OlSxdz+eVXcumly0ZkvaPoH0dzlMD55y/g/PMXBL3mdDrRaDQATJ06jebmJj9hWDIeFksEdns3cdRutxERMXwX/9BCNpiqkKlFgeFuXwIO/SGY/yjnDYU++9QPJGyHtyDUcxiZ3LKpUqlQqTigQ2gkEEjxycjI4c477+Waa1bQ0tIyouuOdrz99oeHbO2DRSEbLA556P3yy38hOjqaSy9dxt69JSQlJQfNGC8snMrzzz/NxRdfTmNjI4IgEhMT088VDyf0VC06kFo0nHC3LwSGnT21MUeaIN4TUoumxBJQq1X+1xgxD3Igik9sbByxsaEJJYwmdHV18dvf3ovdbsPj8bBy5S0UFEw51Lc1aIzW0PuQF3M6Ozt56KF7cTgcaDQabr31TjIzs3juuac47bQ5TJpUwEsvvcDmzf9DFEVWrryVnJzcATfFk0/+kaKibZjNZgAee2xNUO/n6IRkKF0uB5999i+KiiTVE6fTPWyB3/4giUWY0Golz14q/gy+62cw6NnSGMgB7UmOD08h6sjmQL700gtERkayZMklVFVV8MAD9/Dyy28c1HsIRzGndv4JIZ+b+sG3w14vVBxyQzkUhLIprr9+OY8++vhh531+/vlGnnhiNfn5k7jiimUUFBQOeuTuYBA43lb2KGW5s3BPSwRJMk4q0Aw886Zne+dgJyZK6MmBPDILNFarFb1eh8FgZN++MlavfoQ///nlg3oP4TCUNfOOC/nc9E++G/Z6oeKQh95DwZIll6DX6wCpOKDXB3cQCIJATU01q1c/QltbC+edN5/zzw99aNGhxLhx2Tz77F/IzMzyvxL+cRWA4s3JpPXAPGWgnFm4BngFKsBLknEDtzT6fN2FKKOxWwMylBna3fnHI89A9lYUXbXqfvLzJ9PS0sxDD93LjTceOtL4cDBaeZSj3lAOZVM4nQ4WLVrC0qWXIQg+Vq68jokTJ5Gb238/52jAuHHZPV4JrJQPz2AG8hG93r47hHrTxuwZHocKOccqTYR00tbmGHQY3VsHUl8toP3NoTlS0FtRFKCsrJT771/FihU3MX36sYfgzoaPoznKMCNwU8ycGUxS9fl8uFxOzGZZ5v0psrNzmTv3vENxq2GGoPw3VIMZGMIOVUtysONl9XodZrMRmVIU7hyrbIB3797Nl19+yTnnzPPTyI68/GMoKC/fxz333MGDDz56yLprwhF6V86ZFvK5mZ9tG/Z6oeKw3FHl5fu49947uf/+hw8wkgDV1VVcf/1yfD4fXq+XoqLtjB/ft9bc4QW1/0eDKGoUpZzeoNVqiYqyEB0diSgKtLV1YrPZh0TalsNlm82BwaAnNjYKo9FwQP5UnitjMhmw2520t3eOSCHK4/HS2WnFZLJQXl7ORRctZs2aP+J2j2xP+2jFCy88i9vt5qmn/sgNN1zDXXfdeqhvaUgQBDHkn4OJw9KjvOuuWykt3UtKyhhAEul87LE1rFv3OunpY5k161TefPNVPv98I1qtlrlz57FgQf9S74cvglWLfD4vdXW1FBRMVgQjRsJQ9SSFg1TFDkeL4cA4sILd3NzMv/71D+bPv1BhOhzFwUU4PMqyWYUhn5vz5Y5hrxcqDktDORgIgsDjjz9GaeledDodd911ryLAAbBhw3t88MG7aDQali1bzsknzz6Edzt02GxWNmx4n3feWU9+/iRWr/7DiE5MBKnF0GyWKtggdbnY7Y6wtBn2jsAKNhymAZGCgfbm4YZwGMrSmaEL3uR+vXPgk8KEUV/M6YnArp1Q8N///hu3280LL6xl584dPPvsEzz22BoAWlqaefvtdbz44mu43W5+9avlHH/8iej1+pG6/RHDu+++TWlpCQ8/vJqJEyfickkh6EhQi3p6k62tHQABlWkvdvvgBn/1jyNTxae/vTkceL1e1Go1avXh9yA52sIYBnz++UbKy8toa2tjwYJFIVWxi4q2ceKJMwFK5dpoAAAYyUlEQVQoKChkz57dyrHdu3dRWDgVvV6PXq8nLW0sZWV7yc+fPGLvYaRw+eVX9fJqeKlFgZzL3loM5cq00agnMtI8bGm1I5niA/3vzaFi794S1q17jWuvvYGkpORhX+9gY7QaysPikWO323nyyT/wyit/JSEhkYyMDH7/+4dZvfqRAb0Wm82GxdLdkaNWq5Xf6XnMbDbT1dU1Mm/ikEAq+oAaUVT7Dc/gN6LBoCcmJgqz2ai0UfYnk+Z0umlr68ThcCrKRZKqe2iQClSygZR/jiwjCf3vzcFCzqB9/fWXfP75RkpK9ozoaOORgugTQ/45mDgsPMp//esflJfv46GHfq8QsRcsWMwjj9xPScmefoV8LRYLdrtd+X9pTra212N2u53IyMNdcKM3DH6+z4FTEwc/xbA3abX+ZoP3N4fmSER/e3OwUKlU2O02mpqamDXrVL75ZjOTJhUQFxcfrts9KBitohijfjfa7XZ2797F9OnHkpmZhSAI+Hw+9Ho9Wq2O/ftrlHPlGdiBKCycyubNksDnzp07yM7OVY7l50+mqOgHXC4XXV1dVFaWM25czsF5Y4cEsgHqm1qkVquwWEzExkah0Wjo7Oyis7NrWMpBsrRaR4cVjUZNbGwUINDY2IhkvAM9SA2HwbYMC/rbm/1hz57dNDTUAwR5oHV1dUycmM+KFTdTXl5GWdneERQXGRl4vaH/HEyMeo9SFAV+/HGnQu9RqVRKkvr661cqpHKXy4nBIM1QCSz4nHLK6WzZ8g3XXXc1oiiyatX9QTSixYuXsmLFLxEEgWuu+RUGg9S18uijD1JXV4fH42bZsuXMmnWqck/r17/Bhx9+oPSR//rXq4LEhkc/Aqk1kpdZVlbCp5/+HzfddBOiyIhMTQycDb53715uvvkmpk07hksuuZxJk0KnhRwp6G1vhoItW76hoaGO22+/G61Wq+z3nJxccnIkYzt5ciFfffVfJkzIJyoqeiTfRlgxWrMFo54eVFu7n1tvXcmTTz5HSkoKXq83KDzp7OzgP//ZxLfffoMoilx99TVBbYDyjOfB4OOPN1BaupebbrqNzs4OrrzyEt5992Pl+G9/ey9Lllxy2I4kDURNTTVr1vyeffvKuOSSS/j5z5cqcmcjg26Kj91u46OPPuTDD9/nySefIz4+YQTXPXzRrdMpfS9bt37HXXfdRl7eeB57bE1Qukje7zU11fzhD79j+fJrmTIl9G6X4SAc9KAf8kNvDJm+e8+w1wsVo9qjFEWR1NQ0Tj31dP75z0+44oqr0Wq11NRUU1VVyeTJBbz77t/ZtWsHy5YtZ9++Mp566o/ccccq0tLSARQj6fP5FIHggXD66Wdy+ulzlHvQaII/puLi3bz++lpaWlo46aRZfVScDw+43W7OOmsuc+ac7adFBVbKw7dO9+O4m+JjNksKUEuWXBK+hQ4SBoo6wgHZU5T3cH19PSkpKTQ1NZKWlkZGRiaRkZFBzoNarUYURdLTx5KSMoZ//etTcnPzlMhrtGO0epSj2lDKT9DTTz+TP//5GTZv/h8nnDCDTZs+45RTTsPlctLR0c5ll11JYeFUf87nf9TUVJOQkMD/+38vc+KJM5k6dbpiJEPhYcqdHXa7jd/85k5++cvrg47PmXM2F164BIvFwqpVt/PVV/89bInq2dk5ZGcH5mXDJ8IBR65IxaeffkJUVAz33vuQEnWEy1DKe1Tep6Wle3nppReoqqpg4sR8brvtbgoKpnDLLSvo6GgnOjomaF/L/1669DKamxsPGyMJo9dQHhZZ84kT83nyyee49lppjO3tt9/FFVdcTXJyCtXV1Uo3Q3t7O8nJKXR0dLB58/947bW1bNjwLr/4xRVs3Spp1wVupv7Q0FDPypXXcc458zj77LnK66IosmTJJcTExKDT6Zg5cxZ79xaPxNs+xBgetehIp/icfvqZ/PKX1wG9Rx1DgSxlJ+9Rr9fLP/7xMU8++QfOOmsub7zxNpWVlfzrX/8gMTGJCRPyeffdvwMEFTFlD3TcuGyOP37GsO/rYEIQQv85mDgsDKW8gaZOncaVV/6CwsKp6HQ6/5NUUHI0//vff6muriInJ5cPP3yfK6/8Bffe+xAzZ57MunVvsGHDe7z11pt4vV5lMwqCcIDRbG1t4dZbb+D661ceoGNps9m44oqLsNvtiKLI1q1bmDDhSBHc6A2yweyulPdnMHuvYB85BlKG2WzGbLb0GXUMBXLU83//95HyYI+NjaO5uVlR51+wYBFbtnxDZ2cnCxcu5uOPN/D22+twuUIf/zuacbTqPQzIG6hnYSY5OYXp04/j0kt/TnZ2DhqNlnPOORen00lNTTV//OPTAJSWlqBSQU5OHs88s4aGhgZWrrwlKG8pP5HVajWvvroWq9XKK6+8yCuvvAjABRcsxOFwMH/+hVxzza+48cbr0Ol0HHfcCcycOetgfhyHCD3n+0ifl/S8ObLHLPSFhoZ6Vq26g4ULFwdFHaEgcL/J2LJlM2vXvkhm5jg0mp1s3PgpixdfxJw5Z/HDD99zwgkzmDv3PDZt+ozvv/+Wc86Zx1lnzaW4eA9nnjkXk8k0qPbe0YjRGnqP+qp3KLBarWzf/gOFhVOIjo7hjjtuYuLESSxffi1btmxm/fo3WbHiZsaNy+ajjz5g795ibrnl1zz44G+YNKmAM888O2wDpa6++lIlJ5SamhZE+ThSBDgkBKoWwU/JQIIUdaxceS233PJrjjsu9Dkv8p+bbNBcLhdqtRqdTsf7779DbGwcJ588m4ceupfa2v3cc8+DVFaW88UXm7j44svJyxvPK6+8SGVlBXfffR9arXbU9HSHo+r9RVro0dmp+49WvUOCTDCPjIxk1qxTAGhubqazs5OLL74MgH/+8x8UFEwhMzOLurpayspKGTs2g7a2VrZs2YxOp+Pzz//FmDGp3HHHKkwmk3L9wVTKQdr0oijy7LN/OeDYkSTAIeHIKcwMBb1FHY8//jQGg7Hf35MNZGtrC3/+87NUV1dy7LEn8ItfXMe2bVupqanm9ddf4fTTz+SWW+6kqamRSZMK+P77LXz66Sfk5Y3n0kuXIQg+Ze+IojzFcnQYzOFgtHqUh7WhVKlUBxiyhIQEXnhhLSCpoO/dW8KCBYtRq9WUlpbQ0tLMVVf9kr/+9Xlmzz5NIa0/+OBvqKmpDlKHHkylHKTqpNPp5JZbVuDz+bjmmhUUFEhE6iNJgOMo4Oabb+fmm28f8LyeD1u32837779NZWUFWVnjuOyyZVxzzVUcc8xxTJkyjc7ODlavfhKtVsuOHdtZvfoRXnzxNU44YQYGg9Qvr9Fo0Ol0SioqsEJ+uGO0dhId1oayN8gbUxAE/1jb35Gamk5zcxPffvsNubl5aLUaNm/+ihdffJWIiEjUajV1dbU0NjagUqnYvPkrtm3byvnnz+e00+YcQLvoC0ajkYsvvpyf/WwB1dVV3H77jbz55jtotdqfgADHUfQG2UgWFW0jOTmF5OQUHA4HJSV7WLBgERkZWVx99TU8//zTPPPMX/jii8959tknyM3N44MP3mXx4qUYDIYg6pHsOR4JHmRPHOwiTag44j5peWPKmygjIwutVovRaCIpKYkzzzyHjz76AJPJhE6nQ61Ws39/DQ0N9UyYkM/dd9+O0Whk3rwLeOedt3jnnbeUa/dXKQcYOzaDc845F5VKRUZGJtHR0bS0NAM/JQGOwwNtba1ceOF5VFZWhO2avWkN7Nq1k5tuup5nn32Sp59+nKeeepwzzzyH1NQ06uvrAFiy5GKqq6uorKzgwQcfJT4+kW3bfmDZsuXMn39h0PWPdBylBx1iREREcPnlVxEXF8/7779Damoaf/vb62zfvo3Vqx9h0aIlfPHF56SmprJ48VLOOONMbr31TjZt2gjAm2++ynfffauIovbmWX788QaeeeZJAJqbm7DZbEpb3k9PgGP0wuv1snr17w4YczxUyDnCwC4aGd9++zUnn3wKf/nLK9x00+1s2/Y9TU2NJCenUFJSrBjLc8/9GQ8/fB8xMTFcfvmV/OY3DypeZM8C0JGM0Wooj7jQuy/IbV4vvvhnxoxJY9WqB1i79q+88cYrnHPOPObN+xnXXnsVS5deqvzO1q3fERsbR0dHO9988zW7d//IK6+8yIQJ+Vx77YqgQowgCJx//nweeeQBrr9+OSqVirvvvo+3317XrwAHwCeffMgnn3wISDms0tISPvjgU8XjfPLJP1JUtE3pGHrssTUKr+4oBo9nn32SBQsW8dpra8NyPdmAeTweXnttLVVVlZxxxlmccsppbNq0kV//+jcAJCUlM3PmLL74YhMXXLCQl176M8XFe0hJGcNNN91GdXUV0E2DC8xB/lRwtJhziCH3wk6bdgyzZ59GVFQUN910W1Cf7EknzaK9vV35nY0b/8Glly5jw4b3SUhI4PzzFyjh+fbtP3D88Scq58ob+v77Hw7a2IWFU5V/X3DBQi64YOEB9zZv3s+YN+9nADz++O8577wLgsLy4uLdrFnzrKJWdBRDxyeffEhMTAwnnjhzyIayZ5VZEAReffVlrFYrXq+H2bNP4y9/eQ6328WMGSfz0UfvK0W94447geeff5qbbrqN5OQxmM1m5Xpjx2YEXfdIzEEOhKOGcpRAphHJuaRAJaJp047l2Wef4IsvPic6OgaTyUxu7ng2bdrIjBknM3lyIXq9noqKffh8Pux2G99+u5miom3MmXMOkyeHPhipN+zZ8yPl5WXcdtudymuCIFBTU83q1Y/Q1tbCeefNP6Bb6ChCx8cfb0ClUvHdd9/6Zwzdx2OPrQlJuUguFMpV5v37a1Cr1YwZk0ppaQmVlRU8++xfiI6OITIykieeWM3q1U9y7bVXMWvWqURGRvHmm69y3nnS97dixU2KMZQfrj8l77E3HDWUowy9Pa2nTp3GX//6/9i0aSOJiUkUFExhw4b3MJnMZGVlo9fr+frrrwCRGTNO4v777yY3dzwTJ07mmWfWcNZZ57Bo0UVB1xwMF/PVV9dy9dXXBL3mdDpYtGgJS5dehiD4WLnyOiZOnBTSvKCjOBB/+tNflX/fcMM13HHHqgGNpPwdyt+j1Wplw4Z3WbfuDQoLp3L++fP51a9u4pZbVmCxRODxeDjhhBkIgoDb7ebuu+9l69Yt7N69i7POOpcFCxYB3Uo/P3XjGIhwTTkWBIEHHniA4uJi9Ho9Dz/8MJmZmcrxt956i3Xr1qHVarn++us5/fTT+73eT8+37wdyT/npp59JQcEUvF4vbW2tjBuXw9ixGQD89a/Pcc01K/j88438+OMuNm3a6FcRuo+iou34fD62bdtKWVkpEMzF7A9Wq5WqqkqOOea4oNcNBiNLllyM0WjEbLZw7LHHUVpaEu63fhT9QP4Ov/9+C/fccwdPP/04ABs2fMrMmSfz6aefkJycQmbmONavf0PhOBYUFGI0Gpk9+zRWrLiZ559/mQsv/Dnw0yrQDAbhKuZs3LgRt9vN+vXrue2223jssceUY01NTbz22musW7eOl156iTVr1uB2u/u93k/Wo+wNPT0/rVbLsmXLcTqdGI1GfvxxJ3v3lnDeeRfwyCMPcN11N5Cdncvrr7/CN998zYQJ+XR2dvLZZ/+ipGQPsbGxHH/8DBYtWqL8QXz00fvMmHEyCQmJQWtt376V4447/oB7qq6u4v777+bll99AFEWKirYzd+75I/ch/ITQWweVTPEJ3AstLc088cRqdDo9Z5xxFn//+zqMRqmvevLkAnbt2sFXX/2Xq6++hjvvvJnW1hZ27txBWlq6MglRvl5g+H4UByJcoff333/P7NlSi/C0adPYubN7BnhRURHTp09Xmj8yMjLYs2cPU6ZM6fN6/RrKcPRuHs7oFuGIRBAEJk7MZvXq1SQmRjJ58kQqKvaydOkiTjjhCRoaGmhpaaG2dh91ddVccslSTjnlFK677jomTMhm1qxZbN26lbfeepMrr7zsgD+UlpZ68vKylc987dq1ZGRkMGfOHC68cCErVixHp9OxePGFnHjiwVGsHo144YUX+Pzzz/F4PFx88cX8/Oc/H7G1duzYweTJk4mK0lNVVcFpp53G0qWLMBo1/Pvf/0YQ7BxzTAEVFcfy9ddfsHr1aqZMmYLN1skTTzweFOodRWhY6guPZGFXV1cQM0Sj0SiF266urqBiqcViGbD542jo3Q8C85hqtZrk5GQuuOACAAoLC/nuu++4++67WblyJS+88AIZGRmUlJSQm5vL4sWLSUpKQqVS0dXVxXfffceKFSsUhZeerv4vfvELrrzySuX/r7rqKubMmaMce+edd1i3bh0XX3wxANu3b+fyyy8HoLKykosvvphLLrmE+++//wDSs9PpZOXKlVxyySX88pe/pLW1Neyf1cHAN/+/vbMNaXIN4/iPxuzozPkyXzCyuUzUXohNDtJAC4oFEQ0sFUmiRVKEUCmUUUYZzC+JKymrRaaxRX3pBfqifvCFxFoUZBTZcic84Wps0RxmMnc+7Lg285yjHkvP6fl9e3Y93LufZ/C/d1/39dLTw5MnTzCbzTQ3NzM4OPivx7RYLAwPDweuvV4v165dQ6PRUFdXR1VVFb29vZSVldHV1QXAli1bcLvd9PT0IBaLycjI+PNw53e0Wi19fX0sXbp00gB0gR9DZGQkHs/XvvNjY2OBg9uJNo/H84/JH4JQToNgP2NOTg7Nzc2kpqayYcMGSktLefXqFf39/eTk+IuldnZ2EhMTw5o1a0hLSyMyMpK1a9dit9tDYjC7u7sxGAwYDIaQ7J2/4vLlyxw7dixQg1Cv13PgwAFMJhM+n4+2traQ+81mM+np6ZhMJrRaLefPn5+N1/HD6erqIj09nf3797N3717WrVs3o3F8Ph9NTU0UFhbS0NBASUkJTU1NgH/R6ezs5MaNG1y6dImUlBSMRiOrVq0iIiKC1tZWRCIRarWa27dv43Q6WbFiBdXV1SxZsoRNmzaRkJCAzWabNABd4MegVCrp6OgA4OnTp6Snf63hsHr1ah4/fszIyAhutxur1RpinwzhV5wGwdvlsbExwsPDKS0tZevWrSQlJWG1WvF4POTl+TMq7t27R3Z2NomJiVy8eBG1Wk1FRQUJCQmBcd69e0d9fT1JSUkMDQ1RWlrKmzdv/nYeKSkpnDt3LnD9/Plzfv3VX+orNzeXBw8ehNwf7K/Jzc2lu7v7372IOcLlctHb24vBYODkyZNUVFTMKK3v7t27PHz4kJqaGoxGI0ePHmX58uV8/vyZ9+/fMzY2RkREBCKRCK1WS3R0NC9fviQ/P58LFy4AUFhYSH5+PjExMX+myP7C6OgoAEajEblcPpuPLjBNNm7cSFhYGEVFRej1eiorK7l69SptbW3Ex8dTUlJCcXExO3fu5ODBg4Hkj79COMyZIeOhHfBVQLdv345arSYsLAyLxcLr16/ZtWsXDoeD9vb2gLgFHxZIpVKqq6tRKPydI8vLy7Hb7YHrydBoNAwMhPYzH5+DRCLB7Q6tIxrsk5nM/l8hOjoahcIfpqVQKFi4cCFOp5O4uLgpj/Hp0yfu3LmDTqcjNTUVn8+HUqkMvMORkRHS0tKwWq1kZWURGxuLy+VCIpGgVCq5f/8+DocDmUzG5s2bQ8YWi8Wz/cgCM2TBggWcOnUq5LNly76mDBcUFFBQUDD18WZtZj8hweWtxkUzOTkZgNjYWAoKCsjMzKSlpQW5XE5GRsY3J6oSiSQgihaLBalUGiKCUyF4e+fxeIiKigqxB/tkJrP/V1CpVHR2duLz+bDb7QwPD087Wyk8PJy3b98GxDU4T/vjx49IJBJEIhFXrlzB4XDQ0tLC8PAwcrmcuLg4GhsbkcmEtro/G4JQzhITT7EVCgVFRUUBm81mw+HwVxIKdvD7C2T8Rm1tLZWVlSQkJAS2yVMlKyuLnp4eADo6OsjODo3FVCqVtLe3B+wqlWp6DzdPWL9+PZmZmWzbto19+/ZRVVU1rcLK4M/HDvZfBS8ytbW1DA4OsmfPHhYvXsyRI0cwm82UlZWRnJwcWAzH420Ffh7+thWEwMyZmHFht9tJTEwMuefs2bO8ePECu93OypUr0Wg0qNXqKY0/MDDAoUOHuHnzJv39/Rw/fpzR0VEUCgWnT59GJBKh0+loaGjA6/Vy+PBhPnz4gFgs5syZM8THx//zl/xPaWxs5NGjR5w4cSLEX1xeXo5KpaK42N9n3OVyERMTM1fTFJhHCEL5nQkuoDAunkNDQ9TV1XH9+nWqq6u/ayygwLc4nU70ej0ul4uKigpkMhkNDQ309fVRX1/PokWLQha66bYEEfj/IQjlHPHlyxdMJhOtra3IZDJ27NjxzZZZ4PvhdrupqanB4/Fgs9nIy8tDp9MhlUrnemoC8xBBKOeAidvyW7du8ezZM3bv3i1kc/xgnE4nUVFRgWDkiS2RBQRAEMo5Y7KcYoG5w+v1/nRFcgWmjiCU8wDhX4yAwPzmD5HhOmQcrNoZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "std_scale = StandardScaler().fit(X)\n",
    "X_norm = std_scale.transform(X)\n",
    "\n",
    "pca = PCA(n_components=2)\n",
    "X_norm_2d = pca.fit_transform(X_norm)\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "\n",
    "for i in range(meta_info.shape[0]):    \n",
    "    ax.scatter(X_norm_2d[i][0], \n",
    "                X_norm_2d[i][1], \n",
    "                y[i],\n",
    "                c=colors[i],\n",
    "                edgecolors='black',\n",
    "                marker = markers[i],\n",
    "                linewidth='0.2')\n",
    "lr = LinearRegression()\n",
    "lr.fit(X_norm_2d, y)\n",
    "print(lr.score(X_norm_2d, y))\n",
    "\n",
    "x_surf = np.linspace(min([x[0] for x in X_norm_2d])-0.5, max([x[0] for x in X_norm_2d])+0.5, 20)                # generate a mesh\n",
    "y_surf = np.linspace(min([x[1] for x in X_norm_2d])-0.5, max([x[1] for x in X_norm_2d])+0.5, 20)\n",
    "x_surf, y_surf = np.meshgrid(x_surf, y_surf)\n",
    "new_X = [[x_surf[i,j], y_surf[i,j]] for i in range(20) for j in range(20)]\n",
    "new_y = lr.predict(new_X)\n",
    "\n",
    "ax.plot_surface(x_surf, y_surf,\n",
    "                new_y.reshape(x_surf.shape),\n",
    "                rstride=1,\n",
    "                cstride=1,\n",
    "                color='None',\n",
    "                alpha = .1,\n",
    "                linewidth=0.1, edgecolors='black')\n",
    "\n",
    "ax.set_xlabel('Component 1')\n",
    "ax.set_ylabel('Component 2')\n",
    "ax.set_zlabel('Median NRMSE (TS)')\n",
    "ax.view_init(elev=30, azim=45)\n",
    "\n",
    "plt.colorbar(sm)\n",
    "plt.savefig('plots/pca_lr.pdf')#, bbox_inches=\"tight\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-04-01T22:31:43.567828Z",
     "start_time": "2018-04-01T22:31:43.558211Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-3.0568729899430225 6.397333634170555\n",
      "-4.886680447346294 10.699770233275661\n"
     ]
    }
   ],
   "source": [
    "print(min([x[0] for x in X_norm_2d]), max([x[0] for x in X_norm_2d]))\n",
    "print(min([x[1] for x in X_norm_2d]), max([x[1] for x in X_norm_2d]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "## Regression plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-03-25T20:07:41.422502Z",
     "start_time": "2018-03-25T20:07:41.060703Z"
    },
    "hidden": true
   },
   "outputs": [],
   "source": [
    "sns.set(color_codes=True)\n",
    "markers = [\"o\", \"^\"]\n",
    "\n",
    "plot_data = pd.DataFrame(index=meta_info.index.values)\n",
    "\n",
    "plot_data[\"Nature\"] = meta_info['Nature'].values\n",
    "plot_data[\"Component 1\"] = X_norm_2d[:,0]\n",
    "plot_data[\"Component 2\"] = X_norm_2d[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8,8))\n",
    "\n",
    "sns_plot = sns.regplot(x=\"Component 1\", \n",
    "                       y=\"Component 2\", \n",
    "                       data=plot_data);\n",
    "plt.savefig(\"regplot.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "sns_plot = sns.lmplot(x=\"Component 1\", \n",
    "                      y=\"Component 2\", \n",
    "                      hue=\"Nature\", \n",
    "                      col=\"Nature\", \n",
    "                      data=plot_data, \n",
    "                      markers=markers, \n",
    "                      size=8, \n",
    "                      aspect=1, \n",
    "                      palette=\"Set1\");\n",
    "sns_plot.savefig(\"regplot_by_nature.png\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}