{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Frequentist uncertainty vs. Bayesian uncertainty analysis\n",
    "*Mark Bakker, TU Delft & Raoul Collenteur, Hydroconsult, February, 2025 (last updated May 2026)*\n",
    "\n",
    "In this notebook, the fit and uncertainty are compared for Pastas models solved with least squares (frequentist uncertainty) and with MCMC (Bayesian uncertainty). Besides Pastas the following Python Packages have to be installed to run this notebook:\n",
    "\n",
    "- [emcee](https://emcee.readthedocs.io)\n",
    "- [corner](https://corner.readthedocs.io)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import corner\n",
    "import emcee\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import pastas as ps\n",
    "\n",
    "ps.set_log_level(\"ERROR\")\n",
    "ps.show_versions()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. A 'regular' Pastas Model\n",
    "The first step is to create a Pastas Model with a linear `RechargeModel` and a `Gamma` response function to simulate the effect of precipitation and evaporation on the heads. The AR1 noise model is used. We first estimate the model parameters using the standard least-squares approach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head = pd.read_csv(\n",
    "    \"data/B32C0639001.csv\", parse_dates=[\"date\"], index_col=\"date\"\n",
    ").squeeze()\n",
    "head = head[\"1990\":]  # use data from 1990 on for this example\n",
    "\n",
    "evap = pd.read_csv(\"data/evap_260.csv\", index_col=0, parse_dates=[0]).squeeze()\n",
    "rain = pd.read_csv(\"data/rain_260.csv\", index_col=0, parse_dates=[0]).squeeze()\n",
    "\n",
    "ml1 = ps.Model(head)\n",
    "ps.ArNoiseModel(ml1)\n",
    "\n",
    "rm = ps.RechargeModel(\n",
    "    ml1, rain, evap, recharge=ps.rch.Linear(), rfunc=ps.Gamma(), name=\"rch\"\n",
    ")\n",
    "\n",
    "ml1.solve()\n",
    "\n",
    "ax = ml1.plots.results(figsize=(8, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The diagnostics show that the noise meets the statistical requirements for uncertainty analysis reasonably well. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1.plots.diagnostics();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1.parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimated least squares parameters and standard errors are stored for later reference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls_params = ml1.parameters[[\"optimal\", \"stderr\"]].copy()\n",
    "ls_params.rename(columns={\"optimal\": \"ls_opt\", \"stderr\": \"ls_sig\"}, inplace=True)\n",
    "ls_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute prediction interval Pastas\n",
    "pi = ml1.solver.prediction_interval(n=1000)\n",
    "ax = ml1.plot(figsize=(10, 3))\n",
    "ax.fill_between(pi.index, pi.iloc[:, 0], pi.iloc[:, 1], color=\"lightgray\")\n",
    "ax.legend([\"Observations\", \"Simulation\", \"95% Prediction interval\"], ncol=3, loc=2)\n",
    "pi_pasta = np.mean(pi[0.975] - pi[0.025])\n",
    "print(f\"Mean prediction interval width: {pi_pasta:.3f} m\")\n",
    "print(f\"Prediction interval coverage probability: {ps.stats.picp(head, pi): .3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Use the EmceeSolver\n",
    "\n",
    "We will now use MCMC to estimate the model parameters and their uncertainties. The `EmceeSolve` solver wraps the [Emcee](https://lmfit.github.io/lmfit-py/fitting.html#lmfit.minimizer.Minimizer.emcee) package, which implements different versions of MCMC. A good understanding of Emcee helps when using this solver, so it comes recommended to check out their documentation as well.\n",
    "\n",
    "We start by making a `pastas` model with a linear recharge model and a Gamma response function. No noise model is added, as this is taken care of in the likelihood function. The model is solved using the regular solve (least squares) to have a good estimate of the starting values of the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml2 = ps.Model(head)\n",
    "rm = ps.RechargeModel(\n",
    "    ml2, rain, evap, recharge=ps.rch.Linear(), rfunc=ps.Gamma(), name=\"rch\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml2.parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To set up the `EmceeSolve` solver, a number of decisions need to be made:\n",
    "\n",
    "- Select a (log) likelihood function\n",
    "- Select the priors of the parameters\n",
    "- Select the number of steps and thinning\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2a. Create the solver instance\n",
    "\n",
    "The next step is to create an instance of the `EmceeSolve` solver class. At this stage all the settings need to be provided on how the Ensemble Sampler is created (https://emcee.readthedocs.io/en/stable/user/sampler/). Important settings are the `nwalkers`, the `moves`, the `objfunction`. More advanced options are to parallelize the MCMC algorithm (`parallel=True`), and to set a backend to store the results. Here's an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "parameters"
    ]
   },
   "outputs": [],
   "source": [
    "NUM_STEPS = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose the objective function\n",
    "ln_prob = ps.objective.likelihood.GaussianLikelihoodAr1()\n",
    "\n",
    "# Create the EmceeSolver with some settings\n",
    "s = ps.solver.Emcee(\n",
    "    model=ml2,\n",
    "    nwalkers=20,\n",
    "    nsteps=NUM_STEPS,\n",
    "    moves=emcee.moves.DEMove(),\n",
    "    objfunction=ln_prob,\n",
    "    progress_bar=True,\n",
    "    parallel=False,\n",
    ")\n",
    "ml2.parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above code we created an `EmceeSolve` instance with 20 walkers, which take steps according to the `DEMove` move algorithm (see Emcee docs), and a Gaussian likelihood function that assumes AR1 correlated errors. Different objective functions are available, see the Pastas documentation on the different options. \n",
    "\n",
    "Depending on the likelihood function, a number of additional parameters need to be inferred. These parameters are added to the Pastas Model instance. Using the `set_parameter` method of the model, these parameters can be changed. In this example where we use the `GaussianLikelihoodAr1` function, the $\\sigma^2$ and $\\phi$ are estimated; the unknown standard deviation of the errors and the autoregressive parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Priors\n",
    "\n",
    "The first step is to select the priors of the parameters. This is done by using the `ml.set_parameter` method and the `dist` argument (from distribution). Any distribution from `scipy.stats` can be chosen (https://docs.scipy.org/doc/scipy/tutorial/stats/continuous.html), for example `uniform`, `norm`, or `lognorm`. Here, we select normal distributions for the priors. \n",
    "\n",
    "Currently, `pastas` will use the `initial` value of the parameter for the `loc` argument of the distribution (e.g., the mean of a normal distribution), and the `sigma` as the `scale` argument (e.g., the standard deviation of a normal distribution). Only for the parameters with a `uniform` distribution, the `pmin` and `pmax` values are used to determine a uniform prior.\n",
    "\n",
    "<div class=\"alert alert-warning\">\n",
    "<b>Note:</b>\n",
    "This means that either the `pmin` and `pmax` should be set for uniform distributions, or the `sigma` for any other distribution. That is why in this example model was first solved using LeastSquares, in order to obtain estimates for the `stderr` to set the `sigma`. In practice, these could also be set based on expert judgement or information about the parameters.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the initial parameters to optimal values from least squares and a normal distribution\n",
    "# for a good starting point\n",
    "for name in ml1.parameters.index.intersection(ml2.parameters.index):\n",
    "    ml2.set_parameter(name, initial=ml1.parameters.at[name, \"optimal\"])\n",
    "    ml2.set_parameter(name, dist=\"norm\")\n",
    "    ml2.set_parameter(name, sigma=2 * ml1.parameters.at[name, \"stderr\"])\n",
    "\n",
    "ml2.parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2c. Run the solver and solve the model\n",
    "\n",
    "After setting the parameters and creating a EmceeSolve solver instance we are now ready to run the MCMC analysis. We can do this by running `ml.solve`. We can pass the same parameters that we normally provide to this method (e.g., `tmin` or `fit_constant`). Here we use the initial parameters from our least-square solve, and do not fit a noise model, because we take autocorrelated errors into account through the likelihood function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the solver to run MCMC\n",
    "ml2.solve(\n",
    "    initial=False,\n",
    "    tmin=\"1990\",\n",
    "    tune=True,\n",
    "    report=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Posterior parameter distributions\n",
    "\n",
    "The results from the MCMC analysis are stored in the `sampler` object, accessible through `ml.solver.sampler` variable. The object `ml.solver.sampler.flatchain` contains a Pandas DataFrame with $n$ the parameter samples, where $n$ is calculated as follows:\n",
    "\n",
    "$n = \\frac{\\left(\\text{steps}-\\text{burn}\\right)\\cdot\\text{nwalkers}}{\\text{thin}} $\n",
    "\n",
    "## Corner.py\n",
    "Corner is a simple but great python package that makes creating corner graphs easy. A couple of lines of code suffice to create a plot of the parameter distributions and the covariances between the parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Corner plot of the results\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "\n",
    "labels = list(ml2.parameters.index[ml2.parameters.vary])\n",
    "\n",
    "labels = [label.split(\"_\")[1] for label in labels]\n",
    "\n",
    "best = list(ml2.parameters[ml2.parameters.vary].optimal)\n",
    "\n",
    "axes = corner.corner(\n",
    "    ml2.solver.sampler.get_chain(flat=True, discard=(NUM_STEPS // 2)),\n",
    "    quantiles=[0.025, 0.5, 0.975],\n",
    "    labelpad=0.1,\n",
    "    show_titles=True,\n",
    "    title_kwargs={\"fontsize\": 10},\n",
    "    label_kwargs={\"fontsize\": 10},\n",
    "    max_n_ticks=3,\n",
    "    fig=fig,\n",
    "    labels=labels,\n",
    "    truths=best,\n",
    ")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. The trace shows when MCMC converges\n",
    "The walkers take steps in different directions for each step. It is expected that after a number of steps, the direction of the step becomes random, as a sign that an optimum has been found. This can be checked by looking at the autocorrelation, which should be insignificant after a number of steps. Below we just show how to obtain the different chains, the interpretation of which is outside the scope of this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(len(labels), figsize=(10, 7), sharex=True)\n",
    "\n",
    "samples = ml2.solver.sampler.get_chain(flat=True)\n",
    "for i in range(len(labels)):\n",
    "    ax = axes[i]\n",
    "    ax.plot(samples[:, i], \"k\", alpha=0.5)\n",
    "    ax.set_xlim(0, len(samples))\n",
    "    ax.set_ylabel(labels[i])\n",
    "    ax.yaxis.set_label_coords(-0.1, 0.5)\n",
    "\n",
    "axes[-1].set_xlabel(\"step number\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcn_params = pd.DataFrame(index=ls_params.index, columns=[\"mcn_opt\", \"mcn_sig\"])\n",
    "params = ml2.solver.sampler.get_chain(\n",
    "    flat=True, discard=500\n",
    ")  # discard first 500 of every chain\n",
    "for iparam in range(params.shape[1] - 1):\n",
    "    mcn_params.iloc[iparam] = np.median(params[:, iparam]), np.std(params[:, iparam])\n",
    "mean_time_diff = head.index.to_series().diff().mean().total_seconds() / 86400\n",
    "\n",
    "# Translate phi into the value of alpha also used by the noisemodel\n",
    "mcn_params.loc[\"noise_alpha\", \"mcn_opt\"] = -mean_time_diff / np.log(\n",
    "    np.median(params[:, -1])\n",
    ")\n",
    "mcn_params.loc[\"noise_alpha\", \"mcn_sig\"] = -mean_time_diff / np.log(\n",
    "    np.std(params[:, -1])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.concat((ls_params, mcn_params), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Repeat with uniform priors\n",
    "Set more or less uninformative uniform priors. Now also include $\\sigma^2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml3 = ps.Model(head)\n",
    "rm = ps.RechargeModel(\n",
    "    ml3, rain, evap, recharge=ps.rch.Linear(), rfunc=ps.Gamma(), name=\"rch\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Uniform prior selected from 0.25 till 4 times the optimal values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose the objective function\n",
    "ln_prob = ps.objective.likelihood.GaussianLikelihoodAr1()\n",
    "\n",
    "# Create the EmceeSolver with some settings\n",
    "s = ps.solver.Emcee(\n",
    "    model=ml3,\n",
    "    nwalkers=20,\n",
    "    nsteps=NUM_STEPS,\n",
    "    moves=emcee.moves.DEMove(),\n",
    "    objfunction=ln_prob,\n",
    "    progress_bar=True,\n",
    "    parallel=False,\n",
    ")\n",
    "\n",
    "for name in ml1.parameters.index.intersection(ml3.parameters.index):\n",
    "    opt = ml1.parameters.loc[name, \"optimal\"]\n",
    "    if not np.isnan(opt):\n",
    "        sign = np.sign(opt)\n",
    "        ml3.set_parameter(\n",
    "            name,\n",
    "            initial=sign * opt,\n",
    "            dist=\"uniform\",\n",
    "            pmin=sign * 0.25 * opt,\n",
    "            pmax=sign * 4 * opt,\n",
    "            sigma=ml1.parameters.loc[name, \"stderr\"] ** 2,\n",
    "        )\n",
    "\n",
    "ml3.set_parameter(\n",
    "    \"solver_var\",\n",
    "    initial=0.05**2,\n",
    "    pmin=0.05**2 / 4,\n",
    "    pmax=4 * 0.05**2,\n",
    "    vary=True,\n",
    "    dist=\"uniform\",\n",
    ")\n",
    "ml3.set_parameter(\"solver_phi\", dist=\"uniform\")\n",
    "ml3.parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the solver to run MCMC\n",
    "ml3.solve(\n",
    "    initial=False,\n",
    "    tmin=\"1990\",\n",
    "    tune=True,\n",
    "    report=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Corner plot of the results\n",
    "fig = plt.figure(figsize=(8, 8))\n",
    "\n",
    "labels = list(ml3.parameters.index[ml3.parameters.vary])\n",
    "labels = [label.split(\"_\")[1] for label in labels]\n",
    "\n",
    "best = list(ml3.parameters[ml3.parameters.vary].optimal)\n",
    "\n",
    "axes = corner.corner(\n",
    "    ml3.solver.sampler.get_chain(flat=True, discard=(NUM_STEPS // 2)),\n",
    "    quantiles=[0.025, 0.5, 0.975],\n",
    "    labelpad=0.1,\n",
    "    show_titles=True,\n",
    "    title_kwargs={\"fontsize\": 10},\n",
    "    label_kwargs={\"fontsize\": 10},\n",
    "    max_n_ticks=3,\n",
    "    fig=fig,\n",
    "    labels=labels,\n",
    "    truths=best,\n",
    ")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(len(labels), figsize=(10, 7), sharex=True)\n",
    "\n",
    "samples = ml3.solver.sampler.get_chain(flat=True)\n",
    "for i in range(len(labels)):\n",
    "    ax = axes[i]\n",
    "    ax.plot(samples[:, i], \"k\", alpha=0.5)\n",
    "    ax.set_xlim(0, len(samples))\n",
    "    ax.set_ylabel(labels[i])\n",
    "    ax.yaxis.set_label_coords(-0.1, 0.5)\n",
    "\n",
    "axes[-1].set_xlabel(\"step number\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mcu_params = pd.DataFrame(index=ls_params.index, columns=[\"mcu_opt\", \"mcu_sig\"])\n",
    "params = ml3.solver.sampler.get_chain(\n",
    "    flat=True, discard=500\n",
    ")  # discard first 500 of every chain\n",
    "for iparam in range(params.shape[1] - 1):\n",
    "    mcu_params.iloc[iparam] = np.median(params[:, iparam]), np.std(params[:, iparam])\n",
    "mean_time_diff = head.index.to_series().diff().mean().total_seconds() / 86400\n",
    "mcu_params.loc[\"noise_alpha\", \"mcu_opt\"] = -mean_time_diff / np.log(\n",
    "    np.median(params[:, -1])\n",
    ")\n",
    "mcu_params.loc[\"noise_alpha\", \"mcu_sig\"] = -mean_time_diff / np.log(\n",
    "    np.std(params[:, -1])\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.concat((ls_params, mcn_params, mcu_params), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Compute prediction interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nobs = len(head)\n",
    "params = ml3.solver.sampler.get_chain(flat=True, discard=(NUM_STEPS // 2))\n",
    "sim = {}\n",
    "# compute for 1000 random samples of chain\n",
    "np.random.seed(1)\n",
    "for i in np.random.choice(np.arange(NUM_STEPS) * 10, size=NUM_STEPS, replace=False):\n",
    "    h = ml3.simulate(p=params[i, :-2])\n",
    "    res = ml3.residuals(p=params[i, :-2])\n",
    "    h += np.random.normal(loc=0, scale=np.std(res), size=len(h))\n",
    "    sim[i] = h\n",
    "simdf = pd.DataFrame.from_dict(sim, orient=\"columns\", dtype=float)\n",
    "alpha = 0.05\n",
    "q = [alpha / 2, 1 - alpha / 2]\n",
    "pi = simdf.quantile(q, axis=1).transpose()\n",
    "pimean = np.mean(pi[0.975] - pi[0.025])\n",
    "print(f\"prediction interval emcee with uniform priors: {pimean:.3f} m\")\n",
    "print(f\"PICP: {ps.stats.picp(head, pi):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, the prediction interval is dominated by the residuals not by the uncertainty of the parameters. In the code cell below, the parameter uncertainty is not included: the coverage only changes slightly and is mostly affected by the difference in randomly drawing residuals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logprob = ml3.solver.sampler.compute_log_prob(\n",
    "    ml3.solver.sampler.get_chain(flat=True, discard=(NUM_STEPS // 2))\n",
    ")[0]\n",
    "imax = np.argmax(logprob)  # parameter set with larges likelihood\n",
    "nobs = len(head)\n",
    "params = ml3.solver.sampler.get_chain(flat=True, discard=(NUM_STEPS // 2))\n",
    "sim = {}\n",
    "# compute for 1000 random samples of residuals, but one parameter set\n",
    "h = ml3.simulate(p=params[imax, :-2])\n",
    "res = ml3.residuals(p=params[imax, :-2])\n",
    "np.random.seed(1)\n",
    "for i in range(NUM_STEPS):\n",
    "    sim[i] = h + np.random.normal(loc=0, scale=np.std(res), size=len(h))\n",
    "simdf = pd.DataFrame.from_dict(sim, orient=\"columns\", dtype=float)\n",
    "alpha = 0.05\n",
    "q = [alpha / 2, 1 - alpha / 2]\n",
    "pi = simdf.quantile(q, axis=1).transpose()\n",
    "pimean = np.mean(pi[0.975] - pi[0.025])\n",
    "print(f\"prediction interval emcee with uniform priors: {pimean:.3f} m\")\n",
    "print(f\"PICP: {ps.stats.picp(head, pi):.3f}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pastas (3.13.14.final.0)",
   "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.-1"
  },
  "mystnb": {
   "execution_mode": "off"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}