Rename the Rosenblatt method; add the double CDF method. Add a simple example in notebook.

EtaInterpolationFunctions.py

@@ -5,7 +5,7 @@ Functions for eta interpolation
     2. bilinear_xy_lookup(etaX, etaY, X_tab, Y_tab, x_centers, y_centers)
     3. map_xy_no_interp(etaX, etaY, X_tab, Y_tab, x_edges, y_edges)
 '''
-def build_xy_lut_from_hist(H):
+def build_xy_lut_Rosenblatt(H):
     """
     Generate two 2D tables X_tab, Y_tab from a 2D histogram H
     such that (X_tab[i,j], Y_tab[i,j]) is the (x,y) value at the center of bin (i,j).
@@ -44,6 +44,44 @@ def build_xy_lut_from_hist(H):
 
     return (X_tab, Y_tab)
 
+def build_xy_lut_DoubleCDF(H):
+    """
+    Generate two 2D tables X_tab, Y_tab from a 2D histogram H
+    such that (X_tab[i,j], Y_tab[i,j]) is the (x,y) value at the center of bin (i,j).
+    Input:
+        H: 2D histogram of charge centroid (eta_x, eta_y) distribution, shape (nBinX, nBinY), ranging from 0 to 1.
+    Output:
+        X_tab, Y_tab: two 2D tables of interpolated (x,y) values in position space, shape (nBinX, nBinY), ranging from 0 to 1.
+    """
+    H = np.asarray(H, dtype=float)
+    nBinX, nBinY = H.shape
+    etaX_edges = np.linspace(0, 1, nBinX+1)
+    etaY_edges = np.linspace(0, 1, nBinY+1)
+
+    # slight smoothing to avoid zero-count rows/columns causing flat CDF
+    eps=1e-9
+    alpha = max(eps, 1e-12 * max(1.0, H.sum()))
+    Hs = H + alpha
+
+    # marginal 1D CDF F_EtaX
+    pEtaX = Hs.sum(axis=1)                                  # Marginal PDF (nBinX,)
+    F_EtaX_edges = np.concatenate([[0.0], np.cumsum(pEtaX)])# CDF at bin edges (nBinX+1,)
+    F_EtaX_edges = F_EtaX_edges / F_EtaX_edges[-1]                   
+
+    # marginal 1D CDF F_EtaY
+    pEtaY = Hs.sum(axis=0)                                  # Marginal PDF (nBinY,)
+    F_EtaY_edges = np.concatenate([[0.0], np.cumsum(pEtaY)])# CDF at bin edges (nBinY+1,)
+    F_EtaY_edges = F_EtaY_edges / F_EtaY_edges[-1]                   
+
+    # calculate (u,v) at the center of each original bin
+    Xc_1d = 0.5 * (F_EtaX_edges[:-1] + F_EtaX_edges[1:])
+    Yc_1d = 0.5 * (F_EtaY_edges[:-1] + F_EtaY_edges[1:])
+
+    X_tab = np.repeat(Xc_1d[:, None], nBinY, axis=1)        # (nBinX, nBinY)
+    Y_tab = np.repeat(Yc_1d[None, :], nBinX, axis=0)        # (nBinX, nBinY)
+
+    return (X_tab, Y_tab)
+
 def bilinear_xy_lookup(etaX, etaY, X_tab, Y_tab):
     """
     bilinear interpolation, from (eta_x, eta_y) to (u,v):

etaInterpolation.ipynb

@@ -2,14 +2,97 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
     "import numpy as np\n",
     "from glob import glob\n",
     "from matplotlib import pyplot as plt\n",
-    "from EtaInterpolationFunctions import build_xy_lut_from_hist, bilinear_xy_lookup, xy_lookup\n"
+    "from EtaInterpolationFunctions import *\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## A simple example\n",
+    "Mapping a 2D distribution to a uniform distribution using different interpolation methods.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Original x: mean=0.481, std=0.172\n",
+      "Original y: mean=0.522, std=0.174\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAABdIAAAH1CAYAAADyEAw6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOy9d7xdZZX//9l7n3b7zU3vCSmEDqFDgADSQrGAoFgQddRhBNs4iv6+Os6MYhvFEUFmaGNDAUURkN5Db6EkENJ7vb2cuvfvjwwprs+63pMQCfB5v155vWBl7Wc/ZT3rKefkfIIkSRIIIYQQQgghhBBCCCGEEIISvtkVEEIIIYQQQgghhBBCCCF2ZXSRLoQQQgghhBBCCCGEEEL0gy7ShRBCCCGEEEIIIYQQQoh+0EW6EEIIIYQQQgghhBBCCNEPukgXQgghhBBCCCGEEEIIIfpBF+lCCCGEEEIIIYQQQgghRD/oIl0IIYQQQgghhBBCCCGE6AddpAshhBBCCCGEEEIIIYQQ/aCLdCGEEEIIIYQQQgghhBCiH3SRLsR28Pjjj+P9738/Ro4ciUwmgxEjRuCss87CY489VlU5//qv/4ogCLarDg888ACCIMADDzywXc8PlJkzZ2LmzJkD8i2VSpg2bRq++93v7tQ6DYS2tjY0Nzfjj3/845tdFSGEEEIIIYQQQgjxFkcX6UJUyU9/+lMceeSRWLFiBb7//e/jnnvuwQ9/+EOsXLkSM2bMwGWXXTbgsj75yU9Wffn+OtOnT8djjz2G6dOnb9fzO4PLL78cbW1tuPDCC9/sqmDQoEH4whe+gC9/+csoFotvdnWEEEK8QVx33XUIgmDzn1QqhZEjR+IDH/gAXnvttTe7elURBAE++9nP/k2/HfnwfO7cufjXf/1XLFmyxPzdb37zG1x66aVVlbcrfWgOAFdffTVGjx6Nnp6eN7sqQggh/o+/XqtzuRxGjBiBY489FpdccgnWrVv3d6nH32Od7Y/Ozk58+9vfxkEHHYTGxkZks1lMmDABH//4x/Hss89u9tue/nr9S3nsz0DuJLSeC7F96CJdiCqYPXs2Pv/5z2PWrFl4+OGH8ZGPfARHH300PvzhD+Phhx/GrFmz8LnPfQ6zZ8/ut5ze3l4AwJgxY3DYYYdtV10aGxtx2GGHobGxcbuef6Mpl8v4wQ9+gI9//OOoq6t7s6sDAPjMZz6DJUuW4KabbnqzqyKEEOIN5tprr8Vjjz2Ge+65B5/97Gdxyy23YMaMGWhra3uzq7ZLMXfuXHzrW996wy7Sd6UPzQHgvPPOQ11dHb7//e+/2VURQgjxV7y+Vt9999342c9+hv333x/f+973sMcee+Cee+55s6u3U1m4cCEOOOAAfPe738Wxxx6L66+/HnfddRe+9a1vYe3atTjwwAPR0dGxzTPb01933HEHHnvssW3+vP/97/+b9dN6LsT2oYt0IargkksuQRAEuOKKK5BKpbb5u1QqhcsvvxxBEGzzqe7rnxQ/++yzOOusszBo0CBMmjRpm7/bmkKhgC996UsYMWIEamtrcfTRR+OZZ57BhAkT8LGPfWyzH/vU/GMf+xjq6+uxYMECzJo1C/X19Rg7diy+9KUvoVAobPOeb33rWzj00EPR0tKCxsZGTJ8+HVdffTWSJNmuvrnllluwcuVKfOQjH9lsy+fzOOCAAzB58uRtNglr1qzBiBEjMHPmTFQqlarec/zxx2PatGmmnkmSYPLkyTj11FM324YPH44TTjgBP//5z7erTUIIIXZd9t57bxx22GGYOXMmvv71r+OrX/0q1q1bp5/02onsih+ap1IpfPrTn8ZPfvKTzV9UEEIIsWvw+lp91FFH4cwzz8SPf/xjvPDCC6irq8P73vc+rF279s2u4k6hUqngve99LzZs2IDHHnsMP/jBD3DqqafimGOOwXnnnYfbbrsNt99+O9Lp9DbPbU9/HXjggTjssMO2+TN8+PB+66f1XIjtRxfpQgyQSqWC+++/HwcddBDGjBlDfcaOHYsDDzwQ9913n7kgft/73ofJkyfjxhtv7Pdi9/zzz8ell16K888/H3/6059w5pln4r3vfS/a29sHVM9SqYQzzjgDxx9/PP70pz/h4x//OH784x/je9/73jZ+S5Yswac//WnccMMN+MMf/oD3ve99uPDCC/Hv//7vA3rPX3Pbbbdh2LBh2HPPPTfbcrkcbrjhBqxbtw4f//jHAQBxHONDH/oQkiTB9ddfjyiKqnrP5z73Obz66qu49957t7H/5S9/wcKFC/FP//RP29hnzpyJ2bNnD7j/hBBCvDU56KCDAMAcMm+55RYcfvjhqK2tRUNDA0444QTzs2rr16/Hpz71KYwdOxbZbBZDhw7FkUceab79dc899+D4449HY2MjamtrceSRR5r16PUPyV9++WV88IMfRFNTE4YPH46Pf/zj5ptnr3PllVdi6tSpyGaz2HPPPfHb3/72b7b36aefxgc+8AFMmDABNTU1mDBhAj74wQ9i6dKlm32uu+66zd9KO/bYYzf/k+/rrrsOM2fOxG233YalS5du88/B+4N9aA5s+iB/woQJxn97tWDmzJmDIAhw9dVXm7/7y1/+giAIcMstt2y2fehDH0JnZ+eA+k0IIcSby7hx4/Cf//mf6OrqwpVXXrnN3w1kzd6eNWd71llg01p7xhlnoKWlBblcDgcccABuuOGGv/ncH//4R7z44ou4+OKLsffee1OfU045BbW1tX+zrP76a3vRei7E9pP62y5CCADYsGEDent7MXHixH79Jk6ciCeffBIbN27EsGHDNtvPO+88fOtb3+r32blz5+L666/HV77yFVxyySUAgBNOOAHDhw/HBz/4wQHVs1gs4lvf+tbmg/Pxxx+Pp59+Gr/5zW/wjW98Y7Pftddeu/m/4zjGzJkzkSQJfvKTn+D//b//V/VC6f1e+5QpU3DVVVfhnHPOwU9+8hO0trbigQcewB133IGRI0dW9Q4AOO2007Dbbrvhsssuw7ve9a7N9ssuuwyTJk3CKaecso3/9OnTEccxHn/8cZx88slVv08IIcRbg8WLFwMApk6dutn2m9/8Bh/60Idw4okn4vrrr0ehUMD3v/99zJw5E/feey9mzJgBAPjIRz6CZ599Ft/+9rcxdepUtLe349lnn8XGjRs3l/WrX/0KH/3oR/Hud78b//u//4t0Oo0rr7wSJ510Eu68804cf/zx29TnzDPPxDnnnINPfOITmw/TAHDNNdds43fLLbfg/vvvx7/927+hrq4Ol19+OT74wQ8ilUrhrLPOctu7ZMkS7L777vjABz6AlpYWrF69GldccQUOPvhgzJ07F0OGDMGpp56K73znO/ja176Gn/3sZ5vX6UmTJuGQQw7Bpz71KSxcuBA333zzgPqYfWi+M9hvv/1wwAEH4Nprr8UnPvGJbf7uuuuuw7BhwzBr1qzNthEjRmDatGm47bbbNn9wL4QQYtdl1qxZiKIIDz300GbbQNfsatnedfb+++/HySefjEMPPRQ///nP0dTUhN/+9rc455xz0Nvbu82/Fv9r7rrrLgDAe97znu2q81/D+ut1KpUKyuXy5v8PguBvfllN67kQ248u0oV4g3n9J0f++iL6zDPP/JvPPvjggwCAs88+exv7WWedZT4t9giCAKeffvo2tn333Rf33XffNrb77rsP3/nOd/DUU0+hs7Nzm79bt27d3/znYH/NqlWrcPDBB9O/O/vss/HAAw/gy1/+MiqVCr72ta/hhBNOqKr81wnDEJ/97Gfx5S9/GcuWLcO4ceOwcOFC3HHHHfjhD39o+v31DzNWrly5Xe8TQgixa/L6wTGfz2P27Nn4j//4Dxx99NE444wzAGz6kPjLX/4y9tlnH/zlL39BGG76h5izZs3CpEmT8JWvfGWzpsns2bPxyU9+Ev/wD/+wufx3v/vdm/+7t7cXn/vc53Daaadtc+k8a9YsTJ8+HV/72tfwxBNPbFO/T3ziE/jyl78MAHjXu96FBQsW4JprrsHVV1+9zVq1YcMGPPXUU5vX3VmzZmHvvffGxRdf3O8B/6yzztrm7yuVCk477TQMHz4cv/nNb3DRRRdh6NChmDJlCgBgzz333EaXZejQoWhubkY2mx2wXsvfU+T8/PPPx0UXXYT58+dv/nCkra0Nf/rTn/DZz37W/MTe9OnT3/a/tyuEEG8X6urqMGTIEKxatQpAdWt2tWzvOnvBBRdgr732wn333bd5zTnppJOwYcMGfO1rX8NHP/rRzfX8a5YtWwYAf/NLeAPlr/tra0aMGLHN/48ePRorVqzotzyt50JsP/ppFyEGyJAhQ1BbW7v5G28eS5YsQW1tLVpaWraxD+Tb169/8+2vL7FTqRQGDx48oHrW1tYil8ttY8tms8jn85v//8knn8SJJ54IAPif//kfzJ49G0899RS+/vWvAwD6+voG9K6t6evrM+/dmo9//OMolUpIpVK46KKLqi7/r8uqqanZ/BM5P/vZz1BTU0M/tX69TtvTJiGEELsuhx12GNLpNBoaGnDyySdj0KBB+NOf/rT5QPbqq69i1apV+MhHPrLNQbe+vh5nnnkmHn/88c2/wXnIIYfguuuuw3/8x3/g8ccfR6lU2uZdjz76KFpbW3HeeeehXC5v/hPHMU4++WQ89dRT6Onp2eaZ1y/0X2ffffdFPp/HunXrtrEff/zx26z7URThnHPOwYIFC/o9CHd3d+MrX/kKJk+ejFQqhVQqhfr6evT09GDevHlV9OTAWbVq1Tb/2m5n8qEPfQjZbBbXXXfdZtvr31A8//zzjf+wYcOwbt26bb6VJ4QQYtdla82ratbsatmedXbBggV45ZVX8KEPfQgAtln7Z82ahdWrV+PVV1/drvpsL56W2T333IOnnnpq85/bb7/9b5al9VyI7UcX6UIMkCiKcOyxx+Lpp592F9wVK1bgmWeewXHHHWf+OdVAfirl9cvLine truncated
+      "text/plain": [
+       "<Figure size 1500x500 with 3 Axes>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "rng = np.random.default_rng(42)\n",
+    "n = 10_000_000 # number of samples\n",
+    "\n",
+    "# non-uniform 2D distribution: spiral with noise\n",
+    "theta = rng.uniform(0, 4*np.pi, size=n)\n",
+    "r = 0.3 + 0.4 * (theta / (4*np.pi))  # 螺旋半径从 0.3 到 0.7\n",
+    "x_raw = r * np.cos(theta) + 0.1 * rng.normal(size=n)\n",
+    "y_raw = r * np.sin(theta) + 0.1 * rng.normal(size=n)\n",
+    "\n",
+    "# normalize to [0, 1]\n",
+    "x = (x_raw - x_raw.min()) / (x_raw.max() - x_raw.min())\n",
+    "y = (y_raw - y_raw.min()) / (y_raw.max() - y_raw.min())\n",
+    "\n",
+    "# build 2D histogram\n",
+    "nBins = 80\n",
+    "H, _, _ = np.histogram2d(x, y, bins=nBins, range=[[0,1],[0,1]])\n",
+    "H = H.astype(int)\n",
+    "\n",
+    "# build LUT\n",
+    "X_tab, Y_tab = build_xy_lut_Rosenblatt(H)\n",
+    "u, v = bilinear_xy_lookup(x, y, X_tab, Y_tab)\n",
+    "H_uv = np.histogram2d(u, v, bins=nBins, range=[[0,1],[0,1]])[0]\n",
+    "X_tab2, Y_tab2 = build_xy_lut_DoubleCDF(H)\n",
+    "u2, v2 = bilinear_xy_lookup(x, y, X_tab2, Y_tab2)\n",
+    "H_uv2 = np.histogram2d(u2, v2, bins=nBins, range=[[0,1],[0,1]])[0]\n",
+    "\n",
+    "# map to (u, v)\n",
+    "u, v = bilinear_xy_lookup(x, y, X_tab, Y_tab)\n",
+    "\n",
+    "# Visualization\n",
+    "fig, axs = plt.subplots(1, 3, figsize=(15, 5))\n",
+    "\n",
+    "axs[0].imshow(H.T, origin='lower', extent=[0,1,0,1])\n",
+    "axs[0].set_title(\"Original (x, y)\")\n",
+    "axs[0].set_aspect('equal')\n",
+    "\n",
+    "axs[1].imshow(H_uv.T, origin='lower', extent=[0,1,0,1])\n",
+    "axs[1].set_title(\"Rosenblatt (u, v)\")\n",
+    "\n",
+    "axs[2].imshow(H_uv2.T, origin='lower', extent=[0,1,0,1])\n",
+    "axs[2].set_title(\"Double CDF (u, v)\")\n",
+    "axs[2].set_aspect('equal')\n",
+    "axs[2].set_xlim(0,1); axs[2].set_ylim(0,1)\n",
+    "\n",
+    "\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
    ]
   },
   {
@@ -21,7 +104,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [
     {
@@ -102,7 +185,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [
     {
@@ -118,7 +201,7 @@
    ],
    "source": [
     "### build xy lookup table and plot\n",
-    "X_tab, Y_tab = build_xy_lut_from_hist(eta_2dhist)\n",
+    "X_tab, Y_tab = build_xy_lut_Rosenblatt(eta_2dhist)\n",
     "plt.figure(figsize=(8,4))\n",
     "plt.subplot(1,2,1)\n",
     "plt.imshow(X_tab, origin='lower', extent=(0,1,0,1))\n",
@@ -144,7 +227,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
@@ -224,7 +307,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [
     {
@@ -305,16 +388,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 17,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "<matplotlib.colorbar.Colorbar at 0x7faff241da00>"
+       "<matplotlib.colorbar.Colorbar at 0x7f1462d49a60>"
       ]
      },
-     "execution_count": 6,
+     "execution_count": 17,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -351,7 +434,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
@@ -383,7 +466,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [
     {
@@ -397,10 +480,10 @@
     {
      "data": {
       "text/plain": [
-       "<matplotlib.legend.Legend at 0x7fafe3a39d90>"
+       "<matplotlib.legend.Legend at 0x7f14dfb62b80>"
       ]
      },
-     "execution_count": 8,
+     "execution_count": 19,
      "metadata": {},
      "output_type": "execute_result"
     },