{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5f0fc458",
   "metadata": {},
   "source": [
    "# Multi stand level modelling in libcbm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9f64044e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from libcbm.model.cbm import cbm_variables\n",
    "from libcbm.model.cbm import cbm_simulator\n",
    "from libcbm.model.cbm.cbm_output import CBMOutput\n",
    "from libcbm.model.cbm.stand_cbm_factory import StandCBMFactory\n",
    "from libcbm.storage import dataframe\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc413246",
   "metadata": {},
   "source": [
    "define the CBM classifiers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eabb4e64",
   "metadata": {},
   "outputs": [],
   "source": [
    "classifiers = {\n",
    "    \"c1\": [\"c1_v1\"],\n",
    "    \"c2\": [\"c2_v1\"],\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d1e057f",
   "metadata": {},
   "source": [
    "define one or more merchantable volumes.  The relationship between stands and volumes is defined by the classifier sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6e2fefec",
   "metadata": {},
   "outputs": [],
   "source": [
    "merch_volumes = [\n",
    "    {\n",
    "        \"classifier_set\": [\"c1_v1\", \"?\"],\n",
    "        \"merch_volumes\": [\n",
    "            {\n",
    "                \"species\": \"Spruce\",\n",
    "                \"age_volume_pairs\": [\n",
    "                    [0, 0],\n",
    "                    [50, 100],\n",
    "                    [100, 150],\n",
    "                    [150, 200],\n",
    "                ],\n",
    "            }\n",
    "        ],\n",
    "    }\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd4524c8",
   "metadata": {},
   "source": [
    "define the inventory to simulate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4c0f5c82",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_steps = 50\n",
    "n_stands = 1000\n",
    "inventory = dataframe.from_pandas(\n",
    "    pd.DataFrame(\n",
    "        index=list(range(0, n_stands)),\n",
    "        columns=[\n",
    "            \"c1\",\n",
    "            \"c2\",\n",
    "            \"admin_boundary\",\n",
    "            \"eco_boundary\",\n",
    "            \"age\",\n",
    "            \"area\",\n",
    "            \"delay\",\n",
    "            \"land_class\",\n",
    "            \"afforestation_pre_type\",\n",
    "            \"historic_disturbance_type\",\n",
    "            \"last_pass_disturbance_type\",\n",
    "        ],\n",
    "        data=[\n",
    "            [\n",
    "                \"c1_v1\",\n",
    "                \"c2_v1\",\n",
    "                \"Ontario\",\n",
    "                \"Mixedwood Plains\",\n",
    "                15,\n",
    "                1.0,\n",
    "                0,\n",
    "                \"UNFCCC_FL_R_FL\",\n",
    "                \"None\",\n",
    "                \"Wildfire\",\n",
    "                \"Wildfire\",\n",
    "            ]\n",
    "        ],\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cfd0d33",
   "metadata": {},
   "source": [
    "simulate using `StandCBMFactory`\n",
    "\n",
    "note pre-dynamics func argument may be used to inspect and or modify the CBM state, variables and parameters in the simulation loop.  In this example it is used to assign a disturbance type to all stands on timestep = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4edf04ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pre_dynamics(t, cbm_vars):\n",
    "    if t == 10:\n",
    "        parameters = cbm_vars.parameters.to_pandas()\n",
    "        # note type changes will not be allowed in libcbm, so care must be used when using pandas\n",
    "        parameters.loc[:, \"disturbance_type\"] = np.int32(1)\n",
    "        cbm_vars.parameters = dataframe.from_pandas(parameters)\n",
    "    else:\n",
    "        parameters = cbm_vars.parameters.to_pandas()\n",
    "        parameters.loc[:, \"disturbance_type\"] = np.int32(0)\n",
    "        cbm_vars.parameters = dataframe.from_pandas(parameters)\n",
    "    return cbm_vars\n",
    "\n",
    "\n",
    "n_stands = inventory.n_rows\n",
    "\n",
    "cbm_factory = StandCBMFactory(classifiers, merch_volumes)\n",
    "csets, inv = cbm_factory.prepare_inventory(inventory)\n",
    "\n",
    "with cbm_factory.initialize_cbm() as cbm:\n",
    "    cbm_output = CBMOutput(\n",
    "        classifier_map=cbm_factory.classifier_value_names,\n",
    "        disturbance_type_map=cbm_factory.disturbance_types,\n",
    "    )\n",
    "    cbm_simulator.simulate(\n",
    "        cbm,\n",
    "        n_steps=n_steps,\n",
    "        classifiers=csets,\n",
    "        inventory=inv,\n",
    "        pre_dynamics_func=pre_dynamics,\n",
    "        reporting_func=cbm_output.append_simulation_result,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff816825",
   "metadata": {},
   "source": [
    "plot the mean age to confirm the disturbance event had an effect on age"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4233a0fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:xlabel='timestep'>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEGCAYAAAB8Ys7jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArpUlEQVR4nO3dd3xUZfbH8c8JaUBCD0gLoUgXAoQisIrYsKyKvQC6KgGFtfey6s91V1dX7AquLaGpIGLvKIqCJiH03kNLaCEQUiZzfn/MoFEpIZnJ3Jk579eLV2ZuZnLPDZNvnjxz73lEVTHGGBN8IgJdgDHGmMqxADfGmCBlAW6MMUHKAtwYY4KUBbgxxgSpyOrcWaNGjTQpKak6d2mMMUEvMzNzh6om/HF7tQZ4UlISGRkZ1blLY4wJeiKy4VDbbQrFGGOClAW4McYEKQtwY4wJUtU6B34opaWl5OTkUFRUFOhSfCo2NpYWLVoQFRUV6FKMMSEq4AGek5NDfHw8SUlJiEigy/EJVWXnzp3k5OTQunXrQJdjjAlRAZ9CKSoqomHDhiET3gAiQsOGDUPurwpjjLMEPMCBkArvg0LxmIwxzuKIADfGmFBVUFTKwx8sYW9Rqc+/tgW4Mcb4ybod+xn60o+kz93AL+t2+fzrB/xNTGOMCUXfrczj75OzqBEhTLyuLye2bejzfVR4BC4iNURkvoh85L3fWkTmichqEXlbRKJ9Xl01uuCCC+jVqxddunRhwoQJALz22mu0b9+ePn36MHLkSMaOHQtAXl4eF110Eb1796Z3797MmTMnkKUbYxxEVRn/3Rr+9sbPNKtXkw/GDvRLeMOxjcBvBpYBdbz3nwDGqepUEXkFuA54uSrFPPLhEpZu2VuVL/EnnZvV4aG/djnq415//XUaNGjAgQMH6N27N+eccw6PPvooWVlZxMfHM3jwYLp37w7AzTffzK233srAgQPZuHEjZ555JsuWLfNp3caY4FNUWsbd0xcyM3sL55zQlCcv6UataP9NdFToK4tIC+Ac4DHgNvGcYjEYuNL7kLeAh6ligAfSc889x4wZMwDYtGkT6enpnHzyyTRo0ACASy65hJUrVwLw1VdfsXTp0l+fu3fvXvbt20dcXFz1F26McYQtew4wKj2TxVvyufPMDtw4qK3fz0ar6K+GZ4C7gHjv/YbAHlV1ee/nAM0P9UQRSQVSARITE4+4k4qMlP3h22+/5auvvuKnn36iVq1aDBo0iI4dOx52VO12u5k7dy6xsbHVXKkxxol+Wb+LGyZmUlTq5tXhKZzWuUm17Peoc+Aici6Qq6qZldmBqk5Q1RRVTUlI+FM7W0fIz8+nfv361KpVi+XLlzN37lz279/Pd999x+7du3G5XEyfPv3Xx59xxhk8//zzv97Pzs4OQNXGGCeY8vNGrnx1LnExkcy4sX+1hTdU7E3MAcB5IrIemIpn6uRZoJ6IHBzBtwA2+6XCajBkyBBcLhedOnXinnvuoV+/fjRv3pz77ruPPn36MGDAAJKSkqhbty7gmW7JyMigW7dudO7cmVdeeSXAR2CMqW6lZW4efH8x9763iBPbNmLmmIEc3yT+6E/0oaNOoajqvcC9ACIyCLhDVa8SkXeBi/GE+tXATP+V6V8xMTF8+umnf9qekpJCamoqLpeLoUOHcsEFFwDQqFEj3n777Wqu0hjjFDv2FXPjpCx+XreLUSe14a4hHakRUf1XX1flQp678byhuRrPnPhrvinJOR5++GGSk5Pp2rUrrVu3/jXAjTHha8mWfM5/YQ4LNu3hmcuSuffsTgEJbzjGC3lU9VvgW+/ttUAf35fkHE899VSgSzDGOMiHC7Zw57QF1K8VzbujT6Rbi3oBrccRV2Kqasg1f1LVQJdgjPGRMrfy3y9W8NK3a0hpVZ+XhvWkcXzgz0ILeIDHxsayc+fOkGope7AfuJ1maEzw21tUyi1Ts/lmeS5X9GnJI+d1JTrSGW2kAh7gLVq0ICcnh7y8vECX4lMHV+QxxgSvNXn7GJmWwcadhTx6QVeG9U101EAz4AEeFRVlq9YYYxxn1opcbpoyn6gaEUy6vi992/inn0lVBDzAjTHGSVSV8bPX8sRny+l0XB0mjOhFi/q1Al3WIVmAG2OM14ESTzOqDxZs4dxuTXny4u7UjK4R6LIOywLcGGOAzXsOkJqWwdKte7lrSAduONn/zaiqygLcGBP2fl7naUZV4nLz2tUpDO5Yff1MqsIC3BgT1ibN28BDM5eQ2KAWE0ak0K5x8LSFtgA3xoSlEpebRz5cwqR5GxnUIYFnL+9B3ZpRgS7rmFiAG2PCzo59xdwwMZNf1u/mhkFtueOMDgHrZ1IVFuDGmLCyeHM+qWkZ7Cos4dnLkzk/+ZBr0QQFC3BjTNiYmb2Zu6YtpGHtaKaN7k/X5nUDXVKVWIAbY0JemVv5z+fLGf/dWvokNeClYT1pFBcT6LKqzALcGBPS8g+UctOU+Xy3Mo+r+iby0F+7OKYZVVVZgBtjQtbqXE8zqk27CnlsaFeu6tsq0CX51FEDXERigdlAjPfx01T1IRF5EzgZyPc+9BpVzfZTncYYc0y+XradW6ZmExMVwZTUfvROahDoknyuIiPwYmCwqu4TkSjgBxE5uIDknao6zX/lGWPMsVFVXvp2DU99sYIuzeowYXgKzerVDHRZflGRRY0V2Oe9G+X9Z8vNGGMcp7DExZ3TFvLxwq2cn9yMxy/s5uhmVFVVoZl8EakhItlALvClqs7zfuoxEVkoIuNE5JBv6YpIqohkiEhGqC3aYIxxjk27Crno5Z/4ZNFW7jmrI89clhzS4Q0VDHBVLVPVZKAF0EdEugL3Ah2B3kADPKvUH+q5E1Q1RVVTEhISfFO1McaUM3ftTs5/cQ45uwt545rejA6CToK+cEzn0qjqHmAWMERVt6pHMfAGIb5CvTHGeVSV9J/WM+x/86hfK4qZYwYwqEPjQJdVbY4a4CKSICL1vLdrAqcDy0WkqXebABcAi/1XpjHG/F6Jy819Mxbx4MwlnNQ+gRljBtAmIXg6CfpCRc5CaQq8JSI18AT+O6r6kYh8IyIJgADZwGj/lWmMMb/JLSjixolZZGzYzZhT2nLb6cHZjKqqKnIWykKgxyG2D/ZLRcYYcwQLc/aQmpZJ/oFSXriyB+d2axbokgLGrsQ0xgSNGfNzuGf6IhrFxTDthhPp0iy4m1FVlQW4McbxytzKE58tZ8LstfRr04AXr+xJwxBoRlVVFuDGGEfLLyxl7JQsvl+1g6tPbMUD53YmqkZoNKOqKgtwY4xjrdpewMi0DDbvOcDjF57A5X0SA12So1iAG2Mc6Ysl27j17WxqRkcyNbUfvVqFXjOqqrIAN8Y4itutvDBrNU9/uZJuLeoyfngvmtYNzWZUVWUBboxxjP3FLu54dwGfLt7G0B7N+feFJxAbFdr9TKrCAtwY4wibdhUyMi2DldsLeOCcTlw3sHVY9DOpCgtwY0zA/bh6B2MmZ1HmVt78Wx9Oam+N7yrCAtwYEzCqyls/rufRj5fRplFtXh2RQlKj2oEuK2hYgBtjAqLYVcaD7y/mnYwcTu/chHGXJRMXY5F0LOy7ZYypdrl7ixg1MZP5G/dw06nHc8upxxMRhs2oqsoC3BhTrbI37WFUegYFRS5evqonZ53QNNAlBS0LcGNMtZmemcO9MxbROD6G6Tf0p1PTOoEuKahZgBtj/M5V5uZfnyzn9TnrOLFNQ168qicNakcHuqygd9QAF5FYYDYQ4338NFV9SERaA1OBhkAmMFxVS/xZrDEm+OzeX8LYKVnMWb2Ta/oncf85nawZlY9U5LtYDAxW1e5AMjBERPoBTwDjVLUdsBu4zm9VGmOC0optBZz/4hx+Wbeb/1zcjYfP62Lh7UNH/U56Fy7e570b5f2nwGBgmnf7W3jWxTTGGAA+W7yNoS/N4UBpGVNH9ePSlJaBLinkVGgO3LseZibQDngRWAPsUVWX9yE5QHO/VGiMCSput/LcN6t45qtVdG9Zj/HDenFc3dhAlxWSKhTgqloGJHtXp58BdKzoDkQkFUgFSEy0Xr7GhLJ9xS5ufyebz5ds56KeLXhsaFdrRuVHx3QWiqruEZFZwIlAPRGJ9I7CWwCbD/OcCcAEgJSUFK1ivcYYh9q409OMalVuAQ+e25lrByRZMyo/O+ocuIgkeEfeiEhN4HRgGTALuNj7sKuBmX6q0RjjcD+s2sF5L/7Atr1FpF3b1zoJVpOKjMCbAm9558EjgHdU9SMRWQpMFZF/AvOB1/xYpzHGgVSV1+es57GPl9KucRyvjkihVUNrRlVdjhrgqroQ6HGI7WuBPv4oyhjjfEWlZdw/YzHTs3I4s0sT/nupNaOqbvbdNsYcs+17i0hNz2TBpj3cfOrx3GzNqALCAtwYc0yyNu5mdHom+4pdvDKsF0O6HhfoksKWBbgxpsLeydjEAzMWc1zdWNKu60PH46wZVSBZgBtjjspV5uaxT5bxxpz1DGzXiBeu7EG9WtaMKtAswI0xR7R7fwljJmfx45qdXDugNfed3ZFI62fiCBbgxpjDWr5tLyPTMti+t5inLunOxb1aBLokU44FuDHmkD5dtJXb311AXEwkb6f2o0di/UCXZP7AAtwY8ztut/LMVyt57pvV9EisxyvDetGkjjWjciILcGPMrwqKSrn17QV8tWw7l/RqwT+HdiUm0ppROZUFuDEGgPU79jMyLYO1O/bz8F87c3V/a0bldBbgxhhmr8xj7OQsIiKE9Gv70L9do0CXZCrAAtyYMKaqvPbDOv71yTLaN4nn1REptGxQK9BlmQqyADcmTBWVlnHve4uYMX8zZ3U9jqcu6U5ta0YVVOx/y5gwtDX/AKPSM1mYk8/tp7dn7OB2Nt8dhCzAjQkzmRt2MSo9iwMlLiYM78UZXawZVbCyADcmjEz9eSMPzlxM83o1mTyyL+2bxAe6JFMFFVlSraWIzBKRpSKyRERu9m5/WEQ2i0i299/Z/i/XGFMZpWVu/jFzMfe8t4h+bRoyc8xAC+8QUJERuAu4XVWzRCQeyBSRL72fG6eqT/mvPGNMVe3cV8yYyVnMXbuL1JPacNeZHawZVYioyJJqW4Gt3tsFIrIMaO7vwowxVbd0i6cZVd6+YsZd1p2hPawZVSg5pl/DIpKEZ33Med5NY0VkoYi8LiKH7HQjIqkikiEiGXl5eVWr1hhTYR8v3MpFL/9ImVt5d9SJFt4hqMIBLiJxwHTgFlXdC7wMtAWS8YzQ/3uo56nqBFVNUdWUhISEqldsjDkit1t58vPljJmcRedmdfjg7wPo3rJeoMsyflChs1BEJApPeE9S1fcAVHV7uc+/CnzklwqNMRW2t6iUW6dm8/XyXC7v3ZJHzu9izahC2FEDXDxn978GLFPVp8ttb+qdHwcYCiz2T4nGmIpYm7ePkWkZbNhZyKPnd2FYv1Z2cU6Iq8gIfAAwHFgkItnebfcBV4hIMqDAemCUH+oLebNW5NKgVrT9iWuq5NsVufx9ynyiakSQfl1fTmzbMNAlmWpQkbNQfgAO9Wv8E9+XEz7K3Mp/Pl/O+O/WMqhDAm/+rU+gSzJBSFUZP3stT3y2nI7H1eHVEb1oUd+aUYULuxIzAPILS7lp6ny+W5lHTGQEBUWuQJdkgtCBkjLunr6QDxZs4ZxuTXny4m7UirYf6XBi/9vVbHVuASPTMtm0q5DHhnZl9so81u8oDHRZJshs3nOAUekZLNmylzvP7MCNg9rafHcYsgCvRl8v287NU7OJjYpgSmo/eic1IGvDHvYV2wjcVNwv63dxw8RMikrd/G9ECqd2ahLokkyAWIBXA1XlpW/X8NQXK+jSrA4ThqfQrF5NAOJjIy3ATYVNnreRhz5YTIv6tZia2ot2ja2fSTizAPezwhIXd767kI8XbeX85GY8fmE3akb/dl5uXIwnwFXV/gQ2h1XicvPIh0uYNG8jJ7dP4LkrelC3ZlSgyzIBZgHuR5t2FZKansnybXu596yOpJ7U5k8hXTsmkjK3UlTq/l2wG3PQjn3F3Dgxi5/X72L0yW2588wO1IiwX/bGAtxvflqzkxsnZeJyK29c05tBHRof8nFxsZ7/goLiUgtw8yeLN+eTmpbBzv0lPHt5MucnWx858xsLcB9TVdLnbuCRD5eS1LAWr45IoU1C3GEfH+9dg3B/cRnYdKYp54MFW7hr2gIa1Ipm2uj+nNCibqBLMg5jAe5Dxa4yHpq5hKm/bOK0To0Zd1ky8bFHnqc8uIjsPjsX3HiVuZWnvljBy9+uoXdSfV4e1otGcTGBLss4kAW4j+QWFDE6PZOsjXsYe0o7bju9PREVmKeMi/ltCsWYvUWl3DxlPrNW5HFFn0QeOa8L0ZG2+II5NAtwH1iwaQ+j0jPJP1DKi1f25JxuTSv83PjYclMoJqytzt1HaloGG3cV8s8LujKsX6tAl2QczgK8it7LyuGe9xaREBfD9Bv607lZnWN6/sER+D4bgYe1b5Zv5+Yp2URHRjDp+r70bWPNqMzRWYBXkqvMzROfLefV79fRr00DXryyJw0rMU9pc+DhrfxFXp2b1mHCiBSaey/yMuZoLMArIb+wlLFTsvh+1Q6uPrEVD5zbmahKLhJ7cApln02hhJ3CEhd3TVvIRwu38tfuzfjPRd3sVFJzTCzAj9Gq7QWMTMtg854D/PvCE7iiT2KVvl5MZASREWJTKGEmZ3chqWmZLNu2l7uGdOCGk60ZlTl2FuDH4Mul27ll6nxqRkcyZWQ/UpIaVPlrigi1YyJtCiWMzFu7kxsmZVFa5ub1q3tzSsdDX+RlzNEc9e9+EWkpIrNEZKmILBGRm73bG4jIlyKyyvvxkKvShwJV5fmvVzEyLYO2jeP48O8DfBLeB8XFRFJgDa3CQvrcDVz1v3nUqxXF+2MGWHibKqnICNwF3K6qWSISD2SKyJfANcDXqvq4iNwD3APc7b9SA2N/sYs7py3gk0XbGNqjOf++8ARio3w7TxkfG8l+C/CQVuJy8/CHS5g8byOndEjg2St6UOcoF3kZczQVWVJtK7DVe7tARJYBzYHzgUHeh70FfEuIBfimXYWMTMtg5fYC7j+7E9f/pbVf5ilrx1hL2VCWV1DMjZMy+WX9bm4c1Jbbz7BmVMY3jmkOXESSgB7APKBJuVXptwGH7CovIqlAKkBiYtXe8KtOP67ewZjJWZS5lTf+1oeT2yf4bV9xMZHsKSzx29c3gbMoJ5/U9Ax2F5bw3BU9OK97s0CXZEJIhc99E5E4YDpwi6ruLf85VVU8q9P/iapOUNUUVU1JSPBfCPqKqvLmnHUMf/1nGsbF8MHYgX4Nb/B0JLQReOiZmb2Zi1/5kQgRpo3ub+FtfK5CI3ARicIT3pNU9T3v5u0i0lRVt4pIUyDXX0VWl2JXGQ++v5h3MnI4rVMTxl3W/ajNqHwh3qZQQkqZW/nPZ8sZP3stfZIa8NKwntaMyvjFUQNcPJO+rwHLVPXpcp/6ALgaeNz7caZfKqwmuXuLGDUxk/kb93DT4HbcclrFmlH5gp1GGDryC0u5aep8vluZx7B+ifzjXGtGZfynIiPwAcBwYJGIZHu33YcnuN8RkeuADcClfqmwGmRv2sOo9AwKily8fFVPzjqh4s2ofCEuJpL9JWW43VptvzSM763OLWBkWiY5uwv519ATuLJv8LznY4JTRc5C+QE4XKqc6ttyqt/0zBzunbGIxvGeZlSdmh5bMypf+LUjYYmrWqZsjO99tXQ7t7ydTWxUBJNH9qO3D68TMOZwwvZKTFeZm399spzX56yjf9uGvHBlTxrUjg5ILb82tCq2AA825ZtRdWlWhwnDU2hmzahMNQnLAN+9v4SxU7KYs3onfxuQxP1ndyKyks2ofCGufEdCWzUraBSWuLjz3YV8vGgr5yc344mLuvn8Ii9jjiTsAnzFNk8zqm35Rfzn4m5cmtIy0CX9urCxnYkSPMpf5HXf2R0Z+Zc21ozKVLuwCvDPl2zjtrezqR0TydRR/eiZ6Iz2LfExFuDB5Kc1Oxkz2duM6preDOpg/UxMYIRFgLvdynPfrOKZr1bRvWU9JgzvRZM6sYEu61e2qENwUFXSftrA/320lKSGtXh1RAptEuICXZYJYyEf4PuKXdz+TjafL9nORT1b8NjQro6bp/xtYWMLcKcqdpXxj/eX8HbGJk7r1JhxlyXbG84m4EI6wDfu9MxTrsot4MFzO3PtgCRHzlP+trCxBbgT5RYUMTo9k6yNe/j74HbcWo0XeRlzJCEb4D+s2sHYKVmoQtq1fRl4fKNAl3RYNoXiXAs27WFUeib5B0p58cqenNOtei/yMuZIQi7AVZXX56znsY+X0q5xHK+OSKFVw9qBLuuIompEEBMZYW9iOsyM+TncPf23i7w6N6v+i7yMOZKQCvCi0jLun7GY6Vk5nNG5CU9flvzr/LLTxVtHQsdwlbl54rPlvPr9Ovq1acBLV/UK2EVexhxJcKRbBWzL9zSjWrBpD7ecdjw3DT4+qOYp46wjoSPsKSzh71Pm8/2qHVzTP4n7z+lEVAAv8jLmSEIiwDM37Gb0xEwKi12MH96LM7scF+iSjpl1JAy8lds9F3lt2XOAJy46gct6WzMq42xBH+DvZGzigRmLOa5uLBOv60uH4+IDXVKl2Ag8sL5Yso1b386mZnQkU1P70auVNaMyzhe0AV5a5uaxj5fx5o/rGdiuES9c2YN6tYJ3njI+NpKt+UWBLiPsuN3K89+sZtxXK+nWoi7jh/eiaV1rRmWCQ1AG+K79JYyZlMVPa3dy3cDW3HtWx4A2o/IFW9i4+u0vdnH7Owv4bMk2LuzRnH9deILjLvIy5kiCLsCXbd3LyLQMcguK+e8l3bmoV4tAl+QTcTYHXq027iwkNd3TjOqBczpx3cDWjrzIy5gjOeqwVUReF5FcEVlcbtvDIrJZRLK9/872b5keny7ayoUv/UhpmZt3Rp0YMuENtrBxdZqzegfnvfgDW/OLeOvaPlxvnQRNkKrICPxN4AUg7Q/bx6nqUz6v6BDcbuWZr1by3Der6ZFYj/HDetHYQc2ofCEuOpJil5sSl9vWUPQTVeWNOet57JNltGlUm1dHpJDUyNkXeRlzJBVZUm22iCRVQy2H9cDMxUyet5FLU1rw6AVdiYkMvXnKuHL9UKIjg/fNWKcqKi3jgfcXMy0z+C7yMuZwqvIKHisiI4AM4HZV3X2oB4lIKpAKkJhYufNqL0tpSYcm8Yw4sVXI/qkbV64neH276s+ntu8tYlR6Jtmb9nDTqcdzy6nBdZGXMYdT2b/VXwbaAsnAVuC/h3ugqk5Q1RRVTUlISKjUzrq3rMfV/Z3ZSdBX4m1VHr+Yv3E3f33+B1ZuL+CVYT257XTrJGhCR6VG4Kq6/eBtEXkV+MhnFYWp2rYqj8+9m7GJ+2cspkndGNKu60/H46wZlQktlQpwEWmqqlu9d4cCi4/0eHN0cdZS1mdcZW4e+2QZb8xZz4B2DXnhip42LWVC0lEDXESmAIOARiKSAzwEDBKRZECB9cAo/5UYHmwKxTd27y9h7JQs5qzeybUDWnPf2cF/kZcxh1ORs1CuOMTm1/xQS1izKZSqW77Nc5HX9vxinry4G5ektAx0Scb4lZ1H5RA2hVI1ny3eym3vLCAuJpK3R/WjR2L9QJdkjN9ZgDtE7WgbgVeG2608+/Uqnv16Fckt6zF+eC+ahNhFXsYcjgW4Q0REiLWUPUb7il3c9nY2XyzdzsW9WvDPC7paMyoTVizAHaR2TA2bQqmgDTv3MzItgzV5+3nor525JsSvEzDmUCzAHcRG4BXz/ao8xk6ejwikXduHAe0aBbokYwLCAtxB4mKjLMCPQFV57Yd1/OuTZRzfOJ5XR6SQ2LBWoMsyJmAswB0kLqaGBfhhFJWWcd97i3hv/maGdDmO/17a/ddTL40JV/YT4CBxMZHsKCgMdBmOsy2/iFHpGSzIyee209sz9pR21s/EGCzAHSUuxqZQ/ihzw25GT8yksNjFhOG9OKPLcYEuyRjHsAB3kHhbled33v5lIw++v4Sm9WKZdH1f2jeJD3RJxjiKBbiD1PbOgatqWJ8SV1rm5p8fLeWtnzbwl+Mb8fwVPahXy5pRGfNHFuAOEhcTRZlbKXa5w/aClF37S7hxUiZz1+5i5F9ac/cQa0ZlzOFYgDvIwWXVCopcYRngS7fsJTU9g9yCYp6+tDsX9gydRauN8QcLcAeJi/GE9r5iFwnxMQGupnp9vHArd7y7gDo1I3l31Il0b1kv0CUZ43gW4A4SFxMFhFdHQrdbGffVSp7/ZjU9EusxflgvGlszKmMq5KiTiyLyuojkisjictsaiMiXIrLK+9F6d/pAXJj1BC8oKiU1PYPnv1nNpSktmJraz8LbmGNQkXeH3gSG/GHbPcDXqno88LX3vqmicArwdTv2M/SlH5m1Io9HzuvCExd1IyYy/Ob9jamKiqzIM1tEkv6w+Xw8y6wBvAV8C9zty8LCUdyvy6qVBrgS/5q9Mo+xk7OoESGkX9eH/m2tGZUxlVHZOfAm5RY13gY0OdwDRSQVSAVITEys5O7Cw28j8LIAV+Ifqsr/vl/Hvz9dRvsmnmZULRtYMypjKqvKb2KqqoqIHuHzE4AJACkpKYd9nCm3sHEIvolZVFrGve8tYsb8zZx9wnE8ebE1ozKmqir7E7RdRJqq6lYRaQrk+rKocBUTGUGNCAm5KZSt+QdITctk0eZ87jijPWNOaRfWV5oa4yuVDfAPgKuBx70fZ/qsojAm4l1WLYRG4BnrdzF6YhZFpWW8OiKF0zsfdrbNGHOMjhrgIjIFzxuWjUQkB3gIT3C/IyLXARuAS/1ZZDjxrMoTGnPgU3/eyIMzF9O8Xk2mpvalXWNrRmWML1XkLJQrDvOpU31ci+FggAf3FEppmZv/+3Ap6XM3cFL7BJ6/vAd1a0UFuixjQo69i+QwcUHeUnbnvmJumJTFz+t2MeqkNtw1pCM1bPEFY/zCAtxh4mIi2XMgOEfgS7bkk5qWyY59xTxzWTIX9Gge6JKMCWkW4A4TFxtJzu7gW1bto4VbuOPdBdSvFc27o0+kW4t6gS7JmJBnAe4wcdHBNYVS5lae/nIFL85aQ0qr+rw0rCeN462fiTHVwQLcYeJig+c0wr1FpdwyNZtvludyRZ+WPHJeV6IjbfEFY6qLBbjDxMVEsr+kDLdbHb3y+tq8fVyflsHGnYU8ekFXhvVNtItzjKlmFuAOc7Afyv4SF/Gxzjz1btaKXG6aMp+oGhFMvL4v/do0DHRJxoQlC3CH+a0jofMCXFUZP3stT3y2nE7H1WHCiF60qG/NqIwJFAtwh/l1BO6wNzIPlJRx9/SFfLBgC+d0a8qTF3ejVrS9fIwJJPsJdJjyCxs7xeY9B0hNy2Dp1r3ceWYHbhzU1ua7jXEAC3CHcdqqPD+v28UNEzMpcbl57eoUBne0ZlTGOIUFuMP8GuAOGIFPnLuBhz9YQmKDWkwYkUK7xnGBLskYU44FuMM4YQRe4nLz8IdLmDxvI6d0SOCZy3tQt6az3lA1xliAO06gAzyvoJgxk7L4ef0ubhjUljvO6GDNqIxxKAtwh6kdwCmUxZvzSU3LYFdhCc9d0YPzujer9hqMMRVnAe4w0ZERxERGsK+kegN8ZvZm7pq2kIa1o5k2uj9dm9et1v0bY45dlQJcRNYDBUAZ4FLVFF8UFe6qc1m1Mrfyn8+XM/67tfRJasBLw3rSKC6mWvZtjKkaX4zAT1HVHT74OsaruhZ1yD9Qyk1T5vPdyjyG9UvkH+d2sWZUxgQRm0JxoLiYSL9fibk6dx8j0zLYtKuQx4Z25aq+rfy6P2OM71U1wBX4QkQUGK+qE/74ABFJBVIBEhMTq7i78BAXE+nXKzG/XradW6ZmEx0ZweSR/ejTuoHf9mWM8Z+qBvhAVd0sIo2BL0VkuarOLv8Ab6hPAEhJSdEq7i8sxMVEsm1vkc+/rqry0rdreOqLFXRpVofxw1NoXq+mz/djjKkeVQpwVd3s/ZgrIjOAPsDsIz/LHE1cbCT78nw7Ai8scXHntIV8vHAr5yc34/ELu1EzuoZP92GMqV6VDnARqQ1EqGqB9/YZwP/5rLIw5us58E27CklNz2T5tr3ce1ZHUk9qY82ojAkBVRmBNwFmeIMgEpisqp/5pKow58s58J/W7GTM5CxKy9y8fk1vTunQ2Cdf1xgTeJUOcFVdC3T3YS3GKy4mkmKXm9IyN1E1Kndan6qSPncDj3y4lFYNa/G/ESm0SbBmVMaEEjuN0IEO9gTfX+yiXq3oY35+sauMh2YuYeovmzi1Y2PGXZ5MHYet7mOMqToLcAc62NCqoOjYAzy3oIgbJmaRuWE3Y09px22nt3f04sjGmMqzAHegynYkXLBpD6PSM8k/UMqLV/bknG5N/VGeMcYhLMAdqPzCxhU1Y34Od09fREJcDNNv6E/nZnX8VZ4xxiEswB3oWEbgZW7lic+WM2H2Wvq2bsBLV/WkoTWjMiYsWIA7UEWXVcsvLGXslCy+X7WDq09sxQPndq70WSvGmOBjAe5AFZlCWbW9gJFpGWzec4DHLzyBy/tYnxljwo0FuAMdHIEf7mrML5du55ap86kZHcmUkf1ISbJmVMaEIwtwB6od/dtphOW53coLs1bz9Jcr6daiLuOH96JpXWtGZUy4sgB3oIgIoXZ0jd9NoewvdnHHuwv4dPE2LkhuxuMXdSM2yppRGRPOLMAdKi72t2XVNu0qZGRaBiu3F3D/2Z24/i+trRmVMcYC3KniYiLZV+Lix9U7GDM5izK38sbf+nBy+4RAl2aMcQgLcIeKi4lk7pqdfLZ4G60b1eZ/I1JIalQ70GUZYxzEAtyh4mIj2bm/hNM6NWHcZd2Jt2ZUxpg/sAB3qOH9WnFKh8ZcO6C1NaMyxhySBbhDDelqjaiMMUdWpeuuRWSIiKwQkdUico+vijLGGHN0lQ5wEakBvAicBXQGrhCRzr4qzBhjzJFVZQTeB1itqmtVtQSYCpzvm7KMMcYcTVUCvDmwqdz9HO+23xGRVBHJEJGMvLy8KuzOGGNMeX7vPaqqE1Q1RVVTEhLsIhRjjPGVqgT4ZqBlufstvNuMMcZUg6oE+C/A8SLSWkSigcuBD3xTljHGmKOp9HngquoSkbHA50AN4HVVXeKzyowxxhyRqGr17UwkD9hQyac3Anb4sJxgYMccHuyYw0NVjrmVqv7pTcRqDfCqEJEMVU0JdB3VyY45PNgxhwd/HLOtgGuMMUHKAtwYY4JUMAX4hEAXEAB2zOHBjjk8+PyYg2YO3BhjzO8F0wjcGGNMORbgxhgTpIIiwMOh77iIvC4iuSKyuNy2BiLypYis8n6sH8gafUlEWorILBFZKiJLRORm7/ZQPuZYEflZRBZ4j/kR7/bWIjLP+/p+23tlc0gRkRoiMl9EPvLeD+ljFpH1IrJIRLJFJMO7zeevbccHeBj1HX8TGPKHbfcAX6vq8cDX3vuhwgXcrqqdgX7AGO//aygfczEwWFW7A8nAEBHpBzwBjFPVdsBu4LrAleg3NwPLyt0Ph2M+RVWTy5377fPXtuMDnDDpO66qs4Fdf9h8PvCW9/ZbwAXVWZM/qepWVc3y3i7A88PdnNA+ZlXVfd67Ud5/CgwGpnm3h9QxA4hIC+Ac4H/e+0KIH/Nh+Py1HQwBXqG+4yGqiapu9d7eBjQJZDH+IiJJQA9gHiF+zN6phGwgF/gSWAPsUVWX9yGh+Pp+BrgLcHvvNyT0j1mBL0QkU0RSvdt8/tq2RY2DhKqqiITcOZ8iEgdMB25R1b2ewZlHKB6zqpYBySJSD5gBdAxsRf4lIucCuaqaKSKDAlxOdRqoqptFpDHwpYgsL/9JX722g2EEHs59x7eLSFMA78fcANfjUyIShSe8J6nqe97NIX3MB6nqHmAWcCJQT0QODqZC7fU9ADhPRNbjmf4cDDxLaB8zqrrZ+zEXzy/qPvjhtR0MAR7Ofcc/AK723r4amBnAWnzKOw/6GrBMVZ8u96lQPuYE78gbEakJnI5n7n8WcLH3YSF1zKp6r6q2UNUkPD+736jqVYTwMYtIbRGJP3gbOANYjB9e20FxJaaInI1nHu1g3/HHAluR74nIFGAQnpaT24GHgPeBd4BEPG14L1XVP77RGZREZCDwPbCI3+ZG78MzDx6qx9wNz5tXNfAMnt5R1f8TkTZ4RqcNgPnAMFUtDlyl/uGdQrlDVc8N5WP2HtsM791IYLKqPiYiDfHxazsoAtwYY8yfBcMUijHGmEOwADfGmCBlAW6MMUHKAtwYY4KUBbgxxgQpC3ATdESknojc6L3dTESmHe05VdhXsvc0VmMcxwLcBKN6wI0AqrpFVS8+8sOrJBmwADeOZOeBm6AjIgc7Uq4AVgGdVLWriFyDp8NbbeB44CkgGhiOp5Xr2aq6S0Ta4mlRnAAUAiNVdbmIXILnAqoyIB84DVgN1MRzqfe/gY+A54GueLoJPqyqM737HgrUxdOYaaKqPuLf74QJd9bMygSje4Cuqprs7WT4UbnPdcXT2TAWT/jerao9RGQcMALPFb0TgNGqukpE+gIv4enR8Q/gTG8TonqqWiIi/wBSVHUsgIj8C8/l4Nd6L4v/WUS+8u67j3f/hcAvIvKxqmb48ftgwpwFuAk1s7z9xQtEJB/40Lt9EdDN2/2wP/Buuc6HMd6Pc4A3ReQd4D0O7Qw8zZnu8N6PxXNpNMCXqroTQETeAwYCFuDGbyzATagp30/DXe6+G8/rPQJPL+rkPz5RVUd7R+TnAJki0usQX1+Ai1R1xe82ep73x/lIm580fmVvYppgVADEV+aJqroXWOed70Y8untvt1XVear6DyAPTxvjP+7rc+Dv3m6KiEiPcp873bvuYU08c/FzKlOjMRVlAW6CjneaYo54FoB+shJf4irgOhFZACzhtyX6nvQuRLsY+BFYgKftaWfv4rSXAY/iefNyoYgs8d4/6Gc8/c0XAtNt/tv4m52FYowPeM9C+fXNTmOqg43AjTEmSNkI3BhjgpSNwI0xJkhZgBtjTJCyADfGmCBlAW6MMUHKAtwYY4LU/wMtq+v+qc5ByQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "cbm_output.state.to_pandas()[[\"timestep\", \"age\"]].groupby(\n",
    "    \"timestep\"\n",
    ").mean().plot()"
   ]
  }
 ],
 "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.9.7"
  },
  "vscode": {
   "interpreter": {
    "hash": "7036c97a19c395f990150d2191d95cb0b15bafc44a51c61e79499b778f47a5df"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}