{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter estimation using non-linear semi-quantitative data"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we illustrate how to use non-linear semi-quantitatiave data for parameter optimization in pyPESTO. An example model is provided in `pypesto/doc/example/example_censored`.\n",
    "\n",
    "Some measurement techniques do not ensure a linear relationship between the abundance of the biochemical quantities of interest and the measured output . Well-known examples include fluorescence microscopy data like F ̈orster resonance energy transfer (FRET) data (Birtwistle et al. [2011]), optical density (OD) measurement (Stevenson et al. [2016]) and imaging data for certain stainings (Pargett et al. [2014]).\n",
    "\n",
    "Assuming a Gaussian distributed additive noise model, a non-linear semi-quantitative type observable has the following relationship between the biochemical quantities of intereset $(y_k = h(x(t_k, \\theta), \\theta))_{k=1}^{n_t}$ and the measured quantities $(\\widetilde{z}_k)_{k=1}^N$\n",
    "$$\\widetilde{z}_k = g(y_k) + \\varepsilon_k, \\quad k = 1, ..., N$$\n",
    "Where:\n",
    "\n",
    "- $\\{y_k\\}_{k=1}^N$ is the model output at timepoints $\\{t_k\\}_{k=1}^N$, \n",
    "- $\\{\\varepsilon_k\\}_{k=1}^N$ is the normally distributed measurement noise, \n",
    "- and $\\theta$ is the vector of model (unknown) mechanistic model parameters.\n",
    "\n",
    "In pyPESTO, we have implemented an alogorithm which constructs and optimizes a spline approximation $s(y, \\xi)$ (to be exact, a piecewise linear function) of the non-linear monotone function $g(y_k)$. \n",
    "\n",
    "Using this spline appoximation of the function $g$ we can then obtain spline-mapped simulations $\\{z_k = s(y_k, \\xi)\\}_{k=1}^N$ that are comparable to the measured quantities. Thus, we can define a negative log-likelihood objective function and optimize it to obtain maximum likelihood estimates. For better efficiency, the objective function is optimized hierarchically\n",
    "\n",
    "- The dynamical parameters $\\theta$ are optimized in the outer hierarchical loop,\n",
    "- The spline parameters $\\xi$ are optimized in the inner loop for each iteration of the outer one.\n",
    "\n",
    "In the following we will demonstrate how to use the spline approximation approach for integration of nonlinear-monotone data."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem specification & importing model from the petab_problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import petab\n",
    "\n",
    "import pypesto\n",
    "import pypesto.logging\n",
    "import pypesto.optimize as optimize\n",
    "import pypesto.petab\n",
    "from pypesto.C import LIN, InnerParameterType\n",
    "from pypesto.hierarchical.semiquantitative import (\n",
    "    SemiquantInnerSolver,\n",
    "    SemiquantProblem,\n",
    ")\n",
    "from pypesto.hierarchical.semiquantitative.parameter import (\n",
    "    SplineInnerParameter,\n",
    ")\n",
    "from pypesto.visualize import (\n",
    "    plot_splines_from_inner_result,\n",
    "    visualize_estimated_observable_mapping,\n",
    ")\n",
    "from pypesto.visualize.model_fit import visualize_optimized_model_fit"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the spline approach is implemented in the hierarchical manner, it requires us to specify `hierarchical=True` to the petab importer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "petab_folder = \"./example_semiquantitative/\"\n",
    "yaml_file = \"example_semiquantitative.yaml\"\n",
    "\n",
    "petab_problem = petab.Problem.from_yaml(petab_folder + yaml_file)\n",
    "\n",
    "importer = pypesto.petab.PetabImporter(petab_problem, hierarchical=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate the noise parameters of the nonlinear-monotone observables, those noise parameters in the `parameter.tsv` file have to have column entries `estimate` with value 1 and column entry `parameterType` with value `sigma`, as in case of using relative (scaling+offset) data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from pandas import option_context\n",
    "\n",
    "noise_parameter_file = \"parameters_example_semiquantitative_noise.tsv\"\n",
    "# load the csv file\n",
    "noise_parameter_df = pd.read_csv(petab_folder + noise_parameter_file, sep=\"\\t\")\n",
    "with option_context(\"display.max_colwidth\", 400):\n",
    "    display(noise_parameter_df)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only difference in the PEtab formulation is in the `measurementType` column of the measurement file. The non-linear (monotone) semi-quantitative measurements have to be specified file by adding  `semiquantitative` in this column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with option_context(\"display.max_colwidth\", 400):\n",
    "    display(petab_problem.measurement_df)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Note on inclusion of additional data types:\n",
    "It is possible to include observables with different types of data to the same `petab_problem`. Refer to the notebooks on using [relative data](relative_data.ipynb), [ordinal data](ordinal_data.ipynb) and [censored data](censored_data.ipynb) for details on integration of other data types. If the `measurementType` column is left empty for all measurements of an observable, the observable will be treated as quantitative."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constructing the objective and pypesto problem"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spline approximation options that can be specified are the following\n",
    "\n",
    "- `spline_ratio`: float value, determines the number of spline knots as `n_spline_pars` = `spline_ratio` * `n_datapoints`  \n",
    "- `min_diff_factor` : float value, determines the minimal difference between consecutive spline heights as `min_diff_factor` * `measurement_range` / `n_spline_pars`. Set it to 0.0 to disable.\n",
    "- `regularize_spline`: boolean value, determines whether the spline is regularized by a linear function. In most cases, this this option improves method convergence. The goal is to reduce false positive data non-linear splines by pushing them towards a linear line.\n",
    "- `regularization_factor`: float value, determines the strength of the regularization. It is the $\\lambda$ factor in the objective function $J(\\theta, \\xi) = -\\log \\mathcal{L}_{\\mathcal{D}}(\\theta, \\xi) + \\lambda R(\\xi)$, where $-\\log \\mathcal{L}_{\\mathcal{D}}(\\theta, \\xi)$ is the negative log-likelihood and $R(\\xi)$ is the regularization.\n",
    "\n",
    "The `min_diff_factor` is a multiplier of the minimal difference between spline heights. Positive values act as a penalty in the objective function for incorrect orderings; this penalty can improve convergence for most models. However, high `min_diff_factor` values will reduce the spline's ability to approximate functions with flat regions accurately. This issue will be demonstrated in the last section."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now when we construct the `objective`, it will construct all objects of the spline inner optimization for semi-quantitative data:\n",
    "\n",
    "- `SemiquantInnerSolver`\n",
    "- `SemiquantCalculator`\n",
    "- `SemiquantProblem`\n",
    "\n",
    "Specifically, the `SemiquantInnerSolver` and `SemiquantProblem` will be constructed with default settings of \n",
    "\n",
    "- `spline_ratio` = 1/2\n",
    "- `use_minimal_difference` = 0.0\n",
    "- `regularize_spline` = False\n",
    "- `regularization_factor` = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "factory = importer.create_objective_creator()\n",
    "objective = factory.create_objective(verbose=False)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To give non-default options to the `OptimalScalingInnerSolver` and `OptimalScalingProblem`, one can pass them as arguments when constructing the `objective`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = factory.create_objective(\n",
    "    inner_options={\n",
    "        \"spline_ratio\": 1 / 2,\n",
    "        \"min_diff_factor\": 1 / 2,\n",
    "        \"regularize_spline\": True,\n",
    "        \"regularization_factor\": 1e-1,\n",
    "    },\n",
    "    verbose=False,\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, one can even pass them to the importer constructor `pypesto.petab.PetabImporter()`.\n",
    "\n",
    "If changing the `spline_ratio` setting, one has to re-construct the objective object, as this requires a constuction of the new `SemiquantProblem` object with the requested amount of inner parameters."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's construct the pyPESTO problem and optimizer. We're going to use a gradient-based optimizer for a faster optimization, but gradient-free optimizers can be used in the same way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "problem = importer.create_problem(objective)\n",
    "\n",
    "engine = pypesto.engine.MultiProcessEngine(n_procs=3)\n",
    "\n",
    "optimizer = optimize.ScipyOptimizer(\n",
    "    method=\"L-BFGS-B\",\n",
    "    options={\"disp\": None, \"ftol\": 2.220446049250313e-09, \"gtol\": 1e-5},\n",
    ")\n",
    "n_starts = 3"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running optimization using spline approximation"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now running optimization is as simple as running usual pyPESTO miminization:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(n_starts)\n",
    "\n",
    "result = optimize.minimize(\n",
    "    problem, n_starts=n_starts, optimizer=optimizer, engine=engine\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The model optimization has good convergence with a plateu at the optimal point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pypesto.visualize import parameters, waterfall\n",
    "\n",
    "waterfall([result], size=(6, 3))\n",
    "plt.show()\n",
    "parameters([result], size=(6, 3))\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize the estimated spline observable mapping using the `pypesto.visualize.observable_mapping.visualize_estimated_observable_mapping` routine. The routine plots all estimated linear or spline observable mappings of relative or semiquantitative observables, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visualize_estimated_observable_mapping(result, problem, figsize=(10, 8))\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also plot the model fit (observable trajectories) with the spline-mapped simulations included using `visualize_optimized_model_fit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visualize_optimized_model_fit(\n",
    "    petab_problem,\n",
    "    result,\n",
    "    problem,\n",
    ")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Caution when using minimal difference"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate that minimal difference sometimes has a negative effect we will apply it to a very simple synthetic \"model\" -- simulation of the exponential function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timepoints = np.linspace(0, 10, 11)\n",
    "function = np.exp\n",
    "\n",
    "simulation = timepoints\n",
    "sigma = np.full(len(timepoints), 1)\n",
    "\n",
    "# Create synthetic data as the exponential function of timepoints\n",
    "data = function(timepoints)\n",
    "\n",
    "spline_ratio = 1 / 2\n",
    "n_spline_pars = int(np.ceil(spline_ratio * len(timepoints)))\n",
    "\n",
    "\n",
    "par_type = \"spline\"\n",
    "mask = [np.full(len(simulation), True)]\n",
    "\n",
    "inner_parameters = [\n",
    "    SplineInnerParameter(\n",
    "        inner_parameter_id=f\"{par_type}_{1}_{par_index + 1}\",\n",
    "        inner_parameter_type=InnerParameterType.SPLINE,\n",
    "        scale=LIN,\n",
    "        lb=-np.inf,\n",
    "        ub=np.inf,\n",
    "        ixs=mask,\n",
    "        index=par_index + 1,\n",
    "        group=1,\n",
    "        observable_id=\"observable_1\",\n",
    "    )\n",
    "    for par_index in range(n_spline_pars)\n",
    "]\n",
    "\n",
    "inner_problem = SemiquantProblem(\n",
    "    xs=inner_parameters, data=[data], edatas=None, spline_ratio=spline_ratio\n",
    ")\n",
    "\n",
    "options = {\n",
    "    \"minimal_diff_on\": {\n",
    "        \"min_diff_factor\": 1 / 2,\n",
    "    },\n",
    "    \"minimal_diff_off\": {\n",
    "        \"min_diff_factor\": 0.0,\n",
    "    },\n",
    "}\n",
    "inner_solvers = {}\n",
    "results = {}\n",
    "\n",
    "for minimal_diff, option in options.items():\n",
    "    inner_solvers[minimal_diff] = SemiquantInnerSolver(\n",
    "        options=option,\n",
    "    )\n",
    "    print(f\"Using {minimal_diff} options: {option}\")\n",
    "\n",
    "    # Solve the inner problem to obtain the optimal spline\n",
    "    results[minimal_diff] = inner_solvers[minimal_diff].solve(\n",
    "        problem=inner_problem,\n",
    "        sim=[simulation],\n",
    "        amici_sigma=[sigma],\n",
    "    )\n",
    "\n",
    "    plot_splines_from_inner_result(\n",
    "        inner_problem,\n",
    "        inner_solvers[minimal_diff],\n",
    "        results[minimal_diff],\n",
    "        sim=[simulation],\n",
    "    )\n",
    "    plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimized spline for the case with enabled minimal difference is performing much worse even if we use a non-extreme `min_diff_factor` value. This is due to the relative flatness of the data with respect to the true model output.\n",
    "\n",
    "The minimal difference is determined as \n",
    "$$\\text{min diff} = \\text{min diff factor} \\cdot \\frac{\\text{measurement range}}{\\text{n inner pars}}$$ \n",
    "so for nonlinear-monotone functions which are relatively flat on some intervals, it is best to keep the minimal difference disabled.\n",
    "\n",
    "As the true output (e.g. observable simulation of the model with true parameters) is mostly a-priori not known, it's hard to know whether the minimal difference is going to have a bad or good effect on the optimization. So a good heuristic is to run for different values of `min_diff_factor` and compare results."
   ]
  }
 ],
 "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.12.3"
  },
  "vscode": {
   "interpreter": {
    "hash": "b4f64b1cfeae9987d9a74471fe6faf49d769577c41c664ee1b5af662a144b184"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}