aalr.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run ../src/aalr_inference.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = SimulationData.load(\"/scratch/folgert/synthetic.pkl\")\n",
    "dataset.params = ['b', 'mu']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 460,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_sample, output_sample = next(iter(dataset))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 461,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([-0.0460,  0.0065])"
      ]
     },
     "execution_count": 461,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "input_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 462,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([36.0000,  3.5176,  1.2644,  1.1000,  1.0680,  1.0564,  1.0506])"
      ]
     },
     "execution_count": 462,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "output_sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 682,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2,) (7,)\n"
     ]
    }
   ],
   "source": [
    "from hypothesis.nn import ConditionalMLPRatioEstimator as MLPRatioEstimator\n",
    "\n",
    "activation = torch.nn.ELU\n",
    "layers = [64, 64, 64]\n",
    "inputs_shape = (input_sample.shape[0],)\n",
    "outputs_shape = (output_sample.shape[0],)\n",
    "print(inputs_shape, outputs_shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 683,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 100/100 [06:46<00:00,  4.07s/it]\n"
     ]
    }
   ],
   "source": [
    "import tqdm\n",
    "\n",
    "from torch.utils.data import DataLoader\n",
    "\n",
    "from hypothesis.nn.conditional_ratio_estimator import ConditionalRatioEstimatorCriterion\n",
    "from hypothesis.visualization.util import make_square\n",
    "\n",
    "batch_size = 256\n",
    "epochs = 100\n",
    "\n",
    "ratio_estimator = MLPRatioEstimator(\n",
    "    inputs_shape, outputs_shape, layers=layers, activation=activation)\n",
    "ratio_estimator = ratio_estimator.to(hypothesis.accelerator)\n",
    "ratio_estimator.train()\n",
    "criterion = ConditionalRatioEstimatorCriterion(\n",
    "    ratio_estimator, batch_size).to(hypothesis.accelerator)\n",
    "optimizer = torch.optim.Adam(ratio_estimator.parameters())\n",
    "\n",
    "losses = []\n",
    "for epoch in tqdm.trange(epochs):\n",
    "    data_loader = loader = DataLoader(\n",
    "        dataset, num_workers=4, batch_size=batch_size, drop_last=True)\n",
    "    num_batches = len(data_loader)\n",
    "    data_loader = iter(data_loader)\n",
    "    for batch_index in range(num_batches):\n",
    "        optimizer.zero_grad()\n",
    "        inputs, outputs = next(data_loader)\n",
    "        inputs = inputs.to(hypothesis.accelerator)\n",
    "        outputs = outputs.to(hypothesis.accelerator)\n",
    "        loss = criterion(inputs, outputs)\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "        losses.append(loss.item())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 684,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Logarithmic loss')"
      ]
     },
     "execution_count": 684,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEGCAYAAABPdROvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2deZQV1bX/v7vpblrmUWZoEJTBgaGFoGKMioKCAypOOBCVRBNxNiZ5y2f8rfeSEKeARsUBHKJixJGYZ0QQRGVokBmBZm4gzTzTNN29f3/UwLnVVXXr3qp7697q/Vmr1q3h1Dm76ladXeecffYmZoYgCIIgJEtO2AIIgiAI2Y0oEkEQBMEXokgEQRAEX4giEQRBEHwhikQQBEHwRW7YAgRNixYtuLCwMGwxBEEQsoqFCxfuYuaWyZwbOUVSWFiI4uLisMUQBEHIKohoU7LnRqZri4iGE9HE/fv3hy2KIAhCrSIyioSZP2PmMY0bNw5bFEEQhFpFZBSJIAiCEA6iSARBEARfiCIRBEEQfCGKRBAEQfBFZBSJWG0JgiCEQ2QUiVhtuSPhAgRBSBWRUSSCOy+99FLYIgiCEFFCUyRE1IGIZhLRKiJaQUT32aQhIhpPRCVEtJSI+oYhaxQoKysLWwRBECJKmC5SKgE8xMyLiKghgIVE9CUzr1TSDAXQTV8GAHhR/xUEQRAyhNBaJMy8nZkX6esHAawC0M6S7EoAb7LGXABNiKhNmkUVBEEQXMiIMRIiKgTQB8A8y6F2ALYo26WoqWxARGOIqJiIinfu3JkqMQVBEAQbQlckRNQAwFQA9zPzAethm1NqmB8x80RmLmLmopYtk/KCLAiCICRJqIqEiPKgKZG/M/OHNklKAXRQttsD2JYO2QRBEARvhGm1RQBeA7CKmZ9xSPYpgFt1662fANjPzNvTJqQgCIIQlzCtts4FcAuAZUS0WN/3OwAdAYCZXwLwOYDLAJQAOAJgdAhyCoIgCC6EpkiYeQ7sx0DUNAzgV+mRSBAEQUiG0AfbgyIqvraqqqrw9ttvhy2GIAiCZyKjSKLia6u6uhobN24MWwxBEATPREaRCIIgCOEgikQQBEHwhSgSQRAEwReiSARBEARfREaRRMVqSxAEIduIjCKJitWWIAhCthEZRSIIgiCEgygSAQCwdu1afPvtt2GLIQhCFiKKRAAAHDhwALt37w5bDEEQshBRJIIgCIIvRJEIgiAIvoiMIhHzX39o4WEEQRASJzKKRMx/BUEQwiEyikQQBEEIB1EkgiAIgi9EkWQps2fPDlsE31RWVoYtgiAIASCKJEuZMWNG4HlqkY3Tw969ezFx4sS0lScIQuqIjCIRq63sgpmlRSIIESEyikSstgRBEMIhMopEEARBCIdQFQkRvU5EO4houcPxC4hoPxEt1pfH0y1jOnjnnXfCFkHIMsrKysIWQRBMwm6RTAYwJE6ab5i5t748mQaZ0s6aNWvCFkHIMl588cWwRRAEk1AVCTPPBrAnTBkEQRAEf4TdIvHCQCJaQkT/IqJedgmIaAwRFRNR8c6dO9MtnyAIQq0m0xXJIgCdmPksABMAfGyXiJknMnMRMxe1bNkyrQJGBXHaKAhCsmS0ImHmA8x8SF//HEAeEbUIWaxIks7JiLUNubdC1MloRUJErUn/VCai/tDklTB+EaC2VK4HDx7E5MmTwxZD8IBMZk6esM1/3wXwPYDTiKiUiO4gol8S0S/1JNcCWE5ESwCMB3AD15YaKM1I11ZqqKqqwr59+8IWQ/DAs88+G7YIWUtumIUz841xjj8P4Pk0iSMIQpo4cuQIpk2bhpEjR4YtihAAGd21lQjia8s/0tgT0kVlZSW2bNkSthhCQERGkYivLUEQhHCIjCIRBMEbH3zwQdgiCBFDFImQNNXV1dIdloUsX27r2k4QkkYUSQaSLRZUs2fPxoIFC8IWQ8hCsuUZF7whikRBvq4To6KiwldwKqlMBCEaREaRBGG1NX78+AAlEtwQpS0I0SEyiiQIq629e/cGKFH0EWUgCAIQIUWSSqqqqsIWQRAEIWMRReKBv/zlL2GLIAiCkLGIIvFAeXl52CJkJMxsO2AuM5YFL0jXaHQQRSIACNaC6rXXXgssL0EQMp/IKBLxtSUIghAOkVEk4mvLP8l0NchckPjIPRKiTmQUiSAIghAOokhCYPv27Th8+HDYYtQg0S9nGSwVBAEQRRIKs2fPxubNm8MWowaiGARBSAZRJAFQXV2No0ePek4fpQrbT/9/FMYOqqurwxYh4zh27BiOHDkSthhCGhFFEgClpaV45513whYj7URJISbLU089FbYIGccPP/yAWbNmhS1GVlJeXo6PPvoobDESJiFFQkRNiejMVAnjhzDNf4kooUrVaSKfkH3Il3dyyPNvz/Hjx7Fu3bqwxUiYuIqEiL4mokZE1AzAEgCTiOiZ1IuWGEGZ/0pXhSCknkxtzWaqXJmOlxZJY2Y+AGAEgEnM3A/AxakVKzyOHTuW8Dlqi+Tll18OWqRIUpte2Np0rULtxIsiySWiNgBGApgWZOFE9DoR7SAi29ifpDGeiEqIaCkR9Q2yfDv8BGoCNNPeeESlaysq1yGkn9r+3OzYsSNS40heFMmTAL4AUMLMC4ioC4C1AZU/GcAQl+NDAXTTlzEAXgyoXEc2btyY6iISRtzYp5+VK1di9+7dYYuBhQsXhi1CwtR2JeGFgwcPYtOmTWGLERhxFQkz/4OZz2Tme/Tt9cx8TRCFM/NsAHtcklwJ4E3WmAugid46Sgn33ntvUi9u0IPt1rzGjRuXsEyJkug1RJ1ly5Zhx44dYYuBzz77LGwRag2iAJPHy2D7OH2wPY+IviKiXUQ0Kh3CAWgHQPVJXqrvs8o4hoiKiah4586dSRfWrFmzjGyRJDNuky7k5ROSRT5cooOXrq1L9MH2YdAq8lMBPJJSqU5gV0vVePqYeSIzFzFzUcuWLZMvTK8Ujx8/nvB5ib4UUgELUSZblUQ65Y5SHeBFkeTpv5cBeJeZ3bqigqYUQAdluz2AbakuNJG5AfEevF27dqG4uDihc7IFv9cRpRdJEGozXhTJZ0T0I4AiAF8RUUsA6QoZ+CmAW3XrrZ8A2M/M8c2ifLJkyRLPaZ95xn1Kze7du7FmzRq/IglC5EjXh0QmGqtE5WPSwMtg+2MABgIoYubjAA5DGwT3DRG9C+B7AKcRUSkR3UFEvySiX+pJPgewHkAJgFcA3BNEufGYMWOGZ6eKBw8ejMTM9mTlybTrCIogX/So3iO/pKsyFTc2qSc3XgIiygNwC4Dz9RdiFoCXgiicmW+Mc5wB/CqIshLl9ddfx5133on27dv7ykdVGpmoQFQSfbGj9lVV21m3bh1atWqFBg0ahC1KoHh1qJrJ72am46Vr60UA/QD8TV/6Ig3zORIlKF9b6kv06quvei0bVVVV2LVrV41jqvIYP358jX2CkCl8//33nibUxkOe7dqHF0VyNjPfxswz9GU0gLNTLViiBOVra8SIETHbXnxvlZWVYefOnXj++edtjxsv1t69e33Jlmns3LkzkIonahw8eBDz588PWwwhCcRqKzm8KJIqIjrF2NBntmfe6FVAdOnSJWb7ySefjHuO28Q1pwczjIco6Jdk8+bNWLFiRdLn+5nzk8kcOHAAixcvDluMjCZKlajgTZE8AmCm7gV4FoAZAB5KrViZxxNPPOF47Pvvv3c8ZteNlY6vnkQtxeTFFoJ8BuI947V9fC1q1+/FausraL6uxurLacw8M9WCZTKTJ0+O2Y7X/WX3gqa64k4m0FbUHu6wyGalLM+AkAyOVltENMLh0Cm6ueuHKZIpdMaMGYOJEyc6Hre6UXEz/w2rRZLpZMM98KMQsuH6hFiy+QMgbNzMf4e7HGMAGaVIiGg4gOFdu3b1nVfbtm0djyUa+CqZCkUeaCEM1q5dizPOOCMtZWXqMy4fAMnhqEh066ysgZk/A/BZUVHRXakqY9WqVZgyZUrC59m1SJxepLAe5Ex9sbORbL6XEiE0fWTzc2IloZjttR1ViaimvG6Vf9S/cKL0Mqj4/d+M85kZFRUVQYiUNUT1mRCcEUWSJOvXr/eUzmmMRF622sGOHTsScgJam4j6R5YbUbt2USRJ4nUOhFVpzJs3L+45qVIyYpKJtE8UrA33VBC8BLb6FRE1UbabElFanCdmMnPnzvUc333v3r2YOnUqAOBf//qXa1o36y/BP59//rnntH4cWdb2/8vL9Wdaq9yLl+ADBw7IOJINXlokdzHzPmODmfcCSNmAdrIE5WvLoKioKG4aJxPhw4cP49133wWgvVDHjx+Pmf0uXVuiGIXEWblyZcrLiPdcvv322yn3yJCNdYMXRZJDypURUR0A+akTKTmC8rVlMGDAgLhpnFyjrFy5EqtXrzbkqvFg7Nq1C/v27bM7NTTkKzqWLVu2ePYaa6WiogJlZWUBSxQ9En3e3n///RRJkllk43voRZF8AeB9IrqIiC4E8C6A/0utWOHjJ2SvWgEZikR9OI4ePYpVq1Y5nm/3RZKND1c2c+TIkaSVwX/+85+ApRGynYMHD4YtQkrxokh+A82/1t3QYoN8BeDRVAoVJY4fP46tW7fW2B/GPJJkzVDffPPNgCURxSjULp5++umwRUgpcQNbMXM1tPgjGReDJNX06NHDteXghDrHxOlLJNF+0CAqXq/GAVa8mjoLGtnYxx0ktf36vRKl++TYIiGi9/XfZUS01LqkT8TwuP7665M674cffgCgxWs3Bv+tiiCVD1E6rUqi9DKESZTilyxevBjFxcVhiyGkEbeurfv032HQ/G5Zl4wiaKutIFiwYIGjpYkxbmIMuquKxu8Yye7duxOUNHmFkGxLKRu6tlKpJJkZX3zxhbmdiFlypnPo0KGwRUiadJktZ8PznwiOioSZt+u/m5h5E4C9AA4qS0YRtNVWULi1DpgZf/3rX2vsKy8vt03rlXQrBSE53OLY1Aaydca/vCc18TIh8RdEVAZgKYCF+lJr2q2/+93vkj5XfeCsD19lZSUWLVokXUNC5JBnuvbhxWrrYQC9mLmQmTvrS5e4Z3mAiIYQ0WoiKiGix2yO305EO4losb7cGUS5iZCfn2/IkvC5bk38Y8eOYebME/HBDh8+7JpXpn4FZUul4WXWclBkyz0RkkP+35p4USTrAATeBtUnNr4AYCiAngBuJKKeNkmnMHNvfXk1aDm8cO211+Lxxx9P+LxE4pkbprmTJk2yPR5E11adOnU855Fq0q0Yn3rqqYTPWbJkSQokic/06dNTXsauXbtSXobgTpQUkhdF8lsA3xHRy0Q03lgCKLs/gBJmXs/MFQDeA3BlAPkGzumnnw4iQv369ZPOw+7FraysNB8mw8zYabZ8srz22mvmek6O89+dDQ/10aNHceDAgaTPTZR4rcRUMWfOnJSXkaktXCE78aJIXoY2IXEuToyRLAyg7HYAtijbpfo+K9foJscfEFGHAMpNmkceeQRdugTSqwdA694yiGeyG+/F3759u7muKoUtW7bYJU+qjCDYs2dP0ucuW7YM33zzTYDSCEIs6VKwUVPkXhRJJTM/yMyTmPkNYwmgbLtPYOvd/QxAITOfCWA6ANtyiWgMERUTUXGqHaoFjVHp+5378fLLLwchTsoZP15rzP7f//nzspPJL6LX1l2YrcBsaIFmIpn83IWJF0UyU6+o2xBRM2MJoOxSAGoLoz2AbWoCZt7NzMZn+ysA+tllxMwTmbmImYv8+Mjywumnnx5ofl4VSTY/wHayb9682VeeVrNpITFEkQhB4kWR3AR9nATBmv8uANCNiDoTUT6AGwB8qiYgojbK5hUAEvdXEjB9+/b1ZRLshF1la7U0On78eNx8Jk6cmNZKwstcgCArfePaMs178tatW7F8+fIa+7PtA6CkpCRsEYQsJK4iUUx+1cX3QAEzVwL4NTTvwqsAvM/MK4joSSK6Qk82lohWENESAGMB3O633CDIz8/H6NGjA8nLrUUybtw4PPHEEwC0CunPf/5z3Py2bdvmeCys2OHprPRV4wIgfe5iysrKssonmdPHxsaNG9MrSIahjlummii1Cj2F2iWic4joJiK61ViCKJyZP2fmU5n5FGb+H33f48z8qb7+W2buxcxnMfPPmPnHIMoNgk6dOuGMM84ILL9169bVeInVh5qZazhdTPRr1+1rk4gyPvKb0/VWVFRgzZo1AGKNCyorK9PWBZZtLQ+/lVg65+Wkk9LS0rBFyEoF42Vm+1sAngJwHoCz9SV++MBagJ8/3LBeMvLYvn07Nm3a5Jje7qveMCxwcghpBAIy4moY6Zy6yIzKOBNgZs+V8/79+/Hvf//bNo94vtcWL17sqcswaNKheBYvXpzwOV7lSmZeTjZQr169tJTjdp+z7aME8NYiKQJwLjPfw8z36svYVAuWKJnotNGNRLt77PxvGftefNHew7/hMNLa5TJu3LiYbbsHd8GCBY7zKFavXh23gvbLd999Z3pRThYv80C+/PLLtHZnpJOPP/7Y8Zjfr95ko0caZGplmalyZTpeFMlyAK1TLYhfMtVpYzzUF9rt5bbrdjL2GZMY582b51qGoVCsX+Cff/55DWWzYMGCGBcv6vEFCxYEPnHSyuHDh20rK7t7lK4wwfHKMGTLhq6JVMqYDdcfNm73KBu9J7vFI/mMiD4F0ALASiL6gog+NZb0iZi5MDOGD/fnUT/RClANlGU91+iaWrRoke25a9euNdd/+OEHHDx4EBs2bMChQ4fMr/KjR4+aCkr1Tmt0xW3YsCEmcFcqmDdvXox7mcrKyoz4UoySyXGyQc68kAn/VSbwz3/+0/EYM4dm/JIK3FokTwF4GsATAK4C8L/6trHUepo2bQq/LSCnbhVr3O+6desC0EJ2GnMw7Lq7AGD27Nlxy1q4cCH27duHN96IneNZVVWFp59+Gjt27IjpY2dmrF+/HosWLcLu3bsD/epk5pgYKhs2bIhRmNOmTcOKFSscK6jdu3fHnH/48GG88sornspNBLU70to1mY7K8+233w4sLzfrvtpMkP/jggULHI9VVlZmlZVfPNzikcxi5lkALjPW1X3pEzFzufDCC9G1a9fA8psxY4a5bp2wpyqWf/zjHwCAKVOmOLpGUVH321nbqF1DROQ4tuAWu33fvn1Juz+prKyMGef58ccfY7ry3Fok06dPx1dffVUjv61bt8Yt108L47nnnjPXnRR60GTLHI+oWnSFYZSRLXgZIxlss29o0IIIsVgj5qmmwerXutrK8FLRGNY2iVb6RkXuNCN92bJlMV1qpaWl2Lp1a0x640VUIwMmCjOb3V6VlZWYM2eO43jNhAkTXPPat29fjNL88MMPk5Ipld1E6SaIL/JkHWtmAm7X72UeVxDlZCNuYyR3E9EyAKdZ4rVvgBbkKqPINqutRHEK2etmmWOHMYCtzkhXWyTxKkXj/sbr2lq4cCFee+21mIraqPCtTXpmBhFh3rx5NUyg1Upp48aNmD9/vtkic+oaMK4l0S/jpUuTe6zV1ly2E7UKLlHcrj9KHwxB49YieQdabPZPERurvR8zj0qDbAkRttVWt27dQinXz6xxo+JTrUSeffZZ27TWF8zOLT4z491333U8J16+paWlMWbF5eXleOaZZ8w01i6kd955xzbPuXPnupZpHT/Ztm2bGQPES/eF3XUlUwFbFU8QXSeGko1HsgqjuLjYdXwl3r1PFStWrEjbPKgofDAEjZsiYWbeCOBXsMRqD8hpY2QoKioyIylmIk6VhvFFv2nTJsc0xv6FC2MjB6iKZPv27WBmfPvtt1i9erV5ntHSAJxdlfzxj3/E/PnzbV9Ow2X8ihUrMHXqVM/uO6yVmfXarOMnpaWlZgyQeLP7P/nkE/zhD3+wzccg2YrGr0dkwD2YGjPH7Xayym69Hxs2bLDtFi0vLwcz17iGdevW1RjDcpMvWbZv355yk3TBmXgtEuCEk0Y1FkmtidnuhWHDhqFRo0Zhi+GIUwWlDu7bORwENCsxoGbALbV1sH///pi8gJqzqseNG2d+yaoVRkVFBaZPn46Kigps3rw5JV97Ritr/vz55j7VSCEec+bMMU01VaOHKVOmANCuZ8uWLSguLsayZcsSls+4H9XV1SmbQ7Bnzx4cOnQorjWbtTL3OoP9ueees7VAPHLkiGkuvmfPnhiTciD2fk6YMCHp7qOguuRqe9desrhZbQ0j7a3+KTN3CdppY9SI6gPoVLGpFeZ7773neL5ReZeXl+O7774DgBhTXZX9+/fjyy+/TFZUxyBexle4asCgxm+JNzYyd+5c20rSyJeZTQOIqVOnxihDt+fC+No3xoWICKmKp2PEgfEil4oX785e89y/f7/ZYjUw3PcAwQ7SJ+JiJ0ysZv7ZiqvVFmv/xEdpkiWrSXUclGxFHRA3vkzdBsH9VCbJRk9Ux5mqq6vx1ltvxRw/dOgQJk6cCMC+JVNdXY0NGzbY5m0Xf720tDRmPCQbZzJbCbolWVJS4mvM6Ouvv07K11i8rs2dO3fGfHgkq6zUVqhKqif7pgov5r9ziejslEvik7Cttvr1s425JaSRIAZbn376aaxbt67GftXk2oq1MjFaGAcPHqzh5mXTpk149dVXY7zMOlVGhsVbUPMX1HKcynTzoXX8+PG4YzB22CmZeAYL06ZNMz0uqB8XakX7n//8J6ZVqJZz9OjRlM0cVy0R//CHPwTqdy5b/b55USQ/A/A9Ea3TzX+XEVHGmf+GbbUFwPcsdyF8kumjt1aKhrNJu24qw4zbTllZMcalnGa0e5H1wIEDjuNf8VAV89KlS/GnP/3JMa2Tv7Py8nLzq3vGjBnYvHkzqqurMXPmTE8yrFy5EjNnzjTzULsn582b53gf165d69uxJKC1nq0m9tbrdLJ09IJVyWZDd5wdXhTJUACnALgQmvnvMP1XsHD//feHLYIQMGofvhN2LuwBrVIwKrMdO3YkbFpsuK9RFca2bdtMmdTxi9LSUtOyjpnN7sO9e/fWMISwlmO3DcSaV3/44Yc1uiQ//vhj8+vcqWvrn//8p9n62rlzJ6qrq/Hkk08mXMn/7W9/qyGrU/gEQLvueJNumbnGf2J3X5JVxAZLliyxzdvYF4VgYl4iJG5i5k0AjgJgZREsiH159HBy0a9iHUA2mDZtmrleVlYWMzPaMAxYtmwZpk6dau53c0MDaF/odl14e/bsMbvUNm3aFGNtZVSoVVVVZgVu+DMzDCXcuoGsEz+LizWjzcWLF8e4mVENGIATXXJG149q6bdqlRY1m4hqWBUePHjQPMeofFXF895778UoCaO1UlxcHGOSTURma85uXG7t2rXmvTcUYklJCaqqqmLy8dJKmDVrluOxTz75BADw/PPP2+Y5efLkuGVt3LjRd1iFVOIlsNUVRLQWwAYAswBsBPCvFMuVtVx//fVhiyBkCNaB08rKStPazVAkqhJRsVqSqUrpxx+1QKFEhP379+PIkSOYNm2aOY6jKgyV/fv3m60bo4VUVlZmW8mqHgasyk39gjYU0OHDh2sYSqiKw80f2dy5c3H48GFUVVVhwoQJWLJkCZjZtkuPiLB3714cP34cx48fx7FjxzBhwgQwM/bs2YPdu3ebrnqIyHQbpFb0mzZtwsyZM8HMphIyFFNZWRkOHDhgTux0G99RK321q85QtECsO6Ddu3e7tqKsqEYYkydPNpWbmsfGjRvN/2/Pnj2h+Tnz0rX1/wD8BMAaZu4M4CIA36ZUqiymR48eYYsgZCCGUvFqTqv6/dq6datZOc2ZMyfma/nbb7/F0qVLUVFRYVqOOVkeOQ22P/XUU/j6669j0k6aNMmTnF6J54HBmBCqDlxbPTqrbN++HStWrMC6detiFHZVVZXpvt3ODHvVqlWYNGkSZs2aZTuuY90+cuSIY+Xs1FpVlf73338fI8crr7yCQ4cOxSgbFdUzhHUOj5HPhAkTUF1djTlz5mDy5MmmjFOmTEmZ+Xg8vCiS48y8G0AOEeUw80wAvVMsV8KEbbUlCG7YjVNYcYoj45a+srKyRtfQ3//+d3PdqcvKGGfZu3cvjh49WmOiYLqI9wX9zTff2JrEGsrSiLFjbKsxQHJyTlRvRjeX2lJTPSA4dSnZuQJavHgxKioqbGfSu4XLBrRuu+rqas+x4auqqhxNhI2W6caNG/HBBx94yi9V5HpIs4+IGgCYDeDvRLQDQMZ5L2PmzwB8VlRUdFfYsghCOvBi+aQqFa/ejYPqHlEnQTpV1Mb4gdPxvXv3xiiSNWvWICcnp0YFb0QHraysRJ06dQDU7DYqKSlBz549ze3169ejQ4cOWLlyJerXrw9AM1pQx1+c5jXZWbAxc9yWnNdxVKPl+vXXX5tjcOq5JSUlpjJau3Zt6EGyvCiSKwGUA3gAwM0AGgN4MpVCZTvNmjVLOjaHIKQKt7kwKuPGjUsoX6eJoInMf1G9EqiRPO2orq42vSTYYSglp+4ja7lbtmxBq1atzH1vvfUWmjRpAgD46CP7+dh23Yeq4YFbEDI7T95bt25F06ZNzW3DoKKystJs+WzatAnTp08HM8d4k1i6dClOOeUUx/LSgRerrcPMXMXMlcz8BjOP17u6fENEQ4hoNRGVENFjNsfrEtEU/fg8IioMotxUc++994YtgiDE5YUXXrDdn+ikOC9OGV977TXX4+q4T7zuoWQwvuadYuFYzbz37dsX0xLyIpM6DqTGBlJbEtauSEMhzJw50zYsgqr8y8rKHJ2Lrlu3zpOpeqrwYrV1kIgOWJYtRPQRESXtc4uI6gB4Ado8lZ4AbiSinpZkdwDYy8xdATwLILjIMilEzICFbCCdFj7pjOVh11ow3slElKTawgjK+MBqbKGO6arjHMYkSL9zWNKFl8H2ZwA8AqAdgPYAHgbwCoD3ALzuo+z+AEqYeT0zV+j5XWlJcyUAI6j4BwAuoiyppS+99NKwRRAEQSeeD61U4kWJZksYZSe8KJIhzPwyMx9k5gPMPBFaHPcpAJrGO9mFdgBUd62l+j7bNMxcCWA/gObWjIhoDBEVE1FxWOZvVgYOHAgAuOqqq0KWRBCE2bNnhy1CpPGiSKqJaCQR5ejLSOWYn1AcNUoAAB6tSURBVBnudi0La35e0oCZJzJzETMXZZIX3l69eiE314s9gyAIQvbiRZHcDOAWADv05RYAo4joJAC/9lF2KYAOynZ7ANYYnmYaIsqFZjGWNeZQZ5xxBlq3bh1aGF5BEIR04MVqaz0zD2fmFvoynJlLmPkoM8/xUfYCAN2IqDMR5QO4AVp8eJVPAdymr18LYAZnkXvM7t27o0WLFrj55psBAIMHDw5ZIkEQokxY1aMXq632uoXWDiIqI6KpRNTeb8H6mMevAXwBYBWA95l5BRE9SURX6MleA9CciEoAPAigholwNuFkJ5DJYXoFQcgevM4VChovHfiToMVvv07fHqXv8/15zcyfA/jcsu9xZb1cKTey1KtXL9Awo4Ig1E5+/PFHnHrqqWkv18sYSUtmnqRPSKxk5skAMmdEO0to0KABcnJyMGrUqBotkCzqrRMEQaiBF0Wyi4hGEVEdfRkFIJCZ7UGS6U4be/XqBSJC165d8cADD8QcMxTJySefnFCeDRs2DEw+QRCEZPGiSH4OYCSA/wDYDm3Qe3QqhUqGTAi164ba6rCOlRhO5gxat27tKU8xLRYEQSWsD+m4NREzbwZwhbqPiO4H8FyqhIoip556KurVq2d7jIjQoEEDc9trV5d0iQmCoOIUwz7VeGmR2PFgoFLUArp27Yq2bdua23379o05brRSbrrpphqePHv16pV6AQVBEJIkWUWSFf6uMpkrrrjCVB7MDGYGESEnJweXXHJJ3PPHjh3rqmDOOOOMGvtUN9VBYcRxEASh9pKsIpE+lQAoLCxEly5d0KRJEzRt2jQmoluLFi3MdbvKumnTpujQoUON/QZ2c1buvvtunxLX5Lzzzgs8T0EQsgvHMRIiOgh7hUEATkqZRElCRMMBDO/atWvYonjmtttuw8KFC7F161aMHDkSmzdvNhWIqlTq1q1re3737t0TKi+e4+S8vLyEghEJgiAALi0SZm7IzI1slobMnHHmQplutRUPIkKnTp3M1oc6htKsWTN07doVZ511lufuKetA/MCBA+MqkgsuuCAxoQVBEJB815aQYoxuq65du6JBgwYYNWoUAK07zAtdunTBSSedaDh6MSlOJtRLtipuQRCCQxRJyBiD7F6OM3OM5ZdB06ZNYwbeGzZsiD59+tSYZ5KKmGBeFZsgCNFFFEmWYacMioqK0KRJkxr7L7/88hhTYuvERwAxrRY1b7v87HCaGyMIQu0h48Y6ahunnXYaOnXqVGO/0fqwtkjy8/Nt0wLARRddhPLycuzZo4VsqVOnjnlus2bN4sqiKoUwQ5MKgpBdRKZFkum+tpxo2LAhnKI6EhHy8vJMC67WrVujUaNGZnwTg7y8POTl5WHQoEFo0aJFjJWXoUhUU+HRo+093BhKKjc3Ny2z5u0UqCAI2UdkFEm2W205cd1115njEOeeey4KCwvRrVs3U0E89thj6N+/P84///wa5zopg6ZNm5oOH++44w5zv5HnddddF3NuqrqvbrvttviJBEHIeCKjSKKG0aWVm5vrOkheUFBgzoi3kpubazsHRe0us06EVNMYqF6G1fyM8x58MDmPOXblCoKQfcibnKHk5ub6bgl07twZI0aMwJAhQ1zTqd1fF198MerXrx8zs14dpLdTatb4KpkcUjhRV/2CIMRHFEmG0qZNG9x4440Jn9ewYUM0b94cAMyWyk9+8hPH9GrLo379+jjvvPPQvn173H777bbpvZgQ9+nTJ2a7fXvnyMw33XRT3PyCJC8vz3PaVJhLC0IqCcv3nSiSDCaZiqxr164YNGiQaxrr3BU7s2D1uNoyUitiO/nq1atnKidjvOrOO+9Ely5dbGVxCgt6//33x2x37tw5ZtuYYOlkqOCEVeazzz7bMa3XuDCCUNuJjCLJVqutMMjNzTUrydzcXDz88MOu6VUrsVtuucVct5scaS3HDievxar7F1V5FRYWmgPzl19+uXn8nHPOcS1fxVAYdkrTiTFjxnhOKwi1mcgokqhabQVNnTp10KBBg5huM7u5KSrqV7zaAlAtvuKd5wXVv5h67lVXXWWuG8qJiGLGZoYOHYq8vLyYsR1Aa6EBJ5SeVZGcdtppgckvCLWVyCgSwRu/+c1vUpJv7969AQD33ntvXCMBvxW02hIyutGMPE899VRcfPHFjuf26NEjph9Zndk/bNgw23PUNIKQyYT18ROKIiGiZkT0JRGt1X9tXdoSURURLdaXT9MtZxSJ1/pwIt7Av9FqOOmkk2Ie5p49ewLQKnzDu7B1MN7AOu9l7NixAGJfDmsaqyIBTpgV9+/fv8YYit2YiKGY1PGakSNH2sooCJmMNbpqugirRfIYgK+YuRuAr/RtO44yc299ucIhjZBCjHkjbl1AThARfvazn5nbHTt2BKA97Pfcc4+5v3///uZ68+bNzbkzzZo1i/Eu3KpVqxr5A8CZZ57paJasdmUZraYBAwaYlm1OGAoQiO1aU/n973/vmocgpJtErBKDJCxFciWAN/T1NwDYv6lCqDRp0gSPPPKIuX3ttdcCABo0aOA5D2sLwlAM6nyOyy67DIBWwd91110x58Rz1WJ3PN757dq1c3R/r06+NHBSomG9tIKQaYSlSFox83YA0H+dZokVEFExEc0lIkdlQ0Rj9HTFO3fuTIW8tZL77rsvxvLq9NNPB+A8k/2nP/2puW419x08eDDatm1bw6xXpW7duigoKPAkm7UFEpRvMLc+5p49e5rX+MQTTwRSniBEgZQpEiKaTkTLbZYrE8imIzMXAbgJwHNEZNsByMwTmbmImYsSnVdQm4l3r5wqVSfXJmo31q233hpzrG3btqY7FzusA91O6VSFoaaxpn/ooYdsz7emJyLPSmjkyJEx1+gVcQUjRJ2UPeHMfDEzn26zfAKgjIjaAID+u8Mhj23673oAXwOwH6UVkkIdp0gFiViQFBUVmevWil2t9K35Gmmt51i7qFJtzXL11Vc7HnvsMachQEGIBmF9Kn0KwHD9ehuAT6wJiKgpEdXV11sAOBfAyrRJWAtQK1cnSyo/9OnTx7e/MCcrrSFDhpiD96lUEkZ3nhvdunXDWWed5Xg8WUs5QUiUWmX+C+BPAAYT0VoAg/VtEFEREb2qp+kBoJiIlgCYCeBPzCyKJEVceWUiPY5A9+7d46YZOHBgQgPzBvfee69td1C/fv3MlkaLFi3QrFkzXHLJJQC0F8g638PLSxVvrMUwMHDDGh/GaF3JhEYh3aQjjpAdoURIZObdAC6y2V8M4E59/TsAZ6RZNMEjN9xwQ8ryVl2lnHPOOaYJ8tlnn43y8nJzwmF+fj569OiBH3/8EXl5ebj//vtRUVFhW4GrzixTzbBhw1BcXJyWsgQhE5BQu0JGM3DgwJjtgoKCGg4cBw8ejPz8fOTk5MRYfdWpU8ds2ZxyyinmZC31q61Ro0bo37+/rdmvX4x4MoDmG+zIkSOBlyEIKmG1giOjSIhoOIDhhm8lofbg1H3Wv39/MLOrh9/69eub40MHDhwAAHMGvl+YGY8++igA7QV/8MEH8cwzzwSSd0FBAcrLy5M+Pzc3F5WVlYHIIgiRsUsUp42ClZycHNSpU6dG4K28vDzUq1cPQ4cOtT0vKEUCxDqJNFpH6lejk3v9eHidb5Mo//Vf/5WSfGs76Zp3FFbog8goEkHwSqdOnTBixAgMGDDA8zlOSicZnIKGJdK95jSo2qZNm2REMnFy/V9bGD58eNgi+KJfv36hlCuKRKiV2PUlE5GjqW4iSiceHTt2tO2Oe+CBB3znrbptUY0WWrVq5er8UhD8IIpEEHQaNmzoaZLmz3/+c0/5DR061HYSpTFeAsQG+XKaAW/n+l+N3eIUXlW1UmvYsKGpPFQFY0fHjh3NOC5G15yhYLN1ToyXrkBrF6jgHVEkgqDgJYKiMRHSjXvuuQf9+/c3Fcj5558fUwkbSqNfv35xXfTbxUP56U9/6hg/BYCjU0ogNvqkei1DhgwBoDmpZGacfPLJuOaaa9C+fXuzpRMvKqZKMnOIgkS9TlVuJ7nihagGanqgNnB6btT5WUGPP6VqnCwZIqNIJNSuEAZOXV6qd2NAsyBz+pq38y7s5HFY9UDQpEmTmGPW7jpj2xg3MaJH5uXlxbRoWrZsGeNzzJqPtaVEROjRo0fMNdpVzk7OPcNGHUdSFYzaLeikGOyus1evXjH/i9pCVPcHbZp73nnn2a6HQWQUiVhtCWEQ5CA8cKLisZs8mZ+fb7YaVJo3b45Ro0aBiMzKSh0DMRSB4e6lT58+MZZpbdu2jVGIrVq1MitYZrZ1julm2WYM2Dt11amVq5qPF4eY8VoNTtEs1fuhKmC10q9fv76nsSPrdV133XUxrUMnt0BO98NLC9eO/Px80xosbC8KkVEkghA2N9xwg+myJVmMrhBrhaYG/7JDHcMwvA4wc0wFZ1Q2BQUFNcZV1BbJKaecYiodqxxqZWhXeV1zzTVxK2PjGnNzc2MUiRqGQHVNoyrriy464RBDdfRpYFeJq/s6dOgQc8zJAMFYHzp0aIwsahqnaIROH7NGWVZlZx1zs5sLZzdOFpaFlh2iSAQhILp3745zzjknsPwGDx5srhvBv1SYGZ06daqxzw01nPDJJ5+csD+y66+/3jVtjx494uZnYFc5GqgRKq3dh5deeimAE0qtoKAgJi9jLoVxLxo2bGhewx133GGmczMc6NixI5o3b44BAwbUGIRXFYF1nKJr165xwzM8/PDDrsdHjRqFpk2bmqGm69evb5Z58803m4rY2v2WTBTToBBFIggZipOJsqpU8vLyzArzl7/8pW33iZNyadeuXUy+btEmDVmCjArplpdbDBfDbc4ll1xijncQEUaPHg1AC73cqVMnNGrUyLSKs7s26zgWcKI1cPrpp5vHmRmNGzc2uxuNFo6dx+dRo0aZis4JLwYdKtddd5253q1bN0dP3fGMNlKJKBJBCAm3MYFBgwbFDIgDJypDp24uuxaGn75zIkL9+vXRokWLGhVx8+bNbT0jJxIozOCqq5KLtJ2bm4tf/OIX5rbROiMi1K1bF3l5eaZftrp169q26lSY2Rw/qVevXsy9a9y4MUaPHo0BAwaY+7t16+aan1vLxOlcdUJk2OMeiRAZRSJWW0K2oY4JWKlXrx4KCgpqWGZZcapsbrjhBtPCyM7qKt75Bl26dLEdWM/NzY3bheOV3r17+zrfrSUFwDRCsHoOuPTSS2OuId6gd4MGDeLOwVFxm5NkDT0AaPfBaG0kqoytjkzTTWQUiVhtCdmAU3eOU6V83333OeblFna4bdu2pgK59dZb0aZNGxAR+vbta3bZDB482BxUt2tJWLcfeeQRR1nUstXzjMH+QYMGYcSIEeZ+u4o0CKxWZnYRNJkZl19+OTp06ID8/Hycd9556NixY8zYh1OUzkRI9JyrrrrKtkvPyRJNJVmfbUERGUUiCNmA0wCz09erl66qeF+v6qBy9+7dY7pvDMXWrVs3FBYWOuZhdHN5QW1FGZZVF110Ec4880xzv1PXjmqJpN6reFZrVhL5os/JyYmxnLriiivinmNcYzq6n5wmQWYSokgEIY04OUVMpkJq3bp1DRf5yfrQKiwsRPv27ZOWyejbz8nJwdixY5OOU6+OEahf4vHGN1RUuRs1auQaRsAOp7EpFbeWYlBk0xhJ7Xb1KQgZiJMJrbVCa9SoERo1aoTZs2enQyxHHnjggZi5E0SUdvcdTpVu/fr1cdppp2Ht2rW2x61dYYmU59fx5dChQ3Hs2DFHubIJUSSCkGHEm6uRCF4qyXiVVrzjYY9LWuVr3LgxCgsLsWrVqpj9iRgZeE3jhyA9SodNZLq2xGpLqM0EXTFm2xexysknn5wy31Ppui833nhjSsI/p4rIKBKx2hJqA4kMttu5/AiqzDBwckmi4mTJ5nYvvNzTdN+D1q1bJzxxMUwio0gEIepkQishTBluueWWwPP0ez2ZomTDJhRFQkTXEdEKIqomopqe106kG0JEq4mohIiSMwMRhFpGEJVbtkVTTMYs2ivZcP1hE1aLZDmAEQAczU2IqA6AFwAMBdATwI1E1NMpvSBEHb9dHUF9PTu5Sc90grj+k046yfTHJQrmBKFYbTHzKiDuH9sfQAkzr9fTvgfgSgArUy6gIGQgDz/8cNx5KEH17TtVkjk5OXjooYeSyjNVJHPNXiy47PY3atTId6iAKJLJYyTtAGxRtkv1fTUgojFEVExExTt37kyLcIKQbpzcnqd6sN2Km2fesFGvs6CgwGzFeblHiUJEGX0v0knKWiRENB1Aa5tDv2fmT7xkYbPP9l9n5okAJgJAUVGRtDcFwSNevt7r1q0b2pyHRGa0W7n55psTcqvvlsbuPtWrV8+MGVLbSZkiYeaLfWZRCkANZ9YewDafeQpC5MjPzzcrzD59+tjG2WjTpo1rIKd4+ccLcesVL/66jMBNgDcfW2PHjrV1hum1teBn7CTZe+qX2267LZRyncjkme0LAHQjos4AtgK4AcBN4YokCJnHLbfcYlaaamTBxo0bmwPjhpdfrzRp0iQmrK0fVGeQXsZXnAI3OaFO3OvSpYtt5d68eXNPXnTt8NstOGrUKF/n26G6jc+EGfKhKBIiuhrABAAtAfyTiBYz86VE1BbAq8x8GTNXEtGvAXwBoA6A15l5RRjyCkIm42TN5RY4Kx4FBQUxYXmdsIubbuX2228314MaUzj33HNt9zt57lWjFl5yySWORgsGnTt3RosWLRKSSVVgqgNMuxjsVuIFyQLgOFtfjWkfFmFZbX0E4COb/dsAXKZsfw7g8zSKJghCAhjxRtKNGs8+UVSHkueff75tmgYNGpiBwbzy6KOPmut33nlnQuc6xWdRw/lefLHf0YLUERmTA/G1JQjZy1133RVKubm5uXFbJx06dIgZt3HLK2iuvvrqwPNMBZFRJOJrSxASp2HDhoGFzB0zZkzS57ZrZ2vZnxHk5eXFDXlc28nkwXZBEFJM586dA4v33bZt20DyiTJhWXmlmsi0SARBEDIddRwlSogiEQRBSBOpGEfJBESRCIIgCL4QRSIIgiD4IjKKRMx/BUEQwiEyikTMfwVBEMIhMopEEARBCAdRJIIgCIIvRJEIgiAIvhBFIgiCIPgiMrNjiGg4gOEADhDRWpekjQG4mXa1ALAryXNTfTxbZYt3vsgmsols4cvWyeU8d5i5Vi0AJsY5Xuzj3FQfz0rZ4p0vsolsIltmyJbsUhu7tj5L4bmpPh5m2X5ki3e+yJbc+SJbcueLbAFDupYSdIiomJnjh30LAZEtOUS25BDZkqM2ylYbWyTxmBi2AC6IbMkhsiWHyJYctU42aZEIgiAIvpAWiSAIguALUSSCIAiCL0SR6BDRECJaTUQlRPRYGsvdSETLiGgxERXr+5oR0ZdEtFb/barvJyIar8u4lIj6KvncpqdfS0S3JSnL60S0g4iWK/sCk4WI+unXWqKfSz5le4KItur3bjERXaYc+61ezmoiulTZb/s/E1FnIpqnyzyFiDzHRCWiDkQ0k4hWEdEKIrovU+6di2yh3zsiKiCi+US0RJftD275EVFdfbtEP16YrMw+ZJtMRBuU+9Zb35/W90E/vw4R/UBE00K/b6mwKc62BUAdAOsAdAGQD2AJgJ5pKnsjgBaWfeMAPKavPwbgz/r6ZQD+BYAA/ATAPH1/MwDr9d+m+nrTJGQ5H0BfAMtTIQuA+QAG6uf8C8BQn7I9AeBhm7Q99f+wLoDO+n9bx+1/BvA+gBv09ZcA3J2AbG0A9NXXGwJYo8sQ+r1zkS30e6dfSwN9PQ/APP1+2OYH4B4AL+nrNwCYkqzMPmSbDOBam/RpfR/08x8E8A6AaW7/Qzrum7RINPoDKGHm9cxcAeA9AFeGKM+VAN7Q198AcJWy/03WmAugCRG1AXApgC+ZeQ8z7wXwJYAhiRbKzLMB7EmFLPqxRsz8PWtP8ZtKXsnK5sSVAN5j5mPMvAFACbT/2PZ/1r8ELwTwgc11epFtOzMv0tcPAlgFoB0y4N65yOZE2u6dfv2H9M08fWGX/NT7+QGAi/TyE5LZp2xOpPV9IKL2AC4H8Kq+7fY/pPy+iSLRaAdgi7JdCveXLUgYwL+JaCERjdH3tWLm7YBWEQA4OY6cqZQ/KFna6etBy/hrvSvhddK7jpKQrTmAfcxc6Vc2vdugD7Qv2Iy6dxbZgAy4d3r3zGIAO6BVsutc8jNl0I/v18tPyXthlY2Zjfv2P/p9e5aI6lpl8yiD3//0OQCPAqjWt93+h5TfN1EkGnZ9k+myiz6XmfsCGArgV0R0vktaJznDkD9RWVIh44sATgHQG8B2AE+HKRsRNQAwFcD9zHzALWm65bORLSPuHTNXMXNvAO2hfQn3cMkvVNmI6HQAvwXQHcDZ0LqrfpNu2YhoGIAdzLxQ3e2SX8plE0WiUQqgg7LdHsC2dBTMzNv03x0APoL2MpXpTV/ovzviyJlK+YOSpVRfD0xGZi7TX/ZqAK9Au3fJyLYLWldErmW/Z4goD1pF/Xdm/lDfnRH3zk62TLp3ujz7AHwNbXzBKT9TBv14Y2jdnSl9LxTZhuhdhczMxwBMQvL3zc9/ei6AK4hoI7RupwuhtVDCu29uAyi1ZYHmBXk9tAEnY3CpVxrKrQ+gobL+HbSxjb8gdpB2nL5+OWIH9ObziQG9DdAG85rq682SlKkQsQPagckCYIGe1hhcvMynbG2U9Qeg9fcCQC/EDiKuhzaA6Pg/A/gHYgcq70lALoLWx/2cZX/o985FttDvHYCWAJro6ycB+AbAMKf8APwKsYPG7ycrsw/Z2ij39TkAfwrrfdDzuAAnBttDu29pqaizYYFmdbEGWh/t79NUZhf9T1oCYIVRLrT+y68ArNV/jQePALygy7gMQJGS18+hDZaVABidpDzvQuvmOA7tq+SOIGUBUARguX7O89A9K/iQ7S297KUAPkVs5fh7vZzVUKxhnP5n/b+Yr8v8DwB1E5DtPGhN/6UAFuvLZZlw71xkC/3eATgTwA+6DMsBPO6WH4ACfbtEP94lWZl9yDZDv2/LAbyNE5ZdaX0flDwuwAlFEtp9ExcpgiAIgi9kjEQQBEHwhSgSQRAEwReiSARBEARfiCIRBEEQfCGKRBAEQfCFKBIh6yGiVkT0DhGt113NfE9EV/vM8wkielhff5KILk4yn96keNYNEiIqJMUbskuam1JRviAYiCIRshrd+dzHAGYzcxdm7gdt0lV7m7S51n1eYObHmXl6kiL2hmaTHxaFAESRCClFFImQ7VwIoIKZXzJ2MPMmZp4AAER0OxH9g4g+g+YcswERfUVEi/RYEKZXUyL6vR6DYTqA05T9k4noWn29HxHN0ls+XyguUL4moj+TFsNiDREN0uNBPAngetJiV1yvCq7L9ryyPY2ILtDXDxHR07qcXxFRS6X8JUT0PbQZy8a5hUT0jZ5+ERGdox/6E4BBevkP6I4I/0JEC3THg7/Qz29DRLP1dMuJaJDvf0aoNYgiEbKdXgAWxUkzEMBtzHwhgHIAV7PmKPNnAJ4mDaMl0wfACGhO+WLQfVZNgBaPoh+A1wH8j5Ikl5n7A7gfwH+z5oL7cWjxH3oz85QErqs+gEW6nLMA/Le+fxKAscw80JJ+B4DBevrrAYzX9z8G4Bu9/GeheQTYz8xn69d4FxF1htZq+YI1J4VnQZsBLwieSKqpLwiZChG9AM0tSIVeWQJ6PAgjCYD/1b0sV0Nzj90KwCAAHzHzET2fT22yPw3A6QC+1HrUUAea2xYDw1njQmhdSn6oBmAonrcBfEhEjaH5f5ql738LmtdoQIuX8TxpEfuqAJzqkO8lAM40WljQHPh1g+b36XVdWX7MzKJIBM+IIhGynRUArjE2mPlXRNQCQLGS5rCyfjM0h3z9mPm47kG1wDg9TlkEYIVNa8DgmP5bBW/vViViewUKnBLihHtvJxkfAFAGrTWRA63lZQcBuJeZv6hxQFOulwN4i4j+wsxvuosvCBrStSVkOzMAFBDR3cq+ei7pG0OL5XCciH4GoJO+fzaAq4noJCJqCGC4zbmrAbQkooGA1tVFRL3iyHcQWohbOzYC6E1EOUTUASdckgPau2m0Gm4CMIc1d+b7ieg8ff/Nluvazppb+FugtZbsyv8CwN16ywNEdCoR1SeiTtDuyysAXoMW1lgQPCEtEiGrYWYmoqsAPEtEjwLYCa0F8huHU/4O4DMiKoY2DvCjns8iIpqi79sEzW24tawKvUtovN7NlAvNlfgKFxFnAniMtEh7f7SMk3wLza244U1WHes5DKAXES2EFtHOGKgfDa0L6gg0pWDwNwBTieg6vUyjFbYUQCURLYEWb/yv0LrdFukWbzuhhWS9AMAjRHQcwCEAt7pckyDEIN5/BSEDIaJDzNwgbDkEwQvStSUIgiD4QlokgiAIgi+kRSIIgiD4QhSJIAiC4AtRJIIgCIIvRJEIgiAIvhBFIgiCIPji/wMVZTdZ26Q38wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "losses = np.array(losses)\n",
    "plt.plot(np.log(losses), lw=0.5, color=\"black\", alpha=0.5)\n",
    "plt.minorticks_on()\n",
    "plt.xlabel(\"Gradient updates\")\n",
    "plt.ylabel(\"Logarithmic loss\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 685,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ConditionalMLPRatioEstimator(\n",
       "  (mlp): MultiLayerPerceptron(\n",
       "    (mapping): Sequential(\n",
       "      (0): Linear(in_features=9, out_features=64, bias=True)\n",
       "      (1): Sequential(\n",
       "        (0): ELU(alpha=1.0, inplace=True)\n",
       "        (1): Linear(in_features=64, out_features=64, bias=True)\n",
       "      )\n",
       "      (2): Sequential(\n",
       "        (0): ELU(alpha=1.0, inplace=True)\n",
       "        (1): Linear(in_features=64, out_features=64, bias=True)\n",
       "      )\n",
       "      (3): ELU(alpha=1.0, inplace=True)\n",
       "      (4): Linear(in_features=64, out_features=1, bias=True)\n",
       "    )\n",
       "  )\n",
       ")"
      ]
     },
     "execution_count": 685,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ratio_estimator = ratio_estimator.cpu()\n",
    "ratio_estimator.eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 686,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run ../src/models.py\n",
    "%run ../src/summary_statistics.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 687,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 100\n",
    "\n",
    "b = np.linspace(-0.5, 0.5, resolution)\n",
    "mu = np.linspace(0, 0.1, resolution)\n",
    "\n",
    "B, MU = np.meshgrid(b, mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 688,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[-0.5000,  0.0000],\n",
       "        [-0.4899,  0.0000],\n",
       "        [-0.4798,  0.0000],\n",
       "        ...,\n",
       "        [ 0.4798,  0.1000],\n",
       "        [ 0.4899,  0.1000],\n",
       "        [ 0.5000,  0.1000]])"
      ]
     },
     "execution_count": 688,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inputs = torch.FloatTensor(np.stack([B, MU], axis=-1)).view(-1, 2)\n",
    "inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 689,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[552.0000, 143.6310,   9.1077,  ...,   1.9085,   1.7146,   1.6246],\n",
       "        [552.0000, 143.6310,   9.1077,  ...,   1.9085,   1.7146,   1.6246],\n",
       "        [552.0000, 143.6310,   9.1077,  ...,   1.9085,   1.7146,   1.6246],\n",
       "        ...,\n",
       "        [552.0000, 143.6310,   9.1077,  ...,   1.9085,   1.7146,   1.6246],\n",
       "        [552.0000, 143.6310,   9.1077,  ...,   1.9085,   1.7146,   1.6246],\n",
       "        [552.0000, 143.6310,   9.1077,  ...,   1.9085,   1.7146,   1.6246]])"
      ]
     },
     "execution_count": 689,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "obs.repeat(1000, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 690,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([543.0000, 119.5034,   6.2711,   2.0812,   1.6490,   1.5175,   1.4556])"
      ]
     },
     "execution_count": 690,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = EvoModel(N=500, mu=0.1, burn_in=1000, T=10, b=-0.1).run()\n",
    "obs = np.array(hill_numbers(m.trait_fd, q_step=0.5))\n",
    "obs = torch.FloatTensor(obs)\n",
    "obs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 691,
   "metadata": {},
   "outputs": [],
   "source": [
    "outputs = obs.repeat(inputs.size(0))\n",
    "\n",
    "with torch.no_grad():\n",
    "    log_posterior = ratio_estimator.log_ratio(inputs, outputs) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 696,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7fb3d1461610>"
      ]
     },
     "execution_count": 696,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAWl0lEQVR4nO3dbYxc133f8e9/HynZlmzSD3FIuRQgBgXdGGnNSCnQxEUUO9SLigGiNJRbhAEEqEUqtKjhtjLaAokCFJX7oBawCoSoXCjqC0khEJSInaqm1ARBYKui7NQqpTKmlbaiZVRQSdOi+LQP/77YS86Zy5nZS3J2d3j4/QAL3Tv3zMyZC+o/vz333LORmUiS6jW10R2QJK0tC70kVc5CL0mVs9BLUuUs9JJUOQu9JFWuU6GPiN0RcTQijkXEwwOO/0xEfDMiFiPivtaxfRHxneZn37g6LknqJlabRx8R08CfAp8GjgMvAfdn5qtFm+3ALcDngYOZeaB5fDNwGNgFJPAy8MnMPDnuDyJJGqxLor8TOJaZr2fmBeBpYE/ZIDP/V2Z+G1huPffnga9l5ommuH8N2D2GfkuSOprp0GYr8Eaxfxy4q+PrD3ru1najiHgQeBBgatPsJ+e3bQFgdrr/e2Nuaqm3HYvF44v97cpjxXfPTPS/7/SI77mF4nnnstfuzPJ8X7uzS7O9dsU2wMJi8frFdiz1NevbL7enltvter99Rfsrdbk8VvyW1v6Nbbl44nLrWA55Xus1+va8s1qaCO9w8u3M/NCgY10KfQx4rOv/3Z2em5n7gf0AN+/4aP75f/sAAB957zt97bbefOrS9sc2nbi0ffv8W33tPjbbO/ajM6cvbX9kqv/jvndq09COv7XUe96rC++7tP3Ns9v72r3yzrZL20dPfrj/NU70nrd8ovcFMXey/wtm/ge97dniI8+/01/NZ98tvrRO939bTJ/rfblNnVnobZ+70NeOs+d62+fO9x3KCwvFdvG8pf5+5FLx3rlcbI74Z5HtbyZJ43QoD/zvYce6DN0cB24r9rcBb3Z872t5riRpDLok+peAHRFxO/A9YC/w2Y6v/xzwzyLiA83+Z4AvjHrC0uI0J95aScI/PN2fuL//nlsubX/3pi2Xtl/Z1D8a9KFNvVi8ZfbdS9ubZ97ta3fzVH+iLb292Evjf3b2g73tdzb39+lUr09nTt3Ud2zqB72hnPlTvV9u5n7Q14y5IsXPne4l3zLBA8y820vS02f7h6umzvbSeFwojp1f6GvHYpHGW0mdvqReDt1cZRo3xUsTYdVEn5mLwEOsFO3XgGcz80hEPBIR9wJExE9GxHHgl4DfiogjzXNPAL/JypfFS8AjzWOSpHWy6vTK9TZ/22257e/9fQCWbmpdBLypSLQ3FRdc5/vT7c2beuPLN8/1tm+a6U+35cXdtncWemPqJ9+9+dL2u6f7L8YuvzN3aXvmh/3fm7PvFCn+h+Xj/Z9rrtifLVL7zNnWOHyR4uNc/2cuU3wUY+2cb43RF2Pv5Zg8AIu918iFYnupdZ4cl5cmzqE88HJm7hp0zDtjJalyXcbo19XUIsyfXEnCS2f7J+0sbertL89PX9o+NzfX1+7cfG9s/8RML1VOz/YnzCjmL2a23utC7/U529ueOjPd127udO95s/2XAJg9XW73ku/c6f4UPDMkxU+f6zgOD1AkcBYWBm8DWaR2Wkl92Gyazkzw0kQy0UtS5Sz0klS5iRu6iaXehculc/3HluaLoZu5Yrv/hlSW53rfX8vFJ8zp/iGT5UG3czWmF3oHp4shpOlWn2bO9LZnz/Qfm323vMha3Ox0pn+Io7zI2nfj04gLrrQvpJb75RDMcnsKZbF/2V2zQy6stoZkRl6AlTRxTPSSVLnJS/TLvemGS3OtC6TF/U1lil9ufYosFrXpS/RTrQhf7rZCarl8zvT5wdsAM2dz4DbA7Jni2JnyImv/RdCpItFPne94wbWd6MsUP2SaJNCX8K962qSk64qJXpIqN3mJfglmm2UAptuJfnZwUl9uLUuZ5bFyNmR7TL5jop8q7jmaOd/fcPpckdrPtZYsKKdKFttT51tj72WKXyhS9mXLFwyfGtmX4stjo6ZQjuNmJ6dUShPPRC9JlZu8RL+cl2aoLF/oj+DTfWPv5ayb1lIJxVh8Fok+R8yyiVGJfqF3cPpCf8Op88VsmvbY+/nBKb499h5lwh924xP0J/rlEUsHj5pZM2KxsmHj8o7XS9c3E70kVW4iE/3FeeVTi/3fQ8vTZVIvtttj9GWiL15iZKJvDTVPlX/NqUj0U4utMfoytV/oT/TlwmNRpPE4355NUyT3xcGzZ4D+GTPt2TRlUi/Sfbbn0TsuL91wTPSSVDkLvSRVbuKGbsjeMgC50P89NDXT2+8bnhk1dDNiuKbUvhgbxRBNFBc3pxb6hy3K4Zq+qZG0LrqWQy3tIZlhwzWLo1aXbHV42HBN+69IdbwpyguwUj1M9JJUuYlL9LGclxbzypn+76G+hF9ejL1saYPBF2PLx1cOFqm9fX2xSLSx2DsYrRuQyhTfTvRDU/zi8JudygR+2RIFIxL90BQ/agqlF1WlG4KJXpIqN3GJnsze2PZSfwKPqXKMfqp4vPUaZXJvp/hhWtMQo5heWR4r0z3QGl8fldSHL0swNMVftuhYDm4H/dMmrzKpDx2XN/lL1zUTvSRVbiIT/cWx7VhspfHpIsXH4JunmoPd3+uiVmiNMuEvD74Z6bL9dqIvX2PYomOjXr89Dj9qQbIhqfuylG46l244JnpJqtyEJvpmSYB2Mi/XHC6OXZb8rybRXzYvfcjMlfa89L7ZLsMT/dBFx0a1G/Wn/kYk864za0bOlTf5S9Uw0UtS5Sz0klS5CRy6oXdRsz0EUw55TF3FFMpRLlsqYMiwzogpj+3XGDoM0/6rT+XwT98F4u4XUscyXCOpSiZ6SarcBCb67E1FbCf14iYplkctLj/kWPvi5qhjw1J2O7XnkGmYMHSq5GVrxA9L8R3/AtSgtlfFC7BSlUz0klS5yUv05OXj5RdNFY+PGpcfS6IvEvio1N6Xxocn9aHj8O3X7LiMcNf07Zi8JBO9JFVu8hJ9Dlgi4NKx4iapy1Yyu4q36prUS11Te/s1Oqbzq03xV5XcHZOXbggmekmq3AQm+uzNVplqfQ8VSwcnHcfr2689TDuND3neZal91Otf69z2tZgPb4qXbjidEn1E7I6IoxFxLCIeHnB8PiKeaY6/GBHbm8dnI+LJiHglIl6LiC+Mt/uSpNWsWugjYhp4HLgH2AncHxE7W80eAE5m5h3AY8CjzeO/BMxn5o8DnwT+1sUvAUnS+ugydHMncCwzXweIiKeBPcCrRZs9wK832weAL8XKgvEJvCciZoCbgAvAD1d9x4tDIMMuyq50pEPXO74PVzgkU7qKi6decJW0nroM3WwF3ij2jzePDWyTmYvAKWALK0X/XeD7wP8B/mVmnmi/QUQ8GBGHI+LwhTx3xR9CkjRcl0Q/KDq3o+WwNncCS8CPAh8A/igiDl387eBSw8z9wH6AW6e35KXpleNI7aOMSuqlcSfwtb7ZyRQvqdAl0R8Hbiv2twFvDmvTDNPcCpwAPgv858xcyMy3gD8Gdl1rpyVJ3XUp9C8BOyLi9oiYA/YCB1ttDgL7mu37gBcyM1kZrvnZWPEe4KeA/zny3ZKVBL2cK4l7LX+WO/7k8tCfXM6hP/2fq3jeqI8/7Pmr6fj6km48qxb6Zsz9IeA54DXg2cw8EhGPRMS9TbMngC0RcQz4HHBxCubjwHuB/8HKF8Z/yMxvj/kzSJJGiOw6Tr1Obp3akj81t3tlZwzLHIw07rHyq0zTzqaRdK0O5YGXM3Pg0LhLIEhS5SZvCYTSGqfW9UzSY1ku2BQv6SqY6CWpchZ6SarcxA3dJBv0V5HWYFjkmj+HQzWSxsBEL0mVm7hEv3IzU5NkxzG9csypeM1/2zDFSxozE70kVW7yEn1pHdPtul4XMLVLWkcmekmq3GQn+jHbkNk8vTffuPeWdEMz0UtS5SYy0W9o8r4WpnZJE8hEL0mVs9BLUuUmcuhm4jgkI+k6ZqKXpMrd2InepC7pBmCil6TK1ZnoTeqSdImJXpIqN5mJ3kQuSWNjopekylnoJalyFnpJqpyFXpIqZ6GXpMpZ6CWpchZ6SaqchV6SKmehl6TKWeglqXIWekmqnIVekirXqdBHxO6IOBoRxyLi4QHH5yPimeb4ixGxvTj2iYj4ekQciYhXImLT+LovSVrNqoU+IqaBx4F7gJ3A/RGxs9XsAeBkZt4BPAY82jx3BviPwN/OzI8DfxVYGFvvJUmr6pLo7wSOZebrmXkBeBrY02qzB3iy2T4A3B0RAXwG+HZm/neAzPx/mbk0nq5LkrroUui3Am8U+8ebxwa2ycxF4BSwBfgxICPiuYj4ZkT8w0FvEBEPRsThiDi8wPkr/QySpBG6/OGRGPBYdmwzA/wV4CeBM8DzEfFyZj7f1zBzP7Af4JbY3H5tSdI16JLojwO3FfvbgDeHtWnG5W8FTjSP/2Fmvp2ZZ4CvAn/pWjstSequS6F/CdgREbdHxBywFzjYanMQ2Nds3we8kJkJPAd8IiJubr4APgW8Op6uS5K6WHXoJjMXI+IhVor2NPDlzDwSEY8AhzPzIPAE8FREHGMlye9tnnsyIv41K18WCXw1M7+yRp9FkjRArATvyXFLbM674u6N7oYkXVcO5YGXM3PXoGPeGStJlbPQS1LlLPSSVDkLvSRVzkIvSZWz0EtS5Sz0klQ5C70kVc5CL0mVs9BLUuUs9JJUOQu9JFXOQi9JlbPQS1LlLPSSVDkLvSRVzkIvSZWz0EtS5Sz0klQ5C70kVc5CL0mVs9BLUuUs9JJUOQu9JFXOQi9JlbPQS1LlLPSSVDkLvSRVzkIvSZWz0EtS5Sz0klQ5C70kVc5CL0mVs9BLUuUs9JJUuU6FPiJ2R8TRiDgWEQ8POD4fEc80x1+MiO2t4x+LiNMR8fnxdFuS1NWqhT4ipoHHgXuAncD9EbGz1ewB4GRm3gE8BjzaOv4Y8PvX3l1J0pXqkujvBI5l5uuZeQF4GtjTarMHeLLZPgDcHREBEBG/ALwOHBlPlyVJV6JLod8KvFHsH28eG9gmMxeBU8CWiHgP8I+A3xj1BhHxYEQcjojDC5zv2ndJUgddCn0MeCw7tvkN4LHMPD3qDTJzf2buysxds8x36JIkqauZDm2OA7cV+9uAN4e0OR4RM8CtwAngLuC+iPgi8H5gOSLOZeaXrrnnkqROuhT6l4AdEXE78D1gL/DZVpuDwD7g68B9wAuZmcBPX2wQEb8OnLbIS9L6WrXQZ+ZiRDwEPAdMA1/OzCMR8QhwODMPAk8AT0XEMVaS/N617LQkqbtYCd6T45bYnHfF3RvdDUm6rhzKAy9n5q5Bx7wzVpIqZ6GXpMpZ6CWpchZ6SaqchV6SKmehl6TKWeglqXIWekmqnIVekipnoZekylnoJalyFnpJqpyFXpIqZ6GXpMpZ6CWpchZ6SaqchV6SKmehl6TKWeglqXIWekmqnIVekipnoZekylnoJalyFnpJqpyFXpIqZ6GXpMpZ6CWpchZ6SaqchV6SKmehl6TKWeglqXIWekmqnIVekipnoZekynUq9BGxOyKORsSxiHh4wPH5iHimOf5iRGxvHv90RLwcEa80//3Z8XZfkrSaVQt9REwDjwP3ADuB+yNiZ6vZA8DJzLwDeAx4tHn8beCvZeaPA/uAp8bVcUlSN10S/Z3Ascx8PTMvAE8De1pt9gBPNtsHgLsjIjLzW5n5ZvP4EWBTRMyPo+OSpG66FPqtwBvF/vHmsYFtMnMROAVsabX5ReBbmXm+/QYR8WBEHI6IwwtcdliSdA1mOrSJAY/llbSJiI+zMpzzmUFvkJn7gf0At8Tm9mtLkq5Bl0R/HLit2N8GvDmsTUTMALcCJ5r9bcDvAr+Smd+91g5Lkq5Ml0L/ErAjIm6PiDlgL3Cw1eYgKxdbAe4DXsjMjIj3A18BvpCZfzyuTkuSulu10Ddj7g8BzwGvAc9m5pGIeCQi7m2aPQFsiYhjwOeAi1MwHwLuAP5pRPxJ8/PhsX8KSdJQkTlZQ+K3xOa8K+7e6G5I0nXlUB54OTN3DTrmnbGSVDkLvSRVzkIvSZWz0EtS5Sz0klQ5C70kVc5CL0mVs9BLUuUs9JJUOQu9JFXOQi9JlbPQS1LlLPSSVDkLvSRVzkIvSZWz0EtS5Sz0klQ5C70kVc5CL0mVs9BLUuUs9JJUOQu9JFXOQi9JlbPQS1LlLPSSVDkLvSRVzkIvSZWz0EtS5Sz0klQ5C70kVc5CL0mVs9BLUuUs9JJUOQu9JFXOQi9JletU6CNid0QcjYhjEfHwgOPzEfFMc/zFiNheHPtC8/jRiPj58XVdktTFqoU+IqaBx4F7gJ3A/RGxs9XsAeBkZt4BPAY82jx3J7AX+DiwG/h3zetJktZJl0R/J3AsM1/PzAvA08CeVps9wJPN9gHg7oiI5vGnM/N8Zv4ZcKx5PUnSOpnp0GYr8Eaxfxy4a1ibzFyMiFPAlubxb7Seu7X9BhHxIPBgs3v6UB442qn3a+uDwNsb3YkJ4bno8Vz0eC56JuFc/LlhB7oU+hjwWHZs0+W5ZOZ+YH+HvqybiDicmbs2uh+TwHPR47no8Vz0TPq56DJ0cxy4rdjfBrw5rE1EzAC3Aic6PleStIa6FPqXgB0RcXtEzLFycfVgq81BYF+zfR/wQmZm8/jeZlbO7cAO4L+Np+uSpC5WHbppxtwfAp4DpoEvZ+aRiHgEOJyZB4EngKci4hgrSX5v89wjEfEs8CqwCPydzFxao88ybhM1lLTBPBc9nosez0XPRJ+LWAnekqRaeWesJFXOQi9JlbPQNyJic0R8LSK+0/z3AyPa3hIR34uIL61nH9dLl3MRET8REV+PiCMR8e2I+OWN6OtauJYlP2rT4Vx8LiJebf4NPB8RQ+dyX+9WOxdFu/siIiNiYqZbWuh7Hgaez8wdwPPN/jC/CfzhuvRqY3Q5F2eAX8nMi8tb/JuIeP869nFNXMuSH7XpeC6+BezKzE+wclf8F9e3l+uj47kgIt4H/F3gxfXt4WgW+p5yGYcngV8Y1CgiPgl8BPgv69SvjbDqucjMP83M7zTbbwJvAR9atx6unWtZ8qM2q56LzPyvmXmm2f0GK/fK1KjLvwtYCYFfBM6tZ+dWY6Hv+Uhmfh+g+e+H2w0iYgr4V8A/WOe+rbdVz0UpIu4E5oDvrkPf1tqgJT/ay3b0LfkBXFzyozZdzkXpAeD317RHG2fVcxERfxG4LTN/bz071kWXJRCqERGHgB8ZcOgfd3yJXwO+mplvXO8Bbgzn4uLrfBR4CtiXmcvj6NsGu5YlP2rT+XNGxN8EdgGfWtMebZyR56IJgY8Bv7peHboSN1Shz8yfG3YsIv5vRHw0M7/fFK+3BjT7y8BPR8SvAe8F5iLidGaOGs+fSGM4F0TELcBXgH+Smd8Y1OY6dCVLfhxvLflRm05LmETEz7ESED6VmefXqW/rbbVz8T7gLwB/0ITAHwEORsS9mXl43Xo5hEM3PeUyDvuA/9RukJl/IzM/lpnbgc8Dv309FvkOVj0XzXIYv8vKOfiddezbWruWJT9qs+q5aIYrfgu4NzMHBoJKjDwXmXkqMz+Ymdub+vANVs7Jhhd5sNCX/jnw6Yj4DvDpZp+I2BUR/35De7b+upyLvw78DPCrEfEnzc9PbEx3x6cZc7+45MdrwLMXl/yIiHubZk8AW5olPz7H6Bla162O5+JfsPLb7e80/wbaX4pV6HguJpZLIEhS5Uz0klQ5C70kVc5CL0mVs9BLUuUs9JJUOQu9JFXOQi9Jlfv/LVTk0oWbk2gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(B, MU, log_posterior.exp().view(100, 100))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}