{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ratio Simulations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import scipy\n",
    "import skimage\n",
    "import skimage.segmentation\n",
    "import skimage.transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson distribution for 3 channels\n",
    "\n",
    "This simulation compares the arithmetic and geometric averages of three signals, each following a Poisson distribution. \n",
    "\n",
    "While the arithmetic average exhibits a slightly superior reduction in standard deviation, it's noteworthy that the distribution of the arithmetic average tends to be more discernibly separated from the background when signals of different intensities are considered.\n",
    "\n",
    "This insight lays the foundation for more robust segmentation methodologies tailored to multichannel images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "N = 1000000\n",
    "\n",
    "rng = np.random.default_rng(1)\n",
    "\n",
    "s1 = rng.poisson(20, N)\n",
    "s2 = rng.poisson(40, N)\n",
    "s3 = rng.poisson(90, N)\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.hist(s1, histtype=\"step\", bins=20, lw=3, alpha=0.5)\n",
    "plt.hist(s2, histtype=\"step\", bins=20, lw=3, alpha=0.5)\n",
    "plt.hist(s3, histtype=\"step\", bins=20, lw=3, alpha=0.5)\n",
    "plt.hist(\n",
    "    np.power(s1 * s2 * s3, 1 / 3),\n",
    "    histtype=\"step\",\n",
    "    bins=20,\n",
    "    lw=2,\n",
    "    label=\"Arithmetic Mean\",\n",
    ")\n",
    "plt.hist(\n",
    "    np.average([s1, s2, s3], axis=0),\n",
    "    histtype=\"step\",\n",
    "    bins=20,\n",
    "    lw=2,\n",
    "    label=\"Geometric Mean\",\n",
    ")\n",
    "plt.xlabel(\"Signal Intensity\")\n",
    "plt.ylabel(\"Counts\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "\n",
    "# Compute standard deviations\n",
    "s_avg_std = np.std(np.average([s1, s2, s3], axis=0))\n",
    "s_cbrt_std = np.std(np.power(s1 * s2 * s3, 1 / 3))\n",
    "# Print standard deviations\n",
    "print(f\"Standard Deviation of Arithmetic Mean of Signals: {s_avg_std:.3g}\")\n",
    "print(f\"Standard Deviation of Geometric Mean of Signals: {s_cbrt_std:.3g}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ratio calculation\n",
    "\n",
    "We move on exploring the simulation of two channels (signals) and examining their ratios. Our goal is to understand how different averaging methods can influence the resulting distribution and statistical measures of these ratios. We'll also introduce the concept of sliding window analysis to simulate image Gaussian filtering and explore its impact on the ratio calculation.\n",
    "\n",
    "Let's begin by generating two signals, \"cyan\" and \"red,\" each following a Poisson distribution with different intensities. We'll plot histograms of their signal intensities to visualize their distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the number of pixels in the region of interest (ROI)\n",
    "N = 999\n",
    "\n",
    "# Generate Poisson-distributed signals for cyan and red channels\n",
    "cyan = rng.poisson(50, N)\n",
    "red = rng.poisson(50, N)\n",
    "\n",
    "# Plot histograms of signal intensities\n",
    "plt.figure(figsize=(10, 6))  # Set figure size\n",
    "plt.hist(\n",
    "    cyan, bins=30, histtype=\"step\", lw=3, alpha=0.6, label=\"cyan\", color=\"darkcyan\"\n",
    ")\n",
    "plt.hist(red, bins=30, histtype=\"step\", lw=3, alpha=0.6, label=\"red\", color=\"indianred\")\n",
    "plt.xlabel(\"Signal Intensity\")\n",
    "plt.ylabel(\"Counts\")\n",
    "plt.title(\"Histograms of Signal Intensity\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The histograms display the distribution of signal intensities for both the cyan and red channels.\n",
    "\n",
    "Next, we'll calculate the ratio between the cyan and red channels using various methods:\n",
    "\n",
    "1. Simple Ratio: We'll compute the ratio by dividing the cyan signal by the red signal.\n",
    "2. Running Ratio Median: We'll simulate image Mean filtering by applying a sliding window approach with a window size of 9 pixels over both channels and calculate the median of the ratios within each window.\n",
    "3. Ratiorank Filter: This kind of filter employs a sliding window approach. Within each window, the ratiorank filter calculates the ratio between every possible pair of values from the cyan and red signal channels. For a window size of 9 pixels per channel, this results in a total of 81 unique ratio calculations. Finally, aggregating these ratios using the median provides a robust estimate of the ratio central tendency, capturing localized variations in signal ratios."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1. Calculate the ratio of cyan-over-red\n",
    "ratio = cyan / red\n",
    "\n",
    "# 2. Calculate the ratio after running average the signals\n",
    "rr = np.convolve(cyan, np.ones(9) / 9, mode=\"valid\") / np.convolve(\n",
    "    red, np.ones(9) / 9, mode=\"valid\"\n",
    ")\n",
    "\n",
    "\n",
    "# 3. Calculate the ratio using RatioRank\n",
    "def rank(numerators, denominators):\n",
    "    r = [num / den for num in numerators for den in denominators]\n",
    "    return np.median(r)\n",
    "\n",
    "\n",
    "def sliding_window(arr1, arr2, size=9):\n",
    "    \"\"\"Generate sliding window of size `size` over two input arrays `arr1` and `arr2`.\"\"\"\n",
    "    for i in range(len(arr1) - size + 1):\n",
    "        yield arr1[i : i + size], arr2[i : i + size]\n",
    "\n",
    "\n",
    "ratiorank = []\n",
    "for window1, window2 in sliding_window(cyan, red):\n",
    "    ratiorank.append(rank(window1, window2))\n",
    "\n",
    "# Plot histogram of the ratio\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.hist(\n",
    "    ratio,\n",
    "    bins=30,\n",
    "    histtype=\"step\",\n",
    "    lw=3,\n",
    "    alpha=0.6,\n",
    "    color=\"purple\",\n",
    "    label=\"Simple Ratio\",\n",
    ")\n",
    "plt.hist(rr, bins=30, histtype=\"step\", lw=2, alpha=0.6, label=\"Running Average\")\n",
    "plt.hist(\n",
    "    ratiorank,\n",
    "    bins=30,\n",
    "    histtype=\"step\",\n",
    "    lw=2,\n",
    "    alpha=0.6,\n",
    "    color=\"blue\",\n",
    "    label=\"Median Ratiorank\",\n",
    ")\n",
    "\n",
    "plt.xlabel(\"Cyan / Red Ratio\")\n",
    "plt.ylabel(\"Counts\")\n",
    "plt.title(\"Histogram of Cyan-over-Red Ratio\")\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# Plot probability plots for simple ratio, ratiorank, and running average\n",
    "plt.figure(figsize=(10.5, 3.5))\n",
    "# Simple ratio\n",
    "plt.subplot(1, 3, 1)\n",
    "scipy.stats.probplot(ratio, plot=plt, rvalue=True)\n",
    "plt.title(\"Simple Ratio\")\n",
    "plt.xlabel(\"Theoretical Quantiles\")\n",
    "plt.ylabel(\"Ordered Values\")\n",
    "# Running average\n",
    "plt.subplot(1, 3, 2)\n",
    "scipy.stats.probplot(rr, plot=plt, rvalue=True)\n",
    "plt.title(\"Running Average\")\n",
    "plt.xlabel(\"Theoretical Quantiles\")\n",
    "plt.ylabel(\"Ordered Values\")\n",
    "# Ratiorank\n",
    "plt.subplot(1, 3, 3)\n",
    "scipy.stats.probplot(ratiorank, plot=plt, rvalue=True)\n",
    "plt.title(\"Ratiorank\")\n",
    "plt.xlabel(\"Theoretical Quantiles\")\n",
    "plt.ylabel(\"Ordered Values\")\n",
    "# Tight\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simple ratio values exhibit a skewed distribution towards higher values, evident from the deviation observed in the probplot and the notable differences in the tails apparent in the q-q plot. \n",
    "\n",
    "Indeed, ratiorank emerges as a promising alternative to both the simple ratio and running average methods in this context. By computing the median of 81 calculated ratios from sliding windows, ratiorank offers improved robustness and accuracy compared to the running average, which relies on fewer samples. This enhanced fidelity suggests ratiorank as a valuable approach for applications requiring precise ratio calculations,\n",
    "\n",
    "Hence, the ratio of the entire ROI is accurately estimated from the median of the calculated distribution, as the mean may be unreliable for non-normally distributed data.\n",
    "\n",
    "4. Additionally, signal averaging over the ROI prior to ratio calculation can be considered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 4. Calculate the ratio of averaged signals\n",
    "ratio_of_avg = np.average(cyan) / np.average(red)\n",
    "ratio_of_median = np.median(cyan) / np.median(red)\n",
    "\n",
    "# Print statistical measures of the ratios\n",
    "print(\n",
    "    f\"1. Simple Ratio (median vs. mean): {np.median(ratio):.3f} vs. {np.average(ratio):.3f}\"\n",
    ")\n",
    "print(\n",
    "    f\"1. Running Ratio (median vs. mean): {np.median(rr):.3f} vs. {np.average(rr):.3f}\"\n",
    ")\n",
    "print(\n",
    "    f\"1. Ratiorank (median vs. mean): {np.median(ratiorank):.3f} vs. {np.average(ratiorank):.3f}\"\n",
    ")\n",
    "print(f\"Ratio of (median vs.) Averages: ({ratio_of_median:.3f} vs.) {ratio_of_avg:.3f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define variables\n",
    "roi_size = 80\n",
    "signal = 10\n",
    "# Initialize lists to store results\n",
    "ratios_med = []\n",
    "ratios_avg = []\n",
    "rrs_med = []\n",
    "rrs_avg = []\n",
    "rranks_med = []\n",
    "rranks_avg = []\n",
    "r_med = []\n",
    "r_avg = []\n",
    "# Perform simulation\n",
    "for _ in range(199):\n",
    "    # Generate random signals for cyan and red channels\n",
    "    cyan = rng.poisson(signal, roi_size)\n",
    "    red = rng.poisson(signal, roi_size)\n",
    "    # Calculate simple ratios\n",
    "    ratio = cyan / red\n",
    "    # Calculate running averages\n",
    "    rr = np.convolve(cyan, np.ones(9) / 9, mode=\"valid\") / np.convolve(\n",
    "        red, np.ones(9) / 9, mode=\"valid\"\n",
    "    )\n",
    "    # Calculate ratiorank\n",
    "    ratiorank = []\n",
    "    for window1, window2 in sliding_window(cyan, red):\n",
    "        ratiorank.append(rank(window1, window2))\n",
    "    # Store results\n",
    "    ratios_med.append(np.median(ratio))\n",
    "    ratios_avg.append(np.mean(ratio))\n",
    "    rrs_med.append(np.median(rr))\n",
    "    rrs_avg.append(np.mean(rr))\n",
    "    rranks_med.append(np.median(ratiorank))\n",
    "    rranks_avg.append(np.mean(ratiorank))\n",
    "    r_med.append(np.median(cyan) / np.median(red))\n",
    "    r_avg.append(np.mean(cyan) / np.mean(red))\n",
    "\n",
    "# Create DataFrame to organize results and plot boxplot\n",
    "df = pd.DataFrame(\n",
    "    np.column_stack(\n",
    "        (r_avg, r_med, rranks_avg, rranks_med, rrs_avg, rrs_med, ratios_avg, ratios_med)\n",
    "    ),\n",
    "    columns=[\n",
    "        \"mean N / mean D\",\n",
    "        \"median N / median D\",\n",
    "        \"mean ratiorank\",\n",
    "        \"median ratiorank\",\n",
    "        \"mean rr\",\n",
    "        \"median rr\",\n",
    "        \"mean ratio\",\n",
    "        \"median ratio\",\n",
    "    ],\n",
    ")\n",
    "f = df.boxplot(vert=False, showfliers=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import gaussian_kde\n",
    "\n",
    "# Evaluate the KDE on a grid of points\n",
    "x_grid = np.linspace(cyan.min(), cyan.max(), 1000)\n",
    "y_kde = gaussian_kde(cyan)(x_grid)  # Evaluate the KDE at x_grid points\n",
    "\n",
    "# Plot the KDE smoothed line\n",
    "plt.plot(x_grid, y_kde)\n",
    "\n",
    "plt.hist(cyan, density=1, bins=10, histtype=\"step\", lw=2, alpha=0.6, label=\"Cyan\")\n",
    "\n",
    "# Find the peak position of the KDE\n",
    "peak_position = x_grid[np.argmax(y_kde)]\n",
    "peak_value = np.max(y_kde)\n",
    "\n",
    "# Plot the peak position\n",
    "plt.axvline(\n",
    "    x=peak_position,\n",
    "    color=\"r\",\n",
    "    linestyle=\"--\",\n",
    "    label=f\"Peak at {peak_position:.2f}, Value: {peak_value:.2f}\",\n",
    ")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an image.\n",
    "\n",
    "TODO: Complete with Running-ratio and Ratiorank."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "im_cyan = skimage.util.random_noise(25 * np.ones([64, 64]), mode=\"poisson\", clip=0)\n",
    "im_red = skimage.util.random_noise(25 * np.ones([64, 64]), mode=\"poisson\", clip=0)\n",
    "im_ratio = im_cyan / im_red\n",
    "\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.subplot(1, 2, 2)\n",
    "skimage.io.imshow(im_ratio)\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.hist(im_ratio.ravel(), histtype=\"step\", bins=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MAYBE/TODO:\n",
    "\n",
    "- simulare fret e bt:\n",
    "    - $S^{561} = R$\n",
    "    - $S^{488} = G (1 - E) + R (\\alpha + E)$\n",
    "    - $S^{458} = C (1 - E) + R (\\beta + E)$\n",
    "\n",
    "- test ch1, ch2, ch3 and geometric average\n",
    "- test cross-correlation\n",
    "- consider detector gain in simulated error\n",
    "- simulate 3D obj + PSF in generation of synthetic data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage import correlate as ic\n",
    "\n",
    "skimage.io.imshow(ic(im_cyan, im_red))"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}