notebooks/1_getting_started.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f71abee3-f80e-4fd3-8270-8998c1b8da02",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/tee-lab/PyDaddy/blob/colab/notebooks/1_getting_started.ipynb)\n",
    "\n",
    "# Introduction to PyDaddy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7b63b7c-dbaa-4446-94e9-6fadd6d8f253",
   "metadata": {},
   "source": [
    "`pydaddy` is a Python toolbox to derive stochastic differential equations (SDEs) from time-series data. Given samples of a time-series $x(t)$, `pydaddy` attempts to fit functions $f$ and $g^2$ such that\n",
    "\n",
    "$$ \\frac{dx}{dt} = f(x) + g(x) \\cdot \\eta(t) $$\n",
    "\n",
    "where $\\eta(t)$ is uncorrelated white noise. The function $f$ is called the _drift_, and governs the deterministic part of the dynamics. $g^2$ is called the _diffusion_ and governs the stochastic part of the dynamics.\n",
    "\n",
    "PyDaddy estimates the drift function $f$ directly. For diffusion, PyDaddy estimates $g^2$ and not $g$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "258b23f9-9e19-4821-8303-f6df919302ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pydaddy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdd838a7-aed8-4d01-8186-75d8f73dbcf4",
   "metadata": {},
   "source": [
    "## Initializing a `pydaddy` object\n",
    "\n",
    "To start analysis, we need to create a `pydaddy` object with our dataset. This will compute the drift and diffusion parts, and generate a summary plot. To initialize a `pydaddy` object, we need to provide the following arguments:\n",
    " - `data`: the timeseries data, could be either one or two dimensional. This example will deal with 1D data, see [Getting Started with Vector Data](./2_getting_started_vector.ipynb) for a 2D example. `pydaddy` assumes that the samples are evenly spaced. `data` should be a list of Numpy arrays; with one array for the scalar case, and two arrays for the vector case.\n",
    " - `t`: This could either be scalar, denoting the time-interval between samples, or an array denoting the timestamp of each sample.\n",
    " - `bins`: The number of bins to use for computing the average drift and diffusion. Binning is only done for visualization purposes.\n",
    " \n",
    "There are also other optional arguments: see [documentation](https://pydaddy.readthedocs.io/en/latest/api.html) for detailed descriptions of all arguments."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e168e420-3b4e-46c3-a503-da93b380d51e",
   "metadata": {},
   "source": [
    "This example uses a sample dataset, loaded using a helper function. The dataset contains net polarisation from a 1-dimensional agent-based simulation. For this model, the drift and diffusion functions are known -- the drift function is linear while the diffusion is quadratic. In this notebook, we will examine if PyDaddy can estimate the drift and diffusion correctly from this dataset.\n",
    "\n",
    "For more information about the model and the dataset, see [Jhawar and Guttal, 2020](https://doi.org/10.1098/rstb.2019.0381).\n",
    "For details about data formats and loading/saving data, see [Exporting Data](./5_exporting_data.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30c915a2-ff29-4225-bb64-dd755e3421ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "data, t = pydaddy.load_sample_dataset('model-data-scalar-pairwise')\n",
    "ddsde = pydaddy.Characterize(data, t, bins=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "174ff7ca-7663-4275-b2de-7c9a00307faf",
   "metadata": {},
   "source": [
    "`pydaddy.Characterize` initializes a `ddsde` object which can be used for further analysis. It also produces summary plots, showing the time-series, histograms, and the estimated drift and diffusion functions.\n",
    "\n",
    "`pydaddy` can automatically try to fit polynomial functions if called with argument `fit_functions=True`. However, for best results, it is recommended to do the fitting separately, with some level of manual intervention. See [Advanced Function Fitting](./3_advanced_function_fitting.ipynb)  for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "feb39115-17df-466d-98d7-2415141279f7",
   "metadata": {},
   "source": [
    "## Recovering functional forms for $f$ and $g$\n",
    "\n",
    "`pydaddy` has a `fit()` function which can recover functional forms for the drift and diffusion functions, using sparse regression. By default, `pydaddy` fits polynomials (of a specified degree), but it is possible to fit arbitrary functions by specifying a custom library (see [Fitting Non-Polynomial Functions](./6_non_poly_function_fitting.ipynb)).\n",
    "\n",
    "Two parameters need to be specified during fitting:\n",
    " - `order`: The maximum degree of the polynomial to be fitted (see [Advanced Function Fitting](./3_advanced_function_fitting.ipynb) for some tips on how to choose the correct order).\n",
    " - `threshold`: a _sparsification threshold_, that governs the level of sparsity (i.e. the number of terms in the polynomial). For `threshold=theta`, the fitted polynomial will only have terms with coefficients greater than `theta`. \n",
    " \n",
    "We can ask `pydaddy` to try to automatically find an appropriate sparsification threshold by calling `fit()` with argument `tune=True`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c567f3a0-af13-4025-853a-39653da43a35",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fitting with automatic threshold tuning\n",
    "F = ddsde.fit('F', order=3, tune=True)\n",
    "print(F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7adaa8bd-5bd8-4827-acb8-90fef1c5974e",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = ddsde.fit('G', order=3, tune=True)\n",
    "print(G)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9022f125-37f4-4365-8f74-ae057887ffff",
   "metadata": {},
   "source": [
    "In the above example, automatic model selection (`tune=True`) found a reasonable linear function for drift.\n",
    "For diffusion, it looks like there are insignificant $x$ and $x^3$ terms, which can be eliminated by a correct choice of threshold.\n",
    "\n",
    "If the data is too noisy, or if `order` is too high, automatic model selection can give poor results. In such cases, good results can be obtained with some manual intervention: see [Advanced Function Fitting](./3_advanced_function_fitting.ipynb) for more details.\n",
    "\n",
    "In cases like this, we can manually set the appropriate threshold using the `threshold` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1686ec8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsde.fit('F', order=3, threshold=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9caff520-bbc8-40a9-8bac-20d4feecc26e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsde.fit('G', order=2, threshold=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89a113d0-465a-478c-b860-0762f9176df7",
   "metadata": {},
   "source": [
    "The discovered drift and diffusion functions are consistent with the theoretical expectations for the model (see [Jhawar and Guttal, 2020](https://doi.org/10.1098/rstb.2019.0381))."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4895a95-de51-4c6a-9d6e-420cffba6249",
   "metadata": {},
   "source": [
    "The fitted functions can also be printed individually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43c45e8d-989a-4ed6-aa0c-d8b6be7755f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ddsde.F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "429d3087-a876-4392-b83a-c5be3badf6a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ddsde.G)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8888925-f630-4fd7-b755-9d4f0025b1f4",
   "metadata": {},
   "source": [
    "`dddsde.F` and `ddsde.G` are, in fact, callable functions: this means that you can call `ddsde.F` or `ddsde.G` with some desired arguments and evaluate the drift or diffusion for that value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0775b11b-4aa6-40e2-9fa3-637cf212b2fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsde.F(0.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cde5d8d-d2cf-44e5-b2ff-f6b79d475dfb",
   "metadata": {},
   "source": [
    "## Interactive plots for drift and diffusion\n",
    "\n",
    "To get interactive plots for the drift and diffusion functions, use `ddsde.drift()` or `ddsde.diffusion()`. These will be particularly useful for the 2-D case, where the drift and diffusion plots will be 3-D plots (see [Getting Started with Vector Data](./2_getting_started_vector.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcf8daaf-d017-463d-b7c0-5b0a41a0c246",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsde.drift()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5769567",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsde.diffusion()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d97c0051",
   "metadata": {},
   "source": [
    "## Diagnostics\n",
    "\n",
    "For a drift-diffusion model fit to be valid, the data should satisfy some underlying assumptions. `pydaddy.noise_diagnostics()` allows us to verify if the data satisfies these assumptions.\n",
    "The function produces 4 plots:\n",
    "\n",
    "- The distribution of the residuals, which should be a Gaussian.\n",
    "- QQ plot of the residual distribution, against a theoretical Gaussian distribution of the same mean and variance. Ideally (i.e. if the residuals are Gaussian distributed), all points of this plot should fall on a straight line of slope 1.\n",
    "- Autocorrelation plot of the residuals. Ideally, the residuals should be uncorrelated, i.e. autocorrelation time should be close to 0.\n",
    "- The plot of the 2nd and 4th Kramer-Moyal coefficients. Theory dictates that $\\text{KM}(4)$ should equal 3 times $\\text{KM}(2)$, i.e. the plot should be a straight line of slope 1.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7d61c9a-b32c-4403-9472-b591125e21b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsde.noise_diagnostics()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2291a4d1",
   "metadata": {},
   "source": [
    "The `model_diagnostics()` functions checks if the model is self-consistent. \n",
    "\n",
    "To do this, a simulated time series, with the same length and sampling time as the original time series, is generated by integrating the discovered SDE. The drift and diffusion functions are now re-estimated from this simulated time series, with the same fitting parameters as the original fit. If the model is self-consistent, the re-estimated drift and diffusion functions should match the original drift and diffusion.\n",
    "\n",
    "The `oversample` parameter is used to do the SDE simulation at a finer integration time-step than the sampling time of the data. For example, `oversample=5` means that the simulated time-series will be generated with an integration time-step of `t_int / 5`, and then subsampled to `t_int` for the rest of the diagnostics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f39923e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ddsde.model_diagnostics(oversample=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e48c69ee-c511-4e40-90cd-5b525515ae8f",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "- To see PyDaddy in action on a 2-dimensional vector dataset, see [Getting Started with Vector Data](./2_getting_started_vector.ipynb).\n",
    "- To learn more about how to manually fit the drift and diffusion functions, see [Advanced Function Fitting](./3_advanced_function_fitting.ipynb).\n",
    "- To see PyDaddy in action on real world datasets, see [Example analysis: mesoscopic SDEs for fish schools](./7_example_fish_school.ipynb) and [Example analysis: SDEs for cancer cell migration](./8_example_cell_migration.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c37a2937-1ccb-49f1-9215-f9620392f3b2",
   "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.11.0"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}