mirror of
https://github.com/slsdetectorgroup/aare.git
synced 2026-10-31 11:21:28 +01:00
ClusterFinderFrozen: CPU finder for CUDA correctness study
- Diagnostic twin of ClusterFinder: identical decisions, but the pedestal is frozen per frame (snapshot at frame start, updates deferred to frame end), matching the CUDA kernel's update model. - Isolates pedestal-update timing as the sole remaining CPU/GPU mismatch. - Adds class + bindings + factory, validation notebook, and helper utilities.
This commit is contained in:
kferjaoui committed
1 parent
a42d71cf42
commit
3504d96336
7 files changed
+1047
-2
No files matched your search
@@ -0,0 +1,360 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "md-intro",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Frozen-CPU vs CUDA — isolating pedestal-update timing\n",
|
||||
"\n",
|
||||
"`ClusterFinderFrozen` is a diagnostic twin of the serial `ClusterFinder`: **identical**\n",
|
||||
"decision logic (negative skip, Test1, Test3, `value == max` store gate, edge handling,\n",
|
||||
"rounding), differing in **exactly one** respect — *when* the pedestal is updated.\n",
|
||||
"\n",
|
||||
"| finder | pedestal update |\n",
|
||||
"|---|---|\n",
|
||||
"| `ClusterFinder` | per-pixel, **during** the raster scan (scan-order dependent) |\n",
|
||||
"| `ClusterFinderFrozen` | frozen snapshot for all decisions; **deferred** to frame end |\n",
|
||||
"| `ClusterFinderCUDA` | frozen per frame; batch update at frame end |\n",
|
||||
"\n",
|
||||
"`Frozen` and `CUDA` share the same update *model*.\n",
|
||||
"\n",
|
||||
"**Expected result** (CUDA at `n_streams=1`):\n",
|
||||
"- `cpu vs cuda` → the current gap (what we're explaining)\n",
|
||||
"- `frozen vs cuda` → **≈ 0** (timing was the whole story)\n",
|
||||
"- `frozen vs cpu` → ≈ the gap (exactly the flipped frames)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "imports",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import sys; sys.path.append('/home/ferjao_k/aare/build')\n",
|
||||
"\n",
|
||||
"from pathlib import Path\n",
|
||||
"import numpy as np\n",
|
||||
"import time\n",
|
||||
"\n",
|
||||
"from aare import File, ClusterFinder, ClusterFinderFrozen, ClusterFinderCUDA\n",
|
||||
"from helper import (centers, only_sets, train_pedestal,\n",
|
||||
" compare_finders, print_comparison, plot_spectra,\n",
|
||||
" scan_mismatches, plot_masked_mismatch, walkthrough)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "config",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"image (400, 400) pedestal 1000 data 50000 streams 1\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"base = Path('/mnt/sls_det_storage/matterhorn_data/aare_test_data/')\n",
|
||||
"f = File(base / 'Moench03new/cu_half_speed_master_4.json')\n",
|
||||
"\n",
|
||||
"n_frames_pd = 1000\n",
|
||||
"N = 20000\n",
|
||||
"cluster_size = (3, 3)\n",
|
||||
"rows, cols = f.rows, f.cols\n",
|
||||
"image_size = (rows, cols)\n",
|
||||
"capacity = 50_000\n",
|
||||
"\n",
|
||||
"N_STREAMS = 1 # single device pedestal -> apples-to-apples with the CPU finders\n",
|
||||
"N_SIGMA = 5\n",
|
||||
"sx, sy = cluster_size\n",
|
||||
"rx, ry = sx // 2, sy // 2\n",
|
||||
"print(f'image {image_size} pedestal {n_frames_pd} data {N} streams {N_STREAMS}')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "build",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"cf_cpu = ClusterFinder(image_size, cluster_size, n_sigma=N_SIGMA, capacity=capacity)\n",
|
||||
"cf_frozen = ClusterFinderFrozen(image_size, cluster_size, n_sigma=N_SIGMA, capacity=capacity)\n",
|
||||
"cf_cuda = ClusterFinderCUDA(image_size, cluster_size, n_sigma=N_SIGMA,\n",
|
||||
" max_clusters_per_frame=1500, n_streams=N_STREAMS)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"id": "train",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"pedestal train: 0.80s\n",
|
||||
"data: (50000, 400, 400) uint16\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Train ALL three on the SAME pedestal frames, then load the data block.\n",
|
||||
"f.seek(0)\n",
|
||||
"t0 = time.perf_counter()\n",
|
||||
"for _ in range(n_frames_pd):\n",
|
||||
" img = f.read_frame().copy()\n",
|
||||
" cf_cpu.push_pedestal_frame(img)\n",
|
||||
" cf_frozen.push_pedestal_frame(img)\n",
|
||||
" cf_cuda.push_pedestal_frame(img)\n",
|
||||
"print(f'pedestal train: {time.perf_counter()-t0:.2f}s')\n",
|
||||
"\n",
|
||||
"f.seek(n_frames_pd)\n",
|
||||
"data = f.read_n(N)\n",
|
||||
"print('data:', data.shape, data.dtype)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "md-compare",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## The three-way comparison\n",
|
||||
"\n",
|
||||
"All three finders now carry the **same** trained pedestal. `compare_finders` runs each\n",
|
||||
"over the same sampled frames and reports exact (tol=0) pairwise centre mismatches."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"id": "compare",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Scanned 20000 frames\n",
|
||||
"\n",
|
||||
"Total clusters per finder:\n",
|
||||
" cpu 27,666,006\n",
|
||||
" frozen 27,665,954\n",
|
||||
" cuda 27,665,796\n",
|
||||
"\n",
|
||||
"Pairwise exact mismatches (tol=0):\n",
|
||||
" pair A-only B-only total\n",
|
||||
" cpu vs frozen 94 42 136 (0.0005%)\n",
|
||||
" cpu vs cuda 252 42 294 (0.0011%)\n",
|
||||
" frozen vs cuda 158 0 158 (0.0006%)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"SCAN = 20000 # frames sampled across `data` (set to len(data) for the full block)\n",
|
||||
"\n",
|
||||
"totals, pairs, nscan, hists = compare_finders(\n",
|
||||
" # {'frozen': cf_frozen, 'cuda': cf_cuda}, data, scan_count=SCAN)\n",
|
||||
" {'cpu': cf_cpu, 'frozen': cf_frozen, 'cuda': cf_cuda}, data, scan_count=SCAN)\n",
|
||||
"print_comparison(totals, pairs, nscan)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "md-read",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**Reading the table.** If `frozen vs cuda` collapses to ~0 while `cpu vs cuda` and\n",
|
||||
"`frozen vs cpu` are comparable and non-trivial, pedestal-update *timing* is the proven\n",
|
||||
"cause. A small non-zero `frozen vs cuda` residual would be a genuine surprise (FP corner,\n",
|
||||
"edge pixel, or the multi-stream path) worth chasing — not expected background."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "424e2976",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Cluster-energy spectra\n",
|
||||
"\n",
|
||||
"The three finders' cluster-energy distributions, accumulated in the same compare pass\n",
|
||||
"(no extra scan). They should overlap almost perfectly; the ratio panel (each vs `cpu`)\n",
|
||||
"makes any per-bin divergence — e.g. the `frozen vs cuda` residual — visible."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"id": "f057e296",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAAAxYAAAJOCAYAAAAqFJGJAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAyqVJREFUeJzs3QWcVNXfBvDfdrHF9i7L0p3SJSCpgt2C6B8VCwMRxAQUASUMLOwAeVXEApFQURqR7lhYdtmE7Y77fp6Dd5zZ7pnZeb5+xsvcOXPvnXPPzJ7fPXHtNE3ThIiIiIiIqAbsa/JmIiIiIiIiBhZERERERFQr2GJBREREREQ1xsCCiIiIiIhqjIEFERERERHVGAMLIiIiIiKqMQYWRERERERUYwwsiIiIiIioxhhYEBERERFRjTGwIJu2f/9+ueeee6R58+bi6uoqjRo1kssuu0xeffVVuXjxoiHdkCFD1KOuvPPOO/Lpp5/W2fbJethaWbj77rulWbNmdbb95cuXy+uvvy6WrK5/X2ri1KlT4uLiItu2bTOs++qrr+Tyyy+XoKAg9VpoaKiMHTtWtm7dWu397N27V66++mpp2rSpuLm5SePGjaVfv37y5ZdflkiLfT/++ONiLbKysmTmzJnyxx9/mPtQLFZ0dLQ6p4MHDxYfHx+xs7Ozqd/BhoSBBdmsDz74QHr06CG7du2Sp556StauXSurVq2Sm2++Wd577z2ZOHFivR2LrVUmqWwsC7YXWOCc42GJpk6dKiNGjFCVfN2FCxdkwIAB6pjXrVsnixYtkvj4eFXh37RpU7X2k5KSIuHh4fLKK6/ImjVr5PPPP1cB5/jx4+Xll182SfvSSy+pfR87dkysJbCYNWsWA4tynDx5UpYtWybOzs5y1VVX1d/JoVrnWPubJLJ8uPr24IMPqj+Y33//vbrqpsO6J598UgUa1kzTNMnJyVFX/6xBdna21RyrpcjPz1dX9hwd+VNuzeW1Q4cOYomOHDmifh+L/xY+8sgjJdJeeeWVEhAQIB999JG66lwbrTZjxoyRyMhIWbp0qTz33HOG9dh+27ZtZeHCheq1hgaBiLu7u9gSBKWJiYnq33///bdqFSPrxBYLskm4KoYKGf4oGQcVOlw1ueaaa8p8P5q08f7iTdtnzpwp0YR7+vRpue2221R3AewL3QeGDRummv4BV+UOHTqkrvThvXgYdw1JS0tTVw3RXQvHFRYWppqMMzMzTfaN9+EPPlpb2rdvr/b12WeflZsP//d//6euRHp4eKhuYKNGjZI9e/aU6KqC13BFCVeS8G9cWUTwlZuba5I2Ly9PXV1s166d2j8qGuhqpv/B0OHzodLw3XffSffu3VU3NFzRA+TFyJEj1R9WvP/hhx+W1atXm+Q3rliiMn3u3LkSn+l///uf+Pn5qaCqLBWdE+NjRCtWly5d1DG2aNFC3nzzzRLbq+w5Kioqkrfeeku6deumKqVo8u/bt6/8+OOPFZYFvcx98cUXKu+xDxw7zgvy96GHHlIVVJyfwMBAueKKK+Svv/4SS2gxQBnDceGBz47KZ1lK+w7psB5dSnT43Pfff78qj3p5w5X0DRs2qNdRUUXZOXv2rCE/8ajN8lqRypS14pVqfOeMj9f4Yfz5K1vuquvdd9+V4OBgdbGlIp6enipfajvI9ff3L3WbaMlA2UpPT6/WdpHfnTp1Ui3WgwYNUr83+H7PmzdPfU+r4rffflPbw+8OvtfoznXjjTeqAAHlGeUKUGb084hzDDifeP7PP//ITTfdJL6+vtKyZUvDxSG0zOi/F3gNaVCmjK1fv16uvfZaadKkiToHrVq1kkmTJklSUpJJOn1f6AKMlnlvb2/V5WzKlClSUFCgWoBGjx6tziXKPLoEG0O+4PuCoE7//cJv4xtvvCE1ZW/P6mhDwctcZHMKCwvVHwJ0g0KFpK6hMo594kcaf3DwY4++yGj6B1Rc8ccCP/J6dwg92MEfJlydQ//TZ555Rv2Io+L5wgsvyIEDB1QFyriihKuLqEzidVQIUMEsL7jCVUBUpLBEJeu1115Tf2R37txpchUVV8YRaKF7GCq1f/75p6rc45ixL/2PDv64Yf/Tpk2T/v37qwrdiy++qP7o4iqU8RVe/CHFFVHsGxUjBDexsbHq8+LfqNTg+HHlqvgVUvzRnDNnjrz//vsm3SQwLmbFihUqPf7AVvec6FD5Q0UNf5CRn2iqf+yxx1ReoUJX1XOEygT6jCMfZ8+erSqDyAdUPioqC7oZM2aoijoCSPwxRh7pFWHkNY4zIyNDbQv5vnHjRrP130ceoJzccMMNqtzgcx08eFCVi9qAyiXyD2WhTZs26vzhObrqAPIQgQfGCSA/jNVGea2MypY1Y88//7w88MADJuvefvttVXb072VVyh32j0pqRVCejCt4CMpwJbmsSh+2i3yMiYmRuXPnqn3gQkBNYHt4JCcnyzfffCO//vqrLFmypEQ6nKPp06ergBvjO6ojLi5O7rzzTlU2cd5RRvD9QhB41113VWob+O5ibAh+Nz/++GNV2UZ+oJUHvxMhISHq36iw43t/7733qvfpwYYO3xEEoDjvemCI3zkE2I8++qjMnz9f/b7hdwNldd++fSpIBZRv/CZg2/iO4ZjQPW3gwIGqLDg5OZns65ZbbpFx48ap7SMoQdnEbzzKDC5Q4LcNQRvyF0EKjg2QDr+F+A6gXOA9R48eNSnLKAMoF5XBltYGSiOyMXFxcfgLq912222Vfs/gwYPVQ/f777+rbWBpLDIyUq3/5JNP1POkpCT1/PXXXy93+x07djTZvm7u3Lmavb29tmvXLpP13377rdrumjVrDOvw3NvbW7t48WKFnycqKkpzdHTUJk+ebLI+PT1dCw4O1m655RbDugkTJqhtf/311yZpr7rqKq1t27aG51999ZVKt3LlSpN0OHasf+eddwzrIiIiNAcHB+3YsWMmaZ966inNzs5OO3TokMn6UaNGlchvHFdgYKCWm5trWDd//nyVXzgPZansOcEx4lj27t1rsn7EiBGal5eXlpmZWaVz9Oeff6rnzz77bLXKgl7mLr/8cq0iBQUFWn5+vjZs2DDt+uuv18zh9OnT6hzfeeed5abDeURel/UdMob1L774ouF5o0aNtMcff7zc7V999dUm26/N8lqRypa14r8vxeG7h7L4zDPPVOu3AdvGuooeOBe6+Ph4tW7evHllHhe+//p7Q0JCtM2bN2s1NWnSJMM2nZ2dTc6Dsby8PJUn06dPr9Z+9DzZsWOHyfoOHTqo35vK0vO7+O+EscTExBJlV4d1eO2FF14wWb9t2za1fuHChSbrz507p7m5uWnTpk0rdV9FRUXqu3/27Fn1/h9++KHEvopvs1u3bmr9d999Z1iHbQQEBGg33HCDYd2YMWNU2vLov1OVeZT1O61/B0v7DSDLxxYLojqEZmY0a6MlAFdxhg4dKl27dq10s+/PP/+smuvRFI6mah26LOldg9C3WYfuL2gurwiuAmJ7uCpnvF1c5cdV0N9//90kPfZV/KogrpCi5cf4WHG1DumMt4ljx1V0HCvGtRi/H1eZjaELED5v8T7nt99+uzpmY2g5QFcvXNXEVUdc5UQrB64eljfLUFXOSceOHdVrxu644w51lQ9XsHFFsLLn6JdfflHra3pFF10sSoMWDHTtO3z4sEkXNXTzqcwV4uoofoXbGPII+VvTz1ue3r17qyu66IIyfPhw1QpZ/OpsWWqjvNb191//TqBlBleY0TJTnd8GtOxVpssQuh3pzp8/r5bltXquXLlSXV2PiopS5Q/7Q7e+mrSQofUFV94TEhLkp59+Uq2P2IfeQqjDedZbB6oL5xllyBjOs3E3tYog/9HyiJYxXO1HywW6VNX0e43zi/OI8258fnHMKEPG3XCRV2ipQgsTzpvx9xmtbMW79aJbnzF0nUULiPHfErQmoLXCuHUReYV94HOitQ+tJF5eXibb0idEqQy0DFHDw8CCbA7+eKI/LQYF1jX8YUBXFDRfoxkZTe6obKAijEoC+rKWBzOtoA99WZWl4n1o0exeGdgu9OrVq9TXi1d8kF/Fuxahi47xOAZsE03i+CNb3WNFFxZ0MylOb/I3hr7u+COOLiLIT/whRhcAVKJq65zgj3hx+jq9u01lzxG6Kzk4OJS6zaooLd/Q7QGfA90o0PUIZRz7QpcaVCzKg3yo7HiB4tB9xLjPvzG9exb6fdcVjBFCV7gPP/xQfVaM4bj++uvVea0on2ujvNb19x9dm6677jpVzouPS6nKbwMqiJXtCmU8OB3K61KIwFuvcOI48Z1EwI9KanWhuxgeoM8OhO5JEyZMKNF9CMemH2d1ICAtDr9rVdkmAkd0IcL5RRCNIAiBBbovIS8qq3j5wvnFOSvttw/04AVBBMakIaDAd6Bz586qmx7WY/xWaZ8FZdAYvgOl/cZjPcbx6HAesG10yUMgid8YdIlCN62ePXuqNPo4qspgV6iGiYEF2Rz8GGLwJK4go39ydSo++g9w8cHLxSsjEBERYagUHD9+XL7++mtVGUP/W/w4lwcVRPTzRt/dsl43ZjzeoqLtwrfffquOrzZgm/hDXdZsWsUrUaUdK96vBz3F+0KXBn+8MQgRrQfoh40rypUZaFrZc1LafvV1eqWksucIlSJctcb7q1NJLS/f8IceV4nRYmOsMlepcaW1+BXMyirviqNeCcR3rCpjmcr6bumBXPG8xVSyeOCqOa6WP/300+oKbkWzutVGea2M6n7/kW/ol49KNloGigcQVfltwO9dZaaBReVdHzSvv9/4fj4VVRJxDyB8vtqEoAX5hAHLxQMLjMMo/htoDgj88MD3G2NzMEEDxmYhKMC4icooXr7wubAOY4BKm2BEX4cxSwjkcN5w/nQIOmsbzjEGeuOBoBwBFVqY0EqGiTQQnKCcoWWuMnBxry7vYUPmwcCCbBKuvGCu9Pvuu09++OLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 800x600 with 2 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fig = plot_spectra(hists, totals,\n",
|
||||
" title=f'Cluster energy spectra — cluster_size={cluster_size}, n_streams={N_STREAMS}')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "md-resid",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Eyeball any residual `frozen vs cuda` mismatches\n",
|
||||
"\n",
|
||||
"Fresh finders, retrained identically, then the masked side-by-side view of the worst\n",
|
||||
"residual frames. If the table showed 0, this should find nothing to plot."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "resid",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"frozen-only: 158 cuda-only: 0 frames shown: 8\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png":
|
||||
"iVBORw0KGgoAAAANSUhEUgAABEEAABCFCAYAAADTiXqrAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3Qd4FNXXx/FfCiWUhN57EelNEaQpRUWx0ZSmggULUhRUVBQpitIUBbtiwb+CYAEsoCiICqigIk3pTTok1NR9n3N5NyYhgUACKfP9PM+6u7Mzd+/MDs7JmVsCfD6fTwAAAAAAANlcYEZXAAAAAAAA4HwgCQIAAAAAADyBJAgAAAAAAPAEkiAAAAAAAMATSIIAAAAAAABPIAkCAAAAAAA8gSQIAAAAAADwBJIgAAAAAADAE0iCAAAAAAAATyAJgizP5/PprrvuUqFChRQQEKDff/89o6sEZKhNmzbxbwEA/h9xApA2FSpU0PPPP89hRLZBEgRZ3ldffaUpU6Zo9uzZ+vfff1WrVi1lBX369FHlypUVEhKiokWL6vrrr9eaNWtOuuhYYifh45FHHon/fN++fbrqqqtUqlQp5cqVS2XLllXfvn0VERGR7HeuW7dO+fPnV4ECBRItX7RokZo2barChQu7+lx44YWaMGGCzqfvvvtOV199tatDnjx5VKNGDT344IPavn27+/z7779PdBzsmLVr105//PHHaS/Stsw+OxMrV65Ux44d43+D5ModNmzYSb9PiRIlTlrHjmfevHlVsGBBtWnTRkuWLIn/fP/+/br//vtVrVo1t9/lypVTv379FB4erszC/n0lPWcAIKsgTvBGnHCqa5Utt8/9El637fpctWpV3Xbbbfrtt9+S3X7btm3KmTOn2++zcfz4cVd+7dq1FRwcrBtuuOGkdZLGOf5H0thwxowZbt/t97TnTz755KSyJk+erIoVKyp37txq2LChfvjhB2UW3KhBZkASBFne+vXrVbJkSV166aXuD1C7uCQVFRWlzMYuSm+//bZWr16tr7/+2t2puuKKKxQbG5toveHDh7vkjv/x+OOPx38WGBjokieff/65/v77b3eB/+abb3T33Xef9H3R0dHq2rWrmjdvftJnFgBYULRw4UJXH/sOe7z22ms6H1599VWXHLDfzy7uq1at0iuvvOISAePGjUu07tq1a91xmDNnjg4cOOCSQOciYXD06FFVqlRJo0ePPimxkVDNmjUT/T4rVqxI9PkFF1ygl156yS23INKSKvY779mzx32+Y8cO9xg7dqxbx35DC9hvv/12ZTd2bsfFxWV0NQB4DHGCt+KE1LIYzK7bdtNj0qRJOnz4sC655BK9++67J61rx61Lly4uNvjxxx/P6vpnySO7yWH7cSr+OMf/sASN388//6ybbrpJPXv2dDeB7NnqlfDmykcffaQBAwboscce0/Lly93vaTeNtmzZouzGzlngrPiALOzWW2/12Wnsf5QvX94tb9mype++++7zDRw40Fe4cGFfixYt3PJx48b5atWq5cuTJ4+vTJkyvnvuucd36NCh+PLefvttX1hYmG/WrFm+Cy64wBcSEuLr2LGj7/Dhw74pU6a48gsUKODr27evLyYmJn67yMhI3+DBg32lSpVyZTdq1Mj33XffndG+/PHHH24f1q1bF7/Mvm/ChAlnVM4LL7zg9i2phx56yNejR4/4fTydG2+80a1/rm3dutWXM2dO34ABA5L9/MCBA+7ZjqcdH/97s2jRIrfsq6++OuXxsmX+c+NspFTuk08+6atbt+4ZlRUeHu7q/M0336S4zrRp09wxiY6OTnGd2NhY3+jRo32VK1d265YtW9Y3cuRI99nGjRvddyxfvty9T+43/+STT9w6fr///rvvsssu8+XLl8+XP39+X4MGDXy//PJL/HFP+LD9Ts15n/DfU/Xq1X1BQUG+DRs2uHUuvvhit419fumll/o2bdp0RscRAFKDOME7ccKp6m3L7XM/u5bZdTCpW265xV0D9+/fH78sLi7OV6lSJRdrPPzww75evXql+Zy8/vrrT1qeXJyTVJcuXXxXXXVVomVXXnml7+abb45/b9fiu+++O9E6F154oe+RRx45Zb0+++wzX8OGDX25cuVysbP9vsnFQUljDGN1tmX+GMCOX7du3XxFihTx5c6d21elShXfW2+95T5LGlNYzO5n61hdrQ7VqlXzTZo0Kf4z//d+9NFHbhtbx9a3+KF9+/YuPre4okaNGr45c+accl8BWoIgS3vhhRdcS4kyZcq4bPkvv/wS/9k777zjWoVYxt7uIPhbTkycOFF//fWX+3z+/Pl66KGHEpVpWX5b58MPP3R35K15YocOHfTFF1+4x3vvvefufHz88cfx2/Tq1ct9j23z559/qnPnzq6Fwj///JOq/Thy5Ii7I2FNF62pakLPPvusa/pZr149jRo16pStWqxFwcyZM9WyZctEy20/p0+f7u50pIbdOfjpp59OKudcsHrZPiX9HfxO1Q3D7qqczZ0Af5NTa5KZVvYbWzNj++1uvvlmbdiwIcV1bT/t3AkLC1PdunVTXM/ubIWGhibbqslvyJAh7twYOnSouyP2wQcfqHjx4me9H927d3f/juzfkDUHtm5XOXLkcC2srCuQ1cd/V2rQoEGpPu/t39MzzzyjN954w91ts7F7rBmwnVu2jd3VsjF97PcAgPRGnODtOOFMDRw4UIcOHdK8efMSdcOxa5m14LCWF9OmTXPrJGTXsITdbdKifv36roVz69at3XcnZNdMa02a0JVXXul+C2PHya7hSdex9/51kmOtay3Wveaaa9xv++233+qiiy46633wxyZffvmlazn08ssvq0iRIu6zpUuXumdrkWQxhZ2P5vXXX3etVyzWtW2efvppV47F6wk9/PDDrkWNrWP7ft999ykyMtK1UrIWtRYb5cuX76zrDo8gD4SsLrm7/JYhrlev3mm3tTvulu32s7sESVtj9OnTx2WWE7YYsay7LTe2bkBAgG/79u2Jym7durVvyJAhp/x+y3DnzZvXfadlvhN+rxk/frzv+++/d61EXn/9dZdRv/32208qx+4AWKsVK+faa6/1HTt2LP6zvXv3ulYCCxYsOO2dktKlS7u7LYGBgb7hw4f7zgdrjRMaGnra9ZLeIbH9uu6669wdm127dp1RS5AlS5a4Owzbtm1LVR1TKveLL77wffzxx74///zTN2/ePHfeFS9e3NUtIWsJYb+znSfWamLp0qUpfpdtW65cOd9jjz2W4joRERHuDoidE8k5m5YgdhyttVNykts+Nee9/9+TtTLx27dvn1tm5zUAnA/ECd6IE9KjJYgdF/vs2WefjV9mLRoStkKxFqBJr78WU8ycOTNNLUHWrFnje+2113y//fab76effnL7bddZ/+9icuTI4Zs6dWqi7ey9/SbGrslW/x9//DHROqNGjXItnFPSpEkTX/fu3VP8/Exbgtg5llKLmeS2N3YOfvDBB4mWjRgxwtUt4XbPP/98onVq167tGzZsWIp1B5KT8m1GIItLLoNtGXXLLFt22gYFi4mJcYNVWUsM6+9qbLAtG7DUz+6u2zgOCbPKtmz37t3u9bJly9x4Hjb2Q0KWlbYWHKe7+962bVuXCbcxIaxfp91Zt4Gs/Hck/OrUqeMG1uzUqVN86xA/G5zsySefdP1IH330UT3wwANuUCxz5513qlu3bmrRosVpj5kNnGV9YhcvXuxaAlSpUsX1D05pLIzNmzcrtXr06OH67yZlx+5MWgFYawVjv5n1k7U7RMWKFdOZaNSo0UkDjZ0N62PrZ4OdNWnSxJ07dtfCfgO/yy+/3M1atHfvXnenw99/N2m97Zy0uzA20Jn9nimxux92ftldovRi9b3jjjtcSye722WtOhL+O0gqtee9DSRn566ftQSxweHs7o2d+/ZddjzsrhcAnE/ECdkzTkiLE/mREy07zMGDB11LBRvTK2E933rrLXfN9EuPmMIGSLeHn8UUW7dudfFhwt8m6bFI7vikZp2ELEax8yC93HPPPW5weYsVrBWKtQC1lqUpsXHSbF9tPLSE9bA43VrPnurfrbUKse+bO3euiynsexPGHUBySIIg2/InNfzsQmyjittgYCNGjHB/jNlFzf6Hm7A7hXUBSMguGskt8w/waM9BQUGu+aE9J3S65nj2P3Z72B/zjRs3dkkOG+U7pYDC1vGP3p7wD00bKMweNmq5LbdBsKwJof1haU1cbUA0u4j6L4RWZ+tqYV0zevfuHV+Odenw/0G/a9cuN7NJSnWxrkFn0g3FulMkx/6Itu4flghKzR/CFoBZWTY7TNIy7X1yg6RaEJP0Inquzjk7dkm7QdlyCxTtYb+h/d5vvvmm69LiZ01rrSuJnTN2DiQ955LrBpRa1g3MH9j5Jf3t7Le2INiaxFrzVUvCWDeXG2+8MdkyU3veW12TBl7W9cuCFutuZgO42eB61vTYf34DwPlAnJC94gTb3hI0NghpwuuSvbflqYkD7CZDwv20rqZ2s8wGTPXzHx+7oWY3Lc4luy6+//778e8t1tu5c2eideymnL87rHU5sX0/1TppjSsspjAJ44qkv7PdJLK422IK6/ZiN22s24r/HEvKH1PbjaKEx9okjTGS/ru1ZJTdWLHvskSIdcG1wXJt5j0gJSRB4Bm//vqryyjb/xj9/wO3fp3p0XfTLrB2gUluRPUzYRcUu5OeEuunaU4VBPgvSv5yrP9owhlnPvvsM9eSxPqGli5d+qzrUr58eaUHa9lid5Oee+65ZKfbswRGwv6+Fpik1P/XkkAJx4Xxs2UJ766cK3a8LIA63XmQ9NhaCxC7gNt0dxaI+lsCpcSSKBawWJ/dhHeiUmIJI0uyJGzxZHd9kgs07WEtkCyotWSFJUGsNUfSWYvSet7b9vawRJDd7bJAkyQIgIxEnJC14wSLAey6ZLFSwtYC1hrBlqcmDvLine truncated
|
||||
"text/plain": [
|
||||
"<Figure size 1100x4240 with 16 Axes>"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"cf_frozen2 = ClusterFinderFrozen(image_size, cluster_size, n_sigma=N_SIGMA, capacity=capacity)\n",
|
||||
"cf_cuda2 = ClusterFinderCUDA(image_size, cluster_size, n_sigma=N_SIGMA,\n",
|
||||
" max_clusters_per_frame=1500, n_streams=N_STREAMS)\n",
|
||||
"train_pedestal([cf_frozen2, cf_cuda2], f, n_frames_pd, seek=0)\n",
|
||||
"\n",
|
||||
"# scan_mismatches(a, b, ...) snapshots a.pedestal/noise (host) and b.device_*(0).\n",
|
||||
"show, res = scan_mismatches(cf_frozen2, cf_cuda2, data, rx, ry,\n",
|
||||
" scan_count=SCAN, n_show=8, tol=0)\n",
|
||||
"print('frozen-only:', res['cpu_only'], ' cuda-only:', res['cu_only'], ' frames shown:', len(show))\n",
|
||||
"if show:\n",
|
||||
" plot_masked_mismatch(show, data, rx, ry, rows, cols, zoom=30, show_vals=True)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d470cd40",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Walkthrough of the surviving mismatch\n",
|
||||
"\n",
|
||||
"Manual Test1/Test3 recompute of the strongest residual cluster in `show`, under each\n",
|
||||
"finder's decision-time snapshot pedestal. For 7×7 (a single residual) this shows exactly\n",
|
||||
"what the lone mismatched cluster looks like and which test tips it — naming the cause\n",
|
||||
"(FP knife-edge at threshold, a border/tie corner, or a genuine third source)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"id": "c27d894d-64b2-45b2-aaa8-bfddc754b3f9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"frame 35343 centre (x=19, y=349) accepted only by frozen\n",
|
||||
"raw window (ADU):\n",
|
||||
"[[5998 6480 6226]\n",
|
||||
" [5952 9428 6189]\n",
|
||||
" [6250 6202 6378]]\n",
|
||||
"\n",
|
||||
"--- frozen ---\n",
|
||||
" subtracted window:\n",
|
||||
" [[ -54.9 258. -127.8]\n",
|
||||
" [ -56.4 3138.6 46.2]\n",
|
||||
" [ 21.3 -91. -34.8]]\n",
|
||||
" centre value = 3138.63 (local max? True)\n",
|
||||
" max = 3138.63 total = 3099.23\n",
|
||||
" rms(centre) = 54.017\n",
|
||||
" Test1: max > 5*rms = 270.09 -> True\n",
|
||||
" Test3: total > 3*5*rms = 810.26 -> True\n",
|
||||
" ACCEPT = True\n",
|
||||
"\n",
|
||||
"--- cuda ---\n",
|
||||
" subtracted window:\n",
|
||||
" [[ -54.9 258. -127.8]\n",
|
||||
" [ -56.4 3138.6 46.2]\n",
|
||||
" [ 21.3 -91. -34.8]]\n",
|
||||
" centre value = 3138.63 (local max? True)\n",
|
||||
" max = 3138.63 total = 3099.23\n",
|
||||
" rms(centre) = 54.017\n",
|
||||
" Test1: max > 5*rms = 270.09 -> True\n",
|
||||
" Test3: total > 3*5*rms = 810.26 -> True\n",
|
||||
" ACCEPT = True\n",
|
||||
"\n",
|
||||
"RESULT: frozen = ACCEPT , cuda = ACCEPT\n",
|
||||
"pedestal mean gap @centre = 0.0000 ADU; rms gap = 0.0000\n",
|
||||
"NOTE: recompute agrees for this centre pixel — the split is FP rounding right at threshold or a window-neighbour effect; try another pick.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# Dissect the surviving frozen-vs-cuda residual (best seen for 7x7, where it's ~1).\n",
|
||||
"# Recomputes Test1/Test3 under each finder's decision-time snapshot pedestal, so the\n",
|
||||
"# lone mismatched cluster is fully explained (which test flips, and the pedestal gap).\n",
|
||||
"walkthrough(show, data, rx, ry, rows, cols, N_SIGMA, pick=0, labels=('frozen', 'cuda'))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e2f86c8d-5d44-40ad-9fa7-80997aa260e6",
|
||||
"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.11.15"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
@@ -0,0 +1,404 @@
|
||||
"""Helpers for the CPU-vs-CUDA ClusterFinder mismatch analysis.
|
||||
|
||||
Pure, parameterised versions of the functions that used to live inline in the
|
||||
notebook's forensic-view cell, plus the reusable scan + plot routines so the
|
||||
notebook cells stay thin. Import in the notebook with::
|
||||
|
||||
from helper import (centers, only_sets, footprint_mask, shift_dist,
|
||||
train_pedestal, scan_mismatches, plot_masked_mismatch)
|
||||
|
||||
None of these rely on notebook globals — cluster geometry (rx, ry) and the
|
||||
frame shape (rows, cols) are passed in explicitly.
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import boost_histogram as bh
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# Pure helpers
|
||||
# --------------------------------------------------------------------------- #
|
||||
def print_pinning_budget(rows, cols, dtype=np.uint16, headroom_gb=4.0):
|
||||
"""
|
||||
Print system RAM stats and estimate the maximum number of frames
|
||||
that can be safely registered with cudaHostRegister / register_input_buffer.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rows, cols : int Detector frame dimensions.
|
||||
dtype : np.dtype Frame element type (default uint16 = 2 bytes/pixel).
|
||||
headroom_gb: float RAM to keep free for OS + CUDA context (default 4 GB).
|
||||
"""
|
||||
# /proc/meminfo is always available on Linux — no extra dependency needed
|
||||
meminfo = {}
|
||||
with open('/proc/meminfo') as f:
|
||||
for line in f:
|
||||
key, val = line.split(':')
|
||||
meminfo[key.strip()] = int(val.split()[0]) * 1024 # kB → bytes
|
||||
|
||||
total_bytes = meminfo['MemTotal']
|
||||
available_bytes = meminfo['MemAvailable'] # free + reclaimable cache
|
||||
safe_bytes = max(0, available_bytes - int(headroom_gb * 1024**3))
|
||||
|
||||
frame_bytes = rows * cols * np.dtype(dtype).itemsize
|
||||
max_frames = safe_bytes // frame_bytes
|
||||
|
||||
GiB = 1024**3
|
||||
print("── System RAM ──────────────────────────────────────────")
|
||||
print(f" Total RAM : {total_bytes/GiB:.1f} GiB")
|
||||
print(f" Currently available : {available_bytes/GiB:.1f} GiB "
|
||||
f"(free + reclaimable cache)")
|
||||
print(f" Reserved headroom : {headroom_gb:.1f} GiB "
|
||||
f"(OS + CUDA context)")
|
||||
print(f" Safe pinning budget : {safe_bytes/GiB:.1f} GiB")
|
||||
print()
|
||||
print("── Frame layout ────────────────────────────────────────")
|
||||
print(f" Frame size : {rows} × {cols} × "
|
||||
f"{np.dtype(dtype).itemsize} B = {frame_bytes/1024:.1f} kB")
|
||||
print()
|
||||
print("── Pinning estimate ────────────────────────────────────")
|
||||
print(f" Max frames pinnable : {max_frames:,} "
|
||||
f"({max_frames * frame_bytes / GiB:.1f} GiB)")
|
||||
print()
|
||||
print(" Note: no swap on this machine — exceeding available RAM")
|
||||
print(" will trigger the OOM killer. Stay within the budget.")
|
||||
|
||||
def centers(cv):
|
||||
"""Set of ``(x, y)`` integer cluster centres from a ClusterVector."""
|
||||
if cv.size == 0:
|
||||
return set()
|
||||
a = np.asarray(cv)
|
||||
return {(int(x), int(y)) for x, y in zip(a["x"], a["y"])}
|
||||
|
||||
|
||||
def only_sets(cpu_c, cu_c, tol=1):
|
||||
"""CPU-only / CUDA-only centres, ignoring ``<=tol`` px 'shifted' matches.
|
||||
|
||||
``tol=0`` returns the exact set difference (shifted twins stay counted as
|
||||
mismatches); ``tol=1`` drops any mismatch that has a counterpart in the
|
||||
other finder's 8-neighbourhood.
|
||||
"""
|
||||
def near(p, other):
|
||||
x, y = p
|
||||
return any((x + dx, y + dy) in other
|
||||
for dx in range(-tol, tol + 1)
|
||||
for dy in range(-tol, tol + 1))
|
||||
|
||||
cpu_only = {p for p in cpu_c - cu_c if not near(p, cu_c)}
|
||||
cu_only = {p for p in cu_c - cpu_c if not near(p, cpu_c)}
|
||||
return cpu_only, cu_only
|
||||
|
||||
|
||||
def footprint_mask(cs, shape, rx, ry):
|
||||
"""Boolean map: True where a pixel lies in some cluster's footprint.
|
||||
|
||||
Each centre paints a ``(2*ry+1) x (2*rx+1)`` box, clipped at the borders.
|
||||
"""
|
||||
m = np.zeros(shape, bool)
|
||||
for (cx, cy) in cs:
|
||||
m[max(0, cy - ry):cy + ry + 1, max(0, cx - rx):cx + rx + 1] = True
|
||||
return m
|
||||
|
||||
|
||||
def shift_dist(p, other, R=4):
|
||||
"""Chebyshev distance from ``p`` to the nearest member of ``other``.
|
||||
|
||||
Rings are searched inner-to-outer so the first hit is the nearest.
|
||||
Returns ``-1`` if nothing lies within ``R`` px.
|
||||
"""
|
||||
x, y = p
|
||||
for r in range(1, R + 1):
|
||||
for dx in range(-r, r + 1):
|
||||
for dy in range(-r, r + 1):
|
||||
if max(abs(dx), abs(dy)) == r and (x + dx, y + dy) in other:
|
||||
return r
|
||||
return -1
|
||||
|
||||
|
||||
# --------------------------------------------------------------------------- #
|
||||
# Reusable scan + plot
|
||||
# --------------------------------------------------------------------------- #
|
||||
def train_pedestal(finders, f, n_frames, seek=0):
|
||||
"""Push the first ``n_frames`` frames of ``f`` into every finder given."""
|
||||
f.seek(seek)
|
||||
for _ in range(n_frames):
|
||||
img = f.read_frame().copy()
|
||||
for cf in finders:
|
||||
cf.push_pedestal_frame(img)
|
||||
|
||||
|
||||
def scan_mismatches(cf_cpu, cf_cuda, data, rx, ry,
|
||||
scan_count=1000, n_show=8, tol=0):
|
||||
"""Run both finders over ``scan_count`` frames sampled across ``data``.
|
||||
|
||||
Before every frame it snapshots the pedestal each finder will DECIDE with
|
||||
(CPU: host mean/rms; CUDA: device mean/rms on stream 0), so a later
|
||||
recompute uses the exact decision-time baseline.
|
||||
|
||||
``tol`` is forwarded to :func:`only_sets` when scoring a frame; ``tol=0``
|
||||
keeps shifted twins as mismatches. Returns ``(show, totals)`` where
|
||||
``show`` is the ``n_show`` frames with the most mismatches, each a dict
|
||||
with ``fid``, ``score``, ``cpu_c``, ``cu_c``, ``mism`` and the four
|
||||
pedestal arrays ``ped_cpu/noise_cpu/ped_cu/noise_cu``.
|
||||
"""
|
||||
show = []
|
||||
tot_cpu_only = tot_cu_only = 0
|
||||
for fid in np.linspace(0, len(data) - 1, scan_count, dtype=int):
|
||||
snap = dict(ped_cpu=np.asarray(cf_cpu.pedestal).copy(),
|
||||
noise_cpu=np.asarray(cf_cpu.noise).copy(),
|
||||
ped_cu=np.asarray(cf_cuda.device_pedestal(0)).copy(),
|
||||
noise_cu=np.asarray(cf_cuda.device_noise(0)).copy())
|
||||
cf_cpu.find_clusters(data[fid])
|
||||
cpu_c = centers(cf_cpu.steal_clusters(realloc_same_capacity=True))
|
||||
cf_cuda.find_clusters(data[fid])
|
||||
cu_c = centers(cf_cuda.steal_clusters(realloc_same_capacity=True))
|
||||
|
||||
cpu_only, cu_only = only_sets(cpu_c, cu_c, tol=tol)
|
||||
tot_cpu_only += len(cpu_only)
|
||||
tot_cu_only += len(cu_only)
|
||||
mism = cpu_only | cu_only
|
||||
if not mism:
|
||||
continue
|
||||
if len(show) < n_show or len(mism) > show[-1]["score"]:
|
||||
show.append(dict(score=len(mism), fid=int(fid),
|
||||
cpu_c=cpu_c, cu_c=cu_c, mism=mism, **snap))
|
||||
show.sort(key=lambda e: -e["score"])
|
||||
del show[n_show:]
|
||||
return show, dict(cpu_only=tot_cpu_only, cu_only=tot_cu_only)
|
||||
|
||||
|
||||
def compare_finders(finders, data, scan_count=1000, n_bins=200, e_range=(-2, 4000)):
|
||||
"""Run several finders over the same frames and score pairwise agreement.
|
||||
|
||||
``finders`` is a dict ``{name: finder}``; every finder must already be
|
||||
trained on the SAME pedestal frames. Each is run over the same
|
||||
``scan_count`` frames sampled across ``data`` (identical ``find_clusters``/
|
||||
``steal_clusters`` API for CPU, frozen-CPU and CUDA). In the same pass it
|
||||
accumulates a per-finder cluster-energy histogram (from ``cv.sum()``), so no
|
||||
extra scan is needed to draw the spectra.
|
||||
|
||||
Returns ``(totals, pairs, frames_scanned, hists)``:
|
||||
* ``totals[name]`` total clusters found by that finder
|
||||
* ``pairs[(a, b)]`` dict with ``a_only``/``b_only``/``mismatch``
|
||||
summed over frames (exact, tol=0)
|
||||
* ``hists[name]`` boost ``Histogram`` of cluster energies
|
||||
Companions :func:`print_comparison` and :func:`plot_spectra` render these.
|
||||
"""
|
||||
names = list(finders)
|
||||
totals = {n: 0 for n in names}
|
||||
hists = {n: bh.Histogram(bh.axis.Regular(n_bins, *e_range)) for n in names}
|
||||
pairs = {}
|
||||
for i, a in enumerate(names):
|
||||
for b in names[i + 1:]:
|
||||
pairs[(a, b)] = dict(a_only=0, b_only=0, mismatch=0)
|
||||
|
||||
fids = np.linspace(0, len(data) - 1, scan_count, dtype=int)
|
||||
for fid in fids:
|
||||
cs = {}
|
||||
for n, cf in finders.items():
|
||||
cf.find_clusters(data[fid])
|
||||
cv = cf.steal_clusters(realloc_same_capacity=True)
|
||||
cs[n] = centers(cv)
|
||||
totals[n] += len(cs[n])
|
||||
if cv.size:
|
||||
hists[n].fill(np.asarray(cv.sum()).ravel())
|
||||
for (a, b), acc in pairs.items():
|
||||
a_only, b_only = only_sets(cs[a], cs[b], tol=0)
|
||||
acc["a_only"] += len(a_only)
|
||||
acc["b_only"] += len(b_only)
|
||||
acc["mismatch"] += len(a_only) + len(b_only)
|
||||
return totals, pairs, len(fids), hists
|
||||
|
||||
|
||||
def plot_spectra(hists, totals=None, title="Cluster energy spectrum"):
|
||||
"""Overlay per-finder cluster-energy spectra with a ratio panel.
|
||||
|
||||
``hists`` is the ``{name: Histogram}`` returned by :func:`compare_finders`;
|
||||
the first finder is the reference for the ratio panel. Returns the Figure.
|
||||
"""
|
||||
names = list(hists)
|
||||
ref = names[0]
|
||||
edges = hists[ref].axes[0].edges
|
||||
vals = {n: hists[n].values() for n in names}
|
||||
|
||||
fig, (ax_spec, ax_ratio) = plt.subplots(
|
||||
2, 1, figsize=(8, 6), sharex=True,
|
||||
gridspec_kw={"height_ratios": [3, 1]})
|
||||
|
||||
styles = ["-", "--", "-.", ":"]
|
||||
for i, n in enumerate(names):
|
||||
lbl = n if totals is None else f"{n} ({totals[n]:,} clusters)"
|
||||
ax_spec.stairs(vals[n], edges, label=lbl, linestyle=styles[i % len(styles)])
|
||||
ax_spec.set_ylabel("Counts")
|
||||
ax_spec.set_title(title)
|
||||
ax_spec.legend()
|
||||
ax_spec.grid(alpha=0.2)
|
||||
|
||||
with np.errstate(divide="ignore", invalid="ignore"):
|
||||
for i, n in enumerate(names[1:], start=1):
|
||||
ratio = np.where(vals[ref] > 0, vals[n] / vals[ref], np.nan)
|
||||
ax_ratio.stairs(ratio, edges, label=f"{n} / {ref}",
|
||||
color=f"C{i}", linestyle=styles[i % len(styles)])
|
||||
ax_ratio.axhline(1.0, color="gray", linewidth=0.5)
|
||||
ax_ratio.set_ylabel(f"/ {ref}")
|
||||
ax_ratio.set_xlabel("Energy [ADU]")
|
||||
ax_ratio.set_ylim(0.5, 2.0)
|
||||
ax_ratio.legend(fontsize=8)
|
||||
ax_ratio.grid(alpha=0.3)
|
||||
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
return fig
|
||||
|
||||
|
||||
def print_comparison(totals, pairs, frames_scanned):
|
||||
"""Pretty-print the output of :func:`compare_finders`."""
|
||||
print(f"Scanned {frames_scanned} frames\n")
|
||||
print("Total clusters per finder:")
|
||||
for n, t in totals.items():
|
||||
print(f" {n:<16} {t:>12,}")
|
||||
print("\nPairwise exact mismatches (tol=0):")
|
||||
print(f" {'pair':<28} {'A-only':>10} {'B-only':>10} {'total':>10}")
|
||||
for (a, b), acc in pairs.items():
|
||||
ref = max(totals[a], totals[b], 1)
|
||||
pct = 100.0 * acc["mismatch"] / ref
|
||||
print(f" {a+' vs '+b:<28} {acc['a_only']:>10,} "
|
||||
f"{acc['b_only']:>10,} {acc['mismatch']:>10,} ({pct:.4f}%)")
|
||||
|
||||
|
||||
def walkthrough(show, data, rx, ry, rows, cols, n_sigma, pick=0,
|
||||
labels=('A', 'B')):
|
||||
"""Manual Test1/Test3 recompute of the strongest residual mismatch.
|
||||
|
||||
Dissects ``show[pick]`` (from :func:`scan_mismatches`, run as ``(a, b)``)
|
||||
under each finder's decision-time snapshot pedestal — ``a`` = ``ped_cpu`` /
|
||||
``noise_cpu``, ``b`` = ``ped_cu`` / ``noise_cu`` — and prints the raw and
|
||||
pedestal-subtracted window plus the accept/reject each finder reaches. With
|
||||
the double/double build the recompute reproduces the kernel exactly, so a
|
||||
lone surviving mismatch (e.g. the single 7x7 residual) is fully explained:
|
||||
the test that flips (Test1/Test3) and the pedestal gap name the cause.
|
||||
"""
|
||||
if not show:
|
||||
print("No residual mismatches to walk through.")
|
||||
return
|
||||
e = show[pick]
|
||||
sx, sy = 2 * rx + 1, 2 * ry + 1
|
||||
c3 = np.sqrt(sx * sy)
|
||||
frame = data[e['fid']].astype(np.float64)
|
||||
sub_a = frame - e['ped_cpu']
|
||||
sub_b = frame - e['ped_cu']
|
||||
|
||||
# strongest mismatch pixel whose full window stays inside the frame
|
||||
cand = [p for p in e['mism']
|
||||
if rx <= p[0] < cols - rx and ry <= p[1] < rows - ry]
|
||||
if not cand:
|
||||
print(f"frame {e['fid']}: all mismatches on the border — pick another.")
|
||||
return
|
||||
X0, Y0 = max(cand, key=lambda p: max(sub_a[p[1], p[0]], sub_b[p[1], p[0]]))
|
||||
|
||||
def evaluate(mean, rms_img):
|
||||
sig = (frame[Y0 - ry:Y0 + ry + 1, X0 - rx:X0 + rx + 1]
|
||||
- mean[Y0 - ry:Y0 + ry + 1, X0 - rx:X0 + rx + 1])
|
||||
value = frame[Y0, X0] - mean[Y0, X0]
|
||||
rms = rms_img[Y0, X0]
|
||||
thr1, thr3 = n_sigma * rms, c3 * n_sigma * rms
|
||||
mx, total = sig.max(), sig.sum()
|
||||
localmax = bool(value >= mx)
|
||||
t1, t3 = bool(mx > thr1), bool(total > thr3)
|
||||
accept = bool(value >= -thr1) and localmax and (t1 or t3)
|
||||
return dict(sig=sig, value=value, rms=rms, thr1=thr1, thr3=thr3,
|
||||
mx=mx, total=total, localmax=localmax, t1=t1, t3=t3,
|
||||
accept=accept)
|
||||
|
||||
la, lb = labels
|
||||
owner = la if (X0, Y0) in e['cpu_c'] else lb
|
||||
print(f"frame {e['fid']} centre (x={X0}, y={Y0}) accepted only by {owner}")
|
||||
print("raw window (ADU):")
|
||||
print(np.array2string(frame[Y0 - ry:Y0 + ry + 1,
|
||||
X0 - rx:X0 + rx + 1].astype(int)))
|
||||
print()
|
||||
|
||||
ra = evaluate(e['ped_cpu'], e['noise_cpu'])
|
||||
rb = evaluate(e['ped_cu'], e['noise_cu'])
|
||||
|
||||
def detail(tag, r):
|
||||
print(f"--- {tag} ---")
|
||||
print(" subtracted window:")
|
||||
print(" ", np.array2string(r['sig'], precision=1, prefix=' '))
|
||||
print(f" centre value = {float(r['value']):8.2f} "
|
||||
f"(local max? {r['localmax']})")
|
||||
print(f" max = {float(r['mx']):8.2f} "
|
||||
f"total = {float(r['total']):8.2f}")
|
||||
print(f" rms(centre) = {float(r['rms']):8.3f}")
|
||||
print(f" Test1: max > {n_sigma}*rms = {float(r['thr1']):8.2f} -> {r['t1']}")
|
||||
print(f" Test3: total > {c3:.0f}*{n_sigma}*rms = {float(r['thr3']):8.2f} -> {r['t3']}")
|
||||
print(f" ACCEPT = {r['accept']}")
|
||||
print()
|
||||
|
||||
detail(la, ra)
|
||||
detail(lb, rb)
|
||||
print(f"RESULT: {la} = {'ACCEPT' if ra['accept'] else 'reject'} , "
|
||||
f"{lb} = {'ACCEPT' if rb['accept'] else 'reject'}")
|
||||
gap = abs(float(e['ped_cpu'][Y0, X0] - e['ped_cu'][Y0, X0]))
|
||||
rgap = abs(float(e['noise_cpu'][Y0, X0] - e['noise_cu'][Y0, X0]))
|
||||
print(f"pedestal mean gap @centre = {gap:.4f} ADU; rms gap = {rgap:.4f}")
|
||||
if ra['accept'] == rb['accept']:
|
||||
print("NOTE: recompute agrees for this centre pixel — the split is FP "
|
||||
"rounding right at threshold or a window-neighbour effect; "
|
||||
"try another pick.")
|
||||
|
||||
|
||||
def plot_masked_mismatch(show, data, rx, ry, rows, cols,
|
||||
zoom=30, show_vals=True):
|
||||
"""Side-by-side masked view of each frame in ``show``.
|
||||
|
||||
Cluster pixels are coloured by pedestal-subtracted value (value printed
|
||||
when ``show_vals``); non-cluster pixels are white; a red dot marks every
|
||||
cluster centre. The cut is centred on the strongest mismatch pixel of the
|
||||
frame. No box, no tolerance discarding — a shifted twin in the other
|
||||
finder stays visible. Returns the Figure.
|
||||
"""
|
||||
cmap = plt.cm.viridis.copy()
|
||||
cmap.set_bad("white")
|
||||
fig, axes = plt.subplots(len(show), 2, figsize=(11, 5.3 * len(show)),
|
||||
squeeze=False)
|
||||
for r, e in enumerate(show):
|
||||
fid = e["fid"]
|
||||
sub_cpu = data[fid].astype(np.float64) - e["ped_cpu"]
|
||||
sub_cu = data[fid].astype(np.float64) - e["ped_cu"]
|
||||
X0, Y0 = max(e["mism"],
|
||||
key=lambda p: max(sub_cpu[p[1], p[0]], sub_cu[p[1], p[0]]))
|
||||
half = zoom // 2
|
||||
r0, r1 = max(0, Y0 - half), min(rows, Y0 + half)
|
||||
c0, c1 = max(0, X0 - half), min(cols, X0 + half)
|
||||
win_cpu, win_cu = sub_cpu[r0:r1, c0:c1], sub_cu[r0:r1, c0:c1]
|
||||
mask_cpu = footprint_mask(e["cpu_c"], sub_cpu.shape, rx, ry)[r0:r1, c0:c1]
|
||||
mask_cu = footprint_mask(e["cu_c"], sub_cu.shape, rx, ry)[r0:r1, c0:c1]
|
||||
union = mask_cpu | mask_cu
|
||||
vmax = max(np.percentile(np.concatenate([win_cpu[union], win_cu[union]]),
|
||||
99) if union.any() else 50.0, 50.0)
|
||||
|
||||
for ax, win, mask, cs, name in [
|
||||
(axes[r][0], win_cpu, mask_cpu, e["cpu_c"], "CPU"),
|
||||
(axes[r][1], win_cu, mask_cu, e["cu_c"], "CUDA")]:
|
||||
ax.imshow(np.ma.masked_where(~mask, win), cmap=cmap, vmin=0,
|
||||
vmax=vmax, interpolation="nearest")
|
||||
for (cx, cy) in cs:
|
||||
if c0 <= cx < c1 and r0 <= cy < r1:
|
||||
ax.plot(cx - c0, cy - r0, ".", color="red", ms=7)
|
||||
if show_vals:
|
||||
for i in range(win.shape[0]):
|
||||
for j in range(win.shape[1]):
|
||||
if mask[i, j]:
|
||||
ax.text(j, i, f"{win[i, j]:.0f}", ha="center",
|
||||
va="center", fontsize=6,
|
||||
color="white" if win[i, j] < 0.55 * vmax
|
||||
else "black")
|
||||
ax.set_title(f"frame {fid} — {name}: {len(cs)} clusters", fontsize=10)
|
||||
ax.set_xticks([])
|
||||
ax.set_yticks([])
|
||||
axes[r][0].set_ylabel(f"{e['score']} mismatches\nzoom @ ({X0},{Y0})",
|
||||
fontsize=9)
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
return fig
|
||||
Reference in new issue
Block a user