SimulationBasedInference.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "678fdd9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "import scipy.stats as sps\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "175e0d6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_hypo1 = sps.multivariate_normal([0,0],np.array([[2,0.2],[0.2,1]]))\n",
    "X = dist_hypo1.rvs(size = (10,))\n",
    "\n",
    "def simulator(theta, N = 10):\n",
    "    dist_hypo1 = sps.multivariate_normal(theta,np.array([[2,0.2],[0.2,1]]))\n",
    "    X = torch.FloatTensor(dist_hypo1.rvs(size = (N,)))\n",
    "    return X\n",
    "\n",
    "def true_log_likelihood_ratio(X, theta):\n",
    "    dist_hypo2 = sps.multivariate_normal(theta,np.array([[2,0.2],[0.2,1]]))\n",
    "    return dist_hypo2.logpdf(X).sum(axis = -1) - dist_hypo1.logpdf(X).sum(axis = -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "39e55279",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATsAAAEvCAYAAAA6m2ZKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjQUlEQVR4nO3deXhU9dn/8fc9M9kJCZAAyiKoiLKIaEQQ911RrNXiVrtpbZ9qq62trdWni7VPNx9ttbSVWlutttZqcd8X3LUEQWWRfXNBQhJC9mRm7t8fE31of9Ql5yST5Hxe1zXXmJnJfe6Riw/fc873fI+5OyIifV0s2w2IiHQHhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkJLKx0bKyMh81alQ2Ni0ifdiCBQu2uHv5jt7LStiNGjWKysrKbGxaRPowM1v/n97TbqyIRILCTkQiQWEnIpGgsBORSFDYiUgkKOxEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkKOxEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkhBZ2ZhY3s4Vmdn9YNUVEwhLmyO4iYFmI9UREQhNK2JnZcGAGcGMY9UREwhbWyO6XwKVAOqR6IiKhChx2ZnYisNndF3zI5843s0ozq6yqqgq6WRGRjyWMkd10YKaZrQNuB44ws1v//UPuPsfdK9y9ory8PITNioh8dIHDzt0vc/fh7j4KOAN40t0/HbgzEZEQaZ6diERCIsxi7j4PmBdmTRGRMGhkJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgkKOxGJBIWdiESCwk5EIkFhJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgmhrnoiEmXuKUhvgtQW8PrMI5159nQjAGYJIA4W73jOh9ggiJW9/7BYYVa/R1+lsBP5mDxVBckl0P4GntoIqTc7Hu8AyQ/+3Y9S3wohPhISe2CJMZDYAxJjID4MMwvlO0SRwk7kA3i6Adrm4+2vdQTcUkhv/r8PxAZBfDjk7A35x2Px4RAbArFisOKO535gRWSOGqUyD092PDdDuhrSWzIjwvQWPL0FUmuhrRJvuff/tmX98NwKLHcq5E6FxJ6Y6UjUR6WwE9mOu0NyFbQ+hbc+Be2vkhmtxSCxG+ROw3LGQ874TNjEij/mFhKZh+V1/Nwf4kP+5RPbj908XZ/pJ7kCb18CbS/jrfM6PliC503H8k+AvCM6dpHlP9H/HRHAk2vw5rnQ8iCkNmZeTIyHovOw3GmQuw9mBd3el8WKIXcy5E5+PwQ9tQna/om3vQit8/CWBzOjycJZUDAL+7fwlAyFnUSWp7dBy4N48z+gfREQg9zpWNEXIe8wLD402y3ukMWHQsFMrGBm5qRI6zy86S94w/XQ8Bs87yis8CzInapjfNtR2EnkeNsivOkWaHkUaIPEGKz4Usg/qdeNiszikH8kln8knlyPN90OzXfhrY9AzmTo9w0s74Bst9kjKOwkEtyT0PIo3vSnzCjOiqFwFlZwCiQm9IkRkCV2wfp/Gy++CJrvxhtm47Xn4LkHYcVfx3ImZrvFrFLYSZ/m6UZovh1vvAXS70B8JFZ8BRR8Eov1y3Z7XcIsHwrPgIJPQNNteMMNePWpeN6xWPE3scQu2W4xKxR20ie5O7Q8gNf/NDNVJPcArPB7mWNxFs92e93CLB+KzoWC0/HGm6Dpj3j1c1Dycyz/6Gy31+00SUf6HG9fidd+Fq/7RuaKhIG3Exv4Zyz/yMgE3fYs1o9Y8dewsocgsSu+9QK84de4p7PdWrdS2Emf4ekG0tt+ilefDO3LsP4/wAbdheXum+3WegSLD8UG/gXyT8YbrsO3XvT+ZWxRoN1Y6RO85Ul82/cgXQUFp2HFl2Cxgd3bgzs1Tc1UNzbRkkzS0v7eo53WZIrC3BxKCvIpLcinpOORG+/ekaZZHpT8HHL2xOt/gdesg9LfYIkR3dpHNijspFfz9DZ824+hZS4kxmKls7HcSV26zbQ7KzdX89rb77Cueisba7eyvraOjbV1NLa1faxaJfl57DG4jD2HlrPnkMxjTPkg8hJd91fTzDLH8hJ74FsvxqtPhQG/xXL367Jt9gQKO+m1vP11vPaCzGiu6CtYv69gltsl29rS0Mizq9fxzKp1PL9mPXUtrQDkxOMML+3PyAGl7D9yGCMGlFDer4jC3BzyEgkKcnLITyTITcRpam+nrrmFuuYWtnY8v7OtnhWbt3DXwiU0tbcDEDdj6ugRzJo8kSPG7tZloz/LOxgG3YnXfhmv/QIM/CuWM65LttUTmPtHWYchXBUVFV5ZWdnt25W+w5vvxuuugFg5NuBXWM7eoW+jub2dh5au4I5XXmfhm+8AUN6vkIN3G8UBu4xg8oidGV7an3gs+KHvtDsbarbyxrtVLH7nXe5fvJx3ttUzqKiQUyeN51OTJzByYGng7eyIp6rw6k+CFWCD/tGrp+SY2QJ3r9jhewo76U3c03jD1dB4Y+ZyqNJfhn5s7o1NVdyx8HXuff0N6ltbGT1oADMn7sWhu49ir6GDiXXDBORUOs1zq9dz+yuvMW/lWtLuHLTrLlx8+IFM3Dn8y9i8rRKvOQfyj8VKru21k6wVdtInuLfhdZdBy31QeBZWfDlmOaHUTrvz6LKV3PTSAl59axO58TjH7jWG0/edSMXI7K4j9+62Bu5ctJjbKl+lprGJUyaN5xtHTKe8X1Go2/GG3+EN12D9r8QKzwi1dndR2Emv5+kGfOuF0PYC1u8SKDo/lAB6L+RmP/MSK6qqGTVwAGdV7M3Je4+jtCA/hM7D09Daym+efZlbXl5ITiLOVw4+gM9OmUxuSCcz3NN47Reh7WVs0N+xnL1CqdudFHbSq3m6Bq85F5JvYP1/jBV+MpS6T69cy9VPPseKzVsYPWgAFxw8lRPG7xHKMbiutK66lp8+9gxPrVzDyAElXDnjKKaNHhlKbU/X4FtmghVig+ZisXBHj11NYSe9lntz5lhS+3JswPVY3mGBaza0tvL9B5/k/sVvMGrgAC48pHeE3L97bvV6rnrkKd6u28bN55zG5OE7h1LXW1/Ga8/B+n0T63d+KDW7yweFXe/605VIcU/iW78O7Yux0mtDCboVm7dw6o1/5aEly/nqIVO578vncNLEPXtd0AEctNsu/OWzsxjav5gv334Pa6trQ6lreQdA7oF40824f7x5gz1Z7/sTlsjw+quh9Ums/39j+UcFrnfv68uYddNfaWht5U/nnMaFh07r9isYwjawqJDfn3kKZsYX/zqX6samUOpa0Rcy8xdbHg2lXk+gsJMeyZvvhaaboPAcrPDsQLXakkmufOhJvnX3w4zfaQhzv3g2U3YZHlKn2bfLwFJuOONkquob+dLtd9PU1h68aO5BEB+ON/8teK0eQldQSI/j7UvwusshZwpW/J1AtTZtq+drd97Pq29t4gtT9+MbR0wnp4tGc+7OpjdrqatpIJGIE0/EiMVjJBJxEjlxyoaWEI93zfhi0rCduOaTJ3Dh3+/j6/94gNmzZpIIsGtuFoOCWXjDNXhyNZbYLcRus0NhJz2KezO+9asQG4iV/irQPLrapmbOvvkOapqa+dWpMzhu3B4hdgqpVJp1Kzbx+vy1vD5/LYsXrGNb7X9eRWTQ4P4cduIkjjhpMqPHDg197t6RY3fje8cfzg8efJI/vFDJlw6aEqxgwWnQcB3efA9W/I1wmswihZ30KN4wB1JvYgNvxeKDOl0nlU5zydwHebe+kds++ykmDdsptB631jRw9y3P8+DtL1Nf1wzA0OEDOOCwPdlrn5EMGtKfdCpNMpkmlUyTSqZoaW5j/jPLufuW57nrpmcZNWYIh5+0D4efNJnyoSWh9XbmfpN4Yvka/jx/IZ+fum+gOXgWL8MTu2fuldsHKOykx/DkBmj8febGN7nBRiW/fOoFnl+zgatOPDq0oKvaVMddNz3Dw3+fT1trkgOPHs+BR45j4v6jKd+p9EN//4TTD6CutpFnH36dp+5bxB+veYRbf/0EX778JI7/1P6hjfQ+e8BkzvvLXB57YzUzJowNViwxBtoWhNJXtinspMfw+qvAcjJ3+grgkWUrmfPCfE7fdyKfmjwhcF/vbKzhjjnzePyeV0innSNO2odZXzyUEbsO/ti1SgYUceKZUznxzKm8vaGa2Vfew/Xfn8vr89fw1R+cQmFR3ocX+RDTd92Fof37ce/iZYHDzhJj8Jb78HRDr14gABR20kN4y5PQOg8r/nag2xmuqqrmO/c+wj7DduKKYw8L1pM79/3lRf7wi4dwh+NO25/Tzj2EIcMGBKr7np1HDuJHcz7H3+Y8za3XP8bKxW9x+S/PYvTYYCPRmBkzxo/l5pcXUtPYxMCiws4XS4zJPCdXQe4+gfrKNk09kaxzb8frfwyJ3aHwM52u09TWzgV33EdBTg7XnXZioONVLc1tfO9Lf+K3V93H3lN25aZHv8kF3zs5tKB7TywW48wvH85P/ngezY2tXHz6b3j+sSWB686csBfJdJqHlq4IVuj9sFsZuKdsU9hJ9rU+B6mNWL+vBzr7+veFr7OupparP3EcQ/oH2+X63Y/vY8FzK/mvK07iyhs+R9mQ8E4i7MjeU3Zl9tyvMWxUGX+69hGCXsY5dkgZpQX5LN+8JVhjsY6TROmtwer0AAo7yTpveQCsBPIO7XSN9lSKP738ChUjh3HgrsHui/rEPa/wyF2VnH7+Ycw8+8BuW96pdFA/TjprGm+urWL10rcD1TIz8hIJ2lMB7yDmHVNpevnxOggh7MxshJk9ZWZLzWyJmV0URmMSDe7N0PoE5B8TaEn1J5av5u26es6dFuw+CutXvcv1P7ybifuP5tMXHhmoVmccdMwEEjlxnn7otcC1cuIxkulUsCJen3m23h92YZygSAKXuPsrZlYMLDCzx9y9b0zOka7V+jR4I5Y/I1CZ2195jWEl/Tl099GdrtHS1Mb/XPwX8gty+fbVZxBPdP91s8WlhewxcTiLK9cGrpWIxYOP7NINmec+EHaBR3bu/o67v9Lx3/XAMmBY0LoSDd78AMTKIPeATtdYW13Li2s3MmvfiYFWL7l19uNsXFPFpT8/nUGD+3e6TlDjJu/CqqVv09oS7BrX3HiMZODd2I6w027svzKzUcBk4OUw60rf5J6E1nmZ+x5Y50dRDy1dgQGn7jO+0zVSyRT3/PkFjpi5D/tOH9PpOmEYPrqcZHuKtzdUB6rTlgq4CwuQ7ujBshf+YQkt7MysH3AXcLG7b9vB++ebWaWZVVZVVYW1WenNUu8ArVgi2O373tq6jbJ+RYHuyZBKpUm2pxi5W+fn+IWl8tnllAwsYuSu5Z2u8U5dPetqtjJ5RLA5e942P7ML2wcWAggl7CwzX+Au4DZ3/8eOPuPuc9y9wt0ryss7/4cofUhqQ+Y5HmxJ8c0NDQwuDrZ8eBYW7N6hms3beOnJZRw2Y1KgY4bPrl4HwCG7jQrWUNsLkHsAZr3/+oMwzsYa8AdgmbtfE7wliYz3wi4RbKrIu9saGNyv9x9TAph7y/OkU2lmfvrAQHWeWbWOnUuK2b2884speHJjZv5j7rRAvfQUYYzspgPnAEeY2aKOxwkh1JU+zpPrgVyIffxrTLe3uaEx8CTi9ybxZvN2qY31LTz4t5c56JgJ7Dyy8yHVlkrx/Nr1HLzbqGBzBNtezDznBQveniLw2NTdnwN65x11JbtSGyA+IrNQZCe1pVLUNjWHdg/VjWuzczzZ3fnr756kqaGV087r/ORqgBfXbKCprT3QNBwAb3s+8w9RvPcfrwNdQSFZlQQLtspH3IzCnBxqm5oD1ckvyGXqEXux6MVVbK1uCFTr40olU/z6h/dw103PcvQp+zFmfOdnbrW0J/mfR+cxrKQ/B+3W+cMDntoELY9D/tFZvUF4mBR2kj2xgZCuCVQiHoux97ChLNwY7PIqgDO+fDjbapu47PM3srWmewKvubGVKy/8Mw/+7WVOO/cQLr4q2D1xfzXvBdbVbOXHJx1NXoCFELzxJiCNFX4+UD89icJOsscGQLom8EXvk4fvzBvvVtHYFuy2f2MnjuAHv/0sb2+o5rtf+MMHLrEehqpNdVz6mTlUPruCC7//Cc795vHEAkyKXvjm2/zxpQWcvu/EQDfN9tQmaPoLFJyCJcK5+XZPoLCTrLHYQKANPNjt/yaP2ImUO6+//W7gnvaZuhvfn/0Z3ly7hcvO/QP1W8O5NeH2UskUc29+ji/NuIY3123h+7M/w4wzOn8FCWR2Xy+791F2Kinm0qMODlTLG28A0ljRfwWq09Mo7CR7YgMzzwF3ZffpWHY9jF1ZgH2nj+H7s89hw6rNXH7eTaEdw3N3XvvnGr42azZzfvoA4/cbxW/uvogph+0ZuPb1T7/I2uparjrxaPrldf44qKfehqY7oOA0LDEicF89Se+fKSi913srEqfWQ4C/WCUF+Rw4eiR/nr+IL07fP9AtBN+z30F78N/Xf5offfVWPnf0z5l+1HgmVIxmQsUoho8u/8gH7evrmlm97G2WLdrAk/cu5M21VQwa0p/Lf3U2048eH8rB/zsXLebGFyuZNXkC0wMsb+Xu+LYfAob1+3LgvnoahZ1kT+5+QD7e+gSWd1CgUmfutzdfvfN+/vTSK5x3YEUo7U05dE9+e89F3HHj08x/ejlP3rcIgJKBRYzfdxR7ThpBQVEeZobFjJgZZlBb3cCqpW+zeulbbHqz9v16E/YbxawvHsrBx04kv7Dzy1ltb+6rS7jivsc4aLdduOK4w4MVa7oNWp/Ciq/A4juH0l9PYkEPDndGRUWFV1ZWdvt2pedJ134F2l/Hyp8JNMpxd7769/uZt3INf/vCGYzfKdxrXN2dt9ZtYXHlOhYvWMuSBev+Jcj+3U4jB7L7uGEdj53ZfdzO9B8QzlzA9/r5/QuV/O+TzzFt9Ah+d/onyM8JcPa1fQVefWrm0rABv++1003MbIG77/BfO4WdZJU3z8Xrvo0NuhPL2TtQrdqmZk6ecyuFuTn847yzKczt/BLvH0V9XTPJ9iSedtLuuIOnnaLifIqK87tsu6l0mh8/Mo/bKl9lxvix/HTmMYHut+Hemgm6dDU26D4sXhZit93rg8JOJygku/IOB+J4y+OBSw0oLOBnJx/LuupafvLovMD1PkxxSQEDyooZOLg/ZUNKKB9awuCdS7s06Frak3ztzvu5rfJVzp22H1efcnygoAPw+l9AcgVW8rNeHXQfRmEnWWWxUsjdH1oexj34+mvTRo/k3GkV3LFwMbfOXxS4Xk+yeksNZ998B08sX80Vxx7GpUcdQizg7qY33QVNt0DhZ7G8Q0LqtGdS2EnWWcEZkFoHzTtcHexju+jwAzl8zGh+9PBT/M+j80ilA67Wm2UNrW1c+9TzfGLOrby5tY5fz5rJOVMmB67rTX/Ht30Xcqdjxd8KodOeTWdjJfvyj4emWzK7U3lHYPHOr/gBkBuPM3vWTH762DPc/PLC9+efDSnuXctAJdNp7lq0mOvmvciWxiZOHD+Wy445lLIQFj3wptvxbd+D3IOxAbMD3eyot9AJCukRvH0lXv0JyD+OWOn/hlb3rwte5SePPk1uPMGlRx3MaZMnBN7162ruzjOr1/GLx59lZVU1+47YmcuOPpS9hw0Np37jrXj9lZB3GFZ6PRZwMYaeRGdjpVdI118Hjb/GSm/A8gPOGdvOuupa/vuBx/nn+jeZsstwfjTjKEYNGhBa/bC4O8+uXs8fXqzkpXUbGTmghG8deTBH77l7aFNBvPFmvP7HkHckVvqrPjeiU9hJr+DelhndpRuxsgewEO9o5e7cuWgJP3vsGVqTSb5y8AGcVTGJkoKuO3P6UdU1t/DQ0hX8+Z8LWbWlhvJ+RZx3YAVnVUwiNx7O7RzdHRpvwBuugbyjsdJr+1zQgcJOehFvW4jXnAW5+3dMbg13F2tzfQNXPfwUj7yxikQsxrTRIzl2r905cuzuDCwsCHVbH+Strdt4ef1Gnli+mqdXraM9lWLPIeV8fuq+nDB+bGghB+DpGrzuOx13cpuBlfyczG1j+h6FnfQq3nwPXvetjl2t67rkL+aSd97lwSUreOSNlWysrSNuxv67DOeosbsxadhO7DG4LNAVCdtzdzY3NFK5/i1eWreBl9ZtZENtHQDl/Qo5YfxYZk7ci/FDB4d+5YK3voDXXQrprVjxt6Hw07326oiPQmEnvc77B9HzZ3aMRLpmlpS7s2xTFY8sW8nDy1ayriZzCVjcjF3LBjJu6GDGDR3MiAEllBTkUZKfT/+CfEry88nPSZBMp2lua6exrY3Gtnaa2trYXN/Imuoa1mypZc2WGtZU17CtpRWA4rw8puwynKmjRzB11AjGlA/qkvBxb8cbfgmNN0J818xua07w1VV6OoWd9Ere8Fu84VrIPx4r+UWXH2Nyd96q28bSdzazdFMVSzdtZummzVQ17HgRz5xYjPYPmMNX3q+IXcsGsuugAexaNpDJw3di3NDBxENYleWDeHJVZre1/TUoOB3r/13Mum8XPZs+KOw0z056LOv3X2C5eP3P8HQtlM7GYsVdtz0zhpeWMLy0hGP2GvP+61UNjWzaVs+2llbqmluoa2lhW3MrDa1t5OUkKMrNoSg3l8KO54FFBYweNID++d178sO9DRrn4A2/BSvKTCvJP7Zbe+jJFHbSo1nRuRArw+suw6s/BaVXYzkTurWH8n5Fod29rCu4O7Q+jtdfDam1mZMQxVcEnpzd1+hyMenxrOBkbMAfwBvw6ll4w2zck9luq0fwtvl4zZn41gsAw0pvIFZ6rYJuBxR20itY3jSs7AHIPx5v+BVefSre+lK228oKd8dbnyZdfSZeczakNmL9r8TK7g91MnZfo91Y6TUsVoKV/i/ecgy+7Sd47WfwvCOx4kuxRLAbQvcG7qnM6jCNcyC5DGI7YcVXQOGnInMCIgiFnfQ6ln8s5B0KjTfjjb/Dt8zAC8/C+l2AxXreZWBBeXorNN+DN92WWR0mPhrr/xMoOKlPXgXRVRR20iuZ5UO/L0HBaXjDddB0K958N154JlbwCSyxW7ZbDMQ9CW0v4c33QMvDQCvk7I0VX5e53MvCu8IiKhR20qtZfBBW8kO88NOZOXmNv8cbb8ATE7CCkzNnJnvJ6ruebsgEXOsz0Ppo5haT1g8KTsUKZ2E547LdYq+msJM+wXLGYAN+g6eqoOV+vPnezOoe9T/Fcw/A8g6G3GmQ2LPLrsb4uNzTkHwDWp/FW5+F9leAJFgR5B2C5c+AvEP71BJM2aSwkz7F4uVQ9Hms6POZKwma74HWJ/D6n3V8YACeNw3LnQ45kyAxutsuivfUu9D+Gt7+GrS/Cu2LwTtuwJ3YK9N37iGQO1nH4rqAwk76LEvsjhVfAsWXZIKm7QW89QVoexFvebDjUzl4YldI7IEl9oDE7hAfDLEyiA36WKHjngTfBql3IbUBkuvw1HpIbsicWEhv7vhkAhJ7Zq77zZ2UWRY9Pjjsry//RmEnkWDxIVBwClZwSuaKg9RqaF+GJ5dDcjm0VeIt9/1/v+dWArFBECsiMy01DljHM+D1kK4DrwPfwTW0sTKIj8wEWs64zGgyZy/tmmaBwk4ix8wyI7jE7hgnvf+6p7dBci2kt7z/8PQWSG0BbwbS//fwNOAQ3wkSYyFWisX6g5VAvAziu0B8ZKgLkEowCjuRDhbrD7mT/vW1LPUi4esZp6VERLqYwk5EIkFhJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgkKOxGJBIWdiESCwk5EIkFhJyKREErYmdlxZrbczFaZ2XfCqCkiEqbAYWeZxfBnA8cD44AzzUzrR4tIjxLGyG4KsMrd17h7G3A7cHIIdUVEQhNG2A0DNm7385sdr4mI9BjddoLCzM43s0ozq6yqququzYqIAOGE3VvAiO1+Ht7x2r9w9znuXuHuFeXl5SFsVkTkowsj7OYDY8xstGXuTnIGcG8IdUVEQhN4WXZ3T5rZhcAjZO5CcpO7LwncmYhIiEK5B4W7Pwg8+KEfFBHJEl1BISKRoLATkUhQ2IlIJCjsRCQSFHYiEgkKOxGJBIWdiESCwk5EIkFhJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgkKOxGJBIWdiESCwk5EIkFhJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgkKOxGJBIWdiESCwk5EIkFhJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgkKOxGJBIWdiESCwk5EIkFhJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgkKOxGJBIWdiESCwk5EIkFhJyKRoLATkUhQ2IlIJCjsRCQSFHYiEgmBws7MfmFmb5jZa2Y218xKQ+pLRCRUQUd2jwET3H1vYAVwWfCWRETCFyjs3P1Rd092/PgSMDx4SyIi4QvzmN0XgIdCrCciEprEh33AzB4Hhu7grcvd/Z6Oz1wOJIHbPqDO+cD5ACNHjuxUsyIinfWhYefuR33Q+2b2OeBE4Eh39w+oMweYA1BRUfEfPyci0hU+NOw+iJkdB1wKHOruTeG0JCISvqDH7H4NFAOPmdkiM/tdCD2JiIQu0MjO3XcPqxERka6kKyhEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkKOxEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkKOxEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkKOxEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkKOxEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpGgsBORSFDYiUgkKOxEJBIUdiISCQo7EYkEhZ2IRILCTkQiQWEnIpEQStiZ2SVm5mZWFkY9EZGwBQ47MxsBHANsCN6OiEjXCGNkdy1wKeAh1BIR6RKBws7MTgbecvdXQ+pHRKRLJD7sA2b2ODB0B29dDnyXzC7shzKz84HzAUaOHPkxWhQRCc7cO7f3aWYTgSeApo6XhgNvA1PcfdMH/W5FRYVXVlZ2arsiIv+JmS1w94odvfehI7v/xN1fBwZvt5F1QIW7b+lsTRGRrqJ5diISCZ0e2f07dx8VVi0RkbBpZCcikaCwE5FIUNiJSCQo7EQkEhR2IhIJCjsRiQSFnYhEgsJORCJBYScikaCwE5FIUNiJSCQo7EQkEhR2IhIJCjsRiQSFnYhEgsJORCJBYScikaCwE5FIUNiJSCQo7EQkEhR2IhIJCjsRiQRz9+7fqFkVsL4bN1kG9OWbd/fl79eXvxvo+4VtF3cv39EbWQm77mZmle5eke0+ukpf/n59+buBvl930m6siESCwk5EIiEqYTcn2w10sb78/frydwN9v24TiWN2IiJRGdmJSMRFLuzM7BIzczMry3YvYTGzX5jZG2b2mpnNNbPSbPcUBjM7zsyWm9kqM/tOtvsJk5mNMLOnzGypmS0xs4uy3VPYzCxuZgvN7P5s9wIRCzszGwEcA2zIdi8hewyY4O57AyuAy7LcT2BmFgdmA8cD44AzzWxcdrsKVRK4xN3HAVOBC/rY9wO4CFiW7SbeE6mwA64FLgX61IFKd3/U3ZMdP74EDM9mPyGZAqxy9zXu3gbcDpyc5Z5C4+7vuPsrHf9dTyYUhmW3q/CY2XBgBnBjtnt5T2TCzsxOBt5y91ez3UsX+wLwULabCMEwYON2P79JHwqD7ZnZKGAy8HKWWwnTL8kMLNJZ7uN9iWw3ECYzexwYuoO3Lge+S2YXtlf6oO/m7vd0fOZyMrtHt3Vnb9J5ZtYPuAu42N23ZbufMJjZicBmd19gZodluZ339amwc/ejdvS6mU0ERgOvmhlkdvNeMbMp7r6pG1vstP/03d5jZp8DTgSO9L4xn+gtYMR2Pw/veK3PMLMcMkF3m7v/I9v9hGg6MNPMTgDygf5mdqu7fzqbTUVynp2ZrQMq3L1PXIBtZscB1wCHuntVtvsJg5klyJxsOZJMyM0HznL3JVltLCSW+Vf3ZqDG3S/OcjtdpmNk9013PzHLrUTnmF0f92ugGHjMzBaZ2e+y3VBQHSdcLgQeIXPw/o6+EnQdpgPnAEd0/Jkt6hgJSReJ5MhORKJHIzsRiQSFnYhEgsJORCJBYScikaCwE5FIUNiJSCQo7EQkEhR2IhIJ/w8zlbCE6yvxKQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "grid = np.mgrid[-5:5:51j,-5:5:51j]\n",
    "thetas = np.swapaxes(grid,0,-1).reshape(-1,2)\n",
    "\n",
    "X = simulator([2,2], N = 10)\n",
    "nllr = np.array([-true_log_likelihood_ratio(X,t) for t in thetas])\n",
    "grid_z = nllr.reshape(51,51).T\n",
    "plt.contour(grid[0],grid[1], grid_z - grid_z.min(), levels = [0,1,4,9])\n",
    "plt.gcf().set_size_inches(5,5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "556e6815",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hn/y10jq2sn3cs3wnqv65l2rc040000gn/T/ipykernel_69490/577355997.py:16: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  theta_cols = torch.ones(X.shape[0],1)*torch.tensor(theta)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.7307738065719604 tensor([-2.2296,  2.6770])\n",
      "0.42677322030067444 tensor([-0.4010, -2.1199])\n",
      "0.31618884205818176 tensor([ 0.6400, -2.1653])\n",
      "0.11540721356868744 tensor([ 2.6744, -2.6860])\n",
      "0.33912453055381775 tensor([2.9047, 1.0208])\n",
      "0.5983898043632507 tensor([-1.3844,  0.4458])\n",
      "0.5828006863594055 tensor([ 0.2956, -0.9672])\n",
      "0.2861905097961426 tensor([-0.2395,  2.2018])\n",
      "0.2252758890390396 tensor([-1.2478, -2.7299])\n",
      "0.31601226329803467 tensor([ 1.4051, -1.9549])\n",
      "0.17281831800937653 tensor([-0.9609,  2.7089])\n",
      "0.6040951013565063 tensor([-0.5856, -0.9640])\n",
      "0.1972120702266693 tensor([1.6161, 2.6609])\n",
      "0.3866923153400421 tensor([-1.8585,  1.1358])\n",
      "0.3192443549633026 tensor([-0.4490, -2.2246])\n",
      "0.37083008885383606 tensor([1.5957, 1.6919])\n",
      "0.3060976564884186 tensor([-2.1864,  1.4814])\n",
      "0.27830708026885986 tensor([0.7025, 2.5109])\n",
      "0.16201743483543396 tensor([ 2.8468, -2.4704])\n",
      "0.1846662163734436 tensor([-2.1694, -2.5754])\n"
     ]
    }
   ],
   "source": [
    "def label_as(Xa,Xb, label_a, label_b):\n",
    "    ya = label_a * torch.ones((Xa.shape[0],1))\n",
    "    yb = label_b * torch.ones((Xb.shape[0],1))\n",
    "\n",
    "    X = torch.concat([Xa,Xb])\n",
    "    y = torch.concat([ya,yb])\n",
    "    return X,y\n",
    "\n",
    "def prep_data(theta, N = 123):\n",
    "    Xa = simulator([0,0], N = N)\n",
    "    Xb = simulator(theta, N = N)\n",
    "    X,y = label_as(Xa,Xb,0,1)\n",
    "    return X,y\n",
    "\n",
    "def parametrize_input(X, theta):\n",
    "    theta_cols = torch.ones(X.shape[0],1)*torch.tensor(theta)\n",
    "    return torch.cat([X,theta_cols],dim = -1)\n",
    "\n",
    "def train():\n",
    "    model = torch.nn.Sequential(\n",
    "        torch.nn.Linear(2 + 2,128),\n",
    "        torch.nn.ReLU(),\n",
    "        torch.nn.Linear(128,128),\n",
    "        torch.nn.ReLU(),\n",
    "        torch.nn.Linear(128,128),\n",
    "        torch.nn.ReLU(),\n",
    "        torch.nn.Linear(128,1),\n",
    "        torch.nn.Sigmoid()\n",
    "    )\n",
    "    opt = torch.optim.Adam(model.parameters(), lr = 1e-3)\n",
    "\n",
    "    for i in range(20000):\n",
    "        theta = torch.FloatTensor(np.random.uniform(-3,3, size = (2,)))\n",
    "        X,y = prep_data(theta)\n",
    "        X = parametrize_input(X, theta)\n",
    "        pred  = model(X)\n",
    "        loss = torch.nn.functional.binary_cross_entropy(pred,y)\n",
    "        loss.backward()\n",
    "        opt.step()\n",
    "        opt.zero_grad()\n",
    "        if i % 1000 == 0:\n",
    "            print(float(loss),theta)\n",
    "    return model\n",
    "\n",
    "trained_model = train()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "82b63403",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hn/y10jq2sn3cs3wnqv65l2rc040000gn/T/ipykernel_69490/577355997.py:16: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  theta_cols = torch.ones(X.shape[0],1)*torch.tensor(theta)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "tensor(107.3941)"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def eval_data(X, theta):\n",
    "    theta = torch.FloatTensor(theta)\n",
    "    X = parametrize_input(X,theta)\n",
    "    with torch.no_grad():\n",
    "        out = trained_model(X).detach()\n",
    "    a = torch.log(out/(1-out))\n",
    "    return -a.sum()\n",
    "\n",
    "eval_data(X,[2,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "9462c417",
   "metadata": {},
   "outputs": [],
   "source": [
    "def analyze_data(X):\n",
    "    grid = np.mgrid[-4:5:51j,-5:5:51j]\n",
    "    thetas = np.swapaxes(grid,0,-1).reshape(-1,2)\n",
    "\n",
    "    evals = np.array([eval_data(X,t) for t in thetas])\n",
    "    grid_ml = evals.reshape(51,51).T\n",
    "\n",
    "    nllr = np.array([-true_log_likelihood_ratio(X,t) for t in thetas])\n",
    "    grid_truth = nllr.reshape(51,51).T\n",
    "    \n",
    "    ml_result = grid_ml - grid_ml.min()\n",
    "    truth_result = grid_truth - grid_truth.min()\n",
    "    \n",
    "    return ml_result, truth_result\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "f7048f94",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/hn/y10jq2sn3cs3wnqv65l2rc040000gn/T/ipykernel_69490/577355997.py:16: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).\n",
      "  theta_cols = torch.ones(X.shape[0],1)*torch.tensor(theta)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAADQCAYAAAAJQ4B0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABVC0lEQVR4nO3dd3hUZfrG8e+Z9N4rIQkQeocA0qQooIiADRXsvaziqmvdtXfXVVfXgooNrKhYQTpIh9BrCCRAQnrvmXJ+f0T92QiTMJAA9+e6vC4zc+Y9zwnczDxzznlfwzRNRERERERE5MgszV2AiIiIiIjIiUINlIiIiIiIiJPUQImIiIiIiDhJDZSIiIiIiIiT1ECJiIiIiIg4yb05dhoeHm4mJiY2x65FmiQlJaXANM2I5q7jcJQpOdG09EyBciUnFmVKxPUOl6tmaaASExNZv359c+xapEkMw9jf3DU0RJmSE01LzxQoV3JiUaZEXO9wudIlfCIiIiIiIk5SAyUiIiIiIuIkNVAiIiIiIiJOUgMlIiIiIiLiJDVQIiIiIiIiTlIDJSIiIiIi4iQ1UCIiIiIiIk5SAyUiIiIiIuIkNVAiIiIiIiJOcm/uAk4ENruDGqudOpudOqsDm8OBYRi4GQYWi4G3hxu+Xu64u6kfFRERERE5mZ3yDZTdYZJTUsXBggqyi6vIK6smr6SagvIaSqvqKKuqo7LW5tRYXu4WAnw8CQ3wIjzAm7AAb2JD/Wgd5kfrMH8ignxwsxjH+IhERERERORYOaUaqDqbnT3ZpaRll5KaXcrenDIyCyux2h2/buPlbiEyyIfwQB9iQnwJ9PEk0McDXy93PD3c8HS34G6xYJomDhPsjvqzU5U1NiprrZRVWykqryGrqJLNGYW/a768PdxIigmiQ0wQHWKD6B4fRnigd3P8KkREREREpAlO6gbKNE3ScsrYsK+AjekFbD9YRJ2tvlkK9fciKTqQ5HYRtA73Jz7cn9hQPwJ9PDAM15wlMk2T0qo6MgsrOVhYQXpuOanZJXyXsp+6NfV1tIkMoH/7SAZ2iKJTq2CX7VtERERERFzPZQ2UYRhuwHogyzTNca4at7HqbHY2ZxSyKjWX1am5FJbXAvWNyri+CXRPCKVDTDBhAV6NblYcdjs1xcXUlJZiq6nBXleHvbYWh82Gw27HtNsBsLi7Y/HwwM3DA8+AABKDg+nUJQr33vFA/T1V6XnlbEwvYF1aHp+v3MenK/YSHezD8K6xjOjWisTIANf+YuSE01IyJXIyUa5EXEuZklORK89ATQV2AoEuHNMpdofJ5oxClmzPYsWuHCpqbHh7uNG3XQQDO0TRt104of4NXypXVVBA0d69lO7fT0lGBqUHDlCRk0Nlbi4VublU5edTU1JyVHV6BQYS1qEDYR07EtaxI927dWPUkJ64R/dh9Z58Fm/L4rOVe/lkxV66xIVwbnICQzpH4+nudlT7lRNWs2VK5CSmXIm4ljIlpxyXNFCGYcQB5wBPAne6Ykxn5JVW813KfuZvzqSoohYfTzcGdYxmWNcYercJ/8vGo6qggNwtW+r/27qVwl27KNi1i+qiot9t5xUUREBMDP7R0cT06YNvRAQ+oaH4hoXhHRyMu7c3bl5euHt5YXF3x3Bzw7DUz8LnsNlw2GzY6+qoKy+npqSEmpISyg8dojA1lQM//cTWmTN/3ZdnQAAxffow5owzuGz4KHa4hfP9xoM8O3sTb8zz5NzkBC44rS2+Xif1FZfyG82VKZGTmXIl4lrKlJyqXPWJ/CXgHuCw150ZhnEDcANAfHz8Ue1s96ESZq3ax/KdOYBJ/6RIzuwRR//2kXh5/H/TZLdaObRuHZmrV//6X9nBg78+7xsRQUSXLnS+8ELCO3UitF07gtu0ITghAa/AY/tFirWqirxt28jZvJnczZvJXLWKJQ8/DOZD+IaHM/m66/A5/yrm7ytnxrI9fJ9ygCuGd2BMrzjcLJou/RTwEscxUyKniJdQrkRc6SWUKTkFHXUDZRjGOCDPNM0UwzCGH2470zSnAdMAkpOTzcbuxzTrL9P7eEUam9IL8fd254LT2nBucgJRwb6/bpO7dSv75s8nfeFCMpYuxVpZCUBwYiKtBw0iNjmZqJ49ierRA/+oqMYfsIt4+PrSqn9/WvXv/+tjVQUF7J0/n52zZrHquecwXniBPpddxtnXTeWz1Epe/n4rs9emc9vY7nSPD2222uXYOl6ZEjmVKFcirqVMyanMFWegBgPjDcMYC3gDgYZhzDBN8zIXjI1pmmxML+SDJbvZmVVCqL8X15/ZmbF94vH1cqe2vJztn31G2pw57J03j/JDhwAI69iRXlddRZuRI2k9eHCzNkvO8g0Pp/ull9L90ksp2ruX1S++yMbp09ny4YeMvv12zp18C++u2M/d769ibJ94rj2jE/7eHs1dtrjeMc2UyClKuRJxLWVKTlmGabruy4Cfv4G4+0izsCQnJ5vr168/4ngHCip4Y94OUvbmExHozcWDkxjTKw5PdzeyN25k/euvs3XmTKxVVXiHhNBu1CjajRlDu9GjCYyLc81BNbPKvDwW/fOfbHj7bQJiYjjngxksdUTx1Zp0QgO8+deFfejUKqS5yzzpGYaRYppmcjPsdzguzJRIS9Fcmfp538NRruQko0yJuN7hctUiZyWosdqZuWwPX6zeh7eHGzeO6sy45AQ8LAap33/Pyuef58BPP+Hu40P3yZPpecUVtB48GIvbyTdbnV9kJOdOm0af66/nq8su49PRZzLmpZcYds3lPPXFBu5+fzX/mNCTYV1jm7tUEREREZGTnksbKNM0lwBLjmaMlH35vPLDNrKLqxjdM45rz+hEgIfBlg8/YOXzz1OwcyeBrVsz+oUX6HX11fiEHP3ZF4dpkltSTWZhBTklVeSUVJNfWk1ZtZWKmvr/rHYH/HyyzmIx8PNyx9/bAz9vDyICvYkJ8SUmxJe4UD9ahfnjZnHtgrit+vXjhpQUvrzsMubefjunpafz0qNP8vgXG3nqy43klFQxaVA7LcR7knFFpkTk95QrEddSpuRU02LOQFXWWPnf3O0s3JpFXKgfz11+Gl2j/UiZ9iYrn3+essxMonr25LwZM+g6aRJuHk2/96e6zsam9EI2ZRSQllPGvtwyqmptvz7v4WYhMsiHQB8Pgv08iQvzw9O9fuY7AwO7w6Sytr6xyi2pYtuBQipq/v/1Pp5utI8JomNsMAPaR9I1PhSLCxobT39/Jn3xBT/eeSerX3yRkowMnvjgQ16el8r0Rbs5VFzF1HO6u2RfIiIiIiLyZy2igUrPLePxWRvILq5i8pAkLh7Uhp0zZ/DKww9TlplJwumnc+7bb9Nu9Ogmn2Epqaxl2Y5sVqXmsnV/EVa7Ay8PN9pGBXBG91a0jQokPtyfmBBfQvy9Gt2ElFXXkV1cxYH8CnYfKmH3oRJmr83g81X7CA/w5vSuMZzZvRXtooOaVP8vLG5unP3yy4S0bcuPd9yBOWUy98yaRUyILx8vTyPQx5Nrz+h0VPsQEREREZG/1uwN1MItmbz8/Vb8vD147vIBeG1bzbvJF5K3bRut+vdnwnvv0faMM5o0dq3VzvKd2SzadogN+wpwmCatw/wY3y+B/kn1Z4Y83A6/plJtWRlFe/dSvHcvFbm5VBcVUVNcTG1ZGabdjmmamA4Hnv7++ISF4RseTmxMDL179CB09EBq7SardueydPshvlmbwZer0+mREMr5A9oyoEPkUZ0pOm3qVCzu7sz529/49rrruGL6dMqq6/hs5V7iw/0Z1fPkmERDRERERKQlabYGyu5wMG3+TmavzaB7fCg39Qxk5bUXk75wIaFJSVz0+ed0vuCCJp1xyiqq5PuU/fy4KZOKGitRQT5cNKgtI7rG0ibqrxfINU2T/O3bSV+0iIzFizm4ciWVeXl/2s7T3x+vwEAMNzcMiwXDMKgtL6e6qAh+M6Ohu48PUd2702rAAC4aPJgbLhnA6nwHs9dm8Mhn64kN9WXK0PaM7N6qyY1U/1tvpbqwkCUPP4x3aCg3P/9vsgrr14uKDfWla2utFSUiIiIi4krN0kCZJjz6WQpr9uQxsX8iPXYu4JMB9+Dh68vZr7xC3xtvbNI9ToeKKpm+aBc/7czBzWIwqGM045Lj6ZkQdthGrHjfPja++y6b33uPssxMAILbtCHp7LMJ79yZ0KQkQtu1IyA2Fu/gYNw8Pf/6mBwOakpKKD1wgJzNm8ndvJnsDRvY+M47rH3lFQAiu3Xj7gcepHjEED5fnc7zX2/mm3X7uWdiT+LC/Bt9vACn/+tfVBUWsuallwiIjeXBv01l6vQVPPZ5CtNuGkaQ71/XKyIiIiIijdcsDVRWcSVr9+Rx49BEav59Lwu+/56O48cz7s038Y+ObvR4FTVWPl6extdrM7BYDCYPTWJc3wTCArz/cnt7XR27Zs8m5c03SV+0CAyDpDFjGP7YY7QZOZLghIRG12BYLPiEhuITGkp0r17/vy+rldzNmzmwfDkb3n6b2ZMvJbxzZ27810PkjR/M6/N2cstby7llTBfO6h3f+P0aBme9+CKVOTksuPdeonv14qELB3Dr28t5d9Eu7hjXo9FjioiIiIjIX2uWBqq61sZ1PYM4cOVYKnJyOPuVV+h3662NvlzPNE0WbMnirQU7KauqY3SvOK4c3vGwjVPpgQOkTJvGhrffpjI3l6CEBIY/9hi9rrqKoNatXXFof+Lm4UFscjKxyckMuP12dsyaxdJHH+WryZfS9swzeeH1t3l9dTYvfreV9LxybhjVpdFToBsWC+OnTydv+3a+nDKFGzduZEK/RGavTWd8v0TaHuayRRERERERaZxmaaCiPWzsmnQm/lFRXLNiBbHJjV84O6e4ihe/38Km9EI6xwXz1OT+JMX8eYY70+Eg9bvvWP/666T9+CMAHc45h+Sbb6bdmDHHdfFdw2Kh66RJdLnwQlKmTePHO+/k84HJXP3++yyJbMNXa9LJKqrkwQv64OPZuD8aTz8/Js2axbTkZGZdfDEXz5nH/C2ZTJu/k6en9Nf6UCIiIiIiLnD4KeiOofL9GSQMHcoNKSlNap7W7snjlrd+IvVQKbeN7cZ/rhr0p+bJNE12fvUVb/TqxScTJpC3bRun/+tfTE1P59Jvv6X92LHHtXn6LcNiIfmmm7hh/XoCYmP57NxzOS17HbeP7UbK3gKe+XIjdod55IH+ILxTJ8a//TYHV6xg03+eZ/KQJDamF7Anu/QYHIWIiIiIyKmnWRqooLg4Jv/wA77h4Y16nWmazF6bzsOfriM62JfXrx/KuL4Jv5vFzjRN9syZw1v9+vHZ+edjq6nhvA8/ZGp6OiMefbRJ9zcdKxFdunDNihW0GTmSb665ho75O7l5TBdW78nj0xVpTRqz2yWX0OOyy/jpiSfoQSEebhbmb8l0ceUiIiIiIqemZrmEzy8qqtGz7NkdDl7/cQffrt/PwA5R3HdeL7z/cJlbwe7d/HDLLaQvWkRwYiIT3n2XHpddhsXducMsrarjYEEFBwsryC6uoqbOTo3VRq3VgZeHhWBfL4L9PAkN8KZjbDDRwT5HfWmcp78/F3/1FdOHDOHzSZO4ZsUKtneN5cOle+jVJpwucSGNHvOs//6X9EWLWHDT9Qz85zss3naI68/sjKd785xxk1OX3WGSX1ZNVlElh4qqKKqooaSyjtLKWsprrFhtDqx2B3U2B+5uFrzcLXi4W/D18iDU34tQfy/CArxpFepHfLg/wX6euhxVREREmlWzL6TrDLvDwWOfb2B1ai4XDmzLNSM7/W6iBdM0Wffaa8y7807cfXw4+9VX6Xv99Yedcvy3MvLKWbAlk6U7sskrrf71cTeLgY+nO94ebnh6WKizOiiprMX2m0vrQvy86N8+gnOTE2n/F/dfOcvT359Lv/2WtwcM4ONzz+WGtSnszCrm2a828vYtwxtc7Pev+ISEcPYrr/DZBRfQbv8GllXHsj4tn0GdGj/DoUhjFFXUsP1AMdszi9l+sIj03HKsdsevzxuAn5uJr6MWj7pqLHYrFlsdFpsVm+FGnZs7djd36ty9qHb3ocr4fYb9vdzp2CqYvu0i6Ns2goQIfzVUIiIiclydEA3UtPk7WZ2ay81jujCxf5vfPeew2Zhz++2sf/11Opx7Lue+9Rb+UVENjmd3mCzdfogvVu8jLacMN4tBcrsIJvZPJD7cn9Zh/kQE+fzapNmtVgpTU8nbtp2DO3aRmV1Evm842VWRLC6v5MdNmQxJCOCG8clEBfs26RiDWrdm0hdfMH3wYFY/+hB/u+Nh/vXJOpbvzGZEt1aNHq/TeecRP3Qo+f97Bo/rXmVHZrEaKHE5h2myO6uEtXvyWJuWR1pOGQCe7hbivR30tubgfiAV69b1OHZvxrO8GMOsb6g8/Pzw9PfH088Pdx8fME3sVisOqxVrVRXVxcXY7A6sAaFUh7eiOrw1VZHx7Gzbg5R9rYCdBHvAsB5xjO6dSLvoQDVTIiIicsy1+Abq+5T9zF6bwfkD2vypeaotL2fWpEmkzZ3L4Hvv5YynnsKwHP5sjd3hYMm2Q3y0PI3MwkoSIvy5ZUwXhnWNJdjP63fblh44wJ45c0ibM4d9CxZgrawE6ieA8I+Jwds0ibdaibRDep+xrBw4npUvzqNPdQb33j2FwMiIRh9r64EDGXD77az573+5avJkWoX68c26/U1qoAzDYMDtt/P5RRcR7Wlnb25Zo8cQ+SumaZKWU8bibVks3ZFNQVkNFgMSfRwMLNuNx5oF1K1ciMVuA8MgvFMnYvr0IfrisYR16EBwYiLBCQl4BTY8vb5pmthqaqguKqI8K4vSgwcpO3iQvG1r2bNqN/vsvhS178c3Ncl8nZJJhFsdZ/dJYOKI7vh5NX4hbhERERFntOgGalNGAf+bu51+SRFcd2bn3z1XlpXFR+ecQ962bZz71lv0ue66BsfalVXCS99tIT2vnDaRAfzzgj4M7hz9pwkodn/zDWt/vocIICghgR6XX0784MFEdutGWMeOePj4/O41RXv2sHHBMr5Kq2J9SHuuf2Qm953Tmd7njGn0MY984gl2fvklP9x0E+e8+RXTFu5mT3Zpky4R7DhhAgGtWuGVvoO9Rm9M09Q39NJklbVW5m3K5IcNBzhQUIG7xaAt5bTbvBBz7me4V5bi6e9P/JAhxD38EPGDB9Oqf388/f2btD/DMPDw8cGjVSsCW7WiVf/+v3u+rqKCg6tWsXXOPJbvziMjthsf2D35eGU6/fyqufqiYcS3PTbru4mIiMipq8U2UFW1Np79ahOtQv24/7zev7vnyVZTw8fjxlG8bx9TfviBdqNHNzjWqt25PPnFBoL8PHnwgj4M+UPjBJC/Ywff33IL+5cuJSg+nhFPPEHn888nvFOnBpsOwzAI69CBMzt04Ezgq+9X8tbaBJ6ZtZ7/xscS1b17o47b09+fs195hU8nTqR36io83MJZtiO7SQ2Um4cHfa6/nl3zVlMa3ZWKGhsBPvpmXhqnrLqOL1bt49v1+6mstdEmwMKQzJVYP3oNt4oSonv3pv0df6Pd6NHEDRzY6AlimsrT3592o0bRbtQoJlI/iczCWT8wL6ealbHdWP3uGoZVf8TNd0whKC7uuNQkIiIiJ78W20B9uiKNoopaHp7UFz/v338gm/v3v5OzaROXfvfdEZunxduyeG72ZpJiAnny0v4E+v55YomtH33EN9ddh7u3N+PefJPe117b5DWizjtnEFafrbyzwpPXH3yeh76c7vQsgL/oOH48wW3akDrjfaLP/ydZhRVNqgWg3ejRWH5YB9RfwijirBqrna/XpvPpir1U1droGeQgfME7VM3/Gs+QEPrfcA19b7iB8I4dm7tUAMI7duTiBztyMbBt7WZe+noDi4O6sfFf73NBcClj77vriPdHioiIiBxJs6wDdST5ZdV8uSadkd1i6dTq99N4b/3oI1LeeINB99xDh3POaXCcn3Zm8+xXm+jaOoRnLhvwp+bJbrUy9447+HLKFGKTk7l1xw763nDDUS+we8HwrsR42klpP4KVr7/Z6NcbhkHPK64gfdEiwn0s5JRUH/lFhxHTpw9uPx9PUxbnlVOPwzT5cdNBrvnfYqYv2k28WzXD5/wH36nj8dyzmbH/+x93ZmUx5oUXWkzz9Efd+vfkrSeu4pr+0VQkduM93748PPFqtn78MaapHIiIiEjTtcgzUDOX7cE04coRv/9wVpyezrc33ED8kCGMfOKJBsfYm1PKs19tolOrYJ6Y3B9vj983RQ6bjZlnn036woUMmDqVUc8/77JLj9wsFu68eDD/+HA1n8yZxcCbrI0eu8dll7H00UdxO7SfXM/IJtfi7uVFSEI8UP/BWKQhWYWVPP/1JnZmlZDoBx3nvQrL5xLYtStDZsyg28UXN+qMqtXuYH9eOanZpezNKSWvtJqiilpKKmqoqrHWNzOmA0wTH+z4uzkIcnMQ7edOt4QwenVNJCohrsHJYQ7HMAwuHtOXof068dTHq9g29lZeeusjzp49m4nvv4+7t3ejxxQRERFpcQ1Unc3Oku2HGNk9lug/TAm+/JlncNhsnP/RR0dsSD5dsRdPdwuPXtLvT80TwJaZM0lfuJCxr71Gv5tvdukxAPRIDCPAYqfUM5CK7GyC4uMb9frQpCSCEhI4UFQI0Ud32ZEREQOgmcnksEzTZP6WTP43ZzvuFjj9wDKsbz9HWMeOjPz8czqff75TTYxpmhwsrGRdWh5r9+Sx/WDxr+tAeWEnoLoES1EuRl42gbWVGD/39KbFgs3Hn5KAUPICw9kaEsX8onJI2YR/7mwSS/fTP9KdLv16kjBsGIGtnJ+ZMjbUjxdvGsnL321hPpOZ9+M7VI4dy+Tvv//dhDAiLVFljZWtB4pIyykjLbuUjPxyoH6tQg83C5FBPvRpG07fthHEhflpoiARkeOgxTVQG9MLqK6zM7RzzO8er8jJYfP779PzyisJat3wzFqF5TUs35XDhP6JBP3FPU+mw8HK554jqkcPkm+6yaX1/45hwXDYqS4ubnQDBeAfFUWp3Y2IwKP7przY3R/vqlJ8vVrcH7e0ALVWO6/O2ca8zZm0D4CYV/+B40AaZzz5JIPuvtupBanzSqtZsCWThVuyyCyqn/I/yt1Gx+xdmOuW4JW6Ca/iHIITEojp3Zvo3r2J6DII/+hofCMi8A0PxzCM+nWgbDaK8wrZvPsgOzIr2B0RzraYdmwDAhdsJeqZKfT0r6PLeefR+fzzCW3X7oj1ebhZ+Pu5Pam1OVjGtbh/9RK+V1zBhZ9+2qSzWyLHksM02bK/kB83HmT5rhzqbA4MIC7Mj6ToIDzcDGwOE5vdQUZ+OWv25AEQHezDJUOSGN2z9e8mXhIREddqcZ+oV+7KxdfTnZ6JYb97fM1//4u9ro5Bd999xDG+TzmAw2Fybt+Ev3w+9bvvyN+xg/NmzDim39aZhoHhsFNTXNyk1/tFRlLu5kPHoKP7lrzYMwD/8vyjGkNOTrklVTw+awN7skvpX5mG28N3EtKlM+evW0d0z54NvtY0TdbsyePrdRls3FeACSS615C8eznmDx/jVZxLYOvWtD3zTNrcdhltRo4kICamwTF/EdiqFQm9ezD+559ziqtYuOUgP/q5sadNd7IrCtg84z3C772XtiNHMmDqVDqMG9dgnt0sBvdM7EVljZUN503F+90HCbz7bsb85z9O/rZEji3TNFm+M4fpi3dxqKgKPy93RveMY1jXWPy93Nm5Yx/bU7ZTk5+LNSeL2sz9xB3cQxscFLfpQXanwbxUUs0nc1K4flhbhgxuOMMiItI0La6B2pBeQN92EXi6//9ld6ZpsvGdd+g0YQJh7dsfcYzlu7Lp2SaM2FC/v3x+03vvERgXR7eLL3ZZ3X9UXm2l2g6Btjrqfl6Et9F8/KjwCyXmD5cyNobV7qDEK5iEvIwmjyEnp4MFFdzz4Wpq6mwM3/EtdZ+8yYCpUznj6aePeGnbvtwyXv9xO1v2FxHm606/yjT48GXcD6UT3qkTnW+5js7nn090794u+ZIiOsSXKcM6cunpHVi5K4cZy/aw58K7sV54Iz7vP8In48fTaeJExk+fjk9IyGHH8XCz8M8L+zJ1+goyrniIFU9Noc3IkXQYN+6oaxQ5Wh8sTeWjn9JoExnAvRN7MahjFD/tzOa/X67jYIX9561CICAEAjpBe8A0Oa18N8OKdlCbtYj1G63s6n8ejy/YzwUzPuTal59y6iyyiIg4r0U1UKZpUlxRS3Tw7z+8VebmUpmXR+KIEU6N4+nuhqfb4S/LKTt4kMhu3Ro9vXhjfLoiDYdpErFpEWHt7230603TZEN2FfbuXgzs2PR7oDalF2Bz96Sjj7XJY8jJJ6ekivtmrMHhMBmV8hH5n3/I+OnT6X311Q2+rsZq573Fu/l6bTq+HhaGHFqN7a1ncLcYdL/0Uvrf/jmxffses7othsGQzjEM6hTN3I0HeW3udsy73uSCyk2svf9u3h06lCsXL8YvIuKwY/h6uXP7Od25+/1VlI29kgX33kvSWWcd038PRI7ko5/28NFPaZzVqzW3n9OdihorT8xcwbrMCvwOpdF+zxqGDOrOiCsmEdg6nvKKavZv3cG3u0tYbXTC7DmQey9IZpK3B1m793D7p9tYUB2Cz7hxTPriC7wCApr7EEVaFIdpUlheQ05JNT4ebsSF+//lPfMif6VFfWKosdqx2h1/um8pb9s2ACK7dXNqnEBfT0qq6g77fEVuLhFduza90CPIL6vm63UZdKrNIagslxAn7tH4o7xt2ziQ2IcgNzs9EsKO/ILDWJKyF0ttNcmdtZCo1Cssr+G+GWuosdoYtelzcj//kHPeeOOIzdOurGKen72ZzKJKupbtw+/1f+KGnUH3/IMBt9+Of3T0cTqC+kZqbJ94YkN9eeTT9cz07c7fZs9h/kXj+eCMM7hy0SJ8w8MP+/ru8aEkt4tgh9u5BMydweYPPqD3Ndcct/pFfuvb9ft5f0kqZ3Rvxe3ndGf9nhye+2wNlXaDtj99wpSx/Qm++wVWrt/No+8uwf3gHtxSluGXvhX36goSB45nzehrueze6UwhldH33MmFI7vz/hIvdn73Gu8NG8aUH344rhkVaYn255czf3Mmq1NzySmp/nWSo19EBfmQFB3IeQPa0P0oPnvJya9FNVClPzc9f1yvqbENVLCvJwcLKjBN80+XD5mmSWVuLn6RTZ8avCGmaTJ94S5ME9qnfItn165NWldq09xFlLTrw3ldY5p8M7DdYbJ2XyEhe1JIGDehSWPIyaWq1saDH62luKKWsRnzyfnoHca89BLJN9542NeYpsknK/bywZLdBBg2un/0JMHpW0i++WaG3n9/k7LkME2stXVUHcqiOD2dsoMH8Y+OJrxzZ4Jat3Z6YodeieE8d/lpPPjRWl5MtXDHZ1+z4KLxfDhqFNeuWtXgVOVXjejI397Op/zCW1j2xBP0vPLKo14DTqSx9uWW8caP2xnQPpK7xvdg2ca9PPP9TnzyDjG+KIVezz7M8wvSqfpoAwDeNQa17QZhth8KQM9QC7eGlbHtUC4fWlrx9abd7GrXjj6PPInF6EzEM9MovONi3hk0iBs3bMA7OLgZj1akeezKKuZ/c7eTeqgUN4tBr8QwEiICyMvJoyCnAKOqksC8DGrcPFnXqhMrdueSYJbx6JSBxLRLbO7ypQU66gbKMIzWwAdAFGAC00zTfLkpY3n8fNlddZ3td49X5udjWCxOf1DrHBfMwq1ZLN2ezfBusX963i8qipyNG5tSYoNM0+SNeTtYtO0QZ8caFMz9hiH339/ocWy1tXy6JR8jqQ3nDu3S5HpW7sqhzGbQbe86ons93ORx5PhzZa5+6/Uft7M/v5zrYqrYdv9/OP1f/+K0qVMbfM0v92V0suYR+Nzf6DTidMb9mOr0zJIllbWsS8tnw758MvNKyS0opdxm4DAMPMqL8SrNx6skj5DUDwnbsQJvT3daDxpE/9tuo8M55xyxmeoQG8wLVw5k6rsr+aEyjEtmzuSzCy5g2yef0Ouqqw77uvYxQfRtG85Bj6GEzPw32SkptOrf36ljkhPPscrU0VqXlo/NYXLnuT1ws1j4/odVuNV48vCgUHpPfof731tOldVBt7lvMHbS2Qy+fQr4+rN1by7//Hwje0psdLtlEj0Mg8WvLyXgrAnEl24h5anH8XjwUzzj4hjx6qt8ffXV5O/YQetBg5r7kOUk0VIz9UeF5TX86+N1eHm4cePoLliA6Qu2U+swsNTVEFBUSFVkPKVB0Qwt3kJS/mLWVPqxpcsoHr3i71x8Zg+GP/KIlgiQ33HFGSgbcJdpmhsMwwgAUgzDmG+a5o7GDhTq70WQryfpueW/e9zi5obpcPzlGaW/MrZPAgu3ZPG/udvomRhGiL/Xr88ZhkHyzTez6IEHyNu+nUgXXcpXZ7Pz3x+2MX9zJuO6RlB90zgiu3Xj9H/+s9FjzX79Q7I6DGRUjEGrw0yEcSSmafLJT6n4FB3i9L7ttGjoicdlufrF6tRc5m3O5PyeUaRdNZq4005j2EMPNfiaL9ek1zdPpfsIeeF2+t18M2e/8soRz9TU2ews2JLFvE0H2ZVVggl411bilbkHz/Ii2nkZhLdrizU6kor4LuTVdSOt++kcuPAOkqqzqZ37IfvGjyc0KYkBU6fS94YbGrwRPj4igAtPa8sHS1O5/LozCO/cmXWvvdZgAwXQu004KfsKiPcPIX3xYjVQJzeXZ8oViitr8fF0I9jPC1ttLXlpe/GLTaL35MsAKNqUgqdnKLe98C8ST+uPtaqK9a++wrT1+dBpELGzXmJO1hy6/OtJDhRUcMuYLiR4Xs6Gg6XUWu30T4rE11Z/OavFRYvFi/ysRWbqt0zT5MXvtlBjtfPvKwfy0+YMPlx1gMD0rfRIW8FpfTtSXnSI3O3fsbbHuSyO6svuCgvXJftS6eVJ4YiLWPrYFBxWKyOffFJNlPzqqBso0zSzgeyf/7/cMIydQCug0QEyDIM2kQHsyyv73eO/3NxtOhwYTlxi42YxuPPcHtzy1nJe+WEr/7qo7+/+0ve9/nqWPvooa195hXFvvNHYMv+kqKKGxz5PYWdmCVOGtMPx2M2UV1Zy4Wef4eHbuBn0aioq+WRvHT6B5dxy+QVNrmljeiFpeRW0WzaLvv97tMnjSPNwZa6gfjHOV37YRmKEP75vPUZBTQ0T33+/wYkT5m/O5M15O0gsySDkpTsY+cQTDH3ggQbfQCpqrHy7fj9fr82guLKWCEclbdYvwH/DYlp5Oeh15RV0vehaIrr8/syqaZpsP1jM/M2ZLNvhSdXE+xh8lRWPT19kzm23sem997jgo48I69DhsPue0D+RWav38cmKNEbdfDNzb7+dvG3bGrz095flEswhZ5GxaBFD7m38hC9yYnB1plyluKKWYL/6L/m2f/opVquNgLD6mST3zJmD3zfTqbv6KRZVBnLWlq18NGY0OzuP5NCIyZwT50aXMf1Y8eyzfEcCnnHdOb1LLDXuPSjqOABPw6RnYhh7t9ZPInSkBehFGqOlZuq3vkvZz7q0fG4e3YVP56awMKOSsJ2riA3z50DX4WwxQqntNpLgyHT6rJhJVoeB7O56Bp9OfwqPkKWUjLmByBvvYvnTTwOoiZJfuXQFScMwEoHewJqmjtE2OpCMvHIqa/5/1rhfvnmuLS11epz4iACuGN6BFbtzeXzWBnKKq359zjc8nO5TprDpvffY/vnnTS2V6jobX65J58Y3lrEvp4wHJvYk6INn2b9kCee8/joRnTs3eswnn/mA8sgELu8bha9X097sHKbJB0t3411TRqeaLFoNGNCkcaRlcEWuPl6eRlFFDRN9C9j33bec8fTTDTYjadml/OfbLSTYS4h++Q7OeuEFTn/wwQbfOPJKq7npzWW8t3g3rbwd9PrsSZIevYTTjRxuePu/TE3bw/CHH/5T8/TzMdItPpS/n9uDj/5+BhP7J7KyxIN9Nz7LpC+/oiQ9nXcGDqQy//Drmfl7ezA+OYHlO3OIGVt/z1/a3LkN/l7aRQfh7eGGvecgDqxY0eC2cvJwRaZcparWiuPnRXG3ffIJHj4+VBie7D5UwvrXXiO6LIuRXWP4dOVervl6P4uu+S+ZIyYzPCmEv105mo7X3UzuNQ+zL6oznUrSCPH3IiUtl/yeI2hnKcfT3Y3czZsBcD/C8gQiTdWSMvVbHy7dQ9fWIfh7wMKMSqL2phAaH8fWyO7kxnamNiQK98pSSiLbsKHXeELnvIdhraN48Hi6xtd/wZYZGEeva65h+dNPs3/p0mY+ImkpXNZAGYbhD3wB3GGaZtlfPH+DYRjrDcNYn9/Ah6CR3VpRZ3Pw/YYDvz7WZuRIAHbNnt2omi44rS1XDu/A+r35XPf6Ut5ZuIvK2vrG7IynniK6Vy9mTZrE19deS11FhdPjZhdX8ea8HUx+aSFvzttB26hAHujrz+6rJ7LpvfcY9sgj9LziikbVCvD2B3NY6xVPl4oMzj9veKNf/4tv12WwM7OEuB/eZuhdd+rbkhNYQ7lyNlNVtTZ+2HCAoZ1j2P/CY4R17Ej/v/3tsNubpsmb83fg52kh9uWp9Lx4EqfdcUeDdZZW1fHAzDVU1dq4r68/4fddSKuqXP62cyeXfvMN7c8+2+mJIXw83bl5TFduPqsrmzMKKe5yGlctW0ZNaSmLj3DJYbf4UEygzMMPw2KhpqSkwe3dLAbenm4YPn5YKysxTdOpGuXE5ar3KlcZ06s1uaXVfL5qH2EdOxKzcCamaXL7Oyv4acAUcv0imRRr485zezA2zp2obUtoO/sV8h64iZvPv40b313D3vg+dClK5Y4rR/POu9/xyk4rftWl3D5pMAeWL+enJ5+k44QJDX5pItJULS1Tv5UYGUBRRS0x5Vn45O6nuH0fIosOEFCURaBP/ZfUNr8gYkJ8mXLzpWQ8+QWmmzudK/az5rTJuFVXkBxgJbBVKwyLhfBOnY5r/dJyuaSBMgzDg/rwzDRN88u/2sY0zWmmaSabppkc0cAaLe1jgujdJpyv1qRTZ6tfODC2Xz9C27dny4cfNqouN4vB5KHtmX7LcIZ1jeGzlXuZ8tJCnpu9iV0VFq5YupShDz7Ipnff5c0+fchcvfqwY2UWVvD5qr3c9f4qrn51MV+vy6B/UiRPjWvPgB9f5ccxp1ORk8MFn3zC8IcbP2HD/FW7+DzDRnTmNp7+15VNbnpyiqt4Z9FuonJS6VCWQY/LL2/SONL8jpQrZzM1f0smlbU2+hv55G7ZwpD77mvwHqaN6YVs2V9Ez0MpuFeVM/LJJxuss6bOxkOfrCOnpJqbO3mwatI5BMTEcNWyZYR37Oj08f7RuL7xtA7z450Fuwjt1Jl+t97KhmnTyN2y5bCviQis/4a9oLwWN6/6e0qOxDTB+HkCG4dV66WdzFz5XuUqQzrHMLRzDB8sSaWk/1n47t/Jg0m1XDOyI/neIWy77nn+8U0qP2zYT45nCDEXXU72hbez45qn2d/nLBItVbwwvj3PvvB33lm4k88yDaLz9vLGXecQ6e/J55MmEZyYyMT339eXaeJyLTFTv3VWr9ZkF1exZf0Our7/T0L9vVkb3YfTD/7EZ3eeycAnLuDa2g30T4rktbnbsdtsdHn/nxzsMoxDJdV0+uwZRl4xie2ffkri8OFaCkB+5YpZ+AzgHWCnaZr/OfqSYNKgdtw/cw0Lt2Zxdu94DMOgx+WXs+ShhyjJyCA4MbFR44UHevOPCb2Y0C+R71MO8NPObBZuzSLI15OkLuMJeut09nzxMU/97REie/UifuhQApM6UFhp5WBhBRl55eSX1QDQJtyXc9t40y57G/lvvMT3ixZhcXdn6IMPMuS++/D092/08a5Py+U/81IJzErliVvOwrsJY0D9pXsvfrcFbFbiZjzFiP+9qGveT1CuypXDNPl6bQadWgVT8MHz+EZE0O3SSw+7vWmavL9kN2E+bthfe4yBf7u1wbyZpsmTX24k9VAJdwyJZf25pxMUH88VCxcSEBPT1LIBcLNYuO7Mzjz86XqWbD/E8IcfZuuMGSy8/34mf//9X74mIrB+spSCshrcvb2x1dQ4tS+Lpb6htNfVNThZhZy4jsV7lavcNb4HpVW1vJdWROIZU0idPZuL3zuXCf0SeeO/H7Fm+y7KrJXUxLbG8PFjWPc4+rWLoHebcOoqKvj2+xU89tVmSj0C6JSxlqdfuBOLAZ+MH09NSQmXzZ2Ld1BQcx+mnGRacqZ+MbhTNH5e7izKgw6Bfky9ajDXP/EZc7tdwMoXF1B21/us8fTFsS6DgX5VxKz7gp96n0m+6U+Xpe+T3DGWiuxsClNTGXjXXc19ONKCuGIWvsHA5cBWwzA2/fzYA6Zp/tDUAXu3CaNzq2DeXrCTbq1DaR3uT68rr2T5U0/x7fXXM2Xu3Cat19IhNpgOscHcenZX1qXls3J3Dhl55ewoclA98EIA9gKrDgAHUrHY6vAvyyewspCuh1LxXTsfj/wsCoACIKxDBwbfcw99rr+ekDZtGl2PaZp8tSadafO245N/kNt6BdK6V89Gj/OLD5aksimjkA4/TCOpe0e6T57c5LGk2bkkVwfyK8gqquSOc7qz8bb59Lr6aty9vA67fX5ZDbuyShhhTcdqOhj6wAMNjp9ZWMnaPXlcNaIjto9fxGG3c9mPPx518/SLXxaRzimuwqdHHGEdOlB64MBht/9lUUR7XR11FRV4BQQ0OH6t1U5ljRUPaxWGxaJZyk5uLn+vchUfT3eeuLQ/T3+5kVVcyns1Vcz7xxsMGNKHQRPOoo2tiDUvP4PdZsPi7cPB1u3ZFZnIa5EJFEcn4fDwwqe0lInuqVz1/FS2fPA+Pz31FFX5+Uz84AOievRo7kOUk1OLzdQvvDzcOP+0tny41Ebm1f9h7wuf0ZUi9qXuwSciAreSEux1dVhsdWyOSmRVp/rJu1ov+4ywFd8Q//DDzDjrLPwiI+l8QdMn9pKTjytm4VsOuPS6AMMwuPe83kydvoJ/fbKOl68ZTFB8PGf99798d8MNLH/mGU5/8MEmj+/p7sbgTtEM7lR/KtY0TSprbdTU2XFY6ziwbBmZy5ZAaRG2qiqslZV4BgYQeNVkAuPiCE5MpFX//kd1KrfGaufl77awaNshQneuYUprO6ff8nSTx1uy7RAfL0+j1Y5ltD20hUlr1jh9z4m0PK7K1e5DJQC08bGztqqK8CNMbJKWUz9Ri2XdEmL79cM3PNyp7ZPjg/jy44/pfN55BLVufbRl/2r7wSIAOsWFULJ/P5mrVzP8sccOu/2Og8UAhBUfxGG1Ej9kSIPj78oqweYw8cvYjk/79g02l3JiOxbvVa7k5eHGIxcns/NAIe+99y2p5TBrcy6ztuQD7eGWV3+3vcV0EFpXRi+jiDO7J9B/wHh2f/01b/TsSVlmJgnDhjHqm280Nb8cMy09U7+47PT2dIkN5OPpX7GnDNIjOkAgFAP8/BZn2KwE7t9OwsZ3CN+7gT4jBlI9dCgLH3iA6F69uHj2bHzDwprzMKSFccUZqGMiJsSXRy5O5p4PVvPoZ+t55rIB9LnuOjIWL2bJQw+RMHQoCaef7pJ9GYaBv7cH/t4egDeR540j+bxxLhn7r2QWVvDUFxvZl1tK/IIPGdvGh3OefbfJ420/WMS/v9lEWEEG7ea+xeRlS3WdrgD1DZSvlzte+QcBCE1KanD7tOwyLAZULZlDnztuP+L4+/MrsBgGVauXUFNcTM8jrLvUWOvS8vFyt9A9PpQ1z74JQM8G7uvbdrAIDzcLttX1l9YeqYHadqAIA7CsX9LgdOcix0vn+DCefegqyrKymPvPh1m9ZA220GhCk5KI6NKZyE4daRXsQ4SbFazu5G7ZQtqLH7JixQpMu524gQOZ8N57tBk5Uvc8ifysT1IUfZ66iZRp0/jqvmuhxwDajRxOSEw0tTmHCHGrJqCHL2aPIeRs8GXbxx/jGRDAGU8/zYCpU/HQDJbyBy22gQLoEhfCPyb05KkvN/LoZyk8cH5vxr35JofWreOzCy7g/JkzaTd6dHOX6bQ6m50vVqfz8fI0LDYrnWY8Rv/2UYx/a0aT3+j25pTyyKfr8S4vot07D3LxrE+I6t7dxZXLiSo9t5x2UYGUZewBOOKlpul5ZUT6WDBqqog77bQjjp9dXEWgrwepX36BX1TUrzNmukJeaTXLdmTTs0049rISUt58k8Thww97T5bDNEnZW0CHmED2vDuL1oMH4xUY2OA+NqQXEB/mS8WOrUReeJ7Lahc5WoGtWjHp3bcZumEDWz/6iIwlS8iZNY0Ch+NPi+xE9+7NkPvuo/055xB32mlqnEQOo+8NNxDWoQPz7rqLA4/cxQHTBMOon03oZ26engyYOpWhDz6os05yWC26gQIY1jWWqjobr/ywjVve+ol7Jvbi0u++47Pzz2fGmDEMvPtuznjyyRZ947dpmqzanctbC3dyqKiKtlWHCHntfroNGcBFn33W5IkeUg+VcP/MNThKiuk07R4uePUlks46y8XVy4nMZnfg7+2JWVf/5nCkhaj9vT2ocdR/+KouKjri+N0TQlmy/RBZCb2oLf0ah83WpPsT/yi/rJp7PlxNnc3O5YPb8NkFF1CZn89Fs2Yd9jVLth0iI7+c8WRQkJbGmBdfbHAfe7JL2XagiLODKyk1TWVHWqSYPn2I6dMHgJqSEg6tX4/pcODm5YW7lxdB8fEExMY2c5UiJ47E4cO5ISWF6uJi9i9bRvaGDVjc3AhOTCSsQwfCO3fWpCtyRC2+gQI4u3c8CREBPDd7E3e/v4qLBydx9eo1LLznH6z697/JWLyY82fMaJHz829ML+DdRbvZfaiE2AAPBi+bjrngSwbfdx8jH38ci3vT/gg2ZxTy8KfrsJQV0en1O7n4P8/Qy8WXT8mJz93NQp3dgaefHwDWqqoGt48L82dejR2CQhucLvwXZ/eOZ8HmTJbkdKOrhy+Zq1eTOGzYUdVcWF7DvR+uobSqjqcn92f3Y/eyf+lSzpsxg7jDLApdXWdj+qJdJAZ7UnT3XXSdNIkO4xq+DPfzlXvx9XIncM6HGG3aaMFpafG8g4Npe+aZzV2GyEnBJySEThMm0GnChOYuRU5AJ8wsA13iQnjt+qGc2SOOj5encfcnG4m753EmffkVJenpvNGzJ0sefZTasj+t4Xbc2R0OVqfmcu+Hq7lvxhqKKmo4P7CYpIcuwnvjUib/8ANnPv10k5un5TuzefCjtXiW5NPpf1O5+MVn1TzJX/Jwt2Czm3j80kBVVja4fevw+u18+p9OzsaNRxzfzWJw94Se2A030s67g13zFjS5VtM0WbMnl9vfWUFRRQ33j2jNljuuZ8NbbzHk/vvpMWXKYV/73uLd5JfVkDRvOp6+Ppz18ssN7iu7uIqfdmYzqkMomT/+QLdLL9VlTyIiIuKUE+IM1C98vdy5a3xPBrSP5LUft3PvjDUkRYcx7pufKHvlcZY+8girX3yRfrfcQv/bbnPZVMrOME2T/fkVLNqWxfzNmRRV1BLi58lY/xLsrzxCzt404ocO5fyZM5s8S5nDNPlwaSof/ZRGWMkh2r5xNxe8+hK9rrzSxUcjJ4uYYF+W78rBr1skACUZGQ3e25QUHYQBlA4YQ+0zt1Gcnn7E+6biwvy5cXQXXrE5eK04ly33/4dJl51DYlfnFtEtrqhlVWouX61J50BBBRHeBmft+Ja5D0zDcHPjzGefZdDddx/29V+vy2D22gw65++k6rtPmPDeew1OomKaJq/N3YaHmwX/L9+g0GKh9zXXOFWriIiIyAnVQP1iSOcY+rePZOHWLGat2sdLi9IJSL6a5LOvx2/hLH565hlWvfAC3S+7jEF3303EEaZubiq7wyQtp5SVu3JYviuHzMJKLAZ0D/NgcPYmqt94nZJDWbQaMIBxL79E+7Fjm/wtd1FFDc/N3szG9ALi09cT9/FzXPDedLpdfLGLj0pOJm2jA5m76SD21kl4BQWRvmgR3S655LDbRwb5cGaPOJZsN+gdEsWyxx5jwrtHniFyXN8EYnwMXphZzELvjiz8fDex0xcysGMM7bsm0bptPGEB3lTX2SiqqKWwvIZ9uWWsS8tnb279WePw2mK6LPmMwNVzKPH3Y+DddzPgttsIjIs77H7nbDzAa3O30yp/D0H/u48xL754xC8U5m46yNq0fCZG1pL+4XRGPP44oe3aHfEYRUREROAEbaCgfi2ns3vHM6ZXa1L25rNgSxYrd+dQ2/4cwp6fQFRBOovmzWJtv4F0OX0Qp91xB21HjTqqy3TqbHb251ewK6uYjemFbM4opKLGisWAdj52hpfuxPL9TGp2baXE05Oks8+m363v0vbMM49qvxv2FfDs7I1U1Vjp9tMMwlZ9xyXffq1r4eWI+rWLAGBDRhFtRoxg34IjX2J3+bD2LNl+iIobHmHz0zcy5P77CevQ4Yiv69slnplPTGHbtr188dUSNll8+SLXE3IPwKI/L35rOOwEZO4mPnUdwakphFYX0nHcODr9/ROSzjrr1/u2Duf7lP288sM2onJSSXjnAc774P0GL/MD2J9fzutzt9M9NoCye66kVf/+DLnvviMem4iIiMgvmqWBKt63j6K9e13yra/FMOiXFEm/pEiqam2s2JXD6tRcNjkMKib+HYCtFUV8/fF6gj9YTHzrKGIT44jv3J6YNq3xdHfDzWLgZjGw2U2q62xU19morLVRWF5DflkNBWXVZBZWcrCgApujfjazIOpoVXwAv51rcFv2PR5VZZh+frQePpyu999NxwkTjnoWl8paKx8u3cPsNelEeTlo/879BFfkM2Xp0l9nZRJpSGyoH61C/VixO5cLzzyTXbNnk7d9O5Fdux72NVHBvpw/oA2frnTQof/ZzJ06lUu//dape/YMw6B79yS6d0/CVlfH9hVr2bd9DwfSM8nNLcJSU4lvXRU+1kpC3U1iOiYRfsUYIrr8ndi+fZ2aTbOmzsY7i3bxzbr9RBzcRvvPnmHyN18fcUmDsqo6npi1AW9PN9p/+QIF1VVM/OCDJt+LKCIiIqemZvnkUFNewbS+fZn4/vsunf3E18udUT3jGNUz7tfL6zZnFJKeU0LankCyK23s9/CGPCAvE5ZmHnFMN1sd3lWleBdmEZW5B79De/HP2oNXSS4hiYlEdu9O3IP3kjhiBLHJyU2ekvy3TNNkyfZDTJu/k+KKWrpVZOD7+F207tmdixauOew6OCJ/5ew+rXl7wS4mnjMaz4AA5t11F1PmzGnwrOhlw9qzJ6eUDWNvpnrRTPxvuonxb73VqDOp7p6e9BwxhJ4jGl7MtjE2ZxTy4rebyS6pJmbV1/TctZCLF8wnNjm5wddV1lh54KO1ZBdXMSZ9Pjk/fMP4d94hvKNz92mJiIiI/KJZGqiqqARKTjuDTyZOZMi99zLi8cdd0nj8lpvFoGNsMB1jg39+pC+maVJZVcO+TdvYt2kbWWkZ1FRUUltVTW11DdTV4matxd1Wi4e1lgCLDT8vdzz9/PCLiiJoREeCE0YTnJhIRNeux2SdgH25ZbwxbwebMwpJCHSn69yXsa+cz6B77mHk44+36PWupGWa0C+Rb9fv54N12Vz5+BPMu2MqO2bNoutFFx32NZ7ubjx+ST9e+n4r85nCrJR5WK++hvPemubyrDojp6SKtxfs4qed2fiWF9D1839zxtlDGfXZZjz9/Rt8bWWNlQc/Xkt6bhmjDiym6O0XGfnUU5o4QkRERJqkWRoob093NvY8lw7te7P0xTvZN38+46dPJ7pnz2O6X8Mw8PfzocfgfvQY3O+Y7qux8suqmbF0Dz9uOoiftzujrHupvOsu/GKiGT9vHu1GjWruEuUE5enuxjUjO/H0lxspGTue6N7v8eMdd5A0ZgxegYGHfZ27m4W7zu1BZKA3M4EZefs5MPl6bp72Ij4hIcel9qKKGr5Zt59ZK/fisFlpvfgTOh9Yx3lvveHUPYClVXU8+FF983RG5k+UTvs3I598kqH3338cqhcREZGTUbM0UK1C/bhyeAc+XAoBD39KxZy3OdSnD70uv5wRjz1GUHx8c5TVLLIKK5m1eh/zN2dimiZDg2rx+t8/qEzdSa+rrmLMiy/iHRzc3GXKCW5Ylxi+XpvBtIW7ufe5V5g35nS+mDyZS2bPbvAeIMMwuGJ4RzrHhfDcp3Z+jExg7T1vcVb7IM6/5kICw8NcXmuN1c7q1FzmbzpAyr4CTAzCNy8haeXnjLz5Wgb+/a0jnnUC2JVVwhOzUiiprGXwxi8o/2I6I554gqEPPODymkVEROTU0Wx3T08e2p4+bSP435xt7Bh1PVEjL6Fy5nNsa9+ePjfcwNAHHjiu6zgdT6ZpsiOzmC9W7WPl7lzc3Sz0C3Hg/+HzVKxaSljPnlywcCFtRo5s7lLlJGEYBvef35s73l3Bi1uquOrlN1l72/V8d/PNnPvmmxiWhtfU7pcUyQf/GMsn367i69oEPq715/P/LKCDvYCxI3szZFgyPt5Nu7zUYZpkFlaSsjeftak5bD1QhNUBXmWFxG5cQGLWFoZNvoB+b6bgG3bkhs1hmsxem8E7C3YS7GnQ99PHMVM3MfGDD+h5+eVNqlFERETkF4Zpmsd9p8nJyeb69euB+g87i7dmMX3xbgrKakisziFo5gsE5aXT79ZbGXT33Q0uinkiqamzsXj7Ib5bv5+0nDL8vT0YHGzF66OXKVgyn5C2bRnxxBN0u/jiI36glePLMIwU0zQbnqmgGf02Uw05UFDBPz5YhZvF4IKCVWx78iF6Xnkl499+2+nZ6Kx2BwsXrOX7xZtI84jA4emN4bATWl1EmwALbVqF0SouioR2rQkLCfh14gmHaVJWVUdheS0F5TUcKixnZ3ou6YXV1Djqt/EpzCJ49zrC0zdxWu8kel5xBe1GjXK6toKyGv79Tf16aR09Kgl+8kbCI0OZ9MUXRPfq5dQYcny09EyB87kSaQmUKRHXO1yumr2B+kWN1c5Xa9L5YvU+yquttKrKJeCbdwjP2ELPyZeQfNNNxPbte9xrPVo2u4PNGYUs25nNTzuyqay1kRDmS6+6LMwZ/6UgZR1B8fEMvu8++lx7rSaJaKFa+htTY96U0nPL+MeHq/HxcGNc6QZ2PXwPHSdMYOJ77zX6ctGSvALmfDGPrXtz2F9toSgkDoeXj1OvtVhr8c1Jxy97L/6H0mjnVUfXAb1pe8YZJJx+eoP3Z/2R3eHgx02ZvLNwF1arjb77lmG+8ywdx41j4gcfHLd7tsR5LT1ToA97cmJRpkRcr8U3UL+oqrXxfcp+Zq/LoKCsBm97LaGblxCaMp92Qe50v/RSul1yCSFt2x7nqp1XWWNlQ3oBa/fksSo1l/JqKz6ebnT1sxG9fg4ln7yDo7aWyO7dOe3vf6fHlClqnFq4lv7G1Ng3pbTsUh7+dD0llbUM9Sik5p/XEdIqlgnvvkvi8OFNrqMiP5+Mbbs5kJ5JVlYexcXlWKursdfWYa+rxce0EuhuEuxlITw0gMjOnYno3JnwTp0a1TD9wmGaLNuRzQdLUskqqiTOUUbUmw/gX5rHsEceYdBdd+lsbgvV0jMF+rAnJxZlSsT1DperFreCpK+XOxcNasf5p7Vh/d58Fm09xEpPHw71GcPemnJWr15OyPuX0CXUg64Tx9P5vPMI69ChWWu22h2kHiphy/4iNuzLZ/vBYuwOE18PC0lGGeF7lmP9Zib20mJqwsNJvuEGelx+ObHJyY1aV0fEVZJignjthqG8Nnc7i7ebtH72a9w+epr3R4yg97XXMur555t01sY/IoJuIyLoNuIYFP0bdofJqt05fPRTGntzy4h0q6Pn92/gu2YePS+/nDOefprAVq2ObREiIiJySmpxDdQv3CwWBrSPYkD7KCprrazclcuq1FxS/IPI7Xc2qXYbC1J3EXDLE7Q2y+ndvxs9zhxO64ED8fD1PWZ1OUyT3JJq0nJK2XOolN2HStiZVUKt1Q5AlFFDl0O78Fz5A55bV2M4HPi0bUu3yZfQ/pxzaDd6dLOsoyPyR0G+ntx/fm+GdY3hlR+2sXLCvfQ4awprX3mQ1O++Y/S//023Sy5x+v6j46HWamfxtiw+X7WPzMJKwtzt9Fw6A99FX5A4dAhnrlpF3GmnNXeZIiIichJrOZ+MGuDn5cGonnGM6hlHnc3O1v1FpOzLZ3NaAPsSu5KFwWrAfUE+vp+8Rig1hPt5EBMZTOv4GKLjY4hOiCM8KhxfL/e/POtjmiZ1NgcVNVYqaqyUV1spLK+hsLyGgvIaMnNLyCysIK/civXnqx4tpoOg8jwiM7bhu30Ngfu341FVRmhSEq0HDybu5itoe8YZhCYlHd9fmEgjDOoYTff4MN5bvIsfNoBx93u0OrCZnAeepNWDD9L/1lvpc911+ISGNkt9VruDDfvyWbo9m1W7c6mqsxFhLaXbvA8JWDuPmF49OeOH72k3ZozO6IqIiMgxd0I0UL/l6e5G33YR9G0XAaO6UF1nY3dWCXszC9ixbS8ZXjHk2tzJ8PQFuwXSgfQ8IO/XMQy7DTfTgcV04DAsOCwWHIYbHObDl8Vai1dxLj6FWUQUHsKnIAvfnHSi3GoJax1HWKdORF5/EZHdHiW6Vy/8IiKOzy9DxEUCfDy4bWx3LhzYjtlr05nn4caB63uyr7qErfN+5PtXB9K9T2c6n3ceHc8995g3UyWVtaxLy2f93nzWp+VRUWvDy2ElYs9a2qz8jtCcNDpPnEiPh78jacwY3eckIiIix80J10D9kY+nO73ahNOrTTgXDO306+NWu4NDWfns3bGHvEN5FOYXU1xcRmV1LTabHavNjs1uYjjsGA47Focdd9OONza8seNjMQn2Mgj18SDA3xv/dhH4D+mOf/QoAmJjCYqPx93LqxmPXMT1YkJ8uXlMV64c3pEVu3JYvP0QG32DyRp2Mdsqign8OoXglz6gY4QPXQb0ps3wo79s1mGaZBdXsSuzmO2Zxew4WEx6XjkAPtYqglLXE7dxMeGZO0kaOZzOD91J5/PPb9KkEyIiIiJH64RvoA7Hw81CQnwUCfFRzV2KyAnH18v918tmiytqWZWay4Z9+WwICSW/95nsAebUVOD30Xp8Xv2GYB93oiKCiY6LJiIulpCYKIJbxeDt74/V7qDOZsdqc1BWVUdxZR0llbXklVZxIKeErJIa6hz1+/Ww1xGYtYfWezYQkrqeGHcr7UadSYdn7qfd6NF4+vk16+9FRERE5KRtoETENUL8vRjbJ56xfeJxmCZ7c8rYk11K6oECdqf7k1vegVw82G0YUAOkAWm5QO5hxzRsVjzLCvEpOkRY3gF88w4QkJ1GfKgPMb160XrS6bQ98zFC27fXfU0iIiLSoqiBEhGnWQyD9jFBtI8JYmyf+F8ftzsclFbVcehANnn7D1KSk0dpTh5VxSWYNZVQVYVZXYmvm0mglxv+ft74BAcT2DWOwNb9CYyLIzQpSZfFioiISIunBkpEjpqbxUKovzehXdpAlzbNXY6IiIjIMaOpq0RERERERJzkkgbKMIyzDMPYbRhGmmEY97liTJFTnXIl4lrKlIhrKVNyqjrqBsowDDfgf8DZQBfgUsMwuhztuCKnMuVKxLWUKRHXUqbkVOaKM1D9gTTTNPeZplkHfAJMcMG4Iqcy5UrEtZQpEddSpuSU5YoGqhVw8Dc/Z/782O8YhnGDYRjrDcNYn5+f74LdipzUjpgrZUqkUfReJeJaypScso7bJBKmaU4zTTPZNM3kiIiI47VbkZOWMiXiesqViGspU3IyckUDlQW0/s3PcT8/JiJNp1yJuJYyJeJaypScslzRQK0D2huG0cYwDE/gEuAbF4wrcipTrkRcS5kScS1lSk5ZR72QrmmaNsMw/gb8CLgB003T3H7UlYmcwpQrEddSpkRcS5mSU9lRN1AApmn+APzgirFEpJ5yJeJaypSIaylTcqo6bpNIiIiIiIiInOjUQImIiIiIiDhJDZSIiIiIiIiT1ECJiIiIiIg4SQ2UiIiIiIiIk9RAiYiIiIiIOEkNlIiIiIiIiJPUQImIiIiIiDhJDZSIiIiIiIiT1ECJiIiIiIg4SQ2UiIiIiIiIk9RAiYiIiIiIOEkNlIiIiIiIiJPUQImIiIiIiDhJDZSIiIiIiIiT1ECJiIiIiIg4SQ2UiIiIiIiIk9RAiYiIiIiIOEkNlIiIiIiIiJPUQImIiIiIiDhJDZSIiIiIiIiT1ECJiIiIiIg4SQ2UiIiIiIiIk9RAiYiIiIiIOEkNlIiIiIiIiJPUQImIiIiIiDjpqBoowzCeNwxjl2EYWwzD+MowjGAX1SVyylKuRFxLmRJxLWVKTnVHewZqPtDNNM0eQCpw/9GXJHLKU65EXEuZEnEtZUpOaUfVQJmmOc80TdvPP64G4o6+JJFTm3Il4lrKlIhrKVNyqnPlPVDXAHNcOJ6IKFcirqZMibiWMiWnHPcjbWAYxgIg+i+eetA0za9/3uZBwAbMbGCcG4AbAOLj45tUrMjJwhW5UqZE/p/eq0RcS5kSObwjNlCmaZ7Z0POGYVwFjAPOME3TbGCcacA0gOTk5MNuJ3IqcEWulCmR/6f3KhHXUqZEDu+IDVRDDMM4C7gHGGaaZpVrShI5tSlXIq6lTIm4ljIlp7qjvQfqVSAAmG8YxibDMN5wQU0ipzrlSsS1lCkR11Km5JR2VGegTNNMclUhIlJPuRJxLWVKxLWUKTnVuXIWPhERERERkZOaGigREREREREnqYESERERERFxkhooERERERERJ6mBEhERERERcZIaKBERERERESepgRIREREREXGSGigREREREREnqYESERERERFxkhooERERERERJ6mBEhERERERcZIaKBERERERESepgRIREREREXGSGigREREREREnqYESERERERFxkhooERERERERJ6mBEhERERERcZIaKBERERERESepgRIREREREXGSGigREREREREnqYESERERERFxkhooERERERERJ6mBEhERERERcZIaKBERERERESepgRIREREREXGSGigREREREREnqYESERERERFxkksaKMMw7jIMwzQMI9wV44mIciXiasqUiGspU3KqOuoGyjCM1sBo4MDRlyMioFyJuJoyJeJaypScylxxBupF4B7AdMFYIlJPuRJxLWVKxLWUKTllHVUDZRjGBCDLNM3NTmx7g2EY6w3DWJ+fn380uxU5qTmbK2VKxDl6rxJxLWVKTnXuR9rAMIwFQPRfPPUg8AD1p2+PyDTNacA0gOTkZH1bIac0V+RKmRL5f3qvEnEtZUrk8I7YQJmmeeZfPW4YRnegDbDZMAyAOGCDYRj9TdPMcWmVIicZ5UrEtZQpEddSpkQO74gN1OGYprkViPzlZ8MwMoBk0zQLXFCXyClJuRJxLWVKxLWUKRGtAyUiIiIiIuK0Jp+B+iPTNBNdNZaI1FOuRFxLmRJxLWVKTkU6AyUiIiIiIuIkNVAiIiIiIiJOUgMlIiIiIiLiJDVQIiIiIiIiTlIDJSIiIiIi4iQ1UCIiIiIiIk5SAyUiIiIiIuIkNVAiIiIiIiJOUgMlIiIiIiLiJDVQIiIiIiIiTlIDJSIiIiIi4iQ1UCIiIiIiIk4yTNM8/js1jHJg93Hf8ZGFAwXNXcRhtNTaWmpd4NraEkzTjHDRWC7XgjMFLffviOpqvFMmU9Cic9VS/4601Lqg5damTLUMLfXvB7Tc2k6Vuv4yV+4u3EFj7DZNM7mZ9n1YhmGsb4l1QcutraXWBS27tmOgRWYKWu6fg+pqvJZc2zHSInPVUv8cWmpd0HJra6l1HUPKVCO11NpO9bp0CZ+IiIiIiIiT1ECJiIiIiIg4qbkaqGnNtN8jaal1QcutraXWBS27NldrycfaUmtTXY3Xkms7Flrq8aquxmuptbXUuo6Vlnq8LbUuaLm1ndJ1NcskEiIiIiIiIiciXcInIiIiIiLiJDVQIiIiIiIiTmr2BsowjLsMwzANwwhv7loADMN43jCMXYZhbDEM4yvDMIKbuZ6zDMPYbRhGmmEY9zVnLb9lGEZrwzAWG4axwzCM7YZhTG3umn7LMAw3wzA2GobxXXPXcrwpU07V1OJypUy1XC0tU9DycqVMNd6pnCloeblSppyjXNVr1gbKMIzWwGjgQHPW8QfzgW6mafYAUoH7m6sQwzDcgP8BZwNdgEsNw+jSXPX8gQ24yzTNLsBpwK0tqDaAqcDO5i7ieFOmjqwF50qZaoFaaKagBeVKmWqyUzJT0GJzpUw5R7mi+c9AvQjcA7SYmSxM05xnmqbt5x9XA3HNWE5/IM00zX2madYBnwATmrGeX5mmmW2a5oaf/7+c+r+srZq3qnqGYcQB5wBvN3ctzUCZOrIWmStlqsVqcZmCFpcrZaqRTvFMQQvMlTLlHOWqXrM1UIZhTACyTNPc3Fw1OOEaYE4z7r8VcPA3P2fSQv6S/pZhGIlAb2BNM5fyi5eo/4fZ0cx1HFfKlNNafK6UqZbhBMkUNH+ulKnGe4lTMFNwwuRKmXLCqZwr92M5uGEYC4Dov3jqQeAB6k/fHncN1WWa5tc/b/Mg9acpZx7P2k40hmH4A18Ad5imWdYC6hkH5JmmmWIYxvBmLsfllKmTnzJ1fLXUTIFy5SrK1PHXUnOlTLnOqZ6rY9pAmaZ55l89bhhGd6ANsNkwDKg/TbrBMIz+pmnmHMuaGqrrN/VdBYwDzjCbd6GsLKD1b36O+/mxFsEwDA/qwzPTNM0vm7uenw0GxhuGMRbwBgINw5hhmuZlzVyXSyhTLtFic6VMHX8tNVMN1faLFpQrZapxTupMQcvNlTLlGspVC1lI1zCMDCDZNM2CFlDLWcB/gGGmaeY3cy3u1N/IeAb1wVkHTDZNc3tz1gVg1P/L9z5QZJrmHc1czl/6+RuIu03THNfMpRx3ylSD9bTIXClTLVtLyhS0rFwpU013KmcKWlaulCnnKFf1mnsSiZboVSAAmG8YxibDMN5orkJ+vpnxb8CP1N+k91lLCM/PBgOXAyN//j1t+rnrF/mjFpMpaNG5UqakMVpMrpQpOUkoU85RrmghZ6BEREREREROBDoDJSIiIiIi4iQ1UCIiIiIiIk5SAyUiIiIiIuIkNVAiIiIiIiJOUgMlIiIiIiLiJDVQIiIiIiIiTlIDJSIiIiIi4qT/A3yDFlwBorqeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x216 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "f, axarr = plt.subplots(1,4)\n",
    "\n",
    "for ax, N in zip(axarr.ravel(), [5,10,100,200]):\n",
    "    X = simulator([-2,2],N)\n",
    "    ml,tr = analyze_data(X)\n",
    "    ax.contour(grid[0],grid[1], ml, colors = 'maroon', levels = [0,1,4,9])\n",
    "    ax.contour(grid[0],grid[1], tr, colors = 'steelblue', levels = [0,1,4,9])\n",
    "f.set_size_inches(12,3)\n",
    "f.set_tight_layout(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bea054c8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e89476b9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}