Demonstrate optimization of a two-component Gaussian model using a loss function based on Gaussian KDE PDF estimation

demo_kde_loss_double_gaussian.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3d6f3934",
   "metadata": {},
   "source": [
    "# Optimize double-Gaussian model using KDE-based loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cbdb358d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax import random as jran\n",
    "ran_key = jran.PRNGKey(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4325f3b4",
   "metadata": {},
   "source": [
    "## Generate a target dataset from a double Gaussian with $10^5$ points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8759d5d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from kde_experiment import mc_double_gaussian\n",
    "\n",
    "mu1, mu2, sig1, sig2, frac1 = -1.0, 1.0, 0.5, 0.5, 0.35\n",
    "fid_params = np.array((mu1, mu2, sig1, sig2, frac1))\n",
    "\n",
    "NPTS_TARGET_DATA = int(1e5)\n",
    "\n",
    "ran_key, target_data_key = jran.split(ran_key, 2)\n",
    "pop1, pop2, frac1 = mc_double_gaussian(fid_params, target_data_key, NPTS_TARGET_DATA)\n",
    "\n",
    "TARGET_DATA = np.where(np.random.uniform(0, 1, NPTS_TARGET_DATA) < frac1, pop1, pop2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9908a962",
   "metadata": {},
   "source": [
    "## Build loss function and evaluate on $p_{\\rm init}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "601336f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from kde_experiment import build_loss_func\n",
    "\n",
    "NCELLS = 100\n",
    "loss_func = build_loss_func(ran_key, TARGET_DATA, NCELLS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "808fbdea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Array(0.02030281, dtype=float32)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_init = np.array(fid_params) + 0.2\n",
    "ran_key, pred_init_key = jran.split(ran_key, 2)\n",
    "\n",
    "loss_func(p_init, pred_init_key)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f63771",
   "metadata": {},
   "source": [
    "## Test evaluation of gradient of loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "08a0efe1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(Array(0.02030281, dtype=float32),\n",
       " Array([-0.00609772,  0.01659485, -0.027687  ,  0.04816183,  0.09830149],      dtype=float32))"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from jax import value_and_grad\n",
    "\n",
    "loss_and_grad_func = value_and_grad(loss_func)\n",
    "\n",
    "loss, grads = loss_and_grad_func(p_init, pred_init_key)\n",
    "loss, grads"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5064315a",
   "metadata": {},
   "source": [
    "## Walk down the gradient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a54c9cb9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEUCAYAAAAMdcB4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6hUlEQVR4nO3dd3ib9b03/rembXlJ8l6JLWeSrdhAKCQEHHYpw0no85xSWojdPj1tT1cMLeeU/k6Bxu3p7mkt6B5A7EIphQJWGkhYSWxl70jejqcsb1vW+P0h6Y5lyUlsS5ZsvV/X1auRbo2v7kvo7c933SKn0+kEERHRDIlD3QAiIpofGChERBQQDBQiIgoIBgoREQUEA4WIiAJCGuoGhEpycjJyc3ND3Qwiojmlvr4eXV1dfo9FbKDk5uaipqYm1M0gIppTCgoKJj3GLi8iIgoIBgoREQUEA4WIiAKCgUJERAHBQCEiooBgoBARUUAwUIiIKCAYKFPUNTCKEy29oW4GEVHYidhAaW1thUgkwlNPPTWl5z2334QHfvkBeBkZIiJvEbtSPjMzE62trVN+nlohh9XmwJDVjtioiD19REQ+IrZCmS5VrBwAYB60hrglREThhYEyRWqFK1B6hhgoRETjMVCmiBUKEZF/DJQpUseyQiEi8oeBMkWeLi/z4FiIW0JEFF4YKFMUHy2FRCxCD7u8iIi8MFCmSCwWQaWQwcwuLyIiLwyUaVAp5KxQiIgmYKBMgypWzlleREQTMFCmQa2Qc5YXEdEE827vEJ1OBwCora1FWVkZNBpNwN9DFSuHuYGzvIiIxptXFYrBYEBBQQFKSkqwdetWlJaWBuV9VAoZLENWbhBJRDTOvAoUk8mEiooKAEBBQQFqamqC8j7qWDlsDif6R21BeX0iorlo1ru8LBYLdDoduru7sWvXLp/j5eXl0Gg0MJvNAICSkpKrfu3i4mIUFxcDAGpqalBQUBCYRk+g8uznNWhFQrQsKO9BRDTXzGqFotfrodfrYTQaYbFYfI57xjyKi4tRUlICo9GIqqqqab1XRUUFKisrZ9hi/9Tcz4uIyMesVihFRUUAgEOHDvkNFJ1O51W1bN++HWVlZULVodPp/D6vqKgIWq3W53WUSmVA2++h4n5eREQ+wmaWl8Fg8LlPqVRCr9cLt6+m+0uv16OoqAgajUb4d6BxPy8iIl9hMyhvNpuhVqu97pt4+0oMBgO2bt2KLVu2ID8/Xxig98dzCWDP/6ZyKWBVrGvchKvliYguCZsKxV9X1vhjV9N9pdVq0dPTc1XvN91LAANAXJQUUVIxOvpHpvV8IqL5KGwqFKVSKczs8ph4O1yIRCJkKWPQamGgEBF5hE2gqNVqnyrFczsYg+ueLq+pdHWNl6mMQYtlOLCNIiKaw8Kmy0ur1foEh9lsDsqgOjCzLi8AyFLGYO/ZjgC2iIhobgubCgUAtm3b5rXupLq6Omjbp8xUpjIGHf2jGLXZQ90UIqKwMKuBYjAYUF5ejqqqKuj1epSXl3tNF66oqIDJZIJer4dOp0N+fr6wBiXQZt7lFQ0AaOvlOAoRETDLXV5arRZarRY7d+6c9DGXOxZIM+7yUsUAAFp6hrEwKTZQzSIimrPCqstrLslSugOFA/NERAAYKNOWnujq8uLUYSIil4gNlJmOoURJJUiNj0KLZSiwDSMimqPCZtrwbJvpGArgmunFCoWIyCViK5RAcK2W5xgKERHAQJmRLJVrtTwvBUxEFMGBMtMxFMBVoYzaHOgcGA1cw4iI5iiOoczAArUCANBkHkJqfHQgmkVENGdFbIUSCDnuQGk0c6YXEREDZQay3avlm8wcmCciYqDMQLRMgvSEaFYoRESI4EAJxKA8AOSoYxgoRETgoPyMXydHrcCHxu4AtIiIaG6L2AolUBaoFWjrG+F1UYgo4jFQZmiBWgGn07WNPRFRJGOgzNACTh0mIgLAQJmxnHGLG4mIIlnEBkqgZnmlxEUhSipGfTcDhYgiG2d5zZBYLMKy9Hicau0LQKuIiOauiK1QAmlFViJOtPZy12EiimgMlABYlZWI/hEbB+aJKKIxUAJgVVYiAOB4S2+IW0JEFDoMlABYnBYHmUSEEy0cRyGiyMVACYAoqQRL0+NxghUKEUUwBkqArOLAPBFFuIgNlECtQ/FYmZUIy9AYB+aJKGJxHUqAFOaqAQAH68xYmBQbsNclIporIrZCCbRFKXFIjJHhUL051E0hIgoJBkqAiMUiFOaqUFPfE+qmEBGFBAMlgApz1TB1DaKzfzTUTSEimnUMlAAqzHONo9Sw24uIIhADJYBWZiYiWibGgToGChFFHgZKAMmlYlyXl4R95zpD3RQiolnHQAmwW5enwtQ1CFPnQKibQkQ0qxgoAbZ5aSoA4F9nOkLcEiKi2cVACbActQJL0uIYKEQUcSI2UAK99cp4tyxLw8E6M3qHxgL+2kRE4SpiAyUzMxNOpzMogXLP6gzYHE68cKgx4K9NRBSuIjZQgmllViJuWpyMX79Xh5Exe6ibQ0Q0KxgoQfL5m/PR2T+KvxqaQ90UIqJZwUAJkg2aJCxLj8erRwK3ozERUThjoASJSCTC9ZokHGu2wGZ3hLo5RERBx0AJIu1CFUbGHDjT1h/qphARBR0DJYjW5SgBAIcbuaU9Ec1/DJQgylbFICU+CoZGS6ibQkQUdAyUIBKJRFiXo2SFQkQRgYESZOsWqFDfPYTuAV50i4jmt3kXKHq9HgaDAeXl5TAYDKFuDm7ITwIA/PNEW4hbQkQUXPMqUEwmE3bt2gWtVgutVotnn3021E3C6uxErMhMwB8/bIDT6Qx1c4iIgmZeBYpGo0F1dTUAV7gUFhaGuEWucZSHNyzE2fZ+HOSVHIloHpPO9htaLBbodDp0d3dj165dPsfLy8uh0WhgNrt+fEtKSqb8HlVVVTAajX5fPxTuXZOFZ944gz981IDrNEmhbg4RUVDMaoWi1+uh1+thNBphsVh8jpeVlUGj0aC4uBglJSUwGo2oqqqa8vsUFxcjKSkJer0+AK2euRi5BA9os/D2yTaYB62hbg4RUVDMaoVSVFQEADh06JDfQNHpdF5Vxfbt21FWVobi4mLhuL/nFRUVQavVCseUSiWKioqwdetWGI3GgH+O6dhemIPfvl+PVw634NEb80LdHCKigJv1Lq/J+JuRpVQqvaqMK3V/je9KUyqVQrdZOFiWnoA12YnYfagJn/1YLkQiUaibREQUUGEzKG82m6FWq73um3j7SkpKSlBYWAi9Xo9du3ahsrIykE2csW2FOTjb3o/bfrQPz+83hbo5REQBFTYVir+urPHHlErlFV9DqVQK3WOe7rXJeC4B7PHtb387KFdvHK94fTY6+0ex90wHnn7jNDYuScGStPigvicR0WwJmwrFXxdVMLusPJcA9vwv2GECAFFSCf6jaAl++5lrESuX4odvnwv6exIRzZawCRS1Wu1TpYwfZJ9P1LFyPHZTHt482YbahvAZ5yEimomwCRStVusTHGaz+YpdV9Pl6fKajcrEn8du0iBLGYOv7T6KwVFbSNpARBRIYRMoALBt2zavdSfV1dUoLS0Nynt5urxCFShxUVL8z7Y1aDAP4ZuvHOdVHYlozpvVQPFs2lhVVQW9Xu+zgWNFRQVMJhP0ej10Oh3y8/OFQfb56HpNEr5atASvHmnFp397EB39I6FuEhHRtImcEbpjYWZmJi5evDgrs7uuZHdNE5585QTkUjF23rEUD2/IDWl7iIgmU1BQgJqaGr/Hwmba8GzLzMxEa2trqJsBANhWkIOChSo89dop/NerJzFkteNzm/IBAIOjNijkEi6EJKKwF7GBEm40KXH47SOF+I+XjuB7/zyDo00WqGLlqKppxsYlKfjZJ9chRi4JdTOJiCY15UB57rnnYDKZsH37dqxduxaf//znYTKZoNFoUFZWhtzc3CA0MzJIxCL8cNsaZCqjUVXTjN7hMdy8NAV7zrTjnp/tx4rMRDx6Yx7W5ChD3VQiIh9TDhS1Wo1t27YhMTFRCJe33noLAPDyyy/PmUDxTBsOhzGU8WQSMZ64czm+fttSDI/ZkRAtw9sn26DbZ8K75zpR29AD/Vc3sVohorAz5UBRKpVITEwE4LruSFlZmXDMc/9cEE5jKP7IJGLIJK5JeLetSMdtK9JxsM6MbRUf4hd7L+Drty8NcQuJiLxNedrw+MHh6upqFBQU+D1GgXdtnhoPrMvCr9414hd7L8A8aMXImD3UzSIiAjCNCsVoNEKj0eBXv/oViouLkZCQAAB4/vnnodFoAt5A8vbte1dg1ObA9986i++/dRZSsQh3rEzHkrR4jNkdeOwmDRJjZKFuJhFFoGmtQ3nuuecgEonw2GOPobe3F9/73vcAAElJSfj6178e8EYGQzitQ5mO/ec7caFjAA3dQ3jZ0Iy+ERtEIqBwoRp3rkrHHz5sQMWn1nM3YyIKqICuQ3nuuedgNBrx0EMPAXBdttdkMiE/Pz9o26QEQ7iPoVzJTYtTcNPiFADAN+9aDofTiepT7fjSi4dxsN614eRfDjTiqXtXhLKZRBRBZjzLq66uDm+//TaAuTXLaz6RS11DYR9fkwmJWIRRmx1vnmjDP45dxJN3L4dUElZbthHRPBWxs7zmq7tWZQAAYmQSvHWyHc/tr8OHpm6oFDLcsSIdd7qPExEFGmd5zVM3L01FfLQUu948g1OtffjI1I3P/9mAlw3NoW4aEc1TUw4Uo9GI+vp6PP744z6zvOaSUF8PJdiiZRJ8YfMi3LsmE9Vf2Yh9OzfjY4uSsLPqGF4/dhFjdgdeO9oKy5A11E0lonlixrO8+vr68MwzzwAAkpOT58wsr8vNVJiv+kfG8MhvD6G2oQfZqhg09wzjvrWZ+PFD60LdNCKaIy732znt7eufe+45GAwGmM1mbNmyBY899tiMGjnbIjFQAGBkzI4nXj6Ow409WJQahz1nOvDWf2zk9GIiuiqX++2c1vSf2267DQaDARqNRnjxwsJC9PX1zaihFHzRMgl+tH0t3vnGZny/eA1i5VJ88+Xj+MXeCzjb1g8A6BsZww/eOottFR/iYu9wiFtMRHPFlGd5Pf/886isrPSZ0WWxWKDT6eZMlxcBqlg5vrB5EXa9eQY1DT34YfU53JCfhNqGHgxZ7ZBLxXj0dzWo/NwGxEbxSgdEdHlTrlBUKpXf6cFKpRJ5eXkBaRTNns/fnI/zT9+JmieLsHV9Nuq7B3Hfuiz844s3ouJT63GmrQ9feuEw7I6IvLAnEU3BjKYNT+VYuJnvs7ymQiYRIzkuCt97cDX277wFz9y/CiuzErF5aSq+c+8K7DnTge+8dhJWmyPUTSWiMDblfozu7m7U19f7rIg/cuQIjEZjoNoVdHN965XZ8qkNuajrGsJv3q/D3w634Iu3LMaOjdwElIh8TTlQduzYgW3btqGurk7YXdhzxcaXXnop4A2k0PvPe5bjpiXJ+P0H9Xj6jdOIlkvwqesXhrpZRBRmpjXSunv3bhw+fBg1NTWwWCx4/PHHsW4d1zLMVyKRCJuXpuKmRcko/WMt/uvVE1iaFg91rAxP/u0Enrz7GqzM4rY7RJFu2utQ/PnBD34wZ2Z5Reo6lJkastpw24/2IVomgUIuwbHmXixMUuAfX7wR8dG8DgvRfDft7etvv/32q34Tp9OJ2traORMoND0KuRT//YmV+MzvDgEASjdp8Pz+OpT+sRbf/vgKLE3nAkmiSHXZQHE6ndi1axeUSuUVX8jpdOLxxx8PVLsojG1elorHbsyDWCzCE3cuR25SLJ5+/TRu//E+LFArsHV9Nr546+JQN5OIZtllu7wOHz48pbGRurq6ObMWhV1egWUZsqKyphl/O9ICY+cAjj91O2S8DgvRvDPtrVemOtA+V8IE4DqUQFMq5NixUYPSTfkYGXMI27gQUeSI2P00uA4lONZmKwEAR5stnPlFFGHYJ0EBlaOOgTpWjiONFvQOjaGqthlVtc0YGLWFumlEFGQRW6FQcIhEIqzJTsSRJgtK/liDA3VmAMCbJy7iuYcL5tT2PEQ0NQwUCri1OSrsPdsJAHjizmWwOZz4/ltn8Y9jF/HxNZnC4wZHbYiWSSARM2SI5gN2eVHArV2gBADkJcfiszfmoXSjBmuyE/HU309i0N311WoZxsbyvfix/lwIW0pEgcRAoYDTLlBiSVoc/uvj10AmEUMqEeOpe1ege9CKP37UAJvdgS+/eBjdg1Z8aOwOdXOJKEDY5UUBFx8tw9tf2eR137oFKmxckgLdPhMON/bgUL3rEsQnWnthszsg5ZoVojmP/xXTrPnyrYthHrRiz+kOfOuu5fj3zYswMubA+Y4Br8c9/fop/Oa9uhC1koimixUKzZr1C1X40fY1WJwaj5VZiTB2uoLkWLMFyzMSALhW3P/m/XrkJinw2RvnzkJZImKFQrPs/nXZwoLHvKRYxEdJcbTZ1e3ldDqhP90Bu8MJY+cgegatIW4tEU1FxAYKt14JPbFYhFXZidh3rhMby/fi3/9yGG+eaBOmERsae0LcQiKaiojt8uLWK+FhdbYSHxi7ESuX4PXjFwEADxXmoKq2GbUNPbh1eVqIW0hEVytiKxQKD1sLsvGANgv6r23CHSvSAQD3rs3EiswE1DSwQiGaSyK2QqHwkJ8Shx9uWwsA+J9ta3DP2Qxs0CRh/UI1/nygAWN2B2QSMRwOJ/6n+iwkIhG+etvS0DaaiPxihUJhIzZKintWZ0IkEqEwV4VRmwN/rW2G1ebAN6qO4Rd7jfj53gtoMg+FuqlE5AcDhcJS0TVpuHFRMp782wnc8ZN9+KuhGY/ckAuRSIQ/fdQQ6uYRkR8MFApLMokY//tvWixKjcOw1Y5ff7oAT927ArevSMOLh5owbLWHuolENAHHUChsJUTL8Pd/vxEiEYTLCX96Qy7eON6Glw414pGPceEjUThhhUJhTS4Ve12b/to8Na7LU+Pney8IOxcTUXhgoNCcIhKJsPOOZegasOL5/Zf2+3I4nPjzgQa0942EsHVEkY2BQnPO+oUq3LEiHT/Sn8NTfz+JIasNf/yoAd965QSeeeN0qJtHFLHmbaCUlZXBYrGEuhkUJD/cvgaP3JCL331Qjzt/sh/f++cZyKVivHH8Ijr6WaUQhcK8DBSTyYSqqqpQN4OCSCGX4ql7V+DFkuvhcDohk4jwu88UYszuxAsHmmC1OfCXA4149p+n4XA4Q91coogwL2d5mUwmaLXaUDeDZsH1miRUf2UThqx2qGPl2LQkBT/Zcw6/eOcCrDYHAOC+tVnC9vhEFDyzHigWiwU6nQ7d3d3YtWuXz/Hy8nJoNBqYzWYAQElJyZReX6/Xo6ioCBUVFQFpL4W/aJkE0TIJAOA/71mO3TXxcDqdWJaegK9VHsX7F7oYKESzYFa7vPR6PfR6PYxGo9/xjbKyMmg0GhQXF6OkpARGo3FKXVcWiwVqtTqALaa5ZlFqPL5513J86+5r8OD6bOQlx+IDYzdOtvbia7uPClULEQXerFYoRUVFAIBDhw75DRSdTudVtWzfvh1lZWUoLi4Wjvt7XlFREbRaLXQ6HTQaDUwmE0wmE3bv3o1t27ZBqVQG4+PQHHBDfhL+drgFT/7tBA43WvDg+izckJ8c6mYRzUthM4ZiMBh87lMqldDr9cLtK3V/7dy5U/h3RUUFw4TwsUXJ+POBRhxutAAA3r/Q5TdQdPuMEItEeOwmzSy3kGj+CJtZXmaz2ae7arrdV3q9HiaTCTqdLhBNozlsgyYJAJCbpMCaHCXev9Dt8xiHw4lfvWvCb9+vn+XWEc0vYVOhXG7NiMVimVKlUVRUBKPReNnHeC4B7PHtb3+blwOeh1Sxcjx593KszVFi3/ku/Pxf59E7PIaEaCl+8PZZLFTHYk2OEmb39es7+keQGh8d4lYTzU1hEyhKpVKY2eUx8XYg8RLAkcPTjeVwAj/dcx4fmbrRNTCKX+w1Qh0rx/+7OV947LGmXqQljGLIasN17uqGiK5O2ASKWq32qVI8tzkOQoGwNkcJhVyCJ14+joFRG7KUMWixDONX7xqREh8F86AVR5oseOP4RVjtDrxXdkuom0w0p4TNGIpWq/UJDrPZLMwMCzRPlxe7uSKHXCrGD7etwaYlKfhYfhKqPr8ByXFydA1YcdOiZCxJi8eLh5pg6hpEc88wOvtHQ91kojklbCoUANi2bRuqqqqEacLV1dUoLS0Nynuxyysy3bEyA3eszBBu37smC795vw7XadSIkonxwsEm4diRJgu2XJMWimYSzUmzWqEYDAaUl5ejqqoKer0e5eXlXtOFKyoqYDKZoNfrodPpkJ+fL4QLUTA8vGEhrteoccuyNKzJVgIAHtBmQSoW4XBjT2gbRzTHiJxOZ0TunJeZmYmLFy9ydhcJWizDeOQ3B/GL/6vF13YfRXy0FH/Zcb1wvLbBjJOtfXh4Q27oGkkUYgUFBaipqfF7LKy6vGYTu7xooixlDKq/ugmAawD/ZUMz7A4nJGLX9PLvvn4ax5p7Ubw+Gwp5xP6nQzSpsBmUJwon6xYoMWi142RrLwDA2DmAw40W2B1OHG3qDXHriMJTxAYKZ3nR5WzIT0K0TIyHdB/hl+8YUVnTDHehAgPHVoj8iti6nV1edDkZiTH455c34pk3TmPXm2cAADcvTUFzzzAMDQwUIn8itkIhupK85Fg893ABfvLQWmSrYvDojXnQLlCitrEHETqXheiyIrZCIbpan1ibhU+szQIAtFqGsbumGaauQeSnxIW4ZUThhRUK0RRoF6gAAJ/4+fv4+M/eg53XqycSRGygcFCepmNRahw+tykfBbkqHG/pxfGWSzO+WizDGLLaQtg6otCK2C4vDsrTdIhEIjx+5zL0DFqx/rvV2HumA2tzlBiy2nDnj/fhlmWp+PFD60LdTKKQiNgKhWgmVLFyrM1R4p2zHQCAN463oW/EhteOXUSTeSjErSMKDQYK0TRtXpqKYy296BoYRWVNEzISoyEWAf/7jhEfGrvx/oUuXOwdnvT5rZZh1HUNzmKLiYIrYru8iGbq5qWp+J/qc/jvf5zCgTozvnH7UtR3DeKFg4144WAjAEAkAr5113JIxSL86UAjfv/Za5GljAEAfL3yKPpHbHjtizeG8mMQBUzEBopnUJ6bQ9J0rchMwK3LUvHqkVbIJCI8oM2CRCRCjlqBVVmJiJZJ8PsP6vHd108Lz9lzuh0Pb8jFqM2OmoYeREnZSUDzR8QGCgflaabEYhF+/UghugZG0Tc8hoxEV+XxpVsXC4+5Lk+NX75rRHy0FBXvmvCRqRsPb8jF0aZeWG0OWG0O9I2MISFaFtS2Do7a8M1XjuMbty9FtkoR1PeiyBWxgUIUKMlxUUiOi/J7TCwW4QubFwEAjjb1Yu/ZDjgcThys6xYe02oZRkJ6cAPlnbOdePVIK9blKPHIx/KC+l4UuVhvE82S6zVqmAetON8xgAN1Zsglrv/8Wi2TD9wHyvvGLgDAuY6BoL8XRS4GCtEsuV6TBADYf74TtQ092LwsBQDQYhnBz/91Ho/93v9Fi6bCZnfgtaOtcExYwf/+BVegXGhnoFDwMFCIZkmOWoFsVQy++/ppDFntuGd1JmQSEVotw/jniTboT7ejZ9A6o/fYc6YDX3zhMPa7AwQAmsxDaOgeQrRMjHMd/SHZ2NLucGJn1VEca7bM+nvT7InYQOHWKxQK3/74Cuy4KQ9P378St69IR0ZiDOo6B3GuvR8AcKDOLDx2cNTm9ePfZB7C796vu2wgeNa11I7bYt9Tndy/LhuWoTF0zzC0pqO5Zwi7a5rx2lFOhJno70db8fsP6kPdjICI2EDJzMyE0+lkoNCs2nJNGr519zX4v9cthFwqRqYyGu9d6MKY3RUSH5lcg/UtlmFc+7Qef/qoQXhuZU0TnnrtFP7u7tLqHRrzef2Gbtcq/cPjLgK2/0IXUuOjcOfKdADA+RB0e3mC7jzHcHzsPtSE37xfF+pmBETEBgpROMhUxmBg1LWhpCY5VgiUn//rPAatdvztyKW/6FssIwCAZ984g0d/fwjXP7sHXQOjqO8axOf+WIvBURsaul0/3EcaLXA4nLDZHdh/rhOblqRgSVo8AOB8R39QP1PfiG/QCYHCMRwflmErOvtH/R6z2hy4/pk9+MexuVHZMVCIQijbvWpepZDhvnVZONvej9qGHuyuaYZSIYOhsQcd/a4gudg7DHWsHG19I9h7thPDY3acaOnFP0+04c2TbThYbxbGSvpHbTjfMQBDowV9IzZsXpaKtIQoxEdLg/qj3moZhvb/q8bbJ9u87q93B0qLZRiDo5PvyPznAw34YNz4TySwDI1hyGoX/rAYzzxoRVvfCM5cDO4fAYHCQCEKoUx3oKzKVmJDfhKcTuDBX34AmUSEn31yHZxOoPpUOwDXj/UN+Un48fa1+N1nCgEApy724USrawv92voetPYO47ZrXF1bhsYe/OtMB6RiEW5cnAyRSITFqXFBrVBMnYOwOZzYXdPsff+4PcsuTNLt5XA48fTrp/HcftNVvdep1j6f2Wzh6ERLL275wTuwDPkfu/J0XfqrUizDruf0TPLccMNAIQohT6CszkrEuhwlHrsxD18pWoK/feFjuHFRMhYmKfD2yXY4nU609o4gSxmD+9Zl4ealqchSxuD0xX6cdF+T5Y3jF+F0AhuXpEAdK8e/znTgX2faUZirFlbi5yXHob7LdzfkY80W/HTP+RnPAGvvc1VT757r8Brjqe8exKqsRACTj6M09wxjyGrHhc4rV1DGzgHc9dP9ePlwy4zaOxs+MnXD1DWI+m7f826zO9Dvrkw63OduPM85tAz7diMCwNm2fnzhLwaM2uwBbPH0MVCIQmhpejxiZBLctDgZUokYT95zDb5ctBjL0hMgEomwZXkaPjB24WLvCKw2hxBAALA8IwG19WbUdw9BJLpUBeQlK3Df2ixUn2rHufYBYb0LAKQmRKF7cNTnL/vff9CAH1afg2HcYP50tLu758bsTrzl7vYatdnR0jOMTUtSIJeIJ62QzrT1AXAFy8iY/x/I2oYeOJ1OodvurQlda1fSZB6a8WecKs9ECX8VSt/IpW6uzgHfCqXXHSSTVTd7zrTj9WMXcaKlLxBNnTEGClEIpSVE4+R3bsd17kWPE61boMKY3Ym97uuuZCRGC8euyUxAa6/rB3zTkkuhsUAdi/+8Zzl++0ghPrE2E/etyxKOpcRFYczuFH6oPDzrQ377fv2MPk977wjio6VYmKTAa+6B5CbzEBxO19UuNSmxky6u9EyddjpdFchEh+rNePCXH2Df+S40ml3huf9856ThM9GOP9TgpvK9eOB/P4DpKqqgkTF7QK7A2ei+Ps7Ecw54B4X/Li/Xc3oG/Vconl0WTrb2+j0+2yI2ULgOhcKFWCya9NjKrAQAl8ZRxlco12TEC//+5LULAAAKuQTJcXKIRCJsXpaKnzy0Dqnxl0IoJd6159j4v4YHR2240DmAuCgp3jzRhrbeS10vVpsDH5m6fbrC7A4n9KfafX7M2/tGkZ4Qjc1LU1Hb0AOHw4k6dxdbbnIsFqXGTdrldaatHzKJ61z4G2cxuNfWnGrtE/7qHxlz4GVDC7784mG/IeRR3zWI6lPtuHt1BgBg79nOSR/r8eTfTuDR38189wJPoPT5C5Rx93X4CRTPc/yFEQC0umf+nWSFElpch0JzQY5KgfgoKT644JpOnDWhywsAUuOjsGlJCqRiERYmxUIkmjyghEAZ9+N1oqUXTifwjduXwu504i8HLq190e0z4iHdRz4L71490oLH/lCDu3+6HydaLv113N4/grSEaCxOi8OQ1Y7W3mHUdbl+6POSYqFJiUNTzxCsNodP28629WNDfjLEIv+Bcsz9Puc7+tFoHsKy9HjEyiX45ivH8eqRVrximHw8Zd95V4B8/balWJQaJ1xp83KMnQM43tJ71eNKoza7cF47+0fx6/fqYLM70Nzj6fLyDYXx40x+KxT38ckG5T0VyolxFUplTZPwB8hsi9hAIZoLxGIRrslMgNXuQIxMAqXi0q7EOSoF4qKkWJGZgGiZBNdp1Fibo7zs63l2Re4aV6Eca3b9GN21KgMbF6dgd00z7O4xljfdYxT//fppYY0M4Bq7UMfKMTBqwxdfOCz86Lb3jiA1IerSmpf2AZxrH0BynByJChmylNFwOl2D9xd7h/HSIdeFyEZtdtR1DWJVVgIWJsX6DZTj7naebx9AQ/cQFqfF47YV6YiLkiItIQpH3d12b51sg3nCbgD7znUhRx2D3CQFbl6SggMms9Cd1WQe8jsgbh60YmDU5vNa/vxEfx7XP7MHN39/L4asNrx6pAX//Y9TePtUu7Bo1V+V4bkvPlrqt0LxHB+y2v0OvLe4A+Vcez+sNgecTiee/ecZ/Fh/7optDgYGClGYW+meHZWhjPaqPsRiEb734Cr8R9ESAMDvPnMtvnvfysu+1vgKpa5rEK8eacHRZgsyE6OREh+FT16bg7a+Ebx7rgPNPUM40dKHf9+8CAvUCnyj6ihGxuwYGbNj37ku3L0qA4/fuQx1XYP40NgNh8OJjn5Xl9eSVFegnGvvh6GxB2tzVK7P4L5mTKtlGC8caETZX4+jvW9EmG68ND0B+SlxPoFiGbKi0TyEKKlrUL/FMoyFagWeuX8V3ivbjM1LU3G8pRemzgGU/rEWPxn3g2q1OfChsQsbF6cIXYFWu0Oo+j73p1rs/Osxn3NlHnAFib/ZWeONjNnxI/05SMQiDFrtaO4ZRnOP64f+xUNNwuMuN4ayODXOq0I5WGfGyJjdq0ts4s4IfSNj6B+xYWVWAsbsTpzv6EdzzzDMg1acvtjnd72Pze7As/88LSyADTQGClGYW5Hp6toa393lcc/qTKxxVyUyiRiSy4zHAEBCtBRyqRid/aP42Z7z+PKLR/D68YtYne16jVuWpSE5To4XDl7qNnlwfTaevn8lmszD+N+9F/De+S4Mj9mx5Zo03LkyA4kxMvzlYCPMQ1bYHE6kJUQjUSFDanwUDtaZYeocxPqFrkDJVLrGcy72jghjCxc6BnC2zTUgvzQtHotS41DfPQib/VK32HF3d9eWa9IwMuaA3eHEgiQFYuQSKBVyrMpOhGVoTOiae/34ReH5hsYeDFrt2OieuFCQq4JCLsG+851wOp0wdg6gtr7Ha+bbqM0uTOdt6B7EO2c78Ov3/G+P0tE36m6ba/1Pc8+QUDnsd3e1xUdJYRkew8iYHb99v05omycwFqXGodM9Q67FMoxtFR/ilcMtXiHUMyFQLrrHTzzrjk629AnVpsMJHG2y+LT19eMXUfGuSTjfgcZAIQpzngolM9E3UKZKJBIhJS4Knf2jON8xAJVCBqcTuDZPDQCQS8XYWpCD6lPtKH/zLJakxSEvORY35CfjE2sz8ct3jXjqtZOIj5Liek0SomUSPKjNxlsn23Cq1TUwnJbgCo0lafHC7DRPoAgVSu+wECjn2/txsrUXcqkYmpRYLEmLw5jd6bUexfND+YD20oy1hepLV55c4w7EFw42QSoWoWvAio9Mro02PdOEN+S7ZtJFSSVYnpGAc+396BwYxciYay3I+PcbP6uqvnsIv9h7AU+/fgoXe32vXeOZKu35jM09w2hxVyhOJyAVi7A0PR69w2N452wnvvPaKWHzTsvQGOKjpUhPjEH3oBU2u0PYVaC+exC9Q1bhujkTpw57xk9uyE9CrFyCo80WHGu2CBMbxm8Q6mqLExXvmpCfEoui5Wk+nyMQGChEYU6THIvMxGisyk4MyOslx0ehvX8Exs4B3LcuC/t3bsanb8gVjn/51sX45l3LkKWKwb9dv1C4/6mPr8DWghyMjDnwgDYLcqnr5+OT1+ZgzO4UVrinJbi61RanxcHh/kFd7W57bJQUiTEytFqG0Wh2/SBe6BzAiZY+LM9IgEwiFn7433HPxLLaHHj3XCcWJimwfqFaaM/CpFjh30vS4iGXimG1O/Bv1y9EXJRU2Nm4vXcECdFSr8ss5yXHoq5rEE3mS91ZhnE/wN2Dl7qfzrf342hTLxxOYPch1w4AY3YHPvHz9/CvM+3CrLiVWQmQS8WuQLEMC1O8s1UxUMXK0Tc8JgSSZ8fnvuExJMbIkBIfBafTdb+nTS09w+gdHkOO2hXCEysUTxWUrVJg09IU/OPYRRyoM+OajAQsSYtD7YT1NvvPd+HUxT6Ubsy/7MzCmWCgEIU5qUSM98pu8fpxn4mUuCicaOnDkNWORalxyFErvLrKomUSlGzMh/6rm/DwhlzhflWsHM/cvwo1TxbhO5+4NFazOC0eK7MSsP+8aw+u8RUKAKzISkS0TCI8PiMxGsaOQWFiwLn2AZxo7cVKd9deRmIMVmQmYM/pdvSPjOFTvz6Ag3VmfOr6hUiMkSEtIQpRUjFS4y9ddlkuFQuz3u5enYEt16ThzZNtcDpd4zqpCZemTgOuQGnvG8UZd9ePSAQcbrQIxz0D8TEyCfae7YDV7kB8tBQvHWqE3eFEq2UYR5t7sf98l7A7QEZCDLKVMTjT1o/e4THcty4LIpHrOjjKGBl6h8fQ5n6s5/Utw2NQKmRIibs0tuWp3Fosw7AMjyHXHZy9w74VilQsQkp8FD69IRe9w2M40mTB6mwl1i9UwdDg3Y33wsFGJMdF4RPrMhEsDBSiOSCQf1GmxEcJffOLUuIC8pr3rXV1RYlElwb+F6e6Xnv9ApXXYzOVMTjc5PrrOSFaisONPe7B5UsV2K3LXOtYvvXKCRyqN+NH29fgsZs0AIBl6QnIS471OScbNElIT4jGuhwlVmUlond4DD1DY2jvGxGqJo+8ZNeP9P5zXcJzx6+g9/zgr85OxMiYa7zjiTuXo7V3BPvPdwpdWo3dQ2jvG0G0TIyEGCmyVDGorXd1tV2TkYBPXrsAd61yjTNZhsbQ3jshUIasUMbIkepuX0f/CJrcr93cM4y+4TGhEptYobRahpGeGA2JWIRr89RCoK7JUUK7QIW+EZvXmp+zbf0oWKhClFSCYGGgEEWYlHF/2S9KDUyg3LsmE2IRkBQbBZm7z/+azASsX6gSFhN6ZCRGCz/SNy9NFabVrhoXKLcsT4PD6br41MMbcnH/umzh2HfvW4mf/591Pm346pYleOsrGyGViIUFoK2WYVeFEu9boQCui4+lxkdhgyYJ5zsGhKDtds/w0rrHRRanxgnjN8ebe4XupgbzENr6RpGW4JqBl61SYNDqmt6bpYrBM/evwievXYDEGBmGx+xC9TG+QklUyITxoPPtA8JjOvtH4XC6zpdcIvZZi9LaOyJ8TpFIhB035UEsAgoWqoTLTXsurjZmd6DRPARNSiyCiYFCFGE8gaKOlSMpLuoKj746qQnRuGVZKhalXvrBUsil+OvnbxAGqz3Gr/a/dXkqAEAmEWFx2qVwW52ViJT4KKTGR+Frty3xen6OWoFFqfGYSC4VIzHGNU7imRHX3OMJFO/P6elG6h+1YYFagdXumXKeiQXmQSskYhFWu0OuIFeNaJkEGYnRqOseFAKl0TyEtt5hoZsvW3Xps2WP+5yJ7vVDnj3IzOPGUJQxMiTFRWGBWoEjTRY0m4cQM66LMFEhg1IhE6YNO52uXQrOt/cjc9xWPA9os/HRE7ciNzkWOWoF8pJjhVlmjeYh2BxO5AeoIp1MxAYKt16hSJUSJwcQuO4uj59+ch2e/3ThFR/nmTocHy0VwmZperxXV4xYLMKv/m09fveZaxE/bjD9anne4/TFPlhtDp8xlBi5RPgxzlEroHFXLPXu9Rndg1aoFDIsy0iASATcuCgZgCuI6rsGhRlWVpsDp1r7fAJFLhELi0gBCEHnmYpsHrTC6XTCMjQmHFubo8RHpm50D1pRkKvyeq5SIRMqlBcPNeGxP9RALBLh7tXe4yHjP+dNi5PxkcmMUZsdRnfXV7ArFGlQXz2MZWZmorV1blwFjSiQPBVKfoC6uzwU8qv7OfFMHV6gViAzMQaJMTJh2u94EyubqVDHyhElFeOIey3GxAoFAPJSYtHaO4IctQKZyhjIJWJhyq55cBTqWDnykmOx56ubhC6y3GQF3jrZjhi5BGKRa73HoNWONPfrZ6tcXVeZymivMR5PaHiYB60YtNphcziF3Q/WLVDi7+6ZaddrkoRJDsoYGZQKubANy4fGbmQkRmPfzs1C96I/Gxen4A8fNqC2oUfYiVrDCoWIAunSLKzg/rhMJnNcoIjFIlR+bgO+cfvSgL6HSCRCljJG2I7Fb6C4Q2KBe5bbgiSFcKli86AV6lhXJadJiRN2KMhNioV50Iqzbf1eYz7pnmrHXaFkqbzXDI0PlPgoKcyDVmFdiTLG9T7rxk1euC7v0vToRIUMSvegPuDaWXhlVuJlwwQArs9PglQswv7zXTB1DiA5Lson2AKNgUIUYbJVCug+tR7bC3NC8v7pidGQSUTIdf+gL0mLh1IhD/j7ZCijhR/htAldXoDrYmPApRDITYoVdjHuHrQiKdY3hDxt7hqwoiBXDam7CvF0NSXHRUEuFfvsajD+h3x5RoI7UFxt84yvLM+IFxYxalLihBBUxsihUsjRM2TFkNUGU9egsHvC5cRFSXFtnhqvHW3F2fYB5Ae5uwtgoBBFpNtWpF91F1WgyaVi/OnR67DDPQ04WMbvLJCa4BsOm5emYNOSFKxwVxp5yQrUdw/C4XB6VSjj5Y5bTLlArRAqkXR3oIjFIvxo21qUbPT+bOMD07PZp2cNjGcBZJRUghVZCYiLkkKlkAmvnRgjgzLWNYZyrNm1M/SKzKtb5ProjXlo7hnG0SZL0Lu7AAYKEYXAdZokvz/YgeSZTRYXJfUbnpqUOPz+s9ciLsp1bGFSLEZtDteCwqExv+1bmHRpu5csZQwWuKf7jl/ncvfqDJ9ZaAnRrvcQiVwTEADggHv35vFTtx/esBAPb1godNnJpWJEy8S4IT8ZY3Ynfljt2vTyaioUALhlWarwWFYoRETT5Ol28led+OMZUznsHshPivMNFM/UYcAVWJ6A8delNp5UIkZclBQpcVFC+HxU141sVYxX2N2/Lhs771gGwHU5geL12RCJRLhpUTKyVTE4WGeGSiHzunLn5YhEInzp1sUALl0/J5gidpYXEc1vngrF34C8P57xEc9K98kqqNykWFzsHUGWKgYPFS5ARmKM19Yyk0mMkSEpzjUeAgBN5mHcsix10sfftSoDd61yLQoVi0X45LUL8P23zmJFZuJlL6I20e0r0vH2VzYKOxcEEysUIpqXPGtRJq6Sn0xGQjSipGL85WAjpGIRlqX7/4t+SVoc1LFyJMbIsDIrEV/YvOiqXl+TEotl6fFeg/2LpzDTbmtBNuQSMdbkTH2T0CVp8VMKoelihUJE85KnQpm4j9dkxGIRFiYpcK59ALseXDXptjRf2bJkWht1PvdwAcQikdeVFxf7WfE/mdT4aPzjSzd67TQQbhgoRDQvRcsk2PXgKhTmqq/8YLdHb8zDwKgd2wsXTPoYpUI+rWnOnm4xmUQEucS11f5U1wJ5dnAOV/MuUAwGAzQa15Q9s9ks/JuIIs/lgiEQj58OkUgEdawcbX0jAducM1zMuzGUHTt2IC8vDzt27IBaffV/mRARzRZVrNxnhtd8ML8+DYAnnngCxcXFoW4GEdGktixPhfPKD5tzZj1QLBYLdDoduru7sWvXLp/j5eXl0Gg0MJtdU/dKSkqm9Pomkwl6vR4GgwFFRUXQarUBaTcRUaB89bbA7l0WLmY1UPR6PSwWC4xGo9/jZWVlKCwsFCqMsrIyVFVVTani2LlzJwCgoKAAt956K2pra2fecCIiuqJZDZSioiIAwKFDh2CxWHyO63Q6r6pl+/btKCsrEwJFp9P5fZ6nEqmqqsKhQ4ewa9cuKJVKmEymoHwOIiLyFTZjKAaDwec+pVIJvV4v3L5S95dGoxG6uAwGA7Zt2xbYRhIR0aTCZpaX2Wz2mZU11VlaWq0Wer0eVVVVeOmll/yO0Xh4rtjo+R+v3EhENDNhEyj+urKu5thEJSUlKC4uFrq9JpOZmQmn0yn8byqBwvCZOp6zqeH5mjqes6kJxvkKm0BRKpXCzC6PibfDxXe+851QN2HO4TmbGp6vqeM5m5pgnK+wCRS1Wu1TiXhuX67SmC5Plxf/qiEiCoywGZTXarU+wWE2m4WZYYGWmZmJ1tbWoLw2EVEkCptAAYBt27Z5rTuprq5GaWlpUN6rvr4eBQUF03puRkbGtJ8bqXjOpobna+p4zqZmuuervr5+0mMip9M5azsAGAwG6PV6VFRUAABKS0t9VrOXl5dDq9UKa0imulKeiIhCY1YDhYiI5q+wGZQnIqK5jYFCREQBwUAhIqKACKtZXuFuplvrz3c6nQ61tbXYunUrAKCyshJlZWVeV82M5HM400s3ROK5u9w54/fNl+d8Aa5NeLds2TLl79GMzpmTrsrOnTudlZWVk94mp7OiosKpVCqdAJxardZZW1vrdTySz2F1dbWzsrLSWVJS4iwpKfE5fqVzE4nn7krnjN83Xzt37vS6rdFonBUVFV7Hg/k9Y6BcJaVS6XW7trbWWVRUFKLWhKfxX1x/eA5d/4H6+3G80rmJ5HM32Tnj981bT0+Ps7i42Ou+Xbt2OTUajXA72N8zjqFchavZWp8uj+dwclc6Nzx3Uxep50yv13tdB2r8daFm43vGMZSrEIit9SOFTqeDWq326X/lOZzclc4Nz93k+H27RKlUoqenx+u+6upqYfuq2fieMVCuwpW21g/G5pVzUUFBAZRKpTAounXrVqjVahQXF/McXsaVzg3PnX/8vl2exWKBXq/Hnj17hNuXe2wgzhm7vK7CXNpaP5S0Wq3XDJvCwkI8++yzAHgOL+dK54bnzj9+3y5vx44dqKysFLa2mo3vGQPlKsz21vpz1cS+Vo1GI/TL8hxO7krnhufOP37fJldeXi7slegxG98zBspVmO2t9ecik8mELVu2+HwhPX9B8hxO7krnhufOF79vk6uqqoJWqxU+qyd4Z+N7xkC5Sp6t9T2CubX+XKTRaHwuu/zSSy+hrKxMuM1zOLkrnRueO2/8vvmn1+thNptRUFAAi8UCk8nkNXsr2N8z7jY8Bdxa//JMJpPwZezu7kZ+fr7fVbiReA4DcemGSDt3Vzpn/L55s1gsUKlUPvcXFxejsrJSuB3M7xkDhYiIAoJdXkREFBAMFCIiCggGChERBQQDhYiIAoKBQkREAcFAIQoCi8WCsrIy5OfnQyQS+SzAMxgM2LJlC1QqFbZu3eq1QyzRXMVpw0RBVFZWhvLyclRUVPjM5/dcXW/nzp0hah1RYLFCIQoSi8WC/Px8FBcXC4vzxqupqUFxcXEIWkYUHAwUoiCpqalBUVERSktLYTAYfLq1DAaD1265RHMdA4UoSDyBUVRUJOw9RTSfMVCIZkFxcTF2797tdV8kbKNOkYWBQhQk4wOjtLQUFotF2MzQYDCgoKAgRC0jCg4GClEQTAwMT9eXZ3C+pqbGa6dhovmAgUIUBP4Co7S0FHq9/orX7/aoqqq6qscRhQsGCtEs8UwR1ul0V/V4f1ONicIZA4VoFpWUlODZZ5/1ql48F4rS6/XCFQf1ej1MJhN2797tdQU9nU4HvV4PnU4nXIBKpVKhqqoKBoMBZWVlXlfoI5pVTiIKqOrqamdlZaXfY7W1tc6J/9nt3LnTaTQahed6FBUVOXt6eryeW1JS4nXc8/+e5zudTqdGo5nxZyCaDmmoA41oPiktLcXu3buFBYsTV8JrtVqf+7Zv344tW7ZAq9Ve9vrdL730EpRKpVCBmM1m4ZharRb+7XkMB/1ptrHLiyiAKioq0NPTg9ra2km3VRl/fW/ANQPMaDSitLQUZWVlPivqPbctFgsKCwuh1Wqh1WpRW1vr9/XHhwvRbGKgEIWYZ5C+qKgIu3btEgLEs47FM9Nr69atOHTokPC88WMl40PIbDazOqGQYJcXURioqqqCUqmEyWQSdiUuLS2FTqeDRqOBVqtFUVGRMICvVCq9KpGamhqYzWYYDAafCohotnD7eqI5bsuWLaisrORWLhRy7PIiIqKAYKAQzWF6vR41NTV49tlnQ90UInZ5ERFRYLBCISKigGCgEBFRQDBQiIgoIBgoREQUEAwUIiIKiP8fxFd7DTfkWXgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "nsteps = 100\n",
    "\n",
    "p_best = np.copy(p_init)\n",
    "loss_collector = []\n",
    "\n",
    "# Take some big steps\n",
    "for __ in range(nsteps):\n",
    "    ran_key, pred_key = jran.split(ran_key, 2)\n",
    "    loss, grads = loss_and_grad_func(p_best, pred_key)\n",
    "    p_best = p_best - 1.0*grads\n",
    "    loss_collector.append(loss)\n",
    "    \n",
    "# Take some small steps\n",
    "for __ in range(nsteps):\n",
    "    ran_key, pred_key = jran.split(ran_key, 2)\n",
    "    loss, grads = loss_and_grad_func(p_best, pred_key)\n",
    "    p_best = p_best - 0.2*grads\n",
    "    loss_collector.append(loss)\n",
    "    \n",
    "fig, ax = plt.subplots(1, 1)\n",
    "yscale = ax.set_yscale('log')\n",
    "__=ax.plot(loss_collector)\n",
    "xlabel = ax.set_xlabel(r'$N_{\\rm step}$')\n",
    "ylabel = ax.set_ylabel(r'${\\rm loss}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f086e11f",
   "metadata": {},
   "source": [
    "## Generate MC realization of initial and final prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "44b68844",
   "metadata": {},
   "outputs": [],
   "source": [
    "pop1, pop2, frac1 = mc_double_gaussian(p_init, target_data_key, NPTS_TARGET_DATA)\n",
    "pred_init = np.where(np.random.uniform(0, 1, NPTS_TARGET_DATA) < frac1, pop1, pop2)\n",
    "\n",
    "pop1, pop2, frac1 = mc_double_gaussian(p_best, target_data_key, NPTS_TARGET_DATA)\n",
    "pred_best = np.where(np.random.uniform(0, 1, NPTS_TARGET_DATA) < frac1, pop1, pop2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "206e62a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD/CAYAAADllv3BAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfVElEQVR4nO3dXYgj55kv8H8vx3MgrElN4d2d6ckZ21IYZw1hRbVETOAw8ljCTC4WQtQt6KtdxiP5ZiCci1YGBmZYDG3VxcGwN1bhcGYvTkOPKuxdzKLqkyEXJkEtnU4uTCZsV9s+cU9s3NUFMyT2EOhz0a53qvRR+lapqv4/aDwlleS31a1Hbz/1vs+zdHJycgIiIoqMvwp6AERENF0M7EREEcPATkQUMQzsREQRw8BORBQxDOxERBHzX4IewHPPPYcXXngh6GEQEYXKRx99hC+++KLnfYEH9hdeeAG7u7tBD4OIKFTS6XTf+5iKISKKGAZ2IqKIYWAnIooYBnYioohhYCciihgGdiKiiGFgJyKKmMDXsRPR+DRNw9bWljheX19HqVQKcES0CDhjJwqxra0t7O3tAQD29vY8QZ7ii4GdKORSqRTu37+PVCoV9FBoQTCwE4WMpmnIZrPIZrNitk7kxhw7CaZpolarQVVVFAoF5PN52LaN/f19JJNJbGxs+J53dHQEAKhWqwOf8+joCLquo1AoeM6nwZz0SyqVQiqVwvr6urhvb28P2WwWAPPtsXYSsJWVlaCHMBe1Wi3oIQw9BgAnrVbLc1uhUDgplUoDz6vX6yeKogz1nMfHx13PSYNdvnz55PLly94b//faSe3aysnlv/+bk8uXL59885vf7D6HIsUvdg41Y1dVFYlEApZlAYDvLEDTNLRaLayurgIA6vU6KpUKEonEpJ9BoWWaJmzbDvUYyuUy8vk8arWa73mFQgGNRgOqqooZfj+SJGFlZWXsMcWJe/WLM1vHVtFzTulKEqUrSWB9W8zaKZ4G5tidoFwoFFAqlbC/vw9d130fc+/ePeTzeVQqFZTL5VgHdQAzTzUME7Dnme5YXV1FpVLxPcc0TQCI/e/GsNyrXzrTL0SdBs7YNU3zBIVisYhKpYJCodD3McfHx9MZXQQYhgHTNGGaJiRJQiKRQC6XAwC0220Ap0Gu0WigWq1CkiQYhoFyuYxyuQxJklCr1bCzswNJktBut2EYBhKJBBqNBvL5PJrNJqrVKkzThK7rSCQSaDabuHnzpni+fmMYVr1eHzpf69SJNk2zb+A2DAOlUmnkccSZs/pF2NoJbCy02HwDuxN43JxAEaQf//jHga0GSKVSeOedd4Y+P5fLidexMzCurq6iXq+LD8lKpYJarYZcLodyuYzt7W20Wi0Ap687ALz22ms4ODiAJEloNpuwLEt88Obzeezv7wMAFEXB6uoqGo2G7xj6cT4MLMvC/v4+JEkaetbvjLUzsG9vb2N3d1d8INEYOtIvRL34BnbLsiDLsue2zuNeNE2DLMtD5eTjrNVqiSCYSCREesLhBMXO1895DAARyDVNg6IonsdO0pkql8t5nm8UTmqo83elWCxCURSk02l2zRpCz7w60RB8A7tf7ta2bU+AcaTTafHnPnA6K5VluW/q5vDwEEtLS+L49u3buHPnju+gR5kxL7pKpYJMJgPLssQHoSOTyXSdv7a2Bl3XxSy8Xq8DOA3wtm17/ppy7ps35wOq3weDoihd3yt1672scYj0y1YR+PzDmY+PFpdvYJckqesNOOgN2flmzmQy2Nzc7BvYl5eXcXh4OMxYI8FZu23bNlZWVsSsvd1uixUnfh+oKysrUBRFpDQcmUwG7XZ7qJy1M4ZZ2d7eHrgihrn14TCvTuPwXRUjy3JXkHGOe83WAXTl3xOJRM9cfZwkEgmxecexu7sLSZI8+Wjnv50pGbdWq9Xz4mehUIBlWZ6fl6ZpvmPwM+6M2rm4e/PmzbEeTz1sFZ9+EQ3Bd8auKEpXALcsq+9syzRN5PN5HB8fex4X9yVthUIB29vb0DRNvBa5XA7pdFrc5uSedV2HoijY3t4GcPoB6s6xS5KEZDIJ4PRnkU6nUavVkEgkUK/Xsbm5KVI47p9TrzF0clbVABB/PfT6WTu7SZ3zVlZWPDtPnQu+nedubm4in8/zmgvRjC2dnJyc+J3gbExxr9zIZDLi2DRNtNttcdy5MWV1ddX3zcwLacPTdR2WZXley3a7jUql4knLUDQ4m4zul/5u9Me+9QvsfWwjlfk+AJYXiCK/2DlwHbtT58NZ/pZMJj35WV3X0Wg0xG2FQgGqqgIAjo6OOEObItM0u65hjLtyhaJt/fsXxb+dpcF8H8bHwBn7rHHGPhpN02DbtijxYNs2SqVS32seFF6TzNgFV3kBz0VYCr2JZuy0WDjrIqJBWI+diChiOGMnijJuVoolztiJiCKGgZ2IKGIY2ImIIoaBnYgoYkJ78fTa3WYg/9+f/lN3xcVF1263cf36deRyuUAbRxuGgUqlEvg4iKKOM/Y5WF1d9RTkGuf8Xo8f9nkVRUG5XB76/z8rTgMRt1Ffm0leB6K4CO2MPUxG7fvaeX6/RtSjPK8sy6IpxyIZ5XuYxusQV3sf22IHKuvGRB8D+xyMWnu88/xqtSoqOk7yvItolO8hyq/DLLFuTPwwsM9YZ37byTMXi0Uxy2w0GqhUKqJ2fef5vRpR98qb92uOPQxnXOl0WvQjdY/Lr8F2vybawGn9fncpYfdfDb2+B6fMbyaTgSRJkGUZiqKM9DrYti3KE1uWJc4d9NoHwd3+rtO02uGVriRRupL01I2haGNgnzEnv+0ENCfPXK/XRaldy7JQq9VQrVZ7nt+rEXXneUD/5tjDyOVyKBaLaDab4vG5XA4rKyvY39/3bbDdr4k2cNp82/kAAIBm8+lF717fQz6f93SVun79Olqt1kivw2uvveapCe+0Zxz02gfB3f7O4/MPkVo+g/VvPwYwQREwiiVePJ2DXg3A3TNpWZY9XZOGaRje67xWqyXK+PZqjj0M98zVmRl3dmICngZXvybaTjct9/famUpxfw+6rnu6SimKgp2dwa3gOp+jc/ZdLBaxubnp+b7cjx3ndZomp/2d5+vWq7h/69XTmTbRiDhjD8iwwXtUfs2xx5FIJDyz4c4G235NtNvt9kgpDtM0u16XUcsRN5vNns8xzgdn5LBuTGxwxh4yTuu6Tk5z7Js3b6JQKCCdTnvuG5fTXKUfdxs+9xdwOuMeZTasKMrQH0b9XodkMtn1HE79+oXH3qY0JQzsc9DZZHrU8/s1onafN6g59rBjcAdi27ZhmqbvCgq/JtpOgHff12q1PMfuxzrnu8cwTENu93OUSqWu5unb29tsru3iLH3MZrNc/x9RTMXMWLvdRr1eh2maMAwDsix3HddqNc/KEvf9uVyuZyPqzuf1a46dy+W6zvXjzIabzaa4yGgYRt8G235NtHd2dkQTa+exuq73/V6d84dtyN3rdXBWumQyGZimiXK5DEVRev4s3K+9u+VjVHHpYzywNR4Jqqri6OiI2/3nyNO2bp4pGLbMCz2/2MlUDBFRxDCwE4CnqRZd1z0rXIgofJhjJwCnuWz3ph4iCi8GdqI44pr2SGMqhogoYhjYiYgihoGdiChiGNiJ5kzTNLHz09kkRDRNDOxEc+aU6gXwtDQv68PQFA21KkZVVdG0ABhtC3K5XB66JvhIgnojrG+PdLq7ucPGxsaMBhUcNqgej1OqN+iAzpZ50TRwxu50lykUCiiVStjf3+9bWa/XY4OudR00p4HFLCxCASc2qA6v9e9fROp5Cfj8Q+w1P8DWO7eCHhJNycDArmmapzhSsVgcagbeWWGPpqtfY+dFUC6Xh+5D6tegmr1MZ6t0JSkaeqSel4IeDk2Rb2DvFZwlSRpqy/nu7q6o6EfTt8hpj1wuN3T9837fxyjPQURevjl2y7K6us0M031G13Wsra2xaqOLk8JyenkWCgURuPyaQTslam3bRqPRQK1W69vYuRMbVPODgeLJN7D7/alv23bPtmXO7cO2NDs8PMTS0pI4vn37Nu7cuTPUY8PEsixPk+hkMikCYr9m0KqqIpfLiZ6izvWKfo2dO7FB9eL+VUM0S76pGEmSutqMDWpddu/evZFyo8vLyzg5ORFfUQzqALpmj06TaL9m0IlEAtevX4emabBte+wVC2xQTRQvvjN2WZa7Zu3Oca8Zebvd5gWvIbmbRPdrBu3Msmu1GsrlMkql0lSWjrJBNVG0+QZ2RVG63oSWZfUN3pZleYJBs9mEaZpQVdWTU6bT12plZQWyLPf9QDQMA4VCQQT4lZUVmKbZ9TqO2tbNNE2srq72vT+TyfQdk6IookXeMEY5v9/3kUwmRZrFEZoG1UQBGLjccW1tzbNuvdFoeNYtOxfZgNP87cbGhvjK5/OQJAkbGxuxfxN2pgacJtF+zaAbjYbnce5A26+x86D/NxtUE0XfwJ2ntVoNqqqKFQzJZNIzq9J1HY1Go2umpWmaaBysqipKpdLIf4JHhbOCw90k2kltAP2bQSeTSRiGIVbFZDIZ8QHZq7GzHzaopkG4CzU62Mw64tigevFkX/5bAMD9W68GPJKntP+zj60PPgH+9mXs7e09LXlAC8svdrKDEhGhdCWJ0pUksL4tZu0UXgzsEeakWmzbRj6f54qlAGmahq2tLQCnKY+F3cLPlnmRwMAeYWxQvTicUr2pVAqp5yWsf/9i0EOiCGNgJ5qTRSnVS9HHRhtERBHDwE5EFDEM7EREEcPATkQUMQzsREQRw8BORBQxDOxEM6JpGrLZLLLZLPb29oIeDsUIAzvRjDibkoDTNezr6+vBDohigxuUiGaIxbQoCAzsRPPA3aY0R0zFEFEXpzZ7Npv1NEahcOCMnYg83AXKnGsEbLoRLgzsROTB2uzhx1QMEVHEMLATEUUMAzsR9eZ0U/r8Q67qCRkGdiKiiGFgJyKKGAZ2IqKIYWAnIooYrmOnmbl2tyn+/dN/yox9DhGNhoGd5oIBnGh+GNhpYgzaRIuFgZ3mzv1BQETTN1RgV1UViUQClmUB8C8IZNs27t27BwDY398HAFSr1UnHSQuGwTk+9j62kX3rF4CWxfr6OguChcDAwF6pVJDJZFAoFMSxruviuNf51WoVkiQBAFZWVqCqKjY2NqY3aqIFpWkatra2AAB7zQ+Qel4K9a5NVnoMp4HLHTVN8wTxYrGIWq3W9/zd3V0YhiGOE4kEmk3O7miwa3eb4iusPO3wnpc8gTGMSleSuH/rVdy/9SpSy2dYXiAkfGfs7Xa76zZJkjyBu1Or1ep6jkqlMubwKGzCHJSnRbTDYwCkgPgGdsuyIMuy57bOYz+qqiKXy/FPNxoZV9oQjc83FWPb9tj3Oe20ksmk7wAODw+xtLQkvu7cueN7PhER+fOdsUuSJFbCODqP+z3OmaXn83k0m03U6/We5y4vL+Pw8HDY8dKcMbVCFD6+gV2W5a6ZuXPsrHrpvE/TNM8KmHw+zxx7yDCYE4WbbypGUZSuAG5ZFnK5XM/zd3d3UalUfNM0REQ0WwOXO66trUHXdXHcaDRQLpfFsWma4v5cLudZw+6czzXsNIkoLIMkmqeBG5RqtRpUVYVhGDBNE8lk0rOuXdd1NBoNcVuhUICqqgCAo6Mj5PN5BnYiojkaqqSAX2De2Njw3J9IJBjIQ+539/8d5q/+QxwnXnkd38n+sO/tRLRYWASMupi/+g9Yn/we8sVLsD75PQDgO9kf9r19ngGf69uJBmMHJepJvngJP/jJu5AvXup7u/XJ7/Hzt9/EB/+2iT8+ON2l/Pjj3+KrX76LG5/dCmLYgdA0DdlsFtlsVpQTIAoSZ+w0lsQrr4t/n3tJEbP03/zL6z6PiianPkwqlUJq+QzWv/2Y5QQoUAzsBMBbldBJtwDAhScHAIAbn93Cb1z//tfsW8yvu7A+DC0SpmIIgLcqoXzxkmdG/uDhI7xRa+LBw0cBjY5Gsff/bPFF8cQZOwnOrNN9gfJq6rz490vnnxXH/XLob8x2iBQwNt0IBwZ28vWj730LP/ret4IeBi0ANt0IDwb2GPN0+/n64t80OKmbT8+8ybXuczJq2sV9fuq/SUM9pnQlidKV02qt2bd+8bTpxvr2SP9vmj3m2GPM0+0nlcL6+vrEz3k1dR4vnX8WwOlFWPf6diKaD87YY06s5nCZZA26O3Vz9X99McnQBor7ZiVeHKV+GNhjJO6BMI7GSblQ+DGwx4yz/f/820/Xq3dWTbwR0NgmMe8PrVldnyCaBubYY8ap9wJ0r1en4c3i+sSscX17fHDGHkNOvZdZc+9a/de/e2vm/7958+423QG2doIe0liYrokeBnYCMNkFUyJaLAzsFDlRu0jMGTWNioGdZsrZrAS8jqup8/jR974VybTMuBYtaC/aeGg8DOwxNY/Ui7vOjFNAjOUJooN1YxYXAzvNjHuz0umsnYbVb+XKPFe0dP6/3DN41o1ZbAzsFGlRy7cvCtaNWWwM7EQLguvLaVq4QYmIKGI4Y4+4znIBFDzOzGnWGNhjwKkPc+HJAR4/fPR1Wd3nAhuPe0UOlz4STR9TMTHgrg/jbm9HRNHEGXtMyBcv4b1/Dm6WDrg3K0FsVpqnSVfIxL2iIzcvhQdn7BGlaRqy2Sx+/vabYrYeJHdnpQcPH+H9vYcBj2h0YazoSPHEGXtEOUHo2/Jf4cK5Z3D15S8DHU9UNiv16jg1DF4wpXliYI+wVCqFd67+16CH4SuoC6kTb1zaKk5xNOHDD6rFNlRgV1UViUQClmUB8N86bNs2NE0DADSbTeTzeW41poU2aZBn7pkWzcDAXqlUkMlkUCgUxLGu6+K40+bmJqrVqjhOJk+3HTO4U9g9+ONpIbNrd5v46ZmAB0PkY+DFU03TPEG8WCyiVqv1PNe2bZim6bmtXC57Aj1RFAzTZi5OreicSo/ZbFb8xU7B8Z2xt9vtrtskSYJhGH0fYxgGTNNEIpEQ53cGe6KoikMQ7/Tfv/M3ePzlX/D4y7/gP5sfAJ9/iNJf77AgWIB8A7tlWZBl2XNb57GbJEk4Pj723NZoNJDL5SYYIkVR0Gva+2EJhtG5Vzz9+N/+b8CjIWBAYLdt2/c+SZJ8n9y2bRiGgZ2d/k1+Dw8PsbS0JI5v376NO3fu+D4vdQtTeVr3ztfWwTFaB8d4f+8hPj3zJgAg8crr+E72h0ENjybw+Mu/ADj9yyUV7FBizTewS5IkVsI4Oo/9XL9+HfV6HYqi9D1neXkZh4eHQz8nDXbtblNc6AMWb7mje4b3s1//wbNZydlMxcBOND7fwC7Lctes3TkeNFtXVRXlcplpGPLlDvIA8EbtAHhyEOCIiMLPd1WMoihdAdyyrIHBWtd1KIoizvO72EpERNM1cLnj2toadF0Xx41GA+VyWRybpum53zAMWJaFdDotlj/2Wl1DRNF27W5TfNF8DQzstVoNpmnCMAxomoZkMulZ167ruljXbts28vk8yuUyzp49i7NnzyKZTKLZ5A+WiGhehiopsLGx4Xufc78kSTg5OZnOyGhki9ZQg4iCwbK9EcKGGkQEsLpjqLlzlzc+u4XfPDnAhXPP4L3yYq9jJ6LZYmAPOSf98psnB3gg0i80Lc7rC5yusZcvXgp4RESDMRUTcky/zJb79f3uuWfwxstfemrIEy0iztgjYBH6mU5bUA04epEvXsIPfvIuAzqFBmfsREQRwxk7Ec1cmIrURQEDe8jEYRdfZ0lf/GPAA6KROD+/T8+82bNSJ4P87DGw00JxX/x98PC0QuU/MLCHhvvnx0qdwWFgp4XirvbozNopPLp+fk8OcOOzW4FfAI8bXjwlIooYzthp4c176SM3JVHYMbCHEAPPbDmbkuSLlyBfvITEK68HPaTI4oXU2WBgD6HegacV9LAixdmURBRGDOwhJV+8hPfFbtPoBvV5L310iqk5/yYKIwZ2Wlhc+kg0Hgb2EIjDpqReei19nNWFVFbJpCjhcseQufHZLVx4coALX6cLaDpYJZOihDN2oq9FsUpm0AaVF3DjCpnpYWAPCaYKKGzcf/U8/vi3+OrJAW78fYu7UOeAgT0knFTBhXPPxDZVwOJg4cLyEMFhYA+ROKcKZrlCxr3EEYjn60vRwsBOocDZX7ww3z4ZBvYFFdcljqNYpPZ5RIuEgT0EmCogolEwsC8QztKHxwupRP0xsFPosNQAkT8G9gXGteu9dV5IffDwEfAvp6V1r6bOi/uA/rn3G5/dws9+/Qe8v/cQbwB8fedglM1KNJmhAruqqkgkErAsCwBQKpV8z7dtG5qm4ejoCNVqdfJRxhTXrg/Wa/buDuydFRrdwbx1cAwAWHnxLF/fGZtksxJXyIxuYGCvVCrIZDIoFAriWNd1cdzJMAzYto39/f3pjjSm4rx2fRi9Zu+9lkM6M/n39x6K2fnKi2e7Zvg0G1yuOl8DA7umaZ5Zd7FYRKVS6RvYc7kcAKDZbMK27emMkmgI/WbcnTP5l84/i/fKnPlRdPkG9na73XWbJEkwDGNmAyIal3tW6MYZIsWNb9ley7Igy7Lnts5jIiJaLL6B3S+VMq00y+HhIZaWlsTXnTt3pvK8RG5O7t1Jy1A4XbvbFF/Un28qRpIksRLG0Xk8qeXlZRweHk71ORed3y+ls8QRgGhYTZNx5965+oXiwDewy7LcNTN3jiVJmtGQ4s1Z4vjdc8/gwrlncPXlL4MeUuj1y70TRZVvYFcUpSuAW5YlVr7QbHCJIxFNYmDP07W1Nei6Lo4bjQbK5bI4Nk3Tcz8R0SDONY+fv/0mfnf/34MeTuQMDOy1Wg2macIwDGiahmQy6VnDrus6arWaOG6321BVFbquwzAMqKrac9kkEcXT1dR5Ub7h8ce/xVe/fLdrhzBNZunk5OQkyAGk02ns7u4GOYS567x42uuC6ftMxVAMOHsM3itnRq6pH/fyAn6xk0XAFgAvmBKNjjVk+mNgXxC8YEpE08LAPifcUEFE88LAPkMM5kSDsU779DGwB4Q7TIkmq9NO/TGwB4QXTIlYp31WGNgDxAumRNPBFTJeAzcoERFRuHDGTkSRwtk7AzsRRVhcgzwD+5T1W+LoXgUDcCUMUS9c+jgdzLHPibMKBgAuPDnAd889gze4EoZIYHGw6eGMfY7ki5fwg5/wl5WoFy59nB7O2ImIIoaBnYgoYpiKGdOodWAuPDkAAKZhiEbkfs+w1MBwGNhnyL0S5vHDR+LCEBEN5qyQcfv0zJsAMNaKmX6TsSgug2QqZobcK2FeOv+sp+AREfXnXiHTyfrk956lw9SNM/YZ40oYotG5V8h0eqN2AHyd2qTeGNhniHl1IgoCA/uUMa9OFC5RLDvAwD5lTl5dvniJeXWikIlKkGdgH8GwSxyZVyeiIDGwDzBMMPekXz7+LV46/yyDOtEMPXj4CJ++Pf7Sx6hjYJ+Cr3757tN8OtMvRDMl3l9PDvDg4SN89eQAmEFgD3NahoG9h1F3lQKn69TfK4frh08URkEXC+uMD4sY9BnYvzZqMOfqF6LFMI+0zDiTvSDFOrBPEsz/+KANADj3ksLVL0QBmVdaJmyGCuyqqiKRSMCyLABAqVSa6vlh4c6lX3jxLK6mzuNH33sOwHNBD40oljrTMkFcVF3EXPzAwF6pVJDJZFAoFMSxruvieNLz520aKRfm0okWj3v23jo4xh8ftPHVL98V9x3+493gBjdnSycnJyd+J5w9exbHx8fiuN1uo1KpoNFoTOX8dDqN3d3dccbe16T5sF4pl5UXzwLA17P03jUs/Lzb+E+8mf/2ROOi6eLPZDFN4+fys1//Ae/vPQQAtA5O49G5lxQA85vJz3r27hc7fWfs7Xa76zZJkmAYxlTOD1png2mHO5g/TbmMHszdajsmg8iC4c9kMU3j5+JO0Ygg32MmDwCfnnmx6/HTCP7DTDBnFfx9A7tlWZBl2XNb5/Ek50/q2t0mfrX1P0Vp3F4u+FSBcz7Jndm4eMyUgjkRBa9nkHfpjBFO8B+3NHCvD4XOSaR88RJeWf8fYz3/MHwDu23bvvdJkjTR+ZO68dktfPmn3+HBk0djPX6FAZwoVvzKATvcM/xRuf8icP8l4F5FNw++OXbDMLC6uurJmZumiWQyiePj465APer5APCNb3wDf/7zn8Xx+fPnsby8PMG3tJgODw8j+X2FGX8mi4k/l+F89NFH+OKLL3re5ztjl2W5axbuHPcK0qOeDwB/+tOf/IZAREQj8m2NpyhKV0C2LAu5XG4q5xMR0fQN7Hm6trYGXdfFcaPRQLlcFsemaXruH3Q+ERHN1sB17MDpTlJFUWCaJgDvTlJVVdFoNDzr1P3OJyKi2RoqsNP02LYNTdMAAM1mE/l8nh98cxbVkhdhxvfFdMW6CFgQNjc3Ua1WxXEymQTA4DIvi17yIq74vpiugTl2mh7btkV6ylEulz2/0DRbmqZ5gnixWEStVgtwRMT3xfQxsM+ZYRieX2JJkrp+qWk2wlbyIk74vpgupmLmSJIkz+Yt4HTVEJeDzse8S17QcPi+mD4G9gDZtg3DMLCzsxP0UGJh3iUvaDx8X0yOqZgAXb9+HfV6HYoyn/oRcSdJklgJ4+g8puDxfTE5ztgnpOs6tre3fc+RZbnrAp2qqiiXy/xzc47GKXlB88X3xXRwHXsAdF2HJEnil9cwDP4iz0lnIxjDMFCtVvs2gqH54ftiepiKmTPDMGBZFtLptFjm1Wu1Bs0GS14sJr4vposz9jmybRtnz57tur1QKKBerwcwonhiyYvFwvfF9DGwExFFDFMxREQRw8BORBQxDOxERBHDwE5EFDEM7EREEcPATkQUMQzsREQRw8BORBQxDOxERBHz/wHLgKcFilC0NgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "target_hist, target_bins = np.histogram(TARGET_DATA, bins=100, density=True)\n",
    "target_mids = 0.5*(target_bins[:-1] + target_bins[1:])\n",
    "__=ax.plot(target_mids, target_hist, drawstyle='steps', color='k', label=r'${\\rm target\\ PDF}$')\n",
    "\n",
    "__=ax.hist(pred_init, bins=target_bins, density=True, alpha=0.7, label=r'${\\rm initial\\ prediction}$')\n",
    "__=ax.hist(pred_best, bins=target_bins, density=True, alpha=0.7, label=r'${\\rm best\\ prediction}$')\n",
    "\n",
    "leg =ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25284a7a",
   "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

kde_experiment.py

"""
"""

from functools import partial

from jax import jit as jjit
from jax import numpy as jnp
from jax import random as jran
from jax.scipy.stats import gaussian_kde


@jjit
def _mse(pred, target):
    diff = pred - target
    return jnp.mean(diff * diff)


@partial(jjit, static_argnames=["npts"])
def mc_double_gaussian(params, ran_key, npts):
    mu1, mu2, sig1, sig2, frac1 = params
    pop1_key, pop2_key = jran.split(ran_key, 2)
    pop1 = jran.normal(pop1_key, shape=(npts,)) * sig1 + mu1
    pop2 = jran.normal(pop2_key, shape=(npts,)) * sig2 + mu2
    return pop1, pop2, frac1


@partial(jjit, static_argnames=["npts_pred"])
def predict_pdf(params, ran_key, pdf_abscissa, npts_pred):
    """Toy model prediction for the pdf at the input abscissa"""
    pop1, pop2, frac1 = mc_double_gaussian(params, ran_key, npts_pred)
    pred_kde_pop1 = gaussian_kde(pop1.T)
    pred_kde_pop2 = gaussian_kde(pop2.T)
    pred_pdf_pop1 = pred_kde_pop1.pdf(pdf_abscissa)
    pred_pdf_pop2 = pred_kde_pop2.pdf(pdf_abscissa)
    return frac1 * pred_pdf_pop1 + (1 - frac1) * pred_pdf_pop2


def build_loss_func(ran_key, target_data, ncells):
    """

    Parameters
    ----------
    ran_key : jax.random.PRNGKey

    target_data : ndarray of shape (ndim, npts)

    ncells : int
        Number of points at which to evaluate the target PDF

    Returns
    -------
    loss_func accepts (params, ran_key) returns MSE loss

    """
    npts_data = len(target_data)
    # precompute kde model of data
    # returned loss_kern will treat this kde model as fixed data
    data_kde = gaussian_kde(target_data)

    # randomly select points at which to predict pdf
    # do this by drawing samples from a KDE model of the data
    # this eliminates the need to hand-tune bin edges
    target_abscissa = data_kde.resample(ran_key, (ncells,))

    # evaluate pdf of the target data at the randomly generated abscissa
    target_pdf = data_kde.pdf(target_abscissa)

    @jjit
    def loss_func(params, pred_key):
        pred_pdf = predict_pdf(params, pred_key, target_abscissa, npts_data)
        return _mse(pred_pdf, target_pdf)

    return loss_func