{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial about pairwise distance analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pairwise distance distribution p(r) - as derived from a histogram of\n",
    "pairwise distances - represents the probability distribution function\n",
    "to find for a localization at r = 0 another localization at distance r + delta_r."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import locan as lc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc.show_versions(system=False, dependencies=False, verbose=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(seed=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate localization data at two different intensities (localization density) that is (i) homogeneously Poisson distributed (also described as complete spatial randomness, csr) and that (ii) follows a Neyman-Scott distribution (blobs)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_csr_0 = lc.simulate_Poisson(intensity=1e-3, region=((0,1000), (0,1000)), seed=rng)\n",
    "locdata_csr_1 = lc.simulate_Poisson(intensity=1e-2, region=((0,1000), (0,1000)), seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_blob_0 = lc.simulate_Thomas(parent_intensity=1e-4, region=((0, 1000), (0, 1000)), cluster_mu=10, cluster_std=5, seed=rng)\n",
    "locdata_blob_1 = lc.simulate_Thomas(parent_intensity=1e-3, region=((0, 1000), (0, 1000)), cluster_mu=10, cluster_std=5, seed=rng)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of localizations:\")\n",
    "print(\"csr_0:\", len(locdata_csr_0))\n",
    "print(\"csr_1:\", len(locdata_csr_1))\n",
    "print(\"blob_0:\", len(locdata_blob_0))\n",
    "print(\"blob_1:\", len(locdata_blob_1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=2)\n",
    "locdata_csr_0.data.plot.scatter(x='position_x', y='position_y', ax=axes[0, 0], color='Blue', s=1, alpha=0.1, label='locdata_csr')\n",
    "locdata_csr_1.data.plot.scatter(x='position_x', y='position_y', ax=axes[0, 1], color='Blue', s=1, alpha=0.1, label='locdata_csr')\n",
    "locdata_blob_0.data.plot.scatter(x='position_x', y='position_y', ax=axes[1, 0], color='Blue', s=1, alpha=0.1, label='locdata_blobs')\n",
    "locdata_blob_1.data.plot.scatter(x='position_x', y='position_y', ax=axes[1, 1], color='Blue', s=1, alpha=0.1, label='locdata_blobs')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pairwise distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We determine all pairwise distances and plot the pair distance probability distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd_csr_0 = lc.PairDistances().compute(locdata_csr_0)\n",
    "pd_csr_1 = lc.PairDistances().compute(locdata_csr_1)\n",
    "pd_blob_0 = lc.PairDistances().compute(locdata_blob_0)\n",
    "pd_blob_1 = lc.PairDistances().compute(locdata_blob_1)\n",
    "\n",
    "pd_csr_0.results.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd_csr_0.hist(alpha=0.5, label=\"csr_0\")\n",
    "pd_blob_0.hist(alpha=0.5, label=\"blob_0\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd_csr_1.hist(alpha=0.5, label=\"csr_1\")\n",
    "pd_blob_1.hist(alpha=0.5, label=\"blob_1\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relative pairwise distance distribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A pairwise distance distribution relative to the expected distribution for a homogeneous sample (csr) reveals clustering effects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.linspace(0, 100, 100)\n",
    "hist_csr, bin_edges_csr = np.histogram(pd_csr_1.results.pair_distance, bins=bins, density=True)\n",
    "hist_blob, bin_edges_blob = np.histogram(pd_blob_1.results.pair_distance, bins=bins, density=True)\n",
    "bin_widths = np.diff(bin_edges_blob)\n",
    "values = hist_blob / hist_csr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.bar(x=bin_edges_blob[:-1], height=values, align=\"edge\", width=bin_widths, label=\"blob_1\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.12.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}