{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0bf929c5-9741-419c-8b2c-236304df443e",
   "metadata": {},
   "source": [
    "# Compare diffmah fits to Bolshoi-Planck MAHs\n",
    "\n",
    "Change the path in the next cell to the data location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "fec79764-4dd8-4608-9489-a4ac9125121a",
   "metadata": {},
   "outputs": [],
   "source": [
    "drn = \"/Users/aphearin/work/DATA/random_data/0228\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3d6d1a39-5721-4d96-b960-e597781a685b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['bpl_t_table', 'halo_id', 'mpeak_history'])\n",
      "dict_keys(['early_index', 'fit_algo', 'halo_id', 'late_index', 'logm0', 'logtc', 'loss', 'n_points_per_fit', 't_peak'])\n"
     ]
    }
   ],
   "source": [
    "from diffsky.data_loaders import load_flat_hdf5\n",
    "import os\n",
    "mph_bname, fit_bname = \"bpl_mpeak_histories.hdf5\", \"diffmah_fits.hdf5\"\n",
    "\n",
    "bpl_mah_data = load_flat_hdf5(os.path.join(drn, mph_bname))\n",
    "print(bpl_mah_data.keys())\n",
    "diffmah_fit_data = load_flat_hdf5(os.path.join(drn, fit_bname))\n",
    "print(diffmah_fit_data.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "37798e41-5829-44af-b9cc-c91cd24dc039",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "696_053 total halos\n"
     ]
    }
   ],
   "source": [
    "n_halos, n_t  = bpl_mah_data['mpeak_history'].shape\n",
    "print(f\"{n_halos:_} total halos\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6ba4a303-a4a4-4489-bafa-717889fc8f5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.allclose(bpl_mah_data['halo_id'], diffmah_fit_data['halo_id'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "196fe00e-6b81-4fb0-85fc-16c6c587c42f",
   "metadata": {},
   "source": [
    "### Mask halos without a successful diffmah fit and compute the best-fit MAH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "9505ea45-9cb1-4edf-80e9-7c774b6cea5e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99.74% of halos have a diffmah fit\n"
     ]
    }
   ],
   "source": [
    "from diffmah import mah_halopop, DEFAULT_MAH_PARAMS\n",
    "\n",
    "msk_has_fit = diffmah_fit_data['loss'] > 0.0\n",
    "print(f\"{msk_has_fit.mean()*100:.2f}% of halos have a diffmah fit\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "44736a90-38dd-4a13-a5f0-baac782b6ca1",
   "metadata": {},
   "outputs": [],
   "source": [
    "mah_params = DEFAULT_MAH_PARAMS._make([diffmah_fit_data[key][msk_has_fit] for key in DEFAULT_MAH_PARAMS._fields])\n",
    "dmhdt_has_fit, log_mah_has_fit = mah_halopop(mah_params, bpl_mah_data['bpl_t_table'], np.log10(bpl_mah_data['bpl_t_table'][-1]))\n",
    "\n",
    "log_mah_pred = np.zeros((n_halos, n_t))\n",
    "log_mah_pred[msk_has_fit] = log_mah_has_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "802963ef-40b9-4422-8ca6-db9969adf24d",
   "metadata": {},
   "source": [
    "### Plot some example halos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "b390161a-aec1-4a88-962f-41365c18e758",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEnCAYAAABsR64CAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAALEwAACxMBAJqcGAAALb9JREFUeJztnQmUFNX9ti8iboBhU0EQdHBDFhFQ4xpxACOKRmRJTBA1AhGMOyCav7gdYVDRJG4DStS4BAb3aJBFxQ3DpoAiLihukSgwqKgo4Hznvd+5k5qa7p6e6eru6qrnOafPdHd1V92ugvvW/a31KioqKgwAAEAGbJfJlwEAABATAAAIBFYmAACQMYgJAABkDGICAAAZg5gAAEDGICaQNZYuXWp69+5tH/ncR22YMmWKad++vZk5c2ZOjhcFPvjgg3wPAUJAPfJMIJtIDMaOHWvmzJmT1304tJ+NGzea0tLSpJ8ZMWKEFa8BAwaYXCERE6tXr7aT89SpU02TJk2qbRca//Dhw6ts13szZsyw56isrCzh/vUZfUd/x4wZU+0zkyZNstubNWtmX3t/f6rx6VzpvPbo0aPKmLzUNL6ajp/uPiB/bJ/HYwOkRbIJqi4MHjw4p8dLB03UEgeHVkXdu3e3E7ebZP3iIcFzgiix1QSvSTjRKsFN0k5AtH9N/iUlJZWfkSBognbHaNq0qenVq5d9XdP4Fi9enHDlqO+Wl5fXOL6ajp/Ob4T8g5kLYkW3bt3sIyxoYnSTskN35Bs2bKg0tS1atKiawLkVhtDv0XeKiooSHkPC4RUDfda70pHYDBw4sMoxlixZYl+nMz6JkgppeB/6zrx589IaX6rjO2raB+QfxAQgz3gndofuwDVhC03oc+fOrbLdmaxqQnf0iT6nSdntc8KECWbQoEHVtqczPo3D/1133HRFu6bjQ2GAmQtygiYXN/noTttrYhHuLlcTpyYhmThSoQnO2da1X++ddzK0b5mHhNf/4vbl7vbdHX8u0KQpU1CiscoHIXSuZAaSmUrPda7c76gJJ0iJ0DHc79XnJC46B7pWzqxW0/gSCZXfLJaKmo4PhQNiAllHE48mBnenOn36dDtxOMGQmcPrENbEqUks2d2pTCKaVN12TUb6Tk0Oen1e35PZx6Hnit7yOnt1p5xO9Fi6E7pI5fBPNBnr3Ljzpecy+8hPoXOl35nuXb8m/ETi6IRE/g4nou4c6Ds6x8nOp398yUQwHepyfAgniAnkBK8w6LnXiaqViu623WSiSUpik+ju1jlivftzd9Dp3BH7I6D0Hf+dd7qTYW0EIl3027RfiYf3PQmwxumETp9J5+5fv1ef84q3zqH3/Ok8+M+nVgqJTFWJxudF1zFRpFgqanN8CC+ICWSdmswV3jBPr+kl2Z1sohWLVhfJJrhkaIINm21eKyXnuHZoNeWES88VkVZcXGzHXpM5UOi7WtE4nHlQE7V77r9Gel/nxz+ZJxqfQ5P/+vXr0/6t3rGke3wIL4gJZB03YSRDAqJJUmYcbzhoIlL5M1L5BwoBF67r/f2aUP0mN02wyvOQGSgdMRH+1YJb3aW6Nv7rkGh8ftHSNawNqcQcn0lhgZhA3tEE9OGHHyacPPzioclTJh8/CkWtbZa8JuVMchZk108HTdg1mcRkbpMPxju5+iO46jp2v7lIr2XKc8dy+/Ie2+/7SDY+r5glEr6acL60mo4P4QcxgayTaMXgREITWyJfRvPmzSujfLxo4tGk450gnSO5tj4M7UchqYkmxXTusIPKwtbxvJO797xoXBItfza4fBN+n0mylZlzZrv9y+/iPVfjxo2rYkpzPhV3flONL1GgRTKSja+m46ezD8g/lFOBrKFJQROXJiNNGDK1uOxrTTqaQDRZ6rXEw9nwNZHoPd3l6j099+5DyAfgJjdNYv5Q0nr16iUckyZl7UvHdpOx25eLKtLEK3Fy48smGrv8PYmQw92NSefRfc5FPnl/v86rxq3fpnOkz7rfp/ecMPu/69D33UpHfg8Xup3O+Bz6nFe0vL8x1fhSHb82+4D8gpgAAEDGkAEPAAAZg5gAAEDGICYAAJAxiAkAAGRMbEODW7RoYfbee+98DwMAoGBYs2aNWbduXcJtsRUTCYnCPwEAID1SJZJi5gIAgIxBTAAAIGMQEwAAyBjEBAAAMgYxAQCAjEFMAAAgYxATAADImNiKyX/+8x9bpvzqq6/O91AAAAqe2CYt7rnnnlZQAAAgc2K7MgEAgOBATAAAIGMQEwAAyBjEBAAAMgYxAQCAjEFMAAAgYxATAADIGMQEAAAyBjEBAICMQUwAACC+YrJx40YzduxYs3Tp0irvT5kyxT5GjBhhPvjgg7yNDwAgThSsmCxevLiaWEhY1PB++PDhZuDAgVZQAAAg+xSsmPTq1cs0a9asynsSl9LSUvtcoiLBAQCAiFcNlqlKJqn169ebkpKSatsnTZpkioqKzIYNG+xrrThSMWDAAPsQEhIJCgAARHhlMnfuXPtYvXq1FRU/8odISCQOEhF9bubMmWnvXyuUsrKygEcNAAChEhOZqSQUTZo0SbhdKxa3yhCDBw+uNGHVhL6rlU6yfQMAQAx8Jv4ILSFh0EqmJvQZCZVWNel8HgAAItppUT4Sv3Pd/1pCIb+IW31069bNipCiuNxn9Z6EBQAAYigmiXwo3m0SEInEkiVLqmyTeJSXl9eqB7xj/Pjx9IMHAIiSmEgsXASXw/86U+gBDwAQcZ+JzFT+1Yl7jVMdACB8hFJMZK7yi4ZWJkH6P5yZ6+qrrw5snwAAcSWUZi4xaNAgm1fiwoPnzJkTaHkUzFwAABEQE0VeKSLLJSIq210rD61KhHJK9J4+ozIp7du3r5J3AgAA4aFeRUVFhYkhWpl8/vnnRHEBAKRJqpqHoTVzZRvMXAAAwRFbMQGAzJFh47333jPLly+3jy+++ML89NNP9v1EfyH/7L777ubmm28OfL+ICQDUiR9++MEMHTrUTJ8+3b7ebrvtTIsWLexfPRQt6f/rTRSG/NC2bdus7De2YuJCg8l8B6g9X3/9tfnlL39pFixYYPr372/2228/s/3225uddtrJ/OlPf7KfGT16tFm4cGGV7+29997mvvvus89HjRpl3nzzzSrbDzroIHPnnXfa52effXa1Bniy2bu7ahV/Xbt2bZXtxxxzjLn++uvt81NOOcV89dVXVbafcMIJ5oorrrDPe/fubX788ccq23/1q1+Ziy++2GzdutUUFxdX+91nnHGGjSrV7+/Xr1+17eeee64ZMmSIHZfG5+eCCy4wp59+uvnwww/NWWedlbBaet++fc1bb71lRo4cWW27Uhl69uxpFi1aZC677LJq2xW0dPjhh5sXX3zR/N///V+17bfddpvp3LmzyQaxFRN8JhBmvv/+ezvhqPXCmjVrzGeffWa++eabrB1PJqjNmzfb4zZv3ty+99FHH9njfvfdd5WPbdu2md/+9rd2stKEJx599NHK/UhUnJi4FYkX72u3gvHiXbnUZXtt9p9ofP79+8l0u5dU3080ttpsT7b/bBLbaC46MUKY0CQ9e/Zs88gjj9i/n3zySZXtDRo0MLvuumudzETOX6Fj6O+OO+5o9yPx0EPve6cBiYm2b9q0yW53k7J7NGzY0I7lkksuMR07djS77babNW/pexonRBeiuQBCiib3e+65x5on3n//fTtJy3w0bNgwm1ulh0xDmrAT3WlKBL788ku7itEKRn/lx2jVqpV54IEHrFnFX/xUJpIDDzzQPPTQQ/YzWqXr83vssYdp2bKlOfnkk625asuWLdZ0hZ8D0iG2Zi58JpBvPv74Yzvxv/DCC+bQQw+1juxTTz3Vrhy8fPvtt2bFihXW5CWx0GR/wAEH2BXMaaedZs1PXrQviYOESDb+1q1bmzZt2tiHnu+zzz72c9qmRzJYZUBtiK2Y4DOBfKJVgRysMjFNmzbN+iG0slDZIIlAhw4dzMqVK63D+NNPP63yXZmTJCZqACdnsMRBqxf9bdeunWncuLH93BFHHGEfALkAnwlAjsNpFfEj85KiblTQ1DnaJSziqquuMtdcc401T1100UVm//33t499993XCgaVsyFf4DMByBNPP/20TeZTCKx7yE8ybtw4G+YpQenUqZOtO+dEQ6sS0bRp08owWoCwE1szF0BQyAmuMFoVL1X3TzmttbIQinh69913baLYLrvsYoVEwnHDDTfY7a+//joXAiJBKPuZ5AL6mUBd8Tq8lTimsFiZn5SMVlJSYl599dXK7f/85z9t4pxMWQq1VWXsGTNmcPIhcsR2ZYIDHtJddWhlIYFQtrceCuFV509FXSlCShFV3bt3t48uXbrYsFpvEp946aWXrCN94sSJhNpCJImtmAAkEw+V8FBorcxSt9xyi7n00kvttp/97Gc2OmrgwIHWkS4xkYM8HR5++GGz884729BfgCiCmEDsUV8bheQ+99xz5vnnn7f5H08++aStvXTSSSeZRo0amaOPPtom+tWlRIXqPJWVldn9aV8AUQQxgdihEFz5L7TSeOedd6xIuPyN4447zhbbcx0/lc+hRybMmzfPrFu3zvzmN78JZPwAYQQxgVggH8esWbPMM888Y/71r39Zc9Pdd99tQ3EnT55sK7HK35GN4nhTp061Yb4nnnhi4PsGCAtEc119db6vAWSZc845xzYE0spAYqLaVyo1LlR3SiXHu3btmhUhURTXY489ZoYPH16tTApAlIjtyoRormiirPGZM2da05JKlkggVJ5EjnIJiJIE69evn7Px/PWvf7VjOP/883N2TIB8EFsxgeigBkcyYd1///3mqaeesq9VekS9OPbaay9z5ZVX5mVcyi+RKU3RXwohBogyiAkUtCNdq4y5c+daH4jKtJ933nm2050c6Lkqna6mVY8//rj573//W/meeqHfe++91tGvLHiAqIOYQMFlnyuDvLS01GaTX3fddbb9qmpg6W8uyqbLlKZOg+oLom6DCitWmXgvEjKVileOiorjAUQdxAQKAmWhqze47vYVmaVwXvlChAREfbOzjYo0KmBDbWqV3KhVkcxpcuyrX7miwRyuHzpAXEBMIBSoPWyqDtJjxoyxkVhyoquEu5IIdfevnuXZQvv++9//bvuNKCteJiv1CtFYlMyoJlQIBsD/h34mkBe0upg/f741EcnnoeTBsKISKooCUwOq3/3udza5ESCO9OjRwyxevDjhttiuTGjbmx+HuUJl5e9YtWqVfU/1r4499ljbPlamIZVkV0mT9evX20m7f//+leasXKOVjzLiDzvssLwcH6CQiK2YkGeS+/pXCpF95ZVXrHioTa1MVbrrVzKf6lcpAku9zg855BArODJp5TInBADqTmzFBHKHCicWFxdbQZEPQkKiu341ipKJq0+fPnZVokx19QU55ZRTKNMOUGAgJpBVVKpdIbvKu5B/RCsR8fLLL9syJrK/6vlRRx2Vdjl3AAgfsa3NBbnh9ttvt2G906dPt0KiWlWDBg0yxxxzjPVbqce5ExgAKFxYmUDW2LBhg00qPOGEE2xxRflF5C/R+8rXUMvbhg0bcgUAIgBiAoGivh333HOPeeCBB8xHH31kM8MHDx5sI7nkF9FKRP1BWrduzZkHiBCICdQKiYIc6RIKOdaVyCfBeOGFF8xrr71WWZ/qF7/4hc3NkIlLjnWhLPHjjz+eMw4QQRATSAvVvrr11lttaG+irHMl9Kn5k8qcyKS1ZMkSa8aS0FxzzTU2jwQAYigmqj+0cOHCjEM0VXJCiWdQuKiwosSgXbt2ZtiwYaZjx46mbdu29qHWtzJfqfmU+7cycuRIW0dLeSRTpkwxHTp0yPdPAIB8iYnCOCdOnJjxAdRhDjEpTFQr689//rNdYSjiSu1uGzVqlPLzEhQJj/wif/zjH7PSvRAACkhMioqK7F1nplB+u3B9I/JxKMlQvULkUE8mJCrSePnll9vqvTfeeKNdkegBAPEh6W3j6NGjAzlAUPvJVm0uesAnXmFccMEFVkh0ftTDPJmQqECjaldpBbNly5aUlX8BILrE1gFPba7k3HLLLeaOO+6wpdbHjx+f9HMSmaFDh9raWioPLwc8AMSTQMREWc1Lly41TZs2tV3odMePn6QwUXjv2LFjzWmnnZbSZ7Z27VpbY6tz585m5syZttc6AMSXjMXkjTfesH9PP/30Ku8/8sgjtibTrrvumukhIEfITCWBaNOmjW0IlSiS78cffzQ77LCDadmypQ3SkE9MKxMAiDcZh9qsXr3adO3atdr7EhdVhIXCQSsMdRT8y1/+Ypo0aVJtuxIVJR5qnStUnBEhAYCMxWTNmjW2B4XjpptuMs8991zl6yCiwSA3yHE+efJks//++9uWtH4WLVpkM9qV9U4pFAAIVEzkH/G2MP3HP/5hTR+OTBMeIXcos13l4FUG3p8b8uyzz9qOg+qKuGDBAmu+BAAIzGeijnhaibh6S/7ewBIbKAxkumrcuLE588wzq7wvs1e/fv3MQQcdZGbNmmV9JQAAgftMkgmGIryU+AiFkaD45JNPmpNPPrlaSXhdQ1UBnj9/PkICAHUTk6+//trUhBztitxyvpKvvvrKvtZfrVwg/Lz66qvmyy+/tD3XHcozUZiwGDJkCP4vAKi7mEyYMMGkgwRFpq558+bZFYleJ4rwgnCi5ENFZbmkQyUtjho1ytx11135HhoAREFMFCr6/PPPp72z4uLinInIxo0bbXKdkiXTeR+Sm7gkJr169bI+EwnJJZdcYm8Ipk6dymkDgMzFRF3zNMnUr1/flpIfN26cNWclM3+lYxYLCjn75RxO931IzBVXXGFDvM866yxbLl5CMmDAAPPwww/bwo0AABmLie7uZepQP4tBgwbZhkcSF5VN2W+//cx5551n+55oMqqNWSwINI5mzZql/T5UZ/r06WbSpEn2OmolorDuvn37moceegghAYDgQoP32Wcf2wxJvhBF9cyePbtSZCQsmnzOPfdc62x3GdO1ERSZpHQ3vH79elNSUlJtuyY6HXfDhg2VvVEgGJYtW2bb6apUvExbygnSakQlVViRAEBW8kzkCxGK0mrfvr3NetdDQiPkdJewqJdFuqjUisRE5VgSIb+HTGsyubjX8uG411B31HZXhRy1wrzyyittv3b5TVq1amW7JgIAZDXPRKYQlUiRacvrH9EKRqsGf7HHVMgcJWFIVANKaMXiFY7Bgweb0tLS2gwXUpi3dAOgVaQaYKkCMB0RASATan0bKuHQQ6Yvmbe8peZHjBhhgiBRJJZEh8KRwaAcEtXgkhlRqxRdyz322COgvQNAHKmzTcNr+pK5RHkmEpkgkI/E70T3v5awKHLLrWxcwclk78P/Cjbq0aFDB7Nq1Srb112lUgAAMiEjA7lMXfKhaPKWCUpO9L333ttkinwpqbZJKGQmUxCAn2TvJ2vb61BHwTi08H3wwQfNTjvtZM1ad955pz1fAABZFROJhXI23EPOcu9rV7pck7uiriQoQYQHa38ugsvhf50pcW3bu2LFCnPwwQebF1980Ta5AgDIuphoUtfdu/7KhCXBUL0tOcb1XI+gTFt+k5Z/deJeJ3PYQ3pCojpcun4ICQDkTEwkFnKqO7NSNoQjEfJz+EVDK5MgTTLOzBUX85ZWmSrkuHnz5pxdRwCIDynFRHewo0ePts9ff/11c/fdd1cRGtfHxKGQYW90VyYo496bV6I8lqCixeJm5pIpUudObXfFz3/+83wPCQAiRr0KzTR1QHkKLvFQd/jKP1Go6XvvvZd2+K++73JHNNlp5eGNvtL+9Nr5Z4LMgFcvc38zr6gybdo08/vf/96ceuqp5oknnrDXLohACQCIFz1SzJt1FpNE4iAxCNpRni20Mvn8888jb+b64YcfbB01PTp27GgbXX3zzTckKQJAoGISWO0MrSAKqXZWXMxc6lMip7sqPw8dOtTml5DtDgA5K6fiKgHXhokTJwayHwiGBQsWmJ9++sm0adPG1t1auXKlXZ0AAORMTBJV8a0LKtAIuWfhwoW2IvBNN91kX8uX9dlnn5lOnTpxOQAgcJKaueRKUTOsTNA+AnLJBE6UQ4NVb0smrdatW1dGwOl37rLLLubMM8/M9/AAIE5iEtTKJKxE2Wdy1VVX2bpb6j+jKLvly5fbXiW6OaCgIwDkVEw0CUHhobpkkydPtiuS3r172/duvvlm29/d5QwBAOS1nwmEH2W4H3nkkZUry02bNtnKzuoHo+rOAADZILZi4nwmUfOXHHXUUeall16qXFlKSL799lvrQwEAyBax7dEaNZ/Jp59+aqsJyC8iR7tQ8IOy39UmQCIDAJAtYrsyiRoXX3yx9Y2oBa9QmZszzjjDlpr/wx/+UKV3CwBA0CAmEWD+/Pm2KOYVV1xho7UuueQSGxasXu9KJL300kvzPUQAiDixNXNFhW3btpkLL7zQtGzZ0rz77rumS5cutpCj8kkuuugi07Vr13wPEQBiwHbZ7J8RZqLigJdPZNmyZbag45NPPmlF5bnnnjP33nsvQgIAhSUmrmSHl2HDhpmwO+DloC50MVFvkpNOOsmUl5ebBx54wLzyyivmuOOOy/ewACBmBCImN9xwQ2VBxzfeeMOWKVavEsgucrK/+eabtubWvvvua/r27cspB4DC9ZnIRq8a93L26g553rx5ZNBnGXW+VKdLCYq44447KC0PAIUnJl6fiHwPhx56qO2IqBIeSpiT6euyyy4LapzgYevWrfY8S0iUW9KnTx/Trl07zhEAFJ6YqO2rBMRVBZag6LlCVPVXNaIQk+xw7bXXmvXr15sTTzyxoBqSAUB0qbOYlJWVmeLi4qTbZeqC4Pniiy+sOVHdEhXJBQBQ0A74VEIiJ/xXX31lwkwhhgara2L//v3Nli1bzLnnnmvDgAEAIuWAV8lzOd+d2Us2fU18YaUQa3NpJaLQ3549eyZskQwAUNBiIiewHqJZs2b2r5zxECyPPfaY2X///a0JkVpbABC5PBP1ythnn33sw4UFa5UCwUZw6dx26NABIQGAaK5MNmzYYCO7Fi1aZBYuXGj/KmxVeRAQDE888YQtmSK/CQBA2KhX4ZwcAaKGTN26dbMrlbCiLH0lWhYKEmuNV/3cO3funO/hAEAM6ZFi3twuqGxsJc5pwhNFRUWhj+YqpDa8KimvC9ikSROEBACia+ZSguLs2bOtqIhDDjnEVq6FzLnmmmvMLbfcYp8rYg4AIIwEsjLRSkR4I4xczaiwUgh5Ji+//LK58cYbzeGHH25GjRpFODAARHtlsnr1apsRLzOMQoLlgFff8TAT9jyTOXPmmNNOO82WrZk1a5Y9twAAkV6ZqHeJakTJly9zV69evWyGNtQN1TdTjxKt+NSKt1GjRpxKAIhH2175SZyZi1axdUeCfPbZZ9vzOXLkSHPWWWeZ3XbbzfTr1y+oSwUAEE4xkeN94MCBpmnTpvahhEWZvWSigdohX9OmTZtsIujf/vY3ew5pegUAsRCT6dOnm/fff7/Ke/QzqXtVYPHjjz+a+fPnm5KSElO/fv0ArhIAQAGUU/EjMw3UXUwWLFhgdtxxR3POOedwGgEgHmKiqsHqAa/ui3qoBL3ec6/HjRsXxGFiJSZa6Q0aNMi0aNEi30MCAMjMzCUh2HXXXWvciSK3vF0XHfKbuK6LEyZMqHk0UCkmiopr3LgxZwQACl9MJADpiABdF4Nj7dq19u/uu+9uGjRoEOCeAQDyZOZSvsPzzz+fUdfFdLbD/3jrrbdsiLWqLwMAREJM1q1bZxMQFU0kM5Z8H6q5JfNXIpK9H0bCWk5l6dKl9m+nTp3yPRQAgGDMXJrY5s6da7snqkyKSnwoVFWTsLKzJTRqz6ty88qHSNcsFgbCWE7l+++/Nx9//LFNUlSDMQCASKxM1I9EpVJUF0riIaewmjOp9taYMWPM+vXrrfNddbiaN29upkyZkruRR7Qt77Zt22w3RQCAyCUtOp+Hml5JOLQS0UNCIxQGrFWLKtxC3Zk2bZrZbrvtTMeOHTmNABDdPJPTTz/dml8effTRKv4RrWBU6FHboe6MHz/ervxatmzJaQSAaJdTkXDoMW/ePNtNsX///pXbaN6UGW3btq0MCwYAiEVtLq/pS8Udjz/++FD3fA8zSuz89a9/bWbMmGFfszIBgFiVU5GpSz4URXqpPpdKqkDtWbZsmRWSXXbZxdx6662mT58+nEYAiFY5FQmFe6ijove1u6t20V6lpaUFExocJh5++GH799hjjzUXXnhhvocDABCsmEgklFOivzJhSTBUDXjAgAH2uR6YtjJDDvcHH3zQPu/Zs2eGewMACKGYSCzkVJeYKEER4QieV1991Xz22Wf2+ZFHHpmFIwAA5FlMtAIZPXp0ZTfFu+++u4rQyOnuRSHD3uguqBllvKvH++bNm0337t05ZQBQkNSr8NeNTxMlKqrUitrMyhSm/JNJkyaZ9957z+QCHVf+GTn+lUDpUBa+hE4+Ha2m9DwRPXr0MIsXLzb55ocffrArPtXiUoUBAICwkmrerHNosCu14q3jpfIquUI/yAUBOFyQgBIohfrSqzx+WNFqRP6Szz//3Nx77735Hg4AQH47LQqtDtwkngu06lABSn/JfIUq+yvwhpXLL7/cnH/++aZLly62YCYAQKFS55VJIiZOnFhrU5XMUlrRqBqxH5nNZKbasGGDfV2TWGk/ycxaYUPWRa2atDoZMmSINRUCAJi4r0xqi/wtesgsJVHxM3bsWCsMCgKQiOhzWnlEBbUydiXw+/btm+/hAAAUppjITCWhUNhxIrRi0XaHHO1KikyFyuC7VUzYUeSbViNt2rSh5DwAFDx5E5NUJPJ1SHS0kkmFxEcrGIc3yitsJi6tsiQm/fr1w8QFAAVPKMVEqwu/c93/WsKiiK7p06dXio/MYmovrG1a2STyw/jb9rpHLtv3SkyGDh1qs99dwUwAgEImUAd8UCTyoXi3uYx8+R38eE1jYW3bqwZY22+/fWU9LgCAQieUKxOJhd/3USi+kHR46KGHzDPPPGN9Jer3DgBQ6IRSTGTS8q9O3OtkDvva4sxcuTRvCTUUk4nrtddeY1UCAJEhlGYuOc79oqGViUxbQZEvM5dKpmzdutU+P/roo3N+fACA2KxMxKBBg6rklcyZMycSbYGffvpp07BhQ/u8a9eu+R4OAEB+Cz1miiKwFHXlckckFFp5eMN5lQGv164GV5DlWvJR6FHRW61atTItWrSwBTG//fZb06BBg5yOAQAgG/Nm3sQk38jMpQKL48ePz5nfZNWqVaZz5862QvC2bdvM8uXLc3JcAIDQVg0udPLhMznwwANt/bCOHTuaY445JqfHBgCIpc8kqmgh+Omnn9oVCgBAVIitmOQ6NPibb74xxx13nLnvvvvsa8QEAKIEZq4cMX/+fPs45JBD7GvEBACiRGxXJrlGoc0777yz+f77703jxo1N27Zt8z0kAIDAQExyhMKg5XR/++23bTQXzbAAIEogJjnyz6xcudLm0axYsQITFwBEjtiKSS4d8OXl5eaEE06wGe96jr8EAKIGDvgcoLySWbNmmWeffda+RkwAIGrEdmWS60rB4vXXX7d/ERMAiBqISZZZu3atadq0qZk2bZr597//bfbbb79qXSMBAAodxCTLvPjiizbrXaYu9TA5/PDDs31IAICcE1sxyZUDXomKjRo1sh0VtUpBTAAgiuCAzzISkyOPPNJcfvnl9jViAgBRJLYrk1ywbt0689Zbb5mioiJTVlZm3zv44IPzPSwAgMCJ7cokF2y//fbm9ttvN999911lQ7Addtgh38MCAAgcViZZRH3sR44caRthCUVyAQBEEcQki8yePdt88skn9iFhkSMeACCKxFZMsh3NtWXLFnPqqaeayZMnWzHZa6+9snIcAIAwEFufSbbb9i5btsxs3rzZHHHEETaii5LzABBlYrsyyTYLFiywf7t162bef/99065du3wPCQAgayAmWULZ7lr9LF682LbslckLACCqICZZXJnIxKWmWM2bNzfFxcXZOhQAQN5BTLKE/CQ33HCD2bp1q43iql+/frYOBQCQd2LrgM82RG8BQJxgZZIFHn30UXPbbbfZ56oYDAAQdWIrJtnMM7nnnnvMXXfdVflaxwEAiDKxNXNlM89ENbj69OmTlX0DAISR2K5MsoUESn1LlF8CABAXEJMsrEpE9+7d7V98JgAQBxCTgPnoo49MgwYNTNeuXSvfw2cCAFEHMQmYUaNGma+//poKwQAQKxCTLLDTTjtVPsfMBQBxADEJkC+++ML07t3bvPzyy1Xex8wFAFEHMQnY+a5aXK6zIgBAXEBMAmTFihX2b5cuXYLcLQBA6EFMAhaT1q1bm6ZNm1a+h88EAOJAbMUkG+VU3nzzTdOpU6dq7+MzAYCoQzmVgNAKpFWrVraHCQBA3IitmASNVh9PP/10tfcxcwFAHIitmSuXYOYCgKiDmATE9ddfb0uoEBYMAHEEMQkwx2Tz5s205wWAWIKYBBgW3Llz52rv4zMBgDiAmATAd999Z1avXp1QTAQ+EwCIOohJAKxatcquQA466KAgdgcAUHAgJgFVCR4yZIg5+OCDq23DzAUAcYA8kwDQiuT+++9Puh0zFwBEHVYmAaBmWKxAACDORGplMmnSJFNUVGQ++OADM3z4cNOkSZOcHLdnz55mr732Mo8//nhOjgcAEDYiszJRHxFFVA0YMMAKydixY3NyXK1I3n33XdOuXbuk2wEAok5kxGTOnDmmffv29rlWJBKXXFUf3rRpkznggAOSfgafCQBEndCZuTZu3GimTJli1q9fb0pKSpKasjZs2GBfaxUiJCRambh9yNSVC9555x37N5WYAABEnVCtTLSacOYqCYIfma4kJM6Upc/NnDnTbtNrCYi+t3jx4sD8JTX1O6lJTMJu5gqyn0vYiPJvE/y+wubqiP37rFcRwtlOoiFRKC0trfK+OhiWl5dXqYelz8rE5X1PglNcXGyWLFmS9Bg9evSwolMTMlGlOkULFy40Tz31lLn22msTmrMGDx5sli9fbt5++20TRmr6fYVMlH+b4PcVNvUK8N9nqnkzVCuTVEgk/Hh9I1qVjBgxwnTr1s2+N27cuJyM67DDDjPXXXcdfhEAiDUFIybykTRr1qzKe97XWo307t3bCok+K1NYLlCrXuWZAADEmdA54JORyIfi3aZVSm0EZOXKlVVWE2q5u+eee1b7nN7X0i4TGjZsmPE+skUQvy+sRPm3CX5fYdOqAP99rlmzpvDFRGLhIrgc/te1rfQLAAAxM3PJpOVfnbjXucp0BwCAAhcTOdb9oqGVSa9evfI2JgAAKDAxEYMGDarMKxEKCVYEFwAA5JdQ5Zko/FfRWC6/REKhlYdWJd4MeL12Ge4uA76uJMuoD+rzYagmIBYtWmSj3VKNV59Vbs7AgQPt67KysspE0bBRl7EW0rXT79K/f43XH8WYyKwb9mtX18oWhXItN6b4fVH+f1iFihgzZsyYirKysqSvM/18vtH4vBQVFVWUlpYm/by2NWnSRDcXFd26datYsmRJRVip7VgL7drpWum3+R/Jxhzmazdnzhw77uHDh9tH1P4fpvP7ovr/0EusxUQXzIsuWq9evQL7fD4pLy+vGDBgQJX3SkpK7D/kZKT6Bx42ajvWQrp27lql814hXTtNqokm26j8PxyT4PdF/f+hl4LymeQyoz7Tz4cBjc1b8FLjzVUBzDBRiNfObwaRWWfMmDEmavD/MDoUTJ5JrjPqM/18vtFk6a1j5gIWaop+k71Wvysstuggxlpo187vF5HopRO1WEjXzsH/w+hcy9iKSToZ9Zl8PmxojJqU5s2bl/QzysbV73COPjkA9Q86V6VpakNtxlro104OWH/R00K+dl74fxidaxlbM1dtM+qDzsDPNcOGDbOTkjcyzo+2eSNGDj30UDNhwgQTRmoz1kK+dgqF7969e42fK6Rr54X/h9G5lrEVk9pm1BdyBr7s7S7MOhV+H4L+QSeyaYeB2oy1kK+dViTphIQW0rXzwv/D6FzL2IpJbTPqCzUDX3e2GrsbZzKnsxzzin/3T7phjG2v7VgL9dq561WTf6eQrp0f/h9G51rGVkzSyajXhfVuL7QMfFeOXzZY18rYe4fj/X36x6pkK++kO336dJssFTbSGWuhX7uaVk+Feu0Swf/DDyJxLUOVAZ8PUmXUa5smHW8nx6Az8LM5EakzpR858eQ7SfT7vBOUMnnbt28f2t9X01gL+dp5r6H8JcqG9gtKIV27TCtbhP1apvp9GyP+/9BL7MUEAAAyJ9ZmLgAACAbEBAAAMgYxAQCAjEFMAAAgYxATAADIGMQEIOYofFXhrHrkuqq0q87gzQGCwiS2hR4B4H8ojyVZB0TlOjRv3tx+xhUc1OSvXIpMy9GorL7ETLWnwl7IEFKDmADUgLKPNeHVVLk3HRJVKQ5y/0Gh8bjWwd4+Kq4FrYRHyZQADsQEoAYGDx4c2DmaMWNGtWzmIPcfFMXFxVYw/PXLJIQqfxL2UjSQe/CZANSAymKkKt1fG7wlQbKx/yBwjZmSFcKUoGCSAj+ICUCOcOassCNzm0xcqQjjagryC2YuCC3uDlmo+rHXPFTTNleyW5O3q5ysyVx31ePGjbPvaduiRYusOceV5lfRPn3X3XkrusmZdLyrilTHT4Qc1q5ys5zaQt/Rd7371/HrMk5/AUT3u+tSINDtOxXe4zofisas4oXuuyqlrjFMnTq1skGbroOESudBv9cVO4QIoKrBAGFjzJgxFaWlpZWvV69eXfl6wIAB9rWjvLy8olevXva5PuPfpn2JOXPmVHTr1q3Kdu3LbXc0adKkyuslS5ZU7r+msaXCv59k79d1nPqM9uUdZ1lZWY3j8p4jHVPTgnc/6aDfP3z48Crv6djad6LfpfdLSkqqHR8KF8xcEDpcxJD3rlp39qtXr7Z3zbqr9d45u37Z+o7wRkVpmzPJuK5+3u8mugP3d//zRl+lGltd8Ud31WWcrleN1/ei313bCLG6NmHS+VBwQarf5v1det8bJQaFD2ICoUOmHP+kpolHppTFixcnnPDU80GhqprUNLHWq1fPmlk08Xsn2ER5EcqhCGJsQVLbcWpc+o7+uofOQ12SEPX7dJ79SAhckqHOrz/JUQ57l3yYKATa7RuiCWICBUUqB7bzL8gOX15ebn0Peu4NY03UAjdffeBTTfS1Hae749eE7h7ya9RlxeRt3OQ/vlc4dX694iAfj1sJScwSRYPl61xD9kFMIHR4O+j5J0xNUIm2adLUSkSZ1G7S0mfl5A2yREiqsdUWbwvlbI2rLkgs/C2e0xE8N4Zcl2SBcICYQOjQ3a4S45wPxKG7XU1Y2u6d6DSRyyzj/Bj+73nvnrV6yUQIUo0tne96J9pUd+m1HaeEU5O7XwD840wXibBWHolqZsk3ksxc5Wp8JcubKYTQaKgj+Y4AAEiGon0UJaQoIH9UkrbpPT303EUN6bk+7x7arughRScp2kkRUC6KSNuKiopshJE+676v/xaKTNL39HDf80ZspRpbKtz33Hf8+6/rOBPtP91xpYqm0v60zZ1vFy2ncXojtbz7Shax5v1d3u8SzRUN6AEPEHNcocWgggi0mqlNhnzQx4f8gJkLAAJD5j5KrcQTxAQAMkI+EuerwScSXxATAMgIlUdRwEBtzVsQLajNBQA2ykyiIL9FbRMLk1UXTgclQaruGIUjCx8c8AAAkDGYuQAAIGMQEwAAyBjEBAAAMgYxAQCAjEFMAAAgYxATAAAwmfL/AC8GZCuQHkDLAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "icheck = np.random.randint(0, msk_has_fit.sum())\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "yscale = ax.set_yscale('log')\n",
    "__=ax.plot(bpl_mah_data['bpl_t_table'], bpl_mah_data['mpeak_history'][msk_has_fit][icheck, :], color='k')\n",
    "__=ax.plot(bpl_mah_data['bpl_t_table'], 10**log_mah_pred[msk_has_fit][icheck, :], '--', color='k')\n",
    "\n",
    "ymin, ymax = ax.get_ylim()\n",
    "ylim = ax.set_ylim(10**9, ymax)\n",
    "\n",
    "xlabel = ax.set_xlabel(r'${\\rm cosmic\\ time\\ [Gyr]}$')\n",
    "ylabel = ax.set_ylabel(r'$M_{\\rm peak}\\ [M_{\\odot}]$')\n",
    "\n",
    "halo_id = bpl_mah_data['halo_id'][msk_has_fit][icheck]\n",
    "title = ax.set_title(rf'${{\\rm halo\\_id={halo_id}}}$')\n",
    "fig.savefig(f'example_diffmah_fit_bpl_halo_id_{halo_id}.png', bbox_extra_artists=[xlabel, ylabel], bbox_inches='tight', dpi=200) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "729cb1f0-35b1-47e6-a284-49cf0910fc9b",
   "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.13.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}