Demonstration of fitting a toy model for subhalo mass loss to itself

demo_toy_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "annual-helmet",
   "metadata": {},
   "source": [
    "# Demo notebook for toy model of subhalo mass loss\n",
    "\n",
    "Executing this notebook requires the toy_model.py module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "excellent-niagara",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "hydraulic-development",
   "metadata": {},
   "outputs": [],
   "source": [
    "from toy_model import get_adam_opt_funcs\n",
    "\n",
    "opt_state, opt_update, get_params = get_adam_opt_funcs()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "boring-fitness",
   "metadata": {},
   "source": [
    "## Define our target data\n",
    "\n",
    "First choose some random point in parameter space for both the target data and also the initial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "lesser-standard",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Target parameters = [3.07721998 1.72745557 0.08644907]\n"
     ]
    }
   ],
   "source": [
    "from toy_model import get_default_params\n",
    "p_default = get_default_params()\n",
    "\n",
    "rng = np.random.RandomState(seed=42)\n",
    "p_init = rng.normal(loc=p_default, scale=0.3)\n",
    "\n",
    "rng = np.random.RandomState(seed=43)\n",
    "p_target = rng.normal(loc=p_default, scale=0.3)\n",
    "\n",
    "print(\"Target parameters = {}\".format(p_target))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "neural-moses",
   "metadata": {},
   "source": [
    "Now define the time array and the target data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "comparable-collapse",
   "metadata": {},
   "outputs": [],
   "source": [
    "from toy_model import predict_frac_mass_loss\n",
    "tarr = np.linspace(0, 15, 50)\n",
    "mass_loss_target = predict_frac_mass_loss(p_target, tarr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "through-washington",
   "metadata": {},
   "source": [
    "### Calculate the predictions of the initial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "turned-cliff",
   "metadata": {},
   "outputs": [],
   "source": [
    "frac_mass_loss_init = predict_frac_mass_loss(p_init, tarr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "executed-algorithm",
   "metadata": {},
   "source": [
    "### Compare predictions of initial guess to target"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "rubber-optics",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEOCAYAAABIESrBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0nUlEQVR4nO3de3zU9Z3v8dc3d8JtSIBAFIQJiqAVDAHvl0pivZTaUwP0Yru1WxK72z3bnloiu9uW7XaPBttztrdtA7W228sREu1NXTWh1Vq1khDxhlcG8RIESRiQO0m+54/fb8JkMjOZ/JLM5PJ+Ph7zyMzv+/395gPJ5JPv73sz1lpERES8SEt1ACIiMnwpiYiIiGdKIiIi4pmSiIiIeKYkIiIinmWkOoCBNnnyZDtr1qxUhyEiMqxs3bp1n7V2Sl/PG3FJZNasWTQ1NaU6DBGRYcUYs8vLebqdJSIinimJiIiIZ0oiIiLimZKIiIh4piQiIiKeKYmIiIhnSiIiIuKZkoiIiHimJCIiIp4piYiIiGdKIiIi4pmSiIiIeKYkIiIiniV9FV9jjA+oAPKttVUJ1F8NBIA8AGvt+kENUEREEpbUlogxphQoBYoAXwL1q4GAtbbOTR5FxpjywY1SREQSldQkYq1tsNbWAcEET6lw64dsBCoHPDARGRbWr0/9jYihEMNQMmQ3pTLGFEc5HMRpyfRu11Nw+D3IyHYe6dmQPR6mneuUHzsA6VmQkQPGDFTYIjJIAoEAwWBw1Mcw1AzZJILTB9IWcSzydQ8tLS2MPftSfvvhw5SNC3QrOzxuFq+seJRpE3KY/ptPYHY9AWkZkDXOSTCF58PKXziV//pjaD8K4wpOPfL8kJU7QP88kdT61z+8yPaWgyl57/mFE/jGsnP6dE51dTVFRUWDFBEEg0F8Pl9KYxiOhnIS8cUqMMb4rLXBaGWFhYXU3LWJuse3cM+hNo4cPUrHiaNkmXZOtmbw9H8+CcANmYs5L7eI6WPamZp1grzME2Slz2D80ZNMHJMJz/wC9rzQ/eKzL4e/+YPz/KkfwtgpMOVsmHwWZOYMyD9aRHpqaGggEAgQCATw+Xz4/X5KS52bEs3NzYDTSqivr6e6uhqfz0dDQwOVlZVUVlbi8/moqalh8+bN+Hw+mpubaWhowO/3U19fT1lZGY2NjVRXVxMIBKirq8Pv99PY2MiaNWu6rhcrhtFsKCeRIO6IrDCRr6NatqCQZQs+2vX62MkO2g6foPXQCfYdPk5L8Chv7JvNU/uO8EbrYd585wgnOjoBMFsf4expE7hg9nouvCCbxVPaybdBONgC2ROcC3a0wx+/BSePOK/TMmHWJbCkAs6+vl//aJFk6WtLIJVKS0u7kkVFRUW3suXLl1NbW0t5uTPmpqqqipqaGkpLS6msrGTjxo1s3boVoKulsXTpUnbu3InP56OxsZG2tjaqq6sBKCsrY8eOHQAUFxezfPly6uvr48Ywmg3lJNJGz9aIDyBWKySWnMx0Cn1jKPSNiVre0WlpCR5l577DbHsryJadbWxsfIufPdkBgH/yWJbMnsP1503nkk5LWnoGVO2C1tfhvZeg5Rl49WEIvuVc8PA+ePL7MPdaOH0xpKX3JVwR6YOtW7d2JQe/308g0P02tt/vB3r+4g+/dRVKGuvXr6e4uLjbuU1NTYMQ9cgxZJOItbbZGBOMOJwHNAz0e6WnGWbk5TIjL5fLz5oCwMmOTl5sOciWna1s2dnGA8/v5p7Gt5iZl8vHl8ygfNHpTC2YDwXz4dwb4epvQaeTdNi9DZ76ATzxHzB+Olz2FSj+jNPBLyIDrqqqisWLF9PW1kZbW/eu08WLF/eov2LFCurq6rpaF7W1tYCTTILBIA0Np37NhMokuiGVRIwxfqA4bFjvJmNMedjrMqAmGbFkpqexcIaPhTN8VFxexLGTHTz84rv8vy1vsu6hV/g/j7xK2fwCPnnBTC4pmkxamjnV4phTCqsD8Ppm2LIBHrwVnvge3PI4jPElI3yREa2uro7y8nKCwSCLFi3qao00NzdTU+P8iog3imrRokUUFxfT1NREfX191/HFixfT3NycUF9HKIbRLtmTDYvdGejlQKkxZnXEUN5ywuaBWGsrAb8xptQYUwHsiJg3kjQ5mencsPA07qm4iM1fuYKbL5nFXwOtfPquLXzwO4+y+aU9ESdMhHM/Bjc/CDfd5zwPJZB3tkJnZ9L/DSLDmd/vp7W1tduxpqYmfD5f162p0K2sUAd4LFu3bo3aMV5eXk5bW1u3BBQ+LyRaDKOdsdamOoYBVVJSYpN1D/N4ewcPvfAuP/zT67y65xAfWVDIN5bNJ39cnNtWbTvh+4tg6jz40L+D/8qkxCoyEixfvpyysrJuCaCyspJFixbh9/vx+/1dw3CLi4upqqrqqhPeJ1JVVUVdnfP3aFtbGyUlJdTU1HT1qdTU1HTdBisuLu7qV4kVw0hgjNlqrS3p83lKIv13or2T/3z0dX74p9cZn5PJN5bN5yMLCjHRJjF2dsAL98Gf/h2Cu+DadbBkVVLjFRnN6urqaGtr65ZUmpubqaqq6nZra7TxmkS0iu8AyMpI40ulZ3H/P1zGjLxc/vGebXz+503sPnC0Z+W0dDhvOXzhCTjrGqe/5OF/Tn7QIqNUIBDo1rIAuo3Ikr5RS2SAdXRa7n5iJ99+5BUy09L42ofns2LxjOiVOzucBDKhEC75n8kNVGQUW79+PcFgEL/f39UHUlFR0euM9ZFMt7NcqU4iIW+2HqHq3ud4KtDKmmvPpvKKOEslWOus3/XWFpg0C8ZNTVqcIiKg21lDzsz8XH7xt0v48HnTuf2/X+anf9kZu7IxcPIYbPoM/GQp7H05eYGKiPSDksggykhP4/+uXMg150zjm/dv55d/3RW7cmYOfPzXTjK562pnFWIRkSFOSWSQZaan8b1PnE/pvKn8y29fYFPjW7Ern1YMqzbD2MlOq+TQ3uQFKiLigZJIEmRlpPHDTxVzxVlTqLrvOe5rfjt2Zd9MWPlLOH7Qme0uIjKEKYkkSXZGOjWfXsTFRfncWvssf3i2JXblgvnw+Qa4ck3yAhQR8UBJJIlyMtPZ8JkSSmbl8aWN23j4xXdjV572AUhLgwPvQMu2pMUoItIXSiJJlpuVwU8/u5gPnDaRr2x6NvqExBBrnb6Rez4FR3rd1FFEJOmURFJgXHYG3/v4+XR0Wr722xeIOVfHGLjuTji0B373RSepiEhKNDc3s2jRoq71uMShJJIiM/Nz+crVZ9Hw0l7uf2537IqnFUPZN+GVB+DppKyCLzIkLV++vNuKul7qRzs/0esWFxdTWVnZa73RRkkkhT578SzOO30ia3//IvsPn4hd8cIvwNzroP5rsPvZ5AUoMoRUVlb2adXcyPqBQCDqHiN9uW5eXkI7dI8qSiIplJGeRvWN53Hg6Em+9cBLsSsaAzf8EBbdDJNmJy9AkSGktLS0x8KJfakf2kO9v9f1It4GWcNdn3c2NMasAvzARmvtNmPMj9zXAaDaWvvGwIY4ss2bPoFbrijiB396nRsWFnZtz9tDbh5cty65wcnId/f1PY+d81Fne4ITR+BXy3uWL/wknP8pONzqDPyItPhzzpbRB96G+6Lc/rn4izD32j6F2dzczKpVqygtLaW6upqGhgaqqqpYuXJlVwKor6+nqqoKv98ftX5ooyqfz9e1F0hkvdB7gdNyqa+vp7q6uk8LMzY3N9PQ0IDf76e+vp6ysjIaGxspKyujqqqqa++S0PLzfr+/azfGQCBAXV0dfr+fxsZG1qxZ0/Xe69evx+/3EwwGqa+v7zon1vFk8bI9bhuwyVp7IJRQrLUfAjDGfAx4YwDjGxW+eNUcHnxhN//0m+d5+EuXMzY7zrfl3efhoTVw410wviB5QYqkUKg/YseOHYDTeqisrKS2trZrD5C2tjZqamqorq6OWj+UHML3EYmsB04fSW1tbdfWt1VVVX36xbx06VJ27tyJz+ejsbGRtra2rgRVWVnJ1q1bu733xo0bu84tKyvriqW4uJjly5dTX1/PunXrKC0t7VqyPrRrY6zjyeQliQSttQfc5+VAeBvxQJT60ouczHTu+Nh5rKh5iu888ipfXzY/duXMXNj1JDz+HbVMpP9ufiB2WVZu/PKx+fHLJ54ev7yP8vLyuv2yB7q1EPLy8rptKhWtfiLXDe3XDnTtdNhX4XElEgM4LYrwfU38fj+hFcn9fj+rVq2isrKSFStWdCXCWMeTyUufSPg40zKgKUaZ9MGS2XncdOFM7n5yJ8+8uT92xfwiOP8m2Ho3BN9MXoAiQ9BgdXSHts9tamqira1vc7RWrFhBXV0dwWCQ5uZm1qxJbOWJHTt2EAwGaWho6HrU1tYCzt7va9asoba2lkmTJnUNM451PJm8JJEiY8wsY8wdQJ219iCAMebzAxva6FN1zdlMm5DDbfc+z4n2ztgVr1jtfH0sekehiMQX2l89UjAYZNGiRaxZs4by8nJKSkq6lSVi0aJFFBcX09TURH19fdz+lPAEFdrTvbS0tNsDoKGhgfLycurr67HW0tTURCAQiHk8mfqcRKy1G3BaIK9ba1cYYyYaY24HigDtMdkP43My+dZHz+WVPe9zV7z9RyaeDiV/C9t+DfteS16AIikU2oHQa32/309ra2vcek1NTfh8vq5f/KFfyKFO+URi2Lp1a1fHfaTQTooh4bffysvLe1w/NH+lvr6+W3IIXTvW8WTyOjqrCLjHPVSNMzprB6DZcP20dF4BV5w1hbv+EuDmS2aRk5keveJlX4FJZ8DEGFvviowgzc3N1NbWdv31nZeX1+N1TU1Nt9FN4eWlpaWUl5ezcePGrtFM0a5bWlpKSUlJV53i4mJKSkqoq6ujtLS0R91ofD4fRUXOTqZtbW1do7FCiaW2trarJRQasbV+/XoqKiqora3l9ttv79YqASgqKuoa8RUMBlm8eDF+vz/m8WTq8/a4xpgbgYaw0Vnl4aOzrLX3DUKcCRsq2+P2x5Ov7+OTP3ma2z/2AT6xZGaqwxGRBNXV1dHW1tatgzs0lDe81TEUJXN7XI3OGmQXFeVzTuEENjweoLOzlyT/fB3c/+XkBCYicQUCgR4tgfARVyORlyG+kaOzlscoE4+MMVRc7ucf79nGH1/eS+n8OPNB9r8BTT+FBZ+EGYuTFqOI9LR69WrWr1/PunXruvo/gsFg1yirkcjL7axVQD1wC85EwxXu8c8DAWvtHxO4xmqcGe55ANbauKufufWD7kuftTbmBImRcDsL4GRHJ1fe+Sin+caw6ZaLYlc8fgi+txCmzoO/+UPS4hORkSVpt7OijM6a4A73nUMCo7OMMdU4yabOTR5FxpjyOPVXW2vXWWvXu/Ub3KQyomWmp/G5S2ez5Y22+PNGssfBZbfCzj9D4NGkxSciAh4XYHQTiXHXzdqAk1Bus9Z+O4HTK6y14YO0NwLx1ldeGfHezcCouG+zcvEMxudk8JPH4wz3BSi5GSacDpv/TXuOiEhSeUoixphHcFodAZwZ6yXGmEZjzIRezovWUgkC8QY3txljum4oGmMqcBLPiDcuO4NPXXAG//3Cbt5sPRK7YkY2XP8duOR/KomISFJ5mSfyeWB52Ait0HEfUAHEa43k4SzgGK63NQUqga3GmP3A7bi3wvoU9DB28yWzuOsvAe76S4B/veHc2BXnXpO8oEREXF5aIvsjEwiAtTYI9HLfBV+sAjcJ9WCtDeAkjyac4cRxb2W1tLRgjOl6rF27tpeQhraCCTncsPA0NjW9HX/jKnD2YX+02lmCW0QkCfq7AGNfysC5dRW5YlrcFdSMMTU4kxvLcDr0K8Jvb0UqLCzEWtv1GO5JBGDVZX6OnuzgV0/vil/x+Pvw2B3OkF8RkSTwkkTyjTGzIg8aYxbiLIcSTxs9WyM+6GrJRF6zGGdyY7NbpwGYTfw+lBFn7rTxXDl3Cj97chfHTnbErjjpDDjrGtj6c2g/nrwARWTU8jrEd53bkb7RfTQCa6y1d/ZybjOn5nuE5AENMU7JA7qtmOYmm1j1R6yKy/zsO3Sc3z7zTvyKS1bBkX3w4m+TEpeIjG5eh/iuwOlEb8Dpq6iw1q6Mf1aXTRHzQsoIW7jRGOMPlbstj7Lwk92+k+Rv35ViCS+FMvtKyD8TtsSdvykiMiA8JREAa+0z1toN1to7rbXPABhjbk3gvErAb4wpdYfr7ogYbVVO93kjlcaYamNMhVt/hbU2+TuvpFhoKZQd7x3msVffi10xLQ0uqIRxBc4e2SIigyjusifGmIf7ci1gkbU2v99R9cNIWfYkmpMdnSz59wYuO3MK3/vE+akOR0RGEK/LnvQ2T8QAVfTsx4hV946+BiCJy0xP47oPTOe+5nc4fLydsdm9fPtad0BuHoyZlJwARWTU6e12VpV722pnAo8ATsKRQXTDwtM4erKDhpf2xK8YfBO+v8gZqSUiMkjiJpFQX0eirLW9TTaUfio5YxKFE3P43baW+BV9M2HWpdB4F3TGGRYsItIPnjvWJTXS0gzLFhby51ffo623GexLVsGBN+HVvnRtiYgkTklkGLphwWm0d1oefH53/Ipzr4fxhdC4ITmBiciooyQyDM2bPp4zp47j973d0krPgJLPwa4n4VCcYcEiIh4piQxDxhhuWFjIljfaeCd4NH7lJavgyy/CuCnJCU5ERhUlkWHqIwtOA+D+Z3tpjYzxwdjJznPtNSIiA0xJZJiamZ/L+TN9vY/SAmeJ+Luuhm2/GvzARGRU6XMSMcasMsbc7q7aizHmR8aYh92vswY6QInthgWFbN99kNf2vB+/4phJcDQITXcnJS4RGT28tETagDustduMMasAv7X2Q9baL+BsmStJcv15haQZ+H1vt7SMgUWfhXea4N3nkxKbiIwOXpJIMGxnw3Kc3QZDeux4KINnyvhsLpkzmd9tayHeGmgALPg4pGfD1p8lJTYRGR36u7NhGc5S8NHKJAk+sqCQN9uOsO2tYPyKuXlwzkfhuU1w4nAyQhORUcBLEikyxswyxtwB1FlrDwIYYz4/sKFJIj507jSyMtIS62C/8Auw9Os4a2WKiPSf150Ny4DXrbUrjDETjTG342yNqz6RJJuQk8nSs6dy/3O7ae/ojF+58Hxn3khWbnKCE5ERz9PoLJyEEbqNVQ0swtkrvS7GaTKIblhYyL5Dx3kq0Np75ROHnUUZ9748+IGJyIjndXTW7WGjs2Zba6/W6KzUuXLuVMZnZyR2S6v9ODy0Bhp/MviBiciIp9FZI0BOZjrXnDuNh154l2Mne1n2vauDfaM62EWk3zQ6a4T48IJCDh1v58kd+3qvvOizcPwgvPibQY9LREY2jc4aIS705zEuO4P67b3seAgw8yKYPFdzRkSk3/ozOmtH2Ois/wBKUJ9IymRnpHPFWVNoeGkvnZ29NAhDM9gxcPxQMsITkREqw8tJbiIJPT8AfMkYMxFYOlCBSd+Vzp/KA8/v5rl3DrBwhi9+5QtugYv+LilxicjI5SmJuAstluIM6w1XBNzXv5DEqw/OnUp6mqFh+57ek0ia2wg90gaZY5yHiEgfeZknshRYD8xxH5PdxxygakCjkz7x5WZRcsYkGl5KoF8EYN/r8J2z4YV7BzcwERmxvLREzrfWXg1gjJkNYK3d6b5eCGwbqOCk78rmF/CtB17irbYjzMjrZWZ6fhFMmgVbNsDCTzl9JSIifeBldNbO0BM3eYT3g+T1OyLpl7L5BQCJtUaMgQtvgd3bYNcTgxuYiIxInnc2NMZ8zH1aYowZ7z5PaHSWMWa1MabcGFNhjKlIoL7PGFMdqm+M0SiwGM7IH8uZU8clNtQXYMEnIDcfnvz+4AYmIiOSlyG+9xpjvgqsdA+tA3YZY1qB/N7ON8ZUAwFrbZ21dj3OvJPyOPV9QK21tsqtD7Cmr3GPJqXzC3h6ZxsHjpzsvXLmGFhSAa89AgcTWDZFRCSMp5aItfZOa+1K93nAWpsHlFprE/nlXmGtDV+ocSNQGaf+BqAm7PUm1IEfV+m8Ajo6LY++ujexE5ZUwN83woTCwQ1MREYcz7ezIllrn+ltj/UYt6GCOMOFYykHGowxfmNMsbU2aK0NeI905Fs4w8fkcVk0vJRgEsnNg8lznOe97ZAoIhKmP30iEyIf9N5CyMNZBThc5Ovw9wglnZKwY7XuLS6JIT3NcNXZU3n0lb2caO9lj5GQjnbYeBM8esfgBiciI4qXeSI3GmPagK1As/sIPe+tk9wX57rRyvyhJ+5ts2ac218botQFoKWlBWNM12Pt2rW9hDQylc4r4P1j7TS+ETNHd5eeAZ2dsKUGThwZ3OBEZMTwMk/E7/aB9OB2uMcTpOcw4HjDgoPu1/CVggM4t7iiKiwspKVFHcSXnTmF7Iw06rfv4ZI5kxM76eJ/gFcegG2/cnZAFBHphZfbWc2xCqy1d/Zybhs9WyM+99xglPqBKGVBiNlyEdeYrHQuO3My9dv3YBPt55h5IZxWAn/9T+jsZV8SEREGsGMdwBhzVbxy93ZUMOJwHtAQo34ACEYkDB/OxliR15EIpfMKeCd4lJfffT+xE4yBi78IbQF45cHBDU5ERoQ+386y1m42xtzuvmzlVFLw4cwdWdzLJTYZY8rDhvmWETaE1xjjB4rDym8HVuCs14X7HqH3lziumjcVgIbte5g3fUJiJ529DEr/FWZcOIiRichI4aVj/Q6c1XoNpxZeDC3E6OvtfGttJeA3xpS6s9V3RMwbKSds3oi1dh3gc2e5rwZa3WPSi6njc1g4w5f4gozgdLBf+iUYN2XQ4hKRkcNLx3qjtfa2aAXGmMZELhAvCbhl66IcEw/K5hdw58OvsOfgMQom5CR+4mv18HYjfPCfBi84ERn2vPSJBGMVWGu1pvgQUzrPWZBxc6ITD0N2PQmPrXP6R0REYvCSRFpjzUw3xtzav3BkoJ1VMI4ZeWOo3/5u3068oBLSMzX5UETi8pJE/gmoNca8ZoxpNMY87D6a0MKIQ44xhtJ5BTyxo5UjJ9oTP3H8NLjoi/DcRtj11OAFKCLDmpck4gNuwxkxVeE+v819vnnAIpMBUzavgBPtnTz+2r6+nXj5rTDhdHjwq86yKCIiEbx0rFdZa5+JVhA29FeGkMWz8xifk0HD9j186JxpiZ+YNRau/za8/y6YAZ1SJCIjhJd5IlETSG9lkjqZ6WlcOXcqf3x5Lx2dlvS0PmyDO/fawQtMRIY9/Xk5SpTOm0rr4RNseyvo7QJbfw4P//OAxiQiw5+SyChx5VlTSU8zfZt4GK71dXjqB/B2U+91RWTUUBIZJSbmZrJkVh6bvSaRK1bD+OnwwFe0OKOIdFESGUVK5xfw6p5DvNnqYb+Q7PFw9bdg9zZo/vmAxyYiw5OSyChSGlqQ0Wtr5Nwb4YxLYfM34XiCKwOLyIimJDKKnJE/ljOnjvOeRIyB678D5T91WiYiMuoN9H4iWvZkiFs6r4AtO9s4cPSktwtMPRuK3G1jjiS49a6IjFhx54kYYx7uw7UMsAj4dr8ikkFVNn8qP35sB4+9+h4fWVDo/ULbfg0PrYHPPQRT5w1cgCIyrPQ22dAAVcRZuTdCdb+ikUG3cMYk8sdm0bB9T/+SyKxLISMHfnkjfL4BJvTjWiIybPV2O6vKWvuMtXZnvAfO3ungJBwZwtLTDB88eyqPvrKXkx2d3i/kmwmfqoVjB+GX5XDswMAFKSLDRtwkkugyJtbaAziJ5PyBCEoGV+m8Ag4ea6fxjX72aUw/D1b+F+x7BTbepEUaRUYhLwsw4u4nUkrP7XCLgPv6F5IMtsvOnExWRhqbX9rLxUWT+3exoqvghh/C0f3O1roiMqp42WN9KbCe7nurh/Za1+2sYWBsdgYXF+XT8NIerLX9v+CCj8OFX3CeH9zd/+uJyLDhZYjv+dbaq9191quBGmvtbdbaWwD/wIYng6V0XgG7Wo+w471DA3fRPdvhByXw9PqBu6aIDGleksjO0BO3U31pWFlevyOSpFjqzl6v397HvdfjmTIXZl8O//1VeHA1nDw6cNcWkSHJ82RDY8zH3KclxpjQ9OXi/ockyTB94hjOPW2C99nr0aSlQ/ndcMEXYEsN1FwBLdsG7voiMuT0OYlYa+81xnwVWOkeWgfsMsa0AvkDGZwMrtJ5BTS/uZ/WQ8cH7qKZOXDtHfDp38Dxg/CixlmIjGSeWiLW2juttSvd5wFrbR5Qaq1dM6DRyaAqnVeAtfDHlwfwllZI0VXwhSfhg+5GVu9shf1vDPz7iEhKDdjaWdbaZ4wxVw3U9WTwnVM4gekTc6jfPoC3tMLl5kFGNlgLv/sH+NGlsGUDnPCwFL2IDEle54lcRc+RWD6cW1yL+xmTJIkxhg+dM41fb3mT94+dZHxO5mC9EXzyHvjNF+DBW+GP/wYLPwVLKiBv9uC8p4gkhZd5IncAt3Bqnkj4fBFfgtdYbYwpN8ZUGGMq+vj+NX2LWOJZtqCQE+2dg9caCfHNhM/eD599AIqWwpb1sOdFp+zEYc12FxmmvLREGt05Ij0YYxp7O9kYU+1eoy702hhTHnqdwLmaizKAimf6OM03ht8/28LHik8f3Dczxlm4cdal8P4eyHXHYTzxXWj+Bcy/Ac64CGZcCOMLBjcWERkQXvpEgrEKrLX3JnB+RUTC2AhU9naSMUbDhweBMYZlCwr5y2v7aDt8InlvPL7g1DIpMy+EgnNg689g02fgO2c5w4NDs+kPt0JnPxaLFJFB46Ul0mqMmWWtfSOywBhzq7U25n4iMRJBEGcdrt6UAPVoLsqAW7ZgOj9+bAcPvfAun7xgZvIDKLrKebSfgHefg11POhMVjXHKf/5haH3duSU2aZbzmHkRfKDcKX//XcieAFm5yY9dZJTzkkT+CZhtjPHhJIDQUrD5wGzib0qVF1Y/pNelZI0x5cAmnEQiA2z+9An4p4zlD8+2pCaJhGRkweklziPEWrjo72Hfa84Q4f1vwNtNcPzQqSTyvWI4edjZ32RMnjMq7LwVcMk/Oi2YB/4XZOZC5hgn0WTmwumLnfdpPw6vPAjpWc4jLcN55M2Giac75Xu3g0l3JlOaNOcxbiqMmeSUH2xxE5459TU3D7LGOuWhHSBDSREDOROdOTXtx53l9ENCdbInOP8f7cedf2uknAmQngknj8HJKKPdsic4Lb1Y5TkTnX/PyaPRVxbI8UFamjOSrv1Yz/Ixk5xYTxx2YoyU6y5ecfwQdES0cI1xzgc4/j50ROyyadJgjM95fuwgdEb0l6WlO/GDswVBZ0dEeYbz/xOrPD3z1PbOR4NgO+OU7z/VIu4qz4Lscc7zaLt7ZmQ733trnfN7lOc4P4ednXAs2LM8c4zz6OyIvsVCZq7zs9PR7szFipQ11omh46Tz/9ujfJz7s3UCTrg/W8b7QF0vScQH3EbPX/7GPd7buVEZY3zW2mC040DQWhs0XR/C2FpaWgiv941vfIO1a9f2et5oZoxh2XmFfO+Pr7H34DGmTshJdUinGAPn39TzeOgXj7XwoX93PqxH2+CI+zVrrFPefgxevt/5RXniMOD+QrjsVieJHDsItZ/tef3StXDpl+HgO7D+yp7l130blqyC916Bmst6ln/0x7DwE/BOM9x9Tc/ylb+Eectg55/hV+U9yz/9Wyj6oJPgosX3+c1O/M/Xwu+/2LP87/7q7Di59W54KMrH8ksvgG8GPPkD+NO3epZXveH8on/sDqfPKtLXWp0kVf91aPxJ97KMHPgXd6DGA1+B5+7pXp6bD6sDzvPf3OJ8f8L5zoAvPec833gT7Hyse3nBufCFJ5znv/gfzhykcDMuhL91N2W962p47+Xu5XNK4Sb3zvuPLoGDb3cvn38DrPgv5/l3F/T8Rb7wJvjoD53nd84BG5GkLrgFrq12fvbWRRl9eNmtsPRrzs/pnUU9y0M/e8E34XsLe5Z3/ey9BD++tGd56Gfv7aZefvYeO/WzN+G0nvUS5CWJVMXaZ8QYc3sv5wbpub5Wb+ttrbDWJryiX2FhIS0tLYlWF9eyBdP57ubXeOD53dx8yTAYdpvuDkc2Bkpujl0vKxe++rrz3Frnr+aTR06dP2aS8wu344STmDpOOH/5TprllI8rgE/c4/xV2NkOWOcv1+kLnfKJpzsfWtvpllnn64wlTnmeHz78H3Qlr9BftQXnOl+nzHV+KUSafKbzddp5cO2dEYUWJs5wnp6+GK5dF1ZkT8UNcMYlcE2UDUdDf+kXXXXqr+5wGWOcr3Ovg/FRdq0M/eU6/wbIP7N7WVr6qefnrYDCiG2GMsP+SDn/0zArIgmHx7NklRNDuNywhTEu+ns49F738vBBGZf+r56tAd+MU8+vvM394yJMXtjYndK1zl/s4abMPfX8mtt7tlSmud/btMzo//enuXfkM3Ojl4d+dnLzopfPvMj5On56/OtPmhW9vOAc5+vks06VZ40FPtOzbgKM16XAjTGfx9lTPQ+ot9b+pJdTQn0iW621Jt6xiLKgtTbgvi7FSWJlsd6jpKTENjU19fnfI3Dtdx9nTGYa9/3dJakORUSSzBiz1Vrb5y4DTzfCjDGP4PRPBIAmnEUYG40xE+KdZ61tpuforjygIcYpeUC5O69kNc4oLr/7WkN9B9iyBdNpfjPIW22aUS4iienz7Sy3BbLc3RI3/LgPqCB+xzrApoh5IWVA1wRCNzkUW2vrrLUNhCUYd2Ki31q7Dhlwy84rZN1Dr/DA87u55Yoo92pFRCJ4aYnsj0wgAG6n+M6e1XvUC7UmSt2ksCNi3kg5UeaNuHWXc6ol4vMQu8QxIy+XhTN8/H6b+pREJDFeOtbjdaIk1MESryXhlvUodzvXtWXeIFu2oJB/u387r+89xJyp41IdjogMcV5aIvnGmFmRB40xCwHdAxnmPnzedIyB+59Ta0REetfnloi1doMxZpMxZjZOxzo461kFQnuMyPBVMCGHC2bn8YdnW/jHpWeSyNwcERm9vG5KtQKnE70BZ3RWhRLIyLFsQSE73jvMS7ujzHYVEQnjea67tfYZa+0Gd5fDqJMPZXi69tzppKcZ/qBbWiLSiwHb2RC6hv/KMJc3NotL50zmD8+24HUyqoiMDr32iRhjfgTUWGu3ua8fjlUVZwZ7rzPXZehbtqCQW2ufZdtbQc6fOSnV4YjIEJVIx3rknBADVNFz5rkB7hiAmGQIuPqcArJ+k8ampreVREQkpl6TSJRdDCuttVEnFRpjqgYkKkm5CTmZ/I+Fp/GbZ95m9YfmMmlsVqpDEpEhyEufSGv4C2PMbGPMjcaYq2IlFxmePnfpbI6d7OTXW95MdSgiMkR5SSIV4S+stTuttfdaa/9ojPnYAMUlQ8DcaeO57MzJ/NdTb3CiXdvTikhPAzo6i973BpFh5nOXzGbPweM8+PzuVIciIkNQIqOzJgIrcFbbnYizAGK0/Tz8hK3GKyPDFWdNwT9lLD99Yic3LCzUDHYR6SaRjvUDwAZggzGmGmel3mjJIhBtdV8Z3tLSDJ+7ZDb/8tsXaNq1n8Wz1NgUkVP6dDvLWluFs4vhM1EeSiAj1MeKT2PimEzuelzjJkSkuz73iVhr7zXGrArfT90YM1Gd6iNXblYGn7xgJo9sf1e7HopIN31OIu7SJg1AW+iYtfaAtfY+JZKR6zMXnUGaMfzsyTdSHYqIDCFedzbcaa29M0qZbmmNUNMnjuG6D0xnY+NbvH/sZKrDEZEhwksSibci32yvgcjQ97eXzubQ8XZqm95OdSgiMkR43dmwx2q97jFfvyOSIWvBDB8lZ0zi7id30tGp1X1FxFvH+gZgjjGmzRjT6D5agUXW2m8PfIgylHzu0tm81XaU+u17Uh2KiAwBfd4eF5xFGY0xdwBL3UOvAwcHLCoZsq6eX8BpvjH89ImdXHPutFSHIyIp1p+dDYPumln3WmufBVo1Omvky0hP4+ZLZrFlZxvNb+5PdTgikmKeWiLGmFlAKT37QIqA+/oXkgx1H18yk/V/DvD1373Ab//uEjLSB3oJNhEZLrzME1kKrAfmuI/J7mMOzmZVMsKNy87g68vm88I7B/nFX3elOhwRSSEvLZHzrbVXg7OXCDjLwbuvFwLbBio4Gbqu/8B0as96m+888irXnjudaRNzUh2SiKSAl/sQXQsoucljaViZVucbJYwxfPOGczjZ0ck3738x1eGISIp4vpkd1oleYowZ7z4vTvDc1caYcmNMhTGmope6Prf+amNMbW/1JXnOyB/LP1w1hweff5c/vbw31eGISAp4XYDxq8BK99A6YJc7VyS/t/Pd5eQD1to6a+16oMgYUx7nlDXW2nXuYzlQpUQydFRcXsScqeP42u9e4OiJjlSHIyJJ5qklYq2901q70n0esNbmAaXW2jUJnF5hra0Le70RqIxW0Rjjw9nsKlwN6sAfMrIy0vjWR8/l7f1H+f4fX0t1OCKSZF5GZ20yxtwaedxa+0wC50a73RXEGS4cS6kxJjyRBOmZWCSFLvTnU77odNb/OcCre95PdTgikkReWiL1ODsd9mCMmdDLuXmELSHvinzdxZ3QOMlaGwg7XIazFH1ULS0tGGO6HmvXru0lJBkIa649m3E5Gfzzb56nU+tqiYwaXpLIDmBSjLLe+ip8sQrcW1dxuXVKiXM7q7CwEGtt10NJJDnyx2Wz5tqzaXxjP3XNWuVXZLTwMk9kBbDI/YUewLm9BE5iWQTEW4QxSM9hwH0ZFrwBWG6tbe7DOZIkyxfNoG7r2/zvB19iyaw8Zk0em+qQRGSQeUkiJTgtgcjbUAa4rZdz2+jZGvGBc+sq3onGmNVAjbU25q0sSa20NEP1jedx44+e5Ka7nqbulos1CVFkhIt7O8sYc36UTvRV1trN1tpnIh7NwO3RrhPi1glGHM4jTh+HG0c50BxKIMaYeB3xkkL+KeP42c1L2H/4BJ++62n2Hz6R6pBEZBD11idSitPCCBdv98IdCbznpoh5IWU4w3YBMMb4w8vdhJEHNLkTD/0kOKlRUmPBDB8b/qaEXW1H+OzPGjl8vD3VIYnIIEmkY70m4nW84bW9TgK01lYCfmNMqTtpcEfEvJFy3Hkjbr9LvRvDfvexA1icQNySQhcXTeYHnzifF945QMUvmjh2UhMRRUYiY23s4ZjuAou1wEScTnSD0xIJRKuOs7thr7PWB1NJSYltampKZQgS5t6tb/OV2me5en4B//mpYi0bLzJEGWO2WmtL+npe3I51d4HFEmPMRE6NoqqkZ+sEnCRyR18DkJHtxkWnc+DoSb55/3Zuu+951t14HmlpkXdIRWS4Smh0lrX2AHAAwBhTH1r6PZIxRsuRSA+fu3Q2B46e5LubXyMnM42vf/gcsjLUIhEZCfo8xNdauzlOWdTkIvKl0jM5cqKdDY/vpOmN/fzflQuZN723BQ5EZKjTn4OSFMYY/vn6+Wz4TAn7Dp3gIz/4Cz/80+u0d3SmOjQR6QclEUmqsvkFPPLlyymbX8CdD7/C8pqnCLx3KNVhiYhHSiKSdHljs/jhJ4v57scXEnjvMNd973F+9sROLdwoMgwpiUhKGGO4YeFpPPLly7lgdj5r/7CdD3//L2xqfEtzSkSGkbjzRIYjzRMZfqy13Nf8Duv/HOCVPe/jy81kRckMbrrgDGbm56Y6PJFRwes8ESURGTKstTy9s43/euoNHn5xD53WctXcqXz6ojO4dM5kTVQUGUSDMtlQJJmMMVzoz+dCfz67Dxzl/z39Jr/e8iab797L+OwMlszO46Iip3z+9AmatCgyBKglIkPa8fYONr+0l8df28dfA63s3HcYAF9uJhfMzuOC2fnMnTYe/5SxTJuQgzFKLCJeqCUiI1J2RjrXfWA6131gOgC7Dxzlr4FWntrRypM7Wnn4xT1ddcdmpTN7ylj8k8fhnzKWWfljmTI+m8njspk8LotJuVlqvYgMMLVEZFjbc/AYO/YeYse+w+zYe4jAvsME3jvEO8GjRP5op6cZ8sdmMXlcNnljsxiXncHY7AzG52QwLjuDcTnO69zMdLIz08jOSCc7I815ZDrPM9MNGWlpZKQbMtPTyEgzZLhf09MMaSb0FbWKZFhRS0RGpYIJORRMyOHiOZO7HT92soO39x9l36HjzuP947x36Dj73j/Be4eOEzxygvfeP86h4+28f+wkh463M9DTVIyBdGNICyUVnK9pxmCMsxOk6arrPHfyjlNuwq7jHg17fuq8aO8bL6aox3tsG9R/yqGjg5KIjEg5menMmTqOOVPHJVTfWsvRkx0cOtbOsZOdHG/v4Hi787Xr9clOTnZa2js6ae+wnOx0v3Z00tFp6bCWzk5LRyd0WkuntV3Hsc4xa6HTnip33hss1v2K24I6VRZep9uxbvGHjsXJhDGKBuNexEi7wzEa/NnjeUoiIjh/0edmZZCbpY+EjE7f/YS38zTwXkREPFMSERERz5RERETEMyURERHxTElEREQ8UxIRERHPlERERMQzJREREfFMSURERDxTEhEREc9SssaDMWY1EADyAKy16weyvoiIJEfSWyLGmGogYK2tc5NBkTGmfKDqi4hI8qTidlaFtbYu7PVGoHIA64uISJIkNYkYY4qjHA4CpQNRX0REkivZLZE8oC3iWOTr/tSnpaXF2eDHfaxdu7bvUYqISEKS3bHui1VgjPFZa4P9rE9hYSEtLS0ewxMRkb5IdkskiDvCKkzk6/7UFxGRJEp2EmmjZ+vCBxCtVeGhvoiIJFFSk4i1thmndREuD2gYiPoiIpJcqRjiuylinkcZUBN6YYzxR5THrS8iIqmT9CRira0E/MaYUmNMBbAjYh5IOWHzQBKoH5VGZQ1f+t4NX/reDWuFXk4y1tqBDiSlSkpKbFNTE8YYRtq/bbTQ92740vdu+HK/d6av52kBRhER8UxJREREPBtxt7OMMe8Bu3Du72nW4fCk793wpe/d8DXPWju2ryeNuCQiIiLJo9tZIiLimZKIiIh4piQiIiKeKYmIiIhnKdljfTBpP/bhyV2NYBFQ6x5aDlRbawOpi0qiMcb4gAog31pbFaVcn8EhKt73zutncEQlEXc/9sbQsijGmGpjTHkiy6TIkLAC5we8GVilBDL0GGNKcVbSLopRrs/gENXb987V58/giEoiOPuxh2fXjUA1oB/gYcBaOynVMUh81toGAGPMYqJvGqfP4BCVwPfO02dwxPSJaD92kdTSZ3B0GkktkT7vxy5Di3tPtg3dSx+u9Bkc5rx8BkdSEvHFKoi1H7sMKU1AMHQP1hhTa4xp0730YcUXq0CfwWHB02dwxNzOQvuxD2vW2uaITrxGYE2q4hFPgugzOGx5/QyOpCSi/diHMXfkSLgAEO0euwxd+gwOY14/gyMmiWg/9uHLGOMH6t0x7OE0xHcY0Wdw+OrPZ3DEJBGX9mMfhtwmdFXEX6srcYaGyvCiz+Aw1J/P4IhbCt6dLdsM+EEjfIYL9y+h0C+ffGCHvndDjzuMtxSodA/VAA1uKyRUR5/BIai3753Xz+CISyIiIpI8I+12loiIJJGSiIiIeKYkIiIinimJiIiIZ0oiIiLimZKIiIh4piQi4kGUmb2hDZg0sU5GlZG0iq9IMq0AIidibUxFICKppJaIiDdlkQfcVVCbo1UWGamURET6yN1H3JfqOESGAt3OEukDd3FBH+B314gC57ZWHu5Cg9baMnedomqcVW1vd8t9wGJrbVXYstvFQCBy45+w9ad8QJ7Wn5KhSmtnifRRKEFYa8viHXcTRTWwPHy3OJykURV23n5r7aSw17XA7WEL41UDjdrlUYYi3c4SGTjBiNdtgC9it7ho+zO0hUZ7uSupFkf0rWzk1MqrIkOKbmeJDK5glGOtceqXAsGIXeZ8uMuqiww1SiIi/WSM8Ue0NsK1RTkWjHM5H87trsjdAHUrS4Yk3c4S6b+B3Au+azMnkeFASUSk7wJ0/0UfjFM3L8oxX6zKbgukze2k72KMqehDfCJJoyQi0kfuPtQ1xpgKY0y5tbbB7RCvBkrc48XAGsKGArvDg8uBlaE+D7fMD1S718Ad3VUaur77HhriK0OShviKiIhnaomIiIhnSiIiIuKZkoiIiHimJCIiIp4piYiIiGdKIiIi4pmSiIiIeKYkIiIinimJiIiIZ/8fu+pWqXZxMe4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ylim = ax.set_ylim(-0.05, 1.1)\n",
    "xlim = ax.set_xlim(-0.1, 15)\n",
    "__=ax.plot(tarr, mass_loss_target, label=r'${\\rm target}$')\n",
    "__=ax.plot(tarr, frac_mass_loss_init, '--', label=r'${\\rm initial\\ guess}$')\n",
    "\n",
    "xlabel = ax.set_xlabel(r'${\\rm time}$')\n",
    "ylabel = ax.set_ylabel(r'${\\rm fractional\\ mass\\ loss}$')\n",
    "leg = ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unavailable-chair",
   "metadata": {},
   "source": [
    "Now define the `train_step` function based on the previously-retrieved Adam functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "compressed-lesson",
   "metadata": {},
   "outputs": [],
   "source": [
    "from toy_model import _loss_and_grad_func\n",
    "\n",
    "def train_step(step_i, state, train_step_data):\n",
    "    \"\"\"This is unfinished\"\"\"\n",
    "    params = get_params(state)\n",
    "    t_target, m_target = train_step_data\n",
    "    loss_data = t_target, m_target\n",
    "    loss, grads = _loss_and_grad_func(params, loss_data)\n",
    "\n",
    "    return loss, opt_update(step_i, grads, opt_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "textile-conservation",
   "metadata": {},
   "source": [
    "Finally, set up a simple loop where we take one step down the gradient at each step, with a step size that is adaptively set by the Adam algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "living-madagascar",
   "metadata": {},
   "outputs": [],
   "source": [
    "tstep_data = (tarr, mass_loss_target)\n",
    "\n",
    "n_steps = 400\n",
    "loss_history = []\n",
    "for istep in range(n_steps):    \n",
    "    loss, opt_state = train_step(istep, opt_state, tstep_data)\n",
    "    loss_history.append(float(loss))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cross-vehicle",
   "metadata": {},
   "source": [
    "### Inspect the loss curve\n",
    "\n",
    "This should flatten over the course of training, or otherwise it training has not converged and should continue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "stunning-flash",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEOCAYAAABbxmo1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAnwElEQVR4nO3de3QU150n8O+v9QSE1GqBACFASEAwtgELgR/rV0BykrGTSczD2cxsdpO1wT7Jrmd3E7DP7KyZnTljC2eyO5tsbCCb8czueAKI7CbjSYLBjh+xnRghMLEdG5DAPGTzkhrxkNCjf/tHVYtSqx9V/apW9/dzTh911b1d+unG6R+37q17RVVBRESUCh63AyAiouzFJENERCnDJENERCnDJENERCnDJENERCmT73YAmWbSpElaU1PjdhhERGPGvn37zqnq5HBlTDIhampq0Nra6nYYRERjhoh8FKmMt8uIiChlmGSIiChlmGSIiChlmGSIiChlmGSIiChlmGSIiChlmGSIiChlmGSS5OldH+D//CbiVHEiopzEJJMkb7afx8/e6XQ7DCKijMIkkySLqr343ckLGBwKuB0KEVHGYJJJksUzvOgdGMKRs5fcDoWIKGMwySTJohleAMA7J/yuxkFElEmYZJKkpmI8ysYV4ACTDBHRMCaZJBERLJrhxYETF9wOhYgoYzDJJNHi6jIcOn0RV/oH3Q6FiCgjMMkk0aIZXgwFFO919rgdChFRRmCSSaKF1V4AHPwnIgpikkmiyROLMN07DvuZZIiIADDJJN3iGV72ZIiITEwySbZ4hhcnu3tx7tJVt0MhInIdk0ySBR/KPHjS72ocRESZgEkmyW6YXgqPAAeO+90OhYjIdUwySTa+MB/zpkzEgZN8KJOIiEkmBYKD/6rqdihERK5ikkmBxTO8uNA7gI/OX3E7FCIiV+W78UtFZD2ADgA+AFDVLYnUt1nuNw+9qrop4T8iiuDg/4ETftRMmpDKX0VElNHS3pMRkWYAHaraYiaDOhFZFW99G+XrVXWTqm4xy/eYSSdl5laWYHxhHldkJqKc58btsrWq2mI53gZgXQL1Y5U/YL2YqrYBWOosZGfy8zxYVO3Fvo+6U/lriIgyXlqTjIjUhzntB9AYT32b1+sSkR2Wa66FkYhSqn6WF+9/3MMVmYkop6W7J+MD0BVyLvTYSX0711sHoFFEus3bZF0hPZ8ROjs7ISLDr40bN0YJL7Ils8oxFFAc5FRmIsph6U4y3kgFIhKuLFb9mNdT1Q4ATwJoBdCMGLfKqqqqoKrDr3iTzE0zygEAbcd5y4yIcle6k4wf5gwwi9BjJ/VjXk9ENgPYo6pNAJoArLXePkuV8gmFqJ08AW0clyGiHJbuJNOF0b0PLwCoqj+O+lHLzTEbvznYD1XdA2A2IowBJVv9zHK0HedDmUSUu9KaZMwve3/IaR+APfHUt3E9H4DzIdf0R/p9ybZkVjm6LvfjGB/KJKIc5cYU5u0hz8U0AdgcPBCR2pDyqPWjlZs9lybrLzfHajoS+gtsqp9pjsvwlhkR5ai0P/GvqutEZL2INAKoBdAeMttrFYzE0GKnvo3rrTMf2Gy3fGZDyv5Ai7mVJZhYlI+2491YuaQ6Hb+SiCijuLKsTLRlXcyyTWHOxXu9DgBpSSqhPB7B4pl8KJOIchcXyEyx+pnlOHT6Ii72DbgdChFR2jHJpNiSWeUIKPDOCT6USUS5h0kmxRbP9EKED2USUW5ikkmx0uICzKucyHEZIspJTDJpUD/Li/3HuxEI8KFMIsotTDJpUD+zHD19g2g/e8ntUIiI0opJJg2WzDIeymzlLTMiyjFMMmkwe9IETCopxN6j0XY1ICLKPkwyaSAiWFrjw2+ZZIgoxzDJpMnSGh9O+XvR6e91OxQiorRhkkmTZbONbW72HmNvhohyh6MkIyKLUxRH1rtuWikmFuXzlhkR5RSnPZmtIlKakkiyXJ5HUD+rnIP/RJRTnCaZbgCNIrJSRJanIqBstmy2D4fPXELX5X63QyEiSgtHS/2r6j3B9yJSJiIrASiANlU9luTYso51XOYz1091ORoiotSLe+BfVS+o6k4A+wFsEZFdInJ/8kLLPgury1CY7+EtMyLKGU4H/mss7+8XkRcBvAhgN4A1qvoTEVnBZBNeUX4eFs/wcoYZEeUMpz2ZHSLyjIh0AXgAwAZVnauqT6vqBQBQ1ZfMZMNEE8ayGh/e7ezB5auDbodCRJRyTpNMOYB9qupT1QdUdX+4SuZUZ3+CsWWlpbN9GAoo95chopzgNMk0q+oPbdR7OJ5gckH9TC88Ao7LEFFOcDS7DECTiJSp6neiVVJVJpkIJhYX4PqqMj6USUQ5wWlPZjeAreEK+JCmfUtrfDhwwo+rg0Nuh0JElFJOk0w7jHGZcNYmGEvOuLnWh6uDAbxz4oLboRARpZTT22VrACwRES+ADlwb3C8HsARA1NtoZLhldgVEgDfbzw0/oElElI2cJpkGABsAhA4oCIDHkhJRDigbX4Drq0rxVvt5/Emj29EQEaWO0yTzUJRpy08mIZ6ccVvdJDz3xjH09g9hXGGe2+EQEaWEozGZYIIRkQfNhzK3iciD1jKy59a6CvQPBbDvIz4vQ0TZy/HaZeZSMg0wxmRaATSIyF7OLnNmaY0P+R7Bm+3n3A6FiChlHN0uM3stq4NLyFjOe2HMLuPAv00lRflYNMOLN9vPux0KEVHKON5PJjTBAICq+gEcTUpEOeS2ugocPOlHT9+A26EQEaWE0ySjcZZRGLfWVSCgXGKGiLKX0yRTYV3uP8hcELMuGQHlkvqZ5SjM9+At3jIjoizldGfMrSKyXURmwxj4B4BaAB2q+oDd64jIevPzPvO6WxKpb6PcC+BxGCsWAECrqrbZjTdVigvy0DCrnOMyRJS1HM8uU9U1MAb598CYXbbWYYJphpGUWsxkUCciq+Ktb6PcC2CHqm6wJJ/Hbf/BKXZrbQXe/7gH3Zf73Q6FiCjpnO6MWQoYz8So6lYALQBqRWS5g8usVdUWy/E2AOsSqB+rfCuAzZbj7TBWLcgIt82pAAD8poO9GSLKPk57MiMWwVTVo6q6U1VftrMTpojUhzntBxB2cZVY9W1ebxWAPSJSKyL1qupX1Y4wn3PFwmovxhfm8ZYZEWUlx7fLorCz0qMPo9c9iza1Klb9qOWWJNRgObfDvIWWEQryPFg224c3+FAmEWWhqElGRMpE5CFzsH8XgHUisivM6zAAr43fF7FOhC/+WPVjldcGj1W1wxzs34YIe+IAQGdnJ0Rk+LVx48ZIVZPm9jmT0HH2Mjr9vSn/XURE6RR1dpn54OVWAFvNAfajGDm+EdQR7iHNMPwY3eOJ1gOKVd9OOWBMUAjqgHELLayqqip0dnZGCSn57pg7GcDv8dqhs/jysplp/d1ERKlkewqzqm4QkZUJLoTZhdG9D695fb/T+iISq7wjzLX9gNHTifA7027elBJMLS3G64fPMckQUVZxugrzzkR+mXm7yh9y2gdjOrTj+jbKOwD4Q27FeQH4MyXBAICI4I65k/DrI+cwFODCCUSUPZxOYX5IRJ40n/CHudz/LvNnjc3LbA95LqYJlltw5iywVXbr2yh/EsaOnkEPmOcyyp3zJuNC7wDeOel3OxQioqRxOrusC8BTqnpARB4CUKuqn1HVRwCEm048iqqug/FsTaOIrAXQHvKcyypYnnOJVd9G+SYAXhFZb64McN48l1FunzMJIsBrh866HQoRUdI43RnTbxngXwWg2VJmZ+AfwPAXf7SyTWHOxXU9O+WZoHxCIRZOL8Nrh87iTxrnuR0OEVFSJLIKcxNGztriYEKC7pw3GQdO+HGhl0v/E1F2cJpk6kSkRkSeAtCiqj3A8GZmlKA7501GQIE3j/DBTCLKDk5nl22F0YNpV9U15sOaTwGYA5tjMhTZ4hlelBTl47XDHJchouzgdEwmmGiC7y8AeAxgbyYZCvI8uK2uAq8dOgdVhYi4HRIRUUKiJhkReQbAZlU9YB7vilQVwBIAP0xqdDno7k9V4sX3T+PwmUuYN2Wi2+EQESUkVk8mdMaYwFgm3x/m/FNJiimnLZ9fCQB46fdnmGSIaMyLtXbZYyGn1qnq0XB1RSRj9mgZy6aWFWPBtFK8/MFpPHI3d7QmorHN6cB/2ARjmp1gLGRacV0l9n3UDf8V7pZJRGOb44F/cxfM2pDTXhjLtSxNQkw5b/n8Snzv5SN49dBZ/OHi6W6HQ0QUN0dJxpyuXAtjufxQ3mQERMCiai8qJhTipd+fYZIhojHNaU9mb5hxGgCAiOxNQjwEwOMRfHp+JV587xMMDgWQn5fMDUyJiNLH6beXP1JBotsA0EjL51eip28Qbcf9bodCRBQ3p0nmfKQl/UXkW4mHQ0F3zJ2EfI/gpQ9Oux0KEVHcIt4ui/DgpcBYVl9h9Gq6zPMVMGaXfSfZAeaqicUFuLnWh5d/fwaPf+46t8MhIopLtDGZ4IOXdgjM5WUoeZbPn4K/eOF9HDt3GTWTJrgdDhGRY9GSzAZV3W/3QiKScbtNjnX3LDCSzK73PsG6u/hgJhGNPRHHZJwkmHjqU2wzfONxfVUpdr33iduhEBHFhXNjM9xnr5+KtuN+nOnpczsUIiLHmGQy3GdumAoA2PU+Z5kR0djDJJPh5laWYPakCXiRt8yIaAxikslwIoJ7rp+Ct9rP48KVAbfDISJyhElmDPjs9VMxGFA+mElEYw6TzBiwqNqLKaVFnGVGRGOOoyQjIg+JyJMistg8fkZEdpk/a1IRIBkLZt6zYCpePXQWV/oH3Q6HiMg2pz2ZLgBPqeoBEXkIQK2qfkZVHwFQn/zwKOhzN05F30AAL39wxu1QiIhsc7wKs6peMN+vAtBsKbsQpj4lyc2zKzB5YhH+6Z1Ot0MhIrLNaZJRy/smAK0RyijJ8jyCe2+chl99eBYX+zjLjIjGBqdJpk5EaswdMltUtQcAROTB5IdGoT6/aBr6BwPYzQcziWiMcJRkVHUrjB7MEVVdIyJl5sKYdeCYTMrdNKMc073jeMuMiMYMx7PLYCSU4G2yZgBLAHgBtCQ1MhrF4xHct3AaXj98Dv4r/W6HQ0QUUzyzy560zC6brar3cHZZ+ty3sAqDAcUv3+UzM0SU+Ti7bIy5YXopairG44WDH7sdChFRTNE2LQsndHbZ6ghlUYnIegAdAHwAoKpbEqnv5HoisllV19mNNdOICO5bWIUfvHIEZy72oXJisdshERFFlPbZZSLSDKBDVVvMZFAnIqvire/kembdWruxZqov3lSFgAI/O8AJAESU2RKdXVZqJpw5sD8ms1ZVrZMEtgGI1rOIVd/W9UQka8aM5lROxKLqMrTsO+l2KEREUTleINNMNCIizwDYCiPhPKaq34n12Qhf9H4AjfHUd3i9BgC7Y8U4VqxcUo0PPrmI9zo5FEZEmctxkhGRF2H0WjpgTGVuEJG9IlJq4+M+GDPUrEKPndS3dT3z9tl2G/GNGZ9fWIWCPMHOfafcDoWIKCKnz8k8CGC1qj6iqk+br4dh3EJba+MS3ijXDlcWq37M65k//arqtxEfOjs7ISLDr40bN9r5WNqVTyjEivlT8NMDpzAwFHA7HCKisJzOLuu2TGEepqp+ETlq4/N+mDPALEKPndS3c701sWavWVVVVaGzc2wMqK9cUo1fvvcJXv3wLBoXTHE7HCKiURJZINNJWVAXRvc+vICRqOKoH7XcHLPZYyOuMenuT02Gb0IhdrZxAgARZSanPZkKEalR1WPWk+YmZnWxPqyqbSLiDzntQ4REEKu+jev5ADSKSLBsKYBa87maFlXtiBVzJivI8+ALi6rw/G+Po/tyP8onFLodEhHRCPFMYd5kDvRvM197ATyuqk/bvMz2kOdYmgBsDh6ISG1IedT60cpVdY+qbgq+YMwu85vHYzrBBK1pmIH+oQB+sp8TAIgo88QzhXkNjEH+PTBml61V1QccfH4djN5Eo4isBdAe8pzLKliec4lV38b1AABm2Wqz7voIEw3GnAVVpVg0w4vnf/sRVLmlDxFlFknWF5OIfMvOszKZrqGhQVtbW2NXzCDb957A+p0HsW3tLbi5tsLtcIgox4jIPlVtCFcWcUxGRHY5+R0wlvwf80lmLLpv0TT8xQvv4x/fPs4kQ0QZJdrAvwDYAGOacCwC4KlkBETOjS/Mx5fqp+PHe0/gCU4AIKIMEm1MZoOq7lfVozZeHTASErnkKzfPRP9ggNOZiSijREwyqrrfyYVU1c7DmJQi86eWon6mF8+/fZwTAIgoYzieXUaZ649unoWOs5fxZvt5t0MhIgLAJJNV7l04Db4JhfjbN9ipJKLMwCSTRYoL8vDHN8/ESx+cwbFzl90Oh4iISSbb/PEts5DvETz35jG3QyEiYpLJNpWlxfj8wirsaD2Bnr4Bt8MhohwXNcmISI2I3CQi94duSiYiK0RkpYh8W0QeFJHlqQ2V7Prav5iNy/1D2L73hNuhEFGOi5hkRGQIQDMAVdWfqGqPtVxVX1LVnebCmD5k0dbGY92N1WVYVuPDc28ew1CA05mJyD3RejL7VfUBVT0Q6yLmCseOnquh1Pr67TU42d2LX777iduhEFEOi7aszPAqkSJSBqARIzcm6whJQGNrVcks17RgKmZPmoAfvHIEf3DjVFj21CEiSptoPZn24BtVvaCqOwFUAHgcwJ4wPZx2UMbI8wgeuasO73X24NVDZ90Oh4hyVDyblr0UOj5DmemLN01HVVkxfvAr5n8icke0JBNpxPicw/rkksJ8D9beWYu3j3Xh7aNdbodDRDko2pjMl0Uk3OYkjRHOrwL3k8k4Dyydie+9fAT/81dHsGz2MrfDIaIcEy3JeAHUhTl/NMJ5XzICouQaV5iHr98+G0/v+hAHT/qxsNrrdkhElEOiJZnN5jMwtojIt5MQD6XAv7p1Fra81oHv7j6E577G3gwRpU+0/WRsJ5h46lP6lBYX4OG76vDKh2fReoxjM0SUPnGvXWYuObM8dLkZykz/+rZZmFRShE27PuSmZkSUNtGWlXlIRPaar2+JSI15vkxE9gLYB+BhAD8UkWfSEy7Fa3xhPr756Tq8fbQLrx+ONEGQiCi5oo3JbAfgDXMbbAcAqOrwDDMRmS0i31JVzi7LYP/y5pnY+vpRfOfFD3HH3ElcBYCIUi7a7bKHQhOMubzMCgCrredV9SgAfmNluKL8PDy6Yi4OnrzANc2IKC2iJZkLYc6tAeBX1WNhynijfwy4v3465laW4MlffICrg0Nuh0NEWS5akikLc241jNto4bAnMwbk53nwZ/ctwPGuK3jujWNuh0NEWS5akhERWWw5WAljJebNYSouhzERgMaAO+dNxvL5lfjey0dw9uJVt8MhoiwW6zmZh0Vkl4i0AtgCYLV19WVzd8xnESbxUGb703uvQ9/AEL67+5DboRBRFov6nIyqPqyqn4GRXCrM5f4BDE8C6ICxe+Y9MJaboTGibnIJvnprDbbtPY73O7moNhGlRrTnZGqC783ZYyOYe8wcDb4AnE9NiJQqj66YC+/4Qvzn//c7BLhNMxGlQLSezDqH11qbSCCUfmXjC/Cnf3Ad2o778fzbx90Oh4iyULSHMTeIyCoH1/KBS/2POffXT8fOtpNo/uUHuGfBFFSWFrsdEhFlkWhJ5mkYz77sBeCPcR2Bg56MiKyHMZ7jAwBV3ZJI/WjlIuK1xLYUwO5Yvy+XiAj+8os34LN/8zr+6wvv4/tfqXc7JCLKIhGTjKpuAAARuQnGl7eq6suR6psTAWISkWYAe1W1JXgsIquCx07r27je48G/xSxvF5GYiS2X1E4uwTc/PQff3X0IK+vP4NPzK90OiYiyRMxVmFV1v6q+pKovm1OWl1ufn7HU2xnm4+GsDUko2xB9/CdW/YjlZi+mNuR6mwFsAI2w7q5azJtSgsd+chD+K/1uh0NEWcLRUv/BZKOqB8yEsyJcwolERMLdi/HDeMjTcX2b12sUkdqQ8tDEk/OK8vPw3TWLcf5SP/7LT99zOxwiyhJx7ydjJpyXANQFtwOw8TEfgNBds6LtohWrftRyVfWrarmqdljKmwDssRFrzrlhehkeXTEXP3unE//0Tqfb4RBRFogryYjIYhF5VkTOA3gMxnpmW2181BvlmuHKYtV3dD3zXCOi3C7r7OyEiAy/Nm7cGKlqVnrk7josmuHFn/30XZzp6XM7HCIa42wnGTOxPCUiXTD2lGkH0KCqS1X1aVUNt2pzKD/MGWAWocdO6ju93lYYqxe0RapQVVUFVR1+5VqSyc/z4LtrFqG3fwj/YfsBDPEhTSJKQNQkY26x/G0ROQLgJfP0ClWdayaWo5a6drZh7sLo3ocXMG5txVHf9vXMac6bVZW3ymKom1yCP//C9XjjyHl87+XDbodDRGNYtGVlWmGsrFyLa2uXPaaq+yN8pDnWLzN7EP6Q0z5EGCOJVd/u9cyHStuCCUZEwk40oGseWDoD9980HX/z0mG8cYTbNRNRfKL1ZLwwxlv2AJgtIvdbXivNV/D9U7D/MOb2kJUEmmBZxVlEakPKo9a3cb1GGImnVUS85kwzPnEYg4jgL790A+oml+DRH+/n+AwRxUVUw99zF5Fvh26/HPVCIk+p6mM2664H0AZzKnHIE/rrATSpapOd+tHKzYH+7jAhtKjq6jDn0dDQoK2trXb+jJxw+PRFfOH7b+D6qlL8w0M3oyg/z+2QiCjDiMg+VW0IWxYpyeQqJpnR/vngx/jG821YvaQam1YthAg3QSWia6Ilmbifk6Hcce/CaXh0xVzs2HcS/+vX3DaIiOyLtkAm0bBHV8zF4TMX8Vc//z1qJ0/A8vlT3A6JiMYA9mTIFo9H8NerF+P6qjJ84x/2o+14uKEuIqKRmGTItnGFefjRv1mKytIifP25vTh8+qLbIRFRhmOSIUcmTyzC//76zSjI8+CrP3obp/y9bodERBmMSYYcm1kxHn//9WW4dHUQf7T1N/j4AhMNEYXHJENxuW5aKZ772jKcu9SPL2/5DTrZoyGiMJhkKG5LZpXj7//tMnQx0RBRBEwylJD6mUai6b7cj9XPvoUjZy65HRIRZRAmGUrYTTPL8Y9rb8HVwSGsfvZN7Of0ZiIyMclQUtwwvQw7H7kNpeMK8JWtv8XLH5x2OyQiygBMMpQ0syomoOXh2zCnsgQP/l0rtrzWDq6NR5TbmGQoqSZPLMK2dbfgczdMw1/9/AP8x+3voG9gyO2wiMglTDKUdOML8/H9r9yE/9Q0D/93/ymsfvYtfHT+stthEZELmGQoJUQE/27FXGz9agM+On8Z9/2PX+OFg51uh0VEacYkQynVtGAKfv7oHZgzpQTffH4/Htt5EBf7BtwOi4jShEmGUq66fDy2r7sVj9xdh22tJ/DZ//46Xjt01u2wiCgNmGQoLQryPNjw2floefg2FBcYi2uub3kHF3rZqyHKZkwylFZLZpXjn//9HXj4rjq07DuJFX/9CrbvPYFAgFOdibIRkwylXXFBHh773Hz89Bu3Y6ZvPNbvPIgv/uAN7PuIKwUQZRsmGXLNjdVlaHn4Nvy3BxbhdE8fVj7zJr7xfBuOnOFmaETZIt/tACi3eTyCL91UjXsWTMWzr7bjR78+il/87mP84eLpeHTFXNRMmuB2iESUAOGyHyM1NDRoa2ur22HkrK7L/dj8ajv+7q1jGBhS3HvjNDx0Ry1urC5zOzQiikBE9qlqQ9gyJpmRmGQyw5mLfdjyagd+vPcELl0dxC21Pjx4ey0+Pb8SeR5xOzwismCScYBJJrP09A1g29sn8LdvHEXnhT5MKyvG6iXVWN0wAzN8490Oj4jAJOMIk0xmGhgKYM/7p/HjvSfw2mHjQc7b50zCyvpqrLiuEhOLC1yOkCh3Mck4wCST+U75e9HSehLbW0/glL8Xhfke3Dl3Mu5dOBWN101hwiFKMyYZB5hkxo5AQLH/RDdeOPgxfvG7T/BJTx8K8zxYNtuHuz81GXfNm4w5lSUQ4RgOUSoxyTjAJDM2BRPOL373CV45dBZHzlwCAEz3jsMdcydhaY0Py2b7UF0+jkmHKMmYZBxgkskOJ7uv4NVDZ/HKh2fx247z6OkbBABMLS3G0tk+LK0px43Ty3DdtFIUF+S5HC3R2MYk4wCTTPYJBBQfnr6Ivce68PbRLuw91oXTPVcBAHkewdzKEtwwvQw3VJXiU1NLMaeyBJNKCtnjIbKJScYBJpnsp6o45e/Fu6cu4N1TPXi38wLePXUB5y71D9cpG1eAuZUlmGO+aiomYIZvPKrLx2FCERfKILLKuCQjIusBdADwAYCqbkmkfqLlVkwyuUlVcbrnKg6dvogjZy7hyNlLOHLmEtrPXML5y/0j6vomFGJG+ThUm0lnWmkxppQWo7K0CJUTizF5YhFvwVFOyagkIyLNAPaqaku4Y6f1Ez0OxSRDobou9+Oj85dxsrsXJ7qv4ERXL052X8HJ7l6c6u5F/1Bg1GfKxhVgSmkRJk8sQvn4QnjHF5g/C1E+/N74WT6+ECXF+VzJgMasTEsy3apabjmuB9Csqk3x1E+0PBSTDDkRCCi6r/TjdM9VnLnYhzMXr+JMj/HzdE8fzl68Cv+VAXRf6ceF3gFE2zZnXEEeSorzUVJkvCYU5aGkqAATi6+9n1CYh+KCPBQX5qE432O8L8hDcYEH4yzvi/JHns/P44LrlDrRkkxaby6bX/Ch/AAa46mfaDlRojweQUVJESpKirAApVHrBgKKnr4BdJtJx3+lH92XB+DvHcClvkFcujqAS1eHcOnqIC5fHcSlvkGc8vca781X/+DoXpOtOMXYnbQwz4OCfA/yPWIc53tQkCfI9xjnC/OM8/l5194bx4LCPA88HkG+R+ARQV7wvUeQZx6PeIlcqz9cB8jzeJDnATxi/N7h93kCEePaAuOcRwAIrp3zmOdg/BSzjkAgwXojfhp1gtezllk/a/091t8Ps3MZnAMS7GsGJ4VcOw6Wy4hjxCiPeb0smHyS7hFMH4CukHOhx07qJ1pOlDYej8Br3jKbjfi2MBgYCqBvYAh9A8bPq4PG+96BoRHn+waG0DcYQF+/8b5/KID+oQAGhxQDQwEMDAXQP2i8Hwxcex98Xe4fwuDwsaJ/0HgfUMVQQDEYUAQCiiHzeCigUXtplByWvBcz0WFUIhtZHi7Rtf1ZEwrzk9vrTXeS8UYqEBGvqvqd1E+0PMzvQ2dn54h/PTzxxBPYuHFjpMsQpVWwZzGx2O1IRtNgwlFFIAAMBgIIBIAh1RHvh4ZGJicjQZmJSxWqxrUURu8voIAieB4IqFnP/J3GOeN9wPrZ4TLzszDiulZm+Swsnx3+Hdf+rpF/p/kTGnIcvdzaTk4+p9c+aLtuaDlGlYf/XCqGBdOdZPwwZ3hZhB47qZ9o+ShVVVXo7OyMVoWIwhDzlte1LxXOsKP0b7/chdG9Cy8AhOtV2KifaDkREaVQWpOMqrbB6F1Y+QDsiad+ouVERJRabsxr3C4iqyzHTQA2Bw9EpDakPGr9JJQTEVGKuPnEfxuAWmDkE/hmWZP1OZZo9ZNRbsXnZIiInIn2nIw5m4Kv4GvJkiUaryeeeCLuz+YitpczbC/n2GbOxNteAFo1wncqF8gMkUhPRkRGTXekyNhezrC9nGObORNve0XryXCtCSIiShkmGSIiShneLgshImcBfBTnx6sA8ElO+9hezrC9nGObORNve81S1cnhCphkiIgoZXi7jIiIUoZJhoiIUoZJhoiIUoZJhoiIUibdS/1nJXPZmg6Y2wholGVrsp25j89aABWquiFMedS2ysW2tLQZACwFsNtpu+RSu5nttcY8rAOA0P/W2F6RichmVV0Xci517RVpKQC+7L0ANANYFek4l14wtrVeBWMB0s1O2ypX2xJAc8hxO4C1bLeI7bUZgNdyvA/AeraXvf/WYPwjJvRcytrL9T96rL8AdIcc14f+j5hrL/M/wnBJJmpb5WJbwtjfaEfIufUA2tluEdtsX8iX3g5rG7K9IrZbfYQkk9L24phMAkSkPsxpP4x/0ZNFrLbK8bZsFJFay7Ef5orhbLfRVHWJqrZYTtUD2A2wvWJogNlOQeloLyaZxPhg7L5pFXpMhlhtlZNtqap+VS1X1Q7L6SZc21iP7RaFOVawR6+NEbC9wjD31Noepijl7cUkkxhvpAJzcJKu8UYqMNsqVnlOMP/WRgDBgWxvjLqxyrOSiHhFJDhZot1S5I32GRvlWcf8u/wafst5b4zPxSqPibPLEuOHOdvCIvSYDH5Eb6tY5bliK4DVamwdDrDdwjK/MLcAgIjsFpGlqroabK9w1mjk2WB+pLi92JNJTBdGZ3ovMPx/AromVlvlfFuat342q+oey2m2m4XZg1kfcno3jFmNANtrBHNMZU+UKilvL/ZkEqCqbSLiDzntQ/T/UXNSrLbK9bY075m3BROMiDSq6h622ygNAJpFZEu4Lzm21yg+GBNLgsdLAdSaibolHe3FnkzitptfEEFNMObx02ix2ion21JEGmH8H7fV/Jd6LYwZU0FsN5OZhDeEJJgmAJssx2wvk/kPlU3BF4xen988Dk42SWl7can/JDD/VdAGc9pplPufWc3smjcCCD5NvBnGzJ82S52obZVrbWkOnnaHKWoxxxiC9dhuJjMJB7/0KgCcN79ArXXYXiHMiRKrYfQGnwQw3BtMZXsxyRARUcrwdhkREaUMkwwREaUMkwwREaUMkwwREaUMkwwREaUMkwwREaUMkwyRDclcPFFEmkUkLQ//ZeuijzR2MMkQ2bMmdhXbtiF9T5gnM24ix5hkiOxpStaFVLXNugpCiiUtbqJ4MMkQxSAizYiyr0amGqtxU3bhKsxEUZgLA3pxbeVawNjHpBbG3i+tMPaYrwXQFFxvzLKgYC1Grq5cC/NWmao2meu9NcPYt+NJGAtlegEsVdXgxmWxYlwLILjYoRfXdjMcFXeYtaq8AHyquiUkluDtvFoAdXZjIRpFVfnii68oL5h7yIc53whgH4wvYi+A9eb59SH1dgOojXQ963Us53YAaLQR29qQz3kBNMeIeweAestxM4BVlljaAXhD4ht1Hb74svPi7TKi+HXB+DLuUFW/XlsJeGnI0uhtML6og/yRrmM51wFzxVsbgqteQ42eyrZIFYPbCOjIMaFtlmt0AehQy1L6avTCas0tCYgc4e0yosR0hJ7QkUv0B3s53hjX8cfzy9W4zbVDRBTGRlI7NPoy7I0A/CEJw4vYCa0NRs8oWzf3ohRhkiFySERqLb0Of7hyABtg3ALbE65OGF1xxuJV1dXm8zANADaIyBJVXRembjDhdejILZ4BoCWe308UC2+XETlXH6N8H4zdG7eE3AJLxcORjwPGbTI1dkFsQuReST0sG085xF4MxYVJhii20PERv+W911rRnKEFHbk9sNcs88KY+RVOuPPeMOdGMWeXWXVYfo6I2+zBdAXjjHCNBmsyNMeX0vlsD2UR7oxJZIM55dcPoEtVW8wv6cdhjHGEbmXbDOA8jF5DF4wv+2YYs8zazPeNMG6ptVqvo6qbzC/14FTiDWFubYXGZf3y98JICB3h4g739wBAsMwyjbnZcr1aDdnemMguJhkiGhZMMuZtN6KE8XYZERGlDJMMEVlFGjMiiguTDBEBGL5Vtg7GwP/6WPWJ7OCYDBERpQx7MkRElDJMMkRElDJMMkRElDJMMkRElDJMMkRElDL/H9Iz3rwT1dbkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "__=ax.plot(loss_history)\n",
    "xlabel = ax.set_xlabel(r'${\\rm train\\ step}$')\n",
    "ylabel = ax.set_ylabel(r'${\\rm MSE\\ loss\\ history}$')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "beautiful-scratch",
   "metadata": {},
   "source": [
    "### Calculate the predictions of the best-fitting model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "stable-chair",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_best = get_params(opt_state)\n",
    "frac_mass_loss_best_fit = predict_frac_mass_loss(p_best, tarr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "meaningful-skating",
   "metadata": {},
   "source": [
    "### Check the quality of the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "strange-coach",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEOCAYAAABIESrBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8kklEQVR4nO3deXhU1fnA8e/JRhaWScIa9gkIKAokYVFRESYo7tVAXKjVVhKttdqqROwitVVIXGqrVhPqVrU1JFKty09NcF/JAq4oyAAuUbYw7EuW8/vj3hkmk8kkuSQzWd7P88yTmXvOvfNCMnlz7tmU1hohhBDCirBQByCEEKLzkiQihBDCMkkiQgghLJMkIoQQwjJJIkIIISyLCHUAba1v3756xIgRoQ5DCCE6lYqKiu1a636tPa/LJZERI0ZQXl4e6jCEEKJTUUpttnKe3M4SQghhmSQRIYQQlkkSEUIIYZkkESGEEJZJEhFCCGGZJBEhhBCWSRIRQghhmSQRIYQQlkkSEUIIYZkkESGEEJZJEhFCCGGZJBEhhBCWSRIRQghhWdBX8VVK2YAsIFFrndOC+gsBJ5AAoLUuaNcAhRBCtFhQWyJKKQfgAJIBWwvq5wJOrXWxmTySlVIZ7RulEEKIlgpqEtFal2qtiwFXC0/JMuu7FQLZbR6YEKJTKCgI/Y2IjhBDR9JhN6VSSqX4OezCaMk0b/MHsG8bRPQwHuE9oEcvGDjeKD+4C8KjICIalGqrsIUQ7cTpdOJyubp9DB1Nh00iGH0g1T7HfF83UlVVRdzY6Zxv+5retdU8fE4MAM99WcPB6P4cc/tqBvaOJrEokx7ffwBhERDV00gwSZMg80njQh8+DLUHoOeAI48EO0TFtvE/U4jQ+NMLn/NF1e6QvPexSb257dzjWnVObm4uycnJ7RQRuFwubDZbSGPojDpyErE1VaCUsmmtXf7KkpKSyH9kOX+46VcccG3nEp1J3eEDvPn+g0AYsQPfB6Dm31/QNy6BJdln0j/qMK+vWkvSjhrOPa+GPjGRsPpJ2PJZw4uPPBV+9oLx/IMHIa4f9BsLfY+ByOi2+DcLIfwoLS3F6XTidDqx2WzY7XYcDuOmRGVlJWC0EkpKSsjNzcVms1FaWkp2djbZ2dnYbDby8/NZuXIlNpuNyspKSktLsdvtlJSUkJ6eTllZGbm5uTidToqLi7Hb7ZSVlbFo0SLP9ZqKoTvryEnEhTkiy4vva7/OnZDEuSUrPK8P1tTx7bU/o2rHbuqje1HlOsDyvT/HdQjurT+db77fj/O/PyN6WAx/qHuNsQN7s/XlJM6ceTG/vTSdRO2C3VXQo7dxwbpaeP0vULPfeB0WCSNOhilZMPbso/6HCxEMrW0JhJLD4fAki6ysrAZlc+fOpaioiIwMY8xNTk4O+fn5OBwOsrOzKSwspKKiAsDT0pg1axYbN27EZrNRVlZGdXU1ubm5AKSnp7NhwwYAUlJSmDt3LiUlJQFj6M46chKppnFrxAbQVCukKdGR4Ywe0o/RQ/p5jl02dbHneV29ZuOvvuLLb7fh3KP4YN0PvL66DOehnjy9084IWyTOJ+/k2ut/w83JmrDwCMjZDDu+hm1roWo1rHsVXN8aF9y3Hd6/H8bMgSGTISy89f96IUSLVFRUeJKD3W7H6XQ2KLfb7UDjX/zet67cSaOgoICUlJQG55aXl7dD1F1Hh00iWutKpZTL53ACUNrW7xUephg1KJ5Rg+IB+PWs0dRkOfns+12Ubarm9fK1rHLt5t7SDbzkepOzj4llhN5KxnlnET7gWBh/Ecz+C9TXGRf8YQ188AC8dx/0GgSn3Agplxsd/EKINpeTk8PkyZOprq6murph1+nkyZMb1Z83bx7FxcWe1kVRURFgJBOXy0Vp6ZFfM+4y4V+HmrGulLL7zANZ7vM6HcgPRiyR4WFMGhZP1qnJPPPbc6jetJb8RT8nyRZN7v3/5OILz+Oye57jnfXbqK01k4e7xTHKAQudkPEYxI+El2+C+9PggCsYoQvR5RUXGyP/XS4XqampLFq0iIyMDNLS0jx1Ao2iSk1NJSUlhfLyckpKSjytEnfCcTgcDR6BYujugj3ZMMWcgZ4BOJRSC32G8mbgNQ9Ea50N2JVSDqVUFrDBZ95I0ERHhnPBpCE8k3Ui7z6xlCv+vIyvDvbip4+sYviMecy78hq01l4n9IHxF8KVL8P8FcbzGJtR9n0F1NeH4p8hRKdlt9vZsWNHg2Pl5eXYbDZPEnDfynJ3gDeloqLCb8d4RkYG1dXVDRKQ97wQfzF0d6rBL74uIC0tTQfrHuah2jpe+exHrv/Nb9nh2sOVC+/gtnOPJbFngNtW1Rvh/lToPw7OuAPsM4ISqxBdwdy5c0lPT2+QALKzs0lNTcVut2O32z3DcFNSUsjJyfHU8e4TycnJ8bQkqqurSUtLIz8/39Onkp+f72mVpKSkePpVmoqhK1BKVWit05qv6XOeJJGjd7i2ngdeX8dDbzmJ3LeVpE2v8Mw/H6Bv376NK9fXwWcr4I07wLUZ5uTBlAVBjVeI7qy4uJjq6uoGSaWyspKcnBxKSkpCGFloWU0iHapPpLOKigjjt7PH8uJ1pxBd/TVvvPZ/XPdUGT/sOtC4clg4nDAXrnkPjjnT6C959XfBD1qIbsrpdDZoWQANRmSJ1pGWSBurq9c8VPIZD7z7HZFhYUzc9S4P33ELcXFxjSvX1xkJpHcSnPzr4AcrRDdVUFCAy+XCbrd7+kCysrKanbHelcntLFOok4jbNzv2c23+K7z4+7nMvfZWlt//56Yra22s3/XtKogfAT37By1OIYQAuZ3V4QxLjOW/OReQ8afH+Ch2Ko++u7HpykpBzUFYfjn8cxZs/TJ4gQohxFGQJNKOIsLD+Pet85kzfhC3Fb7P7LlXcPDgQf+VI6Ph4n8byeSR2cYqxEII0cFJEmlnkeFh/P2SSRyjv6X0+UL+Whhgwv3gFFiwEuL6Gq2SvVuDF6gQQlggSSQIoiLCePm+m8nIXUHBl4oVld81Xdk2DDKfgkO7YdWy4AUphBAWSBIJkh4R4fzrV7M5KTmRa+54mOmzz6WmpsZ/5QHHwlWlMGNRcIMUQohWkiQSRNGR4Sy7PI3BUQcp//xrXqjc1HTlgcdDWBjs+h6q1gQrRCGEaBVJIkEWGxXB248vZXbOw/zxZaf/CYluWht9I89cBvub3dRRCCGCTpJICPTsEcEDl03m8P49zM68ip07d/qvqBScdRfs3QLP/8pIKkKIkKisrCQ1NdWzHpcwSBIJkWGJsWSOieSz0uXkPfZs0xUHp0D67fDVS/BRUFbBF6JDmjt3boMVda3U93d+S6+bkpJCdnZ2s/W6mw67KVV3cNuV5/LRliJe2xfHwn2HiY+L8l9x2jWw6R0o+QMMPxEGTQhuoEJ0ANnZ2Y3WvGpNfafT6XePkdZcNyEhwbMLojBISySEIsLD+NvPZ7LrQA3X/+N59u/f77+iUnD+g5B6pbHJlRDdkMPhaFUS8a3v3kP9aK9rRaANsjq7VrdElFILADtQqLVeo5R6yHztBHK11pvaNsSubdyg3mSOjebOK88nbPP1PP7A3f4rxibAWXnBDU50fY+d3fjYcRcY2xMc3g9Pz21cPvFSmHQZ7NthDPzwNfnnxpbRu76DFX5u/5z0Kxgzp1VhVlZWsmDBAhwOB7m5uZSWlpKTk0NmZqYnAZSUlJCTk4Pdbvdb371Rlc1m8+wF4lvP/V5gtFxKSkrIzc1t1cKMlZWVlJaWYrfbKSkpIT09nbKyMtLT08nJyfHsXeJeft5ut5Ofn+95z+LiYux2O2VlZSxatMjz3gUFBdjtdlwuFyUlJZ5zmjoeLFZuZ1UDy7XWu9wJRWt9BoBS6kJgUxvG1y384ZLTWPHyjaxNnMa+Q7XE9QjwbfnxU3hlEVz0CPQaELwghQghd3+E+1aSw+EgOzuboqIizx4g1dXV5Ofnk5ub67e+Ozl47yPiWw+MPpKioiIyMoyduXNyclr1i3nWrFls3LgRm81GWVkZ1dXVngSVnZ1NRUVFg/cuLCz0nJuenu6JJSUlhblz51JSUkJeXh4Oh8OzZL1718amjgeTlSTi0lrvMp9nAN5txF1+6otmREeG88SShczL/4C7X/2K3589lvDwcP+VI2Nh8/vwzj3SMhFH78qXmi6Lig1cHpcYuLzPkMDlreSvP8K7hZCQkNBgU6mW9l/41quoqPBc173TYWt5x9XSPpSCgoIG+5rY7XbcK5Lb7XYWLFhAdnY28+bN8yTCpo4Hk5U+Ee9xpulAeRNlohWmjEwgc1J/lt7wU65ftLjpionJMGk+VDwGrm+CFp8QHVFCQkK7XNe9fW55eTnV1a2bozVv3jyKi4txuVxUVlayaFHLVp7YsGEDLpeL0tJSz6OoqAgw9n5ftGgRRUVFxMfHe4YZN3U8mKwkkWSl1Ail1FKgWGu9G0ApdVXbhtb9/P78ifRKGMAbmw5wuLa+6YqnLTS+vuW/o1AIEZh7f3VfLpeL1NRUFi1aREZGBmlpaQ3KWiI1NZWUlBTKy8spKSkJ2J/inaDce7o7HI4GD4DS0lIyMjIoKSlBa015eTlOp7PJ48HU6iSitV6G0QL5Wms9TynVRym1BEgGZI/Jo9ArOpL/PPk4++yn80ig/Uf6DIG0X8Caf8P29cELUIgQcu9AaLW+3W5nx44dAeuVl5djs9k8v/jdv5DdnfItiaGiosLTce/LvZOim/ftt4yMjEbXd89fKSkpaZAc3Ndu6ngwWR2dlQw8Yx7KxRidtQGQ2XBHada4AZw6KpF7H32G05MyGXvMKP8VT7kR4odDn6HBDVCIEKisrKSoqMjz13dCQkKj1/n5+Q1GN3mXOxwOMjIyKCws9Ixm8nddh8NBWlqap05KSgppaWkUFxfjcDga1fXHZrORnJwMGAnKPRrLnViKioo8LSH3iK2CggKysrIoKipiyZIlDVolAMnJyZ4RXy6Xi8mTJ2O325s8Hkyt3h5XKXURUOo1OivDe3SW1npFO8TZYh1le9yj8dJHaznn5BM455Jf8MKTD4c6HCFECxUXF1NdXd2gg9s9lNe71dERBXN7XBmd1c7OmjKWE6/7O3uOz6C+vpkk/2kxvPib4AQmhAjI6XQ2agl4j7jqiqwM8fUdnTW3iTJhkVKKnJ+dx/XPrOH1L7fiODbAfJCdm6D8UZhwKQydHLQYhRCNLVy4kIKCAvLy8jz9Hy6XyzPKqiuycjtrAVACXI0x0XCeefwqwKm1fr0F11iIMcM9AUBrHXD1M7O+y3xp01o3OUGiK9zOAqipq2di9j1se+cZNq95l5iYGP8VD+2Fv0+E/uPgZy8ENUYhRNcRtNtZfkZn9TaH+46iBaOzlFK5GMmm2EweyUqpjAD1F2qt87TWBWb9UjOpdGmR4WGcNWEIrp3VvLbq86Yr9ugJp9wEG98G55tBi08IIcDiAoxmIlHmulnLMBLKLVrrJhZ+aiBLa+09SLsQCLS+cqbPe1cC3eK+zR8XZHDMNQ/x2nfNfJvSroTeQ2Dln2XPESFEUFlKIkqp1zBaHU6MGetpSqkypVTvZs7z11JxAYEGN1crpTw3FJVSWRiJp8vrFR3J/GkjeHnNN3zw8VdNV4zoAWffAyf/WpKIECKorPSJXAUUeY3Qch+3AVcFao0opRxAvtY62ee8nVpr1cQ5dqDCfLkE81ZYU+/RVfpE3LbsPsiIYyfRv08smz+vaP4EIYSwIJhDfHf6JhAArbULCDDNGgBbUwVmMmlEa+3ESB7lGMOJA97KqqqqQinleSxevLiZkDq2Ab2jmX1JFnrChVTvPRS48v5qeDPXWIJbCCGC4GgXYGxNGRi3rnxXTAu4gppSKh9jcmM6Rod+lvftLV9JSUlorT2Pzp5EAPJuvIqwYRP596pmFlw8tAfeWmoM+RVCiCCwkkQSlVIjfA8qpSZiLIcSSDWNWyM28LRkfK+ZgjG5sdKsUwqMJHAfSpczZmAvThoWw9KluXy57uumK8YPh2POhIonoLaZVosQQrQBq0N888yO9ELzUQYs0lrf1cy5lRyZ7+GWAJQ2cUoC0GDFNDPZNFW/y8o8oR/fr3yCpQVPB644ZQHs3w6fPxeUuIQQ3ZvVIb7zgCyMX+blGMN2MwOf5bHcZ15IOl4LNyql7O5ys+WR7n2y2XcS/O27Quy8k8dz+h//w7eDTgu8FMrIGZA4GlYFnL8pRKdSWlpKamoqeXmdbyO2yspK8vLyKC4uJi8vr8vtt24piQBorVdrrZdpre/SWq8GUErd1ILzsgG7UsphDtfd4DPaKoOG80aylVK5Sqkss/48rXXwd14JMaUU159/Ihu27eOtdduarhgWBlOzoecAY49sIboAh8NBZmZL/05tHfdy6+3B6XSSk5PDwoULycjIYMOGDSxfvjxo7x8MAdfOUkq92oprKSAVaHbCYaBlS8yyPK/XTqDbJQ1/zjp+ENfc/EeufHUpm1a/03TFKQuMhxAiIKfT2a4tg9LSUtLTj9xMyc3NbbBJVXu/fzA01xJRwC0Y62S15LGy3SIVRIaHMX54f3bWRuLae6D5E3ZsgAM72z8w0WnNmDGDxx9/HICamhpmzJjBU089BcD+/fuZMWMGhYXG3N5du3YxY8YMVqwwdnvYvn07M2bM4IUXjDXbfvzxR2bMmMErr7wCwLfffsuMGTMoLTW6MJ1OJzNmzOCtt94K5j8xoNzc9t0d1DdB+O5y2N7vHwzNreKb475V1RJKKWkxtLO//O5m5uV/wFtfV3P+xMFNV3R9A/engmMxTL8hWOEJ0a42bNhAcXExNpuNyspKMjIyPEuve29IVVZWxqJFizy/tN2bTLlcLkpKSsjPz6e0tNSzY6HNZmtyN0KrSktLPXuI2Gw2XC4XhYWFOBwOcnNz2/39g8Z7TkVXeKSmpuqurK6uXp94Z6m++G+vNF/5sbO1vne81nW17R+YEO0sNzdXZ2RkNDhmt9v9Pt+wYYN2OBye8yoqKhpcx/u59+u25nv9/Px8vXDhwqC9f2sA5drC71zLHesiNMLCFPb9X/DM9XN496PKwJWnLIBd38C61nRtCdFx+W74ZLfbKSgooKCgoMHmT3a7HffyR3a7nQULFlBQUIDL5Wqw66A4elY2pRIhlj33LF58/QNWb61leqCKY86GXklQtgzGnhWs8IQIGrvdzoYNGwCj/8Hd/wJ4NoLKyDBmFOTn55OdnU1WVhb5+fmNL9aEgoICKioCr1uXk5MT9L3NOwpJIp3Q9OPtTL4om3e+r+O6QBXDIyDt5/DO3bB3G/TsF6wQhQiK6upqUlNTSUhIoLKy0m+fQmlpKRkZGZ5kkpqa6ncb2+LiYk8db8FquTT1/h2d3M7qhJRSnDdhEG+/+x7vVn4WuPKUBfCbzyWBiC7B6XQ2ep2VlUVGRoZnK1o39/yLkpKSBud5Jxq73c6OHQ0WxQiqUL9/W5CWSCd1ur0Xv37mVm47+Dkrix9rumKM7chzrUH5XXFfiA7PPXqpuNiYm1xWVtZg7/KioiKWLFnC5MnGQt/uZJGcnExpaalndNbkyZM9rZCMjAwKCws9o7faUmlpqWd4dEpKCgkJCRQVFeF0OklPT8fhcLTr+wdLq/cT6ei62n4igUy/4X5iBo2mJOfMwBX3V8N/LoaUy2HS/OAEJ4ToVIK2n4hSaoFSaom5ai9KqYeUUq+aX0e09nrCuqvmnsP6nXWs37IncMWYeDjggvIALRYhhLDASp9INbBUa71GKbUAsGutz9BaX4OxZa4IkrNPSOLAV+9y85+WBq6oFKReAd+Xw4+fBiU2IUT3YCWJuPSRnQ0zMHYbdGu046FoP/169SBuy8eUPr+c+vr6wJUnXAzhPaDi8aDEJoToHo52Z8N0jKXg/ZWJILjtzrtJnP9XPv6umfwdmwDHXQCfLIfD+4ISmxCi67OSRJKVUiOUUkuBYq31bgCl1FVtG5poiZ9MG02PyHCeX1PVfOVp18CsP2KsqymEEEfP6s6G6cDXWut5Sqk+SqklGFvjSp9IkPWOjmTUwXXkZZ+Ha9fuwJWTJhnzRqJigxOcEKLLszQ6CyNhuG9j5WLsI2IDips4TbSjWRNGoHv05JXyL5uvfHgflD0CW1tQVwghmmF1dNYSr9FZI7XWs2V0Vuhce/E5jPpZHuXVPZqvXHsIXlkEZf9s/8CEEF2ejM7qAqIjwzlz/EBeXvMNrj3NdJp7OtgLpYNdCHHUZHRWFzGxz0G+vGsueQ8/0Xzl1Cvg0G74/L/tHpcQomuT0VldxIUzUkhIO4cfVGLzlYedCH3HyJwRIcRRO5rRWRu8RmfdB6QhfSIhExMVSea1t/LZoUTq65tpELpnsKPg0N5ghCeE6KIsLQWvtV5mJhO01ru01jcAOYAz4ImiXTmO7U/VN5t45aMWLG0y9Wq4qgR69Gz/wIQQXZalpeDNhRYdGMN6vSUDK44uJGHVtGG9+OHRa1myJZOznvtX4Mph5t8P+6shMsZ4CCFEK1mZJzILKABGmY++5mMURmtEhMigxD6cdNXthJ1wdstO2P413DMWPnu2fQMTQnRZVloik7TWswGUUiMBtNYbzdcTgTVtFZxovSsunctfXlrLt9X7GZrQzMz0xGSIHwGrlsHEy2TDKiFEq1npE9nofmImj1leZQlHHZE4KunHDuCAs4L7Hi9qvrJSMO1q+GENbH6v3WMTQnQ9lvdYV0pdaD5NU0r1Mp+3aHSWUmqhUipDKZWllMpqQX2bUirXXV8pJaPAmjA8MY4DHzzNUwX3t+yECZdAbCK838L6QgjhxcoQ32eVUjcDmeahPGCzUmoH0OwkBaVULuDUWhdrrQsw5p1kBKhvA4q01jlmfYBFrY27O7n6T/cTe97v2bW/pvnKkTEwJQvWvwa7W7ASsBBCeLE6xPcurXWm+dyptU4AHFrrlvxyz9Jaey/UWAhkB6i/DMj3er0c6cAPKHPWZHRYJG+u29qyE6ZkwbVl0DupfQMTQnQ5lm9n+dJar25uj/UmbkO5MIYLNyUDKFVK2ZVSKVprl9Za5qMEMHGoDbX+Te684y8tOyE2AfqOMp5rWblGCNFyR9Mn0tv3QfMthASMVYC9+b72fg930knzOlZk3uISTQgPU/TZs4lPP3iDQzV1LTuprhYK58ObzezXLoQQXqzME7lIKVUNVACV5sP9vLlOcluA6/ors7ufmLfNKjFufy1r6jpVVVUopTyPxYsXNxNS13RH3j30n38P5Zt3tuyE8Aior4dV+XB4f/sGJ4ToMqzME7GbfSCNmB3ugbhoPAw40LBgl/nVe6VgJ8YtLr+SkpKoqpIO4pnHDqZHxOeUfLGFk0f1bdlJJ10HX70Ea542dkAUQohmWLmdVdlUgdb6rmbOraZxa8RmnuvyU9/pp8wFTbZchCkmKpx+37/DvdfNo76+vmUnDZsGg9Pgw39AfQtvgwkhurU261gHUErNDFRu3o5y+RxOAEqbqO8EXD4Jw4axMZbvdYSPSfaB1Eb3oXz99y07QSk46VdQ7YSvXm7f4IQQXUKrb2dprVcqpZaYL3dwJCnYMOaOTG7mEsuVUhlew3zT8RrCq5SyAyle5UuAeRjrdWG+h/v9RQB/uCGLF/Yl89F3B5kypoUnjT0XHH+CodPaNTYhRNdgpWN9KcZqvYojCy+6F2K0NXe+1jobsCulHOZs9Q0+80Yy8Jo3orXOA2zmLPeFwA7zmGhG/17RTBxq45U1m1p+UngETL8BevZrr7CEEF2IlY71Mq31Lf4KlFJlLblAoCRgluX5OSYs6PP9h/wr7xY+OXcdJxwzsuUnri+B78rg9FvbLzghRKdnpU/E1VSB1lrWFO9gLp4zg95pF/DO+u2tO3Hz+/BWntE/IoQQTbCSRHY0NTNdKXXT0YUj2tqZ01M44cJrKN/aypnoU7MhPFImHwohArKSRG4FipRS65VSZUqpV81HObIwYoejlGLW2P68/t4qduza0/ITew2EE38FnxTC5g/aL0AhRKdmJYnYgFswRkxlmc9vMZ+vbLPIRJtJ3LWebx65jgefeq51J556E/QeAi/fbCyLIoQQPqx0rOdorVf7K/Aa+is6kCsvOpPc4t+yJXpo606MioOz74Y9P4Jq0ylFQoguwso8Eb8JpLkyETo9Y2P4ycXzef/r7dTVa8LDWrEN7pg57ReYEKLTkz8vu4lTRvTim/KVvPBOk6vWBFbxBLz6u7YNSgjR6UkS6SZSBsWw/X95PPzYU9YusONr+OAB+K68+bpCiG5Dkkg3MWrEYM74/ePUHX+utQucthB6DYKXbpTFGYUQHpJEupG5Z5zC+m0H+GaHhf1CevSC2X+BH9ZA5RNtHpsQonOSJNKNnJpsw/Xef7j7n/+2doHxF8Hw6bDydjjUijknQoguy8oQX9FJjRpo4/AXr1MSqyHnqtZfQCk4+x7YU2W0TIQQ3V5b7yciy550YEopFj3yMjWpl7HrQI21i/QfC8nmtjH7q9suOCFEpxSwJaKUerUV11JAKnD3UUUk2tVZk4bxyIff89a6bZw3Icn6hdb8G15ZBD9/BfqPa7sAhRCdSnO3sxSQQ4CVe33kHlU0ot1NHBrPvtcf5va1/Tjv+UetX2jEdIiIhqcugqtKofdRJCQhRKfV3O2sHK31aq31xkAPjL3TwUg4ogMLD1MM6h3F5m27qalr4d7r/tiGwWVFcHA3PJUBB3e1XZBCiE4jYBJp6TImWutdGIlkUlsEJdpX3n0PEnfqlZRtOso+jUEnQOa/YPtXUDhfFmkUohuyNDrL3E/EQePtcJOBFUcXkmhvp4zuS1REGK9+8h0nJfc9uoslz4TzH4QDO42tdYUQ3YqVPdZnAQU03Fvdvde63M7qBOJ6RMA7+fz1+kvQupWbVfkz4WKYdo3xfPcPR389IUSnYeVPx0la69kASqmRAGa/CEqpicCatgpOtJ+ZM06jkD58vXUPowf0bpuLbvkCHkmHWbfB1Ky2uaYQokOzMk9ko/uJmTxmeZUlHHVEIihu++3V9DlxHqVrt7XdRfuNgZGnwv/dDC8vhJoDbXdtIUSHZHmyoVLqQvNpmlLKPX055ehDEsEwqE8Mxw6MZcXrH7bdRcPCIeMxmHoNrMqH/NOgak3bXV8I0eG0OolorZ9VSt0MZJqH8oDNSqkdQGJbBifa14F3n2Tlkl/w3TZX2100MhrmLIWf/hcO7YbPZZyFEF2ZpeE0Wuu7vJ47gQSl1CTZ2bBz+VX2z/msdgBvr9/Opf1sbXvx5JlwzfvGFrsA31dAbCLEj2jb9xFChFSbrZ2ltV6tlJrZVtcT7e9Cx8mMmurgza9d7fMGsQkQ0QO0huevg4emw6plcNjCUvRCiA7JUhJRSs1USl3l87gJWfakU1FKccrQaJ5f/jSufQfb843g0mdg0AR4+Sa4d6yx7lb1xubPFUJ0aFbmiSwFrubIPBHv+SK2Fl5joVIqQymVpZRq1VhQpVR+6yIWgSTuWc/Wl+7jgWdeat83sg2DK16EK16C5FmwqgC2fG6UHd4ns92F6KSs9ImUaa1v8VeglCpr7mSlVK55jWL3a6VUhvt1C861tzZg0bTrfjaP/6w9zHo1pP3fTClj4cYR02HPFqOPBOC9v0Hlk3Ds+TD8RBg6DXoNaP94hBBHzcrtLFdTBVrrZ1twfpZPwigEsps7SSklw4fbQWxsLBfPOZX3vt5B9b7DwXvjXgOOLJMybBoMOA4qHofll8M9xxjDg92z6fftgPqjWCxSCNFurLREdiilRmitN/kWKKVu0lo3uZ9IE4nAhbEOV3PSgBJkLkqbO21EDEsWP8Jfh+7lz9fMC34AyTONR+1h+PET2Py+MVFRKaP8iXNgx9fGLbH4EcZj2IlwfIZRvudH6NEbomKDH7sQ3ZyVJHIrMFIpZcNIAO6lYBOBkQTelCrBq75bs0vJKqUygOUYiUS0sUkjB3Dg89d56Y0hoUkibhFRMCTNeLhpDSdeC9vXw85NxuO7cji090gS+XsK1Owz9jeJSTBGhZ0wD06+3mjBvPRbiIyFyBgj0UTGwpDJxvvUHoKvXobwKOMRFmE8EkZCnyFG+dYvQIUbkylVmPHo2R9i4o3y3VVmwlNHvsYmGMObaw8d2QHSnRRREN3HmFNTe8hYTt/NXadHb+P/o/aQ8W/1Fd0bwiOh5iDU+Bnt1qO30dJrqjy6j/HvqTngf2WBaBuEhRkj6Wr9DLqIiTdiPbzPiNFXrLl4xaG9UOfTwlXKOB/g0B6o89llU4VBjM14fnA31Pv0l4WFG/GDsQVBfZ1PeYTx/9NUeXjkke2dD7hA1wco33mkRewpj4IePY3n/nb3jOhhfO+1Ns5vVB5t/BzW18NBV+PyyBjjUV/nf4uFyFjjZ6eu1piL5Ssqzoihrsb4/21U3tP82ToMh82fLWV9oK6VJGIDbqHxL39lHm/uXL+UUjattcvfccCltXYpz4ewaVVVVXjXu+2221i8eHGz53Vn0dHR3FH4Dv94dzNbdx+kf+/oUId0hFIwaX7j4+5fPFrDGXcYH9YD1bDf/Oqen1J7EL580fhFeXgfYP5COOUmI4kc3A1FVzS+vmMxTP8N7P4eCmY0Lj/rbpiyALZ9BfmnNC6/4GGYeAl8XwmPndm4PPMpGHcubHwbns5oXP7T5yD5dCPB+YvvqpVG/J8Wwf9+1bj8lx8aO05WPAav+PlY3vAZ2IbC+w/AG39pXJ6zyfhF/9ZSo8/K1x92GEmq5I9Q9s+GZRHR8PstxvOXboRPnmlYHpsIC53G8/9ebXx/vNmGww2fGM8L58PGtxqWDxgP17xnPH/yJ8YcJG9Dp8EvzE1ZH5kN275sWD7KAfPNO+8PnQy7v2tYfuz5MO9fxvO/TWj8i3zifLjgQeP5XaNA+ySpqVfDnFzjZy9vJI2cchPM+oPxc3pXcuNy98+e6xv4+8TG5Z6fvbXw8PTG5e6fve/Km/nZe+vIz17vwY3rtZCVJJLT1KRCpdSSZs510Xh9rebW25qntS5oYWwkJSVRVVXV0urC9JPJw3nwnc28+EkVP5/eCcYuhEcaX5WCtCubrhcVCzd/bTzX2viruWb/kfNj4o1fuHWHjcRUd9j4y9c9KbLnALjkGeOvwvpaQBt/uQ6aaJT3GWJ8aHW9WaaNr0OnGOUJdjjnPjzJy/1X7YDxxtd+Y4xfCr76jja+DjwB5tzlU6ihz1Dj6ZDJMCfPq0gfiRtg+Mlwpp+R9+6/9JNnHvmr21tEjPF1zFnQy8+ule6/XI89HxJHNywLCz/y/IR5kOSzzVCk1x8pk34KI3ySsHc8UxYYMXiL9VoY48RrYa/P+m/egzKm/7Zxa8A29MjzGbeYf1x4SfD6+XcsNv5i99ZvzJHnZy5p3FIZaH5vwyL9/98PNu/IR8b6L3f/7MQm+C8fdqLxtdegwNePH+G/fMBxxte+xxwpj4oDLm9ctwWU1aXAlVJXYeypngCUaK3/2cwp7j6RCq21CnTMp8xlzopHKeXASGLpTb1HWlqaLi8vb/W/R8DQE88ljHo2f9DOw32FEB2OUqpCa93qLgOrkw1fw+ifcALlGIswlimlAq4prrWupPHorgSgtIlTEoAMc17JQoxRXHbzdSf4c7lzOX70cFyqJ99Wy4xyIUTLtPp2ltkCmWtuiet93AZkEbhjHWC5z7yQdMAzgdBMDila62KtdSleCcacmGjXWuch2lz+fXmckvcGL336A1ef5uderRBC+LDSEtnpm0AAzE7xZtex0Fq7WxMOMyls8Jk3koGfeSNm3bkcaYnYLMQuAhiaEMvEoTaWv/lxqEMRQnQSVjrWA3WitKiDJVBLwixrVG52rre4g11Yoz7+L28+dj+r529kUrKfDlUhhPBipSWSqJQa4XvQ3BpX7oF0ctf9bB7xp1/Ja19sCXUoQohOoNUtEa31MqXUcnN/dXOwN3bAqbXODHCq6ARmnjyF9Iw6Xlu3i4Va05K5OUKI7svS6Cyt9TyMTvRSjNFZWZJAuo4zxyXy6fuv8+Hab0IdihCig7O0syEYm1ABspNhFzSUHWxb8WfuHR5N0b23hjocIUQH1mY7G4Jn+K/o5GZOn0b6jQ/wfd80rE5GFUJ0D822RJRSDwH5Wus15utXm6qKMYO92ZnromNTSpF1yfncVPQxa751MWlYfKhDEkJ0UC25neU7J0QBOTSeea6ApW0Qk+gAZo5JZF/Fc9x+/1f8967fhjocIUQH1eq1s5RSI7XWficVBioLFlk7q21orRmUfBwHeg1m0/svER8XFeqQhBDtKJhrZ+3weeORSqmLlFIzQ51ARNtRSvHSqyXEz7mBf6+SUVpCCP+sJJEs7xda641a62e11q8rpS5so7hEB5A6ejCnjO7L4++s43CtbE8rhGisTUdn0fzeIKKTmRj+HZVLMsl//q3mKwshup2WjM7qA8zDWG23D8YCiP7287DjtRqv6Bp+etap5I2dwn/X/MCvLpQZ7EKIhppNIuaKvcuAZUqpXIyVev0lC6e/1X1F59a/fz/ueehRfv/cZ5Rv3snkEdLYFEIc0arbWVrrHIxdDFf7eUgC6aIuTBlMbO1u/rysuPnKQohupdV9IlrrZ5VSC7z3U1dK9ZFO9a4rNiqC+rce4uX7f8embXtCHY4QogNpdRIxlzYpBardx7TWu7TWKySRdF3/+Ns9DL50CU9+9G2oQxFCdCBWdzbcqLW+y0+Z3NLqomaeNJnzT0ulsOxb9hysCXU4QogOwkoSCTTFfaTVQETH95Nxvdj83L3c+ehzoQ5FCNFBWN3ZsNFqveYx21FHJDqsKaMHUffNGopXfkhdvazuK4Sw1rG+DBillKpWSpWZjx1Aqtb67rYPUXQUcXFxPLOyjJrRMymR7XOFEFjf2fAWjMmFS83HTCCvDeMSHdRZE4Yw2BbDw699HOpQhBAdgOVlT7TWLnPNrGe11h8DO2R0VtcXER5G0jcl/O/Wi3j7882hDkcIEWKWtsdVSo0AHDTuA0kGVhxdSKKj+/VPL+C9td+y5NX1nDR2KBHhbb0EmxCis7AyT2QWUACMMh99zccojM2qRBd3yolTKfhbHmu3HuLJD6U1IkR3ZqUlMklrPRuMvUTAWA7efD0RWNNWwYmO6+zjB/FQVAk3/vpazigtJik+NtQhCSFCwMp9CM/GU2bymOVVJqvzdRNKKWYOOMTudR+x8PGSUIcjhAgRyzezvTrR05RSvcznKS08d6FSKkMplaWUymqmrs2sv1ApVdRcfRE8v/3lApb+p5R3t0TwxpdbQx2OECIErC7AeDOQaR7KAzabc0USmzvfXE7eqbUu1loXAMlKqYwApyzSWueZj7lAjiSSjkEpxa/PnEByvzh+fe9THDhcF+qQhBBBZnWeyF1a60zzuVNrnQA4tNaLWnB6ltbae03xQiDbX0WllA1jPoq3fKQDv8OIighj6uE1fLbsRm56cHmowxFCBFmrO9aVUsuBVb6z07XWq1twrr/bXS6M4cJNcSil7Fprp1d938QiQuiPv7maj3/Yzytbe7Fuyx6OGdCr+ZOEEF2ClZZICcZOh40opXo3c24CXkvIm3xfe5gTGuO9EggY2/SWNnVOVVUVSinPY/Hixc2EJI5WdHQ0/7rzRnrFRHHrs2uol3W1hOg2rCSRDUB8E2XN9VXYmiowb10FZNZxEOB2VlJSElprz0OSSHAk9uxBxtADPH/rXO579o1QhyOECBIrSWQeUKSUWq+UelUpVWg+XgOa6xNx0XgYcGuGBS8D5mqtK1txjgiSq+ZMI37AYPLfWMem7ftCHY4QIgisTDZMw2gJ+N6GUsAtzZxbTePWiA2MW1eBTlRKLQTytdZN3soSoTVw4AA+fPctLnrofS7754cUX30Sg2wxoQ5LCNGOArZElFKTlFI3+RxeoLVeqbVe7fOoBJb4u46bWcflcziBAH0cZhwZQKU7gSilAnXEixCy9+vJY1dM5svn/8HU8y+neu+hUIckhGhHzd3OcmC0MLwF2r1wQwvec7nPvJB0jGG7ACil7N7lZsJIAMrNiYd2WjipUYTGhKE2Zo3ty+79h7ji8TL2HaoNdUhCiHbSkj6RfJ/XgYbXNjsJUGudDdiVUg5z0uAGn3kjGZjzRsyO9BIzhp3mYwMwuQVxixBRSlH46MMUPpbP51W7+cWj73OwRiYiCtEVKa2bHo5pLrBYBPQBnBitkpHm80bVMXY3bHbWentKS0vT5eXloQxBeHm0ZA1XX3I+0877Ka8vu12WjReig1JKVWit01p7XsCOdXOBxTSlVB+OjKLKpnHrBIwksrS1AYiubf5px5I/8QTWHezJLSs+Je+iEwgL871DKoTorFo0OktrvQvYBaCUKnEv/e5LKSXLkYgGoqKi+Kj0Bf5aso6/rVxP3X4XuZedQlSEtEiE6AqsLMC4MkCZ3+QixA2O0ZwW8x33/cLB1Kw7WfvD7lCHJIRoA5a2xxWitZRSPHjjZdRvWcdnfSdw3gPvct2Mkfxy5hjpJxGiE5NPrwianj178uQ/7mFlzhnMGtuXhVdexPjzrsK5bW+oQxNCWCRJRARdQlwUf804nlnTp7I/uj9n/f0dHn9voyzcKEQnJElEhERMTAwvPfMYZY8tZurIRG7Oe5ihU+fw5DvrZE6JEJ2IJBERUgN6R/P4lZNxJNWzv/oHfv/iOqYtWckv71vO2s2y5a4QHZ0kERFySin+89Bd7Fi/msLsE5k6vA8Fv88m7YwL+cXjZbz51VZq6+pDHaYQwg8ZnSU6jLCwMKbZE5k6MoFzkv7LG+t38vZ3Ll574Cu2Pr2Q0362kHkXnsfUkQkcO6g34TKqS4iQkyQiOhylFOemz+DcdDhUW8e/XvmQu1aNo1rH8JeX1nLw28+ofukezr3xXs4+/WQGx2mSeoVznH0oYWGSWIQIJkkiokPrERHOgnNOZsE5xhzXH3Yd4MkXNI84J/BdbU9uf/EL9qx5hepXH2DUdY8zZrSdqK1fsn/zJ1x59fUM6W8j7OAeekVpxiWPkNaLEG1MkojoVAb1iWHh/HNYOP8cALbsPsjK9+N5eVwCyaelsLH6AG+Xr+Kb0n+xdmA6KiycnW//i90fFjMy53n69oxmT/nz7Fz7ARf8Lp+4HhFsWvUaO79dz7xf5hDXI4KvV3/A3p3bOPOCefSICGPjl59Sc+gAJ08/hchwxbcbndTX1TB+/HgiwhTbtvxIuIKhQ4cQphS7XDsJU5CYaKxFunfvXpRSxMXFAXDo0CGUUkRFRQFQV1eHUkpaUaJTkiQiOrUBvaO59MyTuPTMk44c/PkUdu/7O1v31bF97yE+mhTJpzMnM+aUUWzfc5h31sWxO7YXG7fvY++hWta9/xG71pfhHHk+ANtfLODgd1/w7G5j14Nt/8vj8JYNDF5grDu67b93UlP9PUm/eBCArcV/om5vNYOu+BsAWwr/gD58gEGX3024UlQ9nQMohl+eR5iCjY/+hvDoWJIvX4ICvnr4WiJ7JTJq/p9RCj6/P4vovkMYfeltAHx638+JGzyaUfNuRSlYffdP6T1yAqMuugmlFOVLLyF+7DSSL7gegLI759L3hNOxn/tLAD66/ScMmDyHEXOyUAo++MM5JJ1yEcNnXwnA+7+bw9CZlzFs1k+prz3MB7edx/DZVzDktIupPbSfj26/kBFzFjB4+kXU7NvFqiUXYz/nGgZNO4/Du3dQljef5POuY+CUszi480cq7rmSURf+hoGps9m/7Vsq78tizLwc+k2Ywb4fN7L6/l8y9pLf0Xf8dPZ8t46PH7qeY+ffRsK4aeze/AWfFNzIcVf8hfjRqbicH/PZI7cw/hdLsdknsHNdOZ8/8QdOyL6X3sPGUb32Q7546k9MuOZv9BpyDNs/e4cv/3Mnk677B3EDR7J1zRusK8oj9YZlxPQbwpaK11i/4q+k3fQY0fED+WHVy2x4/n4m5zxFj96JVH3wPM4XH2bqrc8QGdeH794pZtMrj3DiH1cQ3iOGb998hs0lT3DSn/5HWEQk36x8im9ef5rpd/wfAJtefYzv31vBybe/AMDGlwv4sewVTrxtBQAb/vcg2z59m2m/KwTg6+f+RvWXq5hyy9MArCu+h10bP2HyzU8A8GXhUvZVfU3qb/4JwNqn/8yBHVWk/PohAL74120c3lPNxGvvB+Czx26l7tBBJlx9r/Gz88+FABx/VR4AHz/8G8J7xDL+yjsAWPPgdUT1SuDYy/9k7QOIJBHRRfWOi6F3HIzq35Np9tmQOftIYcYJDSvfMhOtNQdq6th7sJbtV6exZ99+4mwJHKqtY8Pse9i1ezdDksdRU6/56vjfsW/vPkZPHE9tneaTgb/m4IEDHHfiGOrrNZ/0vJqamhrGnzyKunrNam38sh4/bRhaQ+WeywmLiOTYCUkAJF1wGVExsRw7fiAaSDz/UmJ69WHcuP5oDb3OuZieCf0Zc0xfAKLOysTWfwijRxmvw87KJD5pJKPtiWg0ek4m/YYfQ/JIY+Ht2jmZDBx1PPYRxuvDZ13K4LGTGDHceH1gzqUMO34Kw4fHU19Xy8E5lzJiYhpDh8dTWxPH4TmXMjJlEkOGx1NzMJqaMy8meeIEkobHc2h/BLVnXsyoiScwaHg8BxPDqD/zYkYfP54BQ23s712PPiOTMePH0W+ojX09h8AZmYw9bgx9h9rYEzuMsDMyOfa4Y0gYbGN39AjCz8jkuHGjiB9kwxVlJ/KMTMaPtWMbYGNn+CiizsjkhGOG07ufjR3qGHqckcmEY4bTK9HGtvoxxJyRyaTRQ4mz2dhaO4643ZlMHDWY2N42thweT899mUyyDyK6Zx8GHhxP74OZTBo5gB6xPRmwfwJ9ajKZNKI/kdEx9Js4iXidycThiURERtF3UgoJYQeZNCyesPAIElLSSIysZdJQGwDxaVPoHxfhed0nbSpb4nt6XveechLbBiR6XveccjI7Bw/2vI6bNp1dI0d6XsdMO5V9O8d5XkefdDoH9rg8r6NOOp2aQweYaL6OPNlBXV0NE8zX4dONn/sTzNdhp55BeGQPxpuvOW0OUTFxHDfUxtst+WD5EXA/kc5I9hMRQojWs7qfiNyEFUIIYZkkESGEEJZJEhFCCGGZJBEhhBCWSRIRQghhmSQRIYQQlkkSEUIIYZkkESGEEJZJEhFCCGGZJBEhhBCWhWTtLKXUQsAJJABorQvasr4QQojgCHpLRCmVCzi11sVmMkhWSmW0VX0hhBDBE4rbWVla62Kv14VAdhvWF0IIESRBTSJKqRQ/h12Aoy3qCyGECK5gt0QSgGqfY76vj6Y+VVVVKKU8j8WLF7c+SiGEEC0S7I51W1MFSimb1tp1lPVJSkqiqqrKYnhCCCFaI9gtERfmCCsvvq+Ppr4QQoggCnYSqaZx68IG4K9VYaG+EEKIIApqEtFaV2K0LrwlAKVtUV8IIURwhWKI73KfeR7pQL77hVLK7lMesL4QQojQCXoS0VpnA3allEMplQVs8JkHkoHXPJAW1PdLRmV1XvK967zke9epJVk5SWmt2zqQkEpLS9Pl5eUopehq/7buQr53nZd87zov83unWnueLMAohBDCMkkiQgghLOtyt7OUUtuAzRj392TWYeck37vOS753ndc4rXVca0/qcklECCFE8MjtLCGEEJZJEhFCCGGZJBEhhBCWSRIRQghhWUj2WG9Psh9752SuRpAKFJmH5gK5Wmtn6KIS/iilbEAWkKi1zvFTLp/BDirQ987qZ7BLJRFzP/Yy97IoSqlcpVRGS5ZJER3CPIwf8EpggSSQjkcp5cBYSTu5iXL5DHZQzX3vTK3+DHapJIKxH7t3di0EcgH5Ae4EtNbxoY5BBKa1LgVQSk3G/6Zx8hnsoFrwvbP0GewyfSKyH7sQoSWfwe6pK7VEWr0fu+hYzHuy1ci99M5KPoOdnJXPYFdKIramCpraj110KOWAy30PVilVpJSqlnvpnYqtqQL5DHYKlj6DXeZ2FrIfe6emta706cQrAxaFKh5hiQv5DHZaVj+DXSmJyH7snZg5csSbE/B3j110XPIZ7MSsfga7TBKR/dg7L6WUHSgxx7B7kyG+nYh8Bjuvo/kMdpkkYpL92Dshswmd4/PXaibG0FDRuchnsBM6ms9gl1sK3pwtWwnYQUb4dBbmX0LuXz6JwAb53nU85jBeB5BtHsoHSs1WiLuOfAY7oOa+d1Y/g10uiQghhAiernY7SwghRBBJEhFCCGGZJBEhhBCWSRIRQghhmSQRIYQQlkkSEUIIYZkkESEs8DOz170Bk0ysE91KV1rFV4hgmgf4TsQqDEUgQoSStESEsCbd94C5Cmqlv8pCdFWSRIRoJXMfcVuo4xCiI5DbWUK0grm4oA2wm2tEgXFbKwFzoUGtdbq5TlEuxqq2S8xyGzBZa53jtex2CuD03fjHa/0pG5Ag60+JjkrWzhKildwJQmudHui4mShygbneu8VhJI0cr/N2aq3jvV4XAUu8FsbLBcpkl0fREcntLCHajsvndTVg89ktzt/+DNXu0V7mSqopPn0rhRxZeVWIDkVuZwnRvlx+ju0IUN8BuHx2mbNhLqsuREcjSUSIo6SUsvu0NrxV+znmCnA5G8btLt/dAOVWluiQ5HaWEEevLfeC92zmJERnIElEiNZz0vAXvStA3QQ/x2xNVTZbINVmJ72HUiqrFfEJETSSRIRoJXMf6nylVJZSKkNrXWp2iOcCaebxFGARXkOBzeHBGUCmu8/DLLMDueY1MEd3OdzXN99DhviKDkmG+AohhLBMWiJCCCEskyQihBDCMkkiQgghLJMkIoQQwjJJIkIIISyTJCKEEMIySSJCCCEskyQihBDCMkkiQgghLPt/+3m7GZnJfqwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ylim = ax.set_ylim(-0.05, 1.1)\n",
    "xlim = ax.set_xlim(-0.1, 15)\n",
    "__=ax.plot(tarr, mass_loss_target, label=r'${\\rm target}$')\n",
    "__=ax.plot(tarr, frac_mass_loss_init, '--', label=r'${\\rm initial\\ guess}$')\n",
    "__=ax.plot(tarr, frac_mass_loss_best_fit, ':', color='k', label=r'${\\rm best-fit}$')\n",
    "\n",
    "xlabel = ax.set_xlabel(r'${\\rm time}$')\n",
    "ylabel = ax.set_ylabel(r'${\\rm fractional\\ mass\\ loss}$')\n",
    "leg = ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bibliographic-development",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Target parameters = [3.07721998 1.72745557 0.08644907]\n",
      "Best-fit parameters = [3.075733   1.7490137  0.08698298]\n"
     ]
    }
   ],
   "source": [
    "print(\"Target parameters = {}\".format(p_target))\n",
    "print(\"Best-fit parameters = {}\".format(p_best))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "absolute-period",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ordered-round",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

toy_model.py

"""
"""
from jax import jit as jjit
from jax import numpy as jnp
import numpy as np
from collections import OrderedDict
from jax.experimental import optimizers as jax_opt
from jax import value_and_grad


DEFAULT_PARAMS = OrderedDict(t0=3.0, k=2.0, frac_loss_final=0.2)


def get_default_params(**kwargs):
    pars = OrderedDict(
        [(key, kwargs.get(key, DEFAULT_PARAMS[key])) for key in DEFAULT_PARAMS.keys()]
    )
    p_init = np.array(list(pars.values()))
    return p_init


def get_adam_opt_funcs(step_size=1e-3, **kwargs):
    """Retrieve the three functions used by the JAX implementation of Adam

    Parameters
    ----------
    step_size : float, optional
        Step size parameter defining the Adam configuration.
        Default is 0.01

    **kwargs : floats, optional
        All parameters of the DEFAULT_PARAMS dictionary are accepted

    Returns
    -------
    opt_state : state of the optimizer
        Used to carry the parameters from one step to the next

    opt_update : update function
        Used by the train_step function to take the next step down the gradient

    get_params : function to retrieve the parameter array
        Operates on opt_state, so that get_params(opt_state) returns the parameter array

    """
    opt_init, opt_update, get_params = jax_opt.adam(step_size)
    p_init = get_default_params(**kwargs)
    opt_state = opt_init(p_init)
    return opt_state, opt_update, get_params


@jjit
def jax_sigmoid(x, x0, k, ylo, yhi):
    return ylo + (yhi - ylo) / (1 + jnp.exp(-k * (x - x0)))


@jjit
def predict_frac_mass_loss(params, tarr):
    t0, k, frac_loss_final = params
    frac_mass_loss = jax_sigmoid(tarr, t0, k, 1, frac_loss_final)
    return frac_mass_loss


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


@jjit
def _loss_function(params, loss_data):
    (tarr, target_mass_loss) = loss_data
    pred = predict_frac_mass_loss(params, tarr)
    return _mse(pred, target_mass_loss)


_loss_and_grad_func = jjit(value_and_grad(_loss_function))