MIAMI trial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "\n",
    "os.chdir('C:/Users/rfuchs/Documents/GitHub/M1DGMM')\n",
    "\n",
    "from copy import deepcopy\n",
    "\n",
    "from sklearn.metrics import precision_score\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.preprocessing import LabelEncoder \n",
    "from sklearn.preprocessing import OneHotEncoder\n",
    "\n",
    "import pandas as pd\n",
    "\n",
    "from gower import gower_matrix\n",
    "from sklearn.metrics import silhouette_score\n",
    "\n",
    "\n",
    "from MIAMI import miami\n",
    "from init_params import dim_reduce_init\n",
    "from metrics import misc\n",
    "from data_preprocessing import gen_categ_as_bin_dataset, \\\n",
    "        compute_nj\n",
    "\n",
    "\n",
    "import autograd.numpy as np\n",
    "from autograd.numpy.random import uniform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "###############################################################################\n",
    "###############         Heart    vizualisation          #######################\n",
    "###############################################################################\n",
    "\n",
    "#===========================================#\n",
    "# Importing data\n",
    "#===========================================#\n",
    "os.chdir('C:/Users/rfuchs/Documents/These/Stats/mixed_dgmm/datasets')\n",
    "\n",
    "heart = pd.read_csv('heart_statlog/heart.csv', sep = ' ', header = None)\n",
    "y = heart.iloc[:,:-1]\n",
    "labels = heart.iloc[:,-1]\n",
    "labels = np.where(labels == 1, 0, labels)\n",
    "labels = np.where(labels == 2, 1, labels)\n",
    "\n",
    "y = y.infer_objects()\n",
    "numobs = len(y)\n",
    "\n",
    "# Too many zeros for this \"continuous variable\". Add a little noise to avoid \n",
    "# the correlation matrix for each group to blow up\n",
    "uniform_draws = uniform(0, 1E-12, numobs)\n",
    "y.iloc[:, 9] = np.where(y[9] == 0, uniform_draws, y[9])\n",
    "\n",
    "n_clusters = len(np.unique(labels))\n",
    "p = y.shape[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#===========================================#\n",
    "# Formating the data\n",
    "#===========================================#\n",
    "var_distrib = np.array(['continuous', 'bernoulli', 'categorical', 'continuous',\\\n",
    "                        'continuous', 'bernoulli', 'categorical', 'continuous',\\\n",
    "                        'bernoulli', 'continuous', 'ordinal', 'ordinal',\\\n",
    "                        'categorical']) \n",
    "    \n",
    "# Ordinal data already encoded\n",
    " \n",
    "y_categ_non_enc = deepcopy(y)\n",
    "vd_categ_non_enc = deepcopy(var_distrib)\n",
    "\n",
    "# Encode categorical datas\n",
    "le = LabelEncoder()\n",
    "for col_idx, colname in enumerate(y.columns):\n",
    "    if var_distrib[col_idx] == 'categorical': \n",
    "        y[colname] = le.fit_transform(y[colname])\n",
    "\n",
    "# Encode binary data\n",
    "le = LabelEncoder()\n",
    "for col_idx, colname in enumerate(y.columns):\n",
    "    if var_distrib[col_idx] == 'bernoulli': \n",
    "        y[colname] = le.fit_transform(y[colname])\n",
    "    \n",
    "enc = OneHotEncoder(sparse = False, drop = 'first')\n",
    "labels_oh = enc.fit_transform(np.array(labels).reshape(-1,1)).flatten()\n",
    "\n",
    "nj, nj_bin, nj_ord, nj_categ = compute_nj(y, var_distrib)\n",
    "y_np = y.values\n",
    "nb_cont = np.sum(var_distrib == 'continuous')\n",
    "\n",
    "p_new = y.shape[1]\n",
    "\n",
    "\n",
    "# Feature category (cf)\n",
    "cf_non_enc = np.logical_or(vd_categ_non_enc == 'categorical', vd_categ_non_enc == 'bernoulli')\n",
    "\n",
    "# Non encoded version of the dataset:\n",
    "y_nenc_typed = y_categ_non_enc.astype(np.object)\n",
    "y_np_nenc = y_nenc_typed.values\n",
    "\n",
    "# Defining distances over the non encoded features\n",
    "dm = gower_matrix(y_nenc_typed, cat_features = cf_non_enc) \n",
    "\n",
    "dtype = {y.columns[j]: np.float64 if (var_distrib[j] != 'bernoulli') and \\\n",
    "        (var_distrib[j] != 'categorical') else np.str for j in range(p_new)}\n",
    "\n",
    "y = y.astype(dtype, copy=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#===========================================#\n",
    "# Defining hyperparameters\n",
    "#===========================================# \n",
    "\n",
    "r = np.array([3, 1])\n",
    "numobs = len(y)\n",
    "k = [n_clusters]\n",
    "\n",
    "seed = 1\n",
    "init_seed = 2\n",
    "    \n",
    "eps = 1E-05\n",
    "it = 50\n",
    "maxstep = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#===========================================#\n",
    "# Defining acceptance rule\n",
    "#===========================================# "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The user specifies a range for each variable to lie in: \\[min\\_value, max\\_value\\].\\\n",
    "An observation is accepted if all of its variables lie in the specified range.\\\n",
    "A set of ranges is called a rule.\\\n",
    "Several rules can be specified in order to select different areas in the observations space.\\\n",
    "To select a particular value for a non-continuous variable (let's say 0), the range is set at \\[0,0\\]. \\\n",
    "Example of rule for 3 variables, with the first 2 variables that are continuous, and a categorical variable:\\\n",
    "\\[-1.3, 2.5\\], \\[-5.1, 5.1\\], \\[0,0\\]\\\n",
    "Observations with variable 1 in \\[-1.3, 2.5\\], with variable 2 in \\[-5.1, 5.1\\] and with variable 3 taking modality \\[0,0\\] are selected. \\\n",
    "If one wants observations that can take more than one modality for the categorical variable, let's say 0 and 2, one can specify the following second rule: \\\n",
    "\\[-1.3, 2.5\\], \\[-5.1, 5.1\\], \\[3,3\\]. \\\n",
    "The two rules \\[-1.3, 2.5\\], \\[-5.1, 5.1\\], \\[0,0\\] and \\[-1.3, 2.5\\], \\[-5.1, 5.1\\], \\[3,3\\] are passed to MIAMI to properly select the observations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_rules = 5\n",
    "# The rules are stacked together in a numpy array of shape (2, number of rules, number of variables in the dataset).\n",
    "\n",
    "# Specify 4 impossible rules: All observations have to be equal to 0 and a very permissive rule, all observations\n",
    "# required to be in [-10000, 10000] to test the algorithm\n",
    "\n",
    "authorized_ranges = np.empty((2, num_rules, p_new))\n",
    "authorized_ranges[0] = 0\n",
    "authorized_ranges[1] = 0\n",
    "\n",
    "authorized_ranges[0,0] = -10000\n",
    "authorized_ranges[1,0] = 10000\n",
    "\n",
    "target_nb_pseudo_obs = 500"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "[-3527.350312558943]\n",
      "1\n",
      "[-3527.350312558943, -3494.2783579420393]\n",
      "2\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969]\n",
      "3\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114]\n",
      "4\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816]\n",
      "5\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255]\n",
      "6\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255, -3463.6173028883454]\n",
      "7\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255, -3463.6173028883454, -3466.854975981989]\n",
      "8\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255, -3463.6173028883454, -3466.854975981989, -3466.5504112053186]\n",
      "9\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255, -3463.6173028883454, -3466.854975981989, -3466.5504112053186, -3465.695277648792]\n",
      "10\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255, -3463.6173028883454, -3466.854975981989, -3466.5504112053186, -3465.695277648792, -3465.461631084396]\n",
      "11\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255, -3463.6173028883454, -3466.854975981989, -3466.5504112053186, -3465.695277648792, -3465.461631084396, -3467.800399205382]\n",
      "12\n",
      "[-3527.350312558943, -3494.2783579420393, -3488.103158087969, -3468.711538639114, -3467.064254697816, -3464.2390609269255, -3463.6173028883454, -3466.854975981989, -3466.5504112053186, -3465.695277648792, -3465.461631084396, -3467.800399205382, -3467.868826564085]\n"
     ]
    }
   ],
   "source": [
    "#===========================================#\n",
    "# Running the algorithm\n",
    "#===========================================# \n",
    "\n",
    "import matplotlib .pyplot as plt\n",
    "prince_init = dim_reduce_init(y, n_clusters, k, r, nj, var_distrib, seed = None,\\\n",
    "                              use_famd=True)\n",
    "out = miami(y_np, n_clusters, r, k, prince_init, var_distrib, nj, authorized_ranges,\\\n",
    "            target_nb_pseudo_obs, it, eps, maxstep, seed, perform_selec = False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting utilities\n",
    "varnames = ['age', 'sex', 'chest pain type', 'resting blood pressure',\\\n",
    "            'serum cholesterol in mg/dl', 'fasting blood sugar > 120 mg/dl',\\\n",
    "            'resting electrocardiographic results', 'maximum heart rate achieved',\\\n",
    "            'exercise induced angina', 'ST depression induced by exercise relative to rest',\\\n",
    "            'the slope of the peak exercise ST segment',\\\n",
    "            'number of major vessels (0-3)', 'thal']\n",
    "\n",
    "y_new = out['y'][len(y):]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtuUlEQVR4nO3deXwU9f3H8dcHCKYiIgK2IErUQi2QcMilthqxRcGzFioqKmi11qtaSz1+tuBNa63WW6yKF4hipRSx4AGIBTkNIOBdLAhFQEFAoRyf3x8ziZt1s5lAZjdk38/HI4/sznxn5rOz39nPzHdmvmPujoiI5K462Q5ARESyS4lARCTHKRGIiOQ4JQIRkRynRCAikuOUCEREcpwSQRaZ2XVm9tfdZb47y8y+bWavm9kGM7sj2/GUMrOzzGxStuOojJlNMbOfZzmGEWZ2czXOb5GZFVfX/MJ5LjWzH1XnPKtLTdsmk9XqRGBml5rZHDPbYmYjksYVm9kOM9sY/i03s2fNrGum4nP3W919lzbw8HMsr+75VrMLgTXA3u5+VbaDKeXuT7t7r2zHUdOY2UAze6Ma5/eNJOLu7dx9SnUtoybZTbbJcmp1IgBWADcDj1Y03t33AhoCPYB3gGlmdmx1LNzM6lXHfGqBVsBi192LshvIye3W3Wv9H0EyGJE0rBhYnqLsvcCcCuZTADjBHu4KYCVwVcL4ocAY4CngC+DnQAtgHPAZ8AFwQVL5pxLe9wCmA+uA+UBxwrh9gcfC5X4OjAUaAF8BO4CN4V+LFPM9GVgUzncK8P2EcUuB3wALgPXAaCA/HNcUGB9O9xkwDahTwbo5ApgdzmM2cEQ4fASwFfhfGN+PUkzbBPhHuM5mh9/XGwnjDwVeDmN4F/hZwrgRwH3Ai8AGYCZwSNL3VS+h/BTg5+HrgUnLceAi4P1wHd8HWML484Al4biJQKs0de454L/h+ngdaBcl5nD8jwl2StYT1MeppTGnWE43YE647lYBfw6HvwhcllR2AXBqus8KfB/YDGwPv691EWNO+R0RbCuJ3/8/Eurdj8LXdYHrgA/Dec8FDqjg81ZWl68FFoef6THCuhyOPxEoCaedDhQlTXt1uI62ANcDY5KW/Rfg7vD1oLAubAA+An4RDq8x22SVfiOr44e2pv9RtUTQM/wSG6QYV0CwAY0Kv/BCYHVChR4aVvpTCY62vkWwEd8P5AMdw/LHJpR/Kny9P7AW6BNO++PwfbOEDXs00BjIA46u6HMkzbcNsCmcXx7wW4KEVD+h0s0KK+u+YeW+KBx3G/BgOF0e8EMSfhgTlrcvwYZ3NlAPOCN83yThR+TmNN/PM+HfnkBbYBnhD3S4npcRbHj1gM4EzUztEub9GcEPYj3gaeCZpO+rKolgPLAPcGD4XR0fjjs1XG/fD5dzPTA9zWc6j+BIcw/gLqAkYVy6mJsS/Kj3Ddf5lcA2Kk4EM4Czw9d7AT3C1z8DZiaU60BQn+pH+Kzl1kuEmKN8RzcnzW8pX283g4GFwPcIklEHwrqTNE2Uuvw2cABBnfxX6XLDmD4FuhMknnPD8nskTFsSTvstgqPYLwmaMwmnWZmwfk8ADgnjPTos27mmbJNV/avtTUM7YwXBl7tPmjI3uPsmd19IsNdxRsK4Ge4+1t13EGzUPwCudvfN7l4C/JXgBzPZAGCCu09w9x3u/jLBnl4fM2sO9CaoDJ+7+1Z3nxrx85wOvOjuL7v7VuBPBBX9iIQyd7v7Cnf/jGDPvGM4fCvQnGDPd6u7T/OwNiY5AXjf3Z90923uPopgj/akyoIzs7rAT4Eh7v6luy8GHk8ociKw1N0fC+c9D3ie4Iey1N/cfZa7byP4gerIzhvm7uvc/T/A5IR5/QK4zd2XhMu5FehoZq1SzcTdH3X3De6+heBHoIOZNYoQcx+CZrQx4fd1F8GRRUW2At81s6buvtHd3wyH/x1obWatw/dnA6Pd/X8RPmtFKoo5yneUzs+B6939XQ/Md/e1KcpFqcv3uvuysC7fwtfb5gXAQ+4+0923u/vjBHv+PRKmvTuc9it3/xiYR7ADAMEO4pel69fdX3T3D8N4pwKTCH6Uo8jENlklSgTftD/B3tK6NGWWJbz+mCBzpxrXAvjM3Tckld8/xTxbAf3MbF3pH0ESaU6wl/KZu38e9UMkxfBx6ZswQS1LiiHxh+ZLgj1LgNsJ9lQmmdlHZnZNlGWEKvqcyZoR7EUmrrfE162A7knr5SzgOxHi3xkVzasV8JeEGD4j2GH4xmc0s7pmNszMPjSzLwj28CDYMahsOS1I+PzhRp64PpKdT7CH+Y6ZzTazE8PptgDPAgPMrA7BD+KTET9rRdKtm8q+o3QOIGgWqkyUulzRttkKuCopxgOoeNsFGMnXieTM8D0AZtbbzN40s8/CefWh/Pe7q59jV7fJKlEi+KafAPPcfVOaMgckvD6Q4CiiVGJ2XgHsa2YNk8p/kmKey4An3X2fhL8G7j4sHLevme2TYrrK9gZWEGwEAJiZhfGniqH8jIM92qvc/WCCvftfV3AivdwyQhV9zmSrCZo+WiYMS1y/y4CpSetlL3f/ZYR5l36HeyYMi/rjlGwZQTtwYhzfcvfpKcqeCZwC/AhoRNBEBUHiqMxKEj5/wveVkru/7+5nAPsBfwDGmFmDcPTjBD/IxxLszc6IsHyovE4lq+w7qmx+ywiaWSoTpS5XtG0uA25JinHP8Oi1VHKczwHFZtaS4HdhZLjcPQiOeP4EfNvd9wEm8PX3WxO2ySqp1YnAzOqZWT5B+15dM8tPdUWABfY3syEEh6nXVTLr35nZnmbWjqBddHSqQu6+jOCk1G3hsosI9uCeTlH8KeAkMzsu3KPMDy9Da+nuK4GXgPvNrLGZ5ZnZUeF0q4AmSc0OiZ4FTjCzY80sD7iK4JA41Q9YOWZ2opl9N6yoXxCcQNyeougEoI2ZnRmu89MJ2vrHV7YMd98O/A0YGq7TQ4FzEoqMD+d9dvi588ysq5l9P8K8VxNsXAPCdXoe0X5wUnkQuDb8zjGzRmbWr4KyDQnW8VqCJHRrFZbzItDOzE4L6+rlpEleZjbAzJqFe5XrwsHbAcIf/h3AHXzzaCCdVUBLM6sfsXxl39Eq4OA00/8VuMnMWofbYpGZNUlRLkpdvsTMWprZvgTbcem2+TBwkZl1D5fRwMxOSNpJKyesP1MImn//7e5LwlH1Cc79rAa2mVlvIPEy5JqwTVZJrU4EBCf0vgKuIWiD/yocVqqFmZWe2Z9NcPK32N0ru8loKsHh2avAnyopfwbBHuEK4AWCtvCXkwuFSeMUgsq7mmAPZjBff0dnE7QPvkNw0uuKcLp3CE5efxQe8rZImu+74We/h+AE3knASUltxRVpDbxCsH5mAPd7imu/w/bcEwkq9FqCk18nuvuaCMsAuJRgz/m/BD9Yowg2DMJmtV5Af4J1+F+CPd89Is77AoL1uBZoR4SNLRV3fyFc7jNhc8/bBOdtUnmC4ND/E4IrWN6soFyq5awB+gHDwphbE5z0rMjxwKKwHv8F6O/um5NiKSTY0YjqNYIrWv5rZpV+hxG+o0eAtmH9HJtiFn8m+HGcRPDj9ghBm3nycqLU5ZHhfD4K/24Op51DUBfuJbiQ4QOCk+KVGUlwZFfWLBR+3svDmD8nOAIclzA+69tkVVk1nGfIGWZWAPwbyAtPmO3q/G4EWrr7ebs6r9rEzP4AfMfdz812LLs7MzsHuNDdf5DtWKTmqu1HBDVWeGjXliCx5DQzOzRsDjAz60bQfPZCtuPa3ZnZnsDFwPBsxyI1mxJB9swjOEH6cLYDqQEaEpwn2ERwuH0HweWPspPM7DiCJsZVJDRriKSipiERkRynIwIRkRy323Wu1LRpUy8oKMh2GCIiu5W5c+eucfdmqcbtdomgoKCAOXPmZDsMEZHdipkl3/1fRk1DIiI5TolARCTHKRGIiOS43e4cgUhtsnXrVpYvX87mzZsrLywSQX5+Pi1btiQvLy/yNEoEIlm0fPlyGjZsSEFBAcHN5iI7z91Zu3Yty5cv56CDDoo8nZqGRLJo8+bNNGnSRElAqoWZ0aRJkyofYSoRiGSZkoBUp52pT0oEIiI5TucIRGqQgmterNb5LR12QqVlVq1axZVXXsmbb75J48aNqV+/Pr/97W/5yU9+Uq2xjBs3jsWLF3PNNdXydMWU9tprLzZu3Bjb/KNYunQp06dP58wzzwRgzpw5PPHEE9x9991ZjSudnEoE1b2RVUWUDVIk09ydU089lXPPPZeRI4NOSj/++GPGjRtXyZRVd/LJJ3PyySdX+3zjtG3bNurVq9rP5NKlSxk5cmRZIujSpQtdunSJI7xqo6YhkRz22muvUb9+fS666KKyYa1ateKyyy4DYPv27QwePJiuXbtSVFTEQw89BMCUKVMoLi6mb9++HHrooZx11lmU9mRcUFDAmjXBg83mzJlDcXExACNGjODSSy8FYODAgVx++eUcccQRHHzwwYwZM6Zs+bfffnvZ8oYMGZIy7lGjRlFYWEj79u25+uqry4276qqr6Ny5M8ceeyyrV68G4O6776Zt27YUFRXRv39/ADZt2sR5551H165d6dSpE3//+9/L4uzXrx8nnXQSvXr14vTTT2fChAll8x84cCDPP/88S5cu5Yc//CGdO3emc+fOTJ8ePPzummuuYdq0aXTs2JE777yTKVOmcOKJJwLw2Wefceqpp1JUVESPHj1YsGABAEOHDuW8886juLiYgw8+uOzoYdOmTZxwwgl06NCB9u3bM3p0yqfi7jIlApEctmjRIjp37lzh+EceeYRGjRoxe/ZsZs+ezcMPP8y//x08S+mtt97irrvuYvHixXz00Uf861/pnqj5TStXruSNN95g/PjxZc1FkyZN4v3332fWrFmUlJQwd+5cXn/99XLTrVixgquvvprXXnuNkpISZs+ezdixY4Hgh7Nz587MmzePo48+mhtuuAGAYcOG8dZbb7FgwQIefPBBAG655RZ69uzJ7NmzmTx5MoMHD2bTpk0AzJgxg8cff5zXXnuN/v37l/0A/+9//+PVV1+lT58+7Lfffrz88svMmzeP0aNHc/nll5ct64c//CElJSVceeWV5WIfMmQInTp1YsGCBdx6662cc87Xj+d+5513mDhxIrNmzeKGG25g69at/POf/6RFixbMnz+ft99+m+OPP75K6zgqJQIRKXPJJZfQoUMHunbtCgQ/zE888QQdO3ake/furF27lvfffx+Abt260bJlS+rUqUPHjh1ZunRplZZ16qmnUqdOHdq2bcuqVavKljdp0iQ6depE586deeedd8qWV2r27NkUFxfTrFkz6tWrx1lnnVWWLOrUqcPpp58OwIABA3jjjTcAKCoq4qyzzuKpp54qa+qZNGkSw4YNo2PHjhQXF7N582b+85//APDjH/+YfffdF4DevXvz2muvsWXLFl566SWOOuoovvWtb7F161YuuOACCgsL6devH4sXL670M7/xxhucffbZAPTs2ZO1a9eyfv16AE444QT22GMPmjZtyn777ceqVasoLCzklVde4eqrr2batGk0atSoSus4qpw6RyAi5bVr147nn3++7P19993HmjVrytq03Z177rmH4447rtx0U6ZMYY899ih7X7duXbZtCx7jXa9ePXbs2AGQ9nr2xOlLm5XcnWuvvZZf/OIXFU5XlYdplV5K+eKLL/L6668zbtw4brrpJhYtWoS78/zzz/O9732v3DQzZ86kQYMGZe/z8/MpLi5m4sSJjB49mjPOOAOAO++8k29/+9vMnz+fHTt2kJ+fX2k8qWIvjTHV+mzTpg1z585lwoQJXHvttfTq1Yvf//73kT9/VDoiEMlhPXv2ZPPmzTzwwANlw7788suy18cddxwPPPAAW7duBeC9994raz6pSEFBAXPnzgUol2SiOO6443j00UfLrvz55JNP+PTTT8uV6d69O1OnTmXNmjVs376dUaNGcfTRRwOwY8eOsvMNI0eO5Ac/+AE7duxg2bJlHHPMMfzxj39k3bp1bNy4keOOO4577rmn7Mf5rbfeqjCu/v3789hjjzFt2rSypLh+/XqaN29OnTp1ePLJJ9m+fTsADRs2ZMOGDSnnc9RRR/H0008DQTJt2rQpe++9d4XLXbFiBXvuuScDBgzgN7/5DfPmzat0He4MHRGI1CCZvrrMzBg7dixXXnklf/zjH2nWrBkNGjTgD3/4AwA///nPWbp0KZ07d8bdadasWVl7fEWGDBnC+eefz6233kr37t2rFE+vXr1YsmQJhx9+OBBcDvrUU0+x3377lZVp3rw5t912G8cccwzuTp8+fTjllFMAaNCgAYsWLeKwww6jUaNGjB49mu3btzNgwADWr1+Pu3PllVeyzz778Lvf/Y4rrriCoqIi3J2CggLGjx9fYVznnHMOJ598MvXr1wfg4osv5qc//SnPPfccxxxzTNlRRFFREfXq1aNDhw4MHDiQTp06lc1n6NChDBo0iKKiIvbcc08ef/zxtOtj4cKFDB48mDp16pCXl1cuYVen3e6ZxV26dPGdfTCNLh+VmmbJkiV8//vfz3YYUsukqldmNtfdU17HqqYhEZEcp0QgIpLjlAhERHKcEoGISI7LqauGluafmcWlr8/iskVEKqYjAhGRHJdTRwQiNd7Qau5CYGjlR6K33HILI0eOpG7dutSpU4eHHnqoytf/A4wdO5Y2bdrQtm1bAIqLi/nTn/4UuefNuLpvLigoYM6cOTRt2nSX5rMr1q1bx8iRI7n44ouB4Eaxyy+/vFxne9mkIwKRHDZjxgzGjx/PvHnzWLBgAa+88goHHHDATs1r7NixkfrbqUhp982lunTpUiP78C+9g7gq1q1bx/3331/2vkWLFjUmCYASgUhOW7lyJU2bNi3r56Zp06a0aNGCV199tdyDaV5++WVOO+00ILjb9//+7//o0KEDPXr0YNWqVUyfPp1x48YxePBgOnbsyIcffgjAc889R7du3WjTpg3Tpk0DKu7aOl33zRs3bmTQoEEUFhZSVFSUsuuKV199lU6dOlFYWMh5553Hli1bysbdfvvtdOvWjW7duvHBBx+Uxda+fXs6dOjAUUcdlTa2KVOmcMwxx3DmmWdSWFjI1VdfXe6HfejQodxxxx1s3LiRY489ls6dO1NYWFjWtfU111zDhx9+SMeOHRk8eDBLly6lffv2QNAfU+ln69SpE5MnTwaC7rBPO+00jj/+eFq3bs1vf/vbshgHDhxI+/btKSws5M4779zJb/9rSgQiOaxXr14sW7aMNm3acPHFFzN16lQg6INoyZIlZf35P/bYYwwaNAgIunru0aMH8+fP56ijjuLhhx/miCOO4OSTT+b222+npKSEQw45BAge7DJr1izuuuuusi6hK+raOl33zTfddBONGjVi4cKFLFiwgJ49e5Ybv3nzZgYOHMjo0aNZuHAh27ZtK9cdw957782sWbO49NJLueKKKwC48cYbmThxIvPnzy97EE+6brdnzZrFLbfcwuLFi8t1TQ3w7LPP0q9fP/Lz83nhhReYN28ekydP5qqrrsLdGTZsGIcccgglJSXcfvvt5WK/7777gKA7iVGjRnHuueeWddZXUlJS9plGjx7NsmXLKCkp4ZNPPuHtt99m4cKFZd/LrlAiEMlhe+21F3PnzmX48OE0a9aM008/nREjRmBmnH322Tz11FOsW7eOGTNm0Lt3bwDq169ftqd+2GGHpe1+uvQoIrFcuq6tK/LKK69wySWXlL1v3LhxufHvvvsuBx10EG3atAHg3HPPLfccg9IeQ8844wxmzJgBwJFHHsnAgQN5+OGHy5p7Kut2+6CDDgKgU6dOfPrpp6xYsYL58+fTuHFjDjzwQNyd6667jqKiIn70ox/xySeflHWxXZHErqkPPfRQWrVqxXvvvQfAscceS6NGjcjPz6dt27Z8/PHHHHzwwXz00Udcdtll/POf/0zbaV1UOlkskuPq1q1LcXExxcXFFBYW8vjjjzNw4EAGDRrESSedRH5+Pv369Svrxz8vL6+s6+TE7qdTKW1ySiyXrmvrirh72TIrGp9O4rSlrx988EFmzpzJiy++SMeOHSkpKUkbW2LX1AB9+/ZlzJgx/Pe//y176tnTTz/N6tWrmTt3Lnl5eRQUFKTtiruy2FN1Td24cWPmz5/PxIkTue+++3j22Wd59NFH0y6jMrEdEZjZAWY22cyWmNkiM/tVijLFZrbezErCv+rvaFtEKvTuu++W2xsvKSmhVatWQHBCs0WLFtx8880MHDiw0nml6345UUVdW6ebvlevXtx7771l7z///PNy4w899FCWLl1a1v7/5JNPlnVNDZQ144wePbqsZ9MPP/yQ7t27c+ONN9K0aVOWLVtWpW63+/fvzzPPPMOYMWPo27cvEHRNvd9++5GXl8fkyZP5+OOPK103iV1Tv/fee/znP//5xjMSEq1Zs4YdO3bw05/+lJtuuqlauqaO84hgG3CVu88zs4bAXDN72d2TLyuY5u4nxhiHyO4jwuWe1Wnjxo1cdtllrFu3jnr16vHd736X4cOHl40/66yzWL16ddkloen079+fCy64gLvvvjvtFTEVdW2drvvm66+/nksuuYT27dtTt25dhgwZUtbsBMHDYx577DH69evHtm3b6Nq1a7nnMG/ZsoXu3buzY8cORo0aBcDgwYN5//33cXeOPfZYOnToQFFRUeRut9u1a8eGDRvYf//9ad68edn6Oumkk+jSpQsdO3bk0EMPBaBJkyYceeSRtG/fnt69e5dr5rr44ou56KKLKCwspF69eowYMaLckUCyTz75hEGDBpU9/Oe2226rsGxUGeuG2sz+Dtzr7i8nDCsGflOVRLAr3VBX+zXaVVq27iyWb6rp3VBfeumldOrUifPPPz/boUgV1MhuqM2sAOgEzEwx+nAzm29mL5lZuwqmv9DM5pjZnNKrGEQkXocddhgLFixgwIAB2Q5FYhb7yWIz2wt4HrjC3b9IGj0PaOXuG82sDzAWaJ08D3cfDgyH4Igg3ohFBCh73KTUfrEeEZhZHkESeNrd/5Y83t2/cPeN4esJQJ6ZZe8+cJEs2N2eEig1287UpzivGjLgEWCJu/+5gjLfCcthZt3CeNbGFZNITZOfn8/atWuVDKRauDtr164lPz+/StPF2TR0JHA2sNDMSsJh1wEHArj7g0Bf4Jdmtg34Cujv2iIkh7Rs2ZLly5ejc19SXfLz82nZsmWVpoktEbj7G0DFd4AEZe4F7k1XRqQ2y8vLK7tbVSRb1MWEiEiOUyIQEclxSgQiIjlOiUBEJMcpEYiI5DglAhGRHKdEICKS45QIRERynBKBiEiOUyIQEclxSgQiIjmu0kRgZoeY2R7h62Izu9zM9ok9MhERyYgoRwTPA9vN7LsE3UofBIyMNSoREcmYKIlgh7tvA34C3OXuVwLN4w1LREQyJUoi2GpmZwDnAuPDYXnxhSQiIpkUJREMAg4HbnH3f5vZQcBT8YYlIiKZUumDadx9MXB5wvt/A8PiDEpERDKn0kRgZkcCQ4FWYXkD3N0Pjjc0ERHJhCiPqnwEuBKYC2yPNxwREcm0KIlgvbu/FHskIiKSFVESwWQzux34G7CldKC7z4stKhERyZgoiaB7+L9LwjAHelZ/OCIikmlRrho6JhOBiIhIdkTpa6iRmf3ZzOaEf3eYWaNMBCciIvGLckPZo8AG4Gfh3xfAY3EGJSIimRPlHMEh7v7ThPc3mFlJTPGIiNRsQ7PYIDJ0fSyzjXJE8JWZ/aD0TXiD2VexRCMiIhkX5Yjgl8Dj4XkBAz4DBsYZlIiIZE6Uq4ZKgA5mtnf4/osoMzazA4AngO8AO4Dh7v6XpDIG/AXoA3wJDNT9CSIimVVhIjCzAe7+lJn9Omk4AO7+50rmvQ24yt3nmVlDYK6ZvRx2YleqN9A6/OsOPMDX9y2IiEgGpDsiaBD+b5hinFc2Y3dfCawMX28wsyXA/kBiIjgFeMLdHXjTzPYxs+bhtCIikgEVJgJ3fyh8+Yq7/ytxXHjCODIzKwA6ATOTRu0PLEt4vzwcVi4RmNmFwIUABx54YFUWLSIilYhy1dA9EYelZGZ7ETz3+IoU5xcsxSTfONpw9+Hu3sXduzRr1izqokVEJIJ05wgOB44AmiWdJ9gbqBtl5maWR5AEnnb3v6Uoshw4IOF9S2BFlHmLiEj1SHdEUB/YiyBZNEz4+wLoW9mMwyuCHgGWpDmxPA44xwI9CLq81vkBEZEMSneOYCow1cxGuPvHOzHvI4GzgYUJdyJfBxwYzv9BYALBpaMfEFw+OmgnliMiIrsgyg1lX4bPI2gH5JcOdPe03VC7+xukPgeQWMaBSyLEICIiMYlysvhp4B3gIOAGYCkwO8aYREQkg6Ikgibu/giw1d2nuvt5QI+Y4xIRkQyJ0jS0Nfy/0sxOILiqp2V8IYmISCZFSQQ3hx3OXUVw/8DewJWxRiUiIhkTJRHMdPf1wHpAj60UEallopwjmG5mk8zsfDNrHHtEIiKSUZUmAndvDVxPcPnoXDMbb2YDYo9MREQyIsoRAe4+y91/DXQjeDDN47FGJSIiGVNpIjCzvc3sXDN7CZhO0DNot9gjExGRjIhysng+MBa40d1nxBuOiIhkWtpEYGZ1gRfCZiEREamF0jYNuft2oEOGYhERkSyI0jRUYmbjgOeATaUDK3i+gIiI7GaiJIJ9gbVAYm+jDigRiIjUApUmAnfXMwJERGqxKJePtjGzV83s7fB9kZldH39oIiKSCVFuKHsYuJawF1J3XwD0jzMoERHJnCiJYE93n5U0bFscwYiISOZFSQRrzOwQghPEmFlfgruLRUSkFohy1dAlwHDgUDP7BPg3oE7nRERqiShXDX0E/MjMGgB13H1D/GGJiEimRLlq6FdmtjfwJXCnmc0zs17xhyYiIpkQ5RzBee7+BdAL2A8YBAyLNSoREcmYKInAwv99gMfcfX7CMBER2c1FSQRzzWwSQSKYaGYNgR3xhiUiIpkS5aqh84GOwEfu/qWZNSFoHhIRkVogylVDO8ysABhgZg684e4vxB6ZiIhkRJSrhu4HLgIWAm8DvzCz+yJM96iZfVraR1GK8cVmtt7MSsK/31c1eBER2XVRmoaOBtq7e+mdxY8TJIXKjADuBZ5IU2aau58YYV4iIhKTKCeL3wUOTHh/ALCgsonc/XXgs52MS0REMqTCIwIz+wdB/0KNgCVmVtrxXDdgejUt/3Azmw+sAH7j7osqiOVC4EKAAw88MFURERHZSemahv4U87LnAa3cfaOZ9QHGAq1TFXT34QT9HdGlSxePOS4RkZxSYSJw96mlr83s20DX8O0sd/90Vxcc3q1c+nqCmd1vZk3dfc2uzltERKKLctXQz4BZQD/gZ8DMsCvqXWJm3zEzC193C2NZu6vzFRGRqoly1dD/AV1LjwLMrBnwCjAm3URmNgooBpqa2XJgCJAH4O4PAn2BX5rZNuAroH/plUkiIpI5URJBnaSmoLVEOJJw9zMqGX8vweWlIiKSRVESwT/NbCIwKnx/OjAhvpBERCSTonQxMdjMTgN+QNDr6HB1MSEiUntEOSLA3f8G/C3mWEREJAui3FksIiK1mBKBiEiOqzARmNmr4f8/ZC4cERHJtHTnCJqb2dHAyWb2DEmPp3T3ebFGJiIiGZEuEfweuAZoCfw5aZwDPeMKSkREMiddX0NjgDFm9jt3vymDMYmISAZFuY/gJjM7GTgqHDTF3cfHG5aIiGRKlE7nbgN+BSwO/34VDhMRkVogyg1lJwAd3X0HlD2q8i3g2jgDExGRzIh6H8E+Ca8bxRCHiIhkSZQjgtuAt8xsMsElpEehowERkVojysniUWY2heAJZQZc7e7/jTswERHJjKidzq0ExsUci4iIZIH6GhIRyXFKBCIiOS5tIjCzOmb2dqaCERGRzEubCMJ7B+ab2YEZikdERDIsysni5sAiM5sFbCod6O4nxxaViIhkTJREcEPsUYiISNZEuY9gqpm1Alq7+ytmtidQN/7QREQkE6J0OncBMAZ4KBy0PzA2xphERCSDolw+eglwJPAFgLu/D+wXZ1AiIpI5URLBFnf/X+kbM6tH8IQyERGpBaIkgqlmdh3wLTP7MfAc8I94wxIRkUyJkgiuAVYDC4FfABOA6+MMSkREMifKVUM7wofRzCRoEnrX3SttGjKzR4ETgU/dvX2K8Qb8BegDfAkMdPd5VYxfRER2UZSrhk4APgTuBu4FPjCz3hHmPQI4Ps343kDr8O9C4IEI8xQRkWoW5YayO4Bj3P0DADM7BHgReCndRO7+upkVpClyCvBEeHTxppntY2bNwy6vRUQkQ6KcI/i0NAmEPgI+rYZl7w8sS3i/PBz2DWZ2oZnNMbM5q1evroZFi4hIqQqPCMzstPDlIjObADxLcI6gHzC7GpZtKYalPPfg7sOB4QBdunTRpasiItUoXdPQSQmvVwFHh69XA42rYdnLgQMS3rcEVlTDfEVEpAoqTATuPijmZY8DLjWzZ4DuwHqdHxARybxKTxab2UHAZUBBYvnKuqE2s1FAMdDUzJYDQ4C8cNoHCe5H6AN8QHD5aNyJR0REUohy1dBY4BGCu4l3RJ2xu59RyXgn6MdIRESyKEoi2Ozud8ceiYiIZEWURPAXMxsCTAK2lA7UXcAiIrVDlERQCJwN9OTrpiEP34uIyG4uSiL4CXBwYlfUIiJSe0S5s3g+sE/McYiISJZEOSL4NvCOmc2m/DmCtJePiojI7iFKIhgSexQiIpI1UZ5HMDUTgYiISHZEubN4A193Blef4O7gTe6+d5yBiYhIZkQ5ImiY+N7MTgW6xRWQiIhkVpSrhspx97HoHgIRkVojStPQaQlv6wBdqOC5ASIisvuJctVQ4nMJtgFLCR4zKSIitUCUcwTqHlpEpBZL96jK36eZzt39phjiERGRDEt3RLApxbAGwPlAE0CJQESkFkj3qMo7Sl+bWUPgVwRPEXsGuKOi6UREZPeS9hyBme0L/Bo4C3gc6Ozun2ciMBERyYx05whuB04DhgOF7r4xY1GJiEjGpLuh7CqgBXA9sMLMvgj/NpjZF5kJT0RE4pbuHEGV7zoWEZHdj37sRURynBKBiEiOUyIQEclxSgQiIjlOiUBEJMcpEYiI5DglAhGRHBdrIjCz483sXTP7wMyuSTG+2MzWm1lJ+Jeux1MREYlBlAfT7BQzqwvcB/wYWA7MNrNx7r44qeg0dz8xrjhERCS9OI8IugEfuPtH7v4/gl5L9WQzEZEaJs5EsD+wLOH98nBYssPNbL6ZvWRm7VLNyMwuNLM5ZjZn9erVccQqIpKz4kwElmJY8kPv5wGt3L0DcA8wNtWM3H24u3dx9y7NmjWr3ihFRHJcnIlgOXBAwvuWwIrEAu7+RWn31u4+Acgzs6YxxiQiIkniTASzgdZmdpCZ1Qf6A+MSC5jZd8zMwtfdwnjWxhiTiIgkie2qIXffZmaXAhOBusCj7r7IzC4Kxz8I9AV+aWbbgK+A/u6e3HwkIiIxii0RQFlzz4SkYQ8mvL4XuDfOGEREJD3dWSwikuOUCEREcpwSgYhIjlMiEBHJcUoEIiI5TolARCTHKRGIiOQ4JQIRkRynRCAikuOUCEREcpwSgYhIjou1ryGRnDO0UZaWuz47y5VaQUcEIiI5TolARCTHKRGIiOQ4JQIRkRynRCAikuOUCEREcpwSgYhIjlMiEBHJcUoEIiI5TolARCTHKRGIiOQ4JQIRkRynRCAikuOUCEREcpwSgYhIjlMiEBHJcbEmAjM73szeNbMPzOyaFOPNzO4Oxy8ws85xxiMiIt8UWyIws7rAfUBvoC1whpm1TSrWG2gd/l0IPBBXPCIiklqcRwTdgA/c/SN3/x/wDHBKUplTgCc88Cawj5k1jzEmERFJEuczi/cHliW8Xw50j1Bmf2BlYiEzu5DgiAFgo5m9u5MxNQXW7OS0u+YGSzc2e3FVrqbGprgSpa9foPVVVTUzrhtsV+JqVdGIOBNBqprpO1EGdx8ODN/lgMzmuHuXXZ1PdaupcUHNjU1xVY3iqppciyvOpqHlwAEJ71sCK3aijIiIxCjORDAbaG1mB5lZfaA/MC6pzDjgnPDqoR7AendfmTwjERGJT2xNQ+6+zcwuBSYCdYFH3X2RmV0Ujn8QmAD0AT4AvgQGxRVPaJebl2JSU+OCmhub4qoaxVU1ORWXuX+jSV5ERHKI7iwWEclxSgQiIjmu1iSCXenOorJpY47rrDCeBWY23cw6JIxbamYLzazEzOZkOK5iM1sfLrvEzH4fddqY4xqcENPbZrbdzPYNx8W5vh41s0/N7O0KxmerflUWV7bqV2VxZat+VRZXxuuXmR1gZpPNbImZLTKzX6UoE2/9cvfd/o/gZPSHwMFAfWA+0DapTB/gJYJ7F3oAM6NOG3NcRwCNw9e9S+MK3y8FmmZpfRUD43dm2jjjSip/EvBa3OsrnPdRQGfg7QrGZ7x+RYwr4/UrYlwZr19R4spG/QKaA53D1w2B9zL9+1Vbjgh2pTuLKNPGFpe7T3f3z8O3bxLcSxG3XfnMWV1fSc4ARlXTstNy99eBz9IUyUb9qjSuLNWvKOurIlldX0kyUr/cfaW7zwtfbwCWEPSwkCjW+lVbEkFFXVVEKRNl2jjjSnQ+QdYv5cAkM5trQTcb1SVqXIeb2Xwze8nM2lVx2jjjwsz2BI4Hnk8YHNf6iiIb9auqMlW/osp0/YosW/XLzAqATsDMpFGx1q84u5jIpF3pziJSNxc7KfK8zewYgg31BwmDj3T3FWa2H/Cymb0T7tFkIq55QCt332hmfYCxBL3E1oj1RXDY/i93T9y7i2t9RZGN+hVZhutXFNmoX1WR8fplZnsRJJ4r3P2L5NEpJqm2+lVbjgh2pTuLOLu5iDRvMysC/gqc4u5rS4e7+4rw/6fACwSHgRmJy92/cPeN4esJQJ6ZNY0ybZxxJehP0mF7jOsrimzUr0iyUL8qlaX6VRUZrV9mlkeQBJ5297+lKBJv/aruEx/Z+CM4svkIOIivT5i0SypzAuVPtsyKOm3McR1IcGf1EUnDGwANE15PB47PYFzf4esbDrsB/wnXXVbXV1iuEUE7b4NMrK+EZRRQ8cnPjNeviHFlvH5FjCvj9StKXNmoX+HnfgK4K02ZWOtXrWga8l3ozqKiaTMY1++BJsD9ZgawzYPeBb8NvBAOqweMdPd/ZjCuvsAvzWwb8BXQ34Oal+31BfATYJK7b0qYPLb1BWBmowiudGlqZsuBIUBeQlwZr18R48p4/YoYV8brV8S4IPP160jgbGChmZWEw64jSOIZqV/qYkJEJMfVlnMEIiKyk5QIRERynBKBiEiOUyIQEclxSgQiIjlOiUAkAjPbx8wuDl8Xm9n4Kk4/wsz6xhOdyK5RIhCJZh/g4mwHIRKHWnFDmUgGDAMOCW/42QpsMrMxQHtgLjDA3T3sV/8k4FsEd5/+wnWzjtRwOiIQieYa4EN37wgMJugh8gqgLUFf8EeG5e51967u3p4gGZyY+VBFqkaJQGTnzHL35e6+Aygh6L8G4Bgzm2lmC4GeQLsKphepMdQ0JLJztiS83g7UM7N84H6gi7svM7OhQH42ghOpCh0RiESzgeAxgumU/uivCfuW11VCslvQEYFIBO6+1sz+FT70/CtgVYoy68zsYWAhwfNtZ2c2SpGdo95HRURynJqGRERynBKBiEiOUyIQEclxSgQiIjlOiUBEJMcpEYiI5DglAhGRHPf/VAxkzS6JyVkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA64ElEQVR4nO3deXwU9fnA8c+TAwLIFS6BAAEEEUgIEA5Pgli5PKiCokhBvEWtR1Hpz1ZbtVqttVqr1VqRggfggdQTREEUBYECitwQIIAQCGe4k+f3x0zW3ZhjArO7JHner9e+nt05vvPs7Ox+d74z8x1RVYwxxhiAmGgnYIwx5uRhlYIxxpgAqxSMMcYEWKVgjDEmwCoFY4wxAVYpGGOMCbBK4SQhIr8VkZfLS7nHS0QaicgXIrJPRJ6Kdj4FRGSYiEyPdh6lEZFZInJ9lHN4VUQe8bG8ZSKS4Vd5bpmZInKBn2X65WT7ThZWKSoFEakqIv8WkQ3uj9H/RKR/0PgMEckXkf3uI0tEJotIt0jlqKp/UtUT+rK77yPL73J9diOwA6ilqvdEO5kCqvqaql4Y7TxONiIyUkS+9LG8n1UoqtpBVWf5tYyTSTn5ToaoFJUCEAdsAnoBtYHfAZNFJDlomi2qegpQE+gJrADmiEgfPxIQkTg/yqkAWgA/qF01acqBSvm9VdVK+QCWApe7zzOArCKmeQ5YUMz8yYDi/PPdAmwF7gka/xDwFjAR2AtcDzQBpgE5wBrghkLTTwx63ROYC+wGlgAZQeMSgXHucncBU4EawEEgH9jvPpoUUe4lwDK33FnAGUHjMoHfuOtmDzAJSHDH1Qfed+fLAeYAMcWsm7OAb90yvgXOcoe/ChwFjrj5XVDEvPWA/7rr7FvgEeDLoPHtgBluDiuBK4LGvQr8A/gA2AfMA1oX+rzigqafBVzvPh9ZaDkK3AysdtfxPwAJGj8KWO6O+wRoUcK2NgX40V0fXwAdvOTsjv8Fzh+UPTjb4+yCnItYTndggbvutgF/dYd/ANxexPY/qKT3CpwBHALy3M9rt8eci/yMcL4rwZ//f4O2uwvc57HAb4G1btkLgWbFvN/StuWxwA/uexqHuy274y8CFrvzzgVSC817n7uODgMPAG8VWvYzwLPu82vdbWEfsA64yR1+0nwny/Tb6OcPbXl5AI3cjb2d+zqDoiuF890PtEYR45JxvkxvuB9+CpAdtHE/5H4BBuHskVXD+UI/DyQAae70fYKmn+g+bwrsBAa48/7Cfd0g6Es+CagLxAO9insfhcptC+S65cUD9+JUTlWCNsD57oab6G7oN7vjHgP+6c4XD5xL0I9k0PIScb6Ew3H20K5yX9fTn35QHinhs3nTfVQH2uPs4X0Z9CXbhPMljAO64DRFdQgqOwfnxzEOeA14s9DnVZZK4X2gDtDc/az6ueMGuevtDHc5DwBzS3hPo3D2QKsCfwMWB40rKef6OD/wg911fhdwjOIrha+B4e7zU4Ce7vMrgHlB03XC2Z6qeHivIevFQ85ePqNHCpWXyU/fmzHAd8DpOBVTJ9xtp9A8Xrbl74FmONvkVwXLdXPaDvTAqYRGuNNXDZp3sTtvNZy92wM4TZ6482wNWr8DgdZuvr3cabucLN/Jsj4qS/NRgIjE42zE41V1RSmTb8H5oOuUMM0fVDVXVb/D+TdyVdC4r1V1qqrm43zBzwHuU9VDqroYeBnnx7Owa4APVfVDVc1X1Rk4/wAHiEhjoD/OhrFLVY+q6uzS3rfrSuADVZ2hqkeBv+Bs9GcFTfOsqm5R1Rycf+xp7vCjQGOcf8RHVXWOultmIQOB1ao6QVWPqeobOP90Ly4tORGJBS4HHlTVA6r6AzA+aJKLgExVHeeWvQh4G+dHs8A7qjpfVY/hfM5pHL/HVXW3qm4EPg8q6ybgMVVd7i7nT0CaiLQoqhBVfUVV96nqYZwfhE4iUttDzgNwmtrecj+vv+HscRTnKHCaiNRX1f2q+o07/D2gjYi0cV8PByap6hEP77U4xeXs5TMqyfXAA6q6Uh1LVHVnEdN52ZafU9VN7rb8KD99N28AXlTVeaqap6rjcfYIegbN+6w770FV3QAswvkzAM6fxQMF61dVP1DVtW6+s4HpOD/QXkTiO1kmlapSEJEYYALO7uttHmZpivMvancJ02wKer4Bp0YvalwTIEdV9xWavmkRZbYAhojI7oIHToXSGOffS46q7vKQf2FN3GUC4FZWmwrlEPyjcwDnHyfAkzj/YKaLyDoRud/LMlzFvc/CGvDT8Z8Cwc9bAD0KrZdhwKke8j8exZXVAngmKIccnD8PP3uPIhIrIo+LyFoR2Yvzzw+cPwmlLacJQe/f/cIHr4/CrsP557lCRL4VkYvc+Q4Dk4Fr3O/AVTjfAy/vtTglrZvSPqOSNMNpOiqNl225uO9mC+CeQjk2o/jvLsDr/FSpXO2+BkBE+ovINyKS45Y1gNDP90Tfx4l+J8uk0lQKIiLAv3Gaji53a+XS/BJYpKq5JUzTLOh5c5y9iwLBtfYWIFFEahaafnMRZW4CJqhqnaBHDVV93B2XKCJ1ipivtH8JW3C+EEBgnTQrJofQgp1/uveoaiucf/13F3MQPmQZruLeZ2HZOM0jSUHDgtfvJmB2ofVyiqre4qHsgs+wetAwrz9UhW3CaTcOzqOaqs4tYtqrgUuBC3BOckh2h4uH5Wwl6P0HfV5FUtXVqnoV0BD4M/CWiNRwR4/H+XHug/Mv92sPy4fSt6nCSvuMSitvE05TTGm8bMvFfTc3AY8WyrG6u1dboHCeU4AMEUnC+V143V1uVZw9ob8AjVS1DvAhP32+J8N3skwqTaUAvIDTBnyxqh4sbiJxNBWRB3F2ZX9bSrm/E5HqItIBpx11UlETqeomnANaj4lIgoik4vyze62IyScCF4tIX/efZoJ7aluSqm4FPgKeF5G6IhIvIue5820D6hVqmgg2GRgoIn3cZrR7cHabi/oxCyEiF4nIae5Guxfn4GNeEZN+CLQVkatFJE5ErsQ5NvB+actQ1TzgHeAhd522A34VNMn7btnD3fcdLyLdROQMD2Vn43zRrnHX6Si8/fgU5Z/AWPczR0Rqi8iQYqatibOOd+JUSH8qw3I+ADqIyGXuWTB3UEJFJiLXiEgD99/mbndwHoBbCeQDT/HzvYSSbAOSRKSKx+lL+4y2Aa1KmP9l4GERaeN+F1NFpF4R03nZlkeLSJKIJOJ8jwu+m/8CbhaRHu4yaojIwEJ/2EK4288snCbi9aq63B1VBedYUTZwTJxT3YNPbT4ZvpNlUikqBXHaem/CaYv7UX66HmFY0GRNRKTgDIFvcQ4cZ6hqaRc0zcbZhZsJ/KWU6a/C+ae4BXgXp+18RuGJ3ArkUpwNORvnn80Yfvq8huO0J67AOWB2pzvfCpwD3+vc3eImhcpdiXO84u84B/8uxqkkg9uWi9MG+BRn/XwNPK9FnFvutv9ehLNx78Q5cHaRqu7wsAxwmvVq4+wyT3Dfz2G37H04X7ihOOvwR5x/xFU9ln0DznrcCXTAwxevKKr6rrvcN90moe9xjvMU5T84zQObcc6E+aaY6Ypazg5gCPC4m3MbnAOmxekHLHO342eAoap6qFAuKTh/Orz6DOfMmB9FpNTP0MNn9G+gvbt9Ti2iiL/i/FBOx/mh+zdOG3vh5XjZll93y1nnPh5x512Asy08h3MSxBqcA+qleR1njy/QdOS+3zvcnHfh7BlOCxof9e9kWYkPxyUqJXGucVgPxLsH2060vD8CSao66kTLqkhE5M/Aqao6Itq5lHci8ivgRlU9J9q5mJNXpdhTONm5u3/tcSqZSk1E2rlNBiIi3XGa2N6Ndl7lnYhUB24FXop2LubkZpXCyWERzsHVf0U7kZNATZzjCrk4u+RP4ZxSaY6TiPTFaYbcRlDThzFFseYjY4wxAbanYIwxJqBcd/ZUv359TU5OjnYaxhhTrixcuHCHqjYoaly5rhSSk5NZsGBBtNMwxphyRUQK9zoQYM1HxhhjAqxSMMYYE2CVgjHGmIByfUyhKEePHiUrK4tDhw6VPrExHiUkJJCUlER8fHy0UzEmrCpcpZCVlUXNmjVJTk7GuVDYmBOjquzcuZOsrCxatmwZ7XSMCasK13x06NAh6tWrZxWC8Y2IUK9ePdv7NJVChasUAKsQjO9smzKVRYWsFIwxpiLLyT3Ci7PXkpPrpYftsrFKIQy2bdvG1VdfTatWrejatStnnnkm777rf0ef06ZN4/HHH/e93GCnnHIid7P0R2ZmJq+//lM/bgsWLOCOO+6IYkbGRNeUBZt47KMVTFlQ0t1Zj0+FO9AcbarKoEGDGDFiROCHbMOGDUybNq2UOcvukksu4ZJLLvG93HA6duwYcXFl2+wKKoWrr74agPT0dNLT08ORnjHlwpD0ZiHRT7an4LPPPvuMKlWqcPPNNweGtWjRgttvvx2AvLw8xowZQ7du3UhNTeXFF18EYNasWWRkZDB48GDatWvHsGHDKOjBNjk5mR07nJteLViwgIyMDABeffVVbrvtNgBGjhzJHXfcwVlnnUWrVq146623Ast/8sknA8t78MEHi8z7jTfeICUlhY4dO3LfffeFjLvnnnvo0qULffr0ITs7G4Bnn32W9u3bk5qaytChQwHIzc1l1KhRdOvWjc6dO/Pee+8F8hwyZAgXX3wxF154IVdeeSUffvhhoPyRI0fy9ttvk5mZybnnnkuXLl3o0qULc+c6N0a7//77mTNnDmlpaTz99NPMmjWLiy66CICcnBwGDRpEamoqPXv2ZOnSpQA89NBDjBo1ioyMDFq1asWzzz4byHHgwIF06tSJjh07MmlSkXdPNeakllijCjf1ak1iDa93SS0DVS23j65du2phP/zww8+GRdIzzzyjd955Z7HjX3zxRX344YdVVfXQoUPatWtXXbdunX7++edaq1Yt3bRpk+bl5WnPnj11zpw5qqraokULzc7OVlXVb7/9Vnv16qWqquPGjdPRo0erquqIESN08ODBmpeXp8uWLdPWrVurquonn3yiN9xwg+bn52teXp4OHDhQZ8+eHZLT5s2btVmzZrp9+3Y9evSo9u7dW999911VVQV04sSJqqr6hz/8IbC8xo0b66FDh1RVddeuXaqqOnbsWJ0wYUJgWJs2bXT//v06btw4bdq0qe7cuVNVVd955x391a9+paqqhw8f1qSkJD1w4IDm5ubqwYMHVVV11apVWvD5fv755zpw4MBAvsGvb7vtNn3ooYdUVXXmzJnaqVMnVVV98MEH9cwzz9RDhw5pdna2JiYm6pEjR/Stt97S66+/PlDW7t27i/2sCov2tmWMX4AFWszvqu0pEN6DNqNHj6ZTp05069YNgOnTp/Of//yHtLQ0evTowc6dO1m9ejUA3bt3JykpiZiYGNLS0sjMzCzTsgYNGkRMTAzt27dn27ZtgeVNnz6dzp0706VLF1asWBFYXoFvv/2WjIwMGjRoQFxcHMOGDeOLL74AICYmhiuvvBKAa665hi+//BKA1NRUhg0bxsSJEwPNQdOnT+fxxx8nLS2NjIwMDh06xMaNGwH4xS9+QWJiIgD9+/fns88+4/Dhw3z00Uecd955VKtWjaNHj3LDDTeQkpLCkCFD+OGHH0p9z19++SXDhw8H4Pzzz2fnzp3s2bMHgIEDB1K1alXq169Pw4YN2bZtGykpKXz66afcd999zJkzh9q1i7ufujGVkx1T4KeDNgA39Wp9QmV16NCBt99+O/D6H//4Bzt27Ai0gasqf//73+nbt2/IfLNmzaJq1Z/uPx8bG8uxY86tn+Pi4sjPzwco8Vz54PnVbXpSVcaOHctNN91U7HwF03pRcGrmBx98wBdffMG0adN4+OGHWbZsGarK22+/zemnnx4yz7x586hRo0bgdUJCAhkZGXzyySdMmjSJq666CoCnn36aRo0asWTJEvLz80lISCg1n6JyL8ixqPXZtm1bFi5cyIcffsjYsWO58MIL+f3vf+/5/RtT0dmeAs7BmrH92/ly0Ob888/n0KFDvPDCC4FhBw4cCDzv27cvL7zwAkePHgVg1apV5ObmllhmcnIyCxcuBAipcLzo27cvr7zyCvv37wdg8+bNbN++PWSaHj16MHv2bHbs2EFeXh5vvPEGvXr1AiA/Pz9wfOL111/nnHPOIT8/n02bNtG7d2+eeOIJdu/ezf79++nbty9///vfAz/U//vf/4rNa+jQoYwbN445c+YEKsg9e/bQuHFjYmJimDBhAnl5eQDUrFmTffv2FVnOeeedx2uvvQY4FWv9+vWpVatWscvdsmUL1atX55prruE3v/kNixYtKnUdGlOZ2J4CPx208YOIMHXqVO666y6eeOIJGjRoQI0aNfjzn/8MwPXXX09mZiZdunRBVWnQoAFTp04tscwHH3yQ6667jj/96U/06NGjTPlceOGFLF++nDPPPBNwTjGdOHEiDRs2DEzTuHFjHnvsMXr37o2qMmDAAC699FIAatSowbJly+jatSu1a9dm0qRJ5OXlcc0117Bnzx5Ulbvuuos6derwu9/9jjvvvJPU1FRUleTkZN5///1i8/rVr37FJZdcQpUqzsGyW2+9lcsvv5wpU6bQu3fvwN5FamoqcXFxdOrUiZEjR9K5c+dAOQ899BDXXnstqampVK9enfHjx5e4Pr777jvGjBlDTEwM8fHxIZW3Maac36M5PT1dC99kZ/ny5ZxxxhlRyshUZLZtmYpCRBaqapHndVvzkTHGmACrFIwxxgRYpWCMMSbAKgVjjDEBVikYY4wJsErBGGNMgFUKYfDoo4/SoUMHUlNTSUtLY968ecdVztSpU0O6esjIyKDwKbglCVeX08Ed9EXL7t27ef755wOvt2zZwuDBg6OYkTEVg1UKPvv66695//33WbRoEUuXLuXTTz+lWbPju1K6cKVQVoUrhfT09EBvoSeTgiuXy6JwpdCkSZOQnmGNMcfHKgWfbd26lfr16wf63alfvz5NmjRh5syZ/PKXvwxMN2PGDC677DLAucr4//7v/+jUqRM9e/Zk27ZtzJ07l2nTpjFmzBjS0tJYu3YtAFOmTKF79+60bduWOXPmAMV3x11Sl9P79+/n2muvJSUlhdTU1CK7z5g5cyadO3cmJSWFUaNGcfjw4cC4J598ku7du9O9e3fWrFkTyK1jx4506tSJ8847r8TcZs2aRe/evbn66qtJSUnhvvvuC/mRf+ihh3jqqafYv38/ffr0oUuXLqSkpAS6477//vtZu3YtaWlpjBkzhszMTDp27Ag4/UMVvLfOnTvz+eefA04X3pdddhn9+vWjTZs23HvvvYEcR44cSceOHUlJSeHpp58+zk/fmAqguO5Ty8PjZOw6e9++fdqpUydt06aN3nLLLTpr1ixVVc3Pz9fTTz9dt2/frqqqV111lU6bNk1Vne6pC56PGTMm0LX2iBEjdMqUKYGye/XqpXfffbeqqn7wwQfap08fVS25O+7iupy+99579de//nVgXE5OTsj7OHjwoCYlJenKlStVVXX48OH69NNPq6rTlfcjjzyiqqrjx48PlNmxY0fNyspS1Z+60y4pt+rVq+u6detUVXXRokV63nnnBZZ/xhln6IYNG/To0aO6Z88eVVXNzs7W1q1ba35+vq5fv147dOgQmD749V/+8hcdOXKkqqouX75cmzVrpgcPHtRx48Zpy5Ytdffu3Xrw4EFt3ry5bty4URcsWKAXXHBBoKyC3AuL9rZljF+wrrNLcewIrPzIiSfolFNOYeHChbz00ks0aNCAK6+8kldffRURYfjw4UycOJHdu3fz9ddf079/fwCqVKkS+AfftWvXErvMLti7CJ6upO64i/Ppp58yevTowOu6deuGjF+5ciUtW7akbdu2AIwYMSLQnTYQ6Nn0qquu4uuvvwbg7LPPZuTIkfzrX/8KNAmV1lV4y5YtAejcuTPbt29ny5YtLFmyhLp169K8eXNUld/+9rekpqZywQUXsHnz5kC34MUJ7k67Xbt2tGjRglWrVgHQp08fateuTUJCAu3bt2fDhg20atWKdevWcfvtt/Pxxx+X2KGeMRWddYgHsHYmTBoOV06A0/ufcHGxsbFkZGSQkZFBSkoK48ePZ+TIkVx77bVcfPHFJCQkMGTIkMB9COLj4wPdPQd3mV2Ugmap4Om0hO64i6OqgWUWN74kwfMWPP/nP//JvHnz+OCDD0hLS2Px4sUl5hbcnTbA4MGDeeutt/jxxx8Dd3N77bXXyM7OZuHChcTHx5OcnFxi9+Gl5V5Ud9p169ZlyZIlfPLJJ/zjH/9g8uTJvPLKKyUuw5iKyvYUAFr3cSqE1n1OuKiVK1eG/EtfvHgxLVq0AJyDoU2aNOGRRx5h5MiRpZZVUpfRwYrrjruk+S+88EKee+65wOtdu3aFjG/Xrh2ZmZmB4wUTJkwIdKcNBG5jOWnSpEAPrGvXrqVHjx788Y9/pH79+mzatKlMXYUPHTqUN998k7feeitwJtGePXto2LAh8fHxfP7552zYsKHUdRPcnfaqVavYuHHjz+7xEGzHjh3k5+dz+eWX8/DDD1t32qZSsz0FgLgqvuwhgHMA9/bbb2f37t3ExcVx2mmn8dJLLwXGDxs2jOzsbNq3b19qWUOHDuWGG27g2WefLfHMmuK64y6py+kHHniA0aNH07FjR2JjY3nwwQcDTVPg3Ahn3LhxDBkyhGPHjtGtW7eQ+04fPnyYHj16kJ+fzxtvvAHAmDFjWL16NapKnz596NSpE6mpqZ67Cu/QoQP79u2jadOmNG7cOLC+Lr74YtLT00lLS6Ndu3YA1KtXj7PPPpuOHTvSv3//kKawW2+9lZtvvpmUlBTi4uJ49dVXQ/YQCtu8eTPXXntt4EZGjz32WLHTGlPRWdfZEXbbbbfRuXNnrrvuuminYsroZN+2jPGqpK6zbU8hgrp27UqNGjV46qmnop2KMcYUKayVgohkAvuAPOCYqqaLSCIwCUgGMoErVHWXO/1Y4Dp3+jtU9ZNw5hdpBbfUNMaYk1UkDjT3VtW0oF2V+4GZqtoGmOm+RkTaA0OBDkA/4HkRiT2eBZbnJjFzcrJtylQW0Tj76FKg4Ea644FBQcPfVNXDqroeWAN0L2vhCQkJ7Ny5077Exjeqys6dO0lISIh2KsaEXbiPKSgwXUQUeFFVXwIaqepWAFXdKiIFd5BvCnwTNG+WOyyEiNwI3AjQvHnzny0wKSmJrKwssrOzfX0jpnJLSEggKSkp2mkYE3bhrhTOVtUt7g//DBFZUcK0RV1J9bO/+27F8hI4Zx8VHh8fHx+4StYYY0zZhLX5SFW3uHE78C5Oc9A2EWkM4Mbt7uRZQHB3oknAlnDmZ4wxJlTYKgURqSEiNQueAxcC3wPTgBHuZCOA99zn04ChIlJVRFoCbYD54crPGGPMz4Wz+agR8K7bL04c8Lqqfiwi3wKTReQ6YCMwBEBVl4nIZOAH4BgwWlXL3tG+McaY4xa2SkFV1wGdihi+EyiykyFVfRR4NFw5GWOMKZl1iGeMMSag2D0FEfkvRZz9U0BVLwlLRsYYY6KmpOajv0QsC2OMMSeFYisFVZ0dyUSMMcZEX0nNR99RcvNRalgyMsYYEzUlNR9d5MaCu5dMcOMw4EDYMjLGGBM1JTUfbQAQkbNV9eygUfeLyFfAH8OdnDHGmMjyckpqDRE5p+CFiJwF1ChhemOMMeWUl4vXrgNeEZHaOMcY9gCjwpqVMcaYqCjpQPOZwDequhDoJCK1cO7pvCdi2RljjImokpqPRgALReRNERkJVLcKwRhjKraSDjTfDCAi7YD+wKtuE9LnwMfAV9ZhnTHGVCylHmhW1RWq+rSq9gPOB77E6dl0XriTM8YYE1mlHmgWkcRCg74BZqjq0fCkZIwxJlq8nJK6CMgGVgGr3efrRWSRiHQNZ3LGGGMiy0ul8DEwQFXrq2o9nOMLk4FbgefDmZwxxpjI8lIppKvqJwUvVHU6cJ6qfgNUDVtmxhhjIs7LxWs5InIf8Kb7+kpgl4jEAvlhy8wYY0zEedlTuBpIAqYC7wHN3WGxwBVhy8wYY0zElbqnoKo7gNuLGb3G33SMMcZEk5dTUtOB3wLJwdPb/RSMMabi8XJM4TVgDPAddgzBGGMqNC+VQraqTgt7JsYYY6LOS6XwoIi8DMwEDhcMVNV3wpaVMcaYqPBSKVwLtAPi+an5SAGrFIwxpoLxUil0UtWUsGdijDEm6rxcp/CNiLQPeybGGGOizkulcA6wWERWishSEflORJaGO7GwO3YEVn7kRGOMMYC35qN+Yc8iGtbOhEnD4coJcHr/aGdjjDEnBS9XNG+IRCKRltP4XOZ3eILujc+l8A0jjDGmsvLSfFQhTfnfNm7+thFT/rct2qkYY8xJw0vzUYU0JL1ZSDTGGFOJ9xQSa1Thpl6tSaxRJdqpGGNMmeTkHuHF2WvJyfX/RJlSKwURuUxEVovIHhHZKyL7RGSv75kYY4zxZPzcTB77aAXj52b6XraXPYUngEtUtbaq1lLVmqpay+sCRCRWRP4nIu+7rxNFZIZb0cwQkbpB044VkTXu6a99y/52jDGmMtBC0T9ejilsU9XlJ7CMXwPLgYKK5H5gpqo+LiL3u6/vcy+QGwp0AJoAn4pIW1XNO4FlG2NMhTPirJZUrxIXlmOiXvYUFojIJBG5ym1KukxELvNSuIgkAQOBl4MGXwqMd5+PBwYFDX9TVQ+r6nqcG/h097IcY4ypTMJ5TNTLnkIt4ABwYdAwrx3i/Q24F6gZNKyRqm4FUNWtItLQHd4U+CZouix3WAgRuRG4EaB58+YeUjDGGOOVl4vXrj2egkXkImC7qi4UkQwvsxS1+CLyeQl4CSA9Pf24G9Ryco8wZcEmhqQ3szOQjDHGVWylICL3quoTIvJ3iv5xvqOUss8GLhGRAUACUEtEJgLbRKSxu5fQGNjuTp8FBDeQJQFbyvBeymTKgk089tEKAG7q1TpcizHGmHKlpD2FgoPLC46nYFUdC4wFcPcUfqOq14jIk8AI4HE3vufOMg14XUT+inOguQ0w/3iW7YVdvGaMKbeOHXH6b2vdB+L8bekotlJQ1f+6cXxx0xynx4HJInIdsBEY4i5nmYhMBn4AjgGjw3nmUcGBGmOMKXfC2KFnRLq5UNVZwCz3+U6gTzHTPQo8GomcjDGm3GrWA864yIk+q7TdXBhjTLn1zfOw7F0n+qzSdohnjDHl1tl3hkYflXT2UZFnHRXwcPaRMcaYcIitAk27OtFnJe0pHNdZR8YYY8Js1ScwZQQMGQ/tL/a16JLOPgo560hEajqDdb+vGRhjjCmb/KOAutFfXrrO7igi/wO+B34QkYUi0sH3TIwxxnhz9DBovhN95uXso5eAu1W1hao2B+4B/uV7JsYYY7xZ9VFo9JGXSqGGqn5e8MK95qCG75lE2uH98NkjTjTGmPJkwBPQvKcTfealUlgnIr8TkWT38QCw3vdMIu2rv8EXTzrRGGPKky3/g6yFTvSZl+sURgF/4Keusr8Ajqvn1JNKz1th5xonGmNMedK6j9PFResiO4c4ISVWCiISC0xR1Qt8X3KU7Vv2MTWWTSU3+UJqdrs62ukYY4x3cVV87/OoQInNR26HdAdEpHZYlh5F2+e9haiyfd5b0U7FGGNOGl6ajw4B34nIDCC3YGB5v6K57mVPkjPhcupe9mS0UzHGmJOGlwPNHwC/wzmWsDDoUa4lLBlP4sH1JCzxu2dwY4wJr5zcI7w4ey05uUd8L9vL7TjHi0g1oLmqrvQ9gyj5cu+p9FHhy72nhtx82hhjTnZvz1/H/BmTicm/ght6t/O1bC9XNF8MLAY+dl+nicg0X7OIgvQ+l7G+4QWk97ks2qkYY0yZDE1czUtVn2Fo4mrfy/bSfPQQ0B3YDaCqi4GWvmcSYYlLX+a07BkkLn052qkYY0yZ1OzQl9ihE6nZoa/vZXs50HxMVfeISPCwYrvULjfC2B+5McaEVbROSXV9LyJXA7Ei0sa9z8LcsGQTSWHsj9wYY8IpZ1cOC8fdTc6uHN/L9lIp3A50AA4DrwN7gF/7nkmE7Vsylfw3rmbfkqnRTsUYY8pk/dRH6JL5b9ZPfcT3sr1UCgNV9f9UtZv7eAC4xPdMImz1kq8QzWf1kq+inYoxxpRJy0EPsCj5OloOesD3sr1UCmM9DitXave+g5VVO1C7d7m+Bs8YUwklVoWup+SQWNX/sku6R3N/YADQVESeDRpVCzjmfyqRtfezp0k7vIzFnz0N1z0X7XSMMca7aXfAiv9C3jEYOtHXoks6+2gLzn2aLyH0CuZ9wF2+ZhEFLRKrwSY3GmNMOXKgemOqAQerN6a6z2UX23ykqkvc+zSfpqrj3efTgDWqusvnPCJu/t5a5KsTjTGmPJl7oAn56kS/eTmmMENEaolIIrAEGCcif/U9kwjrsnkyMW40xpjypEuVLcSIE/3mpVKorap7gcuAcaraFSj391eYGjuAfDcaY0x5smnnflA3+sxLpRAnIo2BK4D3fc8gSs6JW0qMG40xpjxp0zgRxI0+81Ip/BH4BFirqt+KSCvA/16YIk0LRWOMKSeq97odad6T6r1u973sUisFVZ2iqqmqeov7ep2qXu57JhE2PfZs8t1ojN9y9u7n47dfIWev/7v3xjDvRdj4jRN95qXr7LYiMlNEvndfp4qI/5fRRVhc487szj+FuMado52KqYDmz5hMn6W/Yf4MO5HB+G9f7TbkI+yr3cb3sr00H/0L5wrmowCquhQY6nsmEXbLxrtIjN3PLRvL/SUX5iTU/RdXMDP1L3T/xRXRTsVUQHOXbUBUmbtsg+9le6kUqqvq/ELDyv0VzQfqtg2JxvgpsXoV+nVsTGJ164XX+K/Tpgkh0U9eKoUdItIa95CsiAwGtvqeSYQt2FEF1I3G+G3VJ/DmMCca47Ov8s4IiX7yUimMBl4E2onIZuBO4JbSZhKRBBGZLyJLRGSZiPzBHZ4oIjNEZLUb6wbNM1ZE1ojIShHx/5ZCQXoc/iIkGmNMedH2lIMh0U9ezj5ap6oXAA2Adqp6jqpmeij7MHC+qnYC0oB+ItITuB+YqaptgJnua0SkPc6xig5AP+B5EYkt+1vyZjFtQ6IxvmrbF4a+5kRjfJbVKAN1o99K6iX17mKGA6CqJXZ1oaoKFJyPF+8+FLgUyHCHjwdmAfe5w99U1cPAehFZg3Nv6K89vZMyWp7XjDNjlrE8rxlnhWMBxhgTJmnxWxBxot9K2lOoWcqjVCISKyKLge3ADFWdBzRS1a0AbmzoTt4U2BQ0e5Y7rHCZN4rIAhFZkJ2d7SWNonNL6kY+giR1O+4yjCnW2pkwabgTjfFZ1trvQd3os2L3FFT1DydauKrmAWkiUgd4V0Q6ljC5FFVEEWW+BLwEkJ6eftzXIzdJ6cXarFdoktLreIswpnit+8CVE5xojM/yNTT6ycvFa0ki8q6IbBeRbSLytogklWUhqrobp5moH7DN7UsJN253J8sCmgXNloRzT4ewaD59FG1iNtN8+qhwLcIYY8KiRavTQdzoMy9nH43DuY9CE5zmnP+6w0okIg3cPQREpBpOz6or3LJGuJONAN5zn08DhopIVRFpCbQBCl8f4Zu2ZIVEY3xlp6SaMDqQnQnqRp+VdOe1Ag1UNbgSeFVE7vQwX2NgvHsGUQwwWVXfF5Gvgckich2wERgCoKrLRGQy8APOxXGj3eanMMkvFI3xUd4xUHWiMT6rWqVKSPSTl0phh4hcA7zhvr4K2FnaTG53GD/rWEhVdwJFNrSq6qPAox5yOmG7pR71dCe7pR71I7FAU7nExoGIE43x2bXbBvNc/gpu2zaY6T6X7aX5aBTOvRR+xLmSebA7rFzLoVZINMZXrXvDuXc70Rif/bXlQtrEbOGvLRf6XraXi9c2quolqtpAVRuq6iBV9b8XpgjL1YSQaIyvMufAl39zojE+m7k1AXWj30rdtxWRBsANQHLw9KparvcWTmdjSDTGV816wBkXOdEYn11waCriRviTr2V7afB8D5gDfAqE8cBvhNVpBrtXONEYnx2Y9TeqLXuXgzWaU33AH6OdjqlgTufHkOgnL5VCdVW9z/clR1nWoXjaqBOt9yPjty/3N6aPCl/ub8yF0U7GVDg/xiTRNG8DP8YkUaaLxjzwcqD5fREZ4PNyo67ZoWUh0Rg/pf/iShYnjyL9F1dGOxVTAe3NjwmJfiq2RBHZJyJ7gV/jVAwHRWRv0PBybSd1Q6IxfkrcPo+um/5D4vZ50U7FVEC1aieGRD+V1PeRp07vyqt1eafSNGYn6/JO9X33yxjr+8iEU409y0Oin7z0ffRLEakd9LqOiAzyPZMIq01uSDTGV3lHYPNCJxrjs3fiLw+JfvLSIPWgqu4peOF2bveg75lE2L7Y2iHRGF999Tf44kknGuOzQdUXhUQ/eTn7qKiKo9xfu59YVeCQG43x29l3hkZjfLRtz34S3VjP57K97CksEJG/ikhrEWklIk8D/l9bHWFND68Iicb4KrYKNO3qRGN8dorb7H1KGJq/vVQKtwNHgEnAZOAgMNr3TCIs1r2nT2yR9/Yx5sTsW/YJeW9ew75l1nW2KV9KbQZS1Vzg/gjkElFHiKOGG43x2+QdyeQeuYgaO5K5LtrJmArnINVCop8q7S/iIRKAfW40xl9X1M+kepX3OVB/MJAS7XRMBVM3bxfEuNFn/l8OV06s0wYh0Rg/HW3ajfX1e3O0abdop2IqoDqxh0OinyptpdBMtodEY/y0+t0/0Xr7DFa/628PlsYA/FizQ0j0k5eus1viHGxOJrTr7Et8zyaC4vNyIcaNxvgsu3pr8t1ojN8a5e8AcaPPvBxTmAr8G/gvFeiGxrHuW4+tvIdVTBj1at8cVsc40Rif7clLoIE6saHPZXv5RTykqs/6vNyo2xNbh4a6jz2xdWgU7WRMhVMzZSBUf52a1veRCYOEvL0h0U9eKoVnRORBYDoQOKqhqv5fXx1BR923ftT2FEw4xFWB0/tHOwtTQcUf3RES/eTlFzEFGA6cz0/NR+q+LrdW0pj2rGcljfH/UI0xxoTPjjpdSNq1gB11uuD3vSO9nH30S6CVqvZS1d7uo1xXCAAX5s0Picb4KSf3CC/OXktOrvWSavw35ci55LvRb14qhSVAHd+XHGXxscdCojF+mrJgE499tIIpCzZFOxVTAe0/sDsk+slL81EjYIWIfEvoMYVyfUrqburSkF3spq4daDa+u6B9I75Zt5ML2tvWZfx3fd5UYmKcCH/ztWwvlUK5v3dCUepwMCQa46fPvs8iZvXHfPZ9LVr3bhftdEwF815sBjfpu7wXm8EtPpftpUO82T4v86SwnRokcYjt1PD9QI0xQxNXM6rqMxxI7AxYpWD8NYIPETf6zcsVzftwzjYCqALEA7mqWsv3bCIov1A0xk81O/SFhIl2nYIJi8004DSy2EwD2vhcdqkHmlW1pqrWch8JwOXAcz7nEXGHOCUkGuOnnMPw4o9tyfG/vzJjaJKQB+JGn5W5QzxVnUo5v0YBQCQ2JBrjJzv7yITTwUMHQN3oMy/NR5cFvYwB0vmpOancOhhfBw670RifDUlvFhKN8dOCKmfS9/DHLKhyJv18LtvL2UcXBz0/BmQCl/qcR8TtPnQQxI3G+C3vCC12zIa8K3AOxRnjn/anVoWNbvRZiZWCOG0rS1X1ad+XHGVtNQvEjcb4bP6MyfRZ+htmAv0uHxXtdEwFU3/jR6Bu9FmJxxRUNQ8o1xepFecQ8SHRGD91P38QS5NH0v38QdFOxVRAx6rWD4l+8nKgea6IPCci54pIl4JHaTOJSDMR+VxElovIMhH5tTs8UURmiMhqN9YNmmesiKwRkZUi0vcE3lep8mJrhURj/JS49Su6bnyVxK1fRTsVUwG9fOyikOgnL8cUznLjH4OGeekl9Rhwj6ouEpGawEIRmQGMBGaq6uMicj9wP3CfiLQHhgIdgCbApyLS1t1b8d2WvBq0jnHiaeFYgKnUtu49SENVtu89SONoJ2MqnMtOOwor3egzL1c09z6eglV1K7DVfb5PRJYDTXEOUme4k40HZgH3ucPfVNXDwHoRWQN0B74+nuWXpkXsNlA3GuOz3/3QFD1yF/JDU17uGe1sTEWTtWM/LdSJyT6XXWrzkYg0EpF/i8hH7uv2InJdWRYiIslAZ2Ae0MitMAoqjoK7yTUFgk/qznKHFS7rRhFZICILsrOzy5JGiISYuJBojJ/GDjiDlKa1GTvgjGinYiqgTTsPhkQ/eTmm8CrwCU6TDsAq4E6vCxCRU4C3gTtVtaR7x0kRw352PYSqvqSq6aqa3qBBA69p/MwPefVDojF+ap3zFXfu+AOtc+yYgvFflWoJIdFPXiqF+qo6GbebIFU9Bnhq5xeReJwK4TVVfccdvE1EGrvjGwPb3eFZENI3XRKwxctyjsdmaRkSjfHTvtx95Gs++3L3RTsVUwG1HXAXmTFJtB1wl+9le6kUckWkHu6/dhHpCewpbSYREeDfwHJV/WvQqGnACPf5COC9oOFDRaSqiLQE2gBhuy1ax2GPMbHqEDoOeyxcizCV2Orv5iOqrP7O7uxn/Fd99gMk52dRffYDvpftpUH9bpwf7NYi8hXQABjsYb6zce7t/J2ILHaH/RZ4HJjsHpfYCAwBUNVlIjIZ+AHnzKXR4TrzCCCteV3SeiZD87qlTmtMWSX/8vcsmhpDy0H+f2mNqVG/Gexwo8+8nH20SER6AafjtPuvVNVSz4NS1S8p+jgBQJH9Cavqo8CjpZXthwPTH6Hawhc4mLuf6hc/HolFmsqkyiksaDWallWsF17jv0UbdtNfnTjA57K9nH00BKimqsuAQcAkLxevnew2r10C6kZjfDZ+biaPfbSC8XMzo52KqYDOrbcHxI0+83JM4XfudQbnAH1xri14wfdMIuzYsbyQaIyfYvKP0CdmITH5R6KdiqmAauZuRNzoNy+VQsGv5kDgBVV9jwrQ7WPuoWMh0Rg/jTp1PS9VfYZRp66PdiqmIjrt/NDoIy8HmjeLyIvABcCfRaQqx3FznpNN9arxcNSNxvjMbsdpwmlt9c4kI2RW70xrn8v28uN+Bc7Fa/1UdTeQCIzxOY+I21etcUg0xpjy4sb5jbjxyN3cOL+R72V7uUfzAVV9R1VXu6+3qup03zOJsLrHskOiMb5aOxMmDXeiMT7rnpwYEv1U7puBjldS81YgbjTGb637wJUTnGiMz849Opt/xT/FuUdn+152pa0U4rctDYnGGFNeZOx9HxEn+q3SVgqL5XRQNxrjN2s+MmE0+5QBqDrRb5W23+j5za7nx6xMNja7nm7RTsZUPNZ8ZMLoWFx1FOFYXHXfy660lcJl9TJpGPct2+tl4tzLxxhjyoeMM5rAKnGizypt89F3G3eCqhON8Zs1H5kwqpkykNirXqdmykDfy660lcKZreshIpzZul60UzEVkTUfmXCKqwKn93eizyptpVCzWjViRKhZrVq0UzEVURi/tMZw7Ais/MiJPqu0lQJt+8LQ15xojDHlSRibJyttpbB21xGunVuPtbusF0tjTDkTxubJSlspPPXe11y29gGeeu/raKdijDFlY8cU/PfnKuO4KG4ef64yLtqpGGNM2RzeD5894kSfVdpKgTMuQhE446JoZ2KMMWUz5yn44kkn+qzSVgpv7u/CDUfu5s395f7OosaYyqbhGYC40V+V9ormy7sm0Swnie5dk6KdijHGlE1cNZAYJ/qs0u4pJG6dQ79l95K4dU60UzHGmLIJ4yn1lXZPgeRz4Zw7nWiMMeVJwdlHYVBp9xRY+znM+asTjTHGAJW5UjDGGPMzlbdSsG4ujDHlVE7uEV6cvZacXOv7yD/WYZkxppyasmATj320gikLNvleduU90GyMMeXUkPRmIdFPlXdPIYxdzxpjTDgl1qjCTb1ak1jD+j7yz6pP4M1hTjTGGANU5krBGGPMz1TeSsHOPjLGlFd257UwsLOPjDHlld15zRhjTEB5vPOaiLwiIttF5PugYYkiMkNEVruxbtC4sSKyRkRWioi16RhjTBSEc0/hVaBfoWH3AzNVtQ0w032NiLQHhgId3HmeF5HYMOZmjDHlV3lsPlLVL4CcQoMvBca7z8cDg4KGv6mqh1V1PbAG6B6u3Iwxplwrj81HxWikqlsB3NjQHd4UCL5eO8sd9jMicqOILBCRBdnZ2WFN1hhjTkphPFHmZDnQLEUM06ImVNWXVDVdVdMbNGgQ5rSMMaZyiXSlsE1EGgO4cbs7PAsI7sQjCdgS4dyMMabSi3SlMA0Y4T4fAbwXNHyoiFQVkZZAG2B+hHMzxphKL2y9pIrIG0AGUF9EsoAHgceBySJyHbARGAKgqstEZDLwA3AMGK2qeeHKzRhjTNHCVimo6lXFjCrycLmqPgo8Gq58jDHGlO5kOdBsjDHmJGCVgjHGmABRLfLMz3JBRLKBDSdQRH1gh0/p+MnyKhvLq2wsr7KpiHm1UNUiz+kv15XCiRKRBaqaHu08CrO8ysbyKhvLq2wqW17WfGSMMSbAKgVjjDEBlb1SeCnaCRTD8ioby6tsLK+yqVR5VepjCsYYY0JV9j0FY4wxQaxSMMYYE1AhKwUR6efe1nONiNxfxHgRkWfd8UtFpIvXecOc1zA3n6UiMldEOgWNyxSR70RksYgsiHBeGSKyx132YhH5vdd5w5zXmKCcvheRPBFJdMeFc3397FazhcZHa/sqLa9obV+l5RWt7au0vCK+fYlIMxH5XESWi8gyEfl1EdOEd/tS1Qr1AGKBtUAroAqwBGhfaJoBwEc493HoCczzOm+Y8zoLqOs+71+Ql/s6E6gfpfWVAbx/PPOGM69C018MfBbu9eWWfR7QBfi+mPER37485hXx7ctjXhHfvrzkFY3tC2gMdHGf1wRWRfr3qyLuKXQH1qjqOlU9AryJc7vPYJcC/1HHN0Adce7v4GXesOWlqnNVdZf78huc+0qE24m856iur0KuAt7wadkl0qJvNRssGttXqXlFafvysr6KE9X1VUhEti9V3aqqi9zn+4Dl/PwulGHdvipipeDl1p7FTeP5tqBhyivYdTj/BgooMF1EForIjT7lVJa8zhSRJSLykYh0KOO84cwLEakO9APeDhocrvXlRTS2r7KK1PblVaS3L8+itX2JSDLQGZhXaFRYt6+wdZ0dRV5u7VncNJ5vC3ocPJctIr1xvrTnBA0+W1W3iEhDYIaIrHD/6UQir0U4faXsF5EBwFScGyGdFOsLZ9f+K1UN/tcXrvXlRTS2L88ivH15EY3tqywivn2JyCk4ldCdqrq38OgiZvFt+6qIewpebu1Z3DThvC2op7JFJBV4GbhUVXcWDFfVLW7cDryLs6sYkbxUda+q7neffwjEi0h9L/OGM68gQym0ax/G9eVFNLYvT6KwfZUqSttXWUR0+xKReJwK4TVVfaeIScK7ffl9oCTaD5y9n3VAS3462NKh0DQDCT1QM9/rvGHOqzmwBjir0PAaQM2g53OBfhHM61R+utCxO85d8yTa68udrjZOu3CNSKyvoGUkU/yB04hvXx7zivj25TGviG9fXvKKxvblvu//AH8rYZqwbl8VrvlIVY+JyG3AJzhH419R53afN7vj/wl8iHMEfw1wALi2pHkjmNfvgXrA8yICcEydXhAbAe+6w+KA11X14wjmNRi4RUSOAQeBoepshdFeXwC/BKaram7Q7GFbX1DsrWbjg/KK+PblMa+Ib18e84r49uUxL4j89nU2MBz4TkQWu8N+i1OhR2T7sm4ujDHGBFTEYwrGGGOOk1UKxhhjAqxSMMYYE2CVgjHGmACrFIwxxgRYpWBMGYlIHRG51X2eISLvl3H+V0VkcHiyM+bEWKVgTNnVAW6NdhLGhEOFu3jNmAh4HGjtXlx0FMgVkbeAjsBC4BpVVfe+ABcD1XCuer1J7cIgc5KzPQVjyu5+YK2qpgFjcHqyvBNoj9OX/dnudM+pajdV7YhTMVwU+VSNKRurFIw5cfNVNUtV84HFOP3pAPQWkXki8h1wPtChmPmNOWlY85ExJ+5w0PM8IE5EEoDngXRV3SQiDwEJ0UjOmLKwPQVjym4fzq0SS1JQAexw+8a3s41MuWB7CsaUkaruFJGv3Bu+HwS2FTHNbhH5F/Adzv18v41slsYcH+sl1RhjTIA1HxljjAmwSsEYY0yAVQrGGGMCrFIwxhgTYJWCMcaYAKsUjDHGBFilYIwxJuD/AX3UKXH7pK6wAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#===========================================#\n",
    "# Visualisation of the resultts\n",
    "#===========================================# \n",
    "\n",
    "# Choose the variables number\n",
    "var1 = -1\n",
    "var2 = 4\n",
    "\n",
    "\n",
    "# Check for 1D distribution\n",
    "plt.hist(y_np[:,var1], density = True)\n",
    "plt.hist(y_new[:,var1], density = True)\n",
    "plt.xlabel(varnames[var1])\n",
    "plt.ylabel('Number of observations')\n",
    "plt.legend(['Genuine observations', 'Synthetic observations'])\n",
    "plt.title('1D projections of genuine and synthetic obervations')\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# Check for 2D distribution\n",
    "plt.scatter(y_np[:,var1], y_np[:,var2],  label = 'Genuine observations', s = 1)\n",
    "plt.scatter(y_new[:,var1], y_new[:,var2], label = 'Synthetic observations', s = 0.5)\n",
    "plt.xlabel(varnames[var1])\n",
    "plt.ylabel(varnames[var2])\n",
    "plt.legend(['Genuine observations', 'Synthetic observations'])\n",
    "plt.title('2D projections of genuine and synthetic obervations')\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}