{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "6b8d7ed0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from load_bacco_data import load_bacco_mahs\n",
    "from astropy.table import Table\n",
    "\n",
    "drn = \"/Users/aphearin/work/DATA/BACCO\"\n",
    "\n",
    "bn = \"nenya_sigma8_0.730_ScaledFrom_nenya_mah.h5\"\n",
    "bn = \"nenya_sigma8_0.770_ScaledFrom_power_mah.h5\"\n",
    "bn = \"nenya_sigma8_0.815_ScaledFrom_nenya_mah.h5\"\n",
    "bn = \"nenya_sigma8_0.860_ScaledFrom_nenya_mah.h5\"\n",
    "mah_data = load_bacco_mahs(os.path.join(drn, bn))\n",
    "halo_ids, log_mahs, tarr, lgm_min = mah_data\n",
    "\n",
    "bn_fits = bn.replace(\".h5\", \".fits.h5\")\n",
    "fit_data = dd.io.load(os.path.join(drn, bn_fits))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09888334",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "id": "0ad7e280",
   "metadata": {},
   "outputs": [],
   "source": [
    "from diffmah.individual_halo_assembly import DEFAULT_MAH_PARAMS\n",
    "FIXED_K = DEFAULT_MAH_PARAMS[\"mah_k\"]\n",
    "logt0 = np.log10(tarr[-1])\n",
    "\n",
    "DIFFMAH_FIT_KEYS = ['logmp_fit', 'mah_logtc', 'early_index', 'late_index']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "id": "ccf293b5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEnCAYAAABCAo+QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzgklEQVR4nO3deVhb550v8O+RxI7hILCNsbGN8J7EjgU4ztLGsUU7k7ZJx5ZN22Qmc29tqDPpbTupYdzeTpw+t7GFZ3rbae4kkKR3muUmNpg2SZumkXCcpY1tQHYWx1sk8IY3kA4yi0BI5/5BUJEFMovQEfD9PA9Pe1b9ZEXnq/Oe97xHkGVZBhER0QAqpQsgIqLow3AgIqIgDAciIgrCcCAioiAMByIiCsJwIFKQ1WpFYWEhCgsLw7I/SZJQVlaG6upq2O12WCyWsOyXph6BXVmJlGW1WlFWVgaz2TzmfRUWFqKiogI6nQ4AIAgCrv+KS5IEURTH/Fo0ufHMgUhh4TxQOxwOfzAAQENDQ9A6e/fuDdvr0eTFcCCaJKxWa0AwAIBerw9aLxxnKDT5MRyIppCysjJIkqR0GTQBaJQugChcysvLodfrIUkSHA4HiouLUV5ejoqKCjgcDjQ2NkIURVitVuTl5cFoNMJkMvl/bVdXVwMA7HY79Ho9DAYDgL9eExBFEdu3b4fD4YAkSairq4PJZPJf9O3/5W40GoO2Kykp8e/bZrPBZDKN+L2EUl1dDbPZDKvVivLycgCAzWZDZWUlbDYbdDodqqurIUkS7Ha7f53i4mJef6DByUSTgNFolBsaGvzTpaWlclVVlSzLsux0OmWdTic7nU5ZlmXZZrP5l/UzmUwB0waDQbbZbP5ps9ks6/X6gHlGo1EuLS0N2E4UxYBps9kc8Nr98wwGg3/aZrMFTId6L6E0NDTIRqMxYN71NTc0NAS8FtFQ2KxEE57dbofVag1oXy8qKkJFRQWAvgu+FRUVWLduHYC+X/T9v+771dXV+c8cgL62+oHdQLVaLSRJCmjTv759f+B6A6d1Ol3Ar3ODwTBkN9MbvReiSGGzEk14FosFoigGHGz7m0/6GQwGGAwGFBYWDnpBtqqqyv//7XY7JEkKapsfrPklPT19VDXr9XpYrVZ/09VI3gtRJDAcaMLr/0V//YH2+rODwsJCWCwWWCyWoHXtdjtMJhPy8vJgMBgGDQKtVhs0L9zt9cN9L+Fit9sHPQMiYjjQhKfX67Fnz56Q69jtdmi1WtTW1iIvLw8NDQ0BB/a8vDz/BevrjUfvHqvViu3btwfNH857CXcdDAcaDK850IRnMBig1WphtVoD5ldWVgLoO7hbLBbo9XqIogiTyYSNGzf61+vfbmAw9AdCf28hAP7/HWg4wVFfXx+wXnV1NfR6/aD3INzovdzIjerR6XQBTVTsqURDUvqKOFG4mEwmuaKiQq6qqvL37iktLZV1Op1cXFzsX89sNssAZIPBELCeyWSSzWaz3NDQIDudTrm4uFiuqqry9wISRdHfq6mqqkrW6XSyXq+XzWaz//UByMXFxf4eQv29g8xms2w2m+WqqqqAnlE2my1o30O9l1AG1lhRUeHfhyiKAfVcv2+ioXBsJaJxFM5xk4giic1KREQUhOFANI4Gu05BNBEwHIjGidVqRUVFBerr6/3DVRBNFLzmQEREQXjmQEREQSbNTXAZGRmYP3++0mUQEU0oTU1NaGlpCZo/acJh/vz5qK+vV7oMIqIJJT8/f9D5bFYiIqIgDAciIgrCcCAioiAMByIiCsJwICKiIAwHIiIKwnAgIqIgDAciIgrCcCAioiAMByIiCsJwICKiIAwHIiIKwnAgIqIgDAciIgrCcCAioiAMByIiCsJwICKiIAwHIiIKwnAgIqIgDAciIgrCcCAioiAMByIiCsJwICKiIAwHIiIKwnAgIqIgDAciIgrCcCAioiAMByIiCqJRugAiIhoeWZbR4+tBr68Xvb5eeHp74LpwDrpF+rC/FsOBiKYkr88Lj8+DHl8PACAlNgUAcMZ1Bl29Xejx9vQt9/ZAjBOxNH0pAOCPjX9Ep6cTPb4e/zoLxAVYk70GAGA6bEK3txsen8d/EL8j6w783cK/g8frwcO1D/vne2Uven29uH/B/Xhg6QNo626D8XXjX5f7vOiVe7F1xVb8t5v/G863n8e9NfcGvI9fXmzB9K3HMC1VG9Z/H4YDEUUVp9sJqVtCZ28n3L1udPd2QxAE3J51OwDgwLkDOOM6g25vd99fbzdS41KxZfkWAMD/bvjfOOk4Cbe3b1u31w1dqg7/vubfAQDf+P03cMJxAl7Z63/N1bNW45kvPQMAKDGX4EL7hYCa1mavxS/X/hIAsPPQTji7nQHLv6b7mj8c3mx6Ez7ZhxhVDDQqDWJUMcgVcwEAgiDA3euGWqVGnCYOGpUGaqghdHlg++gvaLnaiNzuaVB53NB43FB7gRgfEP/Wr/Hpa8+hW3Zja+I1JMq98MlqtMTNR/eib0GjiQnvhwCGAxGNgizLcHvd6PB04FrPNbT3tGNZ+jKoVWocvXIUx1qPodPTia7eLv9B/rHbH4MgCPj1J7/GW01v+Zd19XZBBRXe/ca7AIAnDj2BN5veDHi96QnTsX/TfgBA9alqvHP+HQCAAAHxmnjkpub6w0HqluDqcSFOHYfU+FTMUM1A9rRs/76+qvsq7si6AzHqGMSqYhGjisHs5Nn+5T+67Ufw+DyIUcUgVt23XBuvhaenG20tl7B78XZ0u5zwdUhAhwvqnm6o7G58cPIHEDxd2OWdD8heCF4PBLkXgs8LwfZHWN/6AwT48M+yDyq5F4k9DqR6HdDKEjTCiwCAXAC3fV6HC4lwCSnwCPHwqGLRq4qDWp2KvJhbkHDTV7Bo9d8iPiEp/B/u5xgORFOUx+dBW3cbXN0uSN0SlmiXIDEmER9f/Rhvn3sb7Z52tPe045qn7+C/++7dyEjIwP/95P/iP6z/gV65N2B/7xW9BzFexIFzB/DcJ88BANSCGomaRCTEJMDj8yBWHYs4dRzS4tMwWzMbiTGJSNAkIDkm2b+fosVFuCf7HiRoEhCviUe8Jh4JmgT/8ie+8AQAIF4djxhVDARBCKjj8TseD/m+H1z2IABA9vnQ0d6G9rZWdLla8cl7r8LtvIhY10XEXrsMTddVxHe3INnjQKrPiRhcQwaAjCH265MFuBGLHiEGvdDABxV6BQ18UMMrqCFDgAw1fIIAGSp0xmjhnLYItqQZUE3LRIw4C4na2UiZng3tzDlISUxGyjA+x/EiyLIsK/j6YZOfn4/6+nqlyyAaM9nnQ0+PG93uruFvI8uQIUMlqODsduJIy4dwdDvg/PxXtMvjwkOLHoQuRYcDze/iiSMmdPR2BOzj2S8+jaVpS/Damd9j99GfIykmCckxSUjSJCEpJgn/qv8xMhNnouHqEdRdrQ9YlqxJQt50PeLUcejwdMAre4c8eI/p30aWAZ8XXm8vfD4vfF4v3B1t6JSuwt12FT3tLfC1t0DudELtdiCmR0Kspw2xPjc0vh7EyH1/8XAjWe6EWhj88OeWY+BQaeFSa9EVl46e+Az4kmZAlTwdmuR0xKVMR0LqDCSJGUhInIa4xGTExSVAUE28DqBDHTt55kCkgNbL5/HZgReRbPsDUjwtiJW7EYsexMk9iEcP4gQZcQC8AFrUalzWqDGj14tMrxcXNGpUiKloVavRqlahVa2GQ6XGE1db8OXOLnwaH4efzJoJABBkGSk+H0SfD731LyPF3Y3FsTH4u+SkvvleH1I/X77s+a9gmizjmwAeABB0SK//PQDgns//hqLkr91+PbIaLmEa2lUp6FSnoFMjwqeOg1cVB586DrImHr64VAgJqVAliFAniEhIy8S0jCyIM7KRPE1ElkqFLKXfiIIYDkTjTPb5cPHMKVw+dQjuc0cx7aoVS9wfYZXgw1HNXPxZnIcWjYDpwjTME7SQYlR4SvMJnIIbbXDDh75ft0WqFTCoFuGyfA1vew8gBfFIEeKRgzisQDzasubioCDCLXvwGDqQgngkIw6qz3+592YDBz+v6e5B6jwWmX+OsVOpIajUgKAGBBXUcUmITZmOBHEGksQZSNHOQFJyKjJUqiGbgOjGGA5E40D2+WA/dhhnD74I9dVa9GqcEL0+rHL3wK7Oxr3zF6NV09+Nsq9nzD8s+wcUFWxDV28Xqvd/F7ckZmJm0kzMTJyJzKRMLEpbhMykTADA/Qq+N5oaGA5EY+STfbjSeQXnLn6G3gP7kNTyIZ5OuYCTcSpcidUAszUApuP2lFvxH19+CgsTk1Hw/o8hxomYM20OZifPxqykWchK7mvESNAk4NkvPavsm6Ipj+FANAq/Pf1bHLp0CHbJjsa2Rri9bix19+L/NV/CGc08yJoZWBCfgft0d2HxrJswJ3kOsqdlIz6ur1fOz+76mcLvgCg0hgPRIHyyD2dcZ3Cs9RhOOU7hpPMknG4n9n5tL7o6ruHN46/huOsksnzJuOcakO92IM43Gxe++RJyl+jxlNJvgGiMGA40YbU5W3D6L79Fb+tZyF0OqN1OxHRLiO11QTXg7tdh7Uvlw2exPtzqVkMNAb9J7cYfpnkAABoZyPaoMN8j4NKOeciEhCcB9N+T6kISji/5Lgo2lkKlVof3TRIphOFAUa+r4xps1v1oP/kOVJ1XAABxnZewtMuKfKEvBHpkDdo+77rYpUmBVxUbcp+SyocjcR4cj/PgRGwvrmp8AIBdV1Iwr1eDVW4Bs3tjketRI6tXDQ0E+AQ1zopL0SjOQ0x6DpIzc5GRvQjpM+bgtgnYv50oFIYDRaXGY4dw6YNXIF4+hNyeE7hZ8KJXVkESUiBDgFtIgDVzE8SCjZi7JB8JidMwXaXC9CH2d851DgcvHcTK6SuxIG0B3jv/HiprH4Y2Xou8mXlYnrEct0y/BTdn3Iw4dRxujui7JYo+DAeKClcuNOLsETM8l08i/eK7WNR7CtmyCvaYBWjI+iYSF62BTr8OGQNGnswOsT+Pz4MPmj/Any/8Ge9feB9nr50FAHxf/30sSFuA/Mx8/O7+30GXqgvrHbxEkwXDgRR34tBbyPrjQ8hHJ7yygCZNDg4u2obFhd/Goumzhr2fK51XcKXzCm7OuBmyLOOH7/wQsixj1axVeGDpA7g963bMT5kPoK+7aP9ImUQUjOFAivro7WosPLAVLaoMXL5vD7KX5CE3IQnDPWxfaL8AyxkLLGcsOHr1KBamLUTNfTWIVcfiN3/zG+hEHeLUceP6HogmI4YDKab+D89gxeEynNHMR1rxa8ieOWdE2+86vAsvHX8JALBEuwSP3PoIDPMM/uX9D2chopFTLBwkSUJlZSVaW1thMpkGXQYAdXV1KCwsRHFxsRJlUhjJPh+OH/oTetqdcF/4BKsa/xMnYm/CnH96DSlieshtu73deOfcO3jd9jp+vPrHyEzKxBdnfxEzE2fCMNeA7JRQVyCIaKQUCQeLxQJJkmCz2QZdvnPnzoDAyM3ta2RgQExcjccOoevVR7Gs52P/vA8TV2HRIzVISJo25HZnXWex5+Qe/O6z38HV48KMhBk4d+0cMpMyccfsO3DH7DsiUT7RlKNIOBgMfaf+dXV1kCQpYJkkSbDb7QHzSkpKYDKZGA4TUJuzBcdf3o78y9W4JiTh0LL/ifTFt0OljsEtS/ND3jTmdDtx/+/6hphbN28d1i9cj9syb4NaxRvNiMZbVF5zsFgssNvt0Ol0AABRFIMCg6KX7PPh0Es7kHz+HcztPoVVchfqMu7Hkm+V47b0mUNu1+3txqufvYrjjuN47PbHkBafhie+8ATyZuZhRuKMCL4DIoq6cBBFEU5n4MO7zWaz/2yDot8n77+O1bZfwq6ajxNpa6FdsxW3rbhryPXbe9qx5+QevPDpC2h1t2L59OXo8fYgVh2Lv8352whWTkT9oi4cridJEiwWC2pra0Ou19zcHHAz02OPPYYdO3aMc3U0GO/Bp+BACrK2/QW6GzwA/YPmD/DogUdxzXMNd2Tdgc23bEb+zHzemEaksKgPhy1btqCqqgp6vT7kellZWWhubo5QVTSUC/ZjWN5xEIez/ztWDxEMPd4eXO68jOxp2ViWvgx3Z9+NB5c9iJvSb4pwtUQ0lKgOh/LycpSUlLBJaQI59+YvMQMq5N77P4KW+WQf/mD/A3515FdIiU3B3q/tRWpcKnZ+YacClRJRKFE7lGR1dTX0er0/GCwWi8IV0Y20u5y46fJr+DD1HkzPmh+w7ITjBB7640P40fs/Qlp8Gh7NfxQqIWr/8yOa8qLy22mxWOBwOJCfn+/v2mq1WpUui27g2BtPY5rQhZQ13w2Yf/jiYRT9vghnr53FT+/4KV7+ysu4Pet2haokouFQpFnJarXCYrGguroaQF/zkcFggF6vhyRJKCwsBNB3f0M/o9GoRKk0TD6vF1knn8dJzWIs1q8BALR2tSI9IR36mXo8vOJhfGPJN5Aal6psoUQ0LIIsy7LSRYRDfn4+6uvrlS5jyvrw7SqseGcz6vN3I9ewETsP78Thi4fx6tdfZSAQRbGhjp1RfUGaJoZrbQ4k/LkcV5GGtpt0+PqrX4erx4WS5SVIjElUujwiGgWGA43J1eYmuJ77Oub2nsFPbroPb777AyxKW4TKwkos1i5WujwiGiWGA42au6sDruf+Dpm9F3FyzTNw+95BUVIRthVs4zMUiCY4hgON2oe//h/wxjaj6bafYt09Rvzcdz9iVDFKl0VEYRCVXVkp+ln/9AJOdP8J35k5A5aEJgBgMBBNIjxzoBE7e+Y49pz+Kd5IT8M9c9bgJ6t/onRJRBRmDAcaEUd7C8re/CY+mRaPB7I3oPSef+WdzkSTEL/VNGztLidOPbkBPqETW5IM+Je1OxgMRJMUzxxoWD5tssL34mbkexrxg4wdWL3+e0qXRETjiOFAN1T/6QFs++AR3Ca68c1lT2H12m8oXRIRjTO2CVBIf2n4Pco++Cd4BBn3rPhXrGAwEE0JQ5451NTU4PDhw2N+IldBQQHWr18/pn2QMvbX/heeaCpHl6DC/1r+U6zN5+dINFUMGQ5msxm7du0a8wsUFxczHCYY2efDwZd/hqc7XkSHJhamVeX44s18ljPRVDJkOOh0OqSmjn00zfz8/DHvgyLj4pmTuGL/CN0f1uB26Q1cSSlA5vrHcdv8O5UujYgijEN2EwDg0tnTSHnuTsQI3difmIBp6d/C6n/cDZVarXRpRDSOxnXI7sbGRlitVqSlpcHpdEIQBDYlTTAX93wfKZDx/Vvux7vtR/DivQ8xGIimsDGHw9GjRwEAGzZsCJi/b98+FBYWIiUlZawvQePsw/2vYGXH+/ip7m/xbvsRfGfFd7Bi+gqlyyIiBY25K6vNZsOtt94aNH/Dhg2wWCxj3T2NM3dnO6a/9694NXEO9gnHsW7uOmxdsVXpsohIYWMKh6amJuj1ev/0v/3bv2H//v3+6XBc0KbxdeTlx5CCKzDNSoYuVYef3fUzDolBRGNrVnI6ncjNzfVPv/LKK2htbcXatWsBYMz3SND4OnPyKPLO/hc+Sl2HJ9Z8F7pUHZJikpQui4iiwJjCYeXKldi/f78/DK6/4u10OseyexpHss+H9n2PoDEmAfO/9QtkZGYrXRIRRZExtx8MFQCNjY3Q6XRj3T2Nk7pXn8QlzWfYlK3Fad95pcshoigTMhxcLtcNd7Bhwwbs27fPf62hra0N+/btQ1tbG1auXBmeKimsHFcuIPWT3fjx9OlYol2K/Jm8UZGIAoVsVtq5cyd27tx5w530d2Otra1Fenp6ULdWii6nX96Giunx8GniUX53OWLUfLwnEQUKeeZQXV2Nt99+e9g7W7du3aDdWil69Hp6cLL3A9QlxGPbbaWYlzJP6ZKIKAqFDIeWlhYYDAao1WoUFBRg+/bt2L9//5DNTcNphiJlfWY9AFnwYEVsDowLjUqXQ0RRKmQ4WK1WPP3009i7dy82bdqEhoYGGAwGpKWlYeHChdi6dStqamrQ1NQEAMNqgiJlSR+/iQfaOvDkl59jV2MiGlLIcMjJycGWLVsgiiJ0Oh3eeust+Hw+1NXVobS0FK2trdi8eTNyc3ORnp6OysrKSNVNo2A5Y4FNOoDTsYshaqcrXQ4RRbFh3eewbt06AH3jJeXm5kKv10Ov12PLli0A+rqtms1m7N69e/wqpTFx9bjw+F92YO60dszX3q90OUQU5UZ0n8OGDRuQmpqKmpqagOsLOTk5KC4uZi+lKPbkkSfR1uPCj1sdmH7rV5Quh4ii3IjvkM7JyUFOTg5qa2vR1tYWMDR3SUlJWIuj8DjhOIE9J/fgHreIrJ5WTFt+l9IlEVGUG/XwGQObmtLS0rB27Vrk5OSErTAKD1mW8bODP0NqbCq+13QWtmmrkKcJy2M8iGgSG9PwGS6XC7m5ubDb7SgqKvL3WqLoIUOGcZERm7M2Qedzwqtbq3RJRDQBhPwJ6XK5YLfb/X82my1gGuj7Zdrfm6miooLdWaOMSlDh/gX34+BffgIAyLntawpXREQTQchwEEURgiBAFEXk5ORAp9Nh5cqVMBqN0Ol00Ol0bEqKYq9+9iqkbgkPLn0Q086/A7tqPnRZ85Uui4gmgJDhoNPpUFJSAlEUYTAYGAQTSKenEz9v+DlyUnOwfs7XsLD7EzRkPQCOk0tEwxEyHIxGI7Zt2wYAOHLkCJ599ln/Mp1O53+OQ7+ampqA3kuknN98+hs43A78au2v8NmhN7BS8CLl5r9RuiwimiBChsOuXbv8/3/lypUBQ3A3NjbimWeegSRJEAQBqampKC8vZzhEAVePC88fex7r5q7D8unLcejkv6NDjsfC/HVKl0ZEE8So+zT2D63Rz2q1orW1NSxF0di8dPwltHvasXXFVsg+H+a2/hmnkvKwMi5e6dKIaIII25Pk9Xo9iouLw7U7GoP8mfkoWV6CxdrFOHv6I8zCVfTksAsrEQ3fkOEwmnsWBjZDjWU/NDYFmQV4ZOUjAICL9a8BALIL2IWViIZvyHAwmUxheYGysrKw7IdurNPTiV80/AItXS3+eYnnDuCMag6y5i9WrjAimnCGvOYgyzK2b98+pp3LsgxZlse0Dxq+3372Wzz3yXNYk70GGQkZcHe2Y3HXRziSuQF83hsRjcSQ4RCuMweKDK/Pixc/fRG3Tr8Vt864FQBw4bOPkCt4EDNvtbLFEdGEM2Q4pKamRrIOGqO3z72N8+3n8c/5/+yf57p4GgCQOnuhUmUR0QQVtt5KpKznP30es5NnY232X3sldV9tBABMn7tEqbKIaILi2M2TQLe3G7OSZuHenHuhVqn98wXpDFxIQmpahoLVEdFExHCYBOLUcTB9MfgaUUL7OVxRZyJFgZqIaGJjs9IE19bdhtPO04MuE3ua4UqYHeGKiGgyYDhMcDWna7D+tfU4d+1cwHyf14tM72X0JGcrVBkRTWQMhwnMJ/tQdaoK+hl6ZE8LDIGrF5sQK/RC0HKYdSIaOcXCQZIklJeXD3kH9Y2WE3Do4iGcu3YOmxZvClrWer6vqSlxBp/gQEQjF7YL0rW1taiurgYAFBYWhhy622KxQJIk2Gy2US2nPlWnqiDGiSicVxi0rOPSZwCAtDmLIl0WEU0CYQmHffv2wW63w2g0AugbvvvZZ5/F5s2bB13fYDAAAOrq6iBJ0oiXU984Sh80f4ANCzcgVh0btLy3tRE+WcCMbN4AR0QjF5ZwEEXR/8Q4AFi3bh1qa2vDsWsaQmJMIsxGM3p9vYMu17jO4oqQjkw+w4GIRiEs4SAIQtC8tra2cOyaQkiOTR56Wed5OGJnITOC9RDR5BGWC9JOpxPbt29HTU0NampqsH37djgcjnDsetiam5shCIL/b8eOHRF9/Uj6zPkZNr2+Ccdbjw+5TobnItoT2Y2ViEYnLGcOGzZsgE6nw549ewAARUVFAc+bjoSsrCw0NzdH9DWV8rr9dZxynsKMxBmDLnd3tmM6nLClzo1wZUQ0WYw6HK5/wltaWhq+853v+Ke3b9+OnTt3jrowGpzX58Xv7b/HnbPvRHpC+qDrXD57CvMAxGTwHgciGp1Rh4PBYEBeXl7Aw3z6rz3IsowjR44wHMZB3eU6XOm8gm0F24ZcR2ruC4ekzAWRK4yIJpVRh0NFRQXWrVs35PIjR46MdtcUwhv2N5AUk4Q1c9YMuU7X5b6hujOyeY8DEY3OqMMhVDAcPXoUjY2NQ153sFqtsFgs/pvmysvLYTAYoNfrh7V8Kls9azVyUnMQrwnRRVVqQpcci/QZcyJXGBFNKoIchoc8NzY2oqSkBE6n09/MVFhYGNFmpfz8fNTX10fs9aLZG//xCDLaPsGqnxxQuhQiinJDHTvD0lupoqICFRUVAACtVgsAsNvt4dg1DVB3qQ45qTnISAj98J6XEh+EO8aHfRGqi4gmn7Dc51BUVIScnBzk5OSgtrYWqampcDqd4dg1fa7X14tHDzyK8sPlN1zX0eFBWmLwkBpERMMVljMHh8OBgoIC1NXV4fDhw/4xkdauXXvjjWlY6i/Xw9ntxJfmf+mG6zo7enBzFp//RkSjF5ZwWLduHerq6gAAu3btwr59+3jxOMzeanoLCZoE3DX7rpDrybIMZ2cPtEk8cyCi0RuX5zls2LBh0PGWaHR8sg+1Z2vxhdlfCN1LCUCXx4vuXh9ENisR0RiE5cyhqakJVVVV/vGUZFlGbW2t/2yCxua08zQcbgfWzr1xM52z0wMA0CbFjHdZRDSJhSUcdu3ahcLCwoCmJD6HIXwWaxejdmMtkmOGHoW1n7OjBwB45kBEYxKWcCgpKQm64S0/Pz8cu6bPDTXI3vWcnX3hwGsORDQWo77m4HK5/H9OpxP79+9HU1OT/2/Xrl3hrHPKutRxCY/UPoITjhPDWt/x+ZkDu7IS0ViEPHNwuVxISRm8S6Rerw8aeG8gDrwXHu+efxfvnH8HP8j7wbDWd/rDgdcciGj0QobDzp07hzzAc+C9yHj3/LuYnTwbulTdsNZ3dnogCEBqAsOBiEYvZLNSdXU13n777UGXhQoGABF/2M9k1NXbhYMXD2JN9pphdw12dvYgJT4GGvW49FImoiki5BGkpaUFBoMBarUaBQUF2L59O/bv3w+XyzXo+kPNp9Gpu1SHbm83vjj7i8Pextnp4cVoIhqzkOFgtVrx9NNPY+/evdi0aRMaGhpgMBiQlpaGhQsXYuvWraipqfE/FY7XGMJPP0OPvMy8Ya/v7Ojh9QYiGrNhDdldW1sLSZKwYcMGAH2h0dDQALPZDIvFgra2NoiiCABobW0d14KHwiG7+9z7y/cwKzUez/1jgdKlENEEMKYhu/uvL+zbtw+5ubnQ6/XQ6/XYsmULgL7nOZjNZuzevTuMJU9t7l43BEFAnDpuRNtJnT1YxkH3iGiMRnTVcsOGDUhNTUVNTU3A9YWcnBwUFxf7zyxo7MxnzLjz5Ttx1nV2RNs5OtmsRERjN+I7pAc+t6GtrQ3r16/3LyspKQlrcVPZwYsHkahJxJxpw3/UZ1ePF26PD2m8IE1EYzTmZ0jv27cPaWlpWLt2LXJycsJW2FQmyzIONh/EbbNug0oY/smdf+gM3h1NRGM0ps7wLpcLubm5sNvtKCoq8vdaorGxt9lxpesKVs9aPaLtHBx0j4jC5IbDZ9jtdv+fzWYLmAb6fuWKogidToeKigp2Zw2DgxcPAgBWZ40sHCT/cN0MByIam5DhIIoiBEGAKIrIycmBTqfDypUrYTQaodPpoNPp2JQ0Dm6fdTv+ZdW/YHby7BFt5+jkuEpEFB4hw0Gn06GkpASiKMJgMDAIIkQn6qAThzeW0kBSfzjwzIGIxihkOBiNRmzbtg1A30B6zz77rH+ZTqfD2rWBTyarqakJ6L1EI3ep4xKOtx7H6qzVSNAkjGhb/zUHDrpHRGMUMhwGPpNh5cqVAYPpNTY24plnnoEkSRAEAampqSgvL2c4jNH+s/ux8/BO/GnDn5CQPLJwcHb0ICVew0H3iGjMRt2VNScnx3+HNNA3pIZSQ2dMJg2XGzAraRaykrNGvC0H3SOicAnbT0y9Xo/i4uJw7W5KkmUZDZcboJ+pv/HKg3B29rAbKxGFRVjbH/ho0LE5e+0sWt2tyJs5/FFYB3J29vDMgYjCgo3TUeTIlb6n5+XNGGU4dHj47GgiCotRX3Og8Lsv9z7cnH4zclJH12XYyUH3iChMGA5RRCWosCBtwai2tZ51orPHi1niyHo4ERENhs1KUaK1qxU7/rIDp5ynRrytx+vDj2o+RmZKPIoKssehOiKaahgOUeKjqx9h3+l96PB0jHjbX7/fiBOXruHx+29CchxPBolo7HgkiRIftXwEjaDBUu3SYa1/ztGJX+0/jc4eL2qPX0Hhspn48k2Z41wlEU0VDIco8fHVj7FIuwjxmvgbrvuBrRUPv9SA7l4fMlPjccvsVDx+300RqJKIpgqGQxTw+rz4uOVj3Jd73w3XffHgGex47RjmpSfi2YcKkJORFIEKiWiqYThEgZauFqQnpGP59OVDruPx+vDT1z/FCwfP4J7F0/HLb65ESjy7rRLR+GA4RIGZSTPxxvo3IMvyoMudHT14+CUrPrC3ouSLOpT+zRKoVUKEqySiqYThEEUEIfiAf+ryNWz+TT0utbnx7xtXYEPeHAUqI6Kphl1Zo8DmP21GxYcVQfNrj1/G+v/8C7o8XrxSsprBQEQRw3BQWIenA4cvHYYPPv88WZbx1AEbNj9fj5yMJLz2yJ3Qz01TsEoimmrYrKSwE44TkCHjpvS+rqhujxfbaz7Gb49cwFeXz8Ju4wokxKoVrpKIphqGg8I+bf0UALAsfRkuu9wofqEBH56T8MMvLcI/3bNg0OsQRETjjeGgsOOtxzEjYQYutGhQ/ML7uObuRcXf5/FuZyJSFMNBYQvSFkC6loBNFR8gIzkO+7begaWzUpQui4imOIaDgnw+GS3n78Ab79qwKkfEUw/okZ4cp3RZREQMB6W0d/fiu68cxNvHHfjmqhw8ft9NiNWw8xgRRQeGgwLOtnZi8/N1OOMxQ1z2Bh79m7cYDEQUVXhEirAPbK24//+8j8uubqy5pRupcdOQnpCudFlERAEYDhH0wsEz+PvnDiE9OQ6v/tOdaPHYsTR9KburElHUYbNSBAw2ompcjA82yYa759ytdHlEREEYDuMsYETVu3Uo/XLfiKrHWo7BK3uxNH14T34jIookxcJBkiRUVlaitbUVJpMpaHl5eTl0Oh0cDgcAoLi4ONIljpl/RFWXGz/ftALr9X8dOE8br8X39N/D8oyhn+FARKQURcLBYrFAkiTYbLZBl5eVlaGgoABGo9E/XV1d7Z+eCGqPX8b3XjmKhFg19hSvxsrrBs6blTwLm2/ZrFB1REShKXJB2mAwwGg0QhTFQZdXVlYGBEFRUREqKoKHtI5Gg42oen0wAMCx1mNwuB0KVEhEdGNR11vJarUGzRNFERaLRYFqRsbt8eIHe47C9OYJfHV5FvaW3I5ZqQmDrvuw5WH8ouEXkS2QiGiYou6CtMPhgFarDZh3/XQ0uuxyo/j5enx4vu2GI6q2dLXA4XZgUdqiCFdJRDQ8UXfmIEnSqJY1NzdDEAT/344dO8Je21A+PCfhviffx+kr7aj8+zw8snZhyHsXTjlPAQDDgYiiVtSdOYii6O+h1O/66cFkZWWhubl5vMoa0qtHL2Bb9UeYMS0ONQ/fgSWZNx5R9bTzNABgYdrC8S6PiGhUoi4ctFpt0BlC//RQF7CV4PPJ2P3WSTx1wIZVOdoRjah6ynkK0xOmIy2ej/4kougUdeGg1+uDQsDhcMBgMChT0CDau3vx/VeOwHL8Cr5121zs+NrIRlT99s3fxldyvjKOFRIRjU3UXXMAgE2bNqG6uto/bTabUVJSomBFf3W2tRPr//PPePvkVfz0/pvws6/fPOIRVXWiDnfMvmOcKiQiGjtFzhysVissFos/AMrLy2EwGKDX6wEAFRUVKC8vh8Vigd1uR25ublTcAPeBrRUPv9QAnww8/99X4c4FGSPex9XOq/hz859x95y72axERFFLkGVZVrqIcMjPz0d9ff247f+Fg2fw+GvHMD8jCc89lI956Umj2s9bTW/h0XcexZ6v7sGy9GVhrpKIaGSGOnZG3TWHaOPx+vD468fw4sGzWLtkBn75jVsxLT5m1Puzt9khQEBOak4YqyQiCi+GQwhDjag6FvY2O7KSs5CgGfzOaSKiaMBwGMKpy9fw7d/U4bKrO2hE1bFobGvkWQMRRT2GwyAsn17G9145gsQ4zaAjqo6WT/ahqa0JqzJXhWV/RETjheEwgCzLeOodG3b/6SRuzkpF5T/kDTlw3mgIEPDHDX/EJOkDQESTGMPhc26PF2X7PsKrR5vxtRVZ2G1cjvgYdVhfQxAEZCSMvPsrEVGkMRwAXGpzo+SFvhFVt315MR5ekxty4LzReu/8ezjuOI5v3/xtqFXhDR4ionCKyjukI+no5yOqfvb5iKqhhtoeK8tZC146/hKDgYii3pQ+c/B4ffjuy1bEalR4/tvDG1F1LOySnT2ViGhCmNLhEKNWoeLBfGSmxkObFDuuryXLMuxtdnx5/pfH9XWIiMJhSocDACzLGt+zhX6t7la4elzQpeoi8npERGMx5a85RMrF9ouIUcUwHIhoQpjyZw6Rcsv0W1D3QB1k8B4HIop+DIcIYi8lIpoo2KwUIU8eeRIVH1YoXQYR0bDwzCFCzGfM7MZKRBMGzxwiwOvz4ty1c5g7ba7SpRARDQvDIQIudV6Cx+fB3BSGAxFNDAyHCDjrOgsAmJcyT+FKiIiGh+EQAd3ebsxOns1mJSKaMHhBOgLWZK/Bmuw1SpdBRDRsPHMgIqIgDIcI2GrZisqPKpUug4ho2BgO48zr8+LQxUNo97QrXQoR0bAxHMbZxY6L8Pg8mDeNPZWIaOJgOIyzs9f6urHyHgcimkgYDuOs/x4HdmMloomE4TDOUmJTsHrWasxInKF0KUREw8b7HMbZvbp7ca/uXqXLICIaEZ45jDNZ5sN9iGjimfLhsGPHjnHbtyzLWFe1js9xGIbx/BxoZPhZRAelPwdBniQ/bfPz81FfXz/i7QRBGLdf923dbbjrlbvww/wf4qGbHhqX15gsxvNzoJHhZxEdIvU5DHXsnPJnDuPpfPt5AMDs5NkKV0JENDIMh3HU3N4MgOFARBPPpGlWysjIwPz580e8XXNzM7KyssJfEI0IP4fowc8iOkTqc2hqakJLS0vQ/EkTDkREFD5sViIioiAMByIiCsJwICKiIFN2+Izy8nLodDo4HA4AQHFxscIVTT2VlZVoaGjAxo0bAQBVVVUoKyuDTqdTuLLJT5IkVFZWorW1FSaTKWg5vx+REepzUPr7MSXDoaysDAUFBTAajf7p6upq/zRFzt69e1FZWQm9Xo9nnnmGwRABFosFkiTBZrMNupzfj8i40ecAKPv9mJLNSpWVlQH/oRcVFaGigkNcKMHpdEKWZTQ0NECv1ytdzpRgMBhgNBohiuKgy/n9iIwbfQ6Ast+PKRcOVqs1aJ4oirBYLApUQxRd+P2gflOuWcnhcECr1QbMu36aIqeyshJarZZt21GC34/oouT3Y8qFgyRJIZeFOsWj8MrPz4coiv521I0bN0Kr1bJtW0H8fkQPpb8fU65ZSRRFfwr3u36aIkOv1wdcYCsoKMDOnTsVrIj4/YgeSn8/plw4aLXaoF9H/dP8VRRZ17dj63S6Qdu8KXL4/YgeSn8/plw46PX6oP/IHQ4HDAaDMgVNUXa7HYWFhUEHInZlVRa/H9EhGr4fUy4cAGDTpk2orq72T5vNZpSUlChY0dSj0+lgMpkCDkR79uxBWVmZckURAH4/okE0fD+m7Kis5eXl0Ov1sNvtANhLRgl2u91/EGptbUVubi4/hwiwWq2wWCz+exdKSkpgMBgC+tHz+zH+bvQ5KP39mLLhQEREQ5uSzUpERBQaw4GIiIIwHIiIKAjDgYiIgjAciIgoCMOBaJKRJAklJSUoKSnxd0WNlPLycpSUlATcJ0ET05QbeI9oKhBFccgnvLW2tiI9PR2iKPoHcquurobBYBjzEBmlpaWQJAk7d+7kAIoTHMOBppyysjJIkhSWB9gMNlJpOPcfLpIkYePGjSgpKUFpaWnA/MrKSphMJjQ0NChYIUUbhgNNOUVFRWHb1969e4PuWg3n/sNl3bp1MJlMQWMkiaKITZs2cXgMCsJrDjTl6PX6sD1y0Ww2j+v+w6H/gTFDDZ4niiKbgCgIw4FolPqbj6JdRUUFNm7cGHKdaDzbIWWxWYkipv8XLNA3DPTA5pgbLesfqliSJDgcDuTn56OsrAyiKGL79u1wOByQJAl1dXUwmUz+sfCtVit0Op3/l7Hdbvc3oQz81R/q9QdTXV0NSZJgt9tRXl4OoG9wOofDEbB/q9U6qjr79Q+A1/++RzPwWv++Qxn4uv3XIERRRFVVlX/bwsJCOBwOPPPMMwCALVu2ID8/Hxs3boTdbofZbEZVVdWI66MoJRNFQGlpqVxRUeGfttls/mmj0SjbbDb/MqfTKRsMBlmWZbmioiJoWWlpqSzLsmw2m2W9Xh+w3Gg0+pf3E0UxYLqhocG//xvVFsr1+xlq/mjrNBqNckNDQ0CdVVVVN6xr4L+RzWaTAQTsZzgqKirk4uLigHlVVVWy0+n0Tw98X06nUzaZTEGvTxMXm5Vo3PX3iBn4q7e6uho2mw1WqxV2uz3gl23/c3MrKysBIKDXjyiK/iaQ/qeWDdx2sF/I1z/dbGDvolC1jdb1vZdGU6fdbofVag24dlFUVDTiHlCjfThMcXEx9u7dGzR/4Hsb+L5EUQzoBUUTH8OBxp3FYgk6SJWWlsJkMqG+vn7QA1hubi4aGhpQXFwMu90OQRBQWFiIysrKgAPmYP3y09PTw1JbOI20TovFAlEUYbFY/H92u31UN7XpdDrU19cHzZckyX/TmiAIQTfNGQwG/81sg3XZ7d83TU4MB1JUqAu6/e3zVVVVcDqdKCsrQ1VVVUC3y/7rBAMp9azjUAfukdbZ/4vcYDD4/4xG46jOaIxG46DXAvp/7fcHYVlZWcDBfvv27f4zFYvFMmhvJz5XevJiONC4G/hEsYEkSYLBYBh0mc1mQ2FhIXbu3Amg7yBkMBhgNpvDOiREqNpGKpwPfx+qrtEwmUz+ZqpQrg+w/hoiPQQHRQeGA407nU6HTZs2+a8h9LNYLNDr9dDpdAEHLkmSUF9f778OcP12A3/dOhyOoNcbyYE9VG3D2XbggTPUr+iR1mkwGKDVaoMO6NfXOVxmsxllZWWDjnm0d+/eIZuH+sdoGuq+jYnQlZdGSekr4jR1mEwmuaKiQjabzUG9bkwmk1xVVSVXVVXJJpPJ3yvGZDLJZrPZ/1dVVSXbbDa5oaFBNhqNsiiK/l4yVVVVsk6nk/V6vWw2m/3bA5CLi4tlm80m22w2/3YDeySFqm0476l/m+v3P9o6B9v/cOsK1VvIZDLJpaWl/n/v/t5gDQ0NAT2RBu5rqB5ZA9/XwG3ZW2ly4DOkiSaZ/oHvwnVRvbq6ekR3UIf79UkZbFYioiFZLBYOrTFFMRyIKEBJSYn/WgevKUxdDAciCrBx40Y4HI4RNyfR5MKxlYgmIbvdjo0bN8JkMo34RrWhRm8djvLyctTV1XEgv0mAF6SJiCgIm5WIiCgIw4GIiIIwHIiIKAjDgYiIgjAciIgoCMOBiIiC/H8Ry9CRJL25MwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from diffmah.individual_halo_assembly import _rolling_plaw_vs_logt\n",
    "from bacco_fitter import get_bacco_loss_data\n",
    "\n",
    "tpred = np.linspace(0.1, tarr[-1], 1000)\n",
    "lgtpred = np.log10(tpred)\n",
    "\n",
    "itest = np.random.randint(0, log_mahs.shape[0])\n",
    "# itest = 446\n",
    "# itest = 1841\n",
    "# itest = 1151\n",
    "# itest = 4949\n",
    "# itest = 2605\n",
    "# itest = 2052\n",
    "lgmah = log_mahs[itest, :]\n",
    "p_init, loss_data = get_bacco_loss_data(tarr, lgmah, lgm_min)\n",
    "logt_target, log_mah_target, logt0, u_k = loss_data\n",
    "\n",
    "logmp_fit, logtc, early_index, late_index = [fit_data[key][itest] for key in DIFFMAH_FIT_KEYS]\n",
    "log_mah_pred = _rolling_plaw_vs_logt(lgtpred, logt0, logmp_fit, logtc, FIXED_K, early_index, late_index)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ylim = ax.set_ylim(lgmah[-1]-3, lgmah[-1]+0.5)\n",
    "__=ax.plot(tarr, lgmah, label=r'${\\rm entire\\ MAH}$')\n",
    "__=ax.plot(10**logt_target, log_mah_target, label=r'${\\rm target\\ data}$')\n",
    "__=ax.plot(tpred, log_mah_pred, '--', label=r'${\\rm Diffmah}$')\n",
    "\n",
    "title = ax.set_title(r'${\\rm example\\ fit}$')\n",
    "xlabel = ax.set_xlabel(r'${\\rm cosmic\\ time\\ [Gyr]}$')\n",
    "ylabel = ax.set_ylabel(r'$M_{\\rm halo}\\ [M_{\\odot}]$')\n",
    "# leg = ax.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "a38740d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.0, -1.0)"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "early_index, late_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2ae4b967",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68a9a36b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}