MDN_prediction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "kUhQYDrgiWsH"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as tck\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import math\n",
    "\n",
    "from tensorflow.keras.models import load_model\n",
    "\n",
    "import mdn\n",
    "import joblib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "OHSSb5g_iWsK"
   },
   "outputs": [],
   "source": [
    "input_parameters = [\n",
    "    'Mass', \n",
    "    'Radius',\n",
    "    'H2O_mass_frac', \n",
    "]\n",
    "\n",
    "output_parameters = [\n",
    "    'H2O_radial_frac',\n",
    "    'Mantle_radial_frac',\n",
    "    'Core_radial_frac',\n",
    "    'Core_mass_frac',\n",
    "    'P_CMB', \n",
    "    'T_CMB',\n",
    "]\n",
    "# load trained model and scaler\n",
    "OUTPUT_DIMS = len(output_parameters)\n",
    "N_MIXES = 20\n",
    "model = load_model(\n",
    "    \"RockyExop_MDN_20210412.h5\", \n",
    "    custom_objects={\"MDN\": mdn.MDN(OUTPUT_DIMS, N_MIXES), \"mdn_loss_func\": mdn.get_mixture_loss_func(OUTPUT_DIMS,N_MIXES)},\n",
    "    compile=False\n",
    ")\n",
    "input_scaler = joblib.load(\"MDN_INPUT_scaler20210412.save\")\n",
    "output_scaler = joblib.load(\"MDN_OUTPUT_scaler20210412.save\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DGCfAnYYiWsL"
   },
   "source": [
    "## Predict Earth (for example)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "W0y2fOq0iWsL"
   },
   "outputs": [],
   "source": [
    "# EARTH\n",
    "# input mass and radius for predicting must be earth units\n",
    "input_array_earth = np.array(\n",
    "    [\n",
    "        [1], # mass\n",
    "        [1], # radius\n",
    "        [0.0005], # water mass fraction\n",
    "    ]\n",
    ")\n",
    "scaled_input = input_scaler.transform(input_array_earth.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 998,
     "status": "ok",
     "timestamp": 1617248997141,
     "user": {
      "displayName": "Vector Zhao",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Ghd3FrP7k66YYE7kRBdxDz2xnB_I-RtbyItWFzX=s64",
      "userId": "13898113326268806207"
     },
     "user_tz": -480
    },
    "id": "fIg5F560iWsL",
    "outputId": "4d5338d2-a427-45b3-95f7-e9054abda16f"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 3.03289294e-03  4.25442994e-01  6.33758724e-01  3.85111332e-01\n",
      "   8.76477808e-02  2.56205469e-01  4.70718741e-03  4.19880003e-01\n",
      "   6.38005972e-01  4.00346518e-01  8.59680176e-02  2.51798242e-01\n",
      "   1.24986768e-02  3.21314454e-01  7.47127056e-01  5.77542245e-01\n",
      "   6.96435720e-02  2.23784477e-01  8.11555982e-03  4.38728154e-01\n",
      "   6.17383659e-01  3.65206301e-01  8.61660764e-02  2.52602160e-01\n",
      "   1.05639696e-02  3.81464452e-01  6.78528011e-01  4.63191420e-01\n",
      "   8.18987936e-02  2.46104375e-01  7.25030899e-03  4.53395605e-01\n",
      "   6.00430489e-01  3.39207709e-01  8.59087706e-02  2.51002967e-01\n",
      "   9.46387649e-03  3.75099123e-01  6.80678368e-01  4.69653726e-01\n",
      "   8.20860863e-02  2.47579694e-01  6.38002157e-03  3.98747355e-01\n",
      "   6.62988901e-01  4.36875463e-01  8.07971656e-02  2.43895382e-01\n",
      "   3.01328301e-03  3.91636997e-01  6.67933762e-01  4.49547648e-01\n",
      "   8.38221163e-02  2.45810136e-01  7.49057531e-03  4.03100491e-01\n",
      "   6.56363845e-01  4.24108416e-01  7.85701722e-02  2.42915958e-01\n",
      "   6.12547994e-03  3.71735483e-01  6.89301133e-01  4.83234465e-01\n",
      "   8.27323347e-02  2.47534156e-01 -1.02724731e-02  4.98230428e-01\n",
      "   5.34916043e-01  2.76906878e-01  7.12869614e-02  2.20684618e-01\n",
      "   8.24606419e-03  3.48897398e-01  7.14737833e-01  5.25507331e-01\n",
      "   8.35565925e-02  2.39772782e-01 -2.37859637e-02  3.57452571e-01\n",
      "   6.97689772e-01  5.01353741e-01  7.81995505e-02  2.33624071e-01\n",
      "   3.29214334e-03  4.11195010e-01  6.47837996e-01  4.14613217e-01\n",
      "   9.14572328e-02  2.61547059e-01  3.03068757e-03  3.67722720e-01\n",
      "   6.97700083e-01  4.93314922e-01  7.53150135e-02  2.40133107e-01\n",
      "   7.64885545e-03  3.59394521e-01  7.04508007e-01  5.07576942e-01\n",
      "   7.49671832e-02  2.36628383e-01 -1.57395005e-03  4.82778490e-01\n",
      "   5.70058346e-01  3.08450460e-01  8.20772722e-02  2.47254431e-01\n",
      "   6.07714057e-03  3.55288684e-01  7.12751627e-01  5.23071408e-01\n",
      "   7.59378076e-02  2.37469375e-01  2.35036314e-02  3.21885824e-01\n",
      "   7.45656013e-01  5.81490755e-01  6.57435507e-02  2.15361357e-01\n",
      "   1.62860751e-02  4.13995981e-03  4.29445505e-03  6.34872913e-03\n",
      "   6.74182177e-03  1.31006241e-02  1.55861378e-02  3.35657597e-03\n",
      "   3.59874964e-03  5.20539284e-03  6.62595034e-03  1.24098063e-02\n",
      "   1.88657045e-02  1.06143355e-02  1.12823248e-02  1.82095170e-02\n",
      "   4.71550226e-03  1.08746290e-02  1.67242289e-02  4.45920229e-03\n",
      "   4.82267141e-03  7.70252943e-03  6.51949644e-03  1.33559108e-02\n",
      "   1.56048536e-02  2.88748741e-03  2.69961357e-03  3.95292044e-03\n",
      "   6.42591715e-03  1.16173029e-02  1.66945457e-02  6.62553310e-03\n",
      "   7.27349520e-03  1.08717680e-02  5.62548637e-03  1.22809410e-02\n",
      "   1.57302022e-02  2.70062685e-03  2.65592337e-03  3.92627716e-03\n",
      "   6.10154867e-03  1.09233856e-02  1.53670311e-02  2.66551971e-03\n",
      "   2.69168615e-03  4.04274464e-03  6.75004721e-03  1.23233795e-02\n",
      "   1.48314238e-02  2.68512964e-03  2.61831284e-03  3.90017033e-03\n",
      "   6.53380156e-03  1.15447044e-02  1.59491897e-02  3.03834677e-03\n",
      "   2.90286541e-03  4.46075201e-03  6.77651167e-03  1.20700002e-02\n",
      "   1.54197812e-02  2.74080038e-03  2.80314684e-03  4.24164534e-03\n",
      "   6.25663996e-03  1.17098093e-02  2.48951912e-02  6.64168000e-02\n",
      "   7.11010098e-02  6.96310401e-02  8.94153118e-03  1.85952187e-02\n",
      "   1.70800090e-02  4.73970175e-03  5.08469343e-03  8.40461254e-03\n",
      "   4.90385294e-03  9.72682238e-03  4.37488556e-02  4.75537777e-02\n",
      "   4.95262146e-02  7.18187094e-02  2.32190490e-02  3.51954699e-02\n",
      "   1.55891180e-02  3.24624777e-03  3.26138735e-03  5.06007671e-03\n",
      "   7.10368156e-03  1.27706528e-02  1.60021782e-02  3.16596031e-03\n",
      "   3.27467918e-03  4.97192144e-03  5.30314445e-03  9.72080231e-03\n",
      "   1.58379078e-02  3.27813625e-03  3.37016582e-03  5.67275286e-03\n",
      "   5.35500050e-03  1.02080703e-02  1.86730027e-02  1.14825368e-02\n",
      "   1.25288367e-02  1.83985829e-02  6.28304482e-03  1.27732754e-02\n",
      "   1.84519291e-02  4.40210104e-03  4.30357456e-03  7.99804926e-03\n",
      "   4.05305624e-03  9.11653042e-03  1.88891888e-02  3.37278843e-03\n",
      "   3.75139713e-03  5.80608845e-03  3.60143185e-03  7.78150558e-03\n",
      "  -1.27503145e+00 -1.11920822e+00 -2.23186803e+00 -1.57616317e+00\n",
      "  -5.75977802e-01 -1.89794600e+00 -7.02292800e-01 -7.14424610e-01\n",
      "  -6.41600370e-01 -7.33172059e-01 -7.76785254e-01 -4.23837757e+00\n",
      "  -1.59313595e+00 -5.05068827e+00 -9.21227574e-01 -7.82435775e-01\n",
      "  -1.18445480e+00 -2.88968849e+00 -1.93793118e+00 -1.65921545e+00]]\n",
      "(1, 260)\n"
     ]
    }
   ],
   "source": [
    "pred = model.predict(scaled_input)\n",
    "# pred contains parameters for distributions, not actual points on the graph.\n",
    "print(pred)\n",
    "print(pred.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 910,
     "status": "ok",
     "timestamp": 1617249000269,
     "user": {
      "displayName": "Vector Zhao",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Ghd3FrP7k66YYE7kRBdxDz2xnB_I-RtbyItWFzX=s64",
      "userId": "13898113326268806207"
     },
     "user_tz": -480
    },
    "id": "2KTNlCDGiWsL",
    "outputId": "4aa2ac50-2ee3-41ae-9fae-912909270ccb"
   },
   "outputs": [],
   "source": [
    "# split up the mixture parameters\n",
    "mus = np.apply_along_axis((lambda a: a[:N_MIXES*OUTPUT_DIMS]), 1, pred)\n",
    "sigs = np.apply_along_axis((lambda a: a[N_MIXES*OUTPUT_DIMS:2*N_MIXES*OUTPUT_DIMS]), 1, pred)\n",
    "pis = np.apply_along_axis((lambda a: mdn.softmax(a[-N_MIXES:])), 1, pred)\n",
    "\n",
    "for m in range(OUTPUT_DIMS):\n",
    "    locals()['mus'+str(m)] = []\n",
    "    locals()['sigs'+str(m)] = []\n",
    "#     print(m)\n",
    "    for n in range(20):\n",
    "#         print(n)\n",
    "        locals()['mus'+str(m)].append(mus[0][n*OUTPUT_DIMS + m])\n",
    "        locals()['sigs'+str(m)].append(sigs[0][n*OUTPUT_DIMS + m])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 606
    },
    "executionInfo": {
     "elapsed": 1940,
     "status": "ok",
     "timestamp": 1617250775876,
     "user": {
      "displayName": "Vector Zhao",
      "photoUrl": "https://lh3.googleusercontent.com/a-/AOh14Ghd3FrP7k66YYE7kRBdxDz2xnB_I-RtbyItWFzX=s64",
      "userId": "13898113326268806207"
     },
     "user_tz": -480
    },
    "id": "EE1cZZFriWsM",
    "outputId": "6366b696-ae7a-4bf8-91f3-735125b51579"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkwAAAJNCAYAAADDOCpmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeZxkVX338c+vqnrvno2eYRkYBgURRGWTiCLIE2OI0eCOiUZJjMZ9ze5j4vaoiEv0iTwRV9S4ohhEVCKCIKuDDKsg28wwMzD70t3TSy3n+ePc23WrurqruurW/n2/Xv2qW7furXvuTNfpX/3O755rzjlEREREZH6JZjdAREREpNUpYBIREREpQwGTiIiISBkKmERERETKUMAkIiIiUoYCJhEREZEyUrXsPDo66tauXRtTU0SkE9122207nXMr630c9UciUk4t/VFNAdPatWtZt25dLW8hIh3OzDY24jjqj0SknFr6Iw3JiYiIiJShgElERESkDAVMIiIiImUoYBIREREpo7aAKdnDl6++j+vvfSym5oiISDvL7tvKxC8uIDe+s9lNEYlVTVfJDTzxdL5340MY8K01f8iK4f6YmiUiIu1o3zdfS3bb78juepgl532h2c0RiU1NGaaeVccA4ID1j+yKoz0iItLGstt+B0D6sbub3BKReNUUMCWXrJpdfnjb/pobIyIinSHRN9zsJojEqqaAKdGTH4LbuudAzY0REZHO4HLZZjdBJFY1XyVnwePW3RO1vpWIiHSK7EyzWyASq5oDplOe6G/JsnXPAZxzNTdIREQ6gZXfRKSN1BwwHbVqhIHeJNPpLBPTmTjaJCIi7S5R00XYIi2n5oBp+XAfy4f7ANg9Pl1zg0REpD1FRxkskWxiS0TiV3vANNTH8iEfMO1RwCQi0r0y+b8BzuWa2BCR+NUcMK0Y6WPFsAImEZFu5yIBk4q+pdPUHjANRYbkJhQwiYh0rWx6dtFlVdMqnSWGGqZ+DcmJiAgulw+YlGGSTlNTwDRx988Y7k8xMtALwPhUusweIiLSsQoyTPp7IJ2lpoBp6qGbMDNGBnoAGJvUNwoRkW5VECTlNCQnnaXmITkgEjDpG4WISNeKBElOQ3LSYeIJmPp9wKQhORGRLhbNMKnoWzpMTBkmX8OkDJOISPcqGJLLzuh2WdJRNCQnIiLxKM4qqY5JOkgsAdNgXwoDDsxkyGQ1u6uISDcqmFYANCwnHSWWgClhxvCA6phERLpaUYCkwm/pJLEETJAfllPAJCLSneYESJqLSTpIfAFTvwq/RUS6WlHN0pwhOpE2FlvANKzCbxGR7lacUVKGSTpIjBkmDcmJiHQzV5xhUg2TdJDYAqbBvhTgr5QTEZEuNKeGSX8PpHPEHzBN6wMiItKN3Jyr5DTiIJ0jvoCpVwGTiEhXmzMPk4bkpHMowyQiIrGYm2HS3wPpHAqYREQkHnOukivMMGX3bmbP/zuH6bt/3MBGicSjDgGTxqxFRLpRcc1S8fOJqy8g89hd7P/emxrZLJFYxBYwDfX5aQUmdJWciEh3Kr7ZbnEAlZluYGNE4hV/hmlKAZOISFcqMyRnqf4GNkYkXqphEhGRWJQbkiOXbWBrROIV45CcJq4UEelqc6YVKA6Y9PdB2lfsGaYJZZhERLrS3GkFimqWzPKvOdeIJonEJraAqb8niQHT6SzZXC6utxURkXYRZpgS/gv0nCG66JCcZgGXNhNbwGRmkTomjVOLiHSbMECy3iG/IjP/EJ2umJN2E1vABDDU76cW0FxMIiJdKBiSs95BAFzRVXIFGScFTNJmYg2YdD85EZEulgszTD5gWmiaAWWYpN3EGzCp8FtEpGuFRd/WM1+GKRNZVsAk7aUuAZMyTCIiXSgIkGZrmIoCJjQkJ21MAZOIiMTCzalh0pCcdI76BEyavFJEpPuE0wr0DvjHBYu+i7JPIi1ORd8iIhKLfIZpKHhePNN3dFqBqYa1SyQOGpITEZF4FM/DtMC95TQkJ+0m1oBpIMgwTWpITkSk+4TTCoRXyRUHRSr6ljamDJOIiMQiPyQX1jAVZ5iiRd+qYZL2Up8aJmWYRES6T3GGac5VcpF5mJRhkjajDJOIiMRi9l5yfaXnYdKtUaSdxVvD1KcaJhGRrlVU9L3wPEy6Sk7ai6YVEBGRWLhc4a1RMpt+Q3bvluC1LLhcfuOs/k5Ie9HElSIiEo8go5QYOmh21cTPP1TwWmhO9kmkxamGSURE4hEGTMuPmF018/CvgYWH50TaQczzMCUBmJzO4JyL861FRKTFhUNyif4ls+uSy9fgnCO77b7CbZVhkjYTa8DUm0rSk0yQyTnS2Vz5HUREpHOEQVCih+EXfcwvJ3uZ/PVF7P3yi4u2VYZJ2kusARPks0walhMR6R6zRd1mkEjSc8Spfv30OBNXXzB3ewVM0mZiD5hUxyQi0oUi2SUzw3qDySvTByCXnX97kTZRhwyT5mISEek2Yf0SSf83gESPf5wnMFKGSdqNMkwiIlK7cNLKMFAKAqfZQGqe7UXaRf0CJmWYRES6x2yGyQdMVjbDpIBJ2kv8AZNm+xYR6TrhEJuFQ3LhYy4zG0QVyGhITtpL/DVMs/eTK1HkJyIinSm81UmiMMPkshmw5JzNlWGSdlO3IbmJaX0YRES6RRgAWbKwholcGktEAqaefv+ogEnaTN2G5CanlWESEekauSAASvb6x0QYMGUhEjCFN+bVVXLSblT0LSIiNctnmPzfADPLB02RW2VZ75BfUMAkbaZ+8zCp6FtEpHuEQ2zRAu9wWI65AZNqmKTd1LGGSQGTiEi3cJGZvkOzhd+RuZhstoZJGSZpL/WrYdKQnIhI95idVqBEhik9lV9n/s+OMkzSbjTTt4iI1MxlCyeuhMjklVFm/lEBk7QZ3UtORERqVzytAERqmPIsuIpOV8lJu4k9YBpShklEpPvkShR9J+YGTKQ0D5O0p7rN9K1pBUREuseciSsBKxEwWe9AwfYi7aIOQ3J+grID0xlcZO4NERHpYOEQW8GQ3Nwapr4TXjS7vf5GSDuJPWDqTSXpSSbI5hzpbC7utxcRkRYUFn1bYv4M09K//gF9T3mRH6pzzs8CLtImYg+YQFfKiYh0nZITV0aWU/30rn2mnwE8XK/Cb2kjdQmYosNyIiLS+UrVMEWLvi0ZXdaVctJ+6pRh8h8YBUwiIl2iRIZpvuApn2FS4be0j7pmmDQXk4hId3C5hTNMhdmm4JYpyjBJG6lLwDSk+8mJiHSXUveSK8g2RTNMvYX7iLSBOmWYNNu3iEg3cSWLvufJMKV6C/cRaQO6Sk5ERGo3W/Tdm18XDZJK1jBpSE7aR30yTJrtW0Sku8xmmObWKgFFNUzKMEn7qU8NU68yTCIi3aRs0Xep+ZkyyjBJ+6hrhmlyRrO4ioh0hTLTClhBDZO/Aa9LTzWmbSIxqHMNk9KtIiLdIJwioKJpBcIb8KYPNKRtInGo61VyB6aVYRIR6QolphWY74o56xkEYP/330pmxwNM3vI19l7y5+SmxxvRUpGqpMpvsnhDKvoWEekq4fCa9fTPrpt3SK7HZ5jITLHvkj8nt/8xAKbv+m8GTn11A1orsnj1nYdJRd8iIl3BZaaBfH0SUJhtKhUwwWywBJAb31G/BorUSPMwiYhITdJb1ucDn4IMU+kaJnoHS76Pm9xbj+aJxKIuQ3KDmulbRKQrpDf+hr1ffvHs8+iQHPPcGiXRN1LyvXJT++JvoEhM6pph0r3kREQ62/R9Pyt4Hh2Ss3kKwK1vuOR7uUkFTNK66jwPUwbnXD0OISIirSBXeDV0YYYpemuU5OxyavXTS76Vm9ofb9tEYlSXgKknmaAnmSCbc0xncvU4hIiItAA3PVa4ItWXXy6Ykym/3HP4SSz96x8Wvg7kVMMkLawuARPAcL//IExMafJKEZFOlSvKCoVzLEHxJJbJgu161/4ByZXHFKxzqmGSFlbHgMmnYscmFTCJiHQqN1WYYbL+JfknidK3SZl9uW9JwXPVMEkrq1/ANOA/HOPKMImIdKziITkzyy/PN61A+PpAUcA0M4HL6m+GtKa6BUwj/QqYREQ6XUGhdnFQlOydXbRSAVP/kjnrnG6PIi2qLvMwAYwM+A+KhuRERDpXLsgw9aw9nf4TX1HwWkGQlJz75ybRv3TOOjczAYPL422kSAzqFjCFRd9jyjCJiHSssIZp6Wu+jhXP4J0sfWuUUOkM09icdSKtoO5XyY0rwyQi0pFcZhoyU/4KuMj94WZF7x9XKmCKzvgdBFtuSkNy0ppU9C0iIlUJs0vWN1JQ7B2yMhmmxMCy2eXUyif595xRwCStSUXfIiJSldy0L/i2ErVIQOGQXIkapuToE2eXE0sPC95TQ3LSmupY9B3UME3O1OsQIiLSRGGGKdFf+ma6VmZILnXEKQye/V6Sy9eQfuRG/566Sk5alIq+RUSkKuGUAgW1SFHRDFP0likBM2Po7PcAkNl6p39PBUzSolT0LSIiVQmvaCt1tRsU1jBZiYCpYNsg6FLAJK2qfjVMs0XfmXodQkREmmi26HueIblohsl6+hd8L+sb9u+pGiZpUQ0ImNI45+p1GBERaZLwxrvF94QLFdQtlc0whQGTMkzSmuoWMPWmkvSmEqSzOabT2XodRkREmiQ/JDdfhilya5RUuQyThuSktdUtYIJ8lmm/6phERDpOvui7dIYpOpVA+YBpyL+nAiZpUXUNmJYO+hTsvgOaWkBEpNPkDuwCIDHPvd8KgqRUb8ltZrcNMkyah0laVV0DpuVD/gOyd2K6nocREZEmyI3tACAxsqrk64mB/ISWbubAgu+VUA2TtLi6BkzLhnyGae+EMkwiIp0mNx4ETMMr591m6Hn/TM+Rf0DvUacv+F6qYZJWV7eJKwGWBRmmPcowiYh0nNz4dgASIwfPu83gmW9j8My3lX0vTSsgra5BGSYFTCIi7SY3sWt26oBiLjODO7AbLIENrqj5WLMB08yEpqKRllTngCmsYdKQnIhIO8kd2MPuzz6HfV99ZenXw+zS0CiWSNZ8PEv1+WkIsmnI6Eu2tJ46F30rwyQi0o7Sj9yAm9pH5rG7yI3vnPN6dudDACQPekJsx9TkldLKGjIkt0cZJhGRtpKbyAdJ2T0b57ye2fEAAMmVR8d2TNUxSSur8zxMmlZARKQd5Sb3lVwOZWcDpmNiO2YYMOVmJmJ7T5G4NKyGKaciPhGRtuEiQZKb3Dvn9ezezQAkV6yN7ZiJcGqBKWWYpPXUNWDqTSUZ6kuRc45x3R5FRKRtRIOkUgFTbt9WAJJLD43tmBqSk1ZW14AJ8nVMu8c1LCci0i5yU9EhuRIB0/7HAEgsOSy2Y+Ynr1TAJK2n7gHTyiX+XkI7x6bqfSgREYnJQkNyuakxH9T09GMDy2I7pgX3pMsd2EN6y3rSm2+P7b1FalX/gGnpAADb903W+1AiIhKTaFYpN7mn8LUgu5RcchhmFtsxE8EEmJkt69n7xXPZ++WXkt23Jbb3F6lF3QOmVUt8wLRDAZOISNsoqGE6UJRhCuqXEkviq1+CfMA0fc9PIJeB7AwzD10f6zFEqlX/gGmpH5Lbvl8Bk4hIuyioYZoqnFYgG9YvLY2vfgnAhoJbrGTzc/dlHrs71mOIVKsBNUwakhMRaScuMwMzB/LPi2uYZq+QizdgSpS4J11258OxHkOkWg2rYdqxX0XfIiLtwBVllIonrsxfIVefIbmo7C4FTNIaGjAkl69hyuZy9T6ciIjUKAyQwiE3N7UPF5l8OBvWMMWdYRo5eHbZeofAEuT2bcal9YVbmq/uAVN/T5LRJf1kco7H92pYTkSk1bngqrjE8MGQ6oNsGtL5/ju3PxiSiznDZJEMU2JkFYnla8A5sns2xXockWrUPWACOOIgP3vrozt1B2oRkVYXZphsYBmJ/qXBusg0A/Uq+jaj76nnAtB/6l+SOugJgIblpDU0JmAaHQLg0V0KmEREWl1Yw5QYWIoNLC1Yl5vaj5seh54BLAim4jRy7idZ8hdfZeD0N5A86ChAhd/SGlKNOEiYYdq8U3egFhFpdbkDfkjOBpbNBkVh1ml20sql8U5aGbLeQfqe/Hx/jFFlmKR1NCjD5AOmDTt0fyARkVaXG98BQGJolERw65PwVin1mrSylGRkSM7lcgWF5yKN1pCA6ehD/DeUhx7fTzqrK+VERFqZm9gFQGJ4NDIk52uY6jVpZSnJlU8CIL3xFnZ+6Cj2XvxCclP64i3N0ZCAaWSgh8NXDJHO5nhk2/5GHFJERKo0m2EaXjmbYZodktvr7+0W9xVypSSXHELq8JODRmXIbFnP5I0X1/24IqU0JGACOHa1/9Ddt2VvmS1FRKSZchM7AUgMHYSFQ3JBXVN2r7/EP7n8yIa0Zckr/5OBM97C4Nl/B8DUbf+F05x+0gQNC5iOP2I5AOsf2dmoQ4qISBVmA6bhlbNDb9m9m/3jbh8wJZavaUhbkstWM/z89zH43HeRWHY4ubFtZDbf1pBji0Q1LGA67ehVANz28E5mMtlGHVZERBbBZTPkxrYDkBheNZtJyu7ZCEAumEQyuaIxAVPIzOg7/gUATN/zk4YeWwQaGDCtWjrAEw9ewlQ6y28fVpZJRKQVZfduguwMiSWHYr2DJFcEAdPujeQm95Eb3w7JXhIjhzS8bfmA6QoNy0nDNSxgAjjrKT61e/lvNjTysCIiUqHM5tsBSB18HIAPjFL9uImdzNz3c//aYU/FEsmGty11+Cl+WG7/Y6Q33NTw40t3a2jA9CcnH0FfT5LbHt7JXZt2N/LQIiJShsvlmFr3LQB6jj4LAEsk6H3icwAYu+zd/rUjTm1K+yyRoP9pLwVgat03C15zuazmaZK6ashM36ElA7287JlH8a3rH+TC/17Pp1/3LEaX9DeyCSIiUoLLphm7/B9Ib7wZG1hG/9NfNvta/0mvZOb+/5l93nvs85rRRN+WU1/DgRu/wPTdlzM2uByXniSz9S6yO36P9QzSs+YZ9D31XPqOfwHWO9i0dkrnaWiGCeDPzziaYw5dyra9k7zli9fz/ZseYtOOMbI5fTMQEWm03PQ40/ddxZ4vvIDp278HPf0sOe9iEoPLZ7fpPe5PGDz770gsO4KB0/+GnrWnN629yWWrGfrDfwRg6tZLmL79e2S3/Q5yWdz0GDMP/JKxH76TXZ84kbEfvZfp3/+S7Ng21TxJzayWFOapp57q1q1bt+j99k/O8JFLf8sdG3blG4Kf4HKov4ekGYmEkUwYCTPqcLsiEWmQi9545m3OubqP4dz0ziPdk5/8ZP8k0q85In3cfP1dwfr5tl/ofRa3T8HeDTl26e1dLjM7qzdAYtkRLDnvP+lZfSKtbuaBa0lv+g2J4ZWkDj2B1CFPITe1j5n7fs7UHT8g82jR1ANmWO8wJHvAEpglwMz/EHmMbB95Ms/6+V8r2KrCfYr2qrANC+1Ter1VepwOtOItP6+6P6opYDKzMeD+Re62FNhXxeGq2a+afUaBai7ja1T7dE551ZyXzqn6fard71jn3EgVx1qUKvsj0O9Eo/eB1j6navfTOVV/nGr3a2x/5Jyr+gdYV8U+F1d5rEXvV+U+iz6nBrdP51TDeemcam5fw37/GnUc/U7onGI6ls6p8e1raH/U8Bom4McN3K/aY1WjUe3TOdVG51T9PrXs18r0O9HYfarVib/nOqfa9mtof1TrkNw614DahEbSObWPTjwvnVPrH6fROvG8dE7tQedUqNYMUyfeNlrn1D468bx0Tq1/nEbrxPPSObUHnVNETRkmERERkW7QjBomERERkbZSUcBkZueY2f1m9qCZ/VOJ1/vM7LvB67eY2dq4Gxq3Cs7pPWZ2r5ndaWZXm9mRzWjnYpQ7p8h2LzczZ2YtPzZdyTmZ2SuD/6t7zOxbjW7jYlXwu7fGzK4xs9uD378XNKOdi2FmXzGz7WZ29zyvm5l9LjjnO83s5BqOpf5I/VFTdGJ/BJ3XJ9WtP6rgErwk8BDwBKAXuAM4vmibtwD/GSy/CvhutZftNeKnwnM6GxgMlt/cCecUbDcCXAfcDJza7HbH8P90DHA7sDx4vqrZ7Y7hnC4G3hwsHw9saHa7KzivM4GTgbvnef0FwE/xs+I9E7iljv9+6o/a4JyC7dQftcd5tVWfVK/+qJIM02nAg865h51zM8B3gHOLtjkXuCRYvhT4Q7OWnp+77Dk5565xzh0Int4MHN7gNi5WJf9PAB8GPgFMNbJxVarknN4AfN45twfAObe9wW1crErOyQFLguWlwNYGtq8qzrnrgIXuqH0u8HXn3QwsM7NDqziU+iP1R83Sif0RdGCfVK/+qJKAaTXwaOT55mBdyW2ccxn8zJsHVfDezVLJOUW9Hh+NtrKy52RmJwFHOOeuaGTDalDJ/9OTgCeZ2Q1mdrOZndOw1lWnknP6APAaM9sMXAm8vTFNq6vFfuZqeR/1R82n/qg9+iPozj6pqv4oVcEbl/pmVnxpXSXbtJKK22tmrwFOBc6qa4tqt+A5mVkC+AxwfqMaFINK/p9S+DT4c/Hfuq83sxOcc3vr3LZqVXJOfw58zTn3KTM7HfhGcE7tfPfQuPoI9Ufqj5qlE/sj6M4+qao+opIM02bgiMjzw5mbjpvdxsxS+JTdQumwZqvknDCz5wHvA/7MOTfdoLZVq9w5jQAnANea2Qb8uO3lLV5oWenv3n8759LOuUfw9xI7pkHtq0Yl5/R64HsAzrmbgH78PZ3aWUWfuZjeR/1R86k/ao/+CLqzT6quP6qgeCoFPAwcRb4g7ClF27yVwiLL7zW76CuGczoJXwh3TLPbG9c5FW1/La1fZFnJ/9M5wCXB8ig+zXpQs9te4zn9FDg/WD4u+CBbs9tewbmtZf4iyz+lsMjy1jr++6k/aoNzKtpe/VFrn1fb9Un16I8qPfALgN8HH9j3Bes+hP+mAz7a/D7wIHAr8IRm/2PFcE6/ALYB64Ofy5vd5lrPqWjblu+gKvx/MuDTwL3AXcCrmt3mGM7peOCGoONaDzy/2W2u4Jy+DTwGpPHf3l4PvAl4U+T/6fPBOd9Vy++e+iP1R616Tu3YH1V4Xm3VJ9WrP9JM3yIiIiJlaKZvERERkTIUMImIiIiUoYBJREREpAwFTCIiIiJlKGASERERKUMBU5OZ2WfM7F2R5z83sy9Fnn/KzN6zwP7LzOwt9W7nAscfDx4PM7NLy2y7doG7R18Y3N37whja9C4zG4w8v9LMltX6viKdzsycmX0j8jxlZjvMrOrbl5jZvxQ9H6+ljRUec7avMbNTzexzZbZ/7nznaGbfDu5o/+4Y2lX8b3Fjre8pjaOAqfluBJ4Fs7cLGAWeEnn9Wfj5L+azDH939oqZN+//vZklF/N+AM65rc65ly92v4i/BU52zv19UVsquX1PsXcBswGTc+4FrrVvTSDSKiaAE8xsIHj+R8CWGt/zX8pvUpkq+6Z1zrl3VHm8Q4BnOeee5pz7TNFr1fRNBf8WzrlnVdMuaQ4FTM13A0HAhA+U7gbGzGy5mfXhZ1W93cyGzexqM/utmd1lZuHdpD8OPNHM1ofZGTP7ezP7TfCt6IPBurVm9jszuwj4LYXTwmNmG8zsX83s18ArzOwNwXvcYWY/CDM2ZnaUmd0UvPbhyP7Rb3Rrzez6oK2/NbMFOwUzuxwYAm4xs/PM7Gtm9mkzuwa4wMxOM7Mbzez24PHYYL+kmX0y+Pe408zebmbvAA4Drgn2D89tNFh+j5ndHfy8q+jf5otBluuqyB8MkW7zU/xMyODvIfbt8IUFPovnm9kPzexnZvaAmX0iWP9xYCDon/6r+ECl+qoS24yb2YfM7Bbg9KCf+k3wGb7YzCzY7pSgv7oJP9t7uP9s9mi+9i/gKmBV0P7nmNm1ZvZRM/sV8E4ze5GZ3RK83y/M7ODgOMNm9tVI3/SyUv8Wls/Qm/ks+93BPudF2n6tmV1qZveZ2X+F5ytN0OwZOfXjADYAa/BZljcBH8bPvPps4LpgmxSwJFgexc9ibBRN/w48H7g4eC0BXAGcGWyXA565QBv+IfL8oMjyR4C3B8uXA68Nlt8KjAfLs+3AZ3f6g+VjgHXF25Q4/nhk+WtBu5PB8yVAKlh+HvCDYPnNwA8ir62InMto0bmNAqfgZ3UdAoaBe/C3nFgLZIATg+2/B7ym2b8X+tFPo3+AceBpwKX4GdPX428ke0Xw+nyfxfPxt9dYGuy3ETgifM/iYwSPJfuqEm1ywCsjz1dElr8BvChYvhM4K1i+MNIfVdL+2W2Kjl3QZ+FnJL8o8nw5zE4A/TfAp4LlC4B/j25X5t/iZcD/AEngYGATcGjQrn34e50lgJuAM5r9e9KtP9WkFCV+YZbpWfhp9VcHy/vwQ3bgO5WPmtmZ+MBnNf6DVez5wc/twfNhfNCyCdjonLt5gXZ8N7J8gpl9BD/kNwz8PFj/bPyHG3xndUGJ9+kB/sPMTgSywJMWOOZ8vu+cywbLS4FLzOwYfOfZE6x/Hv6eYRkA51y5G6yeAVzmnJsAMLMfAs/BB4GPOOfWB9vdhu8oRbqOc+5OM1uLzy5dWfTyfJ9FgKudc/sAzOxe4Ej8vdTmM19fdV3Rdln8F6PQ2Wb2D/gvZiuAe8zsOmCZc+5XwTbfAP6kxDEXan+lov3k4cB3zexQ/H3YHgnWPw9/H0MAnHN7yrznGcC3gz5vW5DBegawH3+fs80AZrYe3zf9uop2S40UMLWGsI7pqfghuUeB9+I/LF8Jtnk1sBI4xTmXNn+H7/4S72XAx5xzXyhY6TvAiTLtiL7+NeDFzrk7zOx8/DedULn76bwbf9+rp+O/FU2V2b5cWz4MXOOce0lwHtcG662CtkQtlMqO3v09C2hITrrZ5cAn8Z/7gyLr5/sswtzPULm/LyX7qhKmwi9PZtYPXIS/99ejZvYBfD9YaV+wUPsrFe2b/i/waefc5Wb2XOADwfp69k36u90kqmFqDTcALwR2O+eyQaZkGXA6PgUL/pvR9iBYOhv/7Q1gDBiJvNfPgb82s2EAM1ttZquqaNMI8JiZ9eCDtWhbw29Or56zV76tjznncsBf4tPMtVhKvvD0/Mj6q4A3WVB8aWYrgn0bVHgAACAASURBVPXF/yah64AXm9mgmQ0BLwGur7FtIp3oK8CHnHN3Fa2f77O4kHTQjxSrpq8KvyTuDPZ7OYDzF3XsM7MzgtcX6psW2/6FRN/vdZH1VwFvC5+Y2fJgcb5/i+uA88zXZa7El1HcGkP7JEYKmFrDXfgam5uL1u1zzu0Mnv8XcKqZrcN3BvcBOOd2ATcExYIXOueuAr4F3GRmd+FrEUoFD+W8H7gFP65+X2T9O4G3mtlv8J1FKRcBrzOzm/HDceUyW+V8AviYmd1AYfD1JfxQ451mdgfwF8H6i4GfWlD0HXLO/RafObsVf25fcs7djogUcM5tds59tsRL830WF3Ix/jNaUPRdTV8VBEZfxPePPwJ+E3n5r4DPB0Xfk/O8RTXtX8gHgO+b2fXAzsj6jwDLg375DuDsYH3JfwvgMnwN1h3AL/H1pI/H0D6JUVisJiIiIiLzUIZJREREpAwFTCIiIiJlKGASERERKUMBk4iIiEgZCphEREREylDAJCIiIlKGAiYRERGRMhQwiYiIiJShgElERESkDAVMIiIiImUoYBIREREpQwGTiIiISBkKmERERETKUMAkIiIiUoYCJhEREZEyFDCJiIiIlJGqZefR0VG3du3amJoiIp3otttu2+mcW1nv46g/EpFyaumPagqY1q5dy7p162p5CxHpcGa2sRHHUX8kIuXU0h9pSE5ERESkDAVMIiIiImUoYBIREREpQwGTSOj222Hfvma3QkTKye6GmQea3QrpMgqYRACuuQZOPhnOPbfZLRGRcjb8ATz8JJh5uNktkS6igEkE4HOf84+/+lVz2yEiC8tNQ/pBvzx5U3PbIl1FAZMIwIEDzW6BiFQiuy2yvL157ZCuo4BJBCChj4JIW8jtzy9n9zavHdJ19FdCBBQwibSLbCRgigZPInWmvxIioIBJpF3kFDBJc+ivhAgoYBJpF7mxyLKmAZHG0V8JEVDAJNIuXOQCjWjwJFJn+ishAmDW7BaISCVyU5Hl6ea1Q7qOAiYRUIZJpF24SMDkFDBJ4+ivhAgoYBJpF9EgKRo8idSZ/kqIgAImkXahDJM0if5KiIACJpF2oYBJmkR/JUSgsOjbuea1Q0QWFi361pCcNJACJpFiuVyzWyAi8ymoYVKGSRpHAZMIQDabX06nm9cOEVmYhuSkSRQwiUBhwKQMk0jrcpEvNDkNyUnjKGASgcKAKbosIq0lGjCRAacvONIYCphEQAGTSNvIFD7VsJw0iAImESgchlPAJNK6XFGNoa6UkwZRwCQCqmESaRdOGSZpDgVMIqAhOZG2UZRhUuG3NIgCJhFQwCTSLpRhkiZRwCQCCphE2sWcGiYFTNIYCphEQEXfIm2jOMOkITlpDAVMIqCib5F2EWaYrDd4rgyTNIYCJhHQkJxIu5gNmAaC5zPNa4t0FQVMIqCASaRdhEXficHgue79KI2hgEkEVMMk0jbCDJMCJmksBUwioBomkXZRnGEqnpdJpE4UMImAhuRE2oVThkmaQwGTCChgEmkbYYYpLPpWwCSNoYBJBBQwibQLZZikSRQwiUBh3ZJqmERalyvKMKmGSRpEAZMIKMMk0i6UYZImUcAkAgqYRNpGECCphkkaTAGTCChgEmkX4ZCcMkzSYAqYREATV4q0A5cDgs+q9QfrFDBJYyhgEgFNXCnSFoLsEimwnmBZAZM0hgImEdCQnEg7mC347skHTMowSYMoYBIBBUwi7WC2fkkBkzSeAiYRUMAk0g5mg6MUoIBJGksBkwho4kqRVpWbhB3vg7EfMVvDpAyTNEGq2Q0QaZqf/AQmJuCVr1SGSaRVjf8Idn3ULz9xk380FX1L4ylgku6UycALX+iX/+APFDCJtKrs3vyyir6liTQkJ93poYfyy7ffroBJpFW5mfxybiJYSClgkoZTwCTdaePG/PL4uAImkVaVG48sB9km60FF39JoCpikOz3+eH55fLzwNRV9i7SOaMCUjQRMqmGSBlPAJN1p27b88thY4WvKMIm0jlIZJg3JSRMoYJLutDdSSLp/f+FrCphEWoebzC/n9vlHFX1LEyhgku4UHYZTwCTSXLlp2PZemPifua+56ch2QdG3KcMkjaeASbpTNGAqHpJTDZNIY+35LOz5NDz6fP98aj08+iKY/l3RVXLB57ZU0Xd6C2QiQ+0iMdM8TNKdJibyy9PTha8pwyTSWDP35ZddDjY9xwdH28YhsSz/WqlpBUj7oGrDSZDdA8dsh+TyRrVcuogyTNKdohkmBUwizZU7EFnen88kZbYXDclFMkzRIbn0VsjuADIwta4hTZbuo4BJulM0YJqaKnxNAZNIY2Ueyy9HZ/YmVzgk58KAqaiGKbMl8l4alpP6UMAk3WmhITnVMIk0VnZHfjm8Eg585qnkTN9FNUzR/RUwSZ0oYJLupCE5kdaR3Zlfzu0tfK3ckBzposktt9eliSIKmKQ7KWASaQ0uB9ld+efZSIYJK51hKh6SiwZMyjBJnShgku6kq+REWkNuDxAZBneRAnArDpjmKfouyDApYJL60LQC0p2UYRJpDZmdhc+jwY9zQOTz6UpMK6AMkzSIAibpPjMzkI7MDlx8lZyKvkUaJ1sUMEWH59w0BdmnXOQquWjRt1OGSepPAZN0n8nJwuczM4XPlWESaZw5AVPkuZsCXP75vEXfkc90ZpuvizJVnEi89Bsl3efAgcLnGpITaZ4FM0xThVfJuTAbXFTD5KJZ4ixkd9ejpdLlFDBJ9ykOmDRxpUjzzAZMVvScuQFTqLjou3gbDctJHShgku5TLsOkGiaRxgknnew5MnhelHFyRUPoEARLYUVJtvDWKgCZx+NsoQiggEm6URgw9QTfUDOZwteVYRJpnPAqudSawucLsR4/5UAYNIW1TcmV/jF6M1+RmChgku4TBkxLl5Z+XQGTSOOEGaXZDNOu+bedFQRK4bBcbsw/DpzuHw9cE1vzREIKmKT7hAHTyEjp1xUwiTROccCUq6BgOwyUZuuYggzTyHmAwfgVRTOGi9ROAZN0nzBgWrKk9OsKmEQapzhgqkRxwJQNMky9R8PAmb4IfOKq+NooggIm6UbhPEzzDcmp6FukcWoJmCjKMFkfDDzDL6cfjKV5IiEFTNJ9lGESaQ0uDbl9QBJSqyvfrzjDFNYwWX8+8EpvjK2ZIqCASbqRAiaR1hBml5IHQWJ4ETsWFX2HN+i1PkiuKnxvkZgoYJLu8tBD8N73+mVdJSfSXOEUAslRSAwtsKEVPS3KMIUS/ZBc7peze2NpokhIAZN0lzBYgvkzTKphEmmMbCRgssH5t0sUXdE6X8Bk/ZBYFrz3nnjaKBJQwCTdZTxyV/Ply0tvowyTSGMUBEz9zMkkhYoDJoqKvkPWl88w5ZRhkngpYJLusmJFfvmgg0pvo4BJpDHCgCk16mfunm9YLlGUDbaiGqbZ9X2RDJMCJomXAibpLtEb7Z5zTuFrFny7vewyeNWrYGKice0S6UbRDBMUDcsl84sVDcn1gCUgGdQm5vaC0/C6xEcBk3SX8Aq5b34TDj8ckpFOuSfS+X73u/CpTzW2bSLdZjZgCu4BF80wJSND5pUETIn+/LrEMJDL32NOJAYKmKS7hJNWrglu9JlK5V/r7S3cdqcuSxapq9mAKRgqjwZMiQUCprB2ySKfWeuLbB8My+VU+C3xUcAk3SUMmAYG/GM0YOopqodwrjFtEulW4f3ewgDHohmmSL1hYpCCIbrimb6hKGBaWvj+IjFQwCTdJQyYBoNaiWiQVJxhUsAkUl+5IKAJ644SkRqm6JCc9RYGRKWKvm0gvxxmpMIZwEVioIBJustiMky5HGzcqMBJpF7CgCnMCEUzTGHWCUoETAvUMAEkg6vqFDBJjBQwSXdZKGAqzjB9/euwdi28+90NaZpI18kWBUyJeYbkbLAwYKJEwGSRgEkZJqkDBUzSXYoDpoWG5MJpBT772fq3S6QbVToklxgonWFCAZM0jgIm6S6LGZITkfpxOcjt98vhxJQFRd8H55dtABLlaphKBUz742uvdD0FTNI90mnIZCCRyAdHCw3JRc3M1LdtIt0mNw44HySFAVB0SC4VDZgGi6YQ0JCcNJ4CJuke0exSOKt3pRmm7dvr1y6RblQ8HAeFQVEycuui4iG5UjVM0aLvhIq+JX4KmKR7FA/HQeUZJgVMIvEqvkIOgEx+MTpNgM1Tw6QMkzSQAibpHuUCpoUyTNu21adNIt2q+Ao5KJzCo2DW76Kr5MIhPBV9SwMpYJLO9bnPwSGHwIMP+ufFk1ZCYZDUH+lwiz3+ePztE+lmpYbklv0t9D0dVv8QUofk16cOKQqYIveNK14HKvqWulDAJJ3rne/0maEvfME/L5dhiq4v9thj8bdPpJuVGpLrXQtHrYeRl+RvyAuQXFVU9F0iYEoowyT1lSq/iUgbyuXyy1NT/lEBk0jrKDUkF2UJGPpjSG+A3qOByHBdeNGGapikgRQwSWfavTu/vHevf6wlYNq6Nb62iUjpIblih/8UyIEl/bxNc8wXMOkqOYmfhuSkM+3YkV/evx8++EH4oz/yz5VhEmkO5+DAryA7Btngi0xi+fzbm/lgCYASAVOpYTpQhknqQgGTdKYwqwQwNgYf+ED+uTJMIs0x/iPY9FzY/ELIBZ/RhTJMBUoETAVX0pUImLIq+pb4KGCSzjQ+nl8unkMpGhhFr5IrFTAtDTrzzZv9LOEiUr3xK/zj5HWQ3uyXE8sq27fUkFw0YCrIMA0H+0zMM5QnsngKmKQzhTfOBR/sRM2XYerrY46lS+GwwyCbhS1b4m2jSLdJP5xfnrrJPyYrDZim564LAyMoCp4S+fvS5cYRiYMCJulM0QzTvn2Fr42M5JejAVMyyRzZLKxd65c3bIirdSLdaeaR/HJ2l3+sNMOUeXTuuoKAqagWKqnCb4mXAibpTNEMU7ElS/LL0YApUeLjcOyxcOSRflkBk0j1nINsiRnzK80wrfhH/3jQP+fXRQOmZFHAFNYx7fwAPLgGpm6ruKkipShgks40vkAafmmkyHS+DNPb3gavex186UuwerVfpyvlRKqXGwc3NXf9fPMwFVv2RnjC72H0I/l1ychs4MWBVzi1wL4v+ezUro8trr0iRTQPk3SmajJMPT3wf/4P/Oxn8P73w6pVfv3BB/tH3U9OpHrZHaXXV5phMoPeYwrX9ayGpX8FMw9Az1FF77uy8PnUbys7jsg8lGGSzrRQhikaMEWvkkul4F/+Ba67Lh8sgQImkThkg6tVE5HPnw3mi7OrdehX4MjrC2f9hsJ70QGkN4Kbqe1Y0tUUMElnmi/D9IxnwDnn5J9HM0ypeRKuCphEapcJAqb+0/LrUqvytzmJW+qwohU5HzSJVEkBk3SmUhmmf/xHuPVWGIp8oy0ekitFAZNI7cIMU88RkZV1CpYA+p6aXw7rpGYeLr2tSAUUMElnKpVhig7FhZRhEmmMMMOUXAU9a/3y0Dnzbl6zoefDwOkw/FIYeblfl36ofseTjqeASTpTmGGKXhFXbcA0OuqHDXbt0mzfItUKM0ypVXD4lbD8HXDQP9bveMnlcOSNcPgP8gXhGpKTGihgks4UZpgOi9QxLC1x+XIlQ3KplA+anCu8qa+IVC4byTD1HQcHfxZ6jmzMscOMVnpDY44nHUkBk3SmMMMUDZhKZZiKr5Kbj4blRGqTCeYxK756rRHCwEwBk9RAAZN0plIZpmqH5EABk0itwpvtpg5v/LFnM0wakpPqaeJK6UylMkzVDsmBAiaRaqS3wLY3w9ALIBMGTKsb347UoUCPvzVLbhISA2V3ESmmDJN0plIZpsOK52VBGSaRetr1cRj/sQ+a3JSftDI5Un6/uFkyP52BhuWkSgqYpPM4l88wnXBCfn109u6QAiaR+jnwy8LnzcguhXqP94+Tv4axy/297UQWQUNy0nmmpiCXg74+OOssePWrYc0aSJT4fqCASaQ+XMbf4y2q99jmtAVg4JkwcQU8/kb/fPBsWPPLhfcRiVDAJJ0nHI4bGoJkEr75zfm3jdYtLVTDdOih/vGaa3z2ani49naKdLL0o0AaEiOQG/PrBp/bvPYMnlH4/MA1MPMg9B7dnPZI29GQnHSecDiukqCm0gzTs54FK1bAli3w939fW/tEukH6Qf/YdxKsuRZWfgyWval57Rk4E5a9xT/2HOPXHfj14t9n+n7I7Iy3bdIWFDBJ5wkzTJUETJVmmEZG8pmqK6+svm0i3WImuA1J79EweBYc9E+Q6Gtee8zgkM/Dkb+CZcGw3OSNi3uPyZvgkeNgw9MgNx1/G6WlKWCSzhNmmKI32Z1PNKgqF2A9//nQ3w+bNsH+/dW3T6QbhBmmVhzy6j/RP87ct7j9xn4IOD8J52QV2SlpawqYpPNEa5jKGRkpvVxKMumLxwE2b66ubSLdYiYImHqe2Nx2lNIbDMkVF6WXk34kvzx9b3ztkbaggEk6z1hQYFouAILFZZggP5fT1q2Lb5dINwmDi/DGt60kdQRYH2QfX9z0AulN+eXFZqek7Slgks4TDpeVuhVKsWQyv7xQ0XdodTCPzJYtpV/fvBkuvdTPBSXSrZxr7YDJEtDzBL8cZsIqEb21igKmrqOASTpPtRmmSiyUYXIOTjoJXvEK+NnPFve+Ip0kt9tPJZAYhuRBzW5NaYsdlstNQnZ7/nl0eE66ggIm6TyLCZie/GR4//vhkksqe++FMkzbt8PO4HLj666r7P1EOtFMJLtk1ty2zCcMmNIVZpjC4bhkMIltejO4XPztkpaliSul8ywmYAL40Icqf+9wxu/t2+e+9vvf55c3bZr7ukgnm3kADvwKRl7W2sNxobBtld5bLhMMx/UdB1Npn0XL7oDUwXVpnrQeBUzSecIapkoDpsVYKGCK3jblscfiP7ZIq8pNwqaz/OX2+74Jw+f49S0dMK31j5UGTGH9UupI6NkD07v9bOYKmLqGhuSk84QZpkqKvhcrvIFvqYBpx478sgIm6Sb7v+2DJYDJX8Huz/rlvqc3r03lhAHTTIW1SGHA1HOkv8oOIPNo7M2S1qWASTrDvn35QGmxQ3KLsVDAFF2ngEm6ycRV/jG50j9mH/ePA6c1pz2VCAOmzMbKapGiAVNPEDClFTB1EwVM0v727YMTTvAF2Rs21DdgWrECEgnYtQsymcLXogHTvn0wORn/8UVajXNw4Fq/vPoHQDBVR2I59D65Wa0qLzHkAzw3A5nHy29fKmBShqmrKGCS9nf11X7+o7Ex+M536lvDlEzC6Khf3ll0A87irNPjFXTCIu1u5veQ3eavHhs4A1Z9EpKrYNWnwZLl92+m2cLvCoblSg3JKcPUVRQwSfu744788m9/C7t3++UVK+pzvHBYLlrkDXMDpuLXRTrR5K/84+BZfgqBFe+CY7bBsvOb2qyKVFr47TKQCaYSSR0BqcP9sjJMXUUBk7S/aMB0zz35wGXlyvocb746pvD5kUeWfl2kEx0IA6Yzm9uOalSaYcpsBbKQOhQSfZEaJt1TspsoYJL2d+ed+eV774W9e32dUb0zTMUBUXiV3FOf6h+VYZJO51wkYDqruW2pRqUZpuiUAhDJMG0Bl61Hy6QFKWCS9rZ/PzzyCPT1wRFH5NePjvqgqR5KBUyZjC8EN4PjjvPrFDBJp0tv8EFDYgX0Ht/s1izebMBUJsMUrV8Cn2VKrgKylRWMS0dQwCTtLcwuPeUp/idUr+E4KB0whcHRypVw6KFzXxfpBBP/Azs/DJPrguc/94+Dz/E3tG03lc72nX7IP/Y+MbKvrpTrNm34Gy4SEdYvPf3pcOyx+fWNDpjCK+IOPTQ/G7gyTNJJdl0Ajz4fdv4rbDwNdn4U9n/dvzbykua2rVphxii9yRd2zye8QW/P0fl1ulKu6yhgkvZ2223+8cQT/Y10Q4cdVr9jlgqYwokqDz10/qvoRNrV9L2w432AwfCf+XU73weTN/n5loZf3NTmVS3RD6k1QAZmHpp/u5ngBr29kYBJGaauo4BJ2s83vwnHHw9vfStceaVf95znwKmn5rc54YT6Hb9UQBRmmA45ZOH7zYm0G+dg+7uBLCz7Wzj8v2H1ZZBYAqTg4P+A5NJmt7J6fU/zj9N3lH7dOZi53y/3KsPUzXTzXWkvDz4Ir3sd5HLwu9/5dcuXw9OCTu+kk+D+++GlL61fG0oFTNEM0+rVfnnTJt/ZmtWvLSL1NvETf+uTxFIY/ZBfN3IuDG0Bl4bk8ua2r1b9T4OJK4KA6ZVzX08/BLm9fmLO5CH59T1rgtcfbkgzpfmUYZL28olP+GDpxBN9Ngfgne/0M3Ank3DddfDQQ4X1THE77DAfBG3dCum0X7dpk39cvdpPZzA6ChMTsGVL/dohUm+5Kdj2Hr88+gFIRWoDE8PtHyxB/gbBU3eWfn3yVv84cFrhl5/Z/W6vX9ukpSjDJO1j61a45BLfaX3nOz7T8+ij+XmPAIaH/U89DQzAmjWwcSM8/LAPzu65x78WTilw7LH+1in33w+HH17f9ojUy66PQfoB6D0Olr+12a2pj/4T/ePUrf4mvMVX+01c4R8Hzihc3/skHzRmHoXMNkgdXP+2SlMpwyTt4zOfgZkZeNnLfEASDsU1Y8grLDC//34/7BYGTOHUBuHr69c3vm0icZi4xgdMAIdcDNbT3PbUS88xkFoN2e0wfVfha+ktMPYjvzzy8sLXLAH9z/TLE1fVv53SdAqYpD1s3QoXXeSX/+mfmtsWyGeSbr/dDwHu3w8HHZSvb3re8/zjpZf6gEqkHbgcTN8NOz8Im/8ESMOKv4PBM8ru2rbMYOgcv7z/W/n1U3f6fwM3CcMvgd4nzN13JKiV3PcVfc67gAImaS3O+Zm7L7sMfvELn8HZsQPe8hY4cMAXc59ySrNbCWef7R8vvxw+9Sm//Md/nM92/emfwsgI3HwzfPzj6kylNeWm4MCvYdfH4dEXwgOj8MhTYecHwE3DsjfDygua3cr6W/Y3/nHvRbD3q7D1fNhwis849RwDh3y+9H5LzoPEMjhwLez8N18ELx3LXA0d+amnnurWrVsXY3Ok62Sz/mq366/3P9ddN3+h9JIlft6lo48u/XojTUz4ovPxcf88kfDZpvBqPYCvf91f0Qdwzjlw4YX1ne6gRZnZbc65U8tvWRv1R2W4nL8FyNRtMLUOJm/wj26mcLvUan9fuKV/A0NnN6etjeYcbD0Pxr4fWZmAZW+ClR+D5JL59933LXjs1X45tQaWvdFfRdh7fHvOft7haumPair6fuPWrfBv/1a4slQApnWt155mrpuchD17fObo3nv986jly+G002BqyhdWb9vmr4r7939vjWAJYGgIvvxl+Ku/8vex+/znC4MlgNe+1m/3+tfDz37mf447zhepr1jhM1CJhM9KzfcjFXvjS7fCjn8Dwt+7Mo+u3HbV7FNu+1r2jT7mgpmps/4xujxn3bS/31nmMaB4NmuDvqf6guaBZ/uht9Sa7vvdM4NDL4G+p/ir4vpO8HNOlRqGK7b0LyA1Ctve4edr2vm//U9ixAdNqUMgOQo2AJbyP6TAkkA1/86tvE9nqynDhJnGGaR2a9bAs5/tJ598znP8pJT1unFu3GZm8lMazGf7dvjgB+Eb34Cxsca1rUUYNCTDxH3qj8pKHgL9J/ufgWfBwOmQXNbsVnUGl/P31hv7Pkz8LAhQpdXYcdX3RzUFTGY2Bty/yN2WAvuqOFw1+1Wzzyiwc5H7VHssnVP1+0B156Vzqn6favc71jk3UsWxFqXK/gj0O9HofaC1z6na/XRO1R+n2v0a2x8556r+AdZVsc/FVR5r0ftVuc+iz6nB7dM51XBeOqea29ew379GHUe/EzqnmI6lc2p8+xraHzVj3OPHDdyv2mNVo1Ht0znVRudU/T617NfK9DvR2H2q1Ym/5zqn2vZraH9U65DcOteI2oQG0jm1j048L51T6x+n0TrxvHRO7UHnVKjWDNPFNe7finRO7aMTz0vn1PrHabROPC+dU3vQOUXUdpWciIiISBdok2u3RURERJqnooDJzM4xs/vN7EEzm3MjLzPrM7PvBq/fYmZr425o3Co4p/eY2b1mdqeZXW1mRzajnYtR7pwi273czJyZtfzYdCXnZGavDP6v7jGzb5XappVU8Lu3xsyuMbPbg9+/FzSjnYthZl8xs+1mdvc8r5uZfS445zvN7OQajqX+SP1RU3RifwSd1yfVrT+q4BK8JPAQ8ASgF7gDOL5om7cA/xksvwr4brWX7TXip8JzOhsYDJbf3AnnFGw3AlwH3Ayc2ux2x/D/dAxwO7A8eL6q2e2O4ZwuBt4cLB8PbGh2uys4rzOBk4G753n9BcBP8dMHPxO4pY7/fuqP2uCcgu3UH7XHebVVn1Sv/qiSDNNpwIPOuYedczPAd4Bzi7Y5F7gkWL4U+EOzlp5bv+w5Oeeucc4dCJ7eDBze4DYuViX/TwAfBj4BTDWycVWq5JzeAHzeObcHwDm3vcFtXKxKzskB4c2rlgJbG9i+qjjnrgN2L7DJucDXnXczsMzMDq3iUOqP1B81Syf2R9CBfVK9+qNKAqbVwKOR55uDdSW3cc5l8DNvHlTBezdLJecU9Xp8NNrKyp6TmZ0EHOGcu6KRDatBJf9PTwKeZGY3mNnNZnZOw1pXnUrO6QPAa8xsM3Al8PbGNK2uFvuZq+V91B81n/qj9uiPoDv7pKr6o0puvlvqm1nxpXWVbNNKKm6vmb0GOBU4q64tqt2C52RmCeAzwPmNalAMKvl/SuHT4M/Ff+u+3sxOcM7trXPbqlXJOf058DXn3KfM7HTgG8E55erfvLqJq49Qf6T+qFk6sT+C7uyTquojKskwbQaOiDw/nLnpuNltzCyFT9ktlA5rtkrOCTN7HvA+4M+cc9MNalu1yp3TCHACcK2ZbcCP217e4oWWlf7u/bdzLu2cewR/L7FjGtS+2jfclAAAIABJREFUalRyTq8HvgfgnLsJ6Mff06mdVfSZi+l91B81n/qj9uiPoDv7pOr6owqKp1LAw8BR5AvCnlK0zVspLLL8XrOLvmI4p5PwhXDHNLu9cZ1T0fbX0vpFlpX8P50DXBIsj+LTrAc1u+01ntNPgfOD5eOCD7I1u+0VnNta5i+y/FMKiyxvreO/n/qjNjinou3VH7X2ebVdn1SP/qjSA78A+H3wgX1fsO5D+G864KPN7wMPArcCT2j2P1YM5/QLYBuwPvi5vNltrvWcirZt+Q6qwv8nAz4N3AvcBbyq2W2O4ZyOB24IOq71wPOb3eYKzunbwGNAGv/t7fXAm4A3Rf6fPh+c8121/O6pP1J/1Krn1I79UYXn1VZ9Ur36I830LSIiIlKGZvoWERERKUMBk4iIiEgZCphEREREylDAJCIiIlKGAiYRERGRMhQwtSAzO8TMvmNmDwV3vb7SzJ7U7HaFzGw8eDzMzC4ts+3aBe4YfWFwR+8LY2jTu8xsMPL8SjNbVuv7ikjr90m1MLM+M/uFma03s/NqfK9lZvaWyPOyfaS0D00r0GKCm4TeiJ/87D+DdScCI8656yvYP+mcy1Zx3Ir3M7Nx59xwhduuBa5wzp1Q4rX9wEpXNGuxmaWcvwdYxYLZgk91zu1czH4isrBm9UmNYmbPBC5wzs253cxi275QfyftTxmm1nM2kA47JgDn3Hrn3PXmXWhmd5vZXeG3ITN7rpldY2bfwk/ChZm9xsxuDb41fcHMksUHMrMNZvavZvZr4BVm9gYz+42Z3WFmPwgzNmZ2lJndFLz24cj+s9mjYPl6M/tt8POshU7SzC4HhoBbzOw8M/uamX3azK4BLjCz08zsRjO7PXg8NtgvaWafDM7/TjN7u5m9AzgMuCbYPzy30WD5PcG/2d1m9q5Ie39nZl8MslxXmdlAdf9lIh2t0X3SR4P+Zp2ZnWxmPw8yW28Kthk2s6uDfuYuMzs3WD9kZj8J+q+7I235eJAVu9PMPll0vFXAN4ETg3Y9cRH94sFmdlmw/o6gz/s48MTgvS4s6iP7zeyrQZtvN7Ozg/Xnm9kPzexnZvaAmX0i1v89iU+zZ+TUz5wZSt8BfGae114G/A+QBA4GNgGH4m/0OAEcFWx3HPBjoCd4fhHw2hLvtwH4h8jzgyLLHwHeHixfHu6Pv+3EeLC8lmDqeWAQ6A+WjwHWFW9T4vjjkeWvAVcAyeD5EiAVLD8P+EGw/GbgB5HXVkTOZbTo3EaBU/Ad9hAwDNyDv83EWiADnBhs/z3gNc3+/9ePflrtpwl90puD5c8Ad+LvO7cS2B6sTwFLguVR/IzuFrTli5H3WgqswN/PLRxNWVbimM/FZ4WibaikX/wu8K5gORkcr6C/K+oj3wt8NVh+cvBv1Y+/AfHDwf79wEbgiGb/v+tn7k8KaSdnAN92PkW8zcx+BTwD2I+/F84jwXZ/iA8UfmNmAAPA9nne87uR5RPM7CPAMnxw8fNg/bPxnRHAN4ALSrxPD/Af5lP1WaCa+obvu3z6eylwiZkdg7+LdE+w/nn4+4RlAJxz5W6qegZwmXNuAsDMfgg8Bx8EPuKcWx9sdxu+cxORytWjT7o8eLwLGHbOjQFjZjZlvi5xAviomZ0J5IDV+GDtLuCTZnYBPgC63vzNl6eAL5nZT/BfyipRSb/4v4DXAgTnv8/Mli/wnmcA/zfY/j4z20i+n7zaObcPwMzuBY7E34dOWogCptZzD/DyeV6zBfabKNruEufcP1dwvOh+XwNe7Jy7w8zOx3/zCpUrdns3/l5XT8cP9U5VcOyF2vJh4Brn3EvM1wVcG6y3CtoStdC/WbR2KovvxEWkUKP7pPBzmaPwM5rD/816NT7jdIpzLm2+frHfOfd7MzsFf1+0j5nZVc65D5nZafiA7VXA2/CBTjmV9ouLsZi+SH+bW5BqmFrPL4E+M3tDuMLMnmFmZwHXAecFdTwrgTPxNxctdjXw8mB8HjNbYWZHVnDsEeAxM+vBd0qhG/CdDUXro5YCjznncsBf4lPUtVgKbAmWz4+svwp4U/DNETNbEawfC9pf7DrgxWY2aGZDwEuAsoWqIjKrmX1SKUvxw3PpoA7oyOA9DwMOOOe+CXwSONnMhoGlzrkrgXcBJ1ZxvPn6xavxJQJhbeUS5u+HwP9bvTrY/knAGvxwobQJBUwtxjnn8H/U/ygodLwH+ACwFbgMP6Z/B74T+wfn3OMl3uNe4H8DV5nZnfgag0MrOPz7gVuC7e+LrH8n8FYz+w2+syrlIuB1ZnYzPs08Mc92lfoE/lviDRQGX1/Cj/3faWZ3AH8RrL8Y+KkFRd8h59xv8d8Qb8Wf25ecc7fX2DaRrtHkPqmU/wJONbN1+AAk7KueCtxqZuuB9+HrjUaAK4Jj/gqfCV+shfrFs83sLvyQ/lOcc7uAG4Ki8+LpUi4CksH23wXOd0VXCEtr07QCIiIiImUowyQiIiJShgImERERkTIUMImIiIiUoYBJREREpAwFTCIiIiJlKGASERERKUMBk4iIiEgZCphEREREylDAJCIiIlKGAiYRERGRMhQwiYiIiJShgElERESkDAVMIiIiImUoYBIREREpQwGTiIiISBmpWnYeHR11a9eujakpItKJbrvttp3OuZX1Po76IxEpp5b+qKaAae3ataxbt66WtxCRDmdmGxtxHPVHIlJOLf2RhuREREREylDAJCIiIlKGAiYRERGRMhQwiYiIiJTRUgHT/swBXnHvhXxk46XNboqISFv5ysNpRi+b4OePZ5rdFJGO1FIB0493r+P3k4/x/Z03knO5ZjdHRKRtvHHdDHvT8K93pZvdFJGO1FIB097MxOzytvS+JrZERKR97J52s8sPjOdwzi2wtYhUo6UCph0z+2eXo8GTiIjM79HJfEZ+fxr2KckkEruWCpj2ZMZnlxUwiYhUZvOBwozSlkllmETi1lIB03h2anZ5X+ZAE1siItI+tk4WB0yqARWJW0sFTGORgGl/drKJLRERaR/FGaXiAEpEatdiAVM+SJrOaRBeRKQSj0/5AGlpj3+uITmR+LVUwDQRyTBN5maa2BIRkfaxd8YHSMeO+C5917QCJpG4tVTAFM0qTSlgEhGpyN60D5CeMGwA7Fb3KRK7lgqYZlx+hloFTCIilQmnEVg75Lv03TPKMInErWUCpozLkiP/IZ9SDZOISEX2BRmmtYNhhkkBk0jcWiZgmskV3v9INUwiIpUJM0xHDfsufa8CJpHYtU7A5AoDJmWYREQqE2aYjhpShkmkXlonYMoVB0zKMImIlDOddUxloScBhw3ki751PzmReLVOwKQMk4jIooXDcUt7oD9pDCQh42Ais/B+IrI4LRMwpZ1qmEREFiucUmBpjxU87s8owyQSp5YJmMIhuWTQJA3JiYiUF9YvLQsCpZGUXz+mJL1IrFomYJoOMkxLUgOAhuRERCoRDsktCW6LMqIMk0hdtEzAlA4yTCPJQQAms8owiYiUs29GGSaRRmiZgCks+h5J9gNza5pERGSufUU1TGGGaVwZJpFYtU7AFGSYhoOAqfiqORERmWs86CqHg8xSmGHan1bAJBKnlgmYwozSYLIveJ7VPCIiImUcyPp+cjDlM0tLggzTmL5zisSqZQKmMKPUZz2kgmalXbaZTRIRaXkHgsBoMOkfh4PAaUwZJpFYtUzANB0MyfUmUvQkfE5ZdUwiIgubKMowzQ7JqfsUiVXLBExhcNRjSXrNf+JVxyQisrDJoJscCjJMs0XfyjCJxKplAqaw6Lsv0UOPJQvWiYhIafkMk3++JJxWQN2nSKxaJ2AKskm9lqJ3dkhONUwiIgsJa5iGUoXTCugqOZF4tU7AlAuH5FL0mGqYREQqMRF8rxyYLfr2j2Oah0kkVq0TMLl80XeYYdKQnIjIwg4EgdHQ/2/v3oMkrco7jv+e7rnszmUvLCwXXVjuCCgqF0VRlkgoCy1NLBQTSGLKRNFELS1jWV5SBKsSc/ESxUtRKcUkyEURbwkqERC8LMt9F1EQ2AU3EFZ2WXZ3ZufS0yd/nPP2bfqdfrv7nX7fmfl+qKnp6X7f7nN6Zw5PP+d5z2FZAWBe5Sdgiq6Ss2KlhokpOQCY2/4wTEbLCoyyrAAwL/ITMFUyTP1cJQcACY2FDNNQMaph8vfvY/gEUpWbgGm6pui7sg4TU3IAMKfxkGEaDrVL0cKV7CUHpCs3AdNU2f/VDxT6mJIDgIQaM0wjLCsAzIv8BEyuepUcU3IAkExjhmmwIBVNmi5LUzNkmYC05CdgKk9L8kXfXCUHAK2VnasUfS8LRd9mVtkehTomID35CZhqlhXoY0oOAFqqvUKuYFa5nzomIH35CZgqywpUr5Jj4UoAiBet8h1tixIZIcMEpC43AVOUTWLhSgBIprKPXNHq7ifDBKQvNwFT7V5ylc13yTABQKzqPnL19w9zpRyQutwETJNR0XehWLOXHDVMABCnVYZpjAwTkJrcBEzTtcsKsHAlALQUV8PEat9A+nITMEULVw6yNQoAJDJODRPQM/kJmJrUMDElBwDxKhmmYv39w1wlB6QuNwFTNP3WX7OXHFfJAUC8KMM03EeGCZhvuQmYqgtXFlmHCQASGIvJMFX2k5smYALSkouAyTnXMCXHVXIA0Eq0j9ysou/KVXI9bhCwiOUiYIoCoz4VVLCCBgphHSam5AAgVrRswOyib/+dKTkgPbkImKrTcf5aWK6SA4DWogxT48KV1RqmHjcIWMTyETCVq9NxktTHlBwAtLQ/JsM0TNE3kLpcBEyVRSvDVBwLVwJAa2MxNUwjbI0CpC4XAVNlW5SQWRqorMPEXzsAxBkvNV9WYLSfrVGAtOUiYJoKU29RZqmfGiYAaKmSYYpZVoAaJiA9uQiYpitLCvii78pK32VqmAAgTpRhGmLhSmDe5SJgqhR9hxqmaKVvir4BIN44GSagZ/IRMLn6q+T6qWECgJaiveQalxUYLEhFk6bK0tQMWSYgDfkImGr2kav9ToYJAOJFe8k1LitgZmSZgJTlI2CqLFzZV/edZQUAIF7cXnJSdXsU6piAdOQjYCo3n5LjKjkAiBdlmBqXFZCoYwLSlo+AqTHDxJQcALQ0V4aJK+WAdOUiYJqOK/pmSg4AmpouO5Wc1GfSQHF2hmmYDBOQqlwETNVlBaK95HzAVFJZZVfOrF0AkFfRFXKN26JEogwTq30D6chFwDTZkGEys0qWqUTABACzjMVcIRcZ9esAs58ckJJcBEzTDcsK1N6m8BsAZotbgylCDROQrlwETFFQNFioDZhYvBIA4lQzTM0fp4YJSFc+AqaGZQWkaj3TFIXfADBLtYap+ZQcGSYgXbkImKIsUn/TDBNLCwBAo7h95CKswwSkKxcBU+NeclK1hqlEwAQAs4yX4hetlLhKDkhbLgKmyYZlBaSa1b6ZkgOAWVplmEbDcLp3moAJSEMuAqbKlFxthqnAat8AECfKHMUtK1CtYepZk4BFLRcBU5RFqr1Krro9Cn/tANAoyjDFLSswTNE3kKp8BExNa5jYgBcA4kQZpuWxGSb/nQwTkI58BExNFq7sq+wnx5QcADRqlWGi6BtIVy4Cpsrmu4XZ6zAxJQcAs423qGFiaxQgXbkImCabLFzJOkwAEK9ylRxbowA9kYuAqWmGiaJvAIjVei85/50aJiAduQiYmhZ9V7ZGIcMEAI2qe8k1n5IbLEhFk6bK0tQMWSagW/kImOZYuJIMEwDMVt1LrvnjZlbJMo3xuRPoWj4Cpjm2RqGGCQBmG2+RYZKkUeqYgNTkI2BqsqxAJcPE1igAMEurGiapWse0d3r+2wMsdvkImNzslb77C1wlBwBxkmSYuFIOSE8uAqbpJssKRLdZ6RsAZhuLaphiNt+VqtknrpQDupd5wFR2ZZVUllRd3VuqrWHiLx0AGkUZpmjPuGZY7RtIT+YBU23Bt1n1Dz9aVoCtUQBgtiQZJtZiAtKTfcAUAqLa6TiJzXcBII5zruVK35I02k8NE5CWzAOmSecv3xgs9NfdHwVMJYq+AaDORBgWlxWlgsVPyUU1TOwnB3Qv+4CpHAVM9R+TBliHCQCaiqbYhueYjpO4Sg5IU44CpoYMU1hWYIp1mACgzt4QAEVTbnGoYQLSk3nAVC36bpyS4yo5AGimEjDNUb8kcZUckKbMA6aJkGFaNquGKVqHiSk5AKgVrdw9MseSAlI1oGJKDuhe5gHTVAiYBgrNr5JjaxQAqJd8Si6qYZr3JgGLXuYBUyXD1DAlFwVQXCUHAPWiAKjVlFy0qOXeaTJMQLcyD5gqNUwxGSbWYQKAelEANNpiSo6ibyA9mQdMEzFXyUU/TzIlBwB19oQaptH+uY+j6BtIT+YBU1TDNNg4JReKvqNlBwAAXlTD1KromwwTkJ7MA6aJmJW+o5+nHAETANTam7CGia1RgPRkHjBFC1M2rvQdLTMwQYYJAOrsm25v4Uq2RgG6l3nAFLfSd7WGiYAJAGolXbhysCAVTZoqS1MzZJmAbmQfMLnmNUzRMgMETABQrzIl1yLDZGaVLNMYK7QAXck+YIrJMPVZUSZTSWXWYgKAGtGyAq2KvmuPoY4J6E7mAVN1WYH63LKZVQu/WVoAACqSLlxZewxXygHdyTxg2l+ekiQtLwzOeiyalqPwGwCqkm6NIpFhAtKSo4BpYNZjlcJvlhYAgIrdUz74Wd1i4UpJGo6ulGMYBbqSfcA04wOmoeLsDBNXygFAPeecdvlhU6sHyDABvZJ9wDRnhonVvgGg1r6SVHLSUFEaLLYOmKKg6tkpAiagG5kHTOPlSUnSULOAiaUFAKBOFPgckCC75I+LzpuvFgFLQ+YBUzQl1yzDVFntmxomAJAk7Yrql2YPmU2tCoHVLjJMQFcyD5gqGaY5aphYVgAAvHbql6RqJoopOaA7mQdMc9UwDYQaJpYVAACv/Sm5EDBNEzAB3cg0YCq7ciUYWtZsSo4aJgCo027AtCoMrbuoYQK6kmnANFGelpPTMutX0WY3hXWYAKBeVIu0qs0MEzVMQHcyDZgq03FN6pck1mECgEZRpuiAhEXfq8Nq4LsJmICuZBowRQXfzeqXpOpVcgRMAOB1WsNEhgnoTqYB076ZCUnSaHFZ08cH2EsOAOo802bAtLpmHaayI2gCOpWLgGkkJmBaRg0TANT5vwkf9ByyLFnA1FcwjfZJTtIehlKgYzkJmJY3fZx1mACg3tNtBkySdOCgP/Z3k2SYgE5lHDDtl9Q6wxQVhwPAUuacazvDJEkHhYBpBwET0LFMA6a9LabkhsPVc2Mzkz1rEwDk1XPT0sSMNNInjfQnD5gODsHVjgkCJqBTOZmSax4wDRX8/WNlAiYAiLJLB7eRXao9/mkCJqBjuZ6SizJM4yGwAoCl7PHxsiRp3VB7ARNTckD3cpJhal70PVxgSg4AItv2+YDnyOH2hm6m5IDu5SJgiluHaSjKMDElBwDaOuYDnvXDnU3JkWECOpeLgCl+Si7UMJFhAgBtG/NTcu1mmNYOUsMEdCsXAdNwXMAUTcmRYQIAbesww7SWKTmgaxkvK+CLvuOm5KJ1mCbKU5px5Z61CwDyaGslw9RewHRYCJienHBybI8CdCQXGaa4ou+CFTRUoI4JAH434bRzyq/B1M6ilZK0csC0ol8aK0m7WAcY6EhOAqbmGSapdvFKlhYAsHRtec5nl05aWZBZewGTJB0eliJ4YpxsPdCJzAKmGVeuZI2iWqVmVhSHJEl7Svt70i4AyKMHQsD0wpWdDduHD/nzfjvOlBzQicwCpujKt+HCoAoW34xVfT5g2j0z1pN2AUAeRQHTyR0HTD7D9PgYARPQicwCpuoq383rlyIr+4YlSc+Vxue9TQCQV1u6DJjWMSUHdCXDgKl1/ZIkrQ4B0+4SGSYAS1PZOT24p7spuSPC2k2PMyUHdCSzgCmaYoum3OKsDI+TYQKwVD2012msJD1vuWnNYPsF35J0zIg/7+G9ZJiATmQXMJWigGl4zuNWFckwAVjaNu2ckSSdcUDnQ/YJo/7ch/Y6lcpkmYB25T5gWsmUHIAlbtMunxU6Y03nQ/ZIv+mIIdN0WXp0HwET0K4MAyY/xba6b2TO4w7qXyFJ2jH93Ly3CQDy6M4oYDqg2NXznBjqn6J6KADJZRgw7ZPUOsN0yMAqSdJTU7vnvU0AkDf7S06bd5dVMOnU1d0N2S9Y4c+PrrgDkFxmAdOzCafkooDp6aln2QMJwJLzi51llZx0ysqCRvo7K/iORDVQP39mJo2mAUtKZgHTrmmfYVrdImAaLi7TaHG5Jl2pEmQBwFJxyw4f3Jy9tvvh+qwD/XNs3Fmm8BtoU2YB01NTz0qSDhlY3fLYw8Ix2yefmdc2AUDe3BwCpnPWdle/JEmHLC/omBHTvpJ0726m5YB2ZBIwlV25EjAdliBgOnb5YZKkh/Y/Oa/tAoA8eWp/WZt2ljVQkF51UPcBkyS95mD/PDdsZ1oOaEcmAdPO0j5NuZJW9w1rqBi/8W7kuKEQMI0TMAFYOr61fUZO0nmHFLWiy/qlyIWH90mSrnmipDJ1oUBimQRM2yd3SpIOTZBdkqSTh9ZJku7a98i8tQkA8sQ5p69uLUmS3rKuL7XnPevAgtYPm54Yd7r6CbJMQFKZBEwPjv1WknTM8kMTHX/KyHqtLA5p68QO/Wp8+3w2DQBy4arHS7pvd1mHLDO96fnpTMdJUsFMHz+xX5L0N/dNaus+apmAJLr62FI+YkRbxh6Xc05OUlnO/+f897KcnFN0r8rOqayyrn9moyTpjNFjkjXSinr9mlN11Y7b9fFtV+udh56nNf2jKshUkMnMFP2H3uCdRt6MH3K87tzlMybRTFM04VT53rP7XcfPs2daunvXjL70qM8uXXZyv5YV0/2Lu/iIPl31eEk37yjr9Jv268+P7NOpBxS1ZsC0rCgNFKRCGy/ZbusYP7AQWTdrG51y9wc6Pnlt/0p956QPJ6phkvzWKBf96rPaPrWz05cEkIHNp33mbufcafP9On3XjS26gpz3H9enfzplQGbphxjPTTldfMekbnyKaTksHTMXjnQ8HnUVMJnZXkkPtXnaSkmd7HPSyXmdnHOgpE7WL+hV++hTVSf9ok+dn9Ppecc750Y7eK22dDgeSfxO9PocKd996vQ8+tT563R6Xm/HI+dcx1+S7urgnCs6fK22z+vwnLb71OP20acu+kWfum5fz37/evU6/E7Qp5Reiz71vn09HY+yKPr+Xg/P6/S1OtGr9tGn7tCnzs/p5rw843eit+d0ajH+ntOn7s7r6XjU7ZTcXa4HtQm9RJ8WjsXYL/qU/9fptcXYL/q0MNCnet1mmK7o8vw8ok8Lx2LsF33K/+v02mLsF31aGOhTja4yTAAAAEtBZpvvAgAALBSJAiYze62ZPWRmj5jZh5s8Pmhm14bH7zCz9Wk3NG0J+vQBM3vQzDab2Y/N7Igs2tmOVn2qOe4CM3Nmlvu56SR9MrO3hH+rX5rZ13vdxnYl+N073MxuMbN7w+/f+Vm0sx1m9hUz22FmD8Q8bmb2udDnzWb20i5ei/GI8SgTi3E8khbfmDRv41GCS/CKkh6VdJSkAUn3Szqx4Zh3S/pyuP1WSdd2etleL74S9ukcSUPh9rsWQ5/CcaOSbpO0UdJpWbc7hX+nYyXdK2l1+Hlt1u1OoU9XSHpXuH2ipG1ZtztBv14t6aWSHoh5/HxJN8ov8vxySXfM4/vHeLQA+hSOYzxaGP1aUGPSfI1HSTJMZ0h6xDn3mHNuStI1kt7YcMwbJX0t3P6mpNfYfCxNm56WfXLO3eKcGw8/bpT0/B63sV1J/p0k6ROS/knSRC8b16EkffpLSV9wzj0rSc65HT1uY7uS9MlJWhFur5T0ZA/b1xHn3G2Sds1xyBsl/bvzNkpaZWbJNpOsx3jEeJSVxTgeSYtwTJqv8ShJwPQ8Sb+t+Xl7uK/pMc65kvzKm2sSPHdWkvSp1tvlo9E8a9knM3uJpHXOue/3smFdSPLvdJyk48zsZ2a20cxe27PWdSZJny6VdLGZbZf035Le05umzat2/+a6eR7Go+wxHi2M8UhammNSR+NRks13m30ya7y0LskxeZK4vWZ2saTTJJ09ry3q3px9MrOCpM9IeluvGpSCJP9OffJp8A3yn7pvN7OTnXO757ltnUrSpz+SdKVz7lNmdqak/wh9Wsjbyqc1RjAeMR5lZTGOR9LSHJM6GiOSZJi2S1pX8/PzNTsdVznGzPrkU3ZzpcOylqRPMrNzJX1U0hucc5M9alunWvVpVNLJkm41s23y87bfzXmhZdLfve8456adc1vl9xI7tkft60SSPr1d0nWS5Jz7haRl8ns6LWSJ/uZSeh7Go+wxHi2M8UhammNSZ+NRguKpPkmPSTpS1YKwkxqO+SvVF1lel3XRVwp9eol8IdyxWbc3rT41HH+r8l9kmeTf6bWSvhZuHyifZl2Tddu77NONkt4Wbr8g/CFb1m1P0Lf1ii+yfJ3qiyw3zeP7x3i0APrUcDzjUb77teDGpPkYj5K+8PmSHg5/sB8N910m/0lH8tHmNyQ9ImmTpKOyfrNS6NP/SHpa0n3h67tZt7nbPjUcm/sBKuG/k0n6tKQHJW2R9Nas25xCn06U9LMwcN0n6bys25ygT1dLekrStPynt7dLukTSJTX/Tl8Ifd7Sze8e4xHjUV77tBDHo4T9WlBj0nyNR6z0DQAA0AIrfQMAALRAwAQAANACARMAAEALBEwAAAAtEDABAAC0QMCUMTP7aNjVerOZ3WdmL+vgOdbH7co8xzlXmtkFMY990Mx+bWYPmNn9Zvan4f5bzeyJ2n25zOzbZravph37Qz/uN7Ofm9nx7fYHQPqW4lgTjvvj9nqZDTP7SNZtwNwImDIUlph/vaSXOudeJOlc1e9vk0WbLpH0+5LOcM6dLL/rc+01SJ1pAAAGi0lEQVQy8rslvTIcu0pS44aFjzrnXuycO0V+A9SOBwEzK3Z6bp5fC+i1JTzWrJeUi4ApwRjT9lgZVrJHjxAwZetQSc+4sM2Bc+4Z59yTkmRmp4dPTfeb2SYzGw2flm43s3vC1ysan9DMimb2z2Z2Z/gk+c5wv5nZ5Wb2oJn9l6S1MW36iKR3O+f2hDY955z7Ws3j18ivnixJb5L0rTn6t0LSs03auMHMbjOzG0J7vhz2lpKZ7TOzy8zsDklnmtmpZvYTM7vbzH5oYUdpM3tvOHezmV0T7js7fOK8z8zuDe/ZBjP7fs1rX25mbwu3t5nZ35rZTyW92cyONrMfhNe63cxOmKNvwEKyJMcaSZ+U9KowJrx/jjZvCOPMdWb2sJl90swuCu/HFjM7Ohx3ZRivbg/Hvb7Fe7HBzG4xs6/LL5AYZcruNp/te0e475OSlod2XmUNmTzzmbhLw+1bzezvzewnkt5nZgeZ2fXhte80s1fO8T6hG1mvyLmUvySNyK+a+rCkL0o6O9w/IL9U/enh5xXyy9cPSVoW7jtW0l3h9nqFJeAlvUPSx8LtQUl3yS95/yZJN0kqSjpM/tPbBQ3tGZX07BztvVXSyyRtDs/zo/Da+2rasT/06VH5lVYPb/I8GyRNSDoqPM9NUVvkN0B8S7jdL+nnkg4KP18o6Svh9pOSBsPtVeH79yS9sua97Quv9f2a175c1SX+t0n6UM1jP1bYeiL08+asf0f44iuNryU+1tT+/ce1eUNo56Hh/v+V9HfhuPdJ+my4faWkH8gnG46VX0V6WYvnHZN0ZE0bDgjfl0t6QGHrlKhvje9z+PmDki6teW++WPPY1yWdFW4fLulXWf++LdYv0nkZcs7tM7NTJb1K0jmSrjWzD0u6W9JTzrk7w3F7JMnMhiVdbmYvljQj6bgmT3uepBdZtWZgpfwf9qslXe2cm5H0pJnd3ORcU+sdm2ck/VQ+eFnunNtmVptF92ny0N4LJV0hv79So03OucfCcVdLOkvSN8PzXx+OOV5+g86bwmsU5QdGyQ+kV5nZtyV9O9z3M0mfNrOrJH3LObe9oW3NXBvaMCLpFZK+UXPOYKuTgYVgiY81Sdo8JelO59xT4fkelQ/SJJ8ZOqfmOa5zzpUl/cbMHpN0Qovn3eT8RryR95rZH4bb68JxO1u0u9G1NbfPlXRizXuzwsxGnXN723xOtEDAlLEwqNwqv2v3Fkl/JukeNR9M3i+/n9Qp8p9wJpocY5Le45z7Yd2dZufHPGdtW/aY2ZiZHRUFMzGukXSDpEvnej5J35X01biXi/l5Irwnku/LL51zZzY5/3XyA/MbJH3czE5yzn0yTAGcL2mj+d3dS6qfel7W8Dxj4XtB0u5oAAYWmyU81tQ1T83bvEHSZM1d5Zqfy6r/f2WzsWuu5x1r+PlcSWc658bN7FbNHpOk5OOWwnFnOuf2N3kepIgapgyZ2fFmdmzNXS+W9LikX0s6zMxOD8eNmi/uWyn/abAs6U/kMy6NfijpXWbWH849LnxavE3SW8Nc+6Gq/8RU6x8kfcHMVoTzV0Tz7DVuD8dd3aKLZ8mny5s5w8yONF+7dKH8J8lGD0k6yHzBqsys38xOCuesc87dIulDklZJGjGzo51zW5xz/yifEj9B/v080cwGzWylpNc0a0z4ZL3VzN4cXsvM7JQW/QMWhCU81uyVn/5r1eZ2vNnMCqGu6Sj5cSrp866Un4ocN18j+fKax6aj8+WD1bVmtsbMBuUL9uP8SNJfRz+ErCDmARmmbI1I+rz5K0BK8rurv8M5NxVSzJ83s+Xyc/XnytceXB/+p36L6j9lRP5Nfv77HvM52t9J+gP5T2m/J59efljST2La9KXQrjvNbFp+t+dP1R7gnHOS/iXm/KPN7D75T1xTkv4i5rhfyBdkvlB+gL2h8YDwPlwg6XMh2OmT9NnQ/v8M95mkzzjndpvZJ8zsHPlU/oOSbnTOTZrZdfJTeL+RdG9MeyTpIklfMrOPyddPXSO/Ozew0C3VsWazpJKZ3S9ff/SvMW1ux0PyfTpY0iXOuQkzi3svGv1A0iVmtjk8z8aax66QtNnM7nHOXWRml0m6Q9JW+cA2znvlA8/N8mPkbZIuabNPSMD87yPQOyEt/UHn3FyfmgAgV8zsSvki8m9m3Rb0HlNyAAAALZBhAgAAaIEMEwAAQAsETAAAAC0QMAEAALRAwAQAANACARMAAEALBEwAAAAt/D8RuGR77lnGkAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "colors = [\n",
    "    \"steelblue\",\n",
    "    \"#EA7D1A\",\n",
    "    \"red\", \n",
    "    \"gold\", \n",
    "    \"#2ecc71\",\n",
    "    \"#03a9f4\"\n",
    "]\n",
    "predict_label = [\n",
    "    \"Water radial fraction\", \n",
    "    \"Mantle radial fraction\",\n",
    "    \"Core radial fraction\", \n",
    "    \"Core mass fraction\", \n",
    "    \"Scaled CMB pressure\", \n",
    "    \"Scaled CMB temperature\"\n",
    "]\n",
    "\n",
    "\n",
    "y_label = np.arange(0,1,0.001).reshape(-1,1)\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "fig.subplots_adjust(hspace=0.3, wspace=0.1)\n",
    "for i in range(OUTPUT_DIMS):\n",
    "    ax = fig.add_subplot(3, 2, i+1)\n",
    "    # y_label = y_labels[:, i].reshape(-1,1)\n",
    "    # for i, label in zip(range(OUTPUT_DIMS), predict_label):\n",
    "#     fig, ax = plt.subplots(figsize=(8,4))\n",
    "    # label = \"Water radial fraction\"\n",
    "    ###################Earth###########################\n",
    "    mus_ = np.array(locals()['mus'+str(i)])\n",
    "    sigs_ = np.array(locals()['sigs'+str(i)])\n",
    "    factors = 1 / math.sqrt(2*math.pi) / sigs_\n",
    "    exponent = np.exp(-1/2*np.square((y_label-mus_)/sigs_))\n",
    "    GMM_PDF = np.sum(pis[0]*factors*exponent, axis=1); # 对多个高斯分布求和\n",
    "    plt.plot(\n",
    "        y_label, \n",
    "        GMM_PDF,\n",
    "        color=colors[i],\n",
    "#         label=label,\n",
    "        lw=2,\n",
    "        zorder=10,\n",
    "    )\n",
    "    ax.set_xlim(0, 1)\n",
    "    ax.set_ylim(bottom=0)\n",
    "    ax.set_yticks([])\n",
    "    ax.xaxis.set_minor_locator(tck.AutoMinorLocator())\n",
    "    ax.set_xlabel(predict_label[i])\n",
    "# plt.savefig('512EARTH-water.eps')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "gYMd49hkiWsM"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "MDN_prediction.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}