aphearin/jax_differentiability_pattern_halo_mass_assembly.ipynb

## jax_differentiability_pattern_halo_mass_assembly.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "unsigned-addiction",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "synthetic-palace",
   "metadata": {},
   "source": [
    "This notebook provides a quick self-contained startup guide to implementing your physical model based on JAX. We'll start by demonstrating the basic syntax of differentiation with JAX, and then show how to build a simple differentiable model for the mass assembly history of a Milky Way-like dark matter halo. The [JAX documentation](https://jax.readthedocs.io/en/latest/jax.html) is quite good, so if this demo is interesting to you, the official docs are the next place to look. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "electoral-moisture",
   "metadata": {},
   "source": [
    "## First look at jit and vmap\n",
    "\n",
    "When writing differentiable code with JAX, the typical pattern is you compose your function so that it behaves on scalar arguments, and then you use `vmap` to create a vectorized map of your function over a selection of its arguments. This next cell shows a basic example that uses `vmap` as well as the `jit` decorator that informs JAX to compile the function for us the first time the function is called."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "sensitive-ethiopia",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax import numpy as jnp\n",
    "from jax import vmap\n",
    "from jax import jit as jjit\n",
    "\n",
    "def my_scalar_function(x, y):\n",
    "    return x*jnp.sin(x) + y*x\n",
    "\n",
    "my_vmapped_function = jjit(vmap(my_scalar_function, in_axes=(0, None)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "color-trail",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEOCAYAAACuOOGFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAzqElEQVR4nO3dd3zV9b3H8dc3gbCTkxAghBUSwLAhBByIqCSKYp0EpE5UiLba1norzor2tgrV1qutZbhwM6y4BYILlRUCyB5JmCFAxsne53v/+J5gjAkh4Zz8zjn5PB+P84Azcs7H4+G8891Ka40QQghxtvysLkAIIYRvkEARQgjhEhIoQgghXEICRQghhEtIoAghhHCJVlYXYJXQ0FAdERFhdRlCCOFVNm3alKW17lLXfS02UCIiIkhOTra6DCGE8CpKqYP13SddXkIIIVxCAkUIIYRLSKAIIYRwCQkUIYQQLiGBIoQQwiUkUIQQQriEBIoQQgiXaLHrUIQQwpM4HJp9JwrZetjO4dxicovLqXJAxzb+dAtsS7+uHRne00ZwhwCrS62XBIoQQljoxyN2liYf4dNtx8gpKgfAT4GtfQB+SlFQWkFZpePU44f3DOKywWFcObQ7fUM7WFV2nSRQhBDCApsO5vJ80l7W7MuiTSs/LhscxvgBXRjVJ5iewe1o7W9GJLTW5JVUsPNYPikHc1m16wR/X7GHv6/Yw7j+odx2fgSXRnfFz09Z/F8EqqWe2BgbG6tl6xUhRHPLLSrnb5/tYummI4R2DCDxoiimjulFYNvWZ/wcx/JKWJp8hHfWHyIzv5TosE78IW4Alw/uhlLuDRal1CatdWyd90mgCCFE8/h+fxZ/WLyF3KJyZlwUyX2X9qN9QNM7iiqqHHzyYwYvrN5PelYRg8MDmTUxmosG1Ll3o0tIoNRBAkUI0Vy01rz0dSrPrtxDZGgHXpwWw6DwQJc9f2WVgw+3ZPD86r0czilhXP9QHrlyIAO7u+41qkmg1EECRQjRHCqqHDy+fDvvbTzMNSPCefr6oWfVKjmdssoq3lx7kBe/3E9+aQWTY3rywGXnEBbU1mWvIYFSBwkUIYS7lZRXcfdbm/hm70nuu7Qff4wf4PYxDgB7cTn/+nI/b6w9iJ8fzBgXSeL4KDq2Ofsgk0CpgwSKa82dO5fs7GymTp1KTk4OS5cuZf78+VaXJYRlSiuquGtRMt+nZvG364YybUzvZq/hcE4xc1fs4eOtGYR2DOD3cQOYNroXrfybvqb9dIEiK+WFyyxYsIAJEyYwf/585syZY3U5QlimtKKKGW+YMPn75OGWhAlAr5D2vDhtJMt/O5bI0I48vnw7lz//Ld/vz3LL68k6lEZ48uMd7MzIt+S1B4UH8sSvBlvy2mfCZrORm5trdRlCWK7KofnDe1tYsy+LuTcMY/KonlaXxIheNhYnnseqncd55vPd5JVUuOV1pIXio5YtW0ZUVBTx8fHY7XYAEhISGDVqFCkpKW573ZSUFNLS0tz2/EJ4uqc/28UXOzJ5bNJApozuZXU5pyiluGxwGCvuv4grhoS55TWkhdIIntxCqG3y5MkALF68GJvNBkBiYiKxsbGnrtc2a9asBp936tSpxMTE1HnfsmXLiIuLIykpSbq9RIu06IcDvPxdOrdfEMGdF/a1upw6tT6L8ZOGyKC8jwsODiY9PR2bzcayZctOBY27RUVFMX/+fOLi4prl9YSw2pe7j3PXomQuje7G/FtG4e8BW6G4g08OyiulfjGFSCn1oFJqslJqplJqphV1eZopU6awYMEC7HY7kZGRbnud2t1oMTExrFq1ym2vJ4QnSc8q4vfvbmFQeCAvTBvhs2HSEK/s8lJKzQEi67hto9Z6WfV1pdTk6ust1axZs4iPjycyMrLB1klTu7xSUlKYMGHCzwbl7XY7UVFRTStaCC9SXF7J3W9uwt9f8Z+bRrlt0aI38Lr/cqVU3R34MFNrXfMbcTEwB2jRgRIZGYnNZiMnJ6fBxzZ1zCMmJuYXP5uWlsaUKVOa9HxCeAutNbPe38beEwUsmj6GXiHtrS7JUl4XKEAssAo4FSz1hIwdkA58zGC8u7/cY2NjmTt3LjabjdTUVJYuXVrv4L8QvmLRDwf4eGsGf7r8HLduyOgtvCpQlFKTgSWYUKkpBKj9K3jDv5K3ECEhIW7/co+Jial39pcQvmjXsXz+9tluLo3uyj3jpXsXvGhQXillA+xaa3sdd9sa+LlfyMjIQCl16jJ79mwXVOk5EhMTSUpKIi0tTb7ohXCx0ooq/vDeFgLbtWbu5GEecbiVJ/CmFsoUrfWCeu6zY1opNdW+/jPh4eFkZGS4oi6PlJCQgN1uJyUlpdmmCgvRUjzz+W72HC/g9emjCe3YxupyPIZXBIpzjCTpNA/J4ZetFBtAPS0anyfrP4Rwj6/3nOD1Hw5w+wURXHxOV6vL8SheESiY1kZcjW2fRwORSqkHgWVa6xSllL2OnzldCAkhRKNkF5bxP0t/5JxunXjoimiry/E4XhEoWuskaoSDc9FipNZ6bo2HLam17iQekP3ThRAuYaYI/0h+aQVv3TWGtq39rS7J43jNoHw1Z5gk4GyhVA+6a60TnbfFOR+T2tIXNQohXOft9YdI2nWChyZGEx3m+qN1fYFXtFBqcg7M1zk4X6vFIoQQLrH/RAH/++lOLhrQhdsviLC6HI/ldS0UIYRoTuWVDn7/3hbaB7TiWZkifFpe10IRQojm9NyqPezIyGfhrbF0DWxrdTkeTVooQghRjx/2Z7Hg2zR+fW5v4gd1s7ocjyeBIoQQdbAXl/PHJVvpG9qBxyYNtLocryCBIlwmLS2NhIQEkpJk+Y/wblprHv7vNrKLynjhxpEtekv6xpB3SbhEdYjIefLCFyzddITPt2fy0BXRDOkRZHU5XkMCRbhE9VYvISGn3UJNCI+XnlXE7I92cF5kCDPGue+UU18kgdIYnz8Emdusee2woXDFM2f88GXLljFr1iwiIyNPnU2SkJBAWloaCxculB2IhaiDmSK8mdb+fvxzass9yrepJFB8VPUOw4sXLz51FkpiYiKxsbH1no3S1COAhfAVz63aw49H8ph38yi6B7WzuhyvI4HSGI1oIXiCyZMnM2PGDOx2Ozab7dSf9WnqEcBC+ILv9mUx/xszRXjikDCry/FKMsvLx02ZMoUFCxZgt9uJjJT+YCHqkl1Yxh+XbCGqSwcenzTI6nK8lrRQfNysWbOIj48nMjKywYO2pMtLtETVuwjbiyt4bfpo2gXILsJNJYHi4yIjI7HZbOTk5DT4WOnyEi3RK9+lk7TrBI9fNYjB4TJF+GxIoLQAiYmJTJkyxa2vkZKSQlJSEsnJycyZM4e0tDRmzpzp1tcU4mxtPJDDM5/v5rJB3bhjbITV5Xg9CZQWICQk5LSD8a4QExNDTEwMDz74oFtfRwhXOVlQxm/fTqFncDuenTKcGifCiiaSQXkflZiYSFJSEmlpaTLmIUQtlVUO7ns3hbySCl66aRSBbVtbXZJPkEDxUQkJCdjtdlJSUmR2lxC1PLtyL+vScvjrdUMZFC6nL7qKdHn5qOqtUIQQP/fB5iPM+yaVaWN6M3lUT6vL8SnSQhFCtBjJB3KYtWwb50WG8OTVg60ux+dIoAghWoRD2cXMfHMTPYLbMe/mUQS0kq8/V5N3VAjh83KKyrlj0UaqHJpXbovF1j7A6pJ8kgSKEMKn5ZdWcNurGzicU8z8W0YR2aWj1SX5LAkUIYTPKimv4q7Xk9l1LJ//3BzDeZGdrS7Jp8ksLyGETyqtqOLutzax8WAOL9w4kkuju1ldks+TQBFC+JzCskruWrSR9ek5PHP9UH41PNzqkloECRQhhE/JLSrn9tc3sv1oHv+cMoJrR/awuqQWw2sCRSllA6p3GxwNrNJaL6j1mAeBNCAEoPb9Qgjftv9EIXct2kiGvZT/3BTDZYPloKzm5DWBAjystT51YIdSKlUpdSo0lFJzgI1a62XV15VSk6uvCyF827d7T/Lbd1Jo08qPd2eey6g+IVaX1OJ4xSwvZ+uk9oZU84GaJ0LNrBUei4FEN5cmhLBYRZWDOV/s5rbXNtDD1o7lvx0rYWIRb2qhxCmlIrXWac7rdpwho5SqaztdOyAbWgnhw/afKOCBpT+y9bCdG0f34s+/GkT7AG/6WvMtXtFC0VrbtdbBNcIEIB5Icv49BKh9JOFpjyjMyMhAKXXqMnv2bNcVLIRwq9KKKp5dsYcr/m8NB7KKeOmmGJ65YZiEicW88t13doHFAROcN9lO91ittb327eHh4WRkZLijPCGEm1RUOXh/0xFeWL2PjLxSro/pwSNXDiS0YxurSxN4aaAAC4EErXWK87od58yuGqQTVQgfUVhWyfubjvDq9+kczC5mRC8b/5g6Qla+exivCxTn1OD5WuukGjfn8MtWig1Md1mzFCaEcKnKKgfr03P45McMPt56jMKySkb0svHnqwZxaXRXObLXA3lVoCilJgMp1WGilIrTWidprVOUUvZaDw/hpzEWIYQXOFFQytrUbL7fn8WXu0+QVVhO+wB/Jg4O49YLIhjRy2Z1ieI0vCZQlFJxOEPCOYYSAsTwU2gsqbXuJB4ztVgI4WEcDs3h3GJ2HStgd2Y+u48VsCszn4PZxQAEtm3FuAFduGpody6J7krb1v4WVyzOhFcEijNAVjmv1gyJU+tOtNaJSqkHncETCaTKokYhrFdQWsHuzAJ2H8tnl/PPPZkFFJVXAaAU9A3twJDwIKaN6c3YqFAGhQfi7yddWt7GKwLFOQ7S4KdLaz3X/dUIIepTWlFFyqFcthy2syMjn50Z+aRnFZ26P6hda6LDOpEQ24uB3TsRHRbIgG6daBfQQlogxTlw8Ac4uQvsh6GixCRqhy4QHAG9xkC3IeDnne+HVwSKEMIzaa358UgeX+4+wdq0bLYcslNe5QCgV0g7BncP4vqRPRjcI5CB3QMJC2zb8gbTS/Phx8Xw4xI4suGn2zt0gYCOoKug4DhUlZnb24XAkBsg5lboPsyamptIAkUI0Shaa7YeyWP55qOs3JFJRl4pfgoGhwdx2wV9OD+qMzG9g+WY3eIcWD/PXErzoMtAuORRiBhngiKgw0+P1RryDsPBtbD3c0h5AzYuhH7xcMnD0GOUdf8djaC01lbXYInY2FidnJxsdRlCeI380gr+u+kI7208zO7MAtq08uOiAV24fHAYE6K7EtyhhQdINUcVbHoNVj9lgiT6Krjwj9AjxnRvnYkSOyS/Amv/bYIp5haIexLaW7+8Tim1SWsdW+d9EihCiNM5nl/Kq9+l8/b6QxSWVTKsZxBTR/fi6uHhdGrb2uryPMuJXbD8N5CRYloiE5+BsCFNf77SfPhmDqz7D3TsCtfNh8jxrqu3CU4XKNLlJYSoU1ZhGS+s3sd7Gw5T6XAwaVg4M8dFMrRnkNWleR6tTYtixaNmXOT6hTA04cxbJPVpGwiX/xWGTYFld8Ib18CF98Olj3nkwL0EihDiZ4rLK1n4bToLvk2ltNLBlNhe3D0+kj6dOzT8wy1RWSEsvwd2fQRRl8K186CTi8+v7z4cEr+BLx6G7/4BJ3bCDS9Dm06ufZ2zJIEihDhl5Y5MZn+0g4y8UiYODuNPE88hqktHq8vyXLkH4d1pZhpw/FNw/n3g56ZN3AM6wNUvmHD57E/wymXw6yVg6+We12sCCRQhBBn2Ep74aAerdh7nnG6dWDptJKMjrB8A9mgHvoclt4CjEm5aBv0mNPwzrjD6TgiJhCW3wWtXwK0fQueo5nntBnjFeShCCPfQWvP2+oPE/eMb1uw7yUNXRPPJ7y6UMGnIzg/hzWuhfWeY8VXzhUm1qEvg9o+hohhenQjHdzTv69dDAkWIFupkQRl3LUrm0Q+2M7K3jVX3j+fu8VG09pevhdPa9DosvR3CR8KdK61rHXQfDtO/AL9W8PokyNhsTR01yCdHiBZo5Y5MLn/+W9bsz+LPVw3izTvOpVdIe6vL8mxaw5rn4OPfQ9QEuGU5tAu2tqYuA+COz83g/BvXQuZ2S8uRQBGiBSmtqOKx5duY+eYmwgLb8ul9F3LHhX3xk40YT8/hMFOCVz8FQ6fAtHchwEMCODgCbvsYWrc33XAn91pWigSKEC3E4ZxiEuat5a11h5h5USTLfzuW/t08a9qpR6qqgA9/A+v+DefebRYX+nvYgs7qUEHBG1dDTpolZUigCNECrN51nKte/I4D2UUsuGUUj1w5kIBW8s+/QRUlsPhm2PouXPKYWfnurmnBZyu0n5nxVVkGi64xuxk3Mw99Z4QQrlBZ5WDuF7u5c1EyPYPb8el947hscJjVZXmHEju8eR3sXQGT/gHj/3T2K9/drdsguOUDs4fYG1dDQWazvrwEihA+6mRBGbe8soGXvk5l2phevH/PBfTu7CH9/p6uINPMnDqSDJNfNWs/vEX4CLh5mdkS/41roCi72V5aAkUIH7TpYC5XvbiGzYdzeTZhOE9fP0yO0T1T2almFXpOOty0BIZcb3VFjddrDPz6Pcg9YAbqS+zN8rISKEL4EK01b649wI0L1tK2tT8f/GYsk0f1tLos75GxxYRJeSHc/onZm8tb9b0Ipr5ldkB+O8HsOeZmEihC+IjSiioeWLqVxz/cwbj+Xfjo3gsZ2D3Q6rK8R9o38PpVZvrtHSvN+SXern+86bI7ugnevdFMMnAjCRQhfMDhnGKuf+kHPth8lPvjBvDyrbEEtfOwqa2ebPPb8NYNZqPFO1eYGVO+YtDVcO1/4MB3sORWqCx320vJ5pBCeLmv95zg9+9tQWvNq7eN5pLorlaX5D0cDlj9JHz/PPQdD1MWWb/63R2GTzX7fn3yB3j/Tpj8Gvi7/utfAkUIL+VwaP711X7+mbSX6LBA5t0c47ozS6oqoCTX/OmoMNuOtAk0Bz552qK+pirNM6cr7v4EYu+AK+b6zn9bXWKnmy6vFQ+bLWQunuXyl2hUoCilIoA4IMp5CQLSADuwEUjSWue7tkQhRG05ReU8sGQLX+05yXUje/C364bSLqCRs7iqKsxBTSd2wcndcHKPWWFdeNyESX0COoKtt1mdHRIJYcPMRomd+3nuor/aMraYDR7th2DiHDg30fPXmLjC+b8xLbDoK93y9GcUKEqpG4CpQDaQBCzBhEgOEALYgEhgrlIqGFistf6vG+oVosVbl5bN79/bTG5RBX+5ZjA3n9cHdSZfhhWlcPB7czm8wQzUVhSb+/xam0Do3A8iLoQOXaF9CPgHmN1slYKyAnPGeXG2+SLOPQCpX0Glc6A3oBP0uQAiLzbnnncd5Hlf0tVH9X7xMLQPhemfQe/zrK6qeY2Y5ranPm2gKKWCgIeBlVrrKfU8LM/552bgfefPTVBKPQ08LS0WIVyjyqF58ct9vLB6HxGdO/Dq7aMZHN7A+e7ZqbB/NexfBelrzJe/8ofuwyDmVug52rQwQvo2rbvHUQVZe83W6Uc2Qvq3sG+Fua9jGERPgoFXQcQ467uTcg+anYLTvoL+l5mjejt0trYmH6O01nXfYcIkQWv9cpOfXKkZmNaKx4VKbGysTk5OtroMIc5IZl4p9y/ewtq0bK4b2YO/XDuEjm3q+H2wvNi0QPatMiFSvUlgcF8zhbRfPESMNcfJuov9MKR/Y7Ys2Z9kWkFtbTBgIgz8lVnb0Zw79VaUwob58PUc02KKfxJG3eE93XMeRim1SWsdW+d99QWKN1JKPYgZ0wkB0FovqO+xEijCG2itWb7lKE98uIOKKs1T1wxm8qieP3VxaQ1Z+8wX9/4kEyaVpdCqrWkV9I+HfnHWHQJVUQKpX8Kuj2HP51Bqh1btzAmH0VfBgMtN15o7VFXA9vfhy79C3iETaFc+61FnsHuj0wVKk2Z5KaUitNYHzqoqF1NKzQE2aq2XVV9XSk2uvi6Et8kqLOPRD7axYsdxYnrbeG7KCPqGdjD7TB1aC2lfm+6sPOeusqEDYNR06B8HfcZC63aW1g+YGqInmUtVhVkLsftT5+UT0/0WMdaEyzlXuubLvigLtrwD6+dD/hHTpXfNv8y4jnCrRrdQlFLzgARgRvXAu1LqrrPpGnMFpVSu1jq4xvUYYI7WOr6ux0sLRXgqrTWfbjvGnz/cgaO0kCfOhWu65+J3ZIMJktwD5oEBncyXZL8J5gTB4D6W1t0oWkNGigmWXZ9A1h5ze3BfE4Z9LoCesWYWWUNjL44qM1vt0DoTUulrQFeZFtr5v4X+l0v3lgu5tMvLOeMrRWudXuO26sH7eVa0XJzhsbpWoEQCqVrrOqeZeFWgaG2mdR7fAbnpUF5k/hG1C4agntBtMISe45aFSi1SZTnkH4W8I1BwzKxXKC80eyFVloLy++ni529mQvkHQKs25svPv43z79W3B9S4rfXP/15RYp7fOXsq91g6Kdu34Zd3mOjWx+nuqLH9ePtQMyOp9/nm0n2Y9QPdrpK1D/athIM/mG676mnLfq1MyAT3MZ/3tkGgHVBVDsU5ZrwmN938/wHz2CHXw5AbzL8L4XKuDpQJQFBd04KVUtdbMV1YKRUHzNdaR9W4zQbk1hco4eHh+tixY6euP/HEE8yePdvNlTZSxmbY9Drs+QIKa3yx+AcACqrKfrqtbZAZ7Bx8PZxzhe980bhbVaX5Tfng95C5zVyy95svrV9QZmwCbe7XDhPsuHYcMld3oiqwJyE9z8EvbIg546LrILPuw9Om4bqDw2F+gcrcZlouWXtNuJfYTfj6+Ztpzu1sZj2MrQ/0GAW9zzV/bwnvkYVcPYaSDCxVSi3ErEdZqrX+0nlfZBNrPFu2+u5QStm01vbat4eHh5ORkeHOmpruwPfm7OrD68xGdf0vMwOr4SPN4Gp133h5MdgPmn946d+YmT07PoAOXWDkLaa53yHU2v8WT1ScY7pGdn9m+vTLC8ztQb0hbCgMutb8RhzUEwJ7mLAO6GD+X9T1ZVVVacK9ssyME1SVmVbOqdvKzaWO23Srdmw4VsXrKbnszfcnZugQ/nRVDF0D2zbrW+JR/PxMiHYbZHUlopGaEijPAHOcf48HFiil+mIWOs5wUV2NZcc5s6sGN00dcaOC4/DZ/8Cuj6BTd7OCd8Q084VWl4D20HWguQybYn5b3r8aUhaZvYnWzzMHA134R/fNpPEWVRUmRDa/bdYhOCrNb7NDJ5txiIiLmr4mwb+VuTRyKu7a1Gye+XwXW4/kER0WzlMJgxjbT34BEN6rKYGSpLVe7fz7auAh5xjGKCDFZZU1Tg6/bKXYAOpqnXiknR/BR/eZPvVLHzeti8bO0vHzhwGXmcvJvfDt32Htv2HzW3DpY2YGkF8LO2SpIBM2vgwpb5gtRQJ7wnm/Mf3s3UdY0j2y+VAuL6zex1d7ThIe1JZnE4Zz3cge+PtJV43wbk0dQ9E1urlq3mfJGIrztWvP8ooDZnn8LC9HFXz5F/jun6Yf+Lr5ENrfdc9/fAd8PgsOrIFuQ+FXz5vZM77Ofgi+/z9IedN0L/W/zLTW+sVZEqpaa77ee5J5X6eyPj2HoHat+c3FUdx2QYScpCi8ikvHULTWq5VSI5VS/6O1ftb5AhOAecAywKo9vJbUWncSD8y3qJYzU1YIS28zC9JGTYcr5pjZP67UbTDc9jHsXA4rHoVX4uH8e+GSR6G1D/bT5x2Br5+Gre8BynQZjv2DZQv7KqscfPLjMeZ9k8ruzAK6B7XlsUkDuXFM77pXugvhxVy2Ut65zUpSzenEzc25Uj4F5+QAj14pX5wD70yBoylw1T9g1O3uf83SPFj5uBljCR0A17wEvUa7/3WbQ1mhaZH88KKZfRU7HS64zwysW6CkvIrFGw+xcE06R+0l9O/akcTxUVw9PJyAVrImQnivJk0bVkoFumIPLlc9j6tZGiiFJ+CNayF7nznoZuBVzfv6+1ebTfLyj5qxmkse897WiqPKrIr+8i9mjGRoAkz4s5lOaoHconIWrT3Aoh8OkFtcwag+wdw9PooJ0V3xkzES4QOa1OWltc53tjpWNWWxonPm1w3V3WLCqTjHhEluOty01Gz13dz6TYB7foBVj5vf6PeuhOv+Y8ZwvEna16Yb7/h26DkGpr5tWYvrSG4xL69JZ/HGw5RUVBE3sCt3j48iNqKFz64TLUqDXV7OUInE7Bq8pcEnVGok5uyU/VZvx3I6lrRQygpMmGT+CL9eAlGXNO/r12X/ajO7rCATLrwfxj/o+nEcVzu514Th3i9MSyTuSRh8nSUztvYdL+A/X6fy4dYMFHD1iHDuHh/FgG6dmr0WIZrDWa+Ud26tMgUz0K0xO/pmY9Z/wM9PcFwFLNRa5/3ymTxHswdKVQW8PdnsMzT1TbNZnqcoscOKR2DL29B1MFw3z2zr4WmKsuGbZ2DjK2bNx7gH4Ny7Lemu+/GInX9/tZ8VO47TrrU/08b05s5xfelh84ANGYVwI1dvvRKEabGEOP/MwQRLsqeHSE3NGiham1bA5jfNQPjIm5rndRtrzxfw8e/MiXwXPQjj/ugZW7hUlMKGBfDts2ZV+6jpcPHD0LFLs5eyIT2HF7/cx5p9WQS2bcXtF0Rw+9i+hHQIaPZahLCCq6cN52FOZwSzsFE05IcXTJhc9CfPDROAcyZCr3Xw+YPw9d/MyvJJ/7BuJpjDAduWmgH3vMPmcKjL/mJ2BmhmOzPymbtiN1/vOUloxwBmTYzm5vN606mtBwSuEB5CJsK7286PYNUTZtPGix+xupqGtQ+BG16GgVebYHklDkbeDBNmN1+LQGsztrN6ttmnrPtwuObflpxncSS3mOdW7mX5lqN0atOKh66I5rbzI2gXIIsRhaitoTPlJwAzMSvODzRLRb7kxC74INGsTL/2Je86k2HQ1WbSwLd/h7UvmWA87x6zbUk7m3te0+GAvZ+b18zYbDZrvP5lsxV5M793ZZVVvLwmnRe/3IfWkHhRFPeMjyKovbRIhKhPQy2USMzuwrbqG5RSI85ktleLV1YAi2+BgI4w5U3POD2vsdp0gvinYMTNptvpmzmwbh6cdzfE3gGdwlzzOiW5sHUxJL9qtisPjoBfvQDDp5mzRJrZd/uy+POH20nLKuKKIWE8ftUgwmWwXYgGNRQoUcAGzKyualOBLe4qyCdoDR/eCzlpcNtHENjd6orOTpcBZmZa5jYTKt/MgTXPmSNbh99o1tI0cqddyoshdbVp+ez6yBxc1WMUXLfAtEgsOCyssKySv3y8k8XJh+nTuT2vTx/Nxed0bfY6hPBWp/1Xq7V+SCn1DLBMKZULJAE2pdQlwCZPXAHvEdb9x+ydFf8URFxodTWuEzYUpr4F2amw6TWzFfyuj8yhU33GQo8Ys4OvrRd06AptOprp0pVlkJ8B9gOQuR0Ob4CjySZE2gWblkjsdDNWYpH1adk8sHQrGfYS7rk4it9P6C+bNgrRSGc8bdi5YDEOSHTeFAnkYrrEUoCNmKOBD7i+TNdz27Thg2th0VUwYKL58vXl0+OqKsyRrXs+h/RvzSl7uur0P+PXCsKGmaNsB0w0QWTh0cUVVQ6eXbmHBd+m0TukPc8lDJfV7UKchkumDWutNwOblVKdtdYPOZ94JBCLOQvlESDG2ZJZAszxlnBxmcITsPR2s3r72pd8O0zArFGJHP/T7KvyYjix07RGik5CeZHzDPUACAw370tIpMeMJ50oKOXetzez4UAO08b05rFJA+kgOwAL0WRN+dfzdPVfqkMGWFh9m3MPrxjMwVvJnrz9iktVVcKyO8yOvje/X/8pi74soL3XnLWy8UAOv307hfzSCp6fOoJrR/awuiQhvF5TFzae7v50pVQOpjvMqjPmm9+XfzGHWF07D8KGWF2NOI031x3kyY920CO4HYvuGMPA7oFWlySET3DX5P5YzAB+y7DrE3OG+6jp5kAn4ZEcDs3fPtvF48u3M65/KB/de6GEiRAu5JYOY+eZ8160iu8sZKfC8nsgfCRMfMbqakQ9SiuqeGDJVj7ddoxbzuvDE78aRCv/lvERFaK5yAjk2SgvMosX/fwhYZH3HlLl4+zF5dy1KJnkg7k8cmU0M8ZFonx9woQQFpBAaSqtzamHJ3bCzcsguI/VFYk6ZBWWcfPL60k7WcS/fx3DpGFevshUCA8mgdJUGxaanXAveQz6xVldjahDZl4pN728jqP2El65PZZx/Zt/u3shWhIJlKY4tB5WPGwW5o17wOpqRB0O5xRz08vryS4sY9H0MZwb2dnqkoTweRIojVV4ApbeBkG94Lr53rWDcAuRnlXETQvXUVhWydszzmNEL5vVJQnRIkigNNbWd82RuXetct827qLJjuQWc9PCdZRWOnh35nkMDm+BC0yFsIgESmNd8Ds4ZxKE9rO6ElHLiYJSbn55PQVllbw7Q8JEiOYm/TWNpZSEiQfKLSrnlpc3cKKgjNenj2FIDwkTIZqbtFCE1ysoreD21zaQnl3Ea7ePZlSfYKtLEqJF8opAUUrZMEcRA4wGVmmtF9R6zIOYg8BCAGrfL3xTSXkVdy5KZntGPvNvHsXYfqFWlyREi+UVgQI8rLWeVX1FKZWqlDoVGkqpOcBGrfWy6utKqcnV14VvKq90cPdbm9h4IIfnp44gblA3q0sSokXz+DEUZ+uk9q7F84FZNa7PrBUei/npIDDhgyqrHPz+vc18s/ckT183lGtGyPbzQljN4wPFKU4pVTNU7DhDRikVU8fj7ZjTJYUPcjg0s97fxufbM3ls0kBuHNPb6pKEEHhBl5fW2g7UHmWN56ft8UOAnFr3174ufITWmic/3sH7KUe4P24Ad41rOUfuCOHpvKWFcoqzCyyOn7q8bA08tk4ZGRkopU5dZs+e7cIqhbv8fcUeFq09yIxxffndBJm+LYQn8fgWSh0WAgla6xTndTvOmV011L7+C+Hh4WRkZLi4NOFO//5qPy99ncq0Mb155MqBsgW9EB7GkkBRSk0GpjbwsByt9c8G1p1Tg+drrWueBpnDL1spNjjVXSZ8wKIfDvD3FXu4ZkQ4/3vtEAkTITyQJYHinJHVqCm9zhBKqQ4TpVSc1jpJa52ilLLXengILekIYh/33oZDPPHRDuIHdePZhOH4+0mYCOGJvGIMRSkVhwmJZKWUzTnjq+bsriXOwKkWj5laLLzcB5uP8PAH2xg/oAv/+vVIWsuxvUJ4LI8fQ3EOrK9yXq0ZEqdaOFrrRKXUg87giQRSZVGj9/ts2zEeWLKV8/p2Zv4to2jTyt/qkoQQp+HxgeIcB2mwj0NrPdf91YjmsnrXcX737mZiegfz8m2xtG0tYSKEp5P+A+Fxvt5zgnveSmFQeCCvTh9NhzYe/3uPEAIJFOFhvtieyYw3kunfrSNv3DGGwLatrS5JCHGGJFCEx/hwy1F++04KQ3sE8c6M87C1D7C6JCFEI0hfgvAIizce4qH/buPcviG8cpt0cwnhjeRfrbCU1poXVu/nn0l7GT+gC/NvGSUD8EJ4KQkUYZnySgcP/3cb76cc4YaYnjx9/VACWkkvrBDeSgJFWCK7sIx739nM2rRs/hg/gPsu7SfbqQjh5SRQRLPbdDCXe99JIbuonOenjuDakXI4lhC+QAJFNJsqh2bhmjSeXbGHcFs7/nvPBQzpEWR1WUIIF5FAEc0i9WQh/7N0K5sP2Zk4OIw5k4cR1E7WmAjhSyRQhFsVlFbwr6/289p3B2gX4M//3TiCq4eHy3iJED5IAsWHnCwoIz2riJyicvJKynFoaO3vR/sAf7oFtqV7UFu6dmpDq2bYsbegtIL3Nhxm/repZBWWc0NMT2ZNPIeugW3d/tpCCGtIoHix4vJKVuzI5KvdJ/khNZuswrIGf6a1vyKqS0fOCevEOWGdGBgWSHT3ToQFtj3rVoPWmp3H8lm++SjvbTxMQWkl50d25tXboxnW03ZWzy2E8HwSKF7ocE4x879N5YOUoxSVVxHaMYCx/UIZ1tNGv64dCe0YgK19AP5KUVHloKC0kuP5pRzLK+VQTjF7jxewMT2HD7f8dASyrX1rosM6ER0WyMDu5s/eIe2xtW9db9A4HJqMvBK2HcljfXoO3+3PYv+JQlr5KS4b3I3Ei6IY3svWTO+KEMJqEiheJK+kgudW7uGd9YfwU4pfDQ/nxjG9GNU7GL8GTjEcFB5Y5/PtySxgT2Y+O48VsDsznyXJhykurzr1mLat/ege1I5ObVvRtpU/SkFJRRUFpZUctZdQXuk49bhRfYKZPjaCK4Z0J6SD7MMlREsjgeIlVu7I5NHl28kuLGPamN7ce2k/uge1O6vnDGrXmjF9QxjTN+TUbQ6H5lBOMbszCzhqL+GYvYRj+aUUlVVSWlGFwwEhHQLoFdKeywZ1IyK0AwO6dWJojyBZ5S5ECyeB4uHKKqt4+rPdvP7DAYb0COS120e7de2Gn58iIrQDEaEd3PYaQgjfJIHiwfJKKkh8M5l1aTnceWFfZk2MllaAEMJjSaB4qAx7CdNf20haVqFsTyKE8AoSKB4ow15Cwry15JVU8Pr0MYztF2p1SUII0SAJFA9zsqCMm19eT35JBe/OOI+hPWWvKyGEd5BA8SB5JRXc8sp6juWV8uadYyRMhBBeRUZ4PURllYN730kh9WQhC2+NJTYipOEfEkIIDyItFA/x1892sWZfFnNuGMqF/WXMRAjhfaSF4gGWbDzMa98f4I6xfZk6urfV5QghRJNIoFhsd2Y+j3+4nQv7hfLIldFWlyOEEE3mlV1eSqn5WuvEWrc9CKQBIQBa6wVW1NYYxeWV3PvOZgLbteb5G0c0y7byQgjhLl73DaaUmgNE1nFbmtZ6mTNIopRSky0psBGe/GgnqScL+eeUEYR2bGN1OUIIcVa8KlCUUjH13DVTa72sxvXFQGI9j/UIK3Zksjj5MPeMj5JBeCGET/CqQAFigVU1b6gnZOxAXHMU1BT24nIe/WA7g7oHcn/8AKvLEUIIl/CaQHF2YS2p464QIKfWbbWve5SnPt6JvbicvycMo7WMmwghfIRXfJsppWyAXWttr+NuWwM/V6eMjAyUUqcus2fPPssqz8zqXcf57+aj/OaSfgwOl5XwQgjf4S2zvKacZtaWHefMrhoaXGYeHh5ORkZGQw9zqbySCh75YBvRYZ2495J+zfraQgjhbpYEirP7amoDD8vRWic6x0iSTvc4ftlKsQHU06KxzD9X7eVkQRkLb42Vc02EED7HkkBxzsha1uADjRAgTqlTZ6aPBiKd606Waa1TlFL2On7mdCHU7HZk5PHG2gPcfF4fhvW0WV2OEEK4nMd3eWmtk6gRDkqpmUCk1npujYctUUpNrjF1OB6Y34xlnpbDoXl8+XaC2wfwQPw5VpcjhBBu4VX9Ls4wScDZQqkedHeumo9USsU5H5Naa12KpZalHCHlkJ2HrxxIUPvWVpcjhBBu4fEtlJqcA/N1Ds7XarF4DHtxOc98vpvYPsFcL8f4CiF8mFe1ULzRsyv3YC8u56lrhuDnpxr+ASGE8FISKG60J7OAd9Yf4tbzIxgUHmh1OUII4VYSKG7018920bFNK/4Q19/qUoQQwu0kUNzkm70n+XbvSX43oT+29gFWlyOEEG4ngeIGVQ7N3z7dRe+Q9txyfh+ryxFCiGYhgeIGS5IPs+d4AQ9fEU2bVv5WlyOEEM1CAsXFCssqeW7lHkZHBDNxSJjV5QghRLORQHGxeV+nklVYzqOTBlFjuxghhPB5EigulGEvYeGaNK4ZEc6IXjaryxFCiGYlgeJCz67Ygwb+dLns1yWEaHkkUFzkxyN2/rv5KHde2Jeewe2tLkcIIZqdBIoLaK3530930blDAL+5OMrqcoQQwhISKC6wcudxNqTncH/8ADq1ld2EhRAtkwTKWSqvdPDM57vp17UjN47uZXU5QghhGQmUs/TWuoOkZxXx6KSBtPKXt1MI0XLJN+BZsBeX83+r9zGufygXD+hidTlCCGEpCZSz8MLq/RSUVvDopIGyiFEI0eJJoDRRelYRb647wNTRvYgOk7NOhBBCAqWJnvl8FwH+ftwfP8DqUoQQwiNIoDTBurRsVuw4zj0XR9G1U1uryxFCCI8ggdJIDofmfz/dSXhQW+4aF2l1OUII4TEkUBpp+ZajbD+az4MTo2nbWs46EUKIahIojdS/ayduOrc3Vw8Pt7oUIYTwKK2sLsDbDO0ZxNCeQ60uQwghPI60UIQQQriEBIoQQgiX8JouL6WUDXgYSHXelKy1Tqlx/4NAGhACoLVe0Nw1CiFES+YVLRRnmCzVWs+qERQP17h/DpCmtV7mvD9KKTXZglKFEKLF8opAARYC82tcXwLMqnF9ptZ6WY3ri4FEdxUze/Zsdz21T5L3q/HkPWsceb8ax13vl9Jau+WJXUkppYFgTHeWrVZXVwywWmsdXOO2SCBVa13vjo2xsbE6OTm5qfXgDe+bp5D3q/HkPWsceb8a52zeL6XUJq11bF33eXwLxRkYALE1blvq7AYDEzI5tX6s9vVfyMjIQCl16iK/4QghxNnxhkH5U/ubaK3TAJRSizHdYAmArb4fVErZtNb2uu4LDw8nIyPDpYUKIURL5g2BYnf+WbN/Kg2YXOP+kFo/U/v6L2zatClLKXWwiTWFK6Ukjc6cvF+NJ+9Z48j71Thn8371qe8OSwLFOQNragMPy9FaJ2LCg1otDbvzeWyY7i1brZ+11fEzP6O1liMWhRDChSwJFOeMrGUNPtA8Nk0pZa/VfWUD7M7rKUope60fCwGSXFOtEEKIM+Hxg/JOTwNTalyf6ryt2pJa607i+fk0YyGEEG7mFdOG4dRK+FO01nPruD8F5yC+rJQXQojm5TWBIoQQwrN5wywvj+KcCDAT6Ky1nlXH/bKnWD2UUjOBUcBS500JwJzq6eBCPj+NIZ+n07Piu8pbxlA8glIqDogDoqhj/YvsKXZGpgCrgDnAfPnH/xP5/DSJfJ7qYNV3lQRKI2itk5wz1Oz1PKRZ9xTzRlrrYK210lqPqrmFjgDk89No8nmqm1XfVRIoLlJji5ia7JjfEoQ4Lfn8iObizs+ajKG4TpP2FGtpnP3eOcgYQW3y+WkC+Tw1ids+axIormOr747T7SnWwiRjFqRW78m2VCmVU6vp3VLZ6rtDPj/1ks9T09jqu+NsP2vS5eU6dpqwp1hLorVOqTVoupEaB6W1cHbk89Mo8nlqMjtu+qy16BZKI/cUa0iT9hTzZo19/5RScVrrmlvipAF19ee2RC3u83O25PPUZG77rLXoQGnMnmJn8Fwtbk+xxrx/zkPPVimlgmt9aGWaJy3z83M25PPUdO78rEmXl2vJnmL1cHZNzKr1j38qZv2AMOTzc4bk83TW3PJZk61XGsE53S6On+ZrzweSah1JLHuK1cP5W2X1h7gz5phmeX9qkM/PmZPPU/2s+q6SQBFCCOES0uUlhBDCJSRQhBBCuIQEihBCCJeQQBFCCOESEihCCCFcQgJFCCGES0igCCGEcAkJFCGEEC4hgSKEEMIlJFCEEEK4RIvebVgIT6KUisNsIx4PzAJiq6+f4REKQlhKWihCeIDqc76dRwIsdV5smM37pjg3QhTCo0mgCOEZah8WFau1Xubcpn1GrZMJhfBIstuwEB5GKTUHsEk3l/A20kIRwvNMxnR5CeFVJFCE8ADOAfnqsZTI6u4vpZRNKTXT0uKEOEMSKEJYzBkY1cevxgL2Gnc/LKcQCm8hYyhCWMw5gysR2IiZ1RWDOZY1DUiRAXnhLSRQhBBCuIR0eQkhhHAJCRQhhBAuIYEihBDCJSRQhBBCuIQEihBCCJeQQBFCCOESEihCCCFcQgJFCCGES0igCCGEcIn/B9hfnCR5Z38SAAAAAElFTkSuQmCC\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",
    "\n",
    "x = np.linspace(-10, 10, 500)\n",
    "__=ax.plot(x, my_vmapped_function(x, 5), label=r'${\\rm y=5}$')\n",
    "__=ax.plot(x, my_vmapped_function(x, 1), label=r'${\\rm y=1}$')\n",
    "leg = ax.legend()\n",
    "xlabel = ax.set_xlabel(r'$x$')\n",
    "ylabel = ax.set_ylabel(r'$f(x)$')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "classical-night",
   "metadata": {},
   "source": [
    "## First look at grad\n",
    "\n",
    "The JAX library has a `grad` function that uses autodiff to calculate the derivative of any JAX-implemented function. The first thing to know about `grad` is that is operates on _scalar-valued_ functions only. So if you want vectorized behavior of your derivative, you just need to write your original function as a scalar, use `grad` to take the derivative, and then use `vmap` to operate on the differentiated function. This next cells show the basic pattern."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "polyphonic-armstrong",
   "metadata": {},
   "outputs": [],
   "source": [
    "from jax import grad \n",
    "\n",
    "my_scalar_function_deriv = grad(my_scalar_function, argnums=(0, ))\n",
    "my_vector_function_deriv = jjit(vmap(my_scalar_function_deriv, in_axes=(0, None)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "static-movement",
   "metadata": {},
   "source": [
    "The function `my_scalar_function_deriv` works fine when passed scalar arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "eight-supervisor",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(DeviceArray(1.4593868, dtype=float32),)\n"
     ]
    }
   ],
   "source": [
    "result = my_scalar_function_deriv(5.0, 1.0)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "unlimited-eligibility",
   "metadata": {},
   "source": [
    "But `my_scalar_function_deriv` fails when passed a vector argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "acute-topic",
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Gradient only defined for scalar-output functions. Output had shape: (500,).",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFilteredStackTrace\u001b[0m                        Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-b4b59260f837>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmy_scalar_function_deriv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mFilteredStackTrace\u001b[0m: TypeError: Gradient only defined for scalar-output functions. Output had shape: (500,).\n\nThe stack trace above excludes JAX-internal frames.\nThe following is the original exception that occurred, unmodified.\n\n--------------------",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-b4b59260f837>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmy_scalar_function_deriv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/_src/traceback_util.py\u001b[0m in \u001b[0;36mreraise_with_filtered_traceback\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    137\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0mreraise_with_filtered_traceback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    138\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 139\u001b[0;31m       \u001b[0;32mreturn\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    140\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mis_under_reraiser\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/api.py\u001b[0m in \u001b[0;36mgrad_f\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    758\u001b[0m   \u001b[0;34m@\u001b[0m\u001b[0mapi_boundary\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    759\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0mgrad_f\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 760\u001b[0;31m     \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mvalue_and_grad_f\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    761\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mg\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    762\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/_src/traceback_util.py\u001b[0m in \u001b[0;36mreraise_with_filtered_traceback\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    137\u001b[0m   \u001b[0;32mdef\u001b[0m \u001b[0mreraise_with_filtered_traceback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    138\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 139\u001b[0;31m       \u001b[0;32mreturn\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    140\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mis_under_reraiser\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0me\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/api.py\u001b[0m in \u001b[0;36mvalue_and_grad_f\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    824\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    825\u001b[0m       \u001b[0mans\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvjp_py\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maux\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_vjp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf_partial\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mdyn_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhas_aux\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 826\u001b[0;31m     \u001b[0m_check_scalar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mans\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    827\u001b[0m     \u001b[0mdtype\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdtypes\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mresult_type\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mans\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    828\u001b[0m     \u001b[0mtree_map\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpartial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_check_output_dtype_grad\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mholomorphic\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mans\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/miniconda3/envs/diffhalos/lib/python3.8/site-packages/jax/api.py\u001b[0m in \u001b[0;36m_check_scalar\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    845\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maval\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mShapedArray\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    846\u001b[0m       \u001b[0;32mif\u001b[0m \u001b[0maval\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 847\u001b[0;31m         \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"had shape: {aval.shape}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    848\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    849\u001b[0m       \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf\"had abstract value {aval}\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: Gradient only defined for scalar-output functions. Output had shape: (500,)."
     ]
    }
   ],
   "source": [
    "result = my_scalar_function_deriv(x, 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "threaded-reynolds",
   "metadata": {},
   "source": [
    "The above failure is expected. Let's call the vmapped version and inspect the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "characteristic-score",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'tuple'>\n",
      "1\n",
      "500\n"
     ]
    }
   ],
   "source": [
    "result = my_vector_function_deriv(x, 1.0)\n",
    "print(type(result))\n",
    "print(len(result))\n",
    "print(len(result[0]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sized-privacy",
   "metadata": {},
   "source": [
    "The `my_vector_function_deriv` function returns a tuple with an element for every argument for which we requested a derivative. Let's plot the result for $\\partial f/\\partial x$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "based-measure",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEOCAYAAACaQSCZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABOW0lEQVR4nO3dd3xUVdrA8d+ZVBIgk0ZJL5TQSUjogkgQROwUu64FdNddt6joNtF9dxXcfXdtq8DaXl0VwQYoIkEFkU7oHRISIKGl0kMy5/3jTDAmM8kkmTslOd/PJ58suXfufZy9yTOnPUdIKdE0TdO0hpjcHYCmaZrmHXTC0DRN0xyiE4amaZrmEJ0wNE3TNIfohKFpmqY5xNfdARglIiJCJiQkuDsMTdM0r7Jp06ZTUspIW8dabMJISEhg48aN7g5D0zTNqwgh8uwd86iEIYQwA1OBcCnl9FrHpgIDgPnWH00CZkopc1wapKZpWivlMQlDCJEJmIHkek6bjEoo2cCDOllomqa5jsckDCllFoAQIgOVOGydE+rKmDRN07Qf6VlSmqZpmkM8poXhCOs4RjEQBiClnOPeiDRN01oPb2phbASypJQLrIlijBBior2TCwoKEEJc/poxY4bLAtU0TWuJvKaFIaXMrvWjDcBTwAJb50dFRVFQUGB4XJqmaa2F1yQMIURm9cC4VQ6Q5op77ywoY9HWQk5fuMTAxDDG9+mMn483Nc40TdOazysShhAiCVgmhAiVUpbWOGTotFopJS8tP8C/lu/Dz2QiwM/Ef9fl89YPh5hz9wA6tAs08vaapmkexSs+JlvXW0yvlSymADONvO8/s/bzz6x93NQ/mg1/zGTb01fz4q392XvsNJNeX0PpuQojb69pmuZRPKaFIYRIAzKBidZ/P4Ea5K4eu1hg/RlAOLDMyFlSq/af4qXl+5k4IIZZt/TFZBIA3NA/mmhzG26bu5ZHP9zCW/dmXD6maZrWkomWukVrenq6bGotqXMVlVz19xUEB/iw+JdX0Mbfp8457645xJ8+38kLE/syKT22ueFqmqZ5BCHEJilluq1jXtEl5Wpv/XCIY+UXeP6WvjaTBcAdg+IZEB/K80v2UHbukosj1DRNcz2dMGopPVfB6ysOktmjAxkJYXbPM5kEz97Qi6KzFbyxSpe00jSt5dMJo5YjJeeJbBvAY2O7N3hur6gQxvXqxFurD1F+QbcyNE1r2XTCqKV3dAhZvx1JSqf2Dp3/yFVdOH2hknfX2C0hr2ma1iLohGFDY2Y99Y4OYWhyOO+vy6fK0jInEGiapoFOGE5x1+B4jpae59s9J9wdiqZpmmF0wnCCzJ4d6dg+gPfX57s7FE3TNMPohOEEfj4mbkyNZuW+kxSf1au/NU1rmXTCcJIb+kVTaZF8ub3Q3aFomqYZQicMJ+nRuR1dO7Rl4RZdUl3TtJZJJwwnEUJwQ/8o1h8q5mjpeXeHo2ma5nQ6YTjR9f2iAVi0VbcyNE1reXTCcKK48CB6R7dn2a7j7g5F0zTN6XTCcLLRKR3Jzi+h6MxFd4eiaZrmVDphOFlmj45ICd/uPenuUDRN05xKJwwn6x3dno7tA1i+W3dLaZrWsuiE4WRCCK5K6cjKfSe5WFnl7nA0TdOcRicMA4zp2YGzFVWsyyl2dyiapmlOoxOGAYYmR+Dva2LlPj2OoWlay6EThgEC/XzISAhl1YFT7g5F0zTNaXTCMMiwLhHsOXaaE6cvuDsUTdM0p/B1dwAt1RVdIpnFXlYfKOLG1Gh3h6NpWi1HS8/z0YbD7Dt+Gn9fE0OSwrkpLZoAXx93h+axdMIwSK+o9piD/Ph+/ymdMDTNg1gskjnf5/CPr/dSZZEkRgRz5mIln28p4JVvD/Dq7Wn0izW7O0yPpBOGQUwmwbDkCFYdOImUEiEc3/ZV0zRjWCySpz7ZzryNhxnXqxN/nNCDmNAgAL7ff5KnPtnOrXPW8p970hnWJcLN0XoePYZhoOFdIzhefpEDJ864OxRN04Dnluxm3sbDPDKqC6/dmXY5WQBc0TWST38+jLiwIB56bxMHT+rf29p0wjDQ0ORwANbl6vUYmuZuX24vZO73udw1OJ7fXd3NZqs/sl0A/7knHT8fE7/4bzYVlRY3ROq5dMIwUFxYEB3bB+iEoWluduL0BZ78eBv9Y838aULPeruIY8OCeP7mPuw5dprZKw66MErPpxOGgYQQDEwMZ31uEVJKd4ejaa3WMwt3caHSwj8m98Pft+E/e1f36sS1fTvz8rcH9IZoNeiEYbBBiWEcL79IXtE5d4eiaa3S2pwivtheyC9HdSE5sq3Dr3vqmhSQ8K9l+wyMzrvohOGIgi3w9R9h0aOwdR5UXXL4pYMSwwBYr7ulNM3lpJQ8v2QPndoH8uCIpEa9NiY0iLuHxPNx9hFy9AA4oBNG/aSEb/8Gc0bCujmw63P4dCq8cTWUFzp0iS4d2hIW7K/HMTTNDZbuPM6Ww6X8ZkxXAv0avyBv2shkfH1MzP0+14DovI9HJQwhhFkI8YQQYqad408IISYKIaYKIaYaHtC3f4UVM6Hf7fD4AXgiFya+BSf3wlvj4FzDSUAIwcCEMNblFhkerqZpP7JYJP/4ei9dOrTllrSYJl0jsl0At6TF8HH2EU6e1rtoekzCEEJkAplAMmC2cXwmkCOlXCClnAMkCyEmGhbQwW9g5QuQeifc8CoEtgchoPfNcPdnUF4AH98Ploan3Q1MDONIyXk9eKZpLvTNnhPsP3GGR0Z1wden6X/qHrwikYpKCx+uz3didN7JYxKGlDJLSrkAKLVzylTr8WrzgGmGBFNxFj77BUR0h/F/B1Ottyl2IIx7XiWVLf9t8HKDktQ4xgbdLaVpLjNnZQ7R5jZc27dzs66TFNmWocnhfLTpMBZL657t6DEJoz5CiDQbPy5FtUicb+1rcLoArn8J/NrYPif9PogbAllPN9g1ldKpPe0CfXW3lKa5yOb8EtYfKua+4Yn4NaN1UW1KRiyHi8+zJqd1/w57RcIAwoDaf5Xr/StdUFCAEOLy14wZMxy707li+OFF6D4e4gbbP08IGP8CnCuCtf+u95I+JsGA+FA25ZU4FoOmac3yn1W5tA/0ZUpGrFOuN7ZXJ0La+PHhhsNOuZ638paEYbZ3QAhh81hUVBRSystfDieMsiPQrjNc9aeGz+3UB3pcD+tmw/nSek9Niwtl/4kzlJ13fEqupmmNd+L0BZbuOMbk9FjaBjinvmqgnw83pUazdMcxSs9VOOWa3shbEkYpqpVRU+1/O0fnvvCLddCxp2Pnj3gMLpbDhv/Ue9qA+FCkhC2HS5sfo6Zpds3feIRKi+S2QXFOve7EATFUVFn4ascxp17Xm3hLwiimbivDDCClLHX63RpTirxzP0gcCRvfAkuV3dP6xZoxCcjW3VKaZpgqi+SD9fkMSQpv1KpuR/SKak98eBBfbHdsDVZL5BUJQ0qZTd3ZU2FAluujsSHjASg/AvuW2j2lbYAv3Tq2IztfJwxNM8rK/Sc5UnKeOwY7t3UBak3V+D6dWX2wiOKzrbNbyisShtVHtdZdjAFmuyuYn+g+Xo17bHqr3tMGxIeyJb+01U/N0zSjLNh4hLBgf67u2cmQ61/bpzNVFsnXO1tnt5THJAwhRJoQ4glgIpBpXdV9eTqtlHIakCSEyLSu8j5Ya12G+/j4Qt8pcGA5nD1l97S0uFBOX6xkv95QSdOcruzcJZbtPs71/aIcqkjbFK29W8pjEoaUMltKOUtKmWz9mmXtiqp5zizrAr851tXenqPPJJBVsOszu6cMiA8F0NNrNc0AX2wvpKLS0uQyII6o2S1Vdq71zXj0mITh9Tr2gsgesN1+oyc+PIiwYH89jqFpBvgk+whdOrSld3R7Q++T2aMjVRbJiv0nDb2PJ9IJw1mEgD63QP4aKLW9uEcIQVqcWc+U0jQnyys6y8a8Em5Oi653Nz1n6B9rJizYn2/3nDD0Pp5IJwxn6m0dk9/xsd1T0uJDyTl1lpJWOstC04zw6eajCAE39o82/F4+JsGV3SL5bu8JqlrZBBbnLIPUlLBE6Nwf9nwBw39t85S0ODWOsflwCVeldHRdbB6uotLCp5uP8NnmAnYWlHHhkoXIdgEMSQ7n9kFxl983TatNSsmnm48yNDmcKLOd2m9ONiqlA59sPsqWwyUMiDdmDbEn0i0MZ+t+DRzZAGds92/2izHjYxJ64LuG7PwSxr24kukfb+fE6Qtc3z+K+4Yn0i82hKU7jnHzv1fz8HubOFF+wd2hah5o+9Ey8orOcX2/KJfdc0S3SHxMgm9aWbeUbmE4W7dx8N1zsP9rSL2jzuE2/j707Nye7LxS18fmgT7fcpTH52+jQ/sA3rw3nVHdO/ykD/rsxUreXn2Il5bvZ2NeCa/fmdaqPtFpDftiWyG+JsHYXsasvbAlpI0f6fGhfLPnJI+PTXHZfd1NtzCcrXM/aBcF+5bYPSU1zszWI6Wtrv+zts+3HOXX87aQFm9m0SPDuSqlY50By+AAX34xqguLfjmcYH8f7vjPOlYftL/WRWtdpJQs3lbI8K4RmIP8XXrvEd0i2V1YzqkzrWcnPp0wnE0I6DYWDn4LlbYfpLS4UM5VVLHv+GkXB+c5Nh4q5rH5W8lICOOtewcSGlz/L3u3ju1Y8PBQ4sOCue/tDaw52Lr3JdCUrUfKOFp6nmv7NG+TpKYY1iUCoFU9izphGKH7NVBxBg59b/NwapwZgM35pa6LyYMUnbnIw//NJsrchrl3pdPG38eh10W0DeD9BwcRGxrEtHc3knvqrMGRap7ui20F+PkIw0qB1KdPdAjtAn1bVYtXJwwjJI4AnwDVyrAhLkwt4NvcShfw/XnhTsrOXWL2XQMICfJr1GvD2wbw5r0Z+PqYuP+dDXp/kVZMSskX2wq5omtko58jZ/AxCQYnhbPqgE4YWnP4tVG79dlJGEIIUmPNbG6Fe2N8taOQL7YV8mhmV1I6NW1FbmxYEK/dkUZ+0Tme/HgbUrbusaDWavPhUgrKLrilO6rasORwDhef53DxObfF4Eo6YRgleRSc2Amnj9s8nBpn5kAr24HvzMVK/vjZTnpFtWfqiKRmXWtQUjiPje3Okh3HmNfKt81srb7YVoi/j4nMnu5bz1Q9jvFDK2ll6IRhlKRR6nvOdzYPp1oXorWmHfjmrDjIqTMX+etNffDzaf6jN/WKJIZ1CeeZRbvIOakrALcmFovky+2FjOgWQUibJnZHleTBN/8DH94B838GG9+ES41b69OlQ1s6tAvgh1Yy8K0ThlE69YU2YZBju1uqb0wIQtBqxjGOl19g7ve5TOjbmf6xZqdc02QS/O/k/vj5CJ78eLveZ6QV2Xy4hMKyC1zbtwndURYLfP+/8HIafP8PKDoIh9fD4t/AKxlwZKPDlxJCMDQ5nDUHT7WKrlGdMIxiMkHSSNXCsPEgtQv0o1uHdq1mptS/svZRabHwhJMXOXVsH8gfr+3J+kPFvL8+36nX1jzX4m2F+PuayOzRyO4oiwUW/hKWPwMpE+DX2+EXa+G3O+Huz9Xv7dsT4OA3Dl9yWJcITp2pYG8rmCavE4aRkkbB6UI4udfm4bR4M1sOt/wd+I6Wnmf+xiPcPjCOuPAgp19/UnoMw7qE8/ySPRSWnXf69TXPYrFIvtpxjBFdI2kX2MjuqK//CFveg5HTYdLbEFJj74ykK+H+LAhLgo/uhZP7HLrkkORwANblFDcuFi+kE4aREq9Q3/NW2TycGhtK2flL5Ba17PUEc1fmADBtZLIh1xdC8NxNfamySJ7+fKch99A8x9YjpRSWXWB8n0auvdjxMax9FQZOhSufUotsa2sbCbd/CL7+MP9eqGy4qnRMaBBRIYGsP6QThtYcoYmqTEjeapuHW8MCvlNnLvLB+nxuSo02tJJoXHgQvxrdla93Hee7va2rIFxrs2THMfx8BKMb0x11+hgs+g3EZMDVf7WdLKqZ4+D6l9Usx1X/dOjyGYlhbMgtbvHjGDphGEkIiB8Kh36wOY6RHNmWdgG+LXrg++0fDlFRZeGhK41pXdR03/AEkiKCeWbRLi5WVhl+P831pJQs2VHIsC6NnB315eNQeQFufF21HhrS/RrofQt8/3e7G6LVlJEQxonTF8kratnrMXTCMFr8UDhzDIpz6hwymQT948wttoVx4VIV76/PJ7NHR5Ij2xp+vwBfH56+vhe5p87y5qpDht9Pc72dBeUcLj7PNb0b0R2V+z3sXggjn4CILo6/LvMZQKjq0w0YmKgqKLf0bimdMIyWMFx9t9ctFWtmz7Fyzl6sdGFQrrFkRyHFZyu4Z0iCy+45slskY3p25OVv9nOsTO+f0dIs2VGIj0kwxtHaUVJC1tPQPhqG/KJxNzPHwsAHYcv7cGp/vad2iWyLOciPDbk6YWjNEdENgiIg7webh1PjQrFI2HakzMWBGe//1uSRFBHMUOssElf584SeVFokf/tyt0vvqxlLSsmS7ccYnBRGWAPVjS/bvRCOboJRv1clexpr2K/BNwBWv1zvaSaTID0+jA26haE1S/U4hp2EUb2IbfPhljWOseNoGZvzS7lzcDwmUz0DjAaIDQvi4ZHJLNxawNqclrECt8oiyc4v4Y1VuTy7aBdPf76DV77Zz8p9J7lwqXWM1+w7foacU2e5preDi/UsFvjmrxCZAv1ua9pN20aq12790G6Zn2oDE0M5VHSOE6dbbstW77jnCvHD1Ced0sOqmVtDaLA/SRHBLW4c4901ebTx8+GWATENn2yAh69MZsGmI8xYuJPFvxyOrxNKkbhD+YVLfLAun7dXH6LQ2sUW7O+Dj0lQfkF1Y7YP9OW2gXE8NDK5wX1FvNmX2wsRAsd31tv3FZzaC7e8ASbHSujbNPSXsOktyH5HjYPYkZGgxjE25JY0bQW6F9AJwxXih6rveavBPKXO4f5xZlbuU6UFau84543Kzl/i861HuSk1uul1fpop0M+HP03owUPvZfPfdfncMzTBLXE0lZSShVsL+MviXZw6U8HQ5HCevCaFockRRLYLAFQxxw25xSzIPsLc73P4YH0+z93ct8X+sVqyo5CMhLDL//0NWv0ShMRBzxubd+PwZEgcCZvfhSseU6vBbegdHUIbPx82HCpusf8feOfHLm/TsRcEhtQ7jnHqzEWOlLSMVcpfbCvkwiULt2bEuTWOsb06MaxLOP/4ei/FZxtegOUpyi9c4uH3snn0wy1Em9vw+S+G8f6Dg7mhf/RP/li2DfBlVEoHXr09jSWPjiApsi2/eD+bpz7ZzqUqixv/C5zvwIkz7Dt+hvGOzo46vAHy18CQn4OPEz4Xp90NpfmQu8LuKX4+JlLjzGzMa7njGDphuILJB2IGqgJnNqReHscodV1MBvo4+whdOrSlb0yIW+MQQjDjul6cq6jihaW2y7N4moMnz3DjKz+wbPdxnromhU9+Pox+DhRr7N6pHfMfGsJDI5P5YH0+D727qUWNbXy1oxCAcY6OX6x5RX1IS73LOQGkTIBAs2pl1CMtLpTdhac5V9HyZj2CThiuEzsITu6G86V1DqV0akegn6lFLODLPXWWTXkl3JIW4xHda107tuOeoQl8uCGf7R4+E23r4VImvrZajVs8OJhpI5PxacSEAT8fE09ek8JfbuzNN3tPcP87G1rMAsYlO46RFmemU0hgwyefPg57FqtkEeCk9T9+gdDvVti9CM7Zb0GkxZupskiPf9aaqkkJQwiR4OQ4Wr7YDPX9aN3Syb4+JvrGmMluAQPfn2QfwSTgptRod4dy2aOZXQkPDuDphTs8ttDj6gOnuH3uWoIDfFnw0NDLC8Ga4q7B8fx9Yj9+OFDEY/O3eex/s6Pyi86xs6Dc8dlRm98FSyUM+JlzA+l/O1RVqAks9k6JVfvctITfZVsanTCEEK8Dm4QQN9f42QNOjaolih4AwmS/WyrOzK6CMq/uRrBYJJ9kH2V410jHPgm6SPtAP6aP6052fimfbj7q7nDq2HComJ+9vYGY0CA+fngoCRHBzb7mLQNiePKaFBZtLeCFr72jO86eJZe7oxwYv7BUwaZ31CB1Y1Z1O6JTX1XJdudndk8JC/YnMSK4RfQW2NKUFsYyIF1K+UmNn80XQjxvdMtDCDFVCDFbCJFp/ZothGjeXp+uEtAOOvSqZxwjlEtVkp0F5S4OzHnW5hRxtPQ8t6R5Tuui2i1pMfSPNfPckj2cvuA52+LuLiznvrc3EG1uw38fHETH9s5LtNNGJHHHoDhe++4gX+045rTrutqXO47RNyaE2DAHSuMfWA5l+ZDu5NYFqDVVPW+E3JVw1v76ntRY1VvQEgsRNiVhlAKpNX8gpSyTUj4JpDkjqAZMRiWtmcBsKWXdIk2eKjZDrTq11G1FpF2uXOu9n0w+23KUtgG+js+TdyGTSfDM9b0oOnuRl5bXX+bBVfKKznL3m+sJ9vfl/+4fSERbB6eLOkgIwZ+v60m/WDOPzd9K7invK6N/pOQcWw+XOta6ALXXRVAEdL/WmIB63QSyCvYssntKanzLmvVYU1MSxkbgISFEkRDiNSHEVTWOGf5pX0oZKqUUUsoBUspso+/nVLGD4GI5nNxT51CH9oFEm9t47UypS1UWvt51nMweHQj0a8YiKQP1izUzJT2WN384xI6j7h2UPFF+gbveWM+lKgvv3j+QmFDnbywFqiDjv+9Iw9dH8Ot5W6j0sum2i7aq7qgJfaIaPvl8CexdAn0mOVaRtik69WmwW6p61mO2F3/4s6cpCeN51Kf7yUAZMEcIUSWEKAK859O+O8RYB77rGcfY4qWDZWsOFlF67hLj+3j2gqWnrulBWLA/TyzY5ra1CmXnLnH3m+s5deYib92bQdeO7Qy9X7S5Df9zY2+2Hi7lte8OGnovZ1u4tYD+sWbHdmrc+akalO53q3EB1eyWOm87IaR0akeQv0+Lq94ATUsYWVLK5davJ6WUXYAM4EnA8E/81nGMidbvU42+n1OFJUFQeD0JI5Sjpec5Xu59tWi+3F5IsL8PI7pFujuUeoUE+fGXG3qzq7Ccud+7/vPN+Yoq7n9nAwdPnmH2XQNIjQt1yX0n9I3iun5RvLh8v9tbV446cOI0uwvLub6fA60LUPWeIlOgcz9jA+s+XnVLHVhu87Ca9Rji1d3L9jRpDKNWNxRSymwp5VyMH8PYiEpYC6SUc4AxQoiJtk4sKChACHH5a8aMGQaH5gAhVLfUEfstDPC+HfguVVlYuvMYo3t09NjuqJrG9e7E+D6d+FfWfg6cOO2y+16qsvDz/25iU34J/5qSyhVdXZtc/3JDL8KC/fndR1u9YiX4wq2FmARMcKTMRtFBOLxOtS6MXv8TnaY++O1baveU1LhQdhaUe/WsR1sanTCklMuBEiHEY9U/E0KMFkLsR7U0DGNNTDU/Fm4AnrJ1blRUFFLKy18ekTBAdUsVHbA5y6JXVHv8fbxvAd/anCJKnNkdJaXNHQqdacb1vWgb4Msj7292yS+1xSJ5fP5Wvt17kr/e2McttYbMQf789aY+7D1+mjdW5br8/o0hpWTR1gIGJ4XTwZGZY9vmAQL6TDY8Nkw+0PVqOLDM5gQWUCu+Ky3Sa1pzjmrSwj0p5WYp5d9r/Hs5MAuY46zAbBFCZNb6UQ6umZnlPLGD1PcjG+ocCvD1oWdUe69rYVR3R13ZvRmfmCsvwqa34a1r4blYeDYc/rcnfDIN8tc5LdZqHdoF8o9J/dhz7LTh+2ZIKfnzwh18tqWAx8d25/ZB7quxNaZnR8b07Mi/svZxuNhztxPdcbSc3FNnHeuOklIljMQREOKiKd1dr1ZjGDZ+j+HH3oKWNvDtUMIQQtwshFgqhNhv/ZonhLhfCNG++hwp5VwppWEfW6zrLZYJIcy1DnnXQHtUKph86+2W2na01Cu6DAAqqyws3Xmcq5rTHXV4Pbw2FBY9CudOQf/bYPivVXLdtwTevBo+vAPKC50a+6iUDjwwPJH/W5PH0p3GrFOQUvLckj28tzafh0Ym83MX7G3ekBnX90IgmLFwp8euFVi49Sh+PsKx6bQF2VBySM2OcpXkq0D42O2WimgbQFxYENl5pa6LyQUaTBjWld1XAwv4sRURCsxFdU39u2biMIq1K2q6lLK0xo+noGZseQ//IOjY2+7Ad1pcKBcuWdh7zHV9682xLreY4rMVXNuniWsvtn4Ib42HqktwxwL4+VoY/wKM/jNMegt+uxtGP60GGGdfAXlrnBr/E+NS6BsTwu8+2mrIe/7i8v3MWZnD3UPimT6uu0fU14o2t+E3Y7qyfM8Jvt5V/6ZA7mCxSBZvK2RE10jMQQ5Mj935KZj8oMcE44Or1sasti3Y/7XdU9LizGTnl3hsUm6KehOGECIR2CSlfMjagpgLZAGzgWRgLBAB5AohDJ6aAMACIcQT1q+ZwDLr4Ld3iR0IR7Nt9n+metkCvi+2FxLk78OV3Ts0/sXbPoJPH4L4ITD1O+g6pu6ApX8wXPFbmLYCAtrDuzdCzndOiNx6eV8Ts+8aQJC/D/e9vYGTpy865bpSSv7x9V7+lbWfiQNimHFdL49IFtV+NiyRlE7teGbhTo+rrLo2p4jCsgtc39/B7qidn0HyKGjjmhlnlyWPguM74MxJm4dT40I5cfoiBS1ob/mGWhipwLxaP0uUUn4spcyVUmZJKSejBrvfMLo0iJQyR0o5y/o13SuTBaiB70tn4UTdvvNocxsi2wV4xThGZZWFpTuOMSqlCYv18tbAZw9DwnC4/SMIaqDYXmR3uP9rCEuG96eoefBO0jmkDW/ck0HR2Yvc9ca6Zu+dUWWR/OGzHbz8zQGmpMfy/M19XL5NbUP8fEz85cbeFJRd4JVvDrg7nJ+Yv+kI7QIdrBhwdBOUHVYrsF0t8Ur13c4eGWnWKdPZed7x4c8RDSUMh4b4rd1FmcD0ZkfUGsSkq+82xjGEEKTGmr1ixff63GKKzlZwbWNnR505CR/dDeZ4mPIe+LVx7HXBEXDPQghNhA/vhFPOK/HRJyaE/9ydQe6ps9w+d22TWxrlFy7x0HubeH9dPg9fmczzt/Tx2O1hMxLCuDk1mrnf55Bz8oy7wwHU+/fl9kKu7xfl2IeQ6u6o7uOND662qP4QEGI3YaR0bkeAr4ktXvC77KiGnuSNQHqtloPN8h/WsQXvGoB2l9BENY/7SN1S56CasrmnzlLi4bvEfbG9kDZ+PoxqbHfUl7+DC6UqWbQxN+61wRFwx0fg4wcf3Gp3tW1TDO8awRv3ZHCo6CzXv7KKrY38Rd9+pIzrX17Ft3tOMOO6nkwfl+JR3VC2PDk+hQBfH55ZtMsj+toXby3kYqWFyemxDZ9ssajuqC6jG/8cOYPJBxKvgBzbCcPPx0Sf6Ja1gK/ehCGlLAMEkFOjnPnmmmswar/EmcG1WEKobqkGpuRtPuy5D1qVRbJ05zGuSulAG/9GdEft/Ax2fQ5XPgkdezbt5uY4lWxK8uDzR5y6ZmN41wgWPDQUkxBMen0N//v1Xs5X1L9Oo/hsBc8u2sUNr67i/KUqPpg6mHuHJTotJiN1aBfIrzO7smLfSZZ5wAD4RxsP062jg7s1Ht0I5Ufc0x1VLXEklOapWVo2pMaZ2VFQTkWld8x6bEiDbWXrGovQ6nLm1n93EUJssE6tTQAQQvQHnFyAvgWLyYBT+2x+Qu4bE4JJePaK73W5RZw6U9G4xXoXyuHLx6Fzfxj6aPMCiB+iZlLtWazWbzhR7+gQFv1yOON6d+Klbw5w5d+/5V9Z+9hx9Mf9SsrOX2LV/lP84dPtDH1+OW+tzuX2QXEs++1IMhKavvmRO9wzNIFuHdvy7OJdbl2ZvP/4abYcLmVyeqxjLbOdn4KPP3S/xvjg7Ekaqb7baWWkxoVSUWlhd6H3bltQk0O7o1tbGjX//ZC1jtNcQFr/z80BBjg9wpaquhDh0U3Q5afrEYP8fUnp5NkL+JZsP0agn4lRKY1YrPfDi3D2BNz+Ifg49OjVb8gjcPAb+OopNXge0bX517QKC/bnpdtSuXNwPC9/s59/ZakvAF+ToNK6i52/r4kb+0fxwBVJdDO4iKBR/HxMPHN9b26bu5Z/f3eQ347p5pY4Ptp4GF+T4EZHdmu83B2VqfbudpeIbtC2k5q5N+CeOof7WyvXbs4vcWhvdk/X5N9a6wylOUKIVOu/NzstqtYgOg0QahyjS+0F7Gpv4M82F1BlkY3a19kVqiySJTtUd1SQv4OPUHkBrHkVek9Uuw86g8kEN70Orw6Chb+Ee79UP3OigYlhvHv/IE6UX2BNThF5Ree4WFlFu0A/undqx6DEMMffAw82JDmc6/pF8fqKg0xMi3GsOqwTna+o4qONRxjbu5Nj+4IcWQ+nC6DXM8YHVx8hVCvjwHKVxGo9f51DAunYPqDFDHw3+7fLWiZEJ4vGCmgHHXraH8eIDeXMxUoOesjslZo2HCrm1JmLjeuO+vZvqsLn6D85N5h2nWDs3yB/DWx607nXrqFD+0Bu6B/Nr0Z35fGxKTw0MplR3RuRML3AH8b3wM8keHbxTpff+/MtRyk7f4l7hiQ49oKdn4JPAHQbZ2hcDkkcqSoUnNhV55AQgv5eMuvREXYThhAiUQjxeu3KtJoTxaSrFoal7oDY5Vo0HjiH+8vthQT4mhyfHVWSB1veh/T7IDTB+QH1vx2SRsGyp6HsiPOv30p0CgnkV6O7krX7BN/scd0AuJSSd9bkkdKpHRkJDiy+s1hg10Jrd5ThRSYaljhCfc/7webh1LhQ8orONXt9jyewmzCsdaGmAwOEEF8LIZ4zemFeqxOToaaXFtVdOJUYEUx4sD/rDxW7Pq56VHdHjeregeAABz9dr34ZhAmG/sqYoISA6/4F0qIG1bUm+9mwRJIjg3lmkesGwDccKmF3YTn3DE1wbLD76CbVHdXzeuODc4Q5FkJiIW+1zcPVO/Bt8eBZj45qcFqtlPIFKeXVqBpSD1mLED7mivpRLV71wLeNbikhBIOSwliXU+wR8+OrbTxUzMnTFxnvaHnuMydg87tqnwIjK4mGJqipunu/hH326/to9fP3VQPgeUXnmLvSNcuqZq84SGiQHzc4UgoEYNdnarGeJ3RHVYsborpFbfyu9okJwcckPHoSi6McHsOwlgJ5Uko5FtgMzLJWrb25oddqdkR0UytF7YxjDE4K52jpeY/aTH7JjmME+JoYneJgd9Ta11Tp8mG/NjQuAAY9DOFd4avp6p5akwzvGsH4Pp149bsDHCkxtgT67sJylu85wc+GJTo2HiQl7F4ISVe6Z7GePfFD4MxxKK6bZIP8fenesV2LGPhu6n4Yy60FCacAQgjxkTV56PGOxjCZIGaA3RXfgxLDAVWMzRNYLJIlOwq5snukY91Rly6oNRIp10KEC5bo+PrDNTPVL+2aV42/Xwv2h2t7IhD8z2Jj9wp57buDBPv7OD7YXbgFSvOh5w1GhtV4cUPV93zb1ZT7x5nZkl+KxeI5vQVN4YxZUh9bCxBOBZKt4x2vWRfyaQ2JyYATO+Fi3dlQXTu0JSzYn7U5njGOsSm/hOPljZgdteszOF8MAx80NK6f6DIaUibAyheg7Kjr7tvCRJvb8MhVXfhq5zFW7rNdjbW5ck+dZfG2Au4cHE9IkJ9jL9q1UO1DkXKtITE1WWR3aBNmt/x+aqyZ0xcryTnlebMeG8Npk9at4x1zreMds1D7bS/VA+UNiMlQg7UFdWcmm0yCQYlhHtPC+GJbIf6+Jkb36OjYC9bPVV1EiSONDay2sX9T7+nXf3TtfVuYB65IJCE8iBkLdxoyAD7rqz0E+vnwwBU2y9PVJaUqK5N4RcPVjV1NCOs4hp2B7+rKtV4+jmFIGU3reMcLUsqxUspDRtyjxahexNbAOIa7t9O0WCRf7TjGyG6RtHWkO6pgs6r1k/FA3T0ujBYaD8N/Azs/gdzvXXtvo1iqIH+tmnH25ROw+LewYhbsz4IKY56NAF8fnr2hNzmnzvLC0r1OvfamvBKW7DjGtBHJRLZzYKEeqHUOxQc9rzuqWvwQ1R16uu6U5KSIYNoF+nr9wHfLWXXkrYLCILyL/XGMJPVJal1uMbFhrl19W1N2fgnHyi/wZJ8Ux16w4T/gF6RmR7nDsEdhy39hyRMw7XvnlCJxh/OlsOktWDcbTlu3qA1or7b5PW/tqgwIgbS7YPhvITjcqbcf0S2Su4fE88aqXEZ2i2REt2bs224lpeS5L3cT2S6AB65oRJHGXZ8DQnU5eqLL4xir6xRENJnUAj5vH/huaMe90dbB7AQXxdM6xWSoUgc2puR169CO0CA/t3dLLd6mFutl9nSgO+p8KWxfAH0nu28mi18b1TV1YhdsfMM9MTSHlLB1Hrw8ALJmqD7yiW/C4wfhqcMwPRd+XwB3fqx2Klz7b3gpVb3vTvb78T3o0qEtj83f6pTFZx9tPMzGvBJ+N6ab42t5QI1fxA+Dtk3Y3dEVOvdVH5LqGcfYe6ycsxc9a4fDxmioSyoJtSeGufoHejDbADHpcPakKpNcixrHCGfNwSK3rceoski+3F7IqO4dHOuO2vkJVF6AtLuND64+KRPUCvBv/wpnT7k3lsY4Xwof3gGfTlXrS6Z+B3d/Dr1vUfuBVPMPVqudJ74BD69RSeXj+1XJ96pLTgsn0M+HF2/tT+m5S/x63hYqq5peqruw7Dz/s3g3g5PCHNvzotrJvXByt+cs1rPFx0/9LtuZKZUaF4pFwrYjDu1L55EaShjJwEF+ujHSFOPCaaUuL+Cz3S01rGsER0vPk3vqrAuD+tHGQ8WcOH2Rax1drLflA4hMgag0YwNriBBqmm3FWVju5iJ1jjq5F+aOgv1L4eq/qm1po1Ibfl2HFPjZl6pbavO78OHtTh3b6BUVwrM39GLlvpM8t2RPk65hsUie+mQ7lRbJzFv6Nm7b2l0L1fce1zXp3i4TMxCO71TPXC39L6/4LnVtTE7U0ErvJ4GBQKkQokgIMQ9IE0KM0iu9nahDL9WUtTPwPbKr6jdeYdD0xoYs3lZIoJ+JqxxZrHfqgOpe63eb6we7bYnsDoMegux34Wi2u6Op35FN8OY4NcX63i9g6CNqVzdH+fhB5tMw4V+wfxm8P9mpCxhvHRjHz4Yl8MaqXN5Yldvo17/8zQG+23uS349PIT48uHEv3v25+mPc3sHV4O4SO1AV2bQx6zE02J+E8CCv3oHPkQ2UnpRSmlB7dm9EtTrm8mMSWWqtM3WzHutoIh9f9WncTsKICw8iMSLYsPnw9amyLtYbndLRsf7mrR+oulF9PaghOnI6BEeqOlM2Cj16hIPfwjvXqWJ69y+FuMFNv1b6z+Cm2XDoe/h0mlP/m/8wvgfjenXiL4t38WYjksZHGw/zz6x93JwWzZ2D4xt30+IcOLbdc2dH1VRPuR9Q3VKbD5d6VLmfxmhMaZDNUsoXgAVSyi41ksgCIBT4PWor1yLrwr0EQyJuqWLSoXCbWh1tw8hukazJKXL5jmjrctTOeg51R1kssPVDSL4K2jei9LnRAtvDmGfUNN+tH7g7mrryVqv9ycMS4b6lEObguoT69JsCY/6iyoBnPd3861n5+ph4+fZUxvbqyLOLd/G3L3fXO6YhpeSd1YeY/vE2rugawd9u6tP4fc6ru6M8efyiWlAYhCXDYdsJo3+smZOnL1JQZvv33NM1ZR3Gc9X/w5pE5lrLhKRbk0g6kAU8KYR4wFmBtngxGWC5BMe22Tw8olsEFy5Z2HjItc3ZxdsLCfL3cayU+aGVao/lfrcZH1hj9b1VvcdZT8MFDxp0LNwG709R1U7vXqj293CWob+E9Pth9Utqdzon8fMx8crtadwzJJ45K3O44dUfWJtTd1LGkZJzPPL+Zp5euJPRKR2Yc1c6gX6N6GKrtutzNY5jjnPSf4HBYgeqFoaNVkT1tgXe2i3V6MnptbdrtXE8VwhRDJSgZllpjohJV9+PbFAPXC2Dk8Lx9zGxYt8JhneNqHPcCJVVFr7acYzRPTrSxt+BX/Rt88G/neeVbQBVt2v8CzBnFHw3E8b9zd0RQdFBeO8Wta7i7s+cvoYCIWDcc1C4Vc2c6tjbaTW9/HxMPHNDbwYnhfP0wp3cOmctSZHBpMaG0sbfxP7jZ9hwqBg/HxOPj+3OwyOTGzfIXa0kDwqyYbTzWkmGi0lXLdnSvDr7v6R0ak+Ar4kt+aVM6Ovh4zE2GLLSmx9bGZqj2nWCkDg4vN7m4SB/XzISQ1m5z3XTQ78/cIrisxVMcKQ7quoS7FkMKePVGghPFJWq9l1e97rNQUmXKi+Ed28ESyXc9SmExBhzH98AmPyOGhD/5EGocu4agGv6dGbF46P42019iAsLYtWBkyzeVsj5S1X8/MoufPPYlfxiVJemJQuAHR+r7729qCh2jPUDn41uKX9fE72jQ7x2Bz5Dlr9KKZdjXDJquWLS7SYMgBFdI3luyR4Ky87TOcT4P8qfZB/FHOTnWHdU7gq1GVTPG40Oq3kyn4G9X8Hnv4Sp36o/pK52vgTeuxnOFcM9CyGym7H3C4mBCf+E+ffAqv+FkU849fJt/H24fVActw8yoMtoxyeqK9GInRqN0qEn+AWr3oK+k+ocTo018+7aPCoqLfj7etefSYeitc6AWiqE2G/9mieEuF9PrXWy2IFqDKC8wObhkd3V9Nrv9ho/W6r8wiW+3nmM6/tFOfZQ7/xMdUcle3iF+zZmuPYfcHw7/PCi6+9fcVaNWRQdgFv/+2MtMaP1uhH6TIIVM93funLUyb3q/6feE90dSeP4+EJ0mppebkP/ODMXKy3sOVbu4sCar8G/BEKI14GrUbOhZqF23gtFTa0tEUL8WycOJ2lgAV/3ju2ICwti6c5jhoeyZHshFyst3JzmQFdJdXdU93HgF2h4bM3WY4Kq9bNiJpwwdr+Hn6isgI/uVp88b3lDbQLkSuNfUNOLP31YxeLpti9QU7Rr1WXyCjEZairwpbqbn1VXrvXGQoQN1ZJKBDZZZ0HNlVLORY1NzEatxxgLRAC5Qoh+hkfb0nXqAz7+dudwCyEY17sTPxw4RfkF55V+sOXj7KMkRQbTLyak4ZNzV6puFk/vjqrpmllqsHnB/TZ/qZ3OYoHPfw4HstTCOndMEW0Tqu59cjesecX1928MKWHHAki4Ato5WE7fk8RkqPGpgi11DkWFBNKhXYBXrvhuqIWRCsyr9bNE66ZJuVLKLOvmSRnAG3rtRTP5BkDnfnZbGABje3XkUpXk2z0nDAvjcPE51ucWc0tajGNz5nd9Bv5t1eZF3qJtB7W47cRO4/fNkBK++C1sn69m+wy4x9j71af7OFVja8UsNQPJUxVsVgv2+nhZd1S1y70FdbulhFCVa71xam1DCcOhCetSyhzUIr7pzY6otYvJUL8sdorHpcaGEtkuwNBuqU+y1U51N6ZGN3xy1SXYvRi6jfPc2VH2dM2EIY+oUuy7FxlzDylh6R9UifLhv4UrfmvMfRrjmpmqq+fLx22uFfAIOz4Gk5/n146yp20khCbWu+L7UNE5p1T/daWGEsZGIL1Wy8Hm2gopZSk/LVJoCCHEE0KIiUKIqUKIqUbfz+Vi0qHyvCpgZoPJJLi6Z0e+23vSkFXfVRbJvA35DO8SQbTZgQRwaJXal6HXjU6PxSVGP63Ksnz6sN33vFm+ew7WvgoDp8HoPzv/+k0REgOjnlIFDvcsdnc0dVksanZUl0zVjeatYjLU1FobSbm6EOFWL+uWaqj4YBkgUCU/qidCbxZCPGbvJc4MrjYhxEwgR0q5QEo5B7WHuJe2We1ooBYNwPg+nTlXUcXy3c7vlvpu7wkKyi5wh6NTJHd9pqYQdsl0eiwu4euvZisFtFWzl2zsltYkUqp9LFbMhNQ7YdzznlGMsdqgh9RCviXTbe4n71aHvofTBd7bHVUtdiCcOQZlR+oc6hsTgkl434pvR4oPLgdCpZSf1Ph3FyHEBuvU2gS4vE+Gc5aR2jdVSllzh5h5wDSD7+laIbHQtmO9CWNwUjid2gfySXbdB7G5/rsun8h2AY5tlFRVqbpyuo31vu6omtpHwW0fwrkiePem5u+dYamCRY/Cqn/CgHvhupfUSnNP4uOnpheXH4Xv/+7uaH5q83sQGOKZFQMa43L1hrrjGMEBvnTv1N7r9vh26CmuXQ5ESvkQalrtXOCgEKIKmA84d0VQDUIIW5srlKLGTloOIaw78NlPGD4mwU1p0Xy37yQnTzuvfPWRknN8u/cEt2bE4ufjwKORt0r9kfXW7qiaovrDbR+ogdZ3rmt6S+N8qdqLIvsduOIxNSupMSXKXSlusKr7tfoVVZbeE1wog90L1doLb/4QAqoF5xtodxJLenwom/NLmrUhlas1+WOPlHJOjWKD6VLKrlJKI1eihAHFtX5W+9+XFRQUIIS4/DVjxgwDQ3OymAz1h+us/W1Zb06NpsoiWbjV9iK/pnh3TR4CmJLh4E5oOz9T+3h0GeO0GNwq6Uq4fR6UHII5I+1WHLXryCaYc6WaOjv+7zD6T57VDWXLmGfVH+YlT3jGAPiOj9Vujal3ujuS5vPxU+Vo7CSMjMQwzlZUsavQexbwNbudbK1Y64qlo2Z7B4QQdY5FRUUhpbz85XUJA1Q5bju6dmxH35gQPt7knG6p0xcu8f66fMb36UxMaFDDL6jZHeXvwPneImmk2uXOxw/eugaWP2tz97SfOHsKvnwC/jNabVh075cw8EHXxNtcbTvAqN/DweWeMQC++T1VWsORXQa9QUy6Kv5oYyOrgQlhAKzPtfu51+N4WMdqvUpRrYyaav+7ZYjqD8Kn3m4pgEnpsewqLGdTXvMHzj5Yn8/pi5VMG5Hs2AvyfoBzp7xrsZ6jOvWBqStUwbvv/wEvpcI3f4Wjm9S2p1KqOlAHlsPCX8E/e8H6OSpJ/GItxA1y939B42Q8qP5If/V7p27r2mgndqv3OPVOz2+ZOSo6HaouwrEddQ51CgkkNqwNGw55T8IwpPigQYqp28oww+UpvS2HfzB07NVgwrg5NZoXvtrDm6tyGRDf9OmHFZUW3lx1iCFJ4fRxZGU3qD0K/IKg69VNvq9HCwqDm+eo/SRWvmD9mmU9KLg8IdA3UO0uOOQR44sIGsXHV3WhvT1eDdRf9Qf3xJH9Lph8oc9k99zfCDV7C2Lq1g3LSAhjxd6TSCkbv7GUG3hNwpBSZgshSmv9OIyWWkY9JgO2faRm3NgZNA0O8OX2QfHMWXmQw8XniA1rWtfQvI2HOVZ+gedv6ePYCyxVqjuq69UtqzvKlrhBcOcCOHNCrTkpzlF97IEh0KEHxA1tGe9BwjBVnPCHF6H/bc7Z9a8xKs7Clvegx/Vq0VtLERIN7TqrD3+D6k7oHJgQxifZRzl48ixdOrR1Q4CN401dUgAf1Vp3MQZV16rlicmAitOqYmc97hkaj0kI3ll9qEm3OXuxkhez9jMwMYyR3Rz8Rc1bDWdPtIzZUY5q20F1UY14DK76o9rNrktmy0gW1cb8RY3dfPWU6++9fb6aITWw5a3FJSa93oFvwGu6pbwqYUgppwFJQohM6yrvg7XWZbQc9dSiqalzSBsm9O3M++vzOXG68fsEv7Eql1NnLvLkNSmON4l3fQa+bVpud1Rr1b4zjJwO+75Se4a4ipSwbg507KOm+rY0MRlQkmtzfU9SRDARbf3Z4CUD316VMACklLOsRQ/nWFd7t0zhyRAUAXlrGjz10cxuVFRaeHl54+bSHyu7wJyVOYzt1ZG0OAfHQCxVsGshdB2jxlq0lmXwwxDRHb6aDpca/wGkSfJ+UEUgBz7Ycga7a4quXsBXt5UhhCA9Poz1uoWhNYsQkDBc/TI1MD8+MSKYWwfG8sH6fHY7OKdbSsmfPt/BpSoLT13Tw/G48te0vu6o1sTHD8bPUmtRVr/kmnv+8CIEhasxlJaoetajnWnyGYlhHCk5T0GpC8rsN5NOGJ4sYTiUHVabyTfgt2O6E9LGjycWbHNo5ejH2UdZtus4vx3TjYSIRrQUdn2uZgZ1Hev4azTvknSlmi79/T+ML4FeuA32fw2DHm5Z40E1+QdDx552Zz0OSQoHYM1B+wt1PYVOGJ4sYbj6fmhVg6eGBfvzzA292H60jJlf7an33B1Hy/jDp9sZkhTO/cMTHY/HYvmxOyrA82d0aM0w9q+qBPrS3xt7n1X/VFv7DnzA2Pu4W0wGHM1Wv0O1pHRqR3iwPz8caGYNMxfQCcOTRaaoproDCQNgQt8o7hkSz9zvc3lzVa7Nc3YXlnPXG+sID/bn5dtT8XWkZlS1w2tV9c2WuFhP+6mQGBjxuFr9fcCgmeunDqgJFBn3eXcZc0fEZMDFcji1r84hk0kwtEsEqw6cQnpCeZZ66IThyarHMQ794PBL/jShJ1f37Mizi3fx5MfbLhcnPF9Rxds/5HLzv1cT4OvD+w8OJqJtQOPi2fmZ6o7qprujWoUhv4CwZFX2xIgB8OUz1OLPIY84/9qe5vLAt+1uqWHJ4Zw4fZGDJz2s1HwtOmF4uvjhUJbvcF+yr4+J1+4cwNQRSczfdISBf8ti2PPfkPqXr5mxaBdp8WY++8Wwxo1bgGpK716o1h4EtGvCf4jmdXwDVAn04oNqbw9nyl+nFn8Oe1StcWnpwruoxZ52Br6HdYkAYNV+z+6W8pqV3q1WzXGM0HiHXuJjEvx+fA+mZMSyeGshOafOEB4cwNW9OjIoMaxpJQgOr4PThbo7qrVJHqV2C1z3mhq7csa+7VLCsj9B206qFdMamEwQPcDuAr7YsCDiw4NYdaCIe4c1YlzRxXTC8HQ1xzFS72jUS5Mj2/JoZlfnxLHzE9Ud1X2cc66neY8xz0DuSvjsYXh4DQSHN+962e+oDyDXv9K61vLEZKiaZBfP2Jw0MjQ5gsVbC6issjRubNGFPDMq7UcmE8QPU9tWumtAzFKlptN2vVp3R7VGfm3glrlwvgQ+vl+Vtm+qsiOw9I+QOKJl7HnRGDEZIC1QYHs3iOFdIjh9sZJtR8tsHvcEOmF4g6Qr1XqMIjftipa/Bs4ch143uef+mvt16qPGM3K+ha//2LRrWCyqHLysUtvWtsRV3fWJtlartbceIzkcITx7HEMnDG9Q3W98YLl77r/jEzWbRc+Oat3S7obBP1fjGatfafzrV8xUGzVd/T8Q5rn99IYJClOzzo5usnk4LNifPtEhfLf3hIsDc5xOGN4gNEHNsjBqPnx9qirV7KhuY1tXf7Nm25i/QM8b4Os/wJpXHX9d9ruw4nnofwek32dcfJ4uJl21MOx0L49O6cjmw6WcOlN3hz5PoBOGt0gerQa+XVUQrlreKjh7UndHaYqPL9zyhkoaS3+vdumrb0xDSlj7Oiy0loOf8M/W1xVVU0yG6t4tO2zz8OgeHZASvtnjma0MnTC8RZdMqDwP+atde9+dn4JfMHQZ49r7ap7Lx08ljYHTYO2rMOdKyFlR91NzySH46C5V+bb7eJjynlrb0ZpdHsewPb22V1R7OocEsnz3cRcG5Tg9rdZbJAwDH381jpF8lWvuWVWpakd1v6blFobTmqa6qm3iFWol+P9dr3bpixkIfoFwcp+aLOEbCJkzYOijasZfa9ext3pPjmxUG3LVIoTgqpQOfLr5KBcuVRHoZ3u3TXfR/w96C/9giBvi2oHvg9/A+WKbD7amAdDjOvhVNlz/MkR0U9O/dy+GqgoY+QT8chMM/41OFtV8/aFzf7srvgEye3TkXEUVa3I8r3qtbmF4ky6ZaoVs2VG1V7DRtn4AbcJ0d5RWP782agZV2t3ujsQ7xKTD+rlQWaESSC1DksNpG+DLl9sKGdXds8qm6LTvTbpkqu/7vzb+XudLYc8X0GeizYda07QmikmHqotwfLvNw4F+Pozt1YmvdhzjwqUqFwdXP50wvEmHHhCaqEpOG23X5+qh7ner8ffStNbkcuVa2+sxAG7oH8Xpi5V8t/eki4JyjE4Y3kQI6DFBzUi5YHD5gK0fqj7pqDRj76NprU1IjCq8aGfFN8DQ5HAi2vqzcOtRFwbWMJ0wvE2P68FyCfYZ2C1VckhN3+07pXXPmdc0Iwjx4wI+O3x9TFzbpzNZu09w+sIlFwZXP50wvE10uvp0snuhcffY+qH63neKcffQtNYsJh1KcuGs/ZlQN6RGU1FpYeHWAhcGVj+dMLyNyQQp16oyIRXnnH/9qkrI/j9IGgXmWOdfX9M0teIb6p1emxprpkfn9ry3Nt9jtm7VCcMb9bwBLp2DfV85/9oHlkH50dZd70fTjNa5PwiT3RXfoBbx3TU4nt2F5WTnl7gutnrohOGNEoZD+2i1TsLZNr6pury6X+P8a2uapgS0hQ691EZS9bihfxTtA32ZszLH4UtfqrI0Nzq7dMLwRiYf6DtZrfo+7cSaMyV5sH+ZWoDl4+e862qaVlfcYNXCqKd4Y3CAL/cOS2TpzuPsOVbu0GUfeGcjMxbudFaUP6EThrfqd5vaiGbHAuddc8Nc1UzWK3Y1zXjxQ+DSWTi2td7T7huWQLC/D/9ctq/BS64+cIoV+04SE9rGWVH+hE4Y3iqyu1ojseUD52zdeqEMNr4NvW7Ug92a5gpxQ9X3vDX1nmYO8ufhK5NZuvM4K/fZX8hnsUhmLt1LVEggdw6Od2akl+mE4c1S71DlBQ6vb/61Nr4FFadh6K+afy3No82aNYvp06eTnZ1NVlYW06ZNc3dIrVP7zmpztPz6EwbAgyOSSAgP4k+f76DczrqMN3/IZevhUh4b292wKrc6YXizfrdBoBnWNGG7zJoqK2Dd65A4AqL6OyMyzcPNmTOH0aNHM3v2bGbOnOnucFqvuKEqYTTQSxDg68MLk/pxtOQ8v/lwS52B7Y2Hipn11V7G9OzITanGFSb1ioQhhJgqhJgthMi0fs0WQiS5Oy638w+GAfeq2lIlh5p+nex34HQhDHvUWZFpHsxsNlNSUkJJSQnz58/HbDa7O6TWK34InCuCU/sbPDUjIYynr+vJ8j0nuOfN9RwpOYfFIlm0tYC73lhPTGgbnr+5D8LA6gzeVN58MjAVyAYelFI6Ps+sJRs4VbUw1s2Gcc81/vUXz8CKmRA/XG0D20I8s2gnuwocm1XibD2j2vP0db0cPn/BggVMnz6dpKSky3/AJ02aRE5ODnPnziUtzZh6XtnZ2ZjNZpKS9Gcvt6kex8hfDZHdGjz9riEJBPr58IdPdzB85rf4+5ioqLLQJzqEN+/NILytsTsaek3CkFKGujsGjxQSDb1vgU1vqxZCu06Ne/2aV9We3bd+oOtGucnEiRMBmDdv3uVP+9OmTSM9Pd3up//p06c3eN0pU6bYTTYLFiwgMzOTrKws3S3lTuHJEBypBr4H3OvQSyalxzIkOZzF2wopOVtB3xgzY3t1xNfH+A4j4SlLzusjhJgqpZzTmNekp6fLjRvtr6JsUYoOwqsDIe0emPC/jr+uvABeyYDkUWq/Zc2tQkNDyc3NxWw2s2DBgsuJxGjJycnMnj2bzMxMl9xPq2XeXVC4BX5te38MVxNCbJJSpts65hVjGHB5HGOi9ftUd8fjUcKT1aeTTW/DMQcfOinhi9+BpQrGPGtkdJqDJk+ezJw5cygtLTW0myg7O/sn/05LS2PZsmWG3U9rQPxQKM2HsiPujqRB3tIltREorR63EELMF0IUSyntrlorKCj4yeDP008/zYwZMwwP1G1G/UFtevT5I/DAcvBp4P/aLf+FvV/C1f8DYboP2xNMnz6dMWPGkJSU1GDroqldUtnZ2YwePZqSkh9rE5WWlpKcnNy0oLXmS7hCfc9dCf1vd28sDfCKhCGlzK71ow3AU4DdhBEVFUVBgeeUBTZcUBiM/zvMvweW/RnG/c3+uUezVesicQQM/rnrYtTqlZSUhNlspri4uMFzmzrmkJaWVue1OTk5TJ48uUnX05ygQ081jpHznU4YtgghJgINbbZQLKWcZj0/U0qZVeNYDqC3gqut142QNw3WvqoGw4f8ou45hdvgvZvVA3rLm6ouleYxpk2bZvgf7/T0dGbNmoXZbObgwYN6aq27mUzqw1vOCtVV7MGTT9ySMKxdSQ4VQbKut1gmhAiVUpbWOKSn1doy9m9wugCW/h5O7IKr/qRmTlWcVftcZD2jWiP3LIK2ke6OVqslLCzM8D/eaWlphk3V1Zoo6UrY8TGc3AsdUtwdjV0e3yUlpcwRQkyvlSymAHoeoC0+vjDpHVj+rFqfsfk9VQr9XBFUXlAbI930euOn32qGmTZtGpMmTSIpKUn/IW+tkq5U33O+8+iE4S3TapOA6lHAcOBgQ9NsW9W0WnuKDsKOT6DoAARHQMoEVVLZg5u8rVFWVhalpaUALptKq3mgF/tDhx5wmwH73DRCfdNqPb6FAaqVAcxydxxeJzwZRj7u7ii0Buj1DxqgWhnbF6j9MRqa5egmXrMOQ9M0rUVLulJVjC6oPSnUc+iEoWma5gkSRwBCjWN4KJ0wNE3TPEFQGESlwv6v3R2JXTphaJqmeYru16h9vs+ccHckNumEoWma5im6jQMk7Fvq7khs0glD0zTNU3TqA+1jYN9X7o7EJp0wNE3TPIUQ0G0sHPwGLl1wdzR16IShaZrmSbpfA5fOqeq1HkYnDE1rhXJycpg0aRJZWVkNn6y5VsIVENAedn3m7kjq0AlD01qZrKwscnJyyMnR9Ts9kl8g9LgOdi/yuG4pnTA0rZXJzMwkMzOTsLAwd4ei2dNnIlws97g1GZ5ZsETTmmvJk45vV+tsnfrANc87fPqCBQuYPn06SUlJl/emmDRpEjk5OcydO1dXsG2NEkaoPWu2z4ee17s7mst0wtA0N6uuUDtv3rzLe2FMmzaN9PR0u3tjNHWLVs1L+PhCr5th09twoQwCQ9wdEaAThtZSNeITvieYOHEiDz74IKWlpZjN5svf7WnqFq2aF+k7GdbPVhVsM+53dzSAHsPQNI8xefJk5syZQ2lpKUlJSe4OR3O36AGqe3Pjm2rrVg+gWxia5iGmT5/OmDFjSEpKanAjJd0l1QoIARkPwKJH4fA6tfmZm+mEoWkeIikpCbPZTHFxcYPn6i6pVqLPJFj2Z1j9suMJo7ICfP0NCUd3SWmaB5k2bRqTJ0829B7Z2dnMmjWLjRs3MnPmTObMqXe3Y82d/INh0MOwZzEc2+HYaz6YAl88Zkg4uoWhaR4kLCys3sFuZ0hLSyMtLY0nnnjC0PtoTjJoGqx5Fb57Dm79b/3nHvxW1aHqMsaQUHQLQ9PcbNq0aZdXX+sxB62OoDAY/mvVyjhQTykXiwWWPwMhsZB+nyGh6IShaW42adIkSktLyc7O1rOjNNuG/hLCkuGL38H5UtvnrH0VCjbDVX9S5UUMoLukNM3NMjMz3R2C5ul8A+DGf8Pb18InD8Kt74OP34/H81bD8mchZYJav2EQ3cLQNE3zBnGDYfwLqr7UuzdBySGwVMG2+fDeLRCaANe9pKbjGkS3MDRN07xF+n3gFwQLfwUv9gOTL1gqISoVbp8PweGG3l4nDE3TNG/S71ZIGA47P4Wzp9SK8O7jVf0pg+mEoWma5m1CYtRAuIvpMQxN0zTNITphaJqmaQ7RCUPTNE1ziE4YmqZpmkN0wrBjxowZ7g7Bq+j3q3H0+9U4+v1qHKPeLyE9ZGMOZ0tPT5cbN25s8uuFELTU98YI+v1qHP1+NY5+vxqnOe+XEGKTlDLd1jHdwtA0TdMcohOGpmma5pAW2yUlhDgJ5DXjElFAgZPCaQ30+9U4+v1qHP1+NU5z3q94KWWkrQMtNmFomqZpzqW7pDRN0zSH6IShaZqmOUQnDE3TNM0hOmFomqZpDtHlzW0QQpiBqUC4lHK6jeNPADlAGICUco5LA/RgQoipwABgvvVHk4CZUsoc90XlOfSz4zj9LNXPHX+ndAujFiFEJpAJJANmG8dnAjlSygXW/wOShRATXRulx5sMLANmArP1L7iin50m0c+SDe76O6UTRi1Syiwp5QKg1M4pU63Hq80DphkemBeRUoZKKYWUcoCUMtvd8XgQ/ew0kn6WbHPX3ymdMBpBCJFm48elqEyvaXbpZ0dzFSOfNT2G0ThhQHGtn9X+d6tn7XsuRvfT16SfnSbQz1KTGPas6YTROGZ7B4QQZillqetC8VgbgdLqvmYhxHwhRHGt5nFrZLZ3QD87dulnqWnM9g4091nTXVKNU4r1k04Ntf/dqkkps2sNTG4AnnJXPB6kFP3sNIp+lpqsFIOetRbfwrDODJjSwGnFUkpHBoSKqZu9zQAt9RNiY98/IUSmlDKrxrEcwFafamvT6p6d5tLPUpMZ9qy1+IRhbb46pQkrpcwWQpTW+nEYkGXj9BahMe+fECIJWCaECK31YLb6qZCt8dlpDv0sNZ2Rz5rukmq8j2rNZx4DzHZXMJ7E2n0wvdYv+BTUHHpNPzsO089SsxnyrOny5rVYp6Rl8uOc5dlAVs054NYVlNlAEuiZGzVZPxlWP6jhwEH9/vxIPzuO08+Sfe76O6UThqZpmuYQ3SWlaZqmOUQnDE3TNM0hOmFomqZpDtEJQ9M0TXOIThiapmmaQ3TC0DRN0xyiE4amaZrmEJ0wNE3TNIfohKFpmqY5RCcMTdM0zSEtvlqtpnkKIUQmqsz0GGA6kF79bwfL62uaW+kWhqa5QPU+y9Zy8fOtX2ZUcbjJ1kJ7mubRdMLQNNeovRlQupRygbWM94O1dpbTNI+kq9VqmosJIWYCZt0NpXkb3cLQNNebiOqS0jSvohOGprmAdcC7eiwjqbp7SghhFkJMdWtwmuYgnTA0zWDWhFC9PWY6UFrj8FN6FznNW+gxDE0zmHUG1DRgA2pWVBpq28wcIFsPeGveQicMTdM0zSG6S0rTNE1ziE4YmqZpmkN0wtA0TdMcohOGpmma5hCdMDRN0zSH6IShaZqmOUQnDE3TNM0hOmFomqZpDtEJQ9M0TXPI/wP0lf1t9RF1oQAAAABJRU5ErkJggg==\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",
    "\n",
    "x = np.linspace(-10, 10, 500)\n",
    "__=ax.plot(x, my_vector_function_deriv(x, 5)[0], label=r'${\\rm y=5}$')\n",
    "__=ax.plot(x, my_vector_function_deriv(x, 1)[0], label=r'${\\rm y=1}$')\n",
    "leg = ax.legend()\n",
    "xlabel = ax.set_xlabel(r'$x$')\n",
    "ylabel = ax.set_ylabel(r'$\\partial f/\\partial x$')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "sapphire-microwave",
   "metadata": {},
   "source": [
    "### Now for a more physically interesting example\n",
    "\n",
    "It turns out that the mass assembly history of dark matter halos is reasonably simple to approximate: \n",
    "\n",
    "$$M_{\\rm halo}(t)\\propto (t/t_0)^{\\alpha(t)},$$ \n",
    "where $t_0$ is the present-day cosmic time in Gyr, and $\\alpha(t)$ is a simple sigmoid function that smoothly transitions the power-law growth from a period of rapid early-time growth to a period of slower late-time growth. The next cell gives a JAX-based implementation of this function, and then shows a plot of the typical growth history of a dark matter halo with Milky Way mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "hydraulic-bridges",
   "metadata": {},
   "outputs": [],
   "source": [
    "@jjit\n",
    "def _sigmoid(x, x0, k, ymin, ymax):\n",
    "    \"\"\"Basic sigmoid function with asymptotic values ymin and ymax\"\"\"\n",
    "    height_diff = ymax - ymin\n",
    "    return ymin + height_diff / (1 + jnp.exp(-k * (x - x0)))\n",
    "\n",
    "@jjit\n",
    "def _rolling_plaw_vs_t(t, t0, logmp, x0, k, early, late):\n",
    "    \"\"\"Kernel of the rolling power-law between halo mass and time.\"\"\"\n",
    "    logt = jnp.log10(t)\n",
    "    logt0 = jnp.log10(t0)\n",
    "    rolling_index = _sigmoid(logt, x0, k, early, late)\n",
    "    log_halo_mass = rolling_index * (logt - logt0) + logmp\n",
    "    return log_halo_mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "detected-triple",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEjCAYAAAAYFIcqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuDklEQVR4nO3deXxU9bk/8M+XJexkSCBAQJaAgIAsIWFxq0pQ22pbNbgrLhD02l67/Aza2xbqvdaG1nq73GoQQdzZrK21WBPFHYEkrLJJEtaENZmwBMj2/P4434HJZLbMnOTMzPm8X695wTln5pxnJjPnOee7KhEBERFRqNpYHQAREUU3JhIiIgoLEwkREYWFiYSIiMLCREJERGFhIrERpVS2UqpYKVUY4Hk5SinR/6Z42Z6qlMpTSuV5LC9vobizdNyVSqkst/UOt1jzlFKpHq9brrfltlBcjT4Hq/ZB5uHfIzRMJDYiIvMB5ABI8ZYgPJSIyBwRKQHOJ5dcvZ8iAHPc9luk99siRGSB3n+J/r9rvVNE5gAoArBcx+H+uukA5ojI7BaKq9HnYNU+KHRKKYf7Mv8eoWEisZ8KAMsAZHrbqK/q13vZtBSA+5W908t+W9IyAKmeP3y3Y0/3XKmTZX4Lx2UGp9UB2NhtVgcQC5hI7CkXgK+r9AR4ObGJSJHnFX9rEhEngBJ4/PB1YskDkOHlZalWxkxRYZrVAcQCJhIb0ifXBM86hSiwAk3vPNJ0kZ1TKeX1LovIG6VUDgCH1XHEgnZWB0CWWQDgdhj1CwAApVSGiOQrpRpd3esiIlf9SFBXcPpHmg1ggYjM1pXkc2Dc7Ux3q3vJg3EXNCuIu4elep/uHPrfZfr9rHDb5vSIyZVoUgAUiUi+Xh92bDopJ+h40nXdTcBj+9lfFi4UFya41w0FiCFHv49nPONx+7umwqhvWuHxep8x6nhK9KLDFZOv9UHE6np/KbjwHZwO43v2IoACAMv19mm6vsvn56LX5+jXTQeQhgv1dtNFpMTtOXP0Phww6gtd36kF+s7XFaPfvym5ERE+bPQAkKn/TQVQ6bEtw/UvgGKPbakA8tyWUzyWU2FUeLuWswCkeOwjC0CuZzwAHM2Iv9IVp8f7yXB/P+7P0cvZHst57vGFGpt+38Ue+1ruEWOgY3t+lss9tjvctweIJwNAoZd4cjw/x2A/H8+/pY4nx9f6IGLMdfuuObx8186/B709O5jPRceU7bEfz79plq/vdHP+pnw0frBoy6bEuMKu8Cje8ldh7gx23/rKb5noK3u3Yy6Al8pNcbsKDMIy6OItfYWdr/eRD8Dh9n4cHq9L9yj6KoJbvUq4sXm81xIYJ8Ggju1Ox5/ivj8dg+uKOpAKGMnPM54mz/NouBAoxvN1ajqepQHW+3MbjDuH85+vl++hQ0RKxGiZNz/IzyUXxl2pi8P9Peh9BN34IsDflNwwkdjb+Up3XawVbsW0w1VM5OcEnO86YekTma/n+bIcF074Do/jrEDjE8l5IjJddFGOLqpzoGmyCTU2v88L8tguafB+4i8GMCGMeI6HGqNOsilu/XWyxGh84XV9EPGVwCgycklA0/fsuRzwc9EnfodHgnQvqk3xSA7+OIN8HoGJxO5WwNzmj64y61w/nQCfwYWr2AwJUFfgycudh7s8AJnerjyVUilKqVy3q1enibH5bfoc5LFdHH62JfjZFigef8f0G6NSyiFGHUUPGMVH0/Vzva4PIr5cALN1h9IsAM94ufDwXHb42Z/757ICQJZOhkXw30KxEY++VS3dnD2mMJHYmL46K9FX4Wb8cEr0fucDSPOstNfbihBch0h/8gE8CbeGAtoyGMUPGV5OTIUw7pQWeF6Vuq5gTYrNm4DHdpMP70UoQ2AkypbiM0YAvwbOdwDNF6PBRQqMv4G39YFUwEjaaTCKQOcH8ZpgPxdX8VaqLhorgtH/KBVNvy+eoq0VY8RgIrGfdI/lpQCeNKFYy9MsNO7A6C5XP0I95nLoE4X7Sp08muzTdffikVwcepsDja9oQ4nN252Ca//NObYrmZW433Hp56VJEK2hAsXjTRAxdvRSP+Nq2eZ1fQDT9PHy/RSBNoo32M/FVbyFxp9BPozvuLfiM/fk5B5Lsz5D27O6tp+P1nvAKH4QGCdKh17nQOOWLpkwTtSVMJraOmD82FzrsjyWs6FbbLm26/2k6mMVomkLKgeCbIXk43044KN1ELy0FnN779kwKl9T9T5yoVt9hRKbx/t2tSzKhFFuf/59+zu252fptu9svT3T9XcwMZ5st+9BShAxuta7Hq64va4PIs4Mt5iK9SPH23vwfN/BfC56fYrH55LlI5Zs/Z3JbM5nyEfjh9IfFFGrUkplikc/hkgRybFFO31HkSEexVm6L0eisK9GVGLRFrU63UIsIk/UkRxbjJiNxp1GAZyvV2MdRZRiIqFWoVv5+OrjYalIji0GeR0XTTfM4NDtUYpFW9Qq3FpwOSLtij+SY4tF+vNOQeOhVcDPPnoxkRARUVhYtEVERGGx5ei/PXv2lEGDBlkdBhFRVCksLDwmIr0819sykQwaNAgFBQVWh0FEFFWUUnu9rY/Koi09Rk+O53hLSqks/chtgWEuiIjIi6hMJDDG6GmUKHRSKRBjuITl8D08BxERmSgqE4kYo7J6DjKYggujfBbASDZERNTCLKsj0QOuZcHHsAh6yITz8xZIgAHrdBt0Vzv0NOiJc4iIqGVZckeiOyRlwBgC2uFlew70nNI6gQzxmL0tkNnQs+gREVHLsiSRiDF89Ar4nmwny6OX61IEPzlNoBn6iIjIRBFXR+Jj5jsnfMxx7fHaDAD5IlLibVIlIiIyXyT2I0lA04r0Rss6SaQBcCqlICJFOgEtB1ChlAKMiYmaNY0rEVEsqjxdg892H8Onu47ilzeORHyn9qbuPxITicPXBj1HtFO32prgvk2MGdR6BHOAsrIy6GQDAJg7dy7mzZsXUrBERJGmrr4BG/c78cmuo/h011FsPlgFESC+U3vcPWkAxg8I6lQZtEhMJE40nebS27SXIUtOTkZZWZmZuyQistSBymp8usu46/ii+BhOnq1DGwWMu8iBH08dhquG9cSY/g60baMC76yZIjGRVKDpXYkDaDKnNBGRbdXUNaBgTwVW7zyC1TuPYveRUwCA5PiO+O6lfXHVsF64fEhPxHc2txjLm4hLJLq+w+mxOgGs7yAimztUdRYf7zyC1TuP4PNvjuF0TT3i2rbBpJQE3JF+Ea4e3gtDenVtVHTfGiIukWjLPObNngYOeUJENuOq61i98wg+2nEU28tPADDuOr4/vh+uHZ6Ey4YmonOctadySya20i2sMnChb0gujGa7RW7PyYbR8ioFCNyzvTmSk5OlvLyclexEFHGOnzqHT3YdxeqdRkV51ZlatG2jkDawB64ZkYRrhidhWO/Wv+sAAKVUoYg0GX7KljMkpqWlCYeRJ6JIUXL0FPK2HUb+9sMo3FuJBgF6du2Aq4f3wrUjknD50J6mN9kNha9EEqlFW0REMau+QbBhXyXyth1G3vbDKDl6GgAwsm93/OjaizH1kiSMTo5HmxZoYdUSmEiIiFpBdU0dPvvmGPK3HcZHO47g+OkatG+rMDklEfdfNghTL+mNfo5OVocZEiYSIqIWcuTkWXy4/Qjytx3G57uP4VxdA7p1bIdrRyQh45Le+NbwXuje0foiq3DZMpG4erazsp2IzLa/ohqrtpZj1dZD2LDPCQDo36MT7pw4ANeN7I30wQlo3zbihjkMiy0TCXu2E5GZSo6ewqqth7Bqazm2HjSa6I7u1x0/mzYM00b1xvDe3SxpZdVabJlIiIjCISL45sgp/GtLOd7fegg7Dp0EAIwf4MDPvzMC3x7dFxcldLY4ytbDREJEFAQRwddlJ/D+1kP419ZylBw9DaWA9IEJmHvTSFw/qg+So7SyPFxMJEREfuw+chL/2FiGdzeXo/TYabRRwOSURDxw+WBcP6o3krp1tDpEyzGREBF5OFBZjXc3leMfm8qwvfwElAKmpCQi66oUXDeyNxK7drA6xIhiy0TCVltE5OnoyXP41xYjeRTurQRg1Hn86saRuHFMXyR1552HL7ZMJGy1RUQAUHWmFv/eegjvbi7DF7uPoUGAEX264fHrh+N7Y5NtVWEeDlsmEiKyr7r6Bnz6zVGsLDyIvO2HUVPXgAEJnfEfVw/F98YlY1jvblaHGHWYSIjIFraVncDbRQfwzsYyHDt1Dj06t8ddEwfgB+P7YWz/+Jju59HSmEiIKGYdPXkOf994ECuLDmJ7+Qm0b6tw7Ygk3JraH1cPT0Jcu9jqYW4VJhIiiim19Q34cPthLCs4gE92HUV9g2Bs/3g89f1RuGlMMnp0ibM6xJjDREJEMaH02Gm8tX4fVhYewLFTNejdvQNmXZmCW1P74WLWe7QoWyYSNv8lig1na+vx768P4c11+/BVSQXatjGKru6ceBGuurgX2sXY4IiRypaJhM1/iaLbzkMn8ea6ffjbhoOoOlOLAQmd8fj1w5E5oT96s79Hq7NlIiGi6HO2th7/3FyO19fuxYZ9TsS1bYPrR/fBHekXYUpKYtTMJhiLmEiIKKLtr6jGa2v3Ytn6/aisrsXQpK74xXcvwS2p/ZHAivOIwERCRBGnoUHw+e5jeGXNHny44wjaKIVpl/TGfZcNxJSURPb5iDBMJEQUMarO1GJF4QG89tVelB47jZ5d4/Do1UNx16QBth2iPRowkRCR5UqOnsKiL0qxsvAgztTWI3WAAz++YxxuGN0HHdq1tTo8CoCJhIgsISJYW1qBhZ+V4sMdh9G+TRt8f1wyZlw2CKP7xVsdHjWDLRMJ+5EQWae2vgH/2lKOhZ+VYsvBKiR0icOPrr0Y904eiF7dOM9HNLJlImE/EqLWV3WmFm+t24eXv9yD8qqzSOnVBb+5+VLcktoPHduz+Cqa2TKREFHrKXOewcLPSvHW+n2orqnHlJREPH3zaFw9LIl9P2IEEwkRtYjSY6fxwsfFeHvDATQI8L2xyXjoisGs/4hBTCREZKqvy6rw14+LsWpLOdq3bYM7Jw7ArCtTONtgDGMiISJTrN9Tgb+u3o3VO4+ia4d2mP2tIXjw8sGsQLcBJhIiCpmIYE3xcfxv/jdYt6cCCV3i8Pj1w3HP5IGI79Te6vColTCREFFI1hQfx3P5u7CutAK9u3fAr24ciTsnDkCnOLbAshsmEiJqlrUlRgL5qqQCSd06YN5NI3HHxAFswmtjTCREFJSCPRV4Ln8Xvth9HL26GXcgd01iAiGbJhL2bCcK3sb9Tjz7wU589s0x9OzaAb+8cSTuZgIhN7ZMJOzZThTY7iOn8Pt/78T7Xx9CYpc4/OK7l+DuSQNZB0JN2DKREJFv5VVn8L9532B54X50jmuHn2QMw8wrB6NLB54uyDt+M4gIAOCsrsFfPy7Gy1/uAQS4/7LBePSaIUjsyn4g5B8TCZHNnampx6IvSvHCJ8U4da4ON4/vh59kDGNPdAoaEwmRTTU0CP6xqQw57+9AedVZZFyShMevH4HhfbpZHRpFGSYSIhsq2FOB/35vOzbtd2J0v+547vZxmJySaHVYFKWYSIhsZH9FNX77/g68t7kcvbt3wO+nj8Ut4/txOHcKCxMJkQ2cPFuL/1tdjEWfl6JNG+CxqRdj9rdS0DmOpwAKH79FRDGsoUGwrGA/fvfvnTh+uga3pPZD9vUj0Ce+o9WhUQxhIiGKUZsPOPHLv3+NTfudSBvYA4sfSMeY/g6rw6IYxERCFGMqTtfgd//egbfW70fPrh3w3O1j8YNx/aAU60GoZTCREMWI+gbBm+v24fcf7MTJs3V46PLBeCzjYnTryHlBqGW1sToAK7gGbeSAjRQrivZV4vv/9zl+8c5WjOjTDaseuxK/uHEkkwi1ClvekXDQRooVVdW1+O372/Hmuv3o3b0D/nTneNw0pi+LsahV2TKREEU7EcG7m8vx1LvbUFldg1lXDsZjGcPQlQMrkgX4rSOKMvsrqvHLv2/FxzuPYkz/eCx5MB2jkuOtDotsjImEKErU1Tdg0ReleC7vG7RRwK9uHIkZlw1CW/ZKJ4sxkRBFgU37nXjy7S3YVn4CGZf0xlPfH4VkRyerwyIC4CeRKKVuATARgIR5jPUi8naY+yCypTM19Xj2g51Y9EUpenXrgBfuScX1o/qwMp0iir87kmkAnjDhGAsAMJEQNdO60gpkr9iEPcercc/kAci+YQS6szkvRSB/iaRERKrCPYBSqiDcfRDZSXVNHea/vxNL1uxB/x6d8MasSbhsSE+rwyLyyWciEZHfmXEAs/ZDZAdflRxH9orN2FdRjfsvG4THrx/OudIp4pnyDVVKDQaQCqASQA8AwnoRouCdPleHnPd34JU1ezEwsTPeyprMiaYoaoSdSJRS4wBARFZ6rL8VQJ6InAj3GESx7MviY8hesRkHnWfwwOXGXQjnCaFoYsa3dYhnEgGMxKJbfvHOhMiLs7X1+P2/d2Lh56UYlNgZy2ZPQfqgBKvDImq2sBKJUmoQgCK35f8HoEhEPtKrwq6sJ4pF28pO4MdLN2DX4VO4d/JAPPmdEbwLoagV7je3B4Bit+U7ACQCcCWScPugEMWU+gbBi5+V4NkPdsLROQ6LH0jHNcOTrA6LKCxhJRIR2aCUuhY6cYhImsdTeoSzf6JYsr+iGj9btgnr9lTg26P74OmbL0VClzirwyIKmxn30l6ThW7JVWLC/omimohgReEB/PrdbVAA/nDbWNw8njMWUuzwm0iUUt0DtbrSleq3AqgUkY+UUvEAMgAUi8gGE2MlijrO6ho8+fYWrNp6CBMHJ+APt41F/x6drQ6LyFSBZkh8MpidiMhKnUSmAhislzeGHV0L4QyJ1BrWlVbgO3/8DPnbD+PJb4/Am7MmM4lQTApUtJWplPpARFYHszMR+dCEmFocZ0ikllTfIPjLR7vxxw93YUBCZ6x85DKM6e+wOiyiFhMokfQEkK/LcosA5APIA1DgrcgrmKIwolhWXnUGj721EetKK3Dz+H747x+M5qyFFPMCfcNTYdR3VABIgTEi8BwAopQqwYXEUiQie2AUhQVVHEYUaz74+hCyV25GTV0D/nDbWNyS2t/qkIhahd9EIiKlAF7UdR8lInIdACilUgFMgJFYFgKIV0o59cuYSMhWztbW4zf/2o5X1uzF6H7d8ec7UzG4ZxerwyJqNUHdc7vqPnTrrGIRKYJR1PWiXj8YRlJ5vIXiJIpIxUdP4dHXi7Dj0EnMvGIwHr9hODq0a2t1WEStqlmFt7qp72A9hla+qz5E37ksUEqltESQRJHoH5vK8OTKzejQvi0W35+Oa0awhzrZU7NrAXXSKFVKTVVKxXsMF59rXmhEkelcXT2efs8oypowsAf+ctd49I3n/OlkXyE3J/Eo7qoUkY90kiGKWfsrqvHDN4qw6UAVZl4xGHO+PQLt2wbqjkUU28Id/bc7jEEb05RSswHM0a23iGLOh9sP46fLNqGhQfDCPRNww+g+VodEFBECDpECo9mv6zHEYxkAFAAnjHG1ZoOttijG1NU34Nm8XXj+42KM7Nsdz9+TioGJbJVF5BLojsQJYyh4J4BSGMliA4AV+v8lLM6iWHbkxFn86M0NWFtagTsnXoS5N41Cx/ZslUXkLlAiKYFRge6E0UqLSYNsY/2eCjzyWhFOn6tjB0MiPwIlkhUi8jsAUEqNV0rNdNtW4jYTIvRzbvFoxUUUdUQEr361F0+9uw0XJXTG6zMnYXifblaHRRSxAvVsf8Lt/xtgFGsBMDohKqVmAXDAKP6qApANztFOUexsbT3+629bsbLoAKaOSMIfbh+H+E7trQ6LKKKF0/y3FLpnO3B+2JREM4IissJB5xk8/GohthyswmNTL8ZjUy9GmzacfIooENOGJRWRIqXUArP2R9Saviw+hh++sQG1dQ148b40TBvZ2+qQiKKGz55USqlBzd2Ze1FYOPshai0igoWfleDel9ahR+f2eOeHlzOJEDWTvy65c0w6Ro5J+yEy1Zmaevx46Ub8z3vbkXFJEt559HIM6dXV6rCIoo6/oi2llHomzP0r/SCKKAcqq5H1SiG2HzqBx68fjke+NYT1IUQh8pdIzLojIYooBXsq8PBrhThX24BFMzhqL1G4fCYSEalqzUCIWsPygv34+d+2oJ+jE97KSsfQJBZlEYWLk0mTLdQ3CH67ajte/KwUVwztib/cNR6OznFWh0UUE6Jy/GullEMplaP7rgRcT/Z24mwtHlqyHi9+Vor7LxuElx9IZxIhMlG03pGk4cLow8GsJ5vac+w0HlqyHnuPV+Ppm0fj7kkDrQ6JKOZE5R2JiOQDqAh2PdnTl7uP4fv/9wWOn67Bqw9NYhIhaiGW3ZEopRwAsgAkikiTFmJKqWwYow8nAICIsNc8Be3VNXsw791tGNKrCxbel44BiZ2tDokoZlmSSJRSGTAGexziY3sOgPUissK1rJTKdC0T+VJX34Cn/rkNr6zZi6kjkvC/d4xDt44cdJGoJVlStCUi+TopOH08JcsjaSyFMfsikU+nztVh5isFeGXNXmRdlYIF96UxiRC1AtPuSJRSUwFk6sW8UOcl8dHiygkgI8TQyAbKnGfw4Mvr8c2RU/jNzZfirkkDrA6JyDZMSSRKqVthtJZy3UWkKqVmisjCEHaXgKYV5o2WddFYGgCnUgoiUuRvPcW2rQer8ODL61FdU4/F96fjqmG9rA6JyFbMuiNxumZS1D7UdyihcPjaoJRyiIhTt86a4Lnd13pPZWVlUOrCuEpz587FvHnzQgqWrJW37TD+880NSOgSh5WPcCZDIiuYlUjEy7r4EPflhG6p5cZzOSzJyckoKyszc5fUykQEi77Yg/95bxvG9IvHizPSkNSto9VhEdmSWYmkhx4peL1eTgdQHOK+KtD0rsQBACLiDHGfFEPcW2bdMKoPnrt9HDrFtbU6LCLbMiWRiMhKpVQJgNv1qqV6jvdQ9lWklHJ6rE4AkB9GiBQjTp2rww/fKMLHO49i9lUpmHPDCA7/TmSxkJv/KqUGuT8AVAJ4QT8qw5zLZJlSKtNteRqA3DD214irjoT1ItGlvOoMMp//Ep99cwy/uflSPPmdS5hEiCJAOHck+QAK0XjiKlddiQIwHsCT3l6om/hmQDcX1r3Y812trERktlIqW7fCSgFQbGZnRNaRRJ/t5Sdw/+J1qD7HlllEkSacRDJbRD70tVEpNd7XNp0wigDM9/Mcn9vIXr7YfQwPv1qILh3aYfkjUzCiT3erQyIiNyEXbQVIIuMADA5130Qu72w4iPsXr0OyoxP+9uhlTCJEEcisDomDYdRh9MCFoq48ACH1bicSETz/STHmv78Tk1MSkHtvGuI7cbgTokhkVvPf2bgwFparFzrnBaGQ1DcI5v3ja7z61V7cNDYZv58+Bh3asXkvUaQya9DGpSJSKiKlAKbq+d57mLRv07HVVuQ6W1uPR14rxKtf7cXsq1Lwx9vHMYkQRTiz7kgSlFLrRSQdwESlVDqMToQfmbR/U7HVVmSqOF2DmUvWY8N+J+bdNBL3X85qNqJoYFaHxA9h9GaHiDyhB3HkgIkUtH3HqzFj8TocdJ7BX+9Kxbcv7Wt1SEQUpBaZ2Er3dB/UEvum2LP5gBMPvrwetfWCN2ZOQtogU4dWI6IWZlarrUEApuPC4IoKwFTouxQiX1bvPIJHXy9Cj85xeCtrIoYmdbU6JCJqJrMq25+AMb96vn7kwej1HpFY2R4Zlhfsx8wlBRjcswv+9h+XMYkQRSmzirZyPQdpVEoVmLRv07Gy3Xq5nxTjmVU7cMXQnnjh3gno2qFFSlmJqBWE/OtVSrl3Me6hlLoWxl2Jy2z4GGuL7KuhQfDb93dgwacluHFMXzx721g27yWKcn4TiVKqu4ic8LG5CE0HbXTnc9BGsqfa+gbMWbkZbxcdxIwpAzH3plEcvZcoBgS6I3kSvpNByIM2kv2cqanHo28U4aMdR/DTacPwo2uHNprumIiiV6DK9kyl1DXeNvhLInp7SBNbUexxVtfgnpfWYvXOI/ifH4zGf069mEmEKIYEuiPpCSBf/+iLcKFFVoG3Iq8ARWFkQ+VVZzBj0TrsOVbNjoZEMSrQHUkqgIcB3AZgGYAJMJJJpVLqG6XU80qpW9w6H0ZFnQib/7aO3UdOIfP5NShznsXLD6YziRDFKCUigZ+k1FQADhFZqZdTYSSVaTBmOowH4AQAEUlsqWDNkpaWJgUFEds6OSZs3O/EA4vXoW0bhZcfmIjR/eKtDomIwqSUKhSRNM/1QTX/ddWH6DG0it1mOHxRrx8MI6k8blrEFLU+3XUUD79WiMSucXj1wUkY1LOL1SERUQtqVs92fUdSpYuzurutLxWRBQBWmh0gRZe/bzyIh5asx8DELlj58GVMIkQ20OwOiXrOkVKl1FSlVLyIuM+CmGteaBRtXv6iFL/+5zakD0rAwhlp6N6RMxoS2UHIPds9irsqReQjnWTIZkQEz+Xtwp8+2o3rRvbGn+4cj47t2VudyC7CGrRRF28VA0hRSi3l0PH209Ag+PW72/Cnj3bj9rSL8Ne7U5lEiGwm4BApMOZedz2GeCwDxhApThjjbHF8LRupq2/AE29vwYrCA5h15WD8/DuXsKMhkQ0FKtpyAhD9bymMZLEBwAr9/5JoLM5y9SOZO3cu+5KE6FxdPX781kas2nqIQ54Q2VygRFICowLdCSA/GpOGNxxGPjxnauox+7VCfLrrKH5140g8eAXnVieys0CJZIWI/A4wBmFUSs1021YiIh+5P1kpdYtHKy6KMSfO1uKhl9ejcG8l5t86BrelX2R1SERkMb+JRESecPv/BhjFWgCMTohKqVkAHDCKv6oAZANgIolRFadrcN+itdh56CT+fGcqvjuGQ54QUXjNf0uhe7YD54dNifjhUSg0h6rO4p6X1mJ/RTUW3JeGa4YnWR0SEUUI0+Y3FZEipdQCs/ZHkWPf8Wrc/dJXqDxdiyUPTsTkFF4vENEFpk6U7V4URrHhm8MncffCtaipb8AbsyZhTH+H1SERUYQxNZFQbNlyoAr3LVqL9m3bYGnWFAzv083qkIgoAjGRkFdrS47joSUFcHRuj9dnTsLARA6+SETehTVECsWmj3cewX2L1qF39w5Y/vAUJhEi8suWiYQzJPr2ry3lmPVKAYYmdcWy2VPQN76T1SERUYSzZdEWe7Z7t6xgP55YuRmpA3rgpfvTEd+Jw8ATUWC2TCTU1OIvSvHrd7fhyot7IvfeCegcx68GEQWHZwubExH85aPdeDZvF64fZcwl0qEdh4EnouAxkdiYiOCZVTuw4NMS3JLaD/NvHYN2bW1ZbUZEYWAisan6BsEv3tmKN9ftw31TBmLeTaPQpg2HgSei5mMisaHa+gb8bNkm/GNTGR69Zgj+33XDOZcIEYWMicRmztbW44dvFCF/+xHMuWEEHrl6iNUhEVGUYyKxkVPn6jBrSQG+Kj2O//7BaNw7eaDVIRFRDGAisQlndQ3uX7weWw5W4Q+3jcXN4/tbHRIRxQgmEhs4cvIs7ntpHUqOnsbzd6fiulF9rA6JiGIIE0mMO1BZjXsWrsXhE+ew6P50XHFxT6tDIqIYw0QSw0qOnsI9C9fi5Lk6vDZzIiYMTLA6JCKKQbbsfWaHQRu3lZ3AbblrcK6uAW9lTWYSIaIWY8s7klgftLFwbyUeWLwOXTq0w2szJ2FIr65Wh0REMcyWiSSWfbH7GGa9UoCkbh3w2sxJ6N+js9UhEVGMYyKJIXnbDuPR14swuGcXvDpzIpK6dbQ6JCKyASaSGPH3jQfx02WbMLpfPJY8kA5H5zirQyIim2AiiQGvr92LX7yzFZMGJ2DhjHR07cA/KxG1Hp5xotwLnxTjt6t24NoRSfjr3ano2J5ziRBR62IiiVIigmc/2IW/rN6NG8f0xXO3j0N7ziVCRBZgIolCDQ2Cp/65DS9/uQd3pF+Ep2++FG05lwgRWYSJJMrU1TdgzsotWFl0ADOvGIz/+u4lnEuEiCzFRBJFztXV48dvbcSqrYfwk4xh+M+pQ5lEiMhyTCRRorqmDrNfLcRn3xzDL28ciYeuGGx1SEREAJhIosKJs7V4cPF6FO6rRM6tl+L29AFWh0REdB4TSYQ7fuocZixehx3lJ/HnO8fjxjHJVodERNQIE0kEO1R1Fncv/AoHKs/gxfvScM2IJKtDIiJqgokkQu07Xo27X/oKFadqsOTBiZickmh1SEREXjGRRKBdh0/inoVrUVPfgDdmTcbYixxWh0RE5BMTSYTZfMCJGYvWoV3bNliaNQXD+3SzOiQiIr9sOaZGpM6QuLbkOO56cS06x7XD8tlMIkQUHWx5RxKJMySu3nkED79aiP49OuG1mZPQN76T1SEREQXFlokk0ry3uRw/XroBFyd1wysPTUTPrh2sDomIKGhMJBZbVrAfT6zcjPEDemDR/emI79Te6pCIiJqFicRCiz4vxVP/3IYrL+6J3HsnoHMc/xxEFH145rKAiOBPH+7Gc/m7cN3I3vjzXePRoR0npCKi6MRE0spEBE+/tx0LPy/FLan9MP/WMWjHCamIKIoxkbSi+gbBz9/egqUF+zFjykDMvWkU2nBCKiKKckwkraSmrgE/WboR720px4+uHYqfThvGuUSIKCYwkbSCMzX1eOT1Qny88yh+/p0RyLpqiNUhERGZhomkhZ08W4uHXi7A+r0VeOaWS3HnRM4lQkSxhYmkBVWcrsGMReuwvfwE/njHeHxvLOcSIaLYw0TSQg5VncU9L63F/opqLLhvAq4d0dvqkIiIWgQTSQvYe/w07l64FpWnOZcIEcU+JhKT7Tx0Eve+xLlEiMg+mEhMtGm/EzMWr0Nc2zZYNnsKhvXmMPBEFPuYSEyypvg4Zi5Zj4SucXj9ockYkNjZ6pCIiFoFE4kJPtx+GI+8XoSBCZ3x6kOT0Ce+o9UhERG1GiaSMP1940H8bNkmXNK3O5Y8OBEJXeKsDomIqFUxkYThjbX78F/vbEH6oAS8NCMN3TpyLhEish8mkhDlflKMZ1btwDXDe+H5eyagY3sOA09E9sRE0kwigmc/2IW/rN6NG8f0xR9uG4e4dhwGnojsi4mkGRoaBL9+92ssWbMXd6RfhKdvvhRtOQw8EdlcVF5KK6UcSqkcpVSqx/ospVSG/jfF7OP+/G9bsGTNXsy6cjCeuYVJhIgIiN47kjQAjRKFThxDRGSBXl4OYLqZB718aE/0c3TCD68dyrlEiIi0qEwkIpKvlPJMEpkAit2WU2Gymzh6LxFRE5YlEqWUA0AWgEQRmeNlezaAEgAJAOC60/AjUT+fiIhakSV1JEqpDAAZAIYAcHjZngOgRERW6AQyRCmV2bpREhFRMCxJJCKSLyIrADh9PCVLb3dZCmB2gN0eh757ISKi1hNxrbY8W2JpThh3MP6sgHGH41JkVkxERORbxCUSGHcVFR7rGi3rorE0ALe7Eo+IlABY72r+C6BJvQuFb968eVaHQBQTYum3pETEuoMbdSEOEZntti4TQI6IDHFb5wBQCaCHiDjDPW5ycrKUl5efX547d25M/VFbklIKVn5niGJFNP6WlFKFIpLmuT4S70icaFrXYWrdR3JyMkTk/KM5ScSMhBPqPpr7utZ+X7EuEj+j1ozJ7GPxtxQ7IvGOJBVAoYgof+vCkZaWJgUFBSG91oyriFD30dzXNef5wT43Gq+izBKJ7701YzL7WPwtRd73KRBfdyQR1yFRRIqUUk6P1QkA8s06RmFh4TGl1N4QX56slCoLM4RQ99Hc1zXn+cE+14z3H60i8b23ZkxmH4u/pcj7PgUy0NvKiEsk2jKlVKZbE+BpAHLN2rmI9DJrX0REdmdJ0ZYuqsrAhb4huQDyRaTI7TnZMJrwpgBB9Wwni+hWciUw/lb5ugUdETWTblj0JICl7ufDSGfJHYn+gIoAzPfzHJ/bKHK0xmCZRDbSZEDaaBCJrbYourT4YJlEdiEi+Wjajy7iRWodCbWyMAbR5GCZRG5aYEDaiMc7EuIgmkQmsetviYmEwh1Ek4NlEmktNCBtxGMiIb+CGESTg2USBSGMAWkjHutIKBC/g2iKSIlSar2+pU8BB8sk8qU5A9I6dc/3qLgwYyKhQBy+NiilHCLi9LhVJyLvHL42uP2W8gFMaL2QzMGiLQrEiRYeRJPIJpyI0d8SEwkFUoGmV1IOADBjSH8iG4nZ3xITCfmly2idHqtNHUSTyA5i+bfERELBWObR1t3UQTSJbCQmf0uWzkdCkYGDaBKZw66/JSYSIiIKC4u2iIgoLEwkREQUFiYSIiIKCxMJERGFhYmEiIjCwkRCZGNKKYdSKlc/WnWKV6VUtj5u1M/HYXcctJGInH5m8kuEMeeME0CFiKzQJ/78cIf1EJH5ejbBJ2FMR0BRiomEyA89o51DRMKefMg1wmtL7d8s+uS+HECuiMx3X6+UyoIxVUDUjVBLLYeJhMi/pSbu6zYAnr2Yzdy/WT4EMEcPaX6eiDiVUssQA0N6kLmYSIj8MHlioWnwSCSRNnGRvuOo8EwiLjqZsBiKGmFlO1ErcBVhWR1HEGbDKNbyJxLvoshCvCOhiOS6MtaLCe4D2wWxrUQvOmAM010AIAdGhfEzep0DQLqIzNHTmwJAKoAS14yPuhVTLgCIyLRgju/jvWTq46XoCmzAuDNJcN+/HvCv2XG6Hcc1GKAjmLh8SMWFz88r9+O61Zk4AUwXkRK9Pk/HP0s/9UUYf4flMAYrnCYi00OIjyKRiPDBR0Q9YJxMs9yWU1zL0Ccit20OAHn6/1letuXo/2cAKPTYvty13W1dpcdyqmv/gWIL8J4a7cfP/kONczmAVI84M4OIy/0zSgEg7vsJ8u+VBaNi3n1dJoxGBE3elz5mtufx+YjeB4u2KKLoFkNZ0vhqOhPAEH3FniL6qhc4P7Ncib4yBi4M3+3a5iqGqYBxYnO/2vZ25V2hY3BxBhNbEG/NF6fHcrPj1HdOqdK4vmUp3D6LYHgcszmvWwCjIYHneqfb4vn3Jcbc5PM9n0/Ri0VbFGky4HHidJ10PIqt3BUDmCAis5VSy5VSAmPWueUeJ32nl9ceNyM2kzm9rPMXZwYAp1vRF6CL0kI4dgmANBhFZOe5kiiMpJkFo2guxy355CulMsXoZ+KA9/cQUqKiyMc7EoomDj/bEnQ/jekAesAo2pmulHJvqlrh5XVO88ILXoBe5M2N0wGjziTf7bFCREK5U1oBoEndhdtdhKvjYo7HHcwzuHAHlCHeW305Q4iHogATCUWa8zPHudNXufnetsG4Ss6D0UPaddLLF6OC3MxhP/zF1lypYUdzgde4QiFGD/cUXYzoT6Nkp4vVUlp7mBWKDEwkFFH0Ve4ytzoPlwx9sipxP8npk3iaqwjLy+vcr5oTvBzSYUZsQby8BI1P9k4/z21WnPrqv8Lz5O8lzmBNA5DjYwys2+C7iCpXP3z1jXGEGA9FONaRUMTRdR3ZbnUiDtFNTkVkut7mOimnAJiq/38cRqJxndgdME6IqTDuVlKUUtlijPGUCaOi3KmUKhKRfN18NkW/JkfvIwdAmlIqS0QW+IstwHty6gEKXR3+Vuj3cH7/MJrHNjtOXYE9TceVBn23ICHOBa4Tpmt/Oa7PFUaCy9dxerMARrPeRonG7fNP07EvkDDH6aLIwjnbiWzMNWiieBm0McT9ZQaTWFvq+GQNFm0RkSmUUhnNSSIUO5hIiChkurjOVTfjsDIWsg4TCRGFYzmMptfNKtKi2MLKdiJKUUothzF0fLM6DfroLxIUXfGeDg4CGfVY2U5ERGFh0RYREYWFiYSIiMLCREJERGFhIiEiorAwkRARUViYSIiIKCz/HwcMkQgirqNAAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "TODAY = 13.8\n",
    "tarr = np.linspace(0.5, TODAY, 5000)\n",
    "MILKY_WAY_MASS = 12.0\n",
    "x0, k, early, late = 0.25, 3.0, 1.5, 0.5\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "__=ax.loglog()\n",
    "__=ax.plot(tarr, 10**_rolling_plaw_vs_t(tarr, TODAY, MILKY_WAY_MASS, x0, k, early, late))\n",
    "\n",
    "xlabel = ax.set_xlabel(r'${\\rm cosmic\\ time\\ [Gyr]}$')\n",
    "ylabel = ax.set_ylabel(r'$M_{\\rm halo}\\ [M_{\\odot}]$')\n",
    "title = ax.set_title(r'${\\rm Milky\\ Way\\ halo\\ mass\\ growth}$')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "veterinary-ocean",
   "metadata": {},
   "source": [
    "Since we have implemented our function in JAX, it is now straightforward to calculate $dM_{\\rm halo}/dt,$ the _mass accretion rate_ of our dark matter halo across time, via automatic differentiation. We just need to take care of the appropriate Jacobi factors since we have implemented our function based on logarithmic variables. The next cell shows how to do that using the usual scalar-then-vmap pattern, and plots the result (also converting to the more conventional units of $M_{\\odot}/{\\rm yr}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "appointed-butter",
   "metadata": {},
   "outputs": [],
   "source": [
    "_d_log_mh_dt = jjit(\n",
    "    vmap(grad(_rolling_plaw_vs_t, argnums=0), in_axes=(0, *[None] * 6))\n",
    ")\n",
    "\n",
    "@jjit\n",
    "def _mass_accretion_rate(t, t0, logmp, x0, k, early, late):\n",
    "    d_log_mh_dt = _d_log_mh_dt(t, t0, logmp, x0, k, early, late)\n",
    "    log_mah = _rolling_plaw_vs_t(t, t0, logmp, x0, k, early, late)\n",
    "    dmhdt = d_log_mh_dt * (10.0 ** (log_mah - 9.0)) / jnp.log10(jnp.e)\n",
    "    return dmhdt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "viral-bidding",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEjCAYAAABnxZXbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA8yUlEQVR4nO3deXxU1f3/8ddJQiBsmYQ97An7bgii4lIxLNa6YQK21m5fBVrt11YtoLUVay0GXLqqgNb2Z78qBHCHKsENwYVVdgiETcIamLAlkOX8/pgzOExmJpPJTO4sn+fjMQ+YuTN33jNJ5jPn3HPPUVprhBBCiHATZ3UAIYQQwhMpUEIIIcKSFCghhBBhSQqUEEKIsCQFSgghRFiSAhUDlFJTlFK7lFJrarlfnlJKm3/TPWzPVEotVUotdbueH6LcE03uE0qpiS6321yyLlVKZbo9Lt9smx2iXBe9D1btQ4hoJwUqBmitZwJ5QLqnwuOmSGs9VWtdBBeK1myzn7XAVJf9rjX7DQmt9Ryz/yLzf+ftdq31VGAtkG9yuD4uF5iqtZ4UolwXvQ9W7UPUj1LK5uG2vFB9sQkFT68hmkiBih3HgflAjqeNphWyysOmeYDrH6zdw35DaT6Q6eUP8TiQ636jKcIFIc4VDHarA8S48R5uc/99D3eeXkPUkAIVW2YD3loVqXj4wNRar3VvoTQkrbUdKMLtD9EUrKVAtoeHZVqZWUSMUe43WP37HoAaryGaSIGKIeYPL9X9mE0EWEDNllKW6bq0K6U8tgqF8EYplQfYrM5RH9HwGmqTYHUA0eDmABNwHL8BQCmVrbUuUEpd1BoxXWXO409+fVMzfzRTgDla60lmcMNUHK2zXJdjW0txtNru9uMb6zyzT1c28+9883oWuGyzu2VyFrB0YK3WusDcXu9sptinmjzDzLGxWp/bx/4m8m23aarrsbdaMuSZ1zHDPY/LzzUTx/G8BW6P95rR5CkyV23OTN5u9yOrz/fD1+v3tM289rnAaiDf7HeUOQ6JUmoKjt911+w55nq62Q6Ov4tUvPy+e8tV23vvx/tRW35vv7seX4PpcfD4umvLEpa01nKJgQuQY/7NBE64bct2/gvsctuWCSx1uZ7udj0Tx0AF5/WJQLrbPiYCs93zALY65D/hzOn2erJdX4/rfcz1KW7Xl7rmCzSbed273PaV75axtud2fy/z3bbbXLfXkicbWOMhT577++jv++P+szR58rzd7kfG2t6PPGCi2/sz0Y9tF167yTLF5fVnuu3f9e+gxnvr6fbafi4+3vts9/3X9rNzy1/b++XtNXh93ZF2kS6+GKMdLYLjbt18vgY62P3dt/mWOV+blojLc87Bw8Fcbb7t+Wk+ppvPtAgKzD4KAJvL67G5PW6YWxfgWlyOW9U3m9trLcLxIePXc7sy+dNd92cyFLkOsffhOI6i6p6nxv3cBpzUlvHCMUuTZ14tt/vi9blMpon64m/6OUCGr23O14R57doxwnOmaf27H4uch/djsE521yt+/ly8vfe1jZh1qpHf3O73749L3kBfd1iSLr7Y5BwsMcnZvVfP/dmc3WXae1dCgVIqR2u9wHzg2Ov4HPnmMgnHH7Pr4xfg1m3ppE1XCVz447VRs4gFms3n/fx8bqcsPBeUXcDQeuQp8fUAXxm1ozssXymlcXwhyDc/37VebveplvcjG7fX7/ygNh/SHre5cH/vsnEcn3T9QLfhf9Fw8vfnYq/jft3VeI46/v44Bet1hwVpQcWmBQR3eKqzj3u2j3NIZvDtt7g6F0UPLSVXS4Ecs839mEa6Umq2y7ddexCz+Rxi7+dzO9l8bEutRx5fz+kzo1LKZj4kU3B0E+Wa+3q8vbZwdXw/6sp9XzYcx9sKXC4LtNYZNR96oQB4YvPxnK4/l/qebmF3v6Gu75dLEfP7dYc7aUHFIK11kVKqyNM30wAVmf3OVEqt8dQq01qvNX9w9fkmVwA8RM0TXOfjaBVme/hmvQborr89eHxhg/mgtQcpmyc+n9vtvgU4WoHuMnAU4FDxmhF4DLjPbCvA0dJcivkZeLi9Ps+1Ew/f8s37tNbbNh9dsWvx/H56k4nnvwWrfi5Qt98fcLyGur7usCYtqNgxzO36POAhHfxzPu7G+4mOs80l0Od0Hvx17+6xe9qns7Xl9iFmM9tsXPwNOJBsnlo2zv3X5bmdxwaLXFuI5n5Z/nSf1ZbHEz8yNvFw/Ms50tHj7fV4rpPAfA/7zTY/b4/b3PflZL4guR9rdc3tfozINZfrfvz9udTpva/tvn7+/tR4DX687shi9SgNuYT+gqMbRuP4ALaZ22y4jBLCcdA5H8douSl822/tvG2i2/UpmBF8zu1mP5nmudZQc0SdDT9HpXl5HTa8jBbDw+hBl9c+BceHWabZx2zcRjXVJZvb63aOuMrBcVziwuv29dzu76XLvqeY7TnOn0MQ80xx+T1I9yOj83bnxZnb4+1+/h76/FmY7ROd+3V7fI1t7q/d/f1yeUyOj/0593XR77uH+9b4ufj73vv7s3PNX8f3y+trc98WSRdlXowQDcI5GMHqHJ6EczYhYpEUKNFggjRiMCTCOZsQsUqOQYmQMqOQvJ2jZKlwziaEkBaUCDGX8zFs4dZ9Fs7ZhBBSoIQQQoQp6eITQggRluRE3SBq3bq17tatm9UxhBAiYqxZs+aY1rqNp21SoIKoW7durF692uoYQggRMZRSe71tky4+IYQQYUkKlBBCiLAkBUoIIURYkgIlhBAiLEmBEkIIEZakQAkhhAhLUqCEEEKEJTkPSogwc+LMefYdP8vB0jIOlpZzqrySsooqys5XEacUSYlxJDWKx9Y0kY62JDrYmtAltSlNE+XPWUQX+Y0WwkLnK6tZv9/Oyl3H2PhNKVsOnuRgaXmN+yXGx9G4URxoOFtRRVX1xXNoKgXdWzejb4eWXNLZxuUZrejbviVxcarGvoSIFFKghGhgJ8srWLr5MIs3HuTzohLOnq8iTkFGm+YM755Kv7SWdG/dnA7JTeiQ3ITkpEYkxF/cG3++spqSM+cotpdTbC+j6OgZthws5ev9dt7bcBCA1GaJXNWzNd8d2IFrerWhSaN4K16uEAGTAiVEA6iu1nxSeJTXvtzHx9uPcr6qmo62JHKGdmJEj9Zclt6K5KRGfu8vMSGODslJdEhOYmjXlIu2HSwtY8XOElbuPMaH24/w1vpimjdOYHT/dtwxvAuZXVJQSlpWIvzJchtBlJWVpWUuPuHqVHkF81bt55Uv9rK35CytmzfmpsFp3Di4A0M620JeKCqqqvl8Vwnvbihm8cZDnD5XSZ/2Lbjz8q7cltlJWlXCckqpNVrrLI/bpEAFjxQo4XSyvIJ/r9jDi5/tprSsgmHdUrjz8m6M7d+exARrBs+eOVfJW+uL+c8Xe9ly8CRtWjRm0tXp/GB4FxlgISwjBaqBSIES5RVVvLxiD89/vJOT5ZVk923LL0f2ZHBnm9XRLtBa80XRcf7+USErdpaQ2iyRe67twZ2XdbWseIrYJQWqgUiBil1aa97ffIgnFm9l//EyRvZpy6+zezGwU7LV0Xxas/cEfy7YwfLCY3Rr1ZRp1/dhTP/2coxKNBgpUA1EClRs2ldylofe2MCKnSX0atec33+vP1f2bG11rDr5ePsR/rR4KzsOn+ay9FT+dOtA0ts0tzqWiAFSoBqIFKjYUlWt+ffKPcx6fzvxcYrfjOnNHcO71BgSHikqq6p5fdV+8v67jXOV1dx3XU8mXp1Oowh9PSIySIFqIFKgYsfekjPcP/9r1uw9wbW92/DErQNJsyVZHSsojpwsZ/o7m1m88RB92rfgqdzBDOgY3l2VInL5KlDy1UiIOnpr/QFu+OtnFB4+xTPjB/PPnwyLmuIE0LZlE567Yyhzf5TF8TPnufW5Fcz9tIjqavkyKxqWjC0Vwk9nz1cy/e3NzF/9DVldU/jL9y+hYxQVJnej+rUjq2sK0xZt4InFW/l4xxGezh1C++QmVkcTMUJaUEL4YV/JWW79x0ry13zDvdf24PWJl0V1cXJKaZbICz8cypPjBrJ2r50b/rqclbuOWR1LxAgpUELUYsXOY9z0j884dLKcf//0Uh4c0ztiB0IEQinF7Zd24Z1fXomtaSPufOkr5n5ahBy/FqEWO39lQtSR1pqXV+zmR//8ijbNG/P2vSO4ulcbq2NZpkfb5rx175WM7teOJxZv5d7X1nHmXKXVsUQUkwIlhAdV1ZpH397MY+9s4drebVn0iyvo2qqZ1bEs17xxAs/dkcm06/uwZONBcl/4nEMelgcRIhikQAnhpryiintfXcv/+3wvd1/VnTl3DqVFE/9nGo92SikmX5PByz+9lH3Hz3LLP1awpfik1bFEFJICJYSL0rIKfvTPr1iy6RCP3NCX397QTxb98+KaXm3In3w5SkHuCyv5aPsRqyOJKCMFSgjjyKlyxr/wOev2neAvtw/hrqvSrY4U9vp2aMkbvxhB11bNuOvfq5m/ar/VkUQUkQIlBHD4ZDm3z/mC/SfO8vJPLuXmIR2tjhQx2ic3IX/y5VyR0YopCzfw0me7rY4kooQUKBHziu1lTJj9OYdLy/n3zy6NuIlew0Gzxgm8+OMsrh/Qnsff3cKzS3fIMHRRb1KgREzbf/wsE+Z8Tsnp87xy13CGdUu1OlLEapwQz9++fwm5Qzvxl2WF/OHdLTI9kqgXmepIxKxiexm3z/mCU+UV/Oeu4WG1qGCkSoiPI++2QTRvksDLK/ZwrrKaP948QAaaiIBIgRIx6djpc/zwxS85WVbBq3dfFvYLC0aSuDjF77/XjyaN4nn+413EKXj85gGyCKKoMylQIuaUllXwo5e+ori0jFf+Z7gUpxBQSjFlTG+qtWb2J0XEKcVjN/WXIiXqRAqUiClnz1fys3+tovDIKV788TA55hRCSimmje1DdbVm7vLdxCnFozf2kyIl/CYFSsSM85XVTHplDev2neAfP8jkmhieV6+hKKV4+Lt9qdbw0me7adIonmnX97E6logQUqBETNBaM23hBpYXHmPmbYO4fmAHqyPFDKUUj9zQl/KKKl74ZBcpTRsx6ZoMq2OJCCAFSsSEZwsKWbTuAPeP6sX4YZ2tjhNzlFL84eYBlJZVMGPJNlKaJsrPQdRKCpSIevNX7+evywoZn9WJX47sYXWcmBUfp3hm/BBKyyqYtmgDLZMaMXZAe6tjiTAmJ+qKqLa88CgPL9rIVT1b88StA+UAvcUSE+KYfedQBne28b+vrWPlTlmdV3invE1HopQaB1wK1PdU8FVa60X13IdllFITzX+HAnla6yJv983KytKrV69umGCiVtsPneK251fSKSWJ/MmXy5IZYcR+9jwTZn/BAXsZC35+OX3at7Q6krCIUmqN1jrL4zYfBep5YFoQnn+O1npCEPbT4JRSmQBa67VKqWxgqtZ6lLf7S4EKHyfOnOemf3zGuYpq3rp3BB2Sk6yOJNwcLC3jln+sIE4p3rxnBO1aNrE6krCArwLlq4uvSGtdWt8LEMmf2OnAJPP/1YDHN1GEl8qqau59bS2HS88x+86hUpzCVIfkJP75k2GcLKvgZ/9aJcvHixq8Fiit9SxPtyul5iulHvT3CbztpyEppWxKqSlKqTwv26copXKUUhNduvTQWi/QWjsLVBaRXWxjxhOLt7JiZwl/GjeQS7qkWB1H+NA/LZm/35HJtkOn+OVr66isqrY6kggjgQySWArM9bRBKRV2Hcmmay4byABsHrbn4WgtLtBazwEylFI5HnY1CcgNZVZRf/NX7+flFXv42Yju5AztZHUc4Ydre7flDzf358NtR3jsnS2yTIe4IJACtQvw9rV0opfbLaO1LtBaLwDsXu4y0Wx3mse33XrAhYESU7XW3vYhwsDafSd45I1NjOjRioe/K7MVRJI7hndl0jXpvPLFXv65Yo/VcUSYCOQ8qPHAUKWUDSji2w/+FBwj3Z4KSrIG4BwE4caOo8XlvE82UKC1LlJKZWutCxoqn/DfsdPn+Pl/1tA+uQl//34mCfFyBkWkmTqmD3uPneWJ97bQq11zruopU1HFukAKVBYwFTjudrsiOKP+GlIqNV/HheumgOUDx835M2sBKVBhpqpac9/r67CfreCNX1xKSrNEqyOJAMTFKZ4eP5jbnj/Dva+u4817RtC9dTOrYwkLBfI1826t9TKt9Tq3y1pgRrADhpjN2wallE1rvVZrnaK1zjAXOQYVhv6yrJAVO0v4w8396ZcWdodBRR00a5zA3B9lEafg7v+3mlPlFVZHEhYKpEB5PYKptV5XjyxWsONoRbkKeP2F4uJilFIXLtOnT69PNuGHT3cc5W8fFnJbZifGZ8ncbtGgc2pTnrtjKHuOneFXr6+nSpaNj1mBFKi54ThaL0DHqdmKsgEEMiAiLS0NrfWFixSo0DpYWsav5q2nV9sW/PEWWbE1mlye0YpHb+zHsm1HePqD7VbHERYJ5BjUCSBbOT4NTmitPwxypgZjZoiwu92cihxnCnsVVdXc++o6zlVU8dwPM0lKjLc6kgiyH17WlS0HT/Hcx7vo26ElNw5OszqSaGB1LlBa69HO/yulkpVSt+Ho9lurtd4TxGwNZb5SKsdlqPkoYLaVgUTtnnp/O2v2nuBv37+EjDbNrY4jQkCZZeJ3HD7F1IUb6NuhBT3atrA6lmhA9RqLa6YzWgisA+Yopd43k8yGDaVUplJqCpCDo+U3xXV4uZkpIl0plW3Od9rldl6UCDOfFR5j9qdF/GB4F/lWHeUSE+L4xw8ySWoUz6RX1nBapkOKKV4ni/X6AKW6OVtKphhNBroDc3BMDFuqlLoOSI7kWcwDIZPFht7xM+cZ++dPaZnUiHfuvVK69mLEyp3H+OFLX/LdgR342/cvkeONUSTQyWK9yVdKPa+UOg5MwDHDQk+t9SwzOSxmGPqicGtNicimtWbKgg3Yz1bwl9uHSHGKIVf0aM2DY3rz7oaD/GvlHqvjiAYSSIFKAdZorVO11hO8DS1XSg3B+/RCQtTZ/325j4Kth5kytjf905KtjiMa2OSrM8ju244n3tvKmr3u59eLaBRIgcrTWr/ox/0mB7BvITwqPHyKx9/dwtW92vCzEd2tjiMs4JxpomNKEr/4v7UcO33O6kgixOpcoLTWHmcy93C/yZE8BF2Ej/KKKv739fU0b5zAU7mDiIuT4w+xKjmpEc/fMRT72Qr+97V1chJvlPNaoJRS3ZVSLyilRjZkICHcPf3BdrYePMms3EG0bSGrrsa6fmktefyWAazcVcI/PtppdRwRQr4WLNyNY1LYoUqpD5RSM5RS3RosmRDAqj3HefGz3dwxvAsj+7SzOo4IE7lDO3HLkDT+XLCDr3bL8aho5bOLz5znNMucnDsHmGzOdXowiqY7EmHqzLlKHpj/NZ1Sknj4u32tjiPCiFKKP946kC6pTbnv9XWcOHPe6kgiBPw+BqW13q21nqa1HoPjxNyZSql5MpRchMqTS7ax/8RZnsoZTLPGgczKJaJZ88YJ/P0HmRw7fY7fLPhaVuKNQgHNJGHOc5qstZ4AKKXUfFOs5HiVCIrPCo/xyhd7+ekV3Rme3srqOCJMDeiYzEPX96Vg6xE5PyoK1XvZUa31Qq31eBzLvWeY41XPm/OghKizk+UVTFnwNeltmjFlbG+r44gw99MR3cju25YZi7ex6UCp1XFEEAVtXWxzvGquOV41Exhljld1C9ZziNjwx3e3cOhkOU/nDqZJI5ktQvimlGJWzmBSmyVy76trZb6+KBK0AuXKHK+apbUeE6EznAuLLNt6mPmrv2HyNRlc0iXF6jgiQqQ0S+Qvtw9h3/GzPPLGRqvjiCAJSYGKVc4VdWWhwsCcLK/g4Tc20rtdC+7L7ml1HBFhhqe34r7revHm+mLeWn/A6jgiCGRoVBClpaVRXFxsdYyI9eSSbRw9dY45d2bROEG69kTd3XNtBp8WHuWRNzeR1S2VjrYkqyOJeghKC8rMOnGbUmqk+VeGnos6+aKohFe/3MfPRnRncGeb1XFEhEqIj+PZ8UOortbcP2+9TIUU4epdoMxovWQzmu9D8+8iU6jkZF5Rq/KKKh5atJHOqUncP7qX1XFEhOvSqimP3tSfL3cf58XlRVbHEfUQjBZUhtZ6vfuNZqXd7CDsX0S5vywrZPexM8y4dRBNE6XXWdRf7tBOjO3fnqc+2M7mYhl6HqnqVaDMEPK1LtcfdDtZV34zhE+bi0uZ82kRuUM7cWXP1lbHEVFCKcWfxg0kpWkiv563nvKKKqsjiQDUtwWVApS4XL8dGOVyXTqAhVeVVdVMXbiBlKaJPHJDP6vjiCiT2iyRWbmD2XH4NHn/3WZ1HBGAehUos5pulsv1LK31Qy53kRNZhFcvfbabTQdO8oeb+5PctJHVcUQUuqZXG35yRTdeXrGHT3cctTqOqKNgHIPyWISUUt0BOUIpPNpbcoZnlu5gVL92XD+gvdVxRBSbdn0ferZtzoP5X8us5xHGZ4HyZxSe1nqhc4i5eUyyUuo2HCP71gUpp4giWmt+99ZmGsXH8fjNA1BKVsgVodOkUTzPThjCibPnefiNjTLreQSprQX1UC3bgQsTxn6olLoO6G6ur693OhGVFm88xKc7jvLA6F60T5YVckXoDeiYzAOje7Nk0yHeWCezTESK2gpUjlLqWn93ZpbhWF+/SCKanSqv4A/vbqZ/WkvuvKyr1XFEDLn7qnSyuqbw6NubOVRabnUc4YfaClRroEApVaWUWmWWfR/pretPTswVtXlm6Q6OnDrHE7cOJCFepoIUDSc+TvFU7mAqqzRTF26Qrr4IUNsnRCYwGRgPzAeGAgXACaVUoVn3aZzLkhp+dQmK2LTpQCn/XrmHO4Z3YYhMZyQs0K11M6Zd34dPdhzl9VX7rY4jauGzQJllM+YCdqBIaz1aax0HDMOx5lMr4EVgl1KqBMeihULUUF2teeTNTaQ2S+Q3Y/pYHUfEsDsv68rl6a3447tb2H/8rNVxhA9+9bGYY0vO0XpDtNZrzeKE47XWqUAPHK2n4yFNKyLWa6v2sX6/nd/e0JfkJDnnSVgnLk4xM2cQAFMWbKBaJpQNW3U6CGDm1ys13XotXW7frbWeAywMdkAR+Y6eOkfekm1cnt6KW4Z0tDqOEHRObcoj3+vH50UlvPLFXqvjCC/qfJTaFKNFwDAPy2rMDk4sEU1mLN5KWUUVj98i5zyJ8HH7sM5c06sNTy7Zxp5jZ6yOIzwIeBiV6fZb5HqSrtZ6d/CiiWjwZVEJi9YdYNLVGfRo29zqOEJcoJTiydsGkhCveDD/a1k7KgzVdzbzlsAuIF0pNc9lNJ8QVFZV8+jbm+loS+Kea3tYHUeIGjokJzH9xv6s3nuCf34m36/Djc/Fd0wBSne5ZLhdB1CYUX7AJGSouTBe/Wof2w6d4vk7MklKlCXcRXgal9mRJZsOMeuD7Vzbpw092rawOpIwamtB2YE1wFwcS2m0AtbhGGI+GuihtY7TWqd6mMk85hQXF6OUYvr06VZHsdzxM+d5+oMdjOjRirEyGawIY461owbQLDGeB/M3SFdfGKlt+dIiHAMf7ECBHGPyLS0tjeLiYqtjhIVZ72/n9LlKHr2xvwyMEGGvbYsmTL+pP/e9vp6XV+zmrqvSa3+QCLnaWlALtNazzMm6NqXUXS6Xke539jCqT8SgTQdKeX3VPn58eTd6tZPuEhEZbhqcRnbftsx6f7uM6gsTKtD5qMx6T9mADcfKuaXAFK11z6ClizBZWVl69erVVsewlNaanBc+Z2/JGT588Du0bCIn5YrIcai0nFHPfkK/Di157e7LiIuT1n+oKaXWaK2zPG2rzzDz3WY2iVla66dwHKtqFej+RHR4Y90B1uw9wZQxfaQ4iYjTPrkJv7uhH1/uPs7/fbXP6jgxL2jTSWut1wJzgrU/EXlOn6tkxpJtDO5sI2doJ6vjCBGQ3KxOXNWzNU8u3so3J2SuPit5LVCBnNOktZ4WjP2IyPS3ZYUcPXWOx27qL10jImIppfjTrQPRwEOLZAVeK/lqQU0N0nPkBWk/IoztOnqaf67YTe7QTrKUhoh4nVObMu36PiwvPEb+mm+sjhOzfA0zV0qpGfXcvzIXEeX+8M4WmiTEM2WsLKUhosMPh3fl3a8P8vi7W7imVxvatWxidaSY46tABasFJaLcR9uP8MmOozxyQ1/atGhsdRwhgiIuTpGXM4ixf/6U376xibk/Girn9DUwrwVKa13akEFEZKqsquaJ97bSrVVTfnR5N6vjCBFU3Vs348HRvXli8Vbe2XCQmwanWR0ppgRtFJ+ITa99tY+dR07z0Hf7kpggv04i+vzsyu4M7mxj+tubKTl9zuo4MUU+UUTASssqeGbpDi5LT2V0v3ZWxxEiJOLjFLNyBnG6vJJH395sdZyYIgVKBOzvHxZiL6vgkRv6Sd+8iGq92rXglyN78O6Gg/x30yGr48QMKVAiIHuOneFfK/eQk9mJAR2TrY4jRMhN/k4G/Tq05PdvbaK0rMLqODFBCpQIyJNLttEoPo7fjOltdRQhGkSj+Dhm5gyi5Mx5nlyyzeo4MUEKlKizL4pK+O/mQ/z8mgzayrkhIoYM6JjMXVd257Wv9vFFUYnVcaJenQuUpyU1PC29IaJTdbXmj+9toUNyE1kzR8SkX2X3oktqUx5atJHyiiqr40Q1vwqUUmqcUsq5sM8wD3dZo5R60CwRL6LYwrXfsOnASaaO7SPLuIuYlJQYz4xxA9l97Ax/+7DQ6jhRrdYCZdZ9ehiwK6UKgWyl1K2uxUhrXWqW3MgOXVRhtbPnK5n1/nYGd7bJCYsipo3o0ZrcoZ2Y/UkRW4pPWh0natVaoMy6T1la63jg50AK8FtMwVJKPW9W2O0GSJ9PFHvhkyKOnDrH77/XV2YrFzHvtzf0xda0EdMWbaCyqtrqOFGpTsegtNYFOJaBz9JaxwGTcaykOxlYChQFP6IIBwdLy5jz6S6+N6gDQ7umWh1HCMvZmiYy/ab+bPimlH+t3GN1nKjka7JYb5Y6/6O1XgYsC14cEa6e/mAH1dUwVWYrF+KCGwZ24M2+B3jqg+2M7teeLq2aWh0pqvhsQSmlSpRSq5RSM5wj9UxRcr/fdZ5G94nosO3QSRau/YYfX9GVzqnyByiEk1KKx28ZQEJcHL99UxY3DLbauvjW4FjGPQWYo5Sqci1YzoESzpaUUuquEOcNa8XFxSilmD59utVRgurJJdto0TiBe67tYXUUIcJOh+Qkpo7tzfLCYyxae8DqOFGlti6+PFN85gIopZKB8cAkYBSQqZTaBRTgKGZDgRdDFze8paWlUVxcbHWMoFqx8xgfbz/Kw9/tg61potVxhAhLdwzvypvri3n8vS1c07sNrZvLumjB4LMF5d6dZ9aI0maQhHOgxM9xrJo7FcgPWVLR4KqrNTOWbKWjLUnWehLCh7g4xZPjBnL2XBV/eGeL1XGiRr2nOtJaF2itJwOjg5BHhJF3NhSz6cBJHhjdiyaN5KRcIXzp2a4F91zbg7e/LubDbYetjhMVAilQKZ4GRGitdwO2eicSYeFcZRWz3t9O3w4tuWVIR6vjCBERfv6dDHq1a84jb2zi9LlKq+NEvDoXKK31LGCyUup99xkl8DwNkohAr3y+l29OlPHwd/vISblC+CkxIY4Z4wZx8GQ5T72/3eo4ES+gLj6t9Wgc5z+9BJwwo/tKgFXBDCesUXq2gr99uJOrerbmqp5trI4jREQZ2jWFH1/ejX9/voc1e09YHSeiBXwMSms9U2udCvQAsoB0rfWioCUTlnnuk52cLK9g2vVyUq4QgXhwTG86tGzCtIUbOF8p0yAFKhiDJHZrrdeZEX4iwh2wl/Hyij3cOqQj/dNkpVwhAtG8cQJ/vHUAhUdO8/zHu6yOE7G8FiilVHel1Auy1lNsefoDR7/5/aN7WZxEiMg2sk87bhqcxt8/KqTw8Cmr40QkrwXKjMqbCgxVSn1gZo/o1mDJRIPbXFzKG+sO8NMrutEpRaY0EqK+fn9jP5o1TmDaoo1UV8s0SHVV24m6pVrrWWZQxBy+Hb0nixNGoSeXbKNlk0b84jsypZEQwdC6eWN+d0M/1uw9wf99udfqOBHH72NQ5ljTNK31GGAdMFMpNU8miY0OywuPsrzwGL8c2YPkpo2sjiNE1BiX2ZGrerYm77/bKbaXWR0nogQ6zHyZ1nqy1noCoJRS802xkuNVEai6WvPkkm10tCVx5+VdrY4jRFRRSvGnWwdSVa353ZubZMbzOgjGKL6FWuvxwEQgwxyvel4pNaTe6USDeG/jQTYXO6Y0apwgUxoJEWydU5vywOheLNt2hPc2HrQ6TsSod4FyMser5prjVTOBUeZ4VbdgPYcIvoqqap5ZuoNe7Zpzs0xpJETI/OSKbgzqlMz0tzdjP3ve6jgRIWgFypU5XjVLaz1Ga70nFM8hgiN/9TfsPnaG34zpQ7xMaSREyCTEx/HkuEGcOFvBE+9ttTpORAhJgRKRobyiir8s20FmFxvZfdtaHUeIqNcvrSWTrk4nf803fFZ4zOo4YS9oBcos+/68ucjIvgjw75V7OHzyHFPH9kEpaT0J0RD+97qedG/djIff2EjZ+Sqr44S1oBQopdRtQCawwFwyYn3593BXWlbBcx/v4ppebRie3srqOELEjCaN4pkxbiD7jp/l2YIdVscJa8FqQdnNMadl5jIL2B2kfVtOKWVTSuUppTKtzhIscz7dRWlZBb8Z09vqKELEnMvSW/H9S7vw4vIiNn4j05h6E6wC5WlgfzTNNJoFpFsdIliOnCrnn5/t4cbBaQzoGE0/JiEix7Tr+9C6eWOmLNxARZXMeO5JsApUipmrb5y5zABSg7Rvy2mtC4DjVucIlr9/uJPzVdXcP0omhBXCKslJjXj8lgFsPXiSucuLrI4TlhKCsROt9UKlVBEwwdw0T2u9rrbHKaVswEOAcz761VrrtcHI5OW5JgKttNZTPWyfAhRhCqvWek4oclhtX8lZXv1yHxOGdaZ762ZWxxEipo3p357rB7TnzwWFXD+gg/xNugm4BaWU6uZ6AU4AL5jLCdOK8vV4G5CvtZ7qUgweCjRPLc+VDWQDGYDNw/Y8oEhrvcBkyVBK5YQii9WeWbqd+DjFfdf1tDqKEAJ47Ob+NEmIY9rCDTLjuZv6tKAKgDWA6/hk57urgEvwXXDmArNdrs83+ww600WHUmoYHgoUMNGtVTUPyMMxIjFqbD14kre+LmbS1Rm0a9nE6jhCCKBtiyb89oa+TF24kXmr9/P9S7tYHSls1KdATdJaL/O2USl1SS2PzwHuVkqlAzbTtWf3sb+JOFo5BW6324C5Wutcf4O7Pd7TyDw7jhZXVHnq/e00b5zAz6/JsDqKEMLF+KzOvLmumD8t3srIPm3lC6QRcBdfLcVpCNDdx3ZnUchyuS3fFBtvzzcHmORaUJzdhMDd/ub2IJWaAyAuum66CLOACb6GmhcXF6OUYvr06fWIExqr9xxn2bYjTL4mQ5bTECLMKKWYMW4g5yur+f1bm6yOEzaCMkhCKdUdR3ddCt92+S0FFnl5yIUh21rrIrOPeTi6/by2hLTWuaaQzTAtrrk4WnL2esS3eduglLJpre2m1Ta0th2lpaVRXFxcjyihobUm77/baNOiMT8d0c3qOEIID7q1bsavR/XiySXb+O+mg4wd0MHqSJYLSoECJpkLfNv68HXekN38u9rltiIc3X4+uRQpG47iVN/xmXZqDomPmiHyAB9vP8qqPSd4/Ob+NE0M1o9cCBFsd13ZnXe+LuZ3b23m8vTWMd/bEazzoOaZGcx3A9dprUtxtKa8KQJwa/nY4UK3XUM6Ts1WlA1q5ItI1dWame9vp0tqUyYMk4OvQoSzhPg48m4bxPEz55mxRGY8D1aBSlVKrTL/v9QMMffVVVcE2N2KkQ3HlEl2X0+klMoHZmitRwGzzSCLgHkZnJFKiEYUNrR3NhSz9aBjMcLEBJm8XohwN6BjMndflc7rq/azcldsz3gelE8sM//eMPP/aTi67mbW8rAZwHiX6xPMbV4ppWYDs11O5s3FUaRsgeR2Md/tvKdRXDwEPiI5FyPs074FNw5KszqOEMJPv8ruSddWTXlo0UbKK2J3xvNQLVi4EM/z87neZyZgU0pNMbM4lJjbPDLDzPNdh5mb1lYujsESXimlMs1z5ADZ5jkvjMbTWk8C0pVS2eZ5dmmtI/4cqPzV37C35Cy/GdObOFmMUIiI4ZzxfG9JbM94rrSu/5nLZiaJXL4dXKBwHIsaVu+dR5CsrCy9evXq2u/YAMorqrj2qY/pkNyEhT+/QtZ7EiICTVu4gfw13/DWPSOidmJnpdQarXWWp23BakFNwzHwocBcluKYZUJY5NUv93GwtJwHR/eW4iREhHro+r6kNktk6sINVMbgjOfBKlCztdYLXdaDWgbUmJBVNIyz5yt57uOdXJ7eiit6tLY6jhAiQMlNG/H4zf3ZXHySFz+LmiX2/FafyWJbOi84ltsY6TZ57LSgpRR18q+Vezh2+jwPjpHlNISIdGMHdGBM/3Y8u3QHe46dsTpOg/J51qZSqqXW+qSXzWupOVmsq9omixUhcLK8gtmfFHFt7zYM7RpV5xsLEbP+cPMAsp/5hGmLNvDqXZfFzKCn2qYVeAjvRaa+k8WKEHhp+W5Kyyq4f5Qs5S5EtGjXsgm//W5fpi3ayGur9nHH8K5WR2oQtXXx5SilrvW0wVdxMttrXbBQBNeJM+d56bPdjO3fnoGdonPEjxCxasKwzozo0YoZi7dRbC+zOk6DqK1AtQYKlFJVSqlVZln3kea4Uw3ebhcN44VPd3HmfCX3j5ZjT0JEG6UUT44bRFW15rdvbCQYpwiFu9oKVCYwGceMD/NxzOhdgGPF3EKl1PNKqXFmUATIMSfLHDlVzr9X7uHmwWn0atfC6jhCiBDonNqUKWN789H2o7yx7oDVcULOZ4EyE8DOxTFXXZHWerTWOg4YhmMqo1bAi8AupVQJMDHEeYUXz320i4oqzX3Z0noSIpr9+PJuDO2awmPvbOHIqXKr44SUX8PMzblNC5VStymlhmit12qt52qtx2utU4EeOFpP7gv/iQZQbC/j1S/3kZPZie6tm1kdRwgRQnFxirzbBlFWUcWjb222Ok5I1ek8KDPHXqnp1mvpcvtus+LtwmAHFLX724eFaDS/vK6H1VGEEA2gR9vm/Cq7J0s2HWLxxoNWxwmZOp+oa4rRImCYUmqc2+aInwE80uw5dob5q7/hB5d2oVNKU6vjCCEayMSr0hnYMZnfv7WJE2fOWx0nJAKeScJ0+y0y3X4jzW2xNxeHxf66rJCEOMU910rrSYhY4lzc0H62gsff3WJ1nJCo11x8pptvF46lKua5jOYTDaDw8CneWH+AH1/RjbYtm1gdRwjRwPqlteQX38lg0boDfLTtiNVxgs5ngTJz7Q0xx5weNMPK3zdDzKuAEzimPJoJZACTGiCzMJ4t2EHTRvFMvibD6ihCCIvcM7IHvdo15+E3NnKqvMLqOEFVWwvKjmO+vbnA7TiGla/DUZBGAz201nFa61StdZbWWs6DaiCbDpSyeOMh/ufK7qQ2S7Q6jhDCIo0T4pmZM5jDJ8uZsWSb1XGCqra5+IpwDHywAwVyjCl8PLt0By2bJPA/V6VbHUUIYbEhnW38z5Xdmbt8N98b1IErMqJjmZ3aWlALtNazzMm6NqXUXS6Xke539jCqT4TA2n0nWLbtCJOuySA5qZHVcYQQYeD+Ub3p1qop0xZu5Oz5SqvjBIXPFpTWeprL/9fh6N4DQCnVXSl1N2ADNFAKTAEWhSSpuODpD7bTqlkiP7mim9VRhBBhIinR0dU3Yc7n5C3ZxmM3D7A6Ur3V1sXnlenum+u8rpTKxHGMSoTQyl3HWLGzhEdu6EuzxgH/+IQQUejS7qn85IpuvLxiD2MHdODyjMj+SA7Wku9ordcCc4K1P1GT1ppnPthBu5aN+eFlsbEejBCibqaM6UO3Vk35zYKvOXMusrv6glag4OIuQRF8H+84yuq9J7h3ZE+aNIq3Oo4QIgwlJcYzK3cwB+xlPBnho/pqOw+qxHUdKB/3u04GSISW1pqnP9hOp5QkJmR1tjqOECKMDeuWys9GdOeVL/aycucxq+MErLYW1Boc3XYpwBxvCxea1XWXKaXuCnHemPX+5sNsOnCS+67rSWJCUBu+Qogo9ODo3nRv3YzfLNjA6Qjt6qvtky7PLKsxWWvdA0jFUbBG4ThZ1+5cuBDIxbGgYcwqLi5GKcX06dODut+qas0zS7eT3roZt17SMaj7FkJEp6TEeJ7KHURxaRkzFm+1Ok5AahtmvszteqlSSmuts5y3KaWygRxgKjE+1VFaWhrFxcVB3++7G4rZcfg0f/3+JSTES+tJCOGfoV1TucucwHv9gA5c2TOyTuCt96ed1rpAaz0Zx9RHIsgqq6p5dukO+rRvwfcGdrA6jhAiwjwwujfpbZoxdeGGiJurL5ACleJpQIQ5L8pW70TiIovWHmBPyVnuH9WLuDhldRwhRIRp0iiep3IHc7C0jD8tjqxRfYEsWDgLmGxmNb/VdWVdYFjwoolzlVX8ZVkhgzslM6pfO6vjCCEiVGaXFO6+Kp3XvtrHpzuOWh3HbwF18WmtRwPLgJeAE2Z0XwmwKpjhYt28Vfs5YC/jgdG9UUpaT0KIwP16VC8y2jRj2sINnIyQrr76rKg7U2udCvQAsoB0sxS8CIKy81X87cOdXNotlasi7MCmECL8OLv6Dp0s50/vRcaovmAMktittV6ntS4NRiDh8J8v9nL01DkeGN1LWk9CiKC4pEsKE6/O4PVV+/l4e/ivwCtjlsPQ6XOVPP/JLq7q2Zrh6ZE92aMQIrz8KrsnPds2Z9rCjZSWhXdXnxSoMPTyZ7s5fuY8D4zubXUUIUSUcXb1HT19jsfe3mx1HJ+kQIWZ0rMVzFleRHbfdgzpbLM6jhAiCg3ubOPea3uwaN0Blmw8aHUcr6RAhZm5y4s4VV7J/aN6WR1FCBHF7h3Zg0Gdknn4jY0cOVVudRyPpECFkWOnz/HPFbu5YVAH+qW1rP0BQggRoEbxcTwzfjBnz1cxbeFGtNZWR6pBClQYeeHjXZRXVPHrbGk9CSFCr0fbFkwd24cPtx1h3qr9VsepQQpUmDhUWs4rX+zl1ks60aNtc6vjCCFixE+u6MaIHq14/N0t7Cs5a3Wci0iBChP/+GgnVdWa+67raXUUIUQMiYtTzMoZTFyc4v7566mqDp+uPilQYWD/8bO8vmof44d1pkurplbHEULEmDRbEo/d1J/Ve08wd3mR1XEukAIVBv66rBClFL8c2cPqKEKIGHXrJR25fkB7nvlgB1sPnrQ6DiAFynJnzlXy8Y6j/HB4VzokJ1kdRwgRo5RSPHHrQFomNeLX89ZzrrLK6khSoKzWrHECHz34He7LlmNPQghrpTZLJO+2gWw7dIpnlxZaHUcKVDho3jiB5KRGVscQQgiu69uO24d1Zvanu1i157ilWaRACSGEuMgj3+tHp5Qk7p+/ntPnKi3LIQVKCCHERZo3TuDp3CF8c6KMP767xbIcUqCEEELUcGn3VCZf41g76r+bDlmSQQqUEEIIj36d3YuBHZOZtmgDh0obfkJZKVBCCCE8SkyI48+3D+FcRTUP5K+nuoFnmZACJYQQwquMNs159MZ+rNhZwoufNewsE1KghBBC+DRhWGfG9G/HrPe3s+lAaYM9rxQoIYQQPimleHLcIFKbJXLf6+soO98ws0xIgRJCCFGrlGaJPDN+CLuOnuGJxQ0z9FwKlBBCCL+M6NGaiVen858v9rF0y+GQP58UKCGEEH57YHQv+nVoydSFGzhyMrRDz6VACSGE8FvjhHj++v0hnD1fyQP5X4d06LkUKD8opWxKqTylVKbVWYQQwmo92rbgkRv6sbzwGC+v3BOy55EC5Z8sIN3qEEIIES7uGN6F7L7tyFuyLWQLHEqB8oPWugCwdt55IYQII0op8m4bSHLTRtz3+jrKK4I/9Dwh6HsMkFJqttZ6Ugj3bwMmAq201lM9bJ8CFAGpAFrrOaHKIoQQ0aBV88Y8nTuY5YVHUSr4+w+LAqWUyiOEXWhKqWzABmT4eP5VWusFzutKqRzndSGEEJ5d3asNV/dqE5J9W97F1xADD7TWBabY2L3cZaJbMZoHhKw1J4QQonaWFygcAxCW1nYnpdRE0xJyv92mlMoP9Mm9FEg7UOO5hBBCNBxLC5RSKgeY7899zTGhSa4FxRxXygfurkeMVGoOgLjouimMWcAEXy2+4uJilFJMnz69HnGEEEKAhcegTHGxa63tys+ja1rrXKVUvlJqhtZ6LTAXmKS1ttcjis1XRq213YziG1rbjtLS0iguLq5HFCGEEE5WtqDGmw/+OtFa5wIPKaWWAlO11vVdoMSOGbnnwv26CCFpcQoRPNH092RJgTLdZHUuTiFynJqtKBtAPVtmwk+PPfaY1RGEiBrR9PdkVQsqFchRSk0x5x9NAtLNdZ/Dzc2AiBla61HA7NruXxvTVWj3kK9BC2gwvvUEuo+6Pq4u94+mb3OhEo7vUUNmCvZzyd9S9FBaN+wa8x5DKDURx7Ekn8d5lFKzgXxn16DLIIlcf1o75nwnm/sJwWa/S13Pg8LlvCh/ZWVl6dWrV9flIa4ZqO/PItB91PVxdbm/P/cNxmuPZOH4+hsyU7CfK5b/lgLJYDWl1BqtdZbHbVa/EFOccnGMkpsBzPFUbMz9ityPW5kiNdccm/L2HJk4ho07C9NsoMC0npz3mQKsxZwwHMhMEkqpo8Deuj7OSAPqO8Ii0H3U9XF1ub8/9w3Ga49k4fj6GzJTsJ8rlv+WAslgta5aa49n+lpeoIQQQghPwuFEXSE8cp6cbf6V2eSFCFCkLhkUFnPxCeHOFKQMZ1erGRzjtRtXCOFTRC4ZJC0oEa5ygF0u1yPqm58Q4SRSlwySFpQIqXosc9LK3C6EIDaXDJIWlAgZM4dhNo5lTmwetufhGJm5wPwxZZj5GYUQLmL1b0kKlAiZei5zUoJMOSUEELtLBkmBEpbwY5mTBVy8wOTaGvcWQkT1kkFyDEpYxecyJ1rrIqXUKtO1kQ7U6HMXQgB1WzLIbmaaiIgvfFKghFVs3ja4LHNSp6mmhIhRNm8b6rpkULiRLj5hFTuyzIkQwWAnSv+WpEAJq8gyJ0IER9T+LUmBEpYIl2VOhIh00fy3JAVKWGm+27kao3DMNC+EqJuo/FuS2cxFyDTUMidCRLtY/VuSAiWEECIsSRefEEKIsCQFSgghRFiSAiWEECIsSYESQggRlqRACSGECEtSoIQQQaeUsimlZptLgy41rpSaYp434tdDinUyWawQIlTsPlZ+bYVjzS87cFxrvcAUlIL6Ts+jtZ5pVp99CMeyLSJCSYESwgJmBVSb1rrei8o5Z6wO1f6DxRSNfGC21nqm6+1KqYk4llSJuBm3RehIgRLCGvOCuK/xgPusAcHcf7AsA6aapR8u0FrblVLziYKpeURwSYESwgJBXjBuFG4FKtwWpDMtpOPuxcnJFCnpjhMXkUESQkQwZ1ee1Tn8MAlH954v4djqExaSFpSIKc5v8uZqquuEmn5sKzJXbTiWM1gN5OE40D/D3GYDhmmtp5pltgEygSLnCsFmVNtsAK31KH+e38tryTHPl24GHoCjJZXqun8z0Widc7o8j3MSUps/ubzI5Nv3zyPX53U5JmUHcrXWReb2pSb/3eauc3H8HPJxTJI6SmudG0A+EY601nKRS0xccHxIT3S5nu68jvmAc9lmA5aa/0/0sC3P/D8bWOO2Pd+53eW2E27XM537ry1bLa/pov342H+gOfOBTLecOX7kcn2P0gHtuh8/f14TcQyocL0tB8fgjxqvyzznFPfnl0vkXqSLT8QEM4Jsor74238OkGFaGOnafEuHCyuRFplv8vDtMgfObc7uqOM4PjBdWweeWgrHTQYnuz/Z/Hhp3tjdrtc5p2npZeqLj2fNw+W98Ifbc9blcXNwDABxv93ucvXC69Ja27XL6EAR+aSLT8SKbNw+kJ0fZm7dd652AUO11pOUUvlKKY1jldJ8t2Ji9/DYkmBkCzK7h9t85cwG7C5dgGC6FAN47iIgC0dX4QXO4oyjGE/E0UWZ51LUCpRSOdpxnpQNz68hoAIowp+0oITwPcgg1ZxnlAuk4OjiylVKuQ6JPu7hcfbgxfNfLbM21DWnDccxqQKXywKtdSAtuwVAjWNDLq0e5wm9eW4trhl822LL1p5HAdoDyCMigBQoESsurDTqynwrL/C0Dce3+qU4ZiRwfpgWaMfAhmBO3+MrW11l1jvNtzzmCoR2zCiRbrpTfbmoiJruxfSGni5JhAcpUCImmG/l812OKTllmw/BItcPT1McspxdeR4e5/otP9XDU9qCkc2PhxdxcRGx+7hvnXKa1spx96LiIae/RgF5XubIG4/3rrrZ5uLt3C5bgHlEmJNjUCJmmGNJU1yOOdm0Gdqstc4125wf9unAdeb/JTgKmLNg2HB80GbiaF2lK6WmaMcccDk4BjjYlVJrtdYFZph2unlMntlHHpCllJqotZ7jK1str8luJkZ1ngi7wLyGC/vHMQy7zjnNwINRJlcWpnWjAxtm7izEzv3lOd9XHIWzwOT0ZA6O4eMXFTCX9z/LZJ+j6zmPnwgvSjuGZAohRNA4J2vVHiaLDXB/Of4U7FA9v7CGdPEJIcKaUiq7LsVJRA8pUEKIsGO6LZ3HvmxWZhHWkQIlhAhH+TiG+Nepa09EFxkkIYQIlXSlVD6OJTbqdDKtl/Od/GIGTAxDJp+NeDJIQgghRFiSLj4hhBBhSQqUEEKIsCQFSgghRFiSAiWEECIsSYESQggRlqRACSGECEv/HyPSrfE7XyMeAAAAAElFTkSuQmCC\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.loglog()\n",
    "__=ax.plot(tarr, _mass_accretion_rate(tarr, TODAY, MILKY_WAY_MASS, x0, k, early, late))\n",
    "\n",
    "xlabel = ax.set_xlabel(r'${\\rm cosmic\\ time\\ [Gyr]}$')\n",
    "ylabel = ax.set_ylabel(r'$dM_{\\rm halo}/dt\\ [M_{\\odot}/{\\rm yr}]$')\n",
    "title = ax.set_title(r'${\\rm Milky\\ Way\\ halo\\ mass\\ accretion\\ rate}$')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "defensive-socket",
   "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
}