{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial about local density analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A local density is computed from the number of neighboring localizations\n",
    "within a specified radius."
   ]
  },
  {
   "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": "markdown",
   "metadata": {},
   "source": [
    "## Synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate localization data that is homogeneously Poisson distributed (also described as complete spatial randomness, csr)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(seed=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_csr = lc.simulate_Poisson(intensity=1e-1, region=((0,100), (0,100)), seed=rng)\n",
    "\n",
    "print('Data head:')\n",
    "print(locdata_csr.data.head(), '\\n')\n",
    "print('Summary:')\n",
    "locdata_csr.print_summary()\n",
    "print('Properties:')\n",
    "print(locdata_csr.properties)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also simulate data that follows a Neyman-Scott distribution (blobs): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata_blob = lc.simulate_Thomas(parent_intensity=1e-3, region=((0, 100), (0, 100)), cluster_mu=100, cluster_std=5, seed=rng)\n",
    "\n",
    "print('Data head:')\n",
    "print(locdata_blob.data.head(), '\\n')\n",
    "print('Summary:')\n",
    "locdata_blob.print_summary()\n",
    "print('Properties:')\n",
    "print(locdata_blob.properties)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scatter plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=2)\n",
    "locdata_csr.data.plot.scatter(x='position_x', y='position_y', ax=ax[0], color='Blue', s=1, alpha=0.1, label='locdata_csr')\n",
    "locdata_blob.data.plot.scatter(x='position_x', y='position_y', ax=ax[1], color='Blue', s=1, alpha=0.1, label='locdata_blobs')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local densities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We determine the local density for each localization and plot the probability density function for local densities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lc.LocalDensity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld_csr = lc.LocalDensity(radii=[3, 10, 30])\n",
    "ld_csr.compute(locdata_csr)\n",
    "ld_csr.results.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld_blob = lc.LocalDensity(radii=[3, 10, 30])\n",
    "ld_blob.compute(locdata_blob)\n",
    "ld_blob.results.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "radius = 3\n",
    "fig, ax = plt.subplots(nrows=1, ncols=2)\n",
    "locdata_csr.data.plot.scatter(x='position_x', y='position_y', ax=ax[0], color=ld_csr.results[radius], s=1, colormap='viridis', alpha=1, label='locdata_csr')\n",
    "locdata_blob.data.plot.scatter(x='position_x', y='position_y', ax=ax[1], color=ld_blob.results[radius], s=1, colormap='viridis', alpha=1, label='locdata_blobs')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "radius = 10\n",
    "fig, ax = plt.subplots(nrows=1, ncols=2)\n",
    "locdata_csr.data.plot.scatter(x='position_x', y='position_y', ax=ax[0], color=ld_csr.results[radius], s=1, colormap='viridis', alpha=1, label='locdata_csr')\n",
    "locdata_blob.data.plot.scatter(x='position_x', y='position_y', ax=ax[1], color=ld_blob.results[radius], s=1, colormap='viridis', alpha=1, label='locdata_blobs')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "radius = 30\n",
    "fig, ax = plt.subplots(nrows=1, ncols=2)\n",
    "locdata_csr.data.plot.scatter(x='position_x', y='position_y', ax=ax[0], color=ld_csr.results[radius], s=1, colormap='viridis', alpha=1, label='locdata_csr')\n",
    "locdata_blob.data.plot.scatter(x='position_x', y='position_y', ax=ax[1], color=ld_blob.results[radius], s=1, colormap='viridis', alpha=1, label='locdata_blobs')\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld_csr.hist(alpha=0.5, density=True, bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld_blob.hist(alpha=0.5, density=True, bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.arange(0, 1, 0.05)\n",
    "ld_csr.hist(alpha=0.5, bins=bins)\n",
    "ld_blob.hist(alpha=0.5, bins=bins);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Local densities with boundary correction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We determine the local density for each localization with a boundary correction applied. For the correction local density values are normalized by the relative overlapp of the encircling region and the support region. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "locdata = lc.simulate_uniform(n_samples=5_000, region=lc.Ellipse((0, 0), 100, 100))\n",
    "locdata.region"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "radius = 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld = lc.LocalDensity(radii=[radius]).compute(locdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld.hist(bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld.results.index = locdata.data.index\n",
    "df = locdata.dataframe.assign(local_density=ld.results[radius])\n",
    "locdata = locdata.update(dataframe=df)\n",
    "locdata.data.local_density.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "lc.render_2d_mpl(locdata, ax=ax, loc_properties=[\"position_x\", \"position_y\"], bin_size=1, other_property='local_density', rescale=lc.Trafo.NONE, vmin=0.2, vmax=0.7);\n",
    "if locdata.region:\n",
    "    locdata.region.plot(ax=ax, fill=False, color='Black');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld_2 = lc.LocalDensity(radii=[radius], boundary_correction=True).compute(locdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld_2.hist(bins=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ld_2.results.index = locdata.data.index\n",
    "df = locdata.dataframe.assign(local_density=ld_2.results[radius])\n",
    "locdata = locdata.update(dataframe=df)\n",
    "locdata.data.local_density.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "lc.render_2d_mpl(locdata, ax=ax, loc_properties=[\"position_x\", \"position_y\"], bin_size=1, other_property='local_density', rescale=lc.Trafo.NONE, vmin=0.2, vmax=0.7);\n",
    "if locdata.region:\n",
    "    locdata.region.plot(ax=ax, fill=False, color='White');"
   ]
  }
 ],
 "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
}