dqr-joint-200S.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UIgs75GktB8G"
   },
   "source": [
    "# Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 173
    },
    "colab_type": "code",
    "id": "CJVCWfrMtDbL",
    "outputId": "c567c68f-0a6a-47f3-faed-e2b40426f770"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found GPU at: /device:GPU:0\n",
      "Tensorflow Version .2.3.1\n",
      "Numpy Version .1.16.2\n"
     ]
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from datetime import datetime\n",
    "import tensorflow.keras.backend as K\n",
    "from tensorflow.keras.callbacks import EarlyStopping, ModelCheckpoint\n",
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import BatchNormalization, ConvLSTM2D, Flatten, RepeatVector, Reshape, TimeDistributed, Lambda, Dense, Dropout\n",
    "from scipy import stats \n",
    "from scipy.optimize import curve_fit \n",
    "from scipy.stats import norm, lognorm\n",
    "import tensorflow as tf\n",
    "from pathlib import Path\n",
    "\n",
    "from load import load_links\n",
    "from DQR.models import tilted_loss, convLSTM\n",
    "\n",
    "device_name = tf.test.gpu_device_name()\n",
    "if device_name != '/device:GPU:0':\n",
    "    print('GPU device not found')\n",
    "else:\n",
    "    print('Found GPU at: {}'.format(device_name))\n",
    "print(\"Tensorflow Version .{}\".format(tf.__version__))\n",
    "print(\"Numpy Version .{}\".format(np.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "eEWMQ8y0tEqg"
   },
   "source": [
    "# Data Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def roll(ix, ts, lags, preds):\n",
    "    X = np.stack([np.roll(ts, i, axis=0) for i in range(lags, 0, -1)], axis=1)[lags:-preds,]\n",
    "    Y = np.stack([np.roll(ts, -i, axis=0) for i in range(0, preds, 1)], axis=1)[lags:-preds,]\n",
    "    Y_ix = ix[lags:-preds]\n",
    "    return X, Y, Y_ix\n",
    "\n",
    "def load_data(line, freq, lags, preds):\n",
    "    train = pd.read_csv(f'Data/{line}_{freq}_train_norm.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    train_counts = pd.read_csv(f'Data/{line}_{freq}_train_sample_count.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    X_train, Y_train, Y_train_ix = roll(train.index, train.values, lags, preds)\n",
    "    \n",
    "    val = pd.read_csv(f'Data/{line}_{freq}_val_norm.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    val_counts = pd.read_csv(f'Data/{line}_{freq}_val_sample_count.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    X_val, Y_val, Y_val_ix = roll(val.index, val.values, lags, preds)\n",
    "\n",
    "    test = pd.read_csv(f'Data/{line}_{freq}_test_norm.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    test_counts = pd.read_csv(f'Data/{line}_{freq}_test_sample_count.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    test_mean = pd.read_csv(f'Data/{line}_{freq}_test_mean.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    test_std = pd.read_csv(f'Data/{line}_{freq}_test_std.csv.gz', index_col=0, compression='gzip', parse_dates=True)\n",
    "    X_test, Y_test, Y_test_ix = roll(test.index, test.values, lags, preds)\n",
    "    _, mu_test, _ = roll(test_mean.index, test_mean.values, lags, preds)\n",
    "    _, sigma_test, _ = roll(test_std.index, test_std.values, lags, preds)\n",
    "        \n",
    "    return (X_train, Y_train), (X_val, Y_val), (X_test, Y_test, Y_test_ix, mu_test, sigma_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "WmmJZeJQt9y9"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_train (8701, 32, 29, 1, 1) y_train (8701, 3, 29, 1, 1)\n",
      "X_val (1309, 32, 29, 1, 1) y_val (1309, 3, 29, 1, 1)\n",
      "X_test (1309, 32, 29, 1, 1) y_test (1309, 3, 29, 1, 1)\n"
     ]
    }
   ],
   "source": [
    "line = '200S'\n",
    "freq = '15min'\n",
    "lags = 32\n",
    "preds = 3\n",
    "\n",
    "\n",
    "data_train, data_val, data_test = load_data(line, freq, lags, preds)\n",
    "\n",
    "# Unpack tensors\n",
    "(X_train, y_train) = data_train\n",
    "(X_val, y_val) = data_val\n",
    "(X_test, y_test, y_ix_test, y_mean_test, y_std_test) = data_test\n",
    "\n",
    "X_train = X_train[..., np.newaxis, np.newaxis]\n",
    "y_train = y_train[..., np.newaxis, np.newaxis]\n",
    "X_val = X_val[..., np.newaxis, np.newaxis]\n",
    "y_val = y_val[..., np.newaxis, np.newaxis]\n",
    "X_test = X_test[..., np.newaxis, np.newaxis]\n",
    "y_test = y_test[..., np.newaxis, np.newaxis]\n",
    "\n",
    "print('X_train', X_train.shape, 'y_train', y_train.shape)\n",
    "print('X_val', X_val.shape, 'y_val', y_val.shape)\n",
    "print('X_test', X_test.shape, 'y_test', y_test.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "CYQ321P8E665"
   },
   "outputs": [],
   "source": [
    "quantiles = np.sort([0.025, 0.975, 0.05, 0.95, 0.10, 0.90, 0.20, 0.8, 0.4, 0.6])\n",
    "y_traink = np.zeros((y_train.shape[0], y_train.shape[1], y_train.shape[2], 1, len(quantiles)+1))\n",
    "y_valk  = np.zeros((y_val.shape[0], y_train.shape[1], y_train.shape[2], 1, len(quantiles)+1))\n",
    "for i in range(len(quantiles)+1):\n",
    "    y_traink[:,:,:,:,i] = y_train[:,:,:,:,0]\n",
    "    y_valk[:,:,:,:,i] = y_val[:,:,:,:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.025, 0.05 , 0.1  , 0.2  , 0.4  , 0.6  , 0.8  , 0.9  , 0.95 ,\n",
       "       0.975])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quantiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "QC4IBOREtLAa"
   },
   "source": [
    "# Joint Model (Encoder-Decoder ConvLSTM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "aGCWTNA71_5U"
   },
   "outputs": [],
   "source": [
    "def joint_convLstm(num_filters, kernel_length, input_timesteps, num_links, output_timesteps, quantiles, prob, loss, opt):\n",
    "    model = Sequential()\n",
    "    model.add(BatchNormalization(name = 'batch_norm_0', input_shape = (input_timesteps, num_links, 1, 1)))\n",
    "    model.add(ConvLSTM2D(name ='conv_lstm_1',\n",
    "                         filters = num_filters, kernel_size = (kernel_length, 1), \n",
    "                         padding='same',\n",
    "                         return_sequences=False))\n",
    "\n",
    "    model.add(Dropout(prob, name = 'dropout_1'))\n",
    "    model.add(BatchNormalization(name = 'batch_norm_1'))\n",
    "\n",
    "    model.add(Flatten())\n",
    "    model.add(RepeatVector(output_timesteps))\n",
    "    model.add(Reshape((output_timesteps, num_links, 1, num_filters)))\n",
    "\n",
    "    model.add(ConvLSTM2D(name ='conv_lstm_2',filters = num_filters, kernel_size = (kernel_length, 1), padding='same',return_sequences = True))\n",
    "    model.add(Dropout(prob, name = 'dropout_2'))\n",
    "\n",
    "    model.add(TimeDistributed(Dense(units = len(quantiles) + 1, name = 'dense_1')))\n",
    "    model.compile(loss = loss, optimizer = opt)\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "lbaB6gJHtVdV"
   },
   "source": [
    "# Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tilted loss for both mean and quantiles\n",
    "def joint_tilted_loss(quantiles, y, f):\n",
    "    loss = K.mean(K.square(y[:,:,:,0,0]-f[:,:,:,0,0]), axis = -1)\n",
    "    for k in range(len(quantiles)):\n",
    "        q = quantiles[k]\n",
    "        e = (y[:,:,:,0,k+1]-f[:,:,:,0,k+1])\n",
    "        loss += K.mean(K.maximum(q*e, (q-1)*e))\n",
    "    return loss\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 153
    },
    "colab_type": "code",
    "id": "_IsTBrgFz4P-",
    "outputId": "bff78231-3bf3-49c8-cf0e-f6be42839ff2",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 1/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 3.8147 - val_loss: 4.1298\n",
      "Epoch 2/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 3.5397 - val_loss: 3.9641\n",
      "Epoch 3/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 3.3562 - val_loss: 3.6331\n",
      "Epoch 4/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 3.2106 - val_loss: 3.3449\n",
      "Epoch 5/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 3.1003 - val_loss: 3.2039\n",
      "Epoch 6/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 3.0147 - val_loss: 3.1208\n",
      "Epoch 7/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.9491 - val_loss: 3.0608\n",
      "Epoch 8/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.8970 - val_loss: 3.0144\n",
      "Epoch 9/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.8564 - val_loss: 2.9771\n",
      "Epoch 10/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.8213 - val_loss: 2.9431\n",
      "Epoch 11/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.7926 - val_loss: 2.9182\n",
      "Epoch 12/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.7678 - val_loss: 2.8945\n",
      "Epoch 13/200\n",
      "175/175 [==============================] - 28s 163ms/step - loss: 2.7470 - val_loss: 2.8770\n",
      "Epoch 14/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.7282 - val_loss: 2.8576\n",
      "Epoch 15/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.7111 - val_loss: 2.8436\n",
      "Epoch 16/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.6961 - val_loss: 2.8273\n",
      "Epoch 17/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.6834 - val_loss: 2.8154\n",
      "Epoch 18/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.6701 - val_loss: 2.8041\n",
      "Epoch 19/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.6604 - val_loss: 2.7929\n",
      "Epoch 20/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.6505 - val_loss: 2.7828\n",
      "Epoch 21/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.6404 - val_loss: 2.7736\n",
      "Epoch 22/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.6325 - val_loss: 2.7652\n",
      "Epoch 23/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.6246 - val_loss: 2.7588\n",
      "Epoch 24/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.6173 - val_loss: 2.7504\n",
      "Epoch 25/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.6104 - val_loss: 2.7433\n",
      "Epoch 26/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.6048 - val_loss: 2.7382\n",
      "Epoch 27/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5991 - val_loss: 2.7319\n",
      "Epoch 28/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5930 - val_loss: 2.7265\n",
      "Epoch 29/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5881 - val_loss: 2.7214\n",
      "Epoch 30/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5837 - val_loss: 2.7163\n",
      "Epoch 31/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5784 - val_loss: 2.7113\n",
      "Epoch 32/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5745 - val_loss: 2.7072\n",
      "Epoch 33/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5707 - val_loss: 2.7034\n",
      "Epoch 34/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5665 - val_loss: 2.6999\n",
      "Epoch 35/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 2.5632 - val_loss: 2.6959\n",
      "Epoch 36/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5602 - val_loss: 2.6927\n",
      "Epoch 37/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5560 - val_loss: 2.6879\n",
      "Epoch 38/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5535 - val_loss: 2.6862\n",
      "Epoch 39/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5505 - val_loss: 2.6836\n",
      "Epoch 40/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5477 - val_loss: 2.6805\n",
      "Epoch 41/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5441 - val_loss: 2.6777\n",
      "Epoch 42/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5422 - val_loss: 2.6754\n",
      "Epoch 43/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.5393 - val_loss: 2.6726\n",
      "Epoch 44/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.5377 - val_loss: 2.6701\n",
      "Epoch 45/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5353 - val_loss: 2.6687\n",
      "Epoch 46/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.5326 - val_loss: 2.6670\n",
      "Epoch 47/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5305 - val_loss: 2.6650\n",
      "Epoch 48/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5286 - val_loss: 2.6625\n",
      "Epoch 49/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5263 - val_loss: 2.6612\n",
      "Epoch 50/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5257 - val_loss: 2.6594\n",
      "Epoch 51/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5232 - val_loss: 2.6566\n",
      "Epoch 52/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5216 - val_loss: 2.6551\n",
      "Epoch 53/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 2.5208 - val_loss: 2.6537\n",
      "Epoch 54/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5193 - val_loss: 2.6533\n",
      "Epoch 55/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5175 - val_loss: 2.6504\n",
      "Epoch 56/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5163 - val_loss: 2.6500\n",
      "Epoch 57/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5154 - val_loss: 2.6487\n",
      "Epoch 58/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5135 - val_loss: 2.6477\n",
      "Epoch 59/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5127 - val_loss: 2.6467\n",
      "Epoch 60/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5106 - val_loss: 2.6455\n",
      "Epoch 61/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5100 - val_loss: 2.6437\n",
      "Epoch 62/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5087 - val_loss: 2.6432\n",
      "Epoch 63/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.5077 - val_loss: 2.6419\n",
      "Epoch 64/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5067 - val_loss: 2.6415\n",
      "Epoch 65/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5064 - val_loss: 2.6399\n",
      "Epoch 66/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.5057 - val_loss: 2.6387\n",
      "Epoch 67/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.5046 - val_loss: 2.6383\n",
      "Epoch 68/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5028 - val_loss: 2.6366\n",
      "Epoch 69/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5014 - val_loss: 2.6360\n",
      "Epoch 70/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.5015 - val_loss: 2.6349\n",
      "Epoch 71/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4995 - val_loss: 2.6339\n",
      "Epoch 72/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4988 - val_loss: 2.6333\n",
      "Epoch 73/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4986 - val_loss: 2.6328\n",
      "Epoch 74/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4968 - val_loss: 2.6317\n",
      "Epoch 75/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4965 - val_loss: 2.6308\n",
      "Epoch 76/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4958 - val_loss: 2.6297\n",
      "Epoch 77/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4949 - val_loss: 2.6290\n",
      "Epoch 78/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4940 - val_loss: 2.6287\n",
      "Epoch 79/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4934 - val_loss: 2.6282\n",
      "Epoch 80/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4926 - val_loss: 2.6263\n",
      "Epoch 81/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4912 - val_loss: 2.6257\n",
      "Epoch 82/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4914 - val_loss: 2.6250\n",
      "Epoch 83/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4907 - val_loss: 2.6246\n",
      "Epoch 84/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4893 - val_loss: 2.6237\n",
      "Epoch 85/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4892 - val_loss: 2.6230\n",
      "Epoch 86/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4878 - val_loss: 2.6217\n",
      "Epoch 87/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4871 - val_loss: 2.6222\n",
      "Epoch 88/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4870 - val_loss: 2.6201\n",
      "Epoch 89/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4871 - val_loss: 2.6207\n",
      "Epoch 90/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.4856 - val_loss: 2.6201\n",
      "Epoch 91/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4848 - val_loss: 2.6193\n",
      "Epoch 92/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4837 - val_loss: 2.6184\n",
      "Epoch 93/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.4835 - val_loss: 2.6173\n",
      "Epoch 94/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4829 - val_loss: 2.6177\n",
      "Epoch 95/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4821 - val_loss: 2.6160\n",
      "Epoch 96/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4814 - val_loss: 2.6155\n",
      "Epoch 97/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 2.4814 - val_loss: 2.6150\n",
      "Epoch 98/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4798 - val_loss: 2.6149\n",
      "Epoch 99/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.4799 - val_loss: 2.6141\n",
      "Epoch 100/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4792 - val_loss: 2.6133\n",
      "Epoch 101/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4790 - val_loss: 2.6127\n",
      "Epoch 102/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4781 - val_loss: 2.6125\n",
      "Epoch 103/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4776 - val_loss: 2.6116\n",
      "Epoch 104/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4776 - val_loss: 2.6110\n",
      "Epoch 105/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4757 - val_loss: 2.6104\n",
      "Epoch 106/200\n",
      "175/175 [==============================] - 28s 163ms/step - loss: 2.4760 - val_loss: 2.6102\n",
      "Epoch 107/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4737 - val_loss: 2.6098\n",
      "Epoch 108/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4751 - val_loss: 2.6088\n",
      "Epoch 109/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4744 - val_loss: 2.6084\n",
      "Epoch 110/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4735 - val_loss: 2.6075\n",
      "Epoch 111/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4734 - val_loss: 2.6075\n",
      "Epoch 112/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4722 - val_loss: 2.6069\n",
      "Epoch 113/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4721 - val_loss: 2.6068\n",
      "Epoch 114/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4715 - val_loss: 2.6065\n",
      "Epoch 115/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4715 - val_loss: 2.6058\n",
      "Epoch 116/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4712 - val_loss: 2.6053\n",
      "Epoch 117/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4694 - val_loss: 2.6049\n",
      "Epoch 118/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.4693 - val_loss: 2.6039\n",
      "Epoch 119/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4694 - val_loss: 2.6032\n",
      "Epoch 120/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4694 - val_loss: 2.6037\n",
      "Epoch 121/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4678 - val_loss: 2.6026\n",
      "Epoch 122/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4683 - val_loss: 2.6022\n",
      "Epoch 123/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4676 - val_loss: 2.6023\n",
      "Epoch 124/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.4666 - val_loss: 2.6016\n",
      "Epoch 125/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4662 - val_loss: 2.6007\n",
      "Epoch 126/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.4664 - val_loss: 2.6000\n",
      "Epoch 127/200\n",
      "175/175 [==============================] - 29s 168ms/step - loss: 2.4662 - val_loss: 2.6002\n",
      "Epoch 128/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4646 - val_loss: 2.5996\n",
      "Epoch 129/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4644 - val_loss: 2.5991\n",
      "Epoch 130/200\n",
      "175/175 [==============================] - 28s 163ms/step - loss: 2.4645 - val_loss: 2.5987\n",
      "Epoch 131/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4637 - val_loss: 2.5984\n",
      "Epoch 132/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4642 - val_loss: 2.5981\n",
      "Epoch 133/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4636 - val_loss: 2.5980\n",
      "Epoch 134/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4617 - val_loss: 2.5974\n",
      "Epoch 135/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4622 - val_loss: 2.5964\n",
      "Epoch 136/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4619 - val_loss: 2.5961\n",
      "Epoch 137/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4614 - val_loss: 2.5961\n",
      "Epoch 138/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 2.4610 - val_loss: 2.5960\n",
      "Epoch 139/200\n",
      "175/175 [==============================] - 29s 167ms/step - loss: 2.4611 - val_loss: 2.5952\n",
      "Epoch 140/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4599 - val_loss: 2.5952\n",
      "Epoch 141/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4603 - val_loss: 2.5946\n",
      "Epoch 142/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4594 - val_loss: 2.5947\n",
      "Epoch 143/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4591 - val_loss: 2.5936\n",
      "Epoch 144/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4591 - val_loss: 2.5934\n",
      "Epoch 145/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4586 - val_loss: 2.5928\n",
      "Epoch 146/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4582 - val_loss: 2.5923\n",
      "Epoch 147/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4579 - val_loss: 2.5926\n",
      "Epoch 148/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 2.4575 - val_loss: 2.5922\n",
      "Epoch 149/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 2.4564 - val_loss: 2.5919\n",
      "Epoch 150/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4555 - val_loss: 2.5911\n",
      "Epoch 151/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4562 - val_loss: 2.5910\n",
      "Epoch 152/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4554 - val_loss: 2.5908\n",
      "Epoch 153/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4558 - val_loss: 2.5908\n",
      "Epoch 154/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4551 - val_loss: 2.5904\n",
      "Epoch 155/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4547 - val_loss: 2.5895\n",
      "Epoch 156/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4541 - val_loss: 2.5894\n",
      "Epoch 157/200\n",
      "175/175 [==============================] - 29s 166ms/step - loss: 2.4545 - val_loss: 2.5887\n",
      "Epoch 158/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4542 - val_loss: 2.5887\n",
      "Epoch 159/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4538 - val_loss: 2.5891\n",
      "Epoch 160/200\n",
      "175/175 [==============================] - 28s 162ms/step - loss: 2.4537 - val_loss: 2.5879\n",
      "Epoch 161/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4529 - val_loss: 2.5882\n",
      "Epoch 162/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4529 - val_loss: 2.5876\n",
      "Epoch 163/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4521 - val_loss: 2.5873\n",
      "Epoch 164/200\n",
      "175/175 [==============================] - 29s 163ms/step - loss: 2.4522 - val_loss: 2.5871\n",
      "Epoch 165/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4522 - val_loss: 2.5871\n",
      "Epoch 166/200\n",
      "175/175 [==============================] - 29s 165ms/step - loss: 2.4519 - val_loss: 2.5866\n",
      "Epoch 167/200\n",
      "175/175 [==============================] - 28s 161ms/step - loss: 2.4513 - val_loss: 2.5855\n",
      "Epoch 168/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4512 - val_loss: 2.5856\n",
      "Epoch 169/200\n",
      "175/175 [==============================] - 29s 164ms/step - loss: 2.4513 - val_loss: 2.5857\n",
      "Epoch 00169: early stopping\n"
     ]
    }
   ],
   "source": [
    "# Best parameters [79, 1, 0.10521284879519917]\n",
    "# Validation Loss 2.3823776245117188\n",
    "\n",
    "num_filters = 79\n",
    "kernel_length = 1\n",
    "prob = 0.10521284879519917\n",
    "num_links = X_train.shape[2]\n",
    "lags = X_train.shape[1]\n",
    "opt = tf.keras.optimizers.Adagrad(learning_rate=1e-3) #1e-3\n",
    "loss = lambda y, f: joint_tilted_loss(quantiles, y, f)\n",
    "net = joint_convLstm(num_filters, kernel_length, lags, num_links, preds, quantiles, prob, loss, opt)\n",
    "mc = tf.keras.callbacks.ModelCheckpoint( filepath=\"Model Checkpoint/joint\",\n",
    "                                        save_weights_only=True,\n",
    "                                        monitor='val_loss',\n",
    "                                        mode='min',\n",
    "                                        save_best_only=True)\n",
    "es = EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=2)\n",
    "net.fit(X_train, y_traink, validation_data=(X_val, y_valk), batch_size=50, callbacks=[es, mc], epochs=200)\n",
    "net.load_weights(\"Model Checkpoint/joint\")\n",
    "y_pred = net.predict(X_test)\n",
    "## Sort only the quantiles\n",
    "y_pred[:,:,:,:,1:] = np.sort(y_pred[:,:,:,:,1:], axis=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('DQR/y_test_pred_joint_200s.npy', y_pred)\n",
    "np.save('DQR/y_test_200s.npy', y_test)\n",
    "np.save('DQR/y_test_mean_200s.npy', y_mean_test)\n",
    "np.save('DQR/y_test_std_200s.npy', y_std_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "colab_type": "code",
    "id": "NTlKyyKHcswk",
    "outputId": "a47fa6d6-7302-47f4-8d21-1162d68d6750"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6c1468b0b8>]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAD8CAYAAABdCyJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsnXecXFXZx79ntqSQvgmQtptIUKlSYsQXK4hAULCAgIMgxUiimFflVXEVX9B9LagQVEQQMJBRwEYNiFLEQgtIkSYBkiUkpGzqZvvM8/5x7pl77p07M3fq7mzu7/PZz945c+aeW845T38eJSJEiBAhQoQIxSI22BcQIUKECBFqGxEhiRAhQoQIJSEiJBEiRIgQoSREhCRChAgRIpSEiJBEiBAhQoSSEBGSCBEiRIhQEiJCEiFChAgRSkJoQqKUqlNK/UspdYfzebZS6hGl1Eql1E1KqUanfYTzeaXz/SzrHBc47S8qpY4u981EiBAhQoTqoxCJZDHwvPX5+8ClIjIH2AKc7bSfDWxx2i91+qGU2hc4BdgPOAa4QilVV9rlR4gQIUKEwYYKE9mulJoBLAXagC8BHwY2AnuKyIBS6p3A/4rI0UqpPznHDyml6oE3gCnA1wBE5LvOOdP9so07efJkmTVrVin3FyFChAi7HB5//PFNIjKlWuPVh+x3GfAVYKzzuQnYKiIDzuc1wHTneDrwGoBDZLY5/acDD1vntH8TiFmzZrFixYqQlxghQoQIEQCUUqurOV5e1ZZS6kPABhF5vArXg1JqgVJqhVJqxcaNG6sxZIQIESJEKAFhbCSHA8crpVYBNwJHAEuACY7qCmAG8Lpz/DowE8D5fjzQYbcH/CYNEblKROaKyNwpU6ommUWIECFChCKRl5CIyAUiMkNEZqGN5feJSBy4HzjR6XYGcKtzfJvzGef7+0QbYm4DTnG8umYDewOPlu1OIkSIECHCoCCsjSQIXwVuVEp9B/gXcI3Tfg1wg1JqJbAZTXwQkWeVUjcDzwEDwOdEJFnC+BEiRIgQYQgglNfWYGHu3LkSGdsjRIgQoTAopR4XkbnVGi+KbK8iEgmYNQtiMf0/kRjsK4oQIUKE0lGKaitCAUgkYMEC6OrSn1ev1p8B4vHBu64IESJEKBWRRFIltLa6RMSgq0u3R4gQIUItIyIkVUJ7e2HtESJEiFAriAhJldDcXFh7hAgRItQKIkJSJbS1wejR3rbRo3V7hAgRItQyIkJSJcTjcNVV0NICSun/V10VGdojRIhQ+4i8tqqIeDwiHBEiRBh+iCSSCBEiRIhQEiJCEiFChAgRSkJESCJEiBAhQkmICEkVEaVIiRAhwnBEZGyvEqIUKREiRBiuiCSSKiFKkRIhQoThioiQVAlRipQIESIMV0SEpEqIUqREiBBhuCIvIVFKjVRKPaqUekop9axS6iKn/Uil1BNKqSeVUn9XSs1x2kcopW5SSq1USj2ilJplnesCp/1FpdTRlbqpoYgoRUqECBGGK8JIJL3AESLyNuAg4Bil1GHAz4G4iBwE/Br4htP/bGCLiMwBLgW+D6CU2hdddnc/4BjgCqVUXTlvZigjSpESIUKE4Yq8hEQ0Op2PDc6fOH/jnPbxwFrn+ARgqXP8O+BIpZRy2m8UkV4ReRVYCcwry13UCOJxWLUKUin9vxpEJHI5jhAhQqURyv3XkRweB+YAPxORR5RS5wDLlVLdwHbgMKf7dOA1ABEZUEptA5qc9oet065x2iJUCJHLcYQIEaqBUMZ2EUk6KqwZwDyl1P7AF4H5IjIDuA74cTkuSCm1QCm1Qim1YuPGjeU45S6LbC7HixcPzvVEiBBheKIgry0R2QrcDxwLvE1EHnG+ugn4L+f4dWAmgFKqHq326rDbHcxw2vxjXCUic0Vk7pQpUwq5vAg+ZHMt7uiIVFwRIkQoH8J4bU1RSk1wjkcBRwHPA+OVUm92upk2gNuAM5zjE4H7RESc9lMcr67ZwN7Ao2W7kwgZyOVaHAVCRogQoVwII5FMBe5XSj0NPAb8WUTuAD4D/F4p9RTwKeB/nP7XAE1KqZXAl4CvAYjIs8DNwHPA3cDnRCRZzpuJ4EUu1+IoEDLCUEfkKFI7UFpYGJqYO3eurFixYrAvo6YxebJWZfnR0qI9xyJEGIrwO4qAjruKXObDQSn1uIjMrdZ4UWT7MMeSJVEgZITaQ5SbrrYQEZJhjigQMkItIspNV1uI0sjvAohqxUeoNTQ367inoPYIQw+RRBIhQoQhhyg3XW0hIiQRIkQYcohUsrWFSLUVIUKEIYlIJVs7iCSSCBEiRIhQEiJCEiFChAgRSkJESCJEiBAhQkmICEmECBEiRCgJESGJECFChAglISIkESJEiBChJESEJEKECEMXr74Kf/rTYF9FhDyI4kgiRIgwdPGWt0B/PwzhLOURIokkQoQIQxn9/YN9BRFCICIkESJEiDDY6OuD/faDe+4Z7CspCmFK7Y5USj2qlHpKKfWsUuoip10ppdqUUv9RSj2vlPqC1X65UmqlUupppdQh1rnOUEq95PydkW3MCBEiRNil0N4Ozz0HCxcO9pUUhTA2kl7gCBHpVEo1AH9XSt0F7APMBN4qIiml1O5O/2PR9dj3Bt4B/Bx4h1JqEvAtYC4gwONKqdtEZEt5bylChAgRagwxh6dPpQb3OopEXolENDqdjw3OnwALgYtFJOX02+D0OQG43vndw8AEpdRU4Gh0vffNDvH4M3BMeW+nBvDHP+p0ptu2DfaVRIgQYaigrk7/TyYH9zqKRCgbiVKqTin1JLABTQweAfYCTlZKrVBK3aWU2tvpPh14zfr5GqctW7t/rAXOOVds3Lix8Dsa6rj4Yv3/lVcG9zoiRKgl1CinHhqGkNTofYYiJCKSFJGDgBnAPKXU/sAIoMcpMH81cG05LkhErhKRuSIyd8qUKeU45dBCT4/+P3Lk4F5HhAhDHImEezxndtLzebjhD7forXjd60lmzaLm7rUgry0R2Qrcj1ZJrQH+4Hz1R+BA5/h1tO3EYIbTlq1910J3t/4/YkRVhkskYNYsrYKtxQkaYddEIgELFrif17QnWbBgeM7fRAL+5ysKgBgpVq+m5u41jNfWFKXUBOd4FHAU8AJwC/B+p9t7gf84x7cBpzveW4cB20RkHfAn4INKqYlKqYnAB522XQuGkFQhwMosxtWr9XC1OEEj7JpobYWuLvdzHUm6unT7cENrK3R36/2gDm0j6eqCxYsH86oKQxiJZCpwv1LqaeAxtI3kDuB7wMeVUs8A3wXOcfovB14BVqJVXosARGQz8G3nHI+hDfWby3gvtQGj2qqCLtS/GIFhuxgjVBbVlmzb272fzQbrbx8OaG8HhSYkMdx9oaOjdpg+JUM49cDcuXNlxYoVg30Z5UVjo47WfeEFnf6hgojFggUfpWrWppcbV14JN9wA//jHYF/JsIKRbG2mZPToytZQnzXLkaTRKp+JbGYrE2lpgVWrKjPmYGHWLBhYvYY1zGQb45iA69FZ7P0qpR537NdVQRTZXm2YlA9V2Mmbmwtrr3ksXAj//OdgX8Www2BItm1tmlgZGImks7N2uPSwaGtzJRJznwa1IoFFhGSwUAVC4l+MoD+3tVV86AjDCGYzez/3cTBPZLRXAvG4lngMzAbb0QFnnTW8iEk8DpMmZKq2oHaYvoiQVBOPPuoeV4GQmMXY0qLVWS0tlVVHRBieMJvZfRzJExya0V4p2PPU5tT7+oo3RA9VL8aL/jdTIqklpi8iJNXEtVaozbA0UgwRDGG7Xy1iKEi2fpVPR0fh5xjKXowfOcGVSGqR6YsISTXR2OgeV4GQDOWFU1FERLqsiMfhDCvFal2d/lzKJleoZOAnJMVgqHoxJhLwrsNdieSGG7SBvVaICESEpLqwgxCrkFNnqC6cimNgYLCvYFghkYClS93PyaT+XCxDUgyD827+5vnc1FT4uNlsOoNp0DbPYu1aI5FITTJ7uy4hueoqmDevumNWWSIZigunUrAX3j57D9TcQhzKKDdDUsz5buD09HFDAyxZUvi4Q9GL0TwL47UFtcns7bqE5LOfhcceq9pwiQRc+jOXkNx9V3Xdf/diJcrxCKkVT5Cw8KfTWPvaQE1ydUMV5WZIij2fsR1cd11xap/Qtp7+frjttqrY2sw9+721ao3Z23UJiUEVVExmo9uwzSUkP/huquIbnVk4+/IsK9mbr/L9mvIECQs/h1vPQE1ydUMV2RiPSZPKe758DE4qVZrtILQXY1sbnHAC3HVXcQMVAHPPtkRit9cKIkJicl9VEGaj68MlJH29qYpvdGbhvHPKywAcOeqfNeUJEhZ+7q2egcD2CMWhrU2rk/zYsaM4qW8wvcDicU2MchKl1av1/zfeqPj1mGdhE5JaZPZ2WUIiSqde2H1sd8X9yc2G1otrbI+RqspGF4/DL3+uo+k/cExD9YjI44/rNDBVgJ97M4Sk1ri6oYp4HMaNs1v0ptfXV5zUN+Tjm6pYZMo8ixnTXEJy9ZXJofMsQmKXJCSJBHSLrgcykm5Wr4ZPfQoWLarMeGZDsyWSGKnqbXQmLUt9mMrKZcLcubDPPlUZys/h1jNQk1zdUMZmK73qWHakj4tlhkJJBoOEl17RhGThgurUBonH4d6/uITkk6v+r7IDVgC7JCFpbYUeNCEZjVaui+icf5WYNGajswnJ6BGp6m10g0FIqgh/Oo2WaQNDi8MdBrCZnt3ZENg+HJBIwH1/0+skRrJ6sVe2Yf+55yo8WPmxSxKS1atdQjIK10YiUhkDrdnoxje5iuYvfzFVnY2utxdOd1wngxTdwwTxOGmVxIP3DUREpMywmZ7ZvArUpi4/H1pboXdAzyOjIq2K44ZFSH5144ghlb4lDHY5QpJIaL1sECGByhlo43FYcpk7WY58f5Wir21ngmEqkaQRc6ZzFJBYdtiE+eP8oXp2jZaWCg/gxerVkEQTEjua3tjfK4U7bnf3hl5G1FwWil2OkLS2auK/Ge27uCdez4yKiupWEOKxR6eqw3XYIvMwlkgA10gaEZLKwFkc557dXz27xkEH6f977VXyqcKkZamrgwE0w2UkEtNeSfzkci8hgdqqkhim1O5IpdSjSqmnlFLPKqUu8n1/uVKq0/o8Qil1k1JqpVLqEaXULOu7C5z2F5VSR5fzRsLCSBwvootK7cez6e8qLar/8x/uZFHVqs1sR9BHEkmEUmCYkmrmMjNjluhBFTYtSzIZLJFU2oHrjXWZhARqp0piGImkFzhCRN4GHAQc49RiRyk1F5jo6382sEVE5gCXAt93+u4LnALsBxwDXKGUqjCdz4SROLoZBcD+/BvQe1ClRfVfJ9zJYiJZK65/tSWSsWMrONAQQCSRVBaGgFTBLTaNMhGvsGlZWlpcicQmJJXWsE2bGkxIoDYCa/MSEtEwEkeD8ycOEbgE+IrvJycAJsXb74AjlVLKab9RRHpF5FV0TfcqJ7tyg6tMANC+aA+JSjPriQTs3JlJSKDCgXP2AhzuWXEjiaSyGESJZN3ryZJqiIRNy9LWBsS8xvaGhso7FZz3eXdvsL07oTYCa0PZSJRSdUqpJ4ENwJ9F5BHg88BtIrLO13068BqAiAwA24Amu93BGqfNP9YCpdQKpdSKjRs3Fno/eWGCq8xGPgZNI4sNrgqL1lZv9KpNSCpql7ElEuMGPFxhuIHhfp+DhUEgJM8846q2SimFUEhalqTyqrac2OWKYv6x2SWSWnCxDkVIRCQpIgcBM4B5Sqn3ACcBPyn3BYnIVSIyV0TmTpkypdynB7TeMahGciU9M9rbsxOS+fOLO2eomg42IRkMTv13v6veWMaZoLe3emPuSiiTvSIsEgloX525TotRB4dNy9LaCn1JTUhmsQqoPJMJeNapUbtnu8ahiIK8tkRkK3A/8H5gDrBSKbUKGK2UWul0ex2YCaCUqgfGAx12u4MZTlvVUVcXTEgq6ZnR3OwlHvbx8uWFny90TQeLe7zyJ/3V908/6aTqjWXS9Pf0VG/MXQjd3XrN3HxTdTwObSm+1Oy4YdOytLe7NpJPsYxxbCtqvEJx13KXkNiSUKkFxKqFMF5bU5RSE5zjUcBRwOMisqeIzBKRWUCXY1wHuA0w9dROBO4TEXHaT3G8umYDewOPMghIJt2JWS3PjLa27BJJMZM0rPHwD7+3J+jAoPinV414mcJhkURSdiQS0LlNz9lYlTwObSneXyGxGHVPmLQszc0uIQE3HUyl1UtX/CxzbxApjskcDISRSKYC9yulngYeQ9tI7sjR/xqgyZFQvgR8DUBEngVuBp4D7gY+JyJVdP9w0dISzOlU0jMjHoexu5XPRhLWeHjJ991xGtC2g0p7ivk3l6oRr0giqRj0fPFu6pWcR4mEVtkGERKlKqfuaWuDugZXNVGtvG3r38iUSKA2DO0QzmvraRE5WEQOFJH9ReTigD5jrOMeETlJROaIyDwRecX6rk1E9hKRt4hI5ZP9Z0FbGzTUeSdoNSbLyZ/IJCTFjhvWeLj2dXdMQ0igshM0aHOpSpqJiJBUDLZ0UGmPQ6O2tTUH5r9ScO65lVP3xONw0qmuRFKtvG0TxgWrvWvB0A67YGQ76EnxjnkuIalWuod5b7e5jlRJ44Y1Hs6YHkxIKjlB29thI5MB+BMf9LRXAsbp4OF/aULy/c++XHO5ioY6mpuDCUkl5pGttjVjNtDPHrGN3HADXHFF+cc0SCTg1jvcbfH8/6583rZEAjp3BEskxTriVBu7JCEBmNWsF8OEMcnqpXuwPDNuWJoqW7U30I4Chuu3N9D/+bK76E2m40pLX83NrudJB02e9nLDdjowbpNf5Qc1l6toqKOtLdOuWKl5ZDMchpCMoI83UrsT//D28g/owMylrZvdNfO1rySZPLmy82jxYhAJVnsvXVobc3iXJSTVdmX0jFmmceNxVzIxp/NvoB85wR1zDJ1Vkb7a2iCmvNxrpTYdm3v1B3JF5XbLh3gcdhvtfaejRuX6RfGwGQ5/CVr6+iozKO5cssesZ4CODjjzzMpt6HY4ApTu6jwYiAhJNQlJmaPMEwntHpjTe8sa532HdlZF+orHYdIEr267UpuOzb2mAqZzrRgrBxOh4pEgvWbMO+3oqIzU19bmJinwE5KvnF+5YEgzV2yJwES39/dXNoFiNo9O+7qGMnZdQjKYeYPs8YuEbZAMQnryWWO+8HhnVWwHiQRs3+p1ZqjUppOTe6V2jJWDhdDxSEBvT3mCA/PBZnT87zRxQ+UIiZkr2aSDjo7KjNvUlH1M+7qGMnZZQtLerl9cKlmldO5QVkISFEdiw0y+227VY/YwgjF0VryssNmYxMe9QmmbTjau2eROs8fajk5OWV9fG1HBg4mw8UiJBBkSiUElOGazPDKYg1TlGD+jJvartiqNJUugsT5YIhmWke3DBYkE/OsJsyiE1aul+uU0S5SEci1ee/ItuVRPyh2MZTd2pi+jUmWF/Xrmcmw6ubjmeNyNQzRjdqHd2QYGtDqiFoyVg4X2dvgSP+Lf7JfRbmCef1AQL1SGYzZZJvyEpCFWOYnEOLDUqUzVFmjJoVLjXvA1y7tS6edbteJhZcAuSUhaW0GS7mSpI1n1cpqlJvbKtnjr6ryTb91aPeYWJjKO7cScTaBSZYX9emY/ISlm08nHNXd2esf0qyMi763saG6GH3E++/FcRrtBLuagUhzzggX6v5+QHH5YZVXR8Th88pRMiaS+XksOlcJRH3DHPP9LKUSyR98PReyShMSfQNFsPBU3ahlCsu++8PvfewlLgcgWR7J0qXfyTZ+qF/06plJHislsSn9Xifv165nLsemEjeLPJgXViufLYCDoffjfk3nOQc+3UhzzFVfAkUdmEpInn0hVnClQ4mUyoQoZgMvs0Vlt7JKExA6uAneyVNyoZRS/xx8PL78MO3cWfaqwSegWf0Hf5xvsCcAerE9/V4n79euZzabT1FT8pjNpUnC7uX6jcsimeoHK19yuVXgM21nmkZ85MM+3paUEIiICr+fO2bpyZSZT0NeTLIkpCOOh9qe7MyWS/v4qFqCrwbpBuyQhaWvz6lrrSFbHqGUmi9kZSyAkkD8JXSLh1oJex1TAJSTF3G+YRZjWM8f0uLuNTLFsGWzaVNymk0jA9oAYtMZGyw60xJuXyb/5QOVrbg8HZJtHQcxByevluutgxgx4NHve1tWrMyWSOpJFMwVhPNQSCdi2NdhGUlGNRSSR1B7icTjwAPfFzZqZqrhRK5GA7/6fHvMb39tNN+ZyuyrDeAsWwBvr9KLY5KQsmciWoox4hbiJgpUbKWBTLwStrcF1qhoa3OuPx+H666GhPjhTLNTk2qw+sqhaDXNg3unIxjKsl7//Xf9/5pmsXexyDwYxUkUzBYsX57a1mTluj3k7xzOPR4AqFqCLJJKhD8NVP/WU++KeOfPHFSciCxbAVie24rXNmpDccVNpEkku+A2km9FS0JxJW4oy4hXiJrpgAUhKj9vXkyrJ2J2NC9y503vOeBz2mp1dtVXpmtvDAjk2sHjcTXR6wD7J0teL8dnOUc0ymQyWSIphChKJ7HEgZo4FRbYDfIC/VF5jYROSn/+8JPvpYGCXIiQ2V+1Rf1x8MbzySvYflgjDCZkJ2olOlnz1kspJJH7vKZPzSjZvLipuJqzBO8jDp6tLR+AXQ0xycYE2EUsk4JWXg1VbteKLP+jIt0OXs9RuCEJil3swiDnJTgtFLvuGmWNBke0AL7NXxTIzGMb2g0f5CMfll1dmwAphlyIkQVlF0yhBiZ7LdmBzQmbMnWiJZPsblSMkfgNpN6PoYQQT2VJUQsOwaeuzefgkk8W54eYiAGYskyrGbHB+iaRWqswNOoYYIWlrc21tBmNGJsvq+WfGgeDIdtBzuBJu5DZji38/evnl8g1UBexShCQoq2ga9fUUg3y2A5sT8ksks6aUrtrKFfFtG0hTxOhkDGPQQReFusSGTVtvFmOQB1UxbrjxePZAsOZmb6oY16vIu9HdfHNhY+6yGMgTxV3O/HQhCEk8DrNbvOt01IgUn/pU4VU3szFCTU0uk2GyJPglkkoV8srJ2OZ7F0MMYUrtjlRKPaqUekop9axS6iKnPaGUelEp9W+l1LVKqQanXSmlLldKrVRKPa2UOsQ61xlKqZecvzMqd1vBsCeTf7LYJWkLQT7bgU28zJhGIjnz1O6ixjTIRsQWLXKvy3inCYp+GoquvmYMrodPe5XjuDOrwT6b+28xYxosuUzYe+RrnjZDxOznb49lG/k7OqKAxGywn8uB+yXzJmwEdOGwUouHhSAkAE2TvOtyx7ZkKGcPP9raYMKoXtqZyYe4HdBzyB9kqFSwXcagnJ5bORnbGiMkiEjOP0ABY5zjBuAR4DBgvvOdAn4DLHT6zAfuctoPAx5x2icBrzj/JzrHE3ONfeihh0o5sWyZiFIiIPIXjtAHzt9/TV9V1DnN+YL+li0TaWlxP3+Di0VADuAp3XDjjSXdj31u+8++pnk8LAJyDMulnRnyS85Kf9fSUsSgo0frH+fAsmUiOxklAvIQ7/BcW1FjXnedCMhH9/ynKKXPsWyZ/sq+1xUckv5QT1/p4w5zLFumX6d5SJPYJKNHu8/W7jerOZkx0fz9Chn3x+P/VwTksvHfDDyPWTv2OxWQw/hn0e/1jz9+RQRkK+OkpTmVMa5ZTxfzDc+YZ/HLiswje/0ey53e53vmmSWdG1ghefb2cv6FKbUrIuIkoaDB+RMRWW5d9KPADKfPCcD1zlcPAxOUUlOBo9H13jeLyBbgz8AxxZG/4hCPu4yVnwNY93pxet9cxuAFC3SFM6MSSld7G+mo0UrUNWfzp7eZR5dLVwxQn/aLL9oAHcJlOR6HxoZMiaToMR/R7pd/uPDJjFiHbFJmLabirjb80rRJFWSnSzdSr0lyaqMYm4E538ZtWiLp2tYfGMthJG2F8Chv5xn2T1+jjULe60dO0mOOZzurWq/OkKZt+16SGIeywjNmuZ02bHWx2RsubnAqmdeYRBLKRqKUqlNKPQlsQBODR6zvGoBPAXc7TdMBWw+xxmnL1u4fa4FSaoVSasXGjRsLuZdQMB4f/o2mPlZc6oUg24FBVxcsX+5GoMecyfKtix3DfomEJIx/gGsjUVCvCUlZksHlcU80aSbKEdnOyJH6f3emKjBoMUJtpuKuNvybsGEybFVgNpdYKM5mYM7Xj97UG+jPOI/fdrCOqSxG66BKyt1m20HvuCPja9vYniLG6872VKly3HZ2CrM3vP0bx8CcOTUX+BSKkIhIUkQOQksd85RS+1tfXwE8KCJ/K8cFichVIjJXROZOmTKlHKf0IChVNACpZFEclpkM2dDe7kag/++39JjHf8yZ0CVOljA/N/c5ZmyM2XPqiX9ioDzJ4HJwTIkEJAe8EsnmzfCPfxQ5liEkAXp5T6qYHMWBIvffTPg3YZv4+m18QdkC7O/DwvS3CYn/PH67oqDSRcvsayxYQrATZgUkzzJ7gxlzAL1Ozzq9xDVzwQVZ3XnN3nD7bUZVoli5qp6bfj1QvfIWZUBBXlsishW4H0clpZT6FjAF+JLV7XVgpvV5htOWrb2qcKN0hU7H6A1uvEOxXhmxLE/Ss1gNF19XHokkjD+92QDEkUjKJjLnKHna2upu6ofwL75OW2mp63NIJOAuxhEN7vM8Ji0gez1zIrjwb8JBRuVsLrEGhUp6pr8pi9xIX8Z5/MXKBEUSvWam7p7KmVsuJ/JI0WZvmDBOSyTTZugx335widLB976Xv7yic23f+T9F94DWHBTjpj9YCOO1NUUpNcE5HgUcBbyglDoHbfc4VUTsHfE24HTHe+swYJuIrAP+BHxQKTVRKTUR+KDTVnXE49qr5xHewSe4CXA33EI5rERC13MOogl2PihAd1LKpTolSiRtbfCQeiff5OKsfcwGsG1HrGqExJ8jqY1vACWkrs8hkdjYcw93zN9xEhDsmRNBw+9abeeVMpt5VgkePZXnzy9sTHM+I5EYQmKr0/zqSkHROEKvmV/fkMyaEywv7EWaJZ1vPA4LP5ti5EjF0886DF8F1UzGff8jJ+jn293jtWXWSubqMBLJVOB+pdTTwGNoG8kdwJVUnmljAAAgAElEQVTAHsBDSqknlVIXOv2Xoz2yVgJXA4sARGQz8G3nHI8BFzttg4KRIyRQZC6Uw8qWCwpg7FjfZBfRE7hMEgnAYfIwF/OtrN8bArn7HmWWSHJs6nV15VOFAGlC8sLTvTmTRk4Y6x2zlgoDDRZsImtKDNgqo7TqcKbe6B7l7en+IrpsQSEcczyuA0TNujPF1jo74ayz3GJltrpy1GjF/1xQhk3dlkhy5YUX0ZOsrrKEJCgg0UhfVUsWWSaE8dp6WkQOFpEDRWR/EbnYaa8Xkb1E5CDnz7SLiHzO+e4AEVlhnetaEZnj/F1Xsbvq6dFW7lWrsnZpaRZiMVdkLjajaa6XvNlPJkVIoXj7O/Rjb/1aDr/9EMjGqXhUwc4EPevs8hAScU7+12mnZNXhJpOu8dCPoozeDiF56N7u3EkjfaqLVa9KRETyIB533+nDvDOQ+Mbj8OLzmkg/xDs9vy+GY16+3J2XY9mRbu/rc89l1JX77ysce6zimGOdraoU5ivsb43moMKEJCgg0dhmbDVjLTiKDM/I9h074Ljj4M47M74youTKlwTqYowZqx/BjKnFZTQ1L3kc27iGs5iNm7PLX0fj388IyZRi9et6gm7ZnOK002Dy5OL0oNmImIhbp2TP3fUEPero0lVbiQQkRT+v9/JgVh1uNtuNUkUavR1CMjLZ6Wn2u6pmbBTbthUx2K4HZRn4sqqMHCK9lQkZXxVjcA8iJIHn8kvx1ZJIqkBIggISDSEp2U2/yhiehMSnUzfEQyn41KfcpI19/YruXv0I7rytuIymbW0QI8k2JnAW1/FR/pi174N/DVanFZvH52sTfxHY3tLi1in5TcLZXFXpEklrK9T73GqDONL5xwbr0889tzg109Lr9aKfzasZ3xndeiIBr7ziGzeHHSeChTB+5M4mbNL72CiUY540yfImxMscZJyrnOrgLMZ2f5qh558vn2rLXtN+Cd7vVABeQlJLqtlhT0i8ekhvQKKg6OorbYLG4zCBrenPXbhBJX7V1o7tXg8U245QjIrg/zafm9GWwcGYG46VLpGEyQCcSMANSzOf5RFH6PKphWLRIrjvPn0Ph/EIMZJMZw378my6z+LF+h2/aeAlz2//cHN5g7rCFPaqSVgSSdZ7c+ZRfUOMOziOpLN1FMoxm0JlZu7bhCTDOcWMWy4HFWuN/+XufmIxrQ046yxvmqH7/5Kib8AZU6mixzR7j4Ffgrfv1U9IRtaXyU2/ShiehKS+Xk+Cnp7AXFjgBh0Z6eCwd6SK3hymz3Q5OuOFAi7HYTYgfGOWEqXrR1aXyFT5JJJsnGcs5g1g6+7O5Pzuu6+4Z2sXVAKYxlrWMJNncUOZOjqguev5jN9e8r3yqSRyJeesdQIzIO42kM/+dOInYqwd91Z6GFkUx2ycU8zG+WZeIkaSWAyuvTbgXBWSSCZ0r0NEzx2/4DowIHT1OM+krq7oNePfeyay2cMwesoc+4zttRbZXlzK26EOpbRU0t2dNY1IUKCTWUBQ2OL4ZmsKHOFgJFqdZjg1swF1dbljBkkkUJpRLev6SotgqqRFAQ4HdVpmu0kRD5oY1gUY2o3rb6Eclr+40Qh6A/v5de0A618v32LMlpzzU5/yakyKnUODiZ6+mKWwEkB5NrzWVti+OsVmtGF+wbl1sCSZy5clK/xlBgD24Xmek/2z22YqIJFMzxHCFiNFMuXYUOrqih7Tzxi+hwe5lY8EMoy5jO21gOEpkYAmJD09WdW//kAns6kXo2I66ePuBB1Jr4dT83tmlC1KNyycXe7ue2Lcc189Tz1RWsTsC7w1sN08t+bm4JiDGMmiJK5YzEtwbYkPtDqkqQlGoYMV17N7+ruZ08q3GHM5NvhRK77/Bn3SkD62n68hirZ76k03Kf79Qn3Rm6u/zADAgTydnYnyu+KWSSJpIku5RJz5G4uRSEBXXx0/+kGyqDXjv6dxbM9oN3E8fkIyoq62JJJhT0iyzXe/asue2AVveNYEveiCHo9u0++ZYRMSe8xyF18y6pYPzddjfPs7iu3dpUXMtrbCvzgIgB0BRtf2dqcYUUAMSQP9BUtc5vqySSSKFA0NOhZiwggtCf6Bj6W/7+0aKJuaqdBrrwXff4Nt9W5Uov186+oy3VN7+xX3/rV4yTYowHEWq7MHNvqCeL/4hWTxKkRnnT7LvjTSTwPBzhiN9SkaGxULFsCA1BGztBWFjGnutdeJ4h/PNpTSqbSMKhS8tekFRUrV0zw1IiSDjkQCVm8YybJrskskftWWvakXrGKyuaRer+olV7oHe8xyFl8KCnTqHYiVHDHb3u5esx0wZdDcrInhgs9ksukTRvUVLHG1tupHaz+nz+J6qo2hk507dQ6vLy3UEslVLOBMrgVgx9bi8qcFIVdyziDUgu+/wZjZrhRnCMno0V6hwy6QtnmbM39VqmDXdROQaMcZTWVt9sBGR7V1yx1aC7+toz+7HScPbr9Vz6MdjAXcYMi6OksyUDpPXFePVu8lqfMUtsqX6SToXk39oXFsRwTuvde1tXV0aIIybox+HtOmKQ6ZV8+EkSXWe6kyhh0hMSlLOgdGMpLcEkmQmqkoFZOt3/ARkqB0D0GqrY7sknbBCCr25PdPh+xp6LOhuTk7IbGf25LL9PNor5+d/v6eY35UsMQVpE8/10dIQOfwquvThKSbUWxjfPoay6lmsut2Z8utBrXj+28w5c0T08ej6EmrZu14IHseGUaojiQdHW5EelgsX67PBPAGezCVdXR1wemnB5xn/XpQim//aLRzfW6+tULf7WXOvDSExMyfCRO0VNvQ4GjSSKXXqE1IoPAiacuXu+lgjGrLj/5+6OzU17b4vxUtj9wMK1fCunXhBxpkDDtCsnixfjE9jEwbvm20tMCyZXDg/sJHPxbjGxeaRZEq3m/blkh8qUP8qaKzqdPKiWyBTv00eJ5JoWXq29qgoU5fcwMDmM0gI0W8Q1ibv7vIKaYOe/312oI9m4L06TaMgV0Ebr1JbzA9jExnbTUbQKlqJiPh2cS+rk7bZ/woKV3+EED7S71p1WxQin4/IQFvRHqoMSzJdiNTmIT2k0+lfESpvV1Ti85O/rNGX8houjLOFRbr1wZLJJs3e1MdGYYPMgkJFH+v2QiJGRPgGxfG2DLzAPfHNYJhR0jMYt/JbhnBTiZQLx4HUinaX1P8+DL9CBQpOr3dw8MmJNdck/G1Sfew+LyUM0EVKZRngubicAtFUKBTihjPcADTWZsu2FOovTQeh0MOdqWDNzUnWbYMNm3K4XLsRBA/v3Vq7vQmAbDTegfB9tTq7nAlErPRncpvaKCvZDVTkMdWf7/OpWYyCBgGJeNZ1AJsidpihILqZdjp1YstQWs7ZGxjfNY0Kbzxhv7/oQ8xpVmLg35CUsi7nTY1WCJpbs5kvgyz55fiofh7zeZxCJazT49iUe9lurE3e/+hhmFHSAy2My6DA7AnwLatwmMrFFu2u2qmYsR0IMN1J6Yka1CXzenYG2Shzij+c9ufbbWKrZJY5vjunsKNQLg09H40T3cv9OUXB3Km00ApfjfvBwDcn3qvp0sYtYTZyJomBkck24zCnhP1BqhGjUpvdP/DDzm18Q8lq5mybRybN7sZBCoRPFa1+BRr8n3oqN70eIsW6XfU3g4zpmdKJEHZgsOgrc1LSLKuU1M2YNYsjjmujh5GZBCSQrIPLz7PK5GMZUc6g7G/0qZZp0GExJ/6KBfmz3fXYC5CMhWtxhIUqzY4YmCIaqRDBcOOkBijWdAEtSfLxg3CgMQyDN+FiulAeiH2MAKA0exk9WptqzGLP5GApUtdQpIi5uHo7HTe+eCPmAUvh2+nB7clko3sznPsw2xeLV6Pb1O8bJ47DiF5/AnFGedPYTMTA9WMYTi7eBwu+lYwlZ3GWkBLLR8+Um86P/75KCbv7urszlo4ouQNPtsmWUmDeq4AyHKP89cHXEK95Y2e9Hg//7k7/uuvu/PIr9oKjEjPgXgcxo7W7zTnOjWEZNQoli/XWSP8hETbW8LhuPnGLjMV0Ju3yWA8f762kYBXtbWDsYExSmGg13x+QqJI8QO+CmhC0jBeE5KPH9tVMwGuw46QLFmi9dfbGcd43KR9nsm+fj1zBl7gPTxYHvdfZ3N9gkMAVxfa369tNokEnHtmL2fs+Fmai/ZLJIUgSNWys0sx5TMfSX8+44StvIu/eQtbAe00M5PXPIbjguAjJIFcs0NI7lgeo6tLq5tsI6lB6I04S46kXxNP27Xe9uZuiMX45On1JG5y42zf++7S7VBBHluVNqhnC4AsZ3yKIVY9Pe7zPdmpz+NHNmP7mDFZItLz4ITjgyUSzzo1hGTkSG0uCSAkBa1VZx6tQovizegfd3Vpr8nrrtMM2Eh66GFk+vrsFEgQkNU7C8w7NM/uRH7PS8xJf19Xp9WiJn0/QH29YlOXcSzoqpniVsOOkMTjmgvoHzmOKWyinn6amnyT/bHHAJjKG4EeVAVzmr7MqPbC6OjQxOST/b8CYITjuz6abr7CJZzALel+YWEWz+tM87R/sPtWQE+6s689nL/xnnQpU0NI+miknoGiE0Xam/r7Zq3itNO8XPOZZ8Jvb9ILp2OLHrOL0RmEpKCNOIfeL61S6u7WblV2Og0oS6oJT32MYqvzFYgwec1KhV2P3eTO+gI/CeybzdhebPXJgw/S5+tqGM9udFHnqI8aG3W2gFmz4O9/caTYUaNobg4mJGHXaiIBR7zPVW1tY1xanQTu+tu0CU798E7efPAYWlr0mn4rLzDCkqjDjhnkdTiHl9PHyaReN0ayBlCNDWzt9zoW1EKA67AjJKAn9hfO1y41/bcszzSAWsl1/KqtQsV0IL3RGbdTv6je0eHWpvbjSO4F9AYVdlPP583U2gr78hzg2hEMwbRVaoVO0EQC/vJnd8w/bnt/Rp/+fvj6BXrhKMfQHiSRhN2IEwn4zneCJRIPenpc/9x6K/NPmXIWGYcJ2x5iS2OTJ+u/ctkzqqFOs+uxd5NbRPVHXoO2kRRN2ByGZLuzZsw87ex0mZLEL13VVlsb9CgvQxKWGTEhAdu2ueq5HkZyIE/zCxYQc9bD4sX63d13eyePPDeGOXP0mp7O2nSgayEMUL51arAn2qngB/wP/+makU78ep5F1Ie6A1deQqKUGqmUelQp9ZRS6lml1EVO+2yl1CNKqZVKqZuUUo1O+wjn80rn+1nWuS5w2l9USh1dqZsC4OST9f9//UtHrOEu/EUfXw/Ar464nnHj9SO4ntOZPqmbs8/Wm2tBG0IOicTYKgaypDUzG3whZWjzeTPZky7tIpvFyB92gprFmOx3fzvRJ/IbbNmin0dSXELyMV96/bBEZMEC2LZFj/kxfu/tYFs9jUQCZZdIcl2bHVjW0VE+e0YYdVqpxni7HrtR5WSDYYQGqE/P5bfxVPGEzWG+NvYHM18AsT6XkMTjMHXOaCaN7CpYKjSuvbZ6ro9G3s8DLOBqDuRpQL+/1as1UevoHcO997rBhPO5C6UKy0Bh3pV/nc7hJc7hag7kKcDNHn4dZwJ6zE52YwobUc5vh3qAaxiJpBc4QkTeBhwEHOPUYv8+cKmIzAG2AGc7/c8Gtjjtlzr9UErtC5wC7AccA1yhlCowkiEcEgk49N3OKrzoInjXu/jqp9dz49l/ZvVqeAsv0MUovvDQqVzQ6qTDppsbzrqfpUsLN3Dedot+2X5C0tiobTZNTa7k40cxJTWNqqU+liLlEAiD2XXtHpOCISRBEgmEn6D+xZgLps+ESXpMswnt5nCdYb3FbNULkLnZjbQ+d3e7ny2J5JxPl5ZbLN+1ZUOp6oh86rRyGOPb2vQcVUheicTEXOxgLI/wDgD2b/hP8XYicW0kEExIjPRx0636vQ40jmaUFO7J5FcxGUJi4M/fNobONAEx6zaFQqQw4348DrvtlrlmXuLNXM0CnnLSDRlb7nbG0dQEjaMbmMEaZrAGIVYTAa55CYlTOtf4WTY4fwIcAfzOaV8KGEvvCc5nnO+PVFrHcQJwo4j0isir6Jru88pyFxYM57xha4On/fNL53J77wcB4RCe4HEOZUd3PT+6zN3gb1zaW7CBM5GA/70wU7Vl22WWLIHGWHDQRrG2mXhcu8WmGr2b6+Wpz3k+/y8XAV6JxIxpXB/DIEjfC3A153AaN3jajLvu8ccrRo/W0h64RLOzM9yGZ6teIICQrF3Lb690lNuWRHLn3e47rSsht1iYayu1Ty741WngSiBnnFEeY7yJ5M5HSMYpvQV0x8bwKjpjwSkf7S3eTpRDijcwhGTx10aRSMATL4ymrrerYMJpS16gmSmbkIyglz/zAbY663dfnk+rl/wM4OrVhUl/NiOUDYaQ9I8az5IlmmGY0DKBlKqvmeJWoWwkSqk6pdSTwAbgz8DLwFYRMez0GmC6czwdeA3A+X4b0GS3B/zGHmuBUmqFUmrFxo0bC74hwzn7VUkzWQNo7ngiW9jgZIlds9Z9BB0bgzf7XBtCa6vr9WIvijFj3Jcfj8NZp+xM/8Z29S2ppGYqReeAd3MNyoEF0MsI9t1XL6K0iO+4PpYSZX4O13ADp3Mf7+ed/JPGRvjOt/XzmPcOxVVXwYjR+l0YySSsod+/AdiE5C8cCcDchW/X57FsJD/7qbtwS8ktFubaSu0TFn4JJFswaSHEy64N0uu4rmtkbnyjHV5yW2oMNOhN+OB9SwiYc1RbvSP05h3kYjuSHgaoY/3mBhYvhh1Jr7E9bO6rtjbt2utXbRmMposPcC/j2c7BPAG4DJ5rQ3WfSaFELJ8UP4FtDFDHpVftRjwebI8b6ghFSEQkKSIHATPQUkRwLvEyQESuEpG5IjJ3ypQpBf/eLCST38aPEfQyju3poKQRo9xHkG0TzrUh2CkQjETSRAfjVj/t6XfoWxxCMm0am1xvP09JTShQ551KsTPl5SSzTdq1TOPFFzNTPoTdYNva3JTu/+Sd3MFxnu/fzwMcr+7g2mvh4x9zAxIBtnc3pO+1kHH9tiCba57IFkCX35VzF3okko3r3HGKUR3mQyJB3iwIRdenz4J8qjSDQoiXLWUaVQ7ARXyLj/BHxju6e0WKe9AmzU7G0NNfxwB1pZUydiSSL12os0j7s1CAlkjMO+/oCPbaCpP7Kh7Xrr0TxrkSSb9yCYlNxEydkl/x6XTfIIRdN3ZAYhDeOb2dr+9xDfUTxxE/TWXtN9RRkNeWiGwF7gfeCUxQShm2fwakK8W8DswEcL4fD3TY7QG/KRvMQspFSMayg+2Mo6EBerrdlxxUTCafpGBzHH000sUovsXFPM3buOXSV92OO3dqtsi3m336tCRtbZqz8rvS5uV6Uin6Yl6JJEaKkQExG0nqSSa9EolB2A22vl6fv5cR3MoJGd8fd3xdOv0MwMOPxTjzTOgT18unkHGNnWDS+EzVliEkAKd1Xglr16at0zP3dD3kSnLrttHdDQ8+GJh3KwgirhtrOVRqYd5RoVKtLfHtZDe+ycUAXMi3+SMfYysTmcpaz7M2Ndt7GVFaCg+HkMz/qJaEvhhbwo/5oqeLTUhAExLbVdYgzIYej8Ndd+p59NWvKnrFJSS2Wm2S0sSzb9psFi6ELRP3AqDTIrQG+d6JG5CYXbV1647368SUDhNUq9U2w3htTVFKTXCORwFHAc+jCcqJTrczgFud49uczzjf3yci4rSf4nh1zQb2Bh4t140YGDE2n0QSGzeWcePgFd7E5ZwHZBKSMPrJtjZ3oqSIMdraxH/1gw1ux507dXImX6bEVSsHOPPM4I0pL9eTSrFFeUPi6xmgG6+7zz6OK3BdXXASujAbbGurZkBNHiKvKkTjgL0dX3tnk/jjH5VHzegnJGHGjcfhy1/KVG3dxxEAPMZc3fDyy2kj+1mLx6b7laQ6tLFoEbz3vfziKy+HzlxRzoj0fM+qGF26LfEJik1MzujzFG/z1HjZzjhAM03lkEhMOPm7Un/ji1zGrziDi7iQqaxlFN3pd77bbvAmXqGeJPN4xHOq1asL86689jrlmb+7467T8aIJyYNPT+CKK+CiTdrm+FtOyjhdvlQpWorUY17OedwccI7x2x1tvxPcW41sBpVAGIlkKnC/Uupp4DHgzyJyB/BV4EtKqZVoG4jJVngN0OS0fwn4GoCIPAvcDDwH3A18TkTKXk/SiLEpFexuu9/ULdSR4vNfH+dEqCp+xJcB70bnSfCYZzxl6V5tbHrD2jh37tSrweD88wF48olkOutoEHJyPSI8kHy3p8lE6wL0U8+p/JoX2IfRo/WkjNV5JZKwG6xJOW9SbAe5i173k05iMXjXYfq+N27WRDMbISk0IDFZ7455Hj9hL1ZyTdpZEP5222ZmzYJtM/fn/vPvBHSW4rIYLJ95BoDutVvydMxEOewzuZ6VUsXp0u3aIIIiqTKZryls4r08CMCvOZUOh9j0qRIlEqdg1R/v9KZQPoPruZBvs5bpjKEzLZGMHAkP7KbVqUFSSagN1yEkb2zwGtu/y9fTxy3jHZf28VpNTSzGqtjsokrf+rMcP+l4adlIS1y9vVXJZlAphPHaelpEDhaRA0VkfxG52Gl/RUTmicgcETlJRHqd9h7n8xzn+1esc7WJyF4i8hYRuatSNxWPwzXXB0skV647HoBHXxiX5vL8uYMK1W9P29OVSDzteyTToupNv+pi5brd3Ml+ySWw//6k+nLHOeTkRFMpRo+rZznHppv2Iv24+QKXcyOnopTeSK+4Ag5/Tx2NsWTBvvhGkDLcaxAh+XDvb7lSFrD9da0qSI7Ri9EQEjsos6CIaGcD+NFPXS6yh1G8wl5czWe40PFMG0V3motbe+AxALRdNFAeg6WTnnn6nuWp310o7PxpfhSrsjOqFxxC0ivBa8bAJtpjJo0oWSJJoTj/guxjnsjv0/OsowMeqn8PAH+0JCSDUBuuw5DsvofX2G5jQs8bbGECs+bUp9dqf6ou0H6aL1WKrfbOZmtJE5KenqpkM6gUhmVkO5DVcPUmtN3ishuamD/fqQRnERKl4NxzC9t4Pr8oeLIcPq+fCz+zjt7V6ziZm5kz8EKac0ok4NkX63NyOnmlhVSKd783xiv1bwHgGfb3fL0ZLXvbcSV7v7WOKU2pgj1CjJdQLtXWZDpYwNX8knMA2B4b71EzmsUYi2mX6NBwuNeT4pmLP0VdWu1gPPG6uqD1mzHNEZQYkGgYgYce03PkD+v/i4savl3weQrJGJsNS5aUN+eXHaeTIpZVHWxg1sfChbDbhAbYGhyQGgoiJCXGtu7cY07B9dzcuM3tOybAyyvvhusshLPOiZGsCyYks3pf4D+82aNWitVnqoMhPwHXam9XW7GXlR7FIJ3MtLd3UJKDlgvDlpDkE3M3JCdx882aIx830VW9TJoEhx9e2FjHHq0ny5QpXuJ17509vNw9jZfYO91mXBYXLIDu/mBOB/Rmm1daSKXY74AYTVd/jw/E7uNKzvV8vRHX6y3NrcVihRciwQ0i9Ku2ksT4LSfyLPum+85D5zJ7bfsE3vOeTNVWKpVONhAOIvq6R2QSL4AX2IdzuJoz0uFLWhXXJ/V899vFByTaOuu0G6ikuLD/QvviCHKX9WPHjtzXEMbIWu6cX3acjil8lgvPcABHHKElW15+GW65pXipJBVuTFuKtaWIt/JCRt+8G64jkRz1QcXB84IJydt5jGfZD3ClnKY96mis866ZsAS8TrnaCjuJrMEE03bDDYOSHLRcGJaExAQl5kIHTW4RrB5v6dCCDVwOp3PdUu/jnJTS3NQYdnraOzr0JJ3NqxzH8gzuqrERrr8+/waRSqb4yU8V8bNGcG/q/Rlck/GwAausbl1d4cVP8BpmbULyMnvxCX7LP/mvjN/sNm08DzwQbCMx7s6h4EgknhxaPlzDOWyyCKdSetxSAhJtnXW2zAS3cTwSYhn5yxMsWuTW/VIq02Pv05/OTkzKFWNge23l29SP4h7eYCr33ee7rp07s/4mJ5zaPNlUTAb2nLavbw/We/qF2nCtOjmz9gkOwBzHDv6KWzunvR0mTKpj7kHJggh4utSDuNqKIGcGAI4+Gk47bVCSg5YLw5KQmHK7ubDOqUnQ2go7ur02kq4ubYQMvfGYjdlX5vBX5KZmTU6J0WO429NuNp18HCwibNkeS6+PB3gfAOdzCY8wj2c4IN0/nRQyFqO/J1mwi6GZ5IZ7NRuA+f9CQGjRed+cQDIZTEgKEoqMRKLC+dkrpX9iFyUqxmhpq0qyEZIPc0fB51u0SNf6yIWBAa1iBf1+Jk92ic7kyeXx5DHMQRhCYoJtRbxBgPMO6C7uWkSob1A0jMotkdi2OJvo7MH6dBr20BuuvU4n6jr1D/LujG5mHYFDbOvrmTktWRAB96eQTxHjfH5IBwE6zkMPTR/WYjAiDFNC4nelfYk56RoEADN4jfXsSVOTXtxBpUOTyQK4WLu0bA68lwcYPTrTaBpUqyMfF936dUnXgDd4hgNRCD/ifA7jEXosH3yTFPK5/9TR25PKSP0elpiMaNBjGpWDWeg7HLdQGyeePZ66umBje0H14o1EEhKGsI5jRzq7MhRutLRVJdmMpcWcL6w0ZlLJ+N3DOzoKZHSywDAHIxolI22IH4aQmPGvdZikTa/3cNpphRG3RAKuvFLo7Y/R2Z17ItiExCZ002LrWbq0wA3XkkgMIenAuxjXMpV2Z69ISzl1dQXb2vzpfQTFTsYwmQ6mspbD+TuLG67goc9cC9/8ZkHnHooYloTEjzfzEofyePrz68wA4BOf0Is7qHQoFMDFOhP03vtzP85tsUlcdZU2mjZYjFhQ9cB847/WHuwplgvt7XDvfbEMFZgpwBUG06amUHWxdC6iZ9kPpWCHo0bbMdJVL83au4H3vS/T2A6ZFR6zIZGAX1yZoqs3xqxZ4X5jE6m34WYYKNRo2dbm0i+/RFIAXQO8qpdc0tgUNkbMcbcAACAASURBVKTdycFNY+JHMgmf/WzpwWvxOBywb7C94jwuTx/bhATgLsdT0DBBYVXCRuXTuV2P6XeZ92M749JaTZvQBUXC54WzTu/+c4xv/FgTEr/TiNFUgC/T76OFhbwF5fcCPTfXq6m83nI4865byDuvOtObfJTaDEocloQkyE1yc4BIuXQpzJmT6f5rIxQX60gk11ybuSjSAXPART90c+mMsxj4IIkk3/izmnO7FQahuRk6e4IrM+aL1DaTe017CiHG87H9OYFbWMjPEXHrYPf2uOdevRoeeggOmqslknv5APvHnmPhQsdgmwdm09mxXXPMaTtPHiSTXqNlHQNFGS3jcZeJ9ROSfeXfHBj7d/pzA33UZ6k541e9GEL3ddp4hdks51h+y4k8yLvZwB58gpsBTaxyzb+dOzOLihVVE0UEAgiJna/OpP8xMG6r9twNw3jZnmLZiMjn+CmzeJXLOY/TR/6WSZO0QGBfX32qV8fAFHKvzjr97ncVL23Wm0RWuwW6aiIAK1boBfLUUyEG0fCn90k5WXyNFNXWFlyuolaDEoclIVmyRBusz+cSjk8H3CsW8TMO5+/pfl1d8MAD7iYR5J4Xiot1JugbG7yP80ZO5p08xBaHmzvhZJfzsH3Qc9WEzuY2+p2Lg4Mgs8FspP408mFgT+4YKfqTilQKbuOEdBS9MewHSXWr17gb0jOX3RuKiIBXz2zucyvj+Q2n5PxdU5PeuL856acAHDSjo2ijpfFW8xPsf3MAT6VcG1Q3o3iZvTx9lIJlyzJVL0YaewePMJlNTGEj+/A8c1kBkK7cd+65hUlR/f3F1URZs0ZXSLQ36ls4gV/zSXY677efRpqaXCbNqJz8TFA+Ym/n9wpigo7mbn6uPsdqZvHjlsv56i/3Tq8V+/pG0pOuMBj6Xh2uoKs3xnLm8xW+zw/4iqeLXVc9g7lamxkImQ1pw/kMp6TCxFhapTl5cqZzhVEPLl5co0GJIjJk/w499FApFsuWidTViehXlf/PHBzLnem20aP1efLivvtEQD6xxwPSyej0ucx5zuZqEZAvfKYr/ZOWFnfMazgz63U1NWUZs6dHBOT7E/4v7721tLj38d2Gb0oSldFnzJjst9fS4vZ7ln3kJk7K+P1cHhUB2cZY2c4Yz/0fzBNux1/9KszrExERpfRPfsx/y1bGhX6X6Wd200264ZlnQo/px8KF+hS3cHyowedzh3de5Tjv7Rwnj3Fouu801oiAnMNVsnCh7vd/J66QFl4Nfe9B7z4Xli0TeZr95fd8VA5y3tNL7JX+/ZtYKcdzS3o9LFyo38t/8XcRkA9wj2c8pXKvGTOXLmVx+p2aH7+TfwRetz3/zMEvOavge5XbbhMBOZTHPL97ln3kdo4TAfkPczLfn/PhvvPvlJYWfY/2msqJDRv073/yE1m2TD/DYt6jUiHGsgCsEKneXj0sJRLQHEEYzyCfoxUf4w+AVj+E5mIdieQdhynm8DJrmcrzlhfTNZyDQrj86lFpo6StZnmTFY3uR9boWWfMr3wtlrNQVF2dlyOO1cec3GDi6dfZmV1FYKtXjPuvH6uYBcAXuZSprGOc5TO/x1Sr/7ZMX/pssFPXh5W8wHpmu+sARYooR2BgChll89ry404+lD7O9V6uuAI+dGyKuW+PpbeL152SBldfmeKKK+DGa3Zywe/m8meOKvr686lmW1tdNZN5r11WrrZX2IvbnASdXV36eYgEq7ZAf5eLe/Z7itkw6lH/dQfFVwTZFcMGJPrH3Y/nOJXfAHAFi9LtfhX5D38ohaucnDGJxUJncA7CUA9KHLaEJKxO0R+aYAx6CxYUnsJj+d0x3mBPprOWfXk+sKsxStrIZSPJOoEsV0aT4j0I9liJhOvqHGQnybY47GvItqlvYgoxklzL2exkjMeL693vsohWAYTE3nQKtQUlEnDUJ7Xh//Of2FC0jtlsTmEJCejcaw0Nmjjn1OGnUt4XZyz4zrv95Te1nmhvVqa7nM8lHMU9nrRtuZBvA2pvdzd1Mw9zFblqb9cE0ng7vZu/ZfTJlUTRqHzGjc0kJNutOWNft/mNvbHbKiiDvNkDnHUaNJc6GYsixWVWBuJPfMJ7H36VcCFpWYjFik51UgtBicOWkIT1QvIH5ho9bCElNc1k6eoN9zj9E/AdB/WxbFn4qNZEwjW2f/krMU47TWehbrS8N2MxPEZtY+dIe49ksZMELQ6bI8wmkQBZA/OuevQgHQXd0ABPPBHYJwhmAxk9IrxE0tioa0AsWABPrXM8yDZtLNpg6c/HFgb7TFyPUiHsFX5CYo6dDW/9usx3dAlf4R6O9sQBfo6fshcrM5iiMBtQc7NLqI1n1t8CYivs/m1tsK6hhRd4S2CEOeTm1uNxOPMMYeKkGMuWue1GIgm67ngcxrjxtVk9HXMilc+u6G2/+WY4/XT3s78WCoSQgixCElaqaGqqvaDEYUtI8nkhZYOp31wQ95B3gmaivR24/Xb9oa+P+Kx/sGHkTA6ZsSHnBDJxBVs2e8c0G8uyZXofSia9nlFGrHYrvmWPbvffux1xa6f2LuheTzgBjjoKHn64oN/G43DO2cJuY2OexWVvKjbGjtVMQFeXyzVPYWPRBktDRHcG1KPIhme3TOOYvls9bYHjZyMkqRSJBNT7iH1QnZkx7OCnnMe9HIlSehMqZAMy+aAExYu8lUN4nK/xvcC+ZoM3XoevMZPJbArsG6YEAkp5rq+TsTnn/erVcDy30sGkQEKSK4liIgGf/axRbYXb9jo6vEkgxrMtXezLIK8UZMWYhZUqliypvaDEYUtIwiAoONAQEpFwboWJBHz6jOwiczY0N0Ni24f4bWOcl57r4z/vOpPdNq/he6c/l3MCmbiCoKyi/jQcNgxxMP3/xcE0BqgHzLX5EY/DqleFyWzK6TKZ63zPsS/d67cV7p6aSjFypF6Izc36XrJVKOzosO+1jh5GpDedYlQLhoh2q9H5O1u40CkSZSNj/ByEpLUVYj5C8nMWpo8/zu+YyOZ00akmOujv1wS2kA0oHoeJ41010784hCSZqWj8G/zmzdp1tonsHFu+Egj+YJxpLQ20t2dmdUinGwFu53ie5KBA1VY2jj9djGyjWTMFBgE5+AXnspWJnrZ8OdRsiSQeJ69KcrfdaoNw+DFsCUm2lNsmzYRZGP6MqoaQgOaAzjor+0QxE3TjhsJdcefP15LFjr5GGuljtpOVuO0HDTknpu0+CZnEa/Xq4E3aLDLTfx9e4PP8NPDasnJOO3cygr6CCEljoz7fokXwm+XjGSXd1EtfYf7xInT3xTz+9dmglJdL7KMx/U6LNVjG4zpZYyGwE3UaZIyfw0ayerVX/TiKLj5tJaX8HSdxCx9JExIT81EMsdxjd+Gww1Tgmhk9OtiFublZS3xv5iVmEjxozuctkr73f538PZ5Wb8tqyPYbqXsYGWhXnD8/eCg7dgUK0xwEoc5ycc/FvOnBXGM75De2jyxM2B8yyEtIlFIzlVL3K6WeU0o9q5Ra7LQfpJR6WCn1pFJqhVJqntOulFKXK6VWKqWeVkodYp3rDKXUS87fGZW7rczocdCfb7jBy7EZjtPAz+n09WW3twTl08kHk9V3+XItWfTRyAh6aXAmpwwM5JyYtidTtjFtv3SzGM0is4sC/Yjzmctj6TRWedUhTrH5HY3hCclYxxHnyitJx9OYLKgmE3JepFLs6FShPF7MujXMQS8jGEFv0QbLdJRxgbE3QaqwjI0ui0Ty+ApBKS8h+SHnZ5zvPfyNSU6+NkNIzPwoKDpahFlvirFpkyYaYfTzbW3w70adI+qD3MNsn+dhmBIIhnB+9OGv8jZ50vN1V5c7h/2xKatp4c38J+OdZLNr5mO+DAwdb2nxMqJHWGl2IDOqPifxNq6jzrvNx8zkq3EyZJHPPxhdIfEQ53gs8B9gX+Ae4FinfT7wgHV8F9pydRjwiNM+CXjF+T/ROZ6Ya+xS4khEtJ+38ftuatJ/2XzAT97jftnOGLmd46Sevgw/7iCYOIcPc6sIyME8ntMXvLHRHdf8dgnnyQ52S3c6gr9kHc/cU0ODyBTWi4As5Gc5xzSxMMYX/zrO8HT4FEvD+6g/oeMMHlj8B49vfz7/d9P3NK4XAdmLlzx98vrjn322rGFaQT735p7XME1+M+bscD7/Ac/a+P1fxTkiIO/iQc9g65kSeBFXcG5Gc0acw+GHixx5pPt5+3YRkG9P/KGAyOH8Le/NnsjNIiDr2CN930HxCjljot70JpF4vODnc+e3H3efg9o95/rKwGc+I7LnniLiroWwf+dwlX5OtMpyjvG89yCY+XcyvxEBeSvPCej9IFtciP8Z/p6Ppj9M57Xc79XGCy/oTolE+rwNDdnvLW8sTEgw1OJIRGSdiDzhHO9A12ufDgik/fXGQ5rVPQG43rmfh4EJSqmpwNHoMr2bRWQL8GfgmCJoX2iYTJo33KB1mbYXjV9ldfOG9/E8+/Ah7mQlczLOFcTR+aWDfCKz7VVlfttHoyfNvIkMz+Xxct11MGZUuDGN0dNwTRMcY+Gp/JoUitm8GsglJRK6/Ied6vzdh2hO7L0fGsuqVe70X7Yse+6pSZPcsU2ajQk+g2VeI7gIsbrw6ojmZvfdT589glM+0luU3tlWqcRIsYbp/N3n0bSI4DB95YvTgYCo7ywSyfYt+t0GedbdyvHp4zfYI60S3cJERDKv2yCn8Vuk8ORhwPyPuXqY3WUDmzYVYJ+xxgxT9Mu+PONd9g3aONbKnJ2N27fdyMFNV5LLqO13ObZja+xMFHY1VX9pAKXg0AOdtDmWeiRbFQejBq5FFGQjUUrNAg4GHgH+G7hEKfUa8EPgAqfbdOA162drnLZs7RVFIqFd+Pxuvn6V1aRJbgxJC+3g2wiCFmHQBM2Fzk5X92smjD/jqiEk+TZXJeHVae3t7iIzm/h69mAzk9iDDRmTN5HQagX/MzMi/SmfGevpe+aZpDexIJixzQawgrdzuxW4l1ev39fHmEmNGe7RDQ1e4gwBKpXGxqKLL9nXVUcywwV4Ho/we05kMZcBWOl4gt1TM/bqLDaSSROyExJbZbYn6/kBXwVchsLYkfLdjwfJpCfTZWi1WCkKfRGIxUgkYPv2cN2Nyi02NtNlr6FBa16DUu0borB7k36u06aptMou3712O6YY+93bqi0RfZ5spQFSfZqQ/PWfmpC0tmYPlB47tjYN7VAAIVFKjQF+D/y3iGwHFgJfFJGZwBeBa8pxQUqpBY7NZcXGEiKSwYqdyMIBGBdhM5ntTKAm35FB0CI0E3SPyXoAlS0q0IKfM8xGSHJtrq2t0NMTjniB6/s/ejT8no8D8Dz7sIWJvHWPLRmTNxsRMwvomVXuQs5X+2XzZndsO/Hfh7jTc3050d3N2CmjMor+XHcdXHttHp3+iBHQG+ydlg82pxxESEzg3uUsRiHcbkkLQYRExLdRZZFI3voWSY9p4xAezxosaNudwtyPBz09aaJQUNJAHyF5S3MBtUkcG0m27MZ+tLS40sNNy8d6vmuaJKRS3hpbHR1erUM8Dksu08/1/gdUmojkuldbshuHS+1sQmKyF2QrDWBKJ/wqoQlJrlxkNWsfAULpv4AG4E/Al6y2bYByjhWw3Tn+BXCq1e9FtJ3lVOAXVrunX9BfqTaSMHp8u99yjvF8GaS/DNT93qz11PvzTGgdvhnzAto8X36M3+XVlSolshcviaBtHLnGsnXjCxeKKFLSSI+AyMPMk3tiR2fck9FZT2G9nMcSeTMvSB398mmuFQFpZlW6b757NfexbJnI2yasyni+ft39vvtmnuPpluNEDjkk3ce2feXVx8+YoU/S35+jUzCamtxrWMYn0zmoNqK/mEF7xrXO5mURkD/wkfw68EMPFTnuOPdzb68IyPcntAmIHM1dvueVkiWcJwLyJAem2zfSJNsZk/ddZM3bNm6cyOLFIpJ9zQTNxy+f2eHpNJPVOW0x9nu7ebdPS2fTzFDrxbYtiojc2fYvT4c9J/XmnX8iIrJ0qW5cuTLUvdq2m/t4n6fT3rzoyStmvrqUxdJPnXyZS6SJjemcZB/kT7JsWW57UNb3UwSoso0kfwdNJK4HLvO1Pw+8zzk+EnjcOT4Or7H9Uad9EvAq2tA+0TmelGvsUglJPiOeeXGm39180LdwM39jFoq9KBZN1skBj5z6bKiFYX4HIt/lqyIgL7K3CMjJ/CZvssiWFpF38g8RkKO5y3PeXJusf+HcxdHyMPMynonZQD/P5ekvBojJq+gTTGJT+py57jPjPrZs8XTwX18QEQGRv3CEvLTn4SJShCHZdPrtb0PNGRv2/PkNJ8sLvFlAZD53yJe5JOt9P8z/t3fmcVJU1x7/3u6ZYRsQGHaRQYNEQTQqMXGJy8MtqEETNxyXiIJiTNREX1SMcSNPoyZqIj5xxcw8NW7RuO/7kmAiCq6AQERFFhVkGJjpPu+PW7f6VnVVdfdMz9ZTv8+nPjN961bdOnc7955z7jm7yEvsLrUcIwNZkXXf/dYddxQ55JBMgY2NIiAXcKlbjoDMY5x8xmABkd9xrgjIJVzgvvBy/lsEJEFTZHuEGlV06yby619n0Rz17PTpIj1Y78m0E3OzJ28H/na7g2NlsdpSEoncY9SvBB/X/UNPpl6sy0lzba3IL6tuFwHZY9iiyEndPGOPl3/xHU+mydR52tI4iH2J3d08b/Bd2YvnRECOHvxszoVtqTOSPQAB3gbecq6JTvqbwDy0zmRnJ78CrgcWAe8A4613TQEWOteJucpu7R2J6QQm30Mc7MmQpDG0we1BYaxBLj763ZzePf1WVMejO/dx6NXST9XtcvfsryLpqq0VObLiAREylmL5eCr2D5w6Jmd5O4XMoPglV4mAnMkf3JuNJKWCjW5Z9qrdvhKJgO9JpbyZfAirs1fYVZ5i38g2Dd3BmQyXXhpdOQGwy7qbI2QB24ZOPPbv59nT/XE7x4f2AdlhB5FJkzIFNjWJgFy92cUCGWtAM0GDyO85WwTk1/yP+8KzuFoEcnoJDqyjdFrfvPDCvOu3ttZpY5o8mcyiJohh+d97N0fIu2yT11jxv2cwn3ky9mN1JM2GiZ3CDSIgQ1kuPXuG9117F23a9k129GQ6jPs8bd/LMbw0i7M/coY0kZD9eVwE5KpDX8q5sC3Uw28U2pqR5GO19bKIKBHZXkS+41yPOuk7i8gOIvI9EXnTyS8i8jMR+ZaIjBORuda7bhWRUc51W66yWwo7wp0fVVUZWfoox0jrDK5lPmPdPEGnZ0HLX21ZtLHaev5FHXMgTFViexQ2eoO/cBwD+YIX2ROAWX3O5chpfT0CX781yLHHQv+mLwBYyaC83WH4dRFf0s890GYjldL10y2pZfSzmcZc9JmBi7iITVS4ehQT+8VGRQXccUfA9+ShQwpCDzZQ7+gGwnRHORX2776rdQEFwPYxliQVqo/y+7iyXcicwB3s4XNs6OrJQnQko7cWt0zwKnqNDs0OOmWcHS5hS090RRuhFkHGEMHRdwR52vUbMBgjlbRPZzQQrdP097OxY7N1A93YmBWd0I8gS7Nly2AFQ9iJNzmTP2raCDemMGeq6utx3bmsYoA7fqNoranRMWEAjuJu/siZbj7TxqKbivXrdfMZXdrHbEmStFsnL71enlMX2NE9/EahZE+2Q6Yj+JmJMf0zeP55/fdjtmJ73mYpukXDGIkfhpEs/1Txyivhyn3bWsMo6kdUJ1itBjJkcz0x9Pj6c53hww+BcGuQyrRWrl51c9+83WH4J4k19KcfXwZOPmvWwO8uzUxkxiTahFg1E3dNTUbhDZpZmtO+kYrXEA5fRiOX82vqOIa3Gce+POXx7xU22PzpxhrHxZ13wg/CnREGoaZGh1s1hwODHDcmk9nKYv8EeTI3Zz23bBmhVluLPvJabdlM435+DMCT7M8vuJbDuN/jNdf0Wbt6q6p0GwX2EcNcu3VzafYbNfgXKX4/dg0OvZV84zGHBc1E3n03u1jDSKqqsidzG/4Fgmnnf7OT6+SxPCQypR8DWMXX9KHRMXBZsyY3rbNmafP2L6u25myuctObAtzIpNOZfmK8P5gDwP/5vDyQSRt0Bg+/kWjL7U+hV0tFWwa5lLP+LeY0/tfdAkduvflGKmhwRUDDen6Zc/saKoL6zLtdl7/+VUTCg3NdzG9EQMoSqYIO29kBv4xIpA9fBW/vL71Ul8Em6csaOZnZAulAMUm+uouU0kLxb+gZqCPZAa1I/ZQhIiC38lMRkAeTh+Vdjp1nBEs8itJCAxMZkcyDHCJvsmNk29piGzvhWfYOrt9ttxU54ghvgYmEXMoFAiJHcpcImQN0YZetlN+ML930fJTeOw//XGe+/vrIfHZ92WVvz1syksUiaBGbCchlEPbNz7K3vMAPXH1jWD/39zMTaAwyB1y3YmHgsyNYIn1ZI92plzI2SS3HyEK2Cn13Pn3BPGwHwLOv59lTnmMv2Y8nRED+wJkiIPsPfdtTp5ChOe8gWQWAjqYjac+rWIwkF/yd+ARuEwEZyeLAzvJd3pCd+ac8x16yiv5yA6fIavrlNcmEduCVKz2ZpnBLpI7nas5yLXXyjuTowEy0Ri9jTtPbV22tiFx0kQiIIpV13z9h5Ctbn5M8UQR9Ijzo202kxYk87FFcplCe90QxBP+3HMr97o8oBhQEk/dhJnqiGZoraAKcw3Huj3cYKx+rkbI/j8t2vO0t99vfFjnqKG+ByaT8qc/5AiKTqRNBWwhF9SdjeCEgg/kskj4/I96CpSIgr069xVO3kK37Me/r1cubXo62NvtN8rLQ+jPXXzlcajlGXmFXedLRewV9V9j3221rM9qgBdwq+rs/HuJg+TsHyVx28ry7ttarK/Er9/11ZzKaqJH+67Wy3eW5xH/JjuhT/3dxpC7/9+9Fd7QiI2Yk1tVWjMRe5UC2KwX/Zbs0EbR7ipfZLW9GEqhU81k0/ZxrPc9cxvmykXJ5n9FyF0fKbE6W5QwNnLDzQW2tyJh+n4qAnMMVwUziN3rXE0RDr17e9+Vj7VNdLVLGJlnAtq4FlP/bDxvyqgjI/jwudUx2XzRl4EN505b9LWlZTT/5ggFZ3xdVb0bZui9PioAsJdhc1T8B3shU98e9vY6XpvJu7m8P49t6a5HJk72FlpfLO4ecJz175l5xm2sc89wf+/KkPMxE6ck3gfT5maxhQr+oqs0rFGxVVZCLj7Q0kZCZiQsCd/vd2CA3M8VdnQvIXHbS4W199Z1rx2i3rVkgHMldAtlM3f/xZreQTGaYSJC7EnM/CI9d9JoIyE+4J7B+Fg3+vny63X7y/c2XeW84JsdthbZmJCWtI8kXs3xeLoyMO0xHUsl6vqQvLzgK8iGs4Dn2ybu8QDm/T2Pby3KbAjCOd6igkUbKOYJ76MtXoaFJ80FNDSz4TJ9Q26rPKv7GJMZv/hm1tVZ9pFI0hQR0Wr9e628iafKlL1sGTZTzJju7SmP/t08/Wcu7myjjM4YCUE8PJlwd4to1R5kaits4kZ7UU0YjY5kfWLYfM2boWcCEXx7hccygkUxqPYrt5M/WkWx9wFYkGzP9yKPP8utIABIJthuT1gddq7SOZOjmyVBv1gDfkDkgei1ncBCP8kMey6LPxPSw8Sd+DsD81UM4+eTc3mlXrw46QKjYQA8q0huylONjxsB2zOckbuUsxwMAOG7vE15dknFrE+VmxW5b41Hibo5GqewT4xup4H84l6v4FQB78SLrqWTOHP3usMOQqVS4M9EDj9LeGYYNCNbLrP0yxdDhSV77yOfY1O9BtsQQMxIHdnxto9gdy4KsaGVjxugOeiOnsA/P8QgTeYp9uZxz8yonVKnmYyR+D6NJUvyT8dzEVBIIW/AfDyNplsWHU+apiZuYxEP88yeXewdvKljBbGCf5s1l7VNXl5kzmyjzMBL72/fbR6dXDS5noeOKffEhZ3LMcflHKAyq3w30oDsNXM/PmM84BqONGhKJcKMAMwl/6YtBYSOVgjlztPGG8ZxrW23d9sJW4R/qYyR1dbBhU4Irr9DxSA4/TM+ML7+WDLWOA20kcp3DEIxPKBMnxPYGfOyx2Z+wJR/zID/iWSYUatTmwQZ60IMNWYxqwYLMpHsgj3ECtwMwnE/YSLeCo1ba/WyE5b5eJDtvGU2kSLKIb7lpuw9f4vbxqEVEaGA8p9LXrwoOj53a5Lib6eHzQBAzkq4Bu4MaK5A6jnXvmVXSggXQTTVy0KRyEskEB/MI+/MU661VYTKZYTzTp+cZNjNiR1JdDTuOa0ISGWuQkSxxV6LNtvhwZ3ZnUvePxlSKZEX4BB5khRZEq3FFYfI3UebWcda3O99y0illPD1iCmNZwMHzZhY04dTUZC/0G+hOkjRTuBWA0Xzo0hDmAsRMwqaez+bKwPKMmaox606VdXfoTPLv1Vt48oa5SHHd+YhCoeOR1N7hVJjjB8vfPCKZSJFPsr/+FoeRGJ9qxnW9f6cwgJU8yX7050ueZ+9AuvwICgRn0EB3fsYsruLsrLp86D7d1o2UsxId/riMFBvS3TjjjALc3ePtZ/c77n5e53sBOYUkaS74bRk33JXxDdO3PuP6qFmLL4ch3MJJgbd7lFt+y1atgmHDmlFIJ0RbytEKvdpKR2JgFG+bofUVf+cgVwbuykydQ2Ny8cWFn7KOQjotaUsAfBMned83YYKsGL27HD8oY6HzEAe33OKjrCwjKD79dO+9s86STd3DXW8kk/kV4ZfLX890WcHAYFn0I4+IgPyg2xstqlf/txrLug1ofcXx3J5TV2La91JmSBMJMRZrUbqg6mqRGWhrtzX0ldG878nooWOLLUROPNFTR2uplKv4pYDIqcwSARnM56FWTabLGD3Oq3xfBOR8LhMQmcqNslPfRW7+7Xhb7uYIWcII9yXD+CSSLlM/YS7qQVzPDALZdfn007pNUU+vQAAAH89JREFUecE1phCQGzgl6z252tmvR5mfHBfojiaJ9hIgl1yiXc+cfrr+nUh43hVGd9gp83v/nLGu9PeHnj1Fvhy+nchhh2Ue+OgjkfPP1wc/2xDEOpL2Q02NXuF9TV8W8i3XyaDnYJQRqpaXu+cMzALEyMub5cFTKdTgwe7PXqz37mCamhg0tIw5D2eWhIdM7t2imM51ddCQKnNpuvHPm7yrwlSKDRvDdyQm/Kl5V9jK0i9CMKKtdDrg250dyTcbvTu0QmOu26JKyIibujt6L3+EvSinnJv3Xsc6epOMcGVvVrdLl2bKWkdvdwUeSIe1I7FDIZtzSZlzJMlQj7Hi7FKMXmZXXgf0mYmfcC+zOYVLvvq5m/8AnuBI7nE8XGt86nPCbe8q//IXXYY/EJy/fm1nkll12ZjZkXzMlm6yLQI0iGrnICeL61I96Um2YsfU3Vvzyxg5uoLEn68DYMHEc9w8NTVaYuBHRYX3nJld/pnnZERU5rAh6Po64QTo29vrSZlRo6JPRpcIYkbigxkEDXT3eHD9fKmjMHUGxb/eKSeR0IcFzSA38vJC5b4urL325MM2eplEU5PuoAMsJV5vrxfUQmAGZaNkJuzerPN6QE2laJRwRmKU8kED3I7Q6BchGNFWoGihKfvktkEhBgUzZ3p1Cv5Jy38aOkzMUVMDe35nLetVZehkbsRzdXV6vljuTMx9WOtRhGfRYTESOxSyn5FE6akM/PSdxTXcyxEAHn2UESlWso4RLOVbLPQ8N2FCboW3UYrbzMQ2MMiqS2fMJLuVs5oBLEeLe+wYHzbC2jkwzgrBjKRCaZrveaBM90u0yHCXZy/3jE9z4NBmnmGHN2fMgK83ZBjJEEfPBrrfP/ooWS75uwpiRuKDGQRGeQgwkUdooDsHDnube/5PD4o77ixzV4M2Cl05ezBkSOZ/v+YzldJ6FFtIXZk9SeULMyjtCdsE7Kmvh1NO0WVKInhQ2JNI0AAHrbCcNk3L6W1FfCPllNEUrNdpzFht+VGoTNtunyhGkkvHtPIfH7NEqgPv2btGY+X1khP86iF+5J6ituHSYTESo6MTVLMYiX263Q+bdkN3PT35DyNYbCmihw2Dp5/OWZQLW6/YSLkn3QOnTc+7sFzXF6cA0K1H8PQT1s5BDMYwkm40UGadcE+IZiQbGu26U5m+bSEfazHQCySbTntH4n5fzEhiQGZwGEaiSDMd7aOk32cLOP3UzDY9DIWa4oJezT78dGZV9/lSHyMxO5LevTOa5BbsSIx1jT1h2zEX1q+Hl19I0at3MqfvpSh66+v1Sm32bOjlxGRKk6AHDbz+YoCPJGdHUt7dy0iUCoh5HoFccS5sRmLigyuljW38O8otNn7Eh4zOeodS3onH1MNyhjOIFUzlpuxnSGfosBiJERnZOxKzkMmHkdgiIz9sFyIVbKKRMoKk2itW5CzGA1vMZUfqnDHDJ+J0GuKgQ8tZsgR+fLguu2GDBLovCmPqQQymnp70Yj3r6M2r7MYtTOFCLg50L2PgN13PB2a3accPGsQXHM49CIoRLNXfl0qxeFmyIAOCUkDMSHwwg2O7xLvsw/M8wGFu0KAKNrmDMoqRFLpyNqKhr+sz71zyfoO3AzY16R2JUhlnXs308mYPIjskrB1CFODD91P0qEzm9EeU6zOWLYNXXsn4oRzJEgA+nf139t3Xl9lhJD/6SblnkhEpTGzoZ25+55TTmM3NAZY3DQ06oqZdTpVawxcMysrrp9v+vZJBbApwStifNRk6fOa/NTWakZg22YfneJdtQ0VANoLKMrBpH9BnU1YwNYMw0V0UzGp+j90zaf5AUa88n9Er1tXBAw/qhk2QRiSjPsjlfDTIxNw4Hi2nie8ylyncxsVc5IrzwphwWCCqMJjdZpokp/C/gN6RnMAcAHapmMfMmbB+XYpXXkvmFxishBAzkgDU1ED/tDYkn8RDrtfXAazKi5EUaoprREO2TqYi3eAVkRnRFsBuu+m/3wsye8wNexDZzMPekQAkRG/Tc239o5zRgZ5g7TKvRCs8h/EpzzxDNsMEnn4+W3RYXx9+UCyoTBv+yIIjWcpJ3EoZjYzmA8+9dFrvUA4Y9g533bKe7tLApvJenjxBK+dc9QA68mao9190lE2zIxnMCt5nG8C7bA9zdDiQL0LLNBhdHc5IWiKRWflFgJwX3WZ/uyfDSGbMgI2NmmbDMEUyERCjDEeCFP2rGMDgALqjdiRQONO0FyY3MZUUCQay0p0bunfTbbbuqxQbU96KbJG4u5MgZiQB8K8ezI5kFAvdlU4YIxk2rHArKtNJjdKw3jk851lVG9EW6BizDz+sXas2A/Yg6uaIeNbSO+sQZJnKX97rP39lYCZcu8x57ECKBIPRshTPIHMYydLlwRPA6tX5re5mzvSeAbudn7r/r7dW+LdwEh+wDSf5PPT2Yw1PfLY9laccA8DEIyo9Ho7N5GB/S9BZGj8GOZPe0qUEMpLK3oqKpJ6UNuPrLN2HUroOevTIqMvMiv4r+gbWxXA+IeFMrEs+3BR6OG7vvQOTc6KuDhYtyjCShC9E8NrVGUaybFlGDGbvhvMVB5tFjcEaguMHG68UYYykUKZpL0yEBA105zdcxra8B0BqXb22Ygw5xNsccXdnQsxIfDBipgu41E37DvMAOJDHI3ckY8bA8uWFl2liaRul6AoG050Gb4xte0cyejQcdFDhBTkIGkTL2ZwhrOBXXEWls0sZNig3IzH1ZZ8EDhJV2K9Jk2Q5m7M7rwB6UjWxVn52qmHUwRMA5Le6q6nxWlymSbIfT3IqN3hW5MfzFwBuZirbkvF3bpjqwamHANh5r0p3x2GYYpDYwr97q66G7ZnHUdwFZBYlABvqsxlJj14J9twjTXW1zrux+2YehmF2aatXw4YN2uLITHJNlPMk+2XVRRkp9uMpXQ8bN9GQDt6RLFwYmJwTM2aApDNMIWtB4oyZQcPL6dYN7nMOEtaSOWpfqJTW9Kcw0dW4Iasi7xfKNP1WgAZbO1Zv/fiS+vrwcAOdOdZIPsjJSJRSWyilnlNKvauUWqCUOsO693Ol1PtO+u+t9POUUguVUh8opQ6w0g900hYqpfLzKdLGMGKmmVxAN0vU9D7fZkuWuPEFDCNJJPRgFtGn3luC6dzAhVzMM0zwiLmAjI6E6DMb+cA+/2HwCcMBuIpzOJZaAOrX5WYkQRZbQaIKf5kPMondeBXwikQMow5bSUL+q7tNPl3+0+zHjZyatQiYzVQAqsn498iq/8rKYPPTHGKLmTNhvtqeN5zT12da/qbSjSmefjbb1xYIStL0Zh0N3Tbj2msdhXaAqG/GDG99HMCTrHNMjk/hf9mOd0ij+L5zvqScRjakghlJc1fNy5Z5dxd+Eam9+GpogIVsjUKYzzigeZ4ZDDO/hZOYymyO5k7P/aq03vmF9aPmME27/v0xZvqzBghmJJ0+1kgeyGdH0gT8SkTGoGOw/0wpNUYptQ8wCdhBRMaCjvqilBoDHA2MBQ4EZimlkkqpJDoE7w+BMcBkJ2+Hgj2YbAWmscC5kEuADCNJp1su/1yj+yArGMKlXEg9PenBBjcdcEVbQWc2ClXmzZqlD2LZPGK5dSjNiJwa6nMzknwjFs6apc8oGCxkFN3Y5A5Ag6AIgH60dHWX9IleTLAum3n4nWYasUwQoibgmhrdTkbstKcVLTFBmjfnJSkry7Tf+g0JXns5zZpl60ggLPt6M7e9w8r218c1TiS/Wo5lAdvxH7ZgKxYD2mAkTEfS3HodMcLLSP7DCK7hDKodo4pcesXmHOI1YsNv6M3NTHXPsXyDo8v6QjOSsB1JoUzTbwV4F5OZx/bubz8jsQ8pB4lBSw35hNr9TET+5fy/DngP2ByYDlwuIhude0bjNQm4S0Q2isjH6PjsuzjXQhFZLCKbgLucvB0K/sF0NldyEjczl/F8zEh3IrAHRUvln/4y19Gb3qxjxBbWEsgRbTVnVRyEWbMyLrYAj1WSGRS9e+RmJPlGLAR9RsEc/vrMOZS2M28CUM0SbuOnrtfasEmnkNVdWGRf/6l2rdCGncdmvPT6GcneU7f2ihot5JqAq6vxONicwWX0Yw0J0qRJkErBccfpiebLrxTpVNo97LaCwVpkEtIMI0ZkK/l/y8X0LfuGemdSXcxWORlJS1bNM2dCecLLnM/gOpawJddzGlfy30B4mz76aPPKtNv3CQ7gBk5lsrMzMbqosAVJoUwzaIzbhzD9jMQ+pAylb71VkI5EKTUS2BF4AxgN/EAp9YZS6gWl1HedbJuDx9/2J05aWHqHgt+bwdWcza2cRAM92Jb3OJo7eYBDmct4N09LV8hBIXDLSHHFBZY5blMTHy1ORq5MWwJ7cjGeYzduSPHvd5IoFW53n0+MbxtGh/DXldrt/v48yZ84nSVsyU+ZwwSeBcJXkrvumv/q1X/wzP0+h5H8idM5kMdc3dTid4N3JD2o54Wvv8NXX2XLyYNo9YseR43yxje/jN9wDle6jAT0rmXGDNiU0lZbw/kEyIgcU6nwevYr+UdUJ7j+9oyVmWEk3+Z9qlidxUhymd3mQk0NjBq0NvDeaWTiRIcxkub03Zoab/tuoCencQPvsS2Q2VUHMZJksnCmGTTG7T5qTKzDdCRQ2tZbeTMSpVQlcB9wpoisBcqA/mhx1znAX5VquUMZpdQ0pdRcpdTclStX5n6gyDCiiCBspDt3czQ/5gFWO154E4mWyz/9E0GiSi99j9ovI/ZpWN/Esy8VV9xjr47sVXqQvPeGG4KZST4xvgMxYAAMH05/1nA61wNwLv/DWfyBs7mSppBJ5/nn8ybPFeGF4V4O5wkOdFeW5eLEOyfNRPQyeU9eoMExHU6l9BnQKFpPO03vLmzR4zPP6HvXkLFd/jVXUE6Ty0hAT6iJMn2OxOhrDCMxZRVSz0b88y5jGMZnvM+27MWLHtcekNvsNh/0kbUOXZdn3TuZm9iFNwIPQULz++6cOdnpZme9Odrqpbxb9qTeHDPnoAWTHa/oR/ydU7khkpFACVtv5ePZESgHngB+aaU9Duxj/V4EDATOA86z0p8AdnWuJ6x0T76gq629/xpEhbi1r+7dix9rWURE/vY3EZDfHf6mW9Ya+so1/CLwO5rjcdh4cTUvmcWp7v/zGCe3oEPivs4uGY+qeXr7zRvjxsnjPQ51y+3H6rzqvVD4n9+ZfzrhhbX31n6sFgE5Q+molDX8xc3sj04YGN3SQZQ3WRAZxYfyIIfI9rwlDzBJl8kf3fvV1SJfD9la7k5Olj9zmnxNb0nQlJdH3CAv1CbyZyVr5UNGyTuMlRUMlJuZUvw27dFDBORgHsoivJK1kfXiD9ucD6LG6Cr6y1x2EiE8tnqhEUVFssMML2BbEXTkTPvllzIj9NuaU25zQEcLtYs+DXUHcI0v/VTgEuf/0WixlUIr2ecB3YAtgcVAEr2DWeykVTh5xkaV3V6MJFfI0Wa7is8Xr+pwsxN52C1zLZVyNWcFfk9zvsUMRPOS3XhZDOPwF9CSSTwSe+4pHw79gaRQchEX5sVEmjPxhblgd9uTb0TIhBw+hytE0CGU/fHqoyYCO/Z3PtdwlkkZm1wGVVsrIqNHy8ffP0pe6L6//IPxeYUJCJtUq6tFKir0/xU0SJAb/OZM4oG48EIRkJ2Y6778aP5PqljpllVZWbzJNSy0M4jrSl/QIZuD8kQtCKJg1/UithQB+Qs1npf/jnNDv61V5w0LHZGR7AEI8DbwlnNNdJhBLTAf+BfwX9YzM5wdygfAD630icCHzr0ZucpuL0Yi4o17UFWlr6hY0kXFihUiIL/gGunHarmBU0RAruCcoq1wzEA8gMc8q7brON19+Rr6yl0c2aJJPArLdj5UPmewCMhZXO0O8DFjwieJ5kx8OZmTE7viAi4RELmas0TAmXyzyw+LK14IE7GvsjLrPdtsI3LEESI77SQycWJe9IVNqlGTbVGZiIOqKnHjr3xN77zpb86kHrUjOYbM1lDvPLPztHTcgMhyhoqgF3y1HCPf4zW5kakynn8EljlmTPPKbA46HCNpz6s9GUm7Ip2WzxkkAnI9092eaAIemaslO6OwgXghF4mAvM9oSdAkCZpabeK5pN/VbsFTuNkzyKdP936XUs0vPx9RZSNJuYzzBUTu4zCZz5jQST+sDfKdOCMXI2PGiBx+uA56dcIJOWmrrQ3fcdliGP8VFripuait1fHRhrNMBL2bK6Q+mlNelNTgeG6XF9lDBvF51r1ijRsTJG0M83PSWFHRvPKai5iRWFeXZSQich8/9vTELxggxzHHTQqMLlgAamszYg/72ovnRED+zQ6e9GIyEbPbO5T73QK2JBPFr7lih6jyoiYdEFlPD1lLpVzEhfIJw9zomIVMhLlEW/5og4HYbjuRQw/VCrizz24xXW3FSAztJrrodZweSH/RIopKNBNtFhPPs0zTjualvfk6Z7nF7tO5EDMS6+rKjOTNrX7i9sI32TGrU7ZUvGZWkNmdPi1n8gfZi+daZRDY5U7gKbfglq5Qc5WZa5J/iIPdH2uplNO5Lu/JydRPGHMOukIn0HHjpLFCK65P5NbIiS9fo5DWntj8u8fhLJMkjYHl+kPlFqMf58NME4niiaSnT9ff/zT/ldV3w66WLvwKRcxIrKs9Gcn06ZnVTjJZfLFOGMxAu9fakfyTnd0O6ZGntwCFTELFnNjtCX0XXhcBaaDCU96ECcUrr5BVe5LGLOV6ofXj161FrZiD6nXZzpPcDN/mvUimk0sH0hZtaq/Q27If+b/B1HkQI28N68raWpEerJchfJp3nbe6kY6FmJFYV3sxkgkTgjtCazMTe9I7ijtFQBapreSRmf8ueln5Dv5id3773duyQATkK/pkrd6KhZas2oOuKB1JoWUH7Qr2HvaBm8FmakGTcHNpK9ZipJBvaOsVeVsgbJ4YMyZ8AdFlzX/b82oPRhJ1FqDo5yh88A/K7tS3WucLmwBsBlNVVfzBb5e1BUtFQFYwMOs7ioVCV8tR+adPL0w0k6vsoHZVSqQPX0k1H+dkOs3RkVRWFrdN863fiorSYiRR+hnTP/JdPLQGYkZiXe3BSHKtsFoLuQ6zFdvsOJ9JqDW24rZoyyhmb+KkVmPY+ZztsHVOUWcyCkVUXwrTcxVavmFs/kVAW7RlLhqLUYf5oNh6l3zKixo75jvasg78iBmJdbUHI4kajK21I8l3ZWlOKxdr0NgDsK224n6F9Lf4SMrZ6CmzmCLEfBiJLXYJOyXenHoON2gIp7E55Yf1n0Qi04atNbkWsitqjdV4MdsrX+Rinqa+2/q7bMSMxLo62o6ktXQkhazq/IyuWJ2zLbfiuc4+FBPNEb0Uc4XrtxjLR1xYaPntvfr1W21FTbDFRnvQnu/Or613SjZiRmJd7aUjCVphFdOSyI+WWN8Ua9C09YAMm3yKzaw7guiltdHe8vh86ri1VuPtQXsYvR3JoKCtGUkcateHIG+2tbU6lkZroaVu6IvhUbRQd/AtRVgMiubEpohCEF1h6KyeWQuJCdMayFVvVVUtc1MfhfagPWyszJnTOjR2CrQl1yr06ioHEltyQtmIS4r1HYWKYZpbThQ9rVFePqvmsPMgbS2WKBTtLY8P00MV8xBgGNqL9o7eP4hFW12PkYjktr6pqgr3nlpMRtLagzIX02xtE+sw23+/jqQ9JqeW+Bdrr4ktyqCgrUx+O/qk3h6IGUkXZSQGYSvnKAd8xVrFF9PNdxjy2Rm0NnLtvNpDgRulsG4rrwrNQT4WTDHaHm3NSJQus2Ni/PjxMnfu3Pb+jDZFIqGHYKFoaTOedpqOghgEpSCdbtn77XdFobpaR+wrNurqdJjTZcsycc7D5NlhbVDMevCjrCwT3zvoe8LutTdytWdr1lmMcCil3hSR8blzFgexsr2DoX//9il39uzwe8VUXEaFOW0t5X5dHUyb5g1/O22aN9SwjWIrcP0x3IPKjWIUHXkizhW2tq0U/jHaFzEj6UCoq4O1awt/rjkxqP2ImsiKOblHldNalj0zZkB9vTetvl6nB6GYFmz5MrFitGF7IKo9W9PqL0YHQy7ZF7AF8BzwLrAAOMN3/1eAAAOc3wq4DliIjqq4k5X3BOAj5zohV9ldTUfSXCd8xZChhx0QLLY9fkc6QJYr/noxFLj50hulIyn2Ic1iotAzFbFivG1AR1O2A0MNMwB6o0PljpEMk3kCWGoxkonAYw5D+T7whpPeHx2zvT/Qz/m/X1TZXY2R5DqYWFnZeu7t2+qAYEdyadEWiuBCmFhYiOGOrGwvpD3b20y5K6HDMZKsB+BBYD/n/3uBHYAlFiO5EZhs5f/AYUaTgRutdE++oKurMZKoHUkxXX+Hoa1isHQEJ3ttNYEVwsTa29VJc5Fve3ZW+jojOjQjAUYCy4A+wCTgWifdZiQPA3tYzzwDjAfOBi6w0n8DnB1VXldjJM1x8hcjP7SXSMVE08uHibW3q5PWRqnT15HQ1owkb2W7UqoSuA84E2gCzgcuzPf5AsqZppSaq5Sau3LlymK/vkOjpgb69Am+V2zXIe2NfCyZiomaGm1WnE7rv23hyqKuTrvN0OsmDaXghBOCy29vVyetjVKnrysjL0ailCpHM5E6Ebkf+BawJTBPKbUEGA78Syk1BFiO1p0YDHfSwtI9EJHZIjJeRMYPHDiwcIo6OdasCU7vrH6gglCoOW5nRZC1mEj4oqCt/Z21NUqdvi6NXFsWtNL8DuCaiDxLyIi2DsKrbP+Hk94f+BitaO/n/N8/quyuJtoS6Rpy5K5Ao0j7Wot1VJQ6fR0FdLST7UqpPYCXgHcAczTqfBF51MqzBBgvIquUUgr4M3AgUA+cKCJznXxT0CIxgJkicltU2V3xZHtdHUyZAps2ZdIqKuDWW0vHs2h7nBxvD4wcqXdbfrTW6f0YMQza+mR7Wa4MIvIyencRlWek9b8APwvJdytwa2Gf2PXQ1BT9u7NjxIjgCbbUZOUzZ2qRnS3eikU5MUoR8cn2DoYzzshelafTOr1UMGpUYemdFUGxbVrr9H6MGO2J2GljB0OUE7wO3FQFIcxBYTJZeruvGDHaA7HTxhgljzD/TB3Vw22MGDGiETOSDoaqqsLSOyPCHBR2VseFMWJ0dcSMpIPh2mu1lZaNigqdXiqYNq2w9BgxYnRsxIykg6GmRpv62graUjL9BZg1C6ZPz+xAkkn9e9as9v2uGDFiNA+xsj1GjBgxSgyxsj1GjBgxYnQqxIwkRowYMWK0CDEjiREjRowYLULMSGLEiBEjRosQM5IYMWLEiNEidGirLaXUSnQ8+I6AAcCq9v6INkJXobWr0Aldh9auQidE01otIm0W0KlDM5KOBKXU3LY0p2tPdBVauwqd0HVo7Sp0QseiNRZtxYgRI0aMFiFmJDFixIgRo0WIGUn+mN3eH9CG6Cq0dhU6oevQ2lXohA5Ea6wjiREjRowYLUK8I4kRI0aMGC2DiHTKC9gCeA54F1gAnOGk9weeAj5y/vZz0muAt4F3gFeBHax3HQh8ACwEzo0o83HgK+BhX/rpzrMCDIh4PjAfcA7wlnPNB1JA/xKksx/wgPN9/wC2K4E2rXPKmQ/cCpQ76dsArwEbgbNLmM5Jzre9BcwF9ihhWvcGviYzVi8sUTpD56PQd0Xd7MgXMBTYyfm/N/AhMAb4val84FzgCuf/3axG/CHwhvN/ElgEbAVUAPOAMSFlTgAOCWi4HYGRwJIcDZczn/P+Z0uRTuBK4LfO/9sAz5RAm04ElHPdCUx30gcB3wVmks1ISonOSjIi8u2B90uY1r397yxFOqPmo9B35crQWS7gQWA/NIcdajXuBwF5+wHLnf93BZ6w7p0HnBdRTlRnimy4fPIB/wdMLUU6gUeAH1i/FwGDS4FWJ+9ZwExf2kX4GEkp0ml9z3ulSmvUO0uJTl965HxkrpLQkSilRqK58Bvoiekz59bnwOCAR04CHnP+3xz4j3XvEyetzaGU6one1t4Xcn8knZvOecCPAZRSuwDVwPCgjJ2NVqVUOXAcWtxQyHMj6eR0KqUOU0q9j14oTIl4diSdnFZgV6XUPKXUY0qpsSHPjaTz05lzPrJR1jqf1nZQSlWiCT1TRNYqpdx7IiJKKfHl3wfdcHu06Yfmh0OAV0Rkjf9GidB5OXCtUuottGz432j5qwedlNZZwIsi8lK+D5QKnSLyAPCAUmpP4FJgX/9DJULrv9CuR75RSk0E/gZsbT9QInQahM5HfnTqHYnDSe8D6kTkfid5hVJqqHN/KPCFlX974GZgkoisdpKXoxVlBsOB5Uqp7yml3nKuHzXz+55wnr85z0eORssq/e8pCTpFZK2InCgi3wGOBwYCizs7rUqp3zq0/LKA95QcnSLyIrCVUmpAKdLq9N9vnP8fBcptWkuFTguB81Eg8pGfdcQLrSC6A7jGl34lXuXW753/R6AtGXbz5S9DT2ZbklFujY0od29aQUcCbAasAXqVKp1AX6DC+X8qcEdnpxU4GW110yPk/kVkK9tLhk5gFBll+07oiVCVKK1DLFp3AZZZv0uGTude4HwU+q58MnXEC70VFDKmh2+hrRCqgGfQ5nZPkzGjvRn40so713rXRLSVxSJgRkSZLwErgQ1o2eUBTvovnN9NwKfAzSHPh+YDfgrcVcp0ohWJH6IVkPfjWK10clqbnDI8JqHoSecTYC3aRPMToE8J0vlrtLnrW2hzZ7/5bynRerpD6zzgdSwmUEp0Rs1HYVd8sj1GjBgxYrQInVpHEiNGjBgx2h8xI4kRI0aMGC1CzEhixIgRI0aLEDOSGDFixIjRIsSMJEaMGDFitAgxI4kRI0aMGC1CzEhixIgRI0aLEDOSGDFixIjRIvw/DMkGPeh/hcEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "up_to = 1000\n",
    "plt.plot(y_ix_test[:up_to], np.sum(y_test[:up_to,0,:,0,0]*y_std_test[:up_to,0,:]+y_mean_test[:up_to,0,:], axis=1), 'bo')\n",
    "plt.plot(y_ix_test[:up_to], np.sum(y_pred[:up_to,0,:,0,0]*y_std_test[:up_to,0,:]+y_mean_test[:up_to,0,:], axis=1), 'r')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "BFLwKOhBtYV2"
   },
   "source": [
    "# Quantile Sampling (Example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 296
    },
    "colab_type": "code",
    "id": "pcGCqD48E4NW",
    "outputId": "e95114f8-dca7-4660-aa3c-2d1afaa0d787",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.3215702389568721 0.909438556843736\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6d35b06c50>]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmYFNW5x/HvC4qIuCSAG9uowQVNNDoauG4oKogKGlHEccUbNEi8blEjbtFgFsVE0SgoRqKjBAEVBQF3cGdw4QqKEiLblQBiRARBmPf+cWqGZhimG6Zrqnv693mefrqrurr7nVbmN3VOnXPM3REREQFokHQBIiKSOxQKIiJSSaEgIiKVFAoiIlJJoSAiIpUUCiIiUkmhICIilRQKIiJSSaEgIiKVtkq6gM3VvHlzLyoqSroMEZG8Mm3atKXu3iLdcXkXCkVFRZSVlSVdhohIXjGzuZkcp+YjERGppFAQEZFKCgUREamkUBARkUqxhYKZPWxmi83so008b2Z2j5nNNrPpZnZwXLWIiEhm4jxTeAToWsPzJwLtoltf4P4YaxERkQzEFgruPhlYVsMhPYC/e/A2sJOZ7RZXPSIikl6S4xRaAvNTthdE+75IphwRkZi5w+rVsGoVrFwZ7mt6/N134fiK2ymnwKGHxlpiXgxeM7O+hCYm2rRpk3A1IlKw3GH5cliyBJYuXX9btizsr7h9/XX1j1esCO+xpXbfvV6HwkKgdcp2q2jfRtx9KDAUoLi4uBbfqIhINdasgS++gIULN74tWrRhCKxdu+n3adoUdtwRdthh/a1Vq/WPmzaFJk1g223DreJxdfu23RYaN4Zttgm3rbcGs9i/iiRDYSzQ38xGAD8DvnZ3NR2JSPaVl4df8HPmwD//ueH955/D4sUbv6ZxY2jZEnbbDdq1g44doXnz9bcWLdY//sEPYPvtKR3RkAEDYN4n0KYNDPwVlJTU+U9bK7GFgpk9AXQCmpvZAuBmYGsAd38AGA90A2YDK4EL46pFRArE2rUwezbMnAkzZoTbzJnw6aehTb5Cw4bQti3suSf06BH+mm/ZcsPbD36wWX+Zl5ZC376hSwBg7tywDfkVDOa1ad9KQHFxsWtCPBFh1Sr48EMoK4Np0+C99+CTT0JTUIU99oD27WG//eBHPwohsNde0Lp1aI7JoqKiEARVtW0bTkaSZmbT3L043XF50dEsIgXOHT77DKZMgTffDEEwYwasWxee33lnKC6GE08MIbD//rDvvrDddnVW4rx5m7c/VykURCT3lJeHs4DJk0MQvP46/Pvf4blmzcIVON27hyA45JDQ3FMHnbA1adOm+jOFfLtgUqEgIrlh0SKYNGn9bcmSsL9tWzj+eDjyyHDbd9+MA6C0lNDxOy/q+B0YX/v+wIEb9ilAuJBo4MB4Pi8uCgURSYY7vP8+PPUUPPtsODOAcFXPCSeE2zHHhPb/LVDXHb8V71lXIRQXdTSLSN1Zty70CYwZE8Jg7lxo0ACOOAK6doUuXeCgg8K+Wsr1jt+6po5mEckNFWcEjz0GTzwRmokaNQpNQjfdFKZuaJF26eDNVl86fuuaQkFE4jF/fmjDefTRMFZg663h5JOhV69wldAOO8T68fWl47euKRREJHvWrYPnn4f77w/37nD44WH7zDPhhz+ss1LqS8dvXVMoiEjtLVoEw4bB0KGhfWbXXUOP6wUXhMFiCagvHb91TaEgIltu+nS4887QV7B2LXTuDHfdFcYQZHnE8JYoKVEIbC6FgohsHnd46aUQBhMnhlHD/frBpZfC3nsnXZ3UkkJBRDLjDk8/DbfdFq4m2mWX0B7zy1+GyeOkXohzjWYRqQ/cYexYOPhg+PnPw0IxDz4YLva//voaA6G0NIwXaNAg3JeW1lXRsqV0piAi1XOH8ePh5pvDLKR77QXDh8PZZ8NW6X911JeppAuNzhREZGPTpsGxx4ZxBcuWwd/+FqalPu+8jAIBwlU/qZeDQtgeMCCGeiVrFAoist68eXDuuWH20Y8+gvvug1mzwqWlGYZB6lttzn7JDQoFEYFvvw1/wu+zDzz5JPzmN2EFs379tvjS0k2NHNaI4tymUBApZBVXFLVvD7ffDqefHpauvP32sAB9LQwcGEYQp9KI4tynUBApVHPmhMnoTjstzEM0ZUqYtC5Lf8qXlIQBzm3bhuUP2rYN2+pkzm26+kik0Hz/PdxxRxhvsNVWMGgQ/OpXsYxA1oji/KNQECkk06eHTuP334eePeHPf4ZWrZKuSnKImo9ECsGaNfDb34b1jBcuhNGjQ4eyAkGq0JmCSH334Ydw/vnh/uyz4Z57oFmzpKuSHKUzBZH6qrw8zFh66KHw73+Hq4xKSxUIUiOdKYjUR4sWhbODSZPg1FPDXEXNmyddleQBnSmI1DfjxsFPfhIuMX3gARgzRoEgGVMoiNQXa9bAFVeE+Yp23x3KyuDii8MgAZEMqflIpD5YuDCsgfzmm2HMwZ/+BI0bJ12V5CGFgki+e/llOOusMAXpiBHQq1fSFUkeU/ORSL4qL4c//AGOPz70GUydqkCQWtOZgkg++vbbsLbBmDEhCB56CJo2TboqqQcUCiL5Zv586N49TFkxaFDoXFZnsmRJrM1HZtbVzGaZ2Wwzu66a59uY2Stm9r6ZTTezbnHWI5L33n6bVT8+lOUfzqFb+XMU3XMlpY8rECR7YgsFM2sI3AecCLQHeptZ+yqH3QCMdPefAmcBf42rHpG89+ijrDuqE18s344O/hbPc2LluselpUkXJ/VFnGcKhwGz3X2Ou68BRgA9qhzjwA7R4x2B/4uxHpH85A433ADnnce7DTpyqL/Lx6z/+0rrHks2xRkKLYH5KdsLon2pbgHOMbMFwHjgV9W9kZn1NbMyMytbsmRJHLWK5KY1a8JU1wMHwn//N0evnsQyNp67SOseS7YkfUlqb+ARd28FdAMeNbONanL3oe5e7O7FLVq0qPMiRRKxfHkYnfz3v4cFcYYOZfe21S+Eo3WPJVviDIWFQOuU7VbRvlQXASMB3P0toDGgSVpEvvgCjj46DEx7+OHQfGSmdY8ldnGGwlSgnZntYWaNCB3JY6scMw/oDGBm+xFCQe1DUtg++QQ6doTPPoPnnoMLL6x8SuseS9xiG6fg7mvNrD8wEWgIPOzuM8zsVqDM3ccCVwEPmtkVhE7nC9zd46pJJOdNmwZduoS1k197LayUVoXWPZY4xTp4zd3HEzqQU/fdlPJ4JnB4nDWI5I0pU+Ckk8IiOC++CHvtlXRFUoCS7mgWEYAJE8IZQsuW8PrrCgRJjEJBJGmjR4dpK/bdFyZPDsEgkhCFgkiShg8P6yAcdli40kiXXEvCFAoiSRkyJAxM69wZJk6EnXZKuiIRhYJIIoYMgUsuCYPTnn0Wttsu6YpEAIWCSN1LDYRRo2CbbZKuSKSSQkGkLg0dGgLhpJMUCJKTFAoidWXoULj44hAIo0crECQnKRRE6oICQfKEQkEkbsOHh0Do1k2BIDlPoSASp9GjoU8fOO44BYLkBYWCSFwmTIDevaFDB3j6aWjcOOmKRNJSKIjEYfJkOO00OOAAGDdO4xAkbygURLJt6tQwBqGoSCOVJe8oFESy6aOPoGtXaN48TH+tuYwkzygURLJlzhw4/vjQmfzii5rtVPJSrIvsiBSMxYvDeghr1oTFcvbcM+mKRLaIQkGktlasCGMQFi6El16C9u2TrkhkiykURGpjzRo4/XT44INw2WnHjklXJFIrCgWRLVVeHgamTZoEDz8crjgSyXPqaBbZUtdeC6WlMHAgXHhh0tWIZIVCQWRL3HUX3HknXHop/OY3SVcjkjUKBZHN9fjjcNVV0LMn3H03mCVdkUjWKBRENserr4Z1lTt1gkcfhYYNEy5IJLsUCiKZ+uSTMJ/Rj34ETz2lCe6kXlIoiGRi8eIwFqFRIxg/XvMZSb2lS1JF0lm1Cnr0gEWLQvNRUVHSFYnERqEgUpPycjj3XHjnnbBIzmGHJV2RSKwUCiI1ue66EAaDBoX+BJF6Tn0KIpsyZAjccQf06wdXXJF0NSJ1QqEgUp0JE8LAtG7dNBZBCkqsoWBmXc1slpnNNrPrNnHMmWY208xmmNnjcdYjkpHp0+HMM+HHP4YRI2ArtbJK4YgtFMysIXAfcCLQHuhtZu2rHNMO+A1wuLvvD1weVz0iqUpLw0VEDRqE+9LS6InFi+GUU2D77eG558K9SAGJ80+gw4DZ7j4HwMxGAD2AmSnH/AK4z92/AnD3xTHWIwKEAOjbF1auDNtz54btBmvX0HtYzxAMU6Zo5TQpSHGGQktgfsr2AuBnVY7ZG8DM3gAaAre4+4QYaxJhwID1gVBh5UrH+/eHFVPC3EbFxckUJ5KwjJqPzGwXMxtmZs9H2+3N7KIsfP5WQDugE9AbeNDMNhoqamZ9zazMzMqWLFmShY+VQjZv3sb7+vFXzl7xYJjxtHfvui9KJEdk2qfwCDAR2D3a/pT07f8LgdYp262ifakWAGPd/Xt3/1f0vu2qvpG7D3X3YncvbtGiRYYli1SvTZsNt4/hZe7mf3hh21Pgd79LpiiRHJFpKDR395FAOYC7rwXWpXnNVKCdme1hZo2As4CxVY55mnCWgJk1JzQnzcmwJpEtMnAgNGkSHu/JP3mSM/jU9mHZ3Y+FnmeRApbpv4BvzawZ4ABm1gH4uqYXRMHRn3CG8TEw0t1nmNmtZtY9Omwi8KWZzQReAX7t7l9uwc8hkrGSEhg6FPZvvZyxdKdhA+ezO8fS6xc7JF2aSOLM3dMfZHYwMBg4APgIaAH0dPfp8Za3seLiYi8rK6vrj5X6prwcTj01zHg6cSJ07px0RSKxMrNp7p72CoqMrj5y9/fM7GhgH8CAWe7+fS1rFEnOjTfCs8/C4MEKBJEUNYaCmf18E0/tbWa4+5gYahKJ1xNPwO23h8EJl16adDUiOSXdmcIpNTzngEJB8su0adCnDxx5ZDhL0JxGIhuoMRTc/cK6KkQkdosWhcVydtklTIfdqFHSFYnknHTNR+e4+2NmdmV1z7v7XfGUJZJlq1eH9RC++grefBM03kWkWumaj7aL7qubFSz9ZUsiucAdLrkE3n4bRo2CAw9MuiKRnJWu+WhI9PBFd38j9TkzOzy2qkSy6S9/gUcegVtugdNPT7oakZyW6eC1wRnuE8ktEyfC1VeHMLjxxqSrEcl56foUOgL/BbSo0q+wA2FWU5Hc9emn0KtXWCxn+HBNYSGSgXR9Co2AptFxqf0Ky4GecRUlUmv/+Q907x6uMHrmGdhuu/SvEZG0fQqvAa+Z2SPuPreOahKpnXXrwvTX//wnvPwytG2bdEUieSPTRXa2MbOhQFHqa9z92DiKEqmVa6+FCRPCrHdHHpl0NSJ5JdNQeBJ4AHiI9FNmiyRn+HAYNAj694df/CLpakTyTqahsNbd74+1EpHaevvtMJ/RscfCXRpXKbIlMr0c41kz62dmu5nZDytusVYmsjkWLAgjllu1gpEjYeutk65IJC9leqZwfnT/65R9DuyZ3XJEtsCqVWFthBUr4MUXoVmzpCsSyVuZrqewR9yFiGwRd7joInjvvXDp6f77J12RSF7L9EwBMzsAaA80rtjn7n+PoyiRjP3xj+vXRzilppneRSQTGYWCmd0MdCKEwnjgROB1QKEgyXn2Wbj++jAm4brrkq5GpF7ItKO5J9AZWBStsXAgsGNsVYmkM2MGnH02HHwwDBumxXJEsiTTUFjl7uXAWjPbAVgMtI6vLJEafPllmMKiaVN4+mnYdtukKxKpNzLtUygzs52AB4FpwArgrdiqEtmU77+HM88Ml6C+9lq4BFVEsibTq4/6RQ8fMLMJwA7uPj2+skQ24aqrwnxGjzwCHTokXY1IvZNpR/NR1e1z98nZL0lkEx58EAYPDsFw/vnpjxeRzZZp81HqoLXGwGGEZiRNiCd1Y/Jk6NcPunQJl6GKSCwybT7a4AJwM2sN/CWWikSq+te/wsppe+0FI0ZAQ63vJBKXLV2KagGwXzYLEanWN9+EK43WrQvjEnbaKemKROq1TPsUBhPmOoIQJD8F3ourKBEAysvhnHPg44/D+gjt2iVdkUi9l2mfwiesX5P5S+AJd38jnpJEIjfcAGPHwj33wHHHJV2NSEGoMRTMbGvgDuA84PNo9y7AYOANMzvI3T+ItUIpTI8/Dr//fVgfoX//pKsRKRjpzhQGAU2Atu7+DUA0ovlOM7sf6ApoBlXJrnffhT594KijwiWomsJCpM6kC4VuQDt3r+hPwN2Xm9kvgaWEifFEsmfhwrA2wm67wejR0KhR0hWJFJR0Vx+VpwZCBXdfByxx97drerGZdTWzWWY228w2OY2lmZ1uZm5mxZmVLfVSxWI533wT+hKaN0+6IpGCky4UZprZeVV3mtk5wMc1vdDMGgL3Ec4m2gO9zax9NcdtD/wP8E6mRUv+KS2FoiJo0CDcl5ZWOcA9NBlNmxae/PGPE6hSRNI1H10KjDGzPoQRzADFwLbAaWleexgw293nAJjZCKAHMLPKcbcBf2TDUdNSj5SWhv7ilSvD9ty5YRugpCQ66Pe/DwPTfv/7MC5BRBJR45mCuy90958BtxKuPvocuNXdD3P3hWneuyUwP2V7QbSvkpkdDLR293GbWbfkkQED1gdChZUrw34AxowJGyUlcO21dV6fiKyX6TQXLwMvZ/ODzawBcBdwQQbH9gX6ArRp0yabZUgdmDevhv1Tp4YBah06hAnvdKWRSKK2dJqLTCxkw4V4WkX7KmwPHAC8amafAx2AsdV1Nrv7UHcvdvfiFi1axFiyxGFTOd5h93mhqWiXXeCZZ7RYjkgOiDMUpgLtzGwPM2sEnAWMrXjS3b929+buXuTuRcDbQHd3L4uxJknAwIHQpMmG+3bZdjnjGpwc2pHGjYOdd06mOBHZQGyh4O5rgf7ARMKVSiPdfYaZ3Wpm6kksICUlMHQotG0bWof2bLOWaXufxQ/+byaMGgXtN7ooTUQSYtUMQ8hpxcXFXlamk4m8dtllYaTykCHrL0MSkViZ2TR3TzsWLM7mI5GNDR68fvU0BYJIzlEoSN0ZPx4uvxx69NDqaSI5SqEgdWP6dOjVCw48MIxm0+ppIjlJoSDxW7AATjoJdtwxrJ623XZJVyQim5DpIjsiW+Y//4ETT4Svv4YpU6Bly/SvEZHEKBQkPqtXw2mnwaxZ8PzzoelIRHKaQkHiUV4OF1wAr74Kjz0GnTsnXZGIZEB9ChKPa68Ns57+4Q8pU6GKSK5TKEj23XMP3HknXHopXHNN0tWIyGZQKEh2jR4dxiKcdhrcfbdmPRXJMwoFyZ4pU0JTUceOGosgkqcUCpId778PJ58c1tocO1bTYIvkKYWC1N5nn0HXrmFw2qRJ0KxZ0hWJyBZSKEjtLFgAxx8fLkF94YVNr6gjInlB4xRkyy1dCiecAMuWwSuvwD77JF2RiNSSQkG2zDffQLduMGcOTJwIhxySdEUikgUKBdl8330Hp54K770HY8bA0UcnXZGIZIlCQTbP99/DWWfByy/D3/8O3bWyqkh9oo5mydzatXD22fDMM3DvvXDuuUlXJCJZplCQzKxbB+edB6NGwV13hSksRKTeUShIeuXl0KcPPPFEWEbziiuSrkhEYqJQkJqVl8PFF4f+g9tu0wR3IvWcQkE2zR3694eHHoIbbgg3EanXFApSvfLyEAj33x/ODm69NemKRKQO6JJU2di6daHJaNgw+PWvw0I5mgJbpCDoTEE2tHZtWEZz2DC46abQsaxAECkYOlOQ9b7/PqyH8OSTMHAgXH990hWJSB1TKEiwejX06hUGpg0aBFdemXRFIpIAhYLAt99Cz54wYQLcdx/065d0RSKSEIVCofvyy7Bi2rvvhn6EPn2SrkhEEqRQKGTz50OXLmH669Gjw8ynIlLQFAqF6pNPwgI5X38d1kPQ9NciQsyXpJpZVzObZWazzey6ap6/0sxmmtl0M3vJzNrGWU+hKC2FoiJo0CDcl5ZWOeDdd+GII2DNGnjtNQWCiFSKLRTMrCFwH3Ai0B7obWbtqxz2PlDs7j8BRgF/iqueQlFaCn37wty5YZaKuXPDdmUwPPssHHss7LgjvPEGHHRQovWKSG6J80zhMGC2u89x9zXACKBH6gHu/oq7r4w23wZaxVhPQRgwAFau3HDfypUw4HqHu++GHj1gv/1CIOy1VzJFikjOijMUWgLzU7YXRPs25SLg+eqeMLO+ZlZmZmVLlizJYon1z7x5G+9ryFqunncZXH556Ex+7TXYdde6L05Ecl5OTHNhZucAxcAd1T3v7kPdvdjdi1u0aFG3xeWZNm023G7KNzxDD/pzL1x9dVgkp0mTZIoTkZwXZygsBFqnbLeK9m3AzI4DBgDd3X11jPUUhIED1//Ob8NcpnAkXZjIOxc+AHfcEXqfRUQ2Ic5LUqcC7cxsD0IYnAWcnXqAmf0UGAJ0dffFMdZSMEpKwv1zV73CPf8+k21sDZN/PY5j/9gl2cJEJC/EFgruvtbM+gMTgYbAw+4+w8xuBcrcfSyhuagp8KSFmTjnuXv3uGoqCO6ULLmbkqVXw357w9NPc+zeeyddlYjkiVgHr7n7eGB8lX03pTw+Ls7PLzirVoXrTx97LHQoDx8OO+yQdFUikkfUwFxfzJkTBqQ99lhYJW30aAWCiGw2TXNRH4waBRddFDqRx46FU05JuiIRyVM6U8hn330Hl14KZ5wRBqS9/74CQURqRaGQrz77DDp2hL/+Fa66CiZPDhMdiYjUgpqP8o07PPQQXHEFbLNNmMvo5JOTrkpE6gmdKeSTRYtC81DfvvCzn8EHHygQRCSrFAr5YvRoOOAAeOmlMLHdCy9A69bpXycishkUCrlu6VI499ywhnJREbz3Hlx2maarEJFY6DdLrnIPYw722w9GjICbboK33grbIiIxUUdzLvrXv+CXvwzLZHboAA8+GJqORERipjOFXLJ6NfzpT7D//vDmm3DvvfD66woEEakzOlPIFePGhUVwZs8Oq6Pdey+00kJ0IlK3dKaQtFmzoFu3cGlpw4YwYQKlZzxN0RGtaNAg9C1Xrq8sIhIzhUJSFi8OZwYHHBDWSx40CKZPp3RpF/r2hblzQ1/z3LlhWIKCQUTqgkKhri1fDjffDHvtBYMHwwUXwKefwpVXQqNGDBgAK1du+JKVK2HAgESqFZECoz6FuvLdd3D//WG9zC+/DOMObrsN9t13g8Pmzav+5ZvaLyKSTTpTiNs334S1kYuKwtnAwQfD1Knw5JMbBQJAmzbVv82m9ouIZJNCIS7LlsFvfwtt28I114S+g1degUmToLh4ky8bOBCaNNlwX5MmYb+ISNzUfJRts2eHy0mHDYMVK8IEdgMGhAnsMlBSEu4HDAhNRm3ahECo2C8iEieFQja4w4svwj33hPEGDRvCmWfCtdfCT36y2W9XUqIQEJFkKBRq48sv4fHH4YEHYOZM2HlnuPFGuPhi2H33pKsTEdlsCoXNtW5dmL764YfhqadgzRo45BAYPhx69QoL34iI5Cl1NGfCHT78MDT077EHdOkS1jO45JKwLnJZGZx3XtpAKC0NFyFppLKI5CqdKdRkxgz4xz9g5MgwHUWDBnD88XDnndC9OzRunPFblZaGkckVA9MqRiqD+g9EJHeYuyddw2YpLi72srKyeN587Vp4++3QWTx2bOgnMINOnULH8c9/HvoNtkBRUQiCqtq2hc8/r03RIiLpmdk0d9/09fARnSksXRqagp57DiZMCOMLttoKjjgirGnQsyfsumutP0YjlUUkHxReKHz1FUyeHAaSvfIKTJ8e9rdoEWYqPflkOOEE2HHHrH5smzbVnylopLKI5JLC6WgeOTJcJdSsGZx6KgwZEoLgd78LTUaLFoUriM44I+NA2JyOY41UFpF8UDhnCu6w/fZhhtJjjgkjjGtx+ejmdhxrpLKI5AN1NG8hdRyLSD7JtKO5cJqPskwdxyJSH8UaCmbW1cxmmdlsM7uumue3MbN/RM+/Y2ZFcdaTTZriWkTqo9hCwcwaAvcBJwLtgd5m1r7KYRcBX7n7j4A/A3+Mq55sU8exiNRHcZ4pHAbMdvc57r4GGAH0qHJMD2B49HgU0NnMLMaasqakBIYODX0IZuF+6FB1HItIfovz6qOWwPyU7QVA1UUFKo9x97Vm9jXQDFgaY11ZoymuRaS+yYuOZjPra2ZlZla2ZMmSpMsREam34gyFhUDrlO1W0b5qjzGzrYAdgS+rvpG7D3X3YncvbtGiRUzliohInKEwFWhnZnuYWSPgLGBslWPGAudHj3sCL3u+DZwQEalHYutTiPoI+gMTgYbAw+4+w8xuBcrcfSwwDHjUzGYDywjBISIiCYl1mgt3Hw+Mr7LvppTH3wFnxFmDiIhkLu+muTCzJUA1E0yk1Zw8uaopIfp+aqbvJz19RzVL+vtp6+5pO2XzLhS2lJmVZTLvR6HS91MzfT/p6TuqWb58P3lxSaqIiNQNhYKIiFQqpFAYmnQBOU7fT830/aSn76hmefH9FEyfgoiIpFdIZwoiIpJGwYSCmd1hZp+Y2XQze8rMdkq6plxjZmeY2QwzKzeznL9Koq6kWxek0JnZw2a22Mw+SrqWXGRmrc3sFTObGf37+p+ka6pJwYQC8AJwgLv/BPgU+E3C9eSij4CfA5OTLiRXZLguSKF7BOiadBE5bC1wlbu3BzoAl+by/0MFEwruPsnd10abbxMm6JMU7v6xu89Kuo4ck8m6IAXN3ScTpqmRarj7F+7+XvT4G+BjwrIBOalgQqGKPsDzSRcheaG6dUFy9h+05LZoyeGfAu8kW8mmxTr3UV0zsxeBXat5aoC7PxMdM4BwOldal7Xliky+IxHJPjNrCowGLnf35UnXsyn1KhTc/bianjezC4CTgc6FOkV3uu9INpLJuiAiNTKzrQmBUOruY5KupyYF03xkZl2Ba4Du7r4y6Xokb2SyLojIJkXrzg8DPnb3u5KuJ52CCQXgXmB74AUz+8DMHki6oFxjZqeZ2QKgIzDOzCYmXVPSoosTKtYF+RgY6e4zkq0qt5jZE8BbwD5mtsDMLkq6phxzOHAucGz0u+ck3Ty8AAAE/0lEQVQDM+uWdFGbohHNIiJSqZDOFEREJA2FgoiIVFIoiIhIJYWCiIhUUiiIiEglhYKIiFRSKEhioumEu1TZd7mZ3Z/mdSuyWMPlZnZe9LiDmb0TXUf+sZndEu3vZGb/lcXPvMXMrq7le3xuZs2jx29moaYLzOze6HF/M+tT2/eU/KRQkCQ9QRghnOqsaH/szGwrwuSIj0e7hgN93f0g4ABgZLS/E5C1UNhcUZ2b5O7Zru1h4FdZfk/JEwoFSdIo4KRo+oiKGSR3B6aYWVMze8nM3jOz/zWzjaarjv6Cfy5l+95ofivM7BAze83MppnZRDPbrZrPPxZ4L2VK9Z2BLwDcfZ27z4xqugS4IjqDONLMTonOKN43sxfNbJfoM2+JFpx51czmmNllKbUNMLNPzex1YJ+U/b8ws6lm9qGZjTazJtH+R8zsATN7B/iTmTUzs0nRIi0PAZbyHiui+1tTRswuNLO/RfvPMbN3o/1DojUiMLMLo5reJYy6JfrZVwKfm9lhaf77ST2kUJDEuPsy4F3CAjYQzhJGRpMVfgec5u4HA8cAg6I5ZNKKJh8bDPR090MIf/kOrObQw4FpKdt/BmZZWJnvYjNr7O6fAw8Af3b3g9x9CvA60MHdf0pYX+GalPfYF+hCWIfhZjPb2swOiX62g4BuwKEpx49x90Pd/UDCNBqpU0S0Av7L3a8EbgZed/f9gaeANlV/GHe/KTrL6URY3+BeM9sP6AUcHj23DiiJQvK30XdwBGEBoVRlwJHVfGdSz9WrWVIlL1U0IT0T3Vf8UjTgdjM7CignrGGwC7Aog/fch9D880KUIw2JzgCq2I3wixgAd7/VzEqBE4Czgd6EX7BVtQL+Ef1ibQT8K+W5ce6+GlhtZoujmo8EnqqYiNHMUifUO8DMfgfsBDQlzLFU4Ul3Xxc9PoqwKh7uPs7MvqruB4+C8zHgLnefZmb9gUOAqdF3sS2wGPgZ8Kq7L4le9w9g75S3WkwIOCkwCgVJ2jPAn83sYKCJu1f85V4CtAAOcffvzexzoHGV165lw7PdiucNmOHuHdN89qqq7+nu/wTuN7MHgSVm1qya1w0m/NIda2adgFtSnlud8ngd6f+NPQKc6u4fRk1fnVKe+zbNa6tzC7DA3f8WbRsw3N03WH7WzE5N8z6NCd+PFBg1H0mi3H0F8AqhiSe1g3lHYHEUCMcAbat5+VygvZltY2Y7AZ2j/bOAFmbWEUJzkpntX83rPwZ+VLFhZielNFG1I/xS/w/wDWGG3dTaKtZUOD+DH3MycKqZbWtm2wOnpDy3PfBF1ORVkuY9zo7qPBH4QdUDzOwU4DjgspTdLwE9zWzn6JgfmllbwspfR0d9FVsDZ1R5u70Ja3ZLgVEoSC54AjiQDUOhFCg2s/8FzgM+qfoid59PuELoo+j+/Wj/GqAn8Ecz+xD4gOqvHnqe0CxT4VxCn8IHwKNASdR88yxwWkVHM+Gv8SfNbBqwNN0PF63P+w/gw+gzp6Y8fSPhF/Qb1f2MKX4LHGVmMwjNSPOqOeZKQjNbRafyre4+E7gBmGRm04EXgN3c/Yvo53gr+uyPq7zX4dGxUmA0dbYUNDN7CrjG3T9LupZcYWY/Ba5093OTrkXqnkJBCpqZ7QPs4u6Tk64lV5jZ8cBn0ZVXUmAUCiIiUkl9CiIiUkmhICIilRQKIiJSSaEgIiKVFAoiIlLp/wH3gjo3zcaX7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "S=500\n",
    "t = 60\n",
    "ln = 12\n",
    "(mu, sigma), _ = curve_fit(lambda x, mu, sigma: norm.cdf(x, mu, sigma), y_pred[t,0,ln,0,1:].astype(float), quantiles)\n",
    "print(mu, sigma)\n",
    "x = np.linspace(y_pred[t,0,ln,0,1] - 0.3, y_pred[t,0,ln,0,-1] + 0.3, 1000)\n",
    "plt.plot(y_pred[t,0,ln,0,1:], quantiles, 'bo')\n",
    "plt.xlabel('Value (Standardized)')\n",
    "plt.ylabel('Quantile')\n",
    "plt.plot(x, norm.cdf(x, mu, sigma), 'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  4.,   7.,  30.,  75.,  80., 107., 113.,  53.,  26.,   5.]),\n",
       " array([-3.22647092, -2.6865004 , -2.14652987, -1.60655935, -1.06658883,\n",
       "        -0.52661831,  0.01335221,  0.55332274,  1.09329326,  1.63326378,\n",
       "         2.1732343 ]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAADMNJREFUeJzt3X+onYV9x/H3Z6auq2VT68XZRHeFhg4p3SwX53CMrhZmtVQ3OrGMNeuEULCbXQttrDD/GIVIR39sjEKorhmIs6hDqbLpnMXtD8OilfojdganNRLN7Vrbbh104nd/5Ol2ZxPvzXnOyUm+9/0Cued5znPu832IvH18zjlPUlVIkvr6qXkPIEmaLUMvSc0ZeklqztBLUnOGXpKaM/SS1Jyhl6TmDL0kNWfoJam5DfMeAOC0006rxcXFeY8hSceVhx566NtVtbDadsdE6BcXF9m9e/e8x5Ck40qSZ9eynZduJKk5Qy9JzRl6SWrO0EtSc4Zekpoz9JLUnKGXpOYMvSQ1Z+glqblj4puxkn7S4ra75rLfZ7ZfMpf9anY8o5ek5gy9JDVn6CWpOUMvSc0ZeklqztBLUnOGXpKaM/SS1Jyhl6TmDL0kNectEKTXMK/bEEjT5Bm9JDVn6CWpOUMvSc0ZeklqztBLUnOGXpKaM/SS1Jyhl6TmVg19khuTHEjy2Ip1pya5N8lTw89ThvVJ8udJ9ib5RpJ3zHJ4SdLq1nJG/2Xgolet2wbcV1WbgfuGZYD3AJuHf7YCX5zOmJKkSa0a+qp6APjOq1ZfCuwcHu8ELlux/q/roAeBk5OcMa1hJUlHbtJr9KdX1f7h8QvA6cPjjcBzK7bbN6yTJM3J6Ddjq6qAOtLXJdmaZHeS3cvLy2PHkCQdxqShf/HHl2SGnweG9c8DZ67YbtOw7idU1Y6qWqqqpYWFhQnHkCStZtLQ3wlsGR5vAe5Ysf6Dw6dvzge+t+ISjyRpDla9H32Sm4F3Aqcl2QdcB2wHvpLkSuBZ4PJh87uBi4G9wA+BD81gZknSEVg19FX1gcM8deEhti3gqrFDSZKmx2/GSlJzhl6SmjP0ktScoZek5gy9JDVn6CWpOUMvSc0ZeklqztBLUnOGXpKaM/SS1Jyhl6TmDL0kNWfoJak5Qy9JzRl6SWrO0EtSc4Zekpoz9JLUnKGXpOYMvSQ1Z+glqTlDL0nNGXpJas7QS1Jzhl6SmjP0ktScoZek5kaFPskfJ3k8yWNJbk7y+iRnJ9mVZG+SW5KcOK1hJUlHbuLQJ9kI/BGwVFVvA04ArgCuBz5XVW8BvgtcOY1BJUmTGXvpZgPwM0k2AG8A9gPvAm4dnt8JXDZyH5KkETZM+sKqej7JnwHfAv4LuAd4CHipql4eNtsHbDzU65NsBbYCnHXWWZOOoXVicdtd8x5BOm6NuXRzCnApcDbwZuAk4KK1vr6qdlTVUlUtLSwsTDqGJGkVYy7dvBv4t6parqr/Bm4HLgBOHi7lAGwCnh85oyRphDGh/xZwfpI3JAlwIfAEcD/w/mGbLcAd40aUJI0xceirahcH33R9GHh0+F07gE8CH0uyF3gTcMMU5pQkTWjiN2MBquo64LpXrX4aOG/M75UkTY/fjJWk5gy9JDVn6CWpuVHX6LX++MUl6fjjGb0kNWfoJak5Qy9JzRl6SWrO0EtSc4Zekpoz9JLUnKGXpOYMvSQ1Z+glqTlDL0nNGXpJas7QS1Jzhl6SmjP0ktScoZek5gy9JDVn6CWpOUMvSc0ZeklqztBLUnOGXpKaM/SS1Nyo0Cc5OcmtSZ5MsifJryY5Ncm9SZ4afp4yrWElSUdu7Bn9F4C/q6pfBH4J2ANsA+6rqs3AfcOyJGlOJg59kp8Dfh24AaCqflRVLwGXAjuHzXYCl40dUpI0uTFn9GcDy8BfJfl6ki8lOQk4var2D9u8AJw+dkhJ0uTGhH4D8A7gi1V1LvCfvOoyTVUVUId6cZKtSXYn2b28vDxiDEnSaxkT+n3AvqraNSzfysHwv5jkDIDh54FDvbiqdlTVUlUtLSwsjBhDkvRaJg59Vb0APJfkrcOqC4EngDuBLcO6LcAdoyaUJI2yYeTr/xC4KcmJwNPAhzj4H4+vJLkSeBa4fOQ+JEkjjAp9VT0CLB3iqQvH/F5J0vT4zVhJas7QS1Jzhl6SmjP0ktScoZek5gy9JDVn6CWpOUMvSc2N/WaspGYWt901t30/s/2Sue27M8/oJak5Qy9JzRl6SWrO0EtSc4Zekpoz9JLUnKGXpOYMvSQ1Z+glqTlDL0nNGXpJas7QS1Jzhl6SmjP0ktScoZek5gy9JDVn6CWpOUMvSc0ZeklqbnTok5yQ5OtJvjosn51kV5K9SW5JcuL4MSVJk5rGGf3VwJ4Vy9cDn6uqtwDfBa6cwj4kSRMaFfokm4BLgC8NywHeBdw6bLITuGzMPiRJ44w9o/888AnglWH5TcBLVfXysLwP2HioFybZmmR3kt3Ly8sjx5AkHc7EoU/yXuBAVT00yeurakdVLVXV0sLCwqRjSJJWsWHEay8A3pfkYuD1wM8CXwBOTrJhOKvfBDw/fkxJ0qQmPqOvqmuqalNVLQJXAP9YVb8L3A+8f9hsC3DH6CklSRObxefoPwl8LMleDl6zv2EG+5AkrdGYSzf/q6q+BnxtePw0cN40fq8kaTy/GStJzRl6SWrO0EtSc4Zekpoz9JLUnKGXpOYMvSQ1Z+glqTlDL0nNGXpJas7QS1Jzhl6SmjP0ktTcVO5eqaNrcdtd8x5B0nHEM3pJas7QS1Jzhl6SmjP0ktScoZek5gy9JDXnxyslHTPm9dHhZ7ZfMpf9Hi2e0UtSc4Zekpoz9JLUnKGXpOYMvSQ1Z+glqbmJQ5/kzCT3J3kiyeNJrh7Wn5rk3iRPDT9Pmd64kqQjNeaM/mXg41V1DnA+cFWSc4BtwH1VtRm4b1iWJM3JxKGvqv1V9fDw+AfAHmAjcCmwc9hsJ3DZ2CElSZObyjX6JIvAucAu4PSq2j889QJw+jT2IUmazOjQJ3kjcBvw0ar6/srnqqqAOszrtibZnWT38vLy2DEkSYcxKvRJXsfByN9UVbcPq19Mcsbw/BnAgUO9tqp2VNVSVS0tLCyMGUOS9BrGfOomwA3Anqr67Iqn7gS2DI+3AHdMPp4kaawxd6+8APg94NEkjwzrPgVsB76S5ErgWeDycSNKksaYOPRV9c9ADvP0hZP+XknSdPnNWElqztBLUnOGXpKaM/SS1Jyhl6TmDL0kNWfoJak5Qy9JzRl6SWrO0EtSc4Zekpoz9JLUnKGXpOYMvSQ1Z+glqTlDL0nNGXpJas7QS1Jzhl6Smhvzl4Ove4vb7pr3CJK0Ks/oJak5Qy9JzRl6SWrOa/SS1r15vt/2zPZLZr6P4z70viEqSa/NSzeS1Jyhl6TmDL0kNTeT0Ce5KMk3k+xNsm0W+5Akrc3UQ5/kBOAvgfcA5wAfSHLOtPcjSVqbWZzRnwfsraqnq+pHwN8Al85gP5KkNZhF6DcCz61Y3jeskyTNwdw+R59kK7B1WPyPJN88irs/Dfj2Udzf0db9+KD/MXY/Puh/jGs6vlw/ah+/sJaNZhH654EzVyxvGtb9P1W1A9gxg/2vKsnuqlqax76Phu7HB/2PsfvxQf9jPJaObxaXbv4F2Jzk7CQnAlcAd85gP5KkNZj6GX1VvZzkI8DfAycAN1bV49PejyRpbWZyjb6q7gbunsXvnpK5XDI6irofH/Q/xu7HB/2P8Zg5vlTVvGeQJM2Qt0CQpObWbeiT/GmSbyR5JMk9Sd4875mmKclnkjw5HOPfJjl53jNNW5LfSfJ4kleSHBOfbpiGzrcQSXJjkgNJHpv3LLOQ5Mwk9yd5Yvh38+p5zwTrOPTAZ6rq7VX1y8BXgT+Z90BTdi/wtqp6O/CvwDVznmcWHgN+G3hg3oNMyzq4hciXgYvmPcQMvQx8vKrOAc4HrjoW/vzWbeir6vsrFk8CWr1ZUVX3VNXLw+KDHPw+QytVtaeqjuYX7Y6G1rcQqaoHgO/Me45Zqar9VfXw8PgHwB6OgTsDHPd/w9QYST4NfBD4HvAbcx5nlv4AuGXeQ2hNDnULkV+Z0ywaIckicC6wa76TNA99kn8Afv4QT11bVXdU1bXAtUmuAT4CXHdUBxxpteMbtrmWg/87edPRnG1a1nKM0rEmyRuB24CPvurqwVy0Dn1VvXuNm97Ewc/9H1ehX+34kvw+8F7gwjpOP0d7BH+GXazpFiI6diV5HQcjf1NV3T7veWAdX6NPsnnF4qXAk/OaZRaSXAR8AnhfVf1w3vNozbyFyHEsSYAbgD1V9dl5z/Nj6/YLU0luA94KvAI8C3y4qtqcOSXZC/w08O/Dqger6sNzHGnqkvwW8BfAAvAS8EhV/eZ8pxovycXA5/m/W4h8es4jTU2Sm4F3cvDOji8C11XVDXMdaoqS/BrwT8CjHGwLwKeGuwXMzboNvSStF+v20o0krReGXpKaM/SS1Jyhl6TmDL0kNWfoJak5Qy9JzRl6SWrufwAgNdR9iDolsAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(norm.rvs(loc=mu, scale=sigma, size=S))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "accelerator": "GPU",
  "colab": {
   "collapsed_sections": [
    "QC4IBOREtLAa",
    "lbaB6gJHtVdV",
    "BFLwKOhBtYV2"
   ],
   "name": "Movia Models.ipynb",
   "provenance": []
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}