notebooks/tutorial_structure_factor.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This Notebook demonstrates the package ``structure_factor`` ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "start_all = time.time()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In most examples, we use a sample from the Ginibre point process as a running example pre-saved and loaded from the module \"data\".\n",
    "It is convenient, since its pair correlation function, structure factor, and hyperuniformity class are analytically known.\n",
    "Exceptionally, we used the Poisson point process with the multiscale hyperuniformity test since this test needs more than one sample of the underlying process and the Ginibre ensemble cannot be sampled on arbitrarily large windows on a personal computer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1- Installation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the package is not already installed/updated on your local machine you can simply use the following first-line/second for the installation/update."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install structure-factor\n",
    "# !pip install structure-factor --upgrade"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2- Create a PointPattern object "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create an object of type ``PointPattern`` from a point process, we need a sample of points, the observation window (typically a ``BallWindow`` or  a ``BoxWindow``), and the intensity of the point process (optional). \n",
    "A ``PointPattern`` from the Poisson, Thomas, Ginibre, or KLY point process could be created using the module ``point_process``. \n",
    "For more details, we refer to the documentation of this project https://for-a-few-dpps-more.github.io/structure-factor/.\n",
    "In what follows we will be loading a Ginibre PointPattern of 10^4 points (instead of creating one) using ``load_data`` from the module ``data``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Ginibre PointPattern \n",
    "from structure_factor.point_processes import GinibrePointProcess\n",
    "from structure_factor.data import load_data\n",
    "\n",
    "# Load\n",
    "point_pattern_ball = load_data.load_ginibre()\n",
    "# Plot\n",
    "point_pattern_ball.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Or create a Ginibre PointPattern as follows\n",
    "\n",
    "#from structure_factor.point_pattern import PointPattern\n",
    "#from structure_factor.point_processes import GinibrePointProcess\n",
    "#from structure_factor.spatial_windows import BallWindow\n",
    "\n",
    "#point_process = GinibrePointProcess()\n",
    "\n",
    "#window = BallWindow(center=[0, 0], radius=20)\n",
    "#points = point_process.generate_sample(window=window)\n",
    "#point_pattern_ball = PointPattern(points=points, window=window, intensity=point_process.intensity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exact structure factor and the pair correlation function of the Ginibre point process are known so we will use them to compare with the estimated values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from structure_factor.point_processes import GinibrePointProcess\n",
    "\n",
    "point_process = GinibrePointProcess()\n",
    "exact_sf = point_process.structure_factor\n",
    "exact_pcf = point_process.pair_correlation_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3- Approximate the structure factor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually, the package gathers estimators of the structure factor for:\n",
    "\n",
    "        1- stationary point processes, \n",
    "        2- stationary and isotropic point processes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1- Estimators assuming stationarity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimators that assume the stationarity of the point process, only apply to a ``PointPattern`` with an observation ``BoxWindow``.\n",
    "As ``point_pattern_ball`` has a ``BallWindow`` we will first restrict it to a new ``PointPattern`` with ``BoxWindow``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Restrict to a BoxWindow\n",
    "from structure_factor.spatial_windows import BoxWindow\n",
    "import numpy as np\n",
    "\n",
    "l = point_pattern_ball.window.radius*2/np.sqrt(2) # length side of the BoxWindow\n",
    "box_window = BoxWindow(bounds=[[-l/2, l/2], [-l/2, l/2]]) # BoxWindow\n",
    "point_pattern_box = point_pattern_ball.restrict_to_window(box_window) # Restrict to box_window\n",
    "# Plot\n",
    "point_pattern_box.plot()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we have the adequate ``PointPattern`` we will initialize StructureFactor and then use the estimators progressively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from structure_factor.structure_factor import StructureFactor\n",
    "\n",
    "sf = StructureFactor(point_pattern_box) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.1- Scattering intensity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The scattering intensity $\\widehat{S}_{\\mathrm{SI}}$ is an estimator of the structure factor $S$ of a stationary point process.\n",
    "It can be evaluated on arbitrary wavevectors, or allowed wavevectors to reduce the bias.\n",
    "It has also 2 debiased versions directly/undirectly.\n",
    "It is a particular case of the tapered estimator $\\widehat{S}_{\\mathrm{T}}$ using Bartlett's taper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scattering intensity on arbitrary wavevectors\n",
    "import numpy as np\n",
    "\n",
    "# Construct a grid of wavevectors\n",
    "k_max=8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(0, k_max, 200)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel())) # Wavevectors\n",
    "\n",
    "# Scattering intensity on k \n",
    "k, s_si_k = sf.scattering_intensity(k=k, debiased=False )\n",
    "\n",
    "# Plot\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axis = plt.subplots(figsize=(7,6))\n",
    "sf.plot_non_isotropic_estimator(k, s_si_k, \n",
    "                                axes=axis,\n",
    "                                plot_type=\"radial\", \n",
    "                                exact_sf=exact_sf,\n",
    "                                error_bar=True, bins=60,\n",
    "                                scale=\"log\", \n",
    "                                rasterized=True, \n",
    "                                label=r\"$\\widehat{S}$\", \n",
    "                                # file_name=\"s_si.pdf\" \n",
    "                                      )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The asymptotic bias of $\\widehat{S}_{\\mathrm{SI}}$ at small non allowed wavenumber $k$ is visible.\n",
    "To deal with this bias we can use one of the following three debiased versions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scattering intensity on allowed wavevectors\n",
    "\n",
    "k_max=8 # Threshold on the maximum wavenumber \n",
    "allowed_k, s_si_allowed_k = sf.scattering_intensity(k_max=k_max)\n",
    "\n",
    "# Plot\n",
    "sf.plot_non_isotropic_estimator(allowed_k, s_si_allowed_k, \n",
    "                                       plot_type=\"all\",   \n",
    "                                       exact_sf=exact_sf, \n",
    "                                       error_bar=True, bins=60,\n",
    "                                       label=r\"$\\widehat{S}$\",\n",
    "                                      # file_name=\"ginibre_s_si.pdf\" \n",
    "                                      )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Scattering intensity directly debiased on arbitrary wavevectors\n",
    "\n",
    "# Grid of wavevectors\n",
    "k_max=8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(-k_max, k_max, 17**2)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel())) # Wavevectors\n",
    "# Scattering intensity  \n",
    "k, s_si_dd = sf.scattering_intensity(k=k, debiased=True,  direct=True)\n",
    "\n",
    "# Plot\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sf.plot_non_isotropic_estimator(k, s_si_dd,  \n",
    "                                       plot_type=\"all\", \n",
    "                                       exact_sf=exact_sf, \n",
    "                                       error_bar=True, bins=80, \n",
    "                                       scale=\"log\",\n",
    "                                       rasterized=True, \n",
    "                                       label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Scattering intensity undirectly debiased on arbitrary wavevectors\n",
    "\n",
    "# Grid of wavevectors\n",
    "k_max=8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(0, k_max, 17**2)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel())) # Wavevectors\n",
    "# Scattering intensity \n",
    "k, s_si_ud = sf.scattering_intensity(k=k, debiased=True, direct=False)\n",
    "\n",
    "# Plot\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "sf.plot_non_isotropic_estimator(k, s_si_ud, \n",
    "                                #positive=True, \n",
    "                                plot_type=\"radial\",\n",
    "                                exact_sf=exact_sf, \n",
    "                                error_bar=True, bins=80, \n",
    "                                scale=\"log\", \n",
    "                                rasterized=True,\n",
    "                                label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see a few negative values, resulting in large inaccuracies in our log-log scale.\n",
    "You can add ``positive=True`` to the plot function to consider only positive values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.2- Tapered estimator with a single sine taper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tapered estimator $\\widehat{S}_{\\mathrm{T}}$ is an estimator of the structure factor $S$ of a stationary point process depending on a taper function.\n",
    "It has 2 debiased versions $\\widehat{S}_{\\mathrm{DDT}}$, and $\\widehat{S}_{\\mathrm{UDT}}$ (directly/undirectly).\n",
    "In what follows, we will be using the first taper of the family of sinusoidal tapers available in the module ``tapers``.\n",
    "To use a new taper follow the example in the documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tapered estimator\n",
    "from structure_factor.tapers import SineTaper\n",
    " \n",
    "k_max = 8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(0, k_max, 250)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel()))\n",
    "\n",
    "# First taper of the family of sinusoidal tapers\n",
    "p = [1, 1]\n",
    "tapers = [SineTaper(p)]\n",
    "\n",
    "# Tapered estimator\n",
    "k, s_t = sf.tapered_estimator(k, tapers=tapers, debiased=False)\n",
    "\n",
    "# Plot\n",
    "sf.plot_non_isotropic_estimator(k, s_t, \n",
    "                                plot_type=\"radial\", \n",
    "                                exact_sf=exact_sf, \n",
    "                                error_bar=True, bins=80, \n",
    "                                scale=\"log\",\n",
    "                                rasterized=True, \n",
    "                                label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The asymptotic bias of $\\widehat{S}_{\\mathrm{T}}$ at small wavenumbers $\\|k\\|_2$ is visible. To deal with this bias we can use one of the following two debiased versions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directly debiased tapered estimator\n",
    "\n",
    "from structure_factor.tapers import SineTaper\n",
    " \n",
    "k_max = 8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(0, k_max, 250)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel()))\n",
    "\n",
    "# First taper of the family of sinusoidal tapers\n",
    "p = [1, 1]\n",
    "tapers = [SineTaper(p)]\n",
    "\n",
    "# Tapered estimator\n",
    "k, s_ddt = sf.tapered_estimator(k, tapers=tapers, debiased=True, direct=True)\n",
    "\n",
    "# Plot\n",
    "sf.plot_non_isotropic_estimator(k, s_ddt, \n",
    "                                plot_type=\"radial\", \n",
    "                                exact_sf=exact_sf, \n",
    "                                error_bar=True, bins=80, \n",
    "                                scale=\"log\", \n",
    "                                rasterized=True, \n",
    "                                label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Undirectly debiased tapered estimator\n",
    "from structure_factor.tapers import SineTaper\n",
    "\n",
    "# Wavevectors\n",
    "k_max = 8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(0, k_max, 250)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel()))\n",
    "\n",
    "# First taper of the family of sinusoidal tapers\n",
    "p = [1, 1]\n",
    "tapers = [SineTaper(p)]\n",
    "\n",
    "# Tapered estimator\n",
    "k, s_udt = sf.tapered_estimator(k, tapers=tapers, debiased=True, direct=False)\n",
    "\n",
    "# Plot\n",
    "sf.plot_non_isotropic_estimator(k, s_udt, \n",
    "                                positive=True,\n",
    "                                plot_type=\"radial\", \n",
    "                                exact_sf=exact_sf, \n",
    "                                error_bar=True, bins=80, \n",
    "                                scale=\"log\",\n",
    "                                rasterized=True, \n",
    "                                label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1.3- Multitapered estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The multitapered estimator $\\widehat{S}_{\\mathrm{MT}}$ is an estimator of the structure factor $S$ of a stationary point process depending on a family of tapers. \n",
    "It is a generalization of the tapered estimator $\\widehat{S}_{\\mathrm{T}}$.\n",
    "It has 2 debiased versions $\\widehat{S}_{\\mathrm{DDMT}}$, and $\\widehat{S}_{\\mathrm{UDMT}}$ (directly/undirectly).\n",
    "In what follows, we will be using the first 4 tapers of the family of sinusoidal tapers available by the module ``tapers``.\n",
    "To use a new family of tapers, follow the example present in the documentation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multitapered estimator\n",
    "from structure_factor.tapers import SineTaper\n",
    " \n",
    "# Wavevectors\n",
    "k_max = 8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(0, k_max, 200)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel()))\n",
    "\n",
    "# First 4 tapers of the family of sinusoidal tapers\n",
    "from structure_factor.tapers import multi_sinetaper_grid\n",
    "tapers = multi_sinetaper_grid(d=2, p_component_max=2)\n",
    "\n",
    "# Scaled multitapered periodogram\n",
    "k, s_mt = sf.tapered_estimator(k, tapers=tapers, debiased=False)\n",
    "\n",
    "# Plot\n",
    "sf.plot_non_isotropic_estimator(k, s_mt, \n",
    "                                plot_type=\"radial\", \n",
    "                                exact_sf=exact_sf, \n",
    "                                error_bar=True, bins=80,\n",
    "                                scale=\"log\",\n",
    "                                rasterized=True, \n",
    "                                label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The asymptotic bias of $\\widehat{S}_{\\mathrm{MT}}$ at small wavenumber $\\|k\\|_2$ is visible. \n",
    "To deal with this bias we can use one of the following two debiased versions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Directly debiased multitapered estimator\n",
    "from structure_factor.tapers import SineTaper\n",
    " \n",
    "# Wavevectors\n",
    "k_max = 8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(-k_max, k_max, 16**2)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel()))\n",
    "\n",
    "# First 4 tapers of the family of sinusoidal tapers\n",
    "from structure_factor.tapers import multi_sinetaper_grid\n",
    "tapers = multi_sinetaper_grid(d=2, p_component_max=2)\n",
    "\n",
    "# Scaled multitapered periodogram\n",
    "k, s_ddmt = sf.tapered_estimator(k, tapers=tapers, debiased=True, direct=True)\n",
    "\n",
    "# Plot\n",
    "sf.plot_non_isotropic_estimator(k, s_ddmt, \n",
    "                                plot_type=\"all\", \n",
    "                                exact_sf=exact_sf,  \n",
    "                                error_bar=True, bins=80, \n",
    "                                scale=\"log\", \n",
    "                                rasterized=True,\n",
    "                                label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Undirectly debiased multitapered estimator\n",
    "from structure_factor.tapers import SineTaper\n",
    " \n",
    "# Wavevectors\n",
    "k_max = 8 # Threshold on the maximum wavenumber\n",
    "x = np.linspace(-k_max, k_max, 16**2)\n",
    "x = x[x != 0]\n",
    "X, Y = np.meshgrid(x, x)\n",
    "k = np.column_stack((X.ravel(), Y.ravel()))\n",
    "\n",
    "# First 4 tapers of the family of sinusoidal tapers\n",
    "from structure_factor.tapers import multi_sinetaper_grid\n",
    "tapers = multi_sinetaper_grid(d=2, p_component_max=2)\n",
    "\n",
    "# Scaled multitapered periodogram\n",
    "k, s_udmt = sf.tapered_estimator(k, tapers=tapers, debiased=True, direct=False)\n",
    "\n",
    "# Plot\n",
    "sf.plot_non_isotropic_estimator(k, s_ddmt, \n",
    "                                plot_type=\"imshow\", \n",
    "                                exact_sf=exact_sf, \n",
    "                                error_bar=True, bins=80,\n",
    "                                scale=\"log\", \n",
    "                                rasterized=True, \n",
    "                                label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2- Estimators assuming stationarity and isotropy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimators that assume the stationarity and isotropy of the point process, apply to a ``PointPattern`` having ``BallWindow``.\n",
    "As ``point_pattern_ball`` has a ``BallWindow`` it will be used directly. \n",
    "In the other case, it's recommended to use a restriction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sf_isotropic = StructureFactor(point_pattern_ball)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2.1- Bartlett's isotropic estimator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bartlett's isotropic estimator $\\widehat{S}_{\\mathrm{BI}}$ is an estimator of the structure factor $S$ of a stationary isotropic point process. \n",
    "It is a particular case of the isotropic tapered estimators using Bartlett's taper.\n",
    "In other words, it's the isotropic version of the scattering intensity. \n",
    "It can be evaluated on arbitrary wavenumbers or the allowed wavenumbers to reduce the bias.\n",
    "The complexity of this estimator is quadratic in the number of points of the ``PointPattern``.\n",
    "When dealing with a large number of points, it is recommended to start with a restricted version of the ``PointPattern``, before using the large sample and to increase the window radius respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bartlett's isotropic estimator on arbitrary wavenumbers\n",
    "# Running time = 156 s\n",
    "import time\n",
    "start_time = time.time()\n",
    "k_norm = np.linspace(0.01, 7, 100)\n",
    "k_norm, s_bi = sf_isotropic.bartlett_isotropic_estimator(k_norm=k_norm)\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "\n",
    "#Plot\n",
    "fig, axis = plt.subplots(figsize=(7,6))\n",
    "sf_isotropic.plot_isotropic_estimator(k_norm, s_bi, \n",
    "                                         axis=axis,\n",
    "                                         exact_sf=exact_sf, \n",
    "                                         label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the relevant bias for small wavenumbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Bartlett's isotropic estimator on allowed wavenumbers\n",
    "# Running time = 156 s\n",
    "import time\n",
    "\n",
    "start_time = time.time()\n",
    "k_norm_allowed, s_bi_k_norm_allowed = sf_isotropic.bartlett_isotropic_estimator(nb_values=60)\n",
    "print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "\n",
    "#Plot\n",
    "fig, axis = plt.subplots(figsize=(7,6))\n",
    "sf_isotropic.plot_isotropic_estimator(k_norm_allowed, s_bi_k_norm_allowed, \n",
    "                                         axis=axis,\n",
    "                                         exact_sf=exact_sf, \n",
    "                                         label=r\"$\\widehat{S}$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This estimator seems to give the best results among all the estimators but, it's time-consuming."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2.2- Quadrature estimator isotropic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Warning**: This section relies on an interface between Python and R. Make sure that the R Programming language is available on your machine. For installation, see https://www.r-project.org/\n",
    "\n",
    "The last two isotropic estimators are based on, first estimating the pair correlation function $g$ of the ``PointPattern``,\n",
    "then interpolating/extrapolating the obtained results, and finally using Ogata/Baddour&Chouinard quadrature to estimate the Hankel transform of a function of g, corresponding to the structure factor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by estimating the pair correlation function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.2.1 - Pair correlation function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Warning**: This section requires the R program language https://www.rstudio.com/ to be available on your local machine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The toolbox contains two estimators of the pair correlation function (pcf) of a ``PointPattern``.\n",
    "These two methods are ``pcf.ppp`` and ``pcf.fv`` inherited from the R package ``spatstat`` https://spatstat.org/ ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Approximating the pair correlation function using pcf.ppp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pcf.ppp\n",
    "import structure_factor.pair_correlation_function as pcf\n",
    "\n",
    "r= np.linspace(0, 20, 200)\n",
    "pcf_ppp = pcf.estimate(point_pattern_ball, # PointPattern\n",
    "                        method=\"ppp\", # Estimation method (could be \"ppp\" or \"fv\")\n",
    "                        r=r, # Estimation radii \n",
    "                        correction=\"all\" # Edge correction \n",
    "                              )\n",
    "\n",
    "# Plot\n",
    "pcf.plot(pcf_ppp, \n",
    "         exact_pcf=exact_pcf, \n",
    "         figsize=(7,6),\n",
    "         color=['grey', 'b', 'darkcyan'], \n",
    "         style=[\".\", \"*\", \"^\"], \n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Approximating the pair correlation function using pcf.fv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pcf.fv\n",
    "import structure_factor.pair_correlation_function as pcf\n",
    "\n",
    "pcf_fv = pcf.estimate(point_pattern_ball, # PointPattern \n",
    "                      method=\"fv\", # Estimation method (could be \"ppp\" or \"fv\")\n",
    "                      Kest=dict(rmax=20), # Maximal estimation radius \n",
    "                      fv=dict(method=\"b\", spar=0.2) # Correction method and sparsity\n",
    "                             )\n",
    "\n",
    "# Plot\n",
    "pcf.plot(pcf_fv, \n",
    "         exact_pcf=exact_pcf, \n",
    "         figsize=(7,6), \n",
    "         color=['grey'], \n",
    "         style=[\".\"]\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After Approximating the pair correlation function the second step is to clean, interpolate, and extrapolate the result.\n",
    "We will be working with the approximation result of the method ``pcf.fv``.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcf_fv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = pcf_fv[\"r\"]\n",
    "pcf_r = pcf_fv[\"pcf\"]\n",
    "\n",
    "# Interpolation\n",
    "pcf_fct = pcf.interpolate(r=r, \n",
    "                          pcf_r=pcf_r, \n",
    "                          drop=True, # Drop outliers (nan, neginf, and posinf)\n",
    "                          extrapolate_with_one=True # Extrapolate with 1 after max(r)\n",
    "                                 )\n",
    "\n",
    "\n",
    "# Plot\n",
    "x = np.linspace(0, 100, 200)\n",
    "plt.plot(x, pcf_fct(x), 'b.', label=\"Approximated pcf\") \n",
    "plt.plot(x, exact_pcf(x), 'g', label=\"Exact pcf\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.2.2- Estimate the structure factor using Ogata quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ogata quadrature \n",
    "k_norm = np.linspace(0.5, 30, 1000) # Wavenumbers\n",
    "k_norm, s_ho = sf_isotropic.quadrature_estimator_isotropic(pcf_fct, # Estimated pcf function\n",
    "                                                 k_norm=k_norm, # wavenumbers\n",
    "                                                 method=\"Ogata\", # Quadrature method (\"Ogata\", or \"BaddourChouinard\")\n",
    "                                                 step_size=0.01, # Step size\n",
    "                                                 nb_points=1000 # Number of points\n",
    "                                                )\n",
    "# Plot\n",
    "fig, axis = plt.subplots(figsize=(7,6))\n",
    "fig = sf_isotropic.plot_isotropic_estimator(k_norm, s_ho, \n",
    "                                               axis=axis, \n",
    "                                               error_bar=True, bins=80, \n",
    "                                               exact_sf=exact_sf, \n",
    "                                               label=r\"$\\widehat{S}$\",\n",
    "                                               #file_name=\"ginibre_s_ho.pdf\" \n",
    "                                              )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3.2.2.2- Estimate the structure factor using Baddour & Chouinard quadrature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Baddour & Chouinard quadrature \n",
    "k_norm = np.linspace(0.5, 30, 1000) # Wavenumbers\n",
    "r_max = np.max(r) # Maximum radius for which the pcf has been approximated\n",
    "k_norm, s_hbc = sf_isotropic.quadrature_estimator_isotropic(pcf_fct, # Estimated pcf function\n",
    "                                                 k_norm=k_norm, # wavenumbers\n",
    "                                                 method=\"BaddourChouinard\", # Quadrature method (\"Ogata\", or \"BaddourChouinard\")\n",
    "                                                 r_max=r_max, # Step size\n",
    "                                                 nb_points=1000 # Number of points\n",
    "                                                )\n",
    "# Plot\n",
    "fig, axis = plt.subplots(figsize=(7,6))\n",
    "fig = sf_isotropic.plot_isotropic_estimator(k_norm, s_hbc, \n",
    "                                            axis=axis, \n",
    "                                            exact_sf=exact_sf, \n",
    "                                            label=r\"$\\widehat{S}$\",\n",
    "                                            #file_name=\"ginibre_s_hbc.pdf\" \n",
    "                                              )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.4- Hyperuniformity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The trade-off for studying the hyperuniformity of a point process is by testing the effective hyperuniformity.\n",
    "As the definition of effective hyperuniformity is not standardized, we have contributed a new asymptotic statistical test of hyperuniformity, called the multiscale test.\n",
    "Users can decide which test to be used depending on their study purpose.\n",
    "\n",
    "If the test of hyperuniformity suggests that the point process is hyperuniform, we can study the power decay to zero of the structure factor $\\alpha$ which indicates the possible class of hyperuniformity.\n",
    "If the multiscale test rejects the hypothesis of hyperuniformity then we conclude that the ``PointPattern`` doesn't exhibit a hyperuniform structure.\n",
    "Hence, either the point process is not hyperuniform or the observation window is not big enough to cover hyperuniformity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Ginibre ensemble is a first-class hyperuniform with power decay $\\alpha_{Ginibre} = 2$, and the Poisson point process is not hyperunifrom.\n",
    "In what follows, we will be using the results of Bartlett's isotropic estimator for studying the effective hyperuniformity and the hyperuniformity class of the Ginibre ``PointPattern``.\n",
    "\n",
    "\n",
    "As the multiscale hyperuniformity test uses a list of ``PointPattern`` and sampling from the Ginibre ensemble is costly on a personal computer, we will use the Poisson point process instead, with the scattering intensity estimator to benchmark the multiscale hyperuniformity test."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.1- Effective hyperuniformity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index $H$ is used in the physics literature to study the effective hyperuniformity of a point process.\n",
    "If $H< 10^{-3}$ then the point process is effectively hyperuniform and so maybe hyperuniform. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from structure_factor.hyperuniformity import effective_hyperuniformity\n",
    "# Effective hyperuniformity\n",
    "summary_H = effective_hyperuniformity(k_norm_allowed, s_bi_k_norm_allowed,\n",
    "                                         k_norm_stop=0.2) # H index\n",
    "\n",
    "# Visualization of the results\n",
    "import matplotlib.pyplot as plt\n",
    "import structure_factor.utils as utils\n",
    "fitted_line = summary_H[\"fitted_line\"] # Fitted line to s_bi\n",
    "x = np.linspace(0, 2, 300)\n",
    "fig, axis =plt.subplots(figsize=(7,5))\n",
    "axis.plot(k_norm_allowed, s_bi_k_norm_allowed, 'b', marker=\".\", label=\"Approximated structure factor\")\n",
    "axis.plot(x, fitted_line(x), 'r--', label= \"Fitted line\")\n",
    "axis.plot(k_norm_allowed, GinibrePointProcess.structure_factor(k_norm_allowed), 'g', label=r\"$S(k)$\")\n",
    "axis.annotate('H={}'.format(summary_H[\"H\"]), xy=(0, 0), xytext=(0.01,0.3),\n",
    "            arrowprops=dict(facecolor='black', shrink=0.01))\n",
    "axis.legend()\n",
    "axis.set_xlabel('wavelength (k)')\n",
    "axis.set_ylabel(r\"Structure factor ($\\mathsf{S}(k)$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As $H<10^{-3}$, we conclude that the Ginibre ensemble is effectively hyperuniformity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.2- Multiscale hyperuniformity test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A point process is hyperuniform iff $\\mathbb{E}[Z]=0$.\n",
    "The multiscale test is an asymptotic hyperuniformity test based on testing whether zero belongs to the confidence interval of $\\mathbb{E}[Z]$ or not.\n",
    "In what follows, we chose a confidence level of $99.7 \\%$, i.e., $3$ standard deviations to test the hyperuniformity of a Poisson point process of intensity $1/(2\\pi)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from structure_factor.point_processes import HomogeneousPoissonPointProcess\n",
    "from structure_factor.spatial_windows import BallWindow\n",
    "from structure_factor.hyperuniformity import subwindows, multiscale_test\n",
    "\n",
    "# PointPattern list\n",
    "poisson = HomogeneousPoissonPointProcess(intensity=1 / (2*np.pi))\n",
    "N = 100\n",
    "R = 80\n",
    "window = BallWindow(center=[0,0], radius=R)\n",
    "point_patterns = [\n",
    "    poisson.generate_point_pattern(window=window) for _ in range(N)\n",
    "]\n",
    "# subwindows and k\n",
    "r_0 = 5\n",
    "subwindows_list, k = subwindows(window, subwindows_type=\"BallWindow\", param_0=r_0)\n",
    "\n",
    "# multiscale hyperuniformity test\n",
    "mean_poisson = 60\n",
    "summary = multiscale_test(point_patterns,\n",
    "    estimator=\"bartlett_isotropic_estimator\", k_list=k,\n",
    "    subwindows_list=subwindows_list, mean_poisson=mean_poisson,\n",
    "    verbose=True)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.hist(summary[\"Z\"], color=\"grey\", label=\"Z\")\n",
    "plt.axvline(summary[\"mean_Z\"], color=\"b\", linewidth=1, label=r\"$\\bar{Z}$\")\n",
    "plt.axvline(summary[\"mean_Z\"] - 3*summary[\"std_mean_Z\"], linestyle='dashed',\n",
    "            color=\"k\", linewidth=1, label=r\"$\\bar{Z} \\pm 3 \\bar{\\sigma}/ \\sqrt{N}$\")\n",
    "plt.axvline(summary[\"mean_Z\"] + 3*summary[\"std_mean_Z\"], linestyle='dashed',\n",
    "            color=\"k\", linewidth=1)\n",
    "plt.xlabel(r\"$Z$\")\n",
    "plt.title(r\"Histogram of the $N$ obtained values of $Z$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(summary[\"mean_Z\"] - 3*summary[\"std_mean_Z\"], summary[\"mean_Z\"] + 3*summary[\"std_mean_Z\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, zero is not in the confidence interval of $\\mathbb{E}[Z]$, hinting that the Poisson point process is not hyperuniform."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3.2- Hyperuniformity class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from structure_factor.hyperuniformity import hyperuniformity_class\n",
    "# Hyperuniformity class\n",
    "summary_alpha = hyperuniformity_class(k_norm_allowed, s_bi_k_norm_allowed,\n",
    "                                         k_norm_stop=0.4)\n",
    "\n",
    "# Visualization of the results\n",
    "import matplotlib.pyplot as plt\n",
    "import structure_factor.utils as utils\n",
    "fitted_poly = summary_alpha[\"fitted_poly\"] # Fitted polynomial to s_bi\n",
    "fig, axis =plt.subplots(figsize=(7,5))\n",
    "axis.plot(k_norm_allowed, s_bi_k_norm_allowed, 'b', marker=\".\", label=\"Approximated structure factor\")\n",
    "axis.plot(k_norm_allowed, GinibrePointProcess.structure_factor(k_norm_allowed), 'g', label=r\"$S(k)$\")\n",
    "axis.plot(k_norm_allowed, fitted_poly(k_norm_allowed), 'r--', label= \"Fitted line\")\n",
    "axis.annotate(r\" $\\alpha$ ={}\".format(summary_alpha[\"alpha\"]), xy=(0, 0), xytext=(0.01,0.4),\n",
    "            arrowprops=dict(facecolor='black', shrink=0.1))\n",
    "axis.legend()\n",
    "axis.set_xlabel('wavelength (k)')\n",
    "axis.set_ylabel(r\"Structure factor ($\\mathsf{S}(k)$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\alpha \\approx 1.99$ is a good approximation of the exact value $\\alpha_{Ginibre}=2$.\n",
    "Thus this test successfully predicts the hyperuniformity class of the Ginibre Ensemble."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"--- %s seconds ---\" % (time.time() - start_all))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "55b4d108f296a679a2dff1927b0e4c2a151de922cfd09e10911d1814a00e3c5c"
  },
  "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}