diff --git a/files/cdf_audit.md b/files/cdf_audit.md deleted file mode 100644 index 6811fc548..000000000 --- a/files/cdf_audit.md +++ /dev/null @@ -1,394 +0,0 @@ -# Adaptive Rectangular CDF — Phase 1 Audit - -Reference document for issue [#322](https://github.com/PyAutoLabs/PyAutoArray/issues/322). -Records exactly how the current JAX-compatible adaptive rectangular source-plane -mesh works, so that Phase 2+ can extend the density signal without touching the -spline machinery. - -This is a developer note. It lives under `PyAutoArray/files/` to keep package -namespaces clean. All paths are relative to the PyAutoArray repo root unless -explicitly absolute. - ---- - -## 1. Class hierarchy - -Two adaptive variants live in `autoarray/inversion/mesh/mesh/`: - -``` -RectangularAdaptDensity (rectangular_adapt_density.py:64) -└── RectangularAdaptImage (rectangular_adapt_image.py:9) - └── RectangularSplineAdaptImage (rectangular_spline_adapt_image.py:16) -└── RectangularSplineAdaptDensity (rectangular_spline_adapt_density.py:16) -``` - -Naming caveat: `RectangularAdaptDensity` is *not* an adaptive-coordinate mesh. -Its `interpolator_from` returns a plain `InterpolatorRectangular` whose -`_mappings_sizes_weights` runs the linear-CDF helper -(`adaptive_rectangular_mappings_weights_via_interpolation_from`, in -`autoarray/inversion/mesh/interpolator/rectangular.py:113`). It is "adaptive" -only in the implicit sense — the linear empirical CDF *is* the warp; there is -no separate uniform path. The `RectangularSpline*` subclasses swap that linear -CDF for the polynomial+Hermite-spline implementation in `rectangular_spline.py`. - -The class hierarchy carries one switching flag downstream: - -- `interpolator_cls` (property) — `InterpolatorRectangular` vs - `InterpolatorRectangularSpline`. -- `MeshGeometryRectangular.spline_deg` (`mesh_geometry/abstract.py:22`) — when - non-None, the geometry's `areas_transformed` and `edges_transformed` - properties route to the spline helpers instead of the linear helpers - (`mesh_geometry/rectangular.py:478,526`). - ---- - -## 2. `mesh_weight_map` data flow - -`mesh_weight_map` is the single per-point weight signal that drives both -the linear and spline CDF transforms. It is shape `(N_base,)` — one scalar per -point in the (non-oversampled) source-plane data grid. - -### 2.1 Where it is produced - -| Mesh class | `mesh_weight_map_from(adapt_data, xp)` returns | -|----------------------------------|------------------------------------------------| -| `RectangularAdaptDensity` | `None` (rectangular_adapt_density.py:159) | -| `RectangularAdaptImage` | normalised `adapt_data` (rectangular_adapt_image.py:82) | -| `RectangularSplineAdaptDensity` | inherited → `None` | -| `RectangularSplineAdaptImage` | inherited from `RectangularAdaptImage` | - -The `RectangularAdaptImage` recipe at `rectangular_adapt_image.py:82-107`: - -```python -mesh_weight_map = adapt_data.array -mesh_weight_map = xp.clip(mesh_weight_map, 1e-12, None) -mesh_weight_map = mesh_weight_map ** self.weight_power -mesh_weight_map = xp.where( - mesh_weight_map < self.weight_floor, - self.weight_floor, - mesh_weight_map, -) -mesh_weight_map = mesh_weight_map / xp.sum(mesh_weight_map) # final normalisation -``` - -Key observations: - -- `weight_floor` is applied *after* the power, *before* normalisation. -- The output sums to 1.0. Downstream uses interpret it as a per-point - probability mass (CDF construction via `cumsum` — see §3). -- `adapt_data` is the user-supplied adapt image, sampled at the *base* data - grid (one value per base grid point), not the oversampled grid. - -### 2.2 Where it is consumed - -All four entry points (linear / spline × areas / mappings) take -`mesh_weight_map` as a keyword and pass it unchanged into `create_transforms*`: - -| Call site | Function | Line | -|-----------------------------------------------------|---------------------------------------------------|------| -| `InterpolatorRectangular._mappings_sizes_weights` | `adaptive_rectangular_mappings_weights_via_interpolation_from` | rectangular.py:297 | -| `InterpolatorRectangularSpline._mappings_sizes_weights` | `adaptive_rectangular_mappings_weights_via_interpolation_from_spline` | rectangular_spline.py:618 | -| `MeshGeometryRectangular.areas_transformed` | `adaptive_rectangular_areas_from[_spline]` | mesh_geometry/rectangular.py:478-496 | -| `MeshGeometryRectangular.edges_transformed` | `adaptive_rectangular_transformed_grid_from[_spline]` | mesh_geometry/rectangular.py:526-548 | - -The interpolator and the mesh geometry independently build their own -transforms (no caching across the two). The CDF is re-fit per `mesh_weight_map`. - ---- - -## 3. The CDF transform itself - -### 3.1 Linear CDF (`rectangular.py:70`) - -`create_transforms(traced_points, mesh_weight_map=None, xp=np)`: - -```python -N = traced_points.shape[0] - -if mesh_weight_map is None: - t = xp.arange(1, N + 1) / (N + 1) # uniform rank probabilities - t = xp.stack([t, t], axis=1) # same t for y and x axes - sort_points = xp.sort(traced_points, axis=0) -else: - sdx = xp.argsort(traced_points, axis=0) # per-axis sort permutation - sort_points = xp.take_along_axis(traced_points, sdx, axis=0) - t = xp.stack([mesh_weight_map, mesh_weight_map], axis=1) # broadcast w to both axes - t = xp.take_along_axis(t, sdx, axis=0) # weights aligned to per-axis sort - t = xp.cumsum(t, axis=0) # weighted empirical CDF -``` - -This is the **separability assumption** in code form. `mesh_weight_map` is -the *same* (N,) weight vector applied independently to each axis, after -sorting traced points by that axis. The resulting CDF is the marginal CDF -along y for axis 0 and the marginal CDF along x for axis 1. - -If two source-plane points have identical x but different y, they will be -adjacent in the y-CDF and far apart in the x-CDF — same weight contribution -in each, but interpreted as different cumulative positions per axis. - -The forward / inverse transforms are then `np.interp` (numpy path) or a -vmapped `jnp.interp` (JAX path), producing piecewise-linear maps -`source plane ↔ [0,1]²`. - -### 3.2 Spline CDF (`rectangular_spline.py`) - -The spline variant replaces `np.interp` with: - -1. **Empirical CDF** (same as §3.1): per-axis weighted cumsum after sort. -2. **Chebyshev-node resample** (rectangular_spline.py:316-342): the empirical - `(t, sort_points)` is resampled at Chebyshev nodes via linear interp - (`_interp1d_numpy` / `_interp1d_jax`) to control Runge's phenomenon. Resample - width is `3 * deg` (default `33`). -3. **Weighted polynomial fit** (rectangular_spline.py:349-355): degree-`deg` - polynomial fit of `cx = f(cy)` with weights `1/dcx/dcy`. This gives a - closed-form **inverse CDF**: `polyval(coefs, y) = x`. -4. **Hermite-spline forward inverse** (rectangular_spline.py:122-247): - `InvertPolySpline` evaluates the polynomial on a `low_res = 20*deg` grid in - `[0,1]`, enforces strict monotonicity (cummax + per-position eps jitter), - then for any new `x` does binary search + cubic-Hermite interpolation to - recover `y`. - -The pair `(fwd_transform, rev_transform)` are C¹ continuous and differentiable. -`InvertPolySpline` is registered as a JAX pytree lazily so it can flow through -`jax.jit` / `jax.vmap` / autograd. - -**The weight signal enters exactly the same way** as in the linear variant — -through the `t = take_along_axis(stack([w, w]), sdx) → cumsum` step in -`_build_inv_poly_jax_impl` (rectangular_spline.py:390-399) and -`_build_inv_poly_numpy` (rectangular_spline.py:325-330). Everything after that -is a smoothing/interpolation pipeline that takes the empirical `(t, -sort_points)` and turns it into a differentiable map. - -### 3.3 Pre-CDF normalisation (the µ/σ step) - -Before CDF construction, every entry point normalises the source-plane grid: - -```python -mu = data_grid.mean(axis=0) # (2,) per-axis mean -scale = data_grid.std(axis=0).min() # scalar — min std over axes -source_grid_scaled = (data_grid - mu) / scale # zero-mean, isotropic-scale -``` - -This is a *single shared scale*, not per-axis. Anisotropic source-plane -distributions therefore retain their aspect ratio. The same `(mu, scale)` pair -is used to scale the over-sampled grid before pushing it through the -forward CDF (rectangular_spline.py:519-528, rectangular.py:174-186). - -The transform is fit on the *base* data grid; the oversampled grid is only -*evaluated* through the transform — never used in CDF construction. - ---- - -## 4. Index-space scatter and the `(N-3) * t + 1` rescale - -`adaptive_rectangular_mappings_weights_via_interpolation_from_spline` -(rectangular_spline.py:505-562) and its linear sibling do: - -```python -grid_over_sampled_transformed = transform(grid_over_sampled_scaled) - # shape (N_over, 2), values in [0, 1] -grid_over_index = (source_grid_size - 3) * grid_over_sampled_transformed + 1 - # values in [1, N - 2] -``` - -Then floor/ceil and bilinear weights: - -```python -ix_down = floor(grid_over_index[:, 0]) -ix_up = ceil(grid_over_index[:, 0]) -iy_down = floor(grid_over_index[:, 1]) -iy_up = ceil(grid_over_index[:, 1]) - -w_tl = (1 - t_row) * (1 - t_col) # t_row = (ix_down - ix_down) / 1 -w_tr = (1 - t_row) * t_col -w_bl = t_row * (1 - t_col) -w_br = t_row * t_col -``` - -Row-major flatten with a flip: - -```python -flat = (source_grid_size - idx[:, 0]) * source_grid_size + idx[:, 1] -``` - -Consequences: - -- Interior pixel indices live in `[1, N-2]`. With `N = source_grid_size`, the - outermost rows/cols (indices 0 and N-1 along each axis) are never the - bottom-left of any sample's bilinear cell. -- `zeroed_pixels` (`rectangular_adapt_density.py:126`) zeros the entire perimeter - ring via `rectangular_edge_pixel_list_from`. Combined with the `(N-3)+1` - buffer, **effective resolution is `(N-2) × (N-2)`**. - - N=33 → 31×31 = 961 effective pixels. - - N=50 → 48×48 = 2304 effective pixels. - - N=70 → 68×68 = 4624 effective pixels (≈ the "~4000+" the issue cites). -- Bilinear weights sum to 1 by construction; `sizes = 4` for every sample - (`rectangular_spline.py:628`). -- The `+ 1e-12` epsilon in the denominators of `t_row` / `t_col` - (rectangular.py:217-218, rectangular_spline.py:553-554) protects against - exact knot hits where `ix_down == ix_up`; under JAX this guards autograd too. - ---- - -## 5. Mesh-geometry side (areas / edges) - -`MeshGeometryRectangular.areas_transformed` computes per-pixel source-plane -areas by transforming uniform unit-square edges *through the inverse CDF* -(mesh_geometry/rectangular.py:466-496, rectangular_spline.py:476-502): - -```python -edges_y = xp.linspace(1, 0, N + 1) # (N+1,) unit-square y-edges -edges_x = xp.linspace(0, 1, N + 1) # (N+1,) unit-square x-edges -pixel_edges = inv_transform([edges_y, edges_x].T) * scale + mu -pixel_lengths = diff(pixel_edges, axis=0).squeeze() -dy, dx = pixel_lengths[:, 0], pixel_lengths[:, 1] -return abs(outer(dy, dx).flatten()) -``` - -This is **also separable**: pixel area is the outer product of per-axis edge -spacings. A non-separable density would not produce the correct per-pixel -area through this code path — the area would still be `dy_i * dx_j` even if -the density at `(y_i, x_j)` was uncorrelated with that product. Phase 4's -"is separability sufficient?" question lands here directly. - -`edges_transformed` is the pcolormesh path — it pushes the unit-square edges -through the inverse CDF and returns the warped (y, x) edge coordinates for -plotting. - ---- - -## 6. Where Phase 2 plugs in - -The composable density signal must arrive at `create_transforms_spline` as a -`mesh_weight_map` argument of shape `(N_base,)`, summed to 1.0, strictly -positive. **No change is needed below `mesh_weight_map_from`** — that method is -the only producer of the weight vector that every downstream call site reads -from. - -Two clean integration shapes: - -### 6.1 Subclass approach (preferred) - -Add `RectangularMultiComponentAdapt` next to `RectangularSplineAdaptImage`, -overriding only `mesh_weight_map_from(adapt_data, xp) -> (N_base,)`. The -override calls `compose_density(components, weights, floor, xp, context)` and -ends with the same `clip → normalise` finish as the parent. This keeps -`interpolator_from` and the spline interpolator class untouched. - -The `context` (magnification map, current reconstruction, residuals, etc.) -is the new wrinkle — it has to arrive at `mesh_weight_map_from`. Two -candidates: - -- Stash it on the mesh instance at construction time - (`RectangularMultiComponentAdapt(shape=..., context=...)`), making the mesh - carry the auxiliary fields. Cleanest at the call site; awkward if `context` - changes between fits (which it does — residuals update each iteration). -- Override `interpolator_from` to accept `context` as a kwarg and pass it - through. Slightly more plumbing but matches how `adapt_data` already flows - in the parent. - -Recommend Option B: keep the mesh stateless, pass `context` through -`interpolator_from`. The PyAutoGalaxy/Lens caller already has the auxiliary -fields at hand when it constructs the interpolator. - -### 6.2 Helper-injection approach (NOT recommended) - -Change `mesh_weight_map_from` to accept the multi-component arguments -directly. This pollutes the base-class signature and forces every existing -caller to know about the new fields. Rejected — keep the subclass. - ---- - -## 7. Open questions for Phase 4 (separability & low-rank) - -The current architecture commits to **separable per-axis CDFs**. For a single -density signal that's well-approximated by `ρ(x,y) = a(x) · b(y)` (Gaussian-like -sources, point-density of regularly traced rays), separability is fine. - -For multi-component densities the worry is: - -1. **Caustic proximity** is anisotropic along curves — the level sets are - 1-D loci in the source plane. The marginal x-CDF and y-CDF of a caustic - density will concentrate pixels along the x and y projections of the - caustic, *not* along the caustic itself. Pixels get wasted on the convex - hull of the projections rather than placed on the curve. -2. **Magnification density** is similarly anisotropic near critical curves. -3. **Source brightness × residual gradient** products may produce diagonal - or curved support that no separable factorisation captures. - -Phase 4 should measure how badly separability hurts vs a full 2-D CDF and -decide whether an outer-product factorisation -`ρ(x,y) ≈ a(x)·b(y) + Σ_k u_k(x)·v_k(y)` (sum of rank-1 outer products) is -worth the engineering cost. The current code has no obvious slot for a -non-separable CDF — would need a new transform class, not a new -`mesh_weight_map`. - -A separate concern: the **areas calculation** (§5) is fundamentally -separable in this code. Even if Phase 2 ships a non-separable weight signal, -the areas reported to magnification/regularisation downstream would still be -`dy_i · dx_j`. Investigating non-separability requires touching -`areas_transformed` *and* `edges_transformed` *and* the CDF construction -itself — a much bigger change than the Phase 2 subclass. - ---- - -## 8. Deferred to Phase 5: quantitative baseline - -The issue calls for baseline reconstruction χ², log-evidence, and peak -per-pixel residual at 500 / 1000 / 4000 source pixels on a reference dataset. -That requires a full lens-model fit (image → tracer → mapper → inversion → -fit), which is `autolens_workspace` / `autolens_workspace_developer` territory -— this issue scope ends at the library prototype + unit tests. - -Followup prompt: `autolens_workspace_developer/rectangular_adapt_cdf_benchmark.md`, -to be authored at the end of Phase 3 once the multi-component API is stable. -The benchmark script should: - -- Use one of the existing `autolens_workspace_test/scripts/imaging/` - pixelization examples as the reference dataset. -- Sweep N ∈ {25, 32, 40, 50, 70} (effective pixel counts 529, 900, 1444, - 2304, 4624) for both single-signal (current) and multi-component (Phase 3) - meshes. -- Hold sampler budget and regularisation strategy fixed. -- Compare against a Delaunay run at the matched compute budget. - ---- - -## 9. JAX compatibility checklist (what NOT to break) - -The current spline path is fully JAX/JIT-compatible. Phase 2/3 components -must preserve every property: - -- No dynamic shapes — every component returns `(N_base,)` regardless of - trigger conditions. Branches that "skip" a component should multiply its - contribution by 0, never resize. -- No scipy callbacks, no `pure_callback`, no kNN / Delaunay / RBF. -- No Python-level control flow on traced values. `jnp.where` for masking. -- Strictly positive output (`xp.clip(w, eps, None)` before normalisation). -- Finite under autograd — the `1e-12` epsilon trick in the bilinear weights is - load-bearing; replicate when introducing new divisions. -- `InvertPolySpline` pytree registration is *lazy* (rectangular_spline.py:163-204); - composing new components that wrap it must trigger - `InvertPolySpline._register_pytree()` themselves or rely on the existing - trigger in `_build_inv_poly_jax`. - -Library unit tests stay numpy-only per -`feedback_no_jax_in_unit_tests` — cross-xp parity goes in -`autolens_workspace_test`. - ---- - -## 10. Files touched in subsequent phases - -- New: `autoarray/inversion/mesh/interpolator/density_components.py` — - `compose_density` + per-component callables. -- New: `test_autoarray/inversion/pixelization/interpolator/test_density_components.py` - — numpy-only unit tests. -- (Phase 6, conditional) New: - `autoarray/inversion/mesh/mesh/rectangular_multi_component_adapt.py` and a - parity test in `test_rectangular_spline.py`. - -No edits to `rectangular_spline.py`, `rectangular.py`, or the existing mesh -classes are required for Phases 2–3. Phase 4's non-separable investigation -might force a deeper change, but only if separability is shown to be the -binding constraint. diff --git a/files/ghost_peak.png b/files/ghost_peak.png deleted file mode 100644 index fc6611139..000000000 Binary files a/files/ghost_peak.png and /dev/null differ diff --git a/files/ghost_peak_experiment.py b/files/ghost_peak_experiment.py deleted file mode 100644 index 4076f87e9..000000000 --- a/files/ghost_peak_experiment.py +++ /dev/null @@ -1,286 +0,0 @@ -""" -Ghost-peak diagnostic for the separable per-axis CDF in -``adaptive_rectangular_*`` (rectangular_spline.py). - -Hypothesis under test ---------------------- -A separable per-axis CDF cannot resolve a multi-modal adapt image (e.g. a -clumpy source with K bright peaks) without producing K**2 hot zones in the -warped mesh — K**2 - K of which are "ghost peaks" at the cross-product of -marginal modes where no source brightness exists. - -Method ------- -1. Sample N traced points uniformly across the source plane (proxy for a - uniformly back-projected image-plane mask — no lensing physics needed - for the geometric diagnostic). -2. Build a synthetic adapt image as a sum of K isotropic Gaussian peaks of - matched amplitude. -3. Compute ``mesh_weight_map`` exactly as ``RectangularRTUAdaptImage`` does: - ``clip(b, eps, None) ** power`` then normalise. -4. Build the spline CDF transforms via ``create_transforms_spline`` (the same - function ``RectangularSplineAdaptImage`` uses). -5. Evaluate the inverse CDF on the rectangular unit-square edge grid to - recover the warped mesh edges in source-plane coordinates. Because the - CDF is separable, the warped 2D mesh is the outer product of warped y - and warped x edges. -6. Count pixel centres (mid-edges) inside each of the K**2 candidate hot - zones (each peak location and each ghost cross-term). -7. Report the per-zone pixel-density ratios. - -Outputs -------- -- ``files/ghost_peak.png`` — overlay of synthetic brightness, traced points, - and the warped mesh edges, for K = 1, 2, 3 peaks side by side. -- console log: per-zone pixel-density ratios and ghost-to-real ratio. - -Run from the workspace root (or a task-bundle root) with its ``activate.sh`` -sourced:: - - source ./activate.sh - python PyAutoArray/files/ghost_peak_experiment.py -""" -import os -import sys - -import numpy as np -import matplotlib - -matplotlib.use("Agg") -import matplotlib.pyplot as plt - -from autoarray.inversion.mesh.interpolator.rectangular_spline import ( - SPLINE_CDF_DEFAULT_DEG, - create_transforms_spline, -) - - -# --------------------------------------------------------------------------- -# Synthetic source -# --------------------------------------------------------------------------- - - -def sample_traced_points(n: int, seed: int = 0) -> np.ndarray: - """N points uniformly in [-1, 1] x [-1, 1].""" - rng = np.random.default_rng(seed) - return rng.uniform(-1.0, 1.0, size=(n, 2)) - - -def brightness_at(points: np.ndarray, peaks, sigma: float = 0.1) -> np.ndarray: - """Sum of isotropic Gaussian peaks, amplitude 1 each. ``peaks`` is - a list of (y, x) tuples in the source-plane coordinate system.""" - b = np.zeros(points.shape[0]) - for (py, px) in peaks: - dy = points[:, 0] - py - dx = points[:, 1] - px - b += np.exp(-(dy ** 2 + dx ** 2) / (2 * sigma ** 2)) - return b - - -# --------------------------------------------------------------------------- -# Run the existing pipeline geometry -# --------------------------------------------------------------------------- - - -def warped_mesh_edges( - traced: np.ndarray, - mesh_weight_map: np.ndarray, - mesh_size: int, - deg: int = SPLINE_CDF_DEFAULT_DEG, -): - """Replicates the geometry path used by the live mesh. - - Mirrors ``adaptive_rectangular_transformed_grid_from_spline`` — - same (mu, scale) normalisation, same per-axis spline CDF. - - Returns (warped_y_edges, warped_x_edges), each shape (mesh_size + 1,), - in source-plane coordinates. - """ - mu = traced.mean(axis=0) - scale = traced.std(axis=0).min() - traced_scaled = (traced - mu) / scale - - _, inv_transform = create_transforms_spline( - traced_scaled, mesh_weight_map=mesh_weight_map, deg=deg, xp=np - ) - - # Per-axis unit-square edges — same as MeshGeometryRectangular.edges_transformed - edges_y_unit = np.linspace(1.0, 0.0, mesh_size + 1) - edges_x_unit = np.linspace(0.0, 1.0, mesh_size + 1) - - # The inverse transform takes (K, 2) -> (K, 2). Each column is processed - # independently (column 0 = y, column 1 = x). To recover the warped 2D - # mesh under the separability assumption, we feed the per-axis edges as - # rows: row i is (y_unit[i], x_unit[i]) and we take col 0 of the output - # for the warped-y edges and col 1 for the warped-x edges. - grid = np.stack([edges_y_unit, edges_x_unit], axis=1) - warped = inv_transform(grid) * scale + mu # (mesh_size + 1, 2) - - return warped[:, 0], warped[:, 1] - - -# --------------------------------------------------------------------------- -# Quantify per-zone pixel density -# --------------------------------------------------------------------------- - - -def pixel_centres(warped_y_edges, warped_x_edges): - """Cell centres are the outer product of consecutive-edge midpoints.""" - y_centres = 0.5 * (warped_y_edges[:-1] + warped_y_edges[1:]) - x_centres = 0.5 * (warped_x_edges[:-1] + warped_x_edges[1:]) - yy, xx = np.meshgrid(y_centres, x_centres, indexing="ij") - return np.stack([yy.ravel(), xx.ravel()], axis=1) - - -def count_in_disc(centres, py, px, radius): - d2 = (centres[:, 0] - py) ** 2 + (centres[:, 1] - px) ** 2 - return int(np.sum(d2 < radius ** 2)) - - -# --------------------------------------------------------------------------- -# Experiment runner -# --------------------------------------------------------------------------- - - -def run_experiment(peaks_label, peaks, real_zones, ghost_zones, - mesh_size=40, n_traced=2000, sigma=0.1, deg=11, - probe_radius=0.15): - traced = sample_traced_points(n_traced, seed=0) - brightness = brightness_at(traced, peaks=peaks, sigma=sigma) - # Mimic RectangularRTUAdaptImage.mesh_weight_map_from with power=1, floor=0 - weight_map = np.clip(brightness, 1e-12, None) - weight_map = weight_map / weight_map.sum() - - y_edges, x_edges = warped_mesh_edges( - traced, weight_map, mesh_size=mesh_size, deg=deg - ) - centres = pixel_centres(y_edges, x_edges) - - # Background pixel density per unit area (reference baseline) - bbox_area = (y_edges.max() - y_edges.min()) * (x_edges.max() - x_edges.min()) - bg_density = centres.shape[0] / bbox_area - probe_area = np.pi * probe_radius ** 2 - - print(f"\n=== {peaks_label}: K = {len(peaks)} peaks ===") - print( - f"Total pixel centres: {centres.shape[0]}, " - f"bbox area: {bbox_area:.3f}, baseline density: {bg_density:.2f} px/unit²" - ) - - real_counts = [] - for (py, px) in real_zones: - n = count_in_disc(centres, py, px, probe_radius) - d = n / probe_area - ratio = d / bg_density - real_counts.append(n) - print( - f" REAL peak at ({py:+.2f}, {px:+.2f}): " - f"{n:3d} centres ({d:6.1f} px/unit², {ratio:.2f}x baseline)" - ) - - ghost_counts = [] - for (py, px) in ghost_zones: - n = count_in_disc(centres, py, px, probe_radius) - d = n / probe_area - ratio = d / bg_density - ghost_counts.append(n) - print( - f" GHOST cross at ({py:+.2f}, {px:+.2f}): " - f"{n:3d} centres ({d:6.1f} px/unit², {ratio:.2f}x baseline)" - ) - - if real_counts and ghost_counts: - mean_real = np.mean(real_counts) - mean_ghost = np.mean(ghost_counts) - ghost_to_real = mean_ghost / max(mean_real, 1e-12) - print( - f" mean(real) = {mean_real:.1f}, mean(ghost) = {mean_ghost:.1f}, " - f"ghost / real = {ghost_to_real:.2f}" - ) - verdict = ( - "CONFIRMED ghost-peak failure (ghost density ≈ real density)" - if ghost_to_real > 0.5 - else "no ghost-peak failure (real ≫ ghost)" - ) - print(f" verdict: {verdict}") - - return traced, brightness, weight_map, y_edges, x_edges, centres - - -def plot_experiment(ax, label, traced, brightness, weight_map, y_edges, x_edges, - centres, real_zones, ghost_zones): - # Source brightness — scatter coloured by adapt image - sc = ax.scatter( - traced[:, 1], traced[:, 0], c=brightness, cmap="hot", - s=4, alpha=0.45, marker="." - ) - - # Warped mesh as gridlines - for ye in y_edges: - ax.plot([x_edges.min(), x_edges.max()], [ye, ye], color="cyan", - lw=0.25, alpha=0.5) - for xe in x_edges: - ax.plot([xe, xe], [y_edges.min(), y_edges.max()], color="cyan", - lw=0.25, alpha=0.5) - - # Pixel centres - ax.scatter(centres[:, 1], centres[:, 0], c="white", s=2, - edgecolors="black", linewidths=0.2) - - # Real zones in green, ghost zones in red - for (py, px) in real_zones: - ax.plot(px, py, "o", color="lime", markersize=12, - markeredgecolor="black", markeredgewidth=1.5) - for (py, px) in ghost_zones: - ax.plot(px, py, "x", color="red", markersize=12, markeredgewidth=2) - - ax.set_xlim(-1.05, 1.05) - ax.set_ylim(-1.05, 1.05) - ax.set_aspect("equal") - ax.set_title(label) - ax.set_xlabel("source x") - ax.set_ylabel("source y") - return sc - - -def main(): - # K = 1 control (no separability problem possible) - peaks_1 = [(0.0, 0.0)] - real_1 = [(0.0, 0.0)] - ghost_1 = [] # no cross-terms with K = 1 - - # K = 2 diagonal — the textbook ghost-peak case - peaks_2 = [(-0.5, -0.5), (0.5, 0.5)] - real_2 = [(-0.5, -0.5), (0.5, 0.5)] - ghost_2 = [(-0.5, 0.5), (0.5, -0.5)] - - # K = 3 — 9 candidate zones (3 real, 6 ghost) - peaks_3 = [(-0.5, -0.5), (0.0, 0.5), (0.5, -0.2)] - real_3 = list(peaks_3) - ghost_3 = [ - (py, px) for py in (-0.5, 0.0, 0.5) for px in (-0.5, 0.5, -0.2) - if (py, px) not in real_3 - ] - - cases = [ - ("K=1 (control)", peaks_1, real_1, ghost_1), - ("K=2 (diagonal)", peaks_2, real_2, ghost_2), - ("K=3 (offset)", peaks_3, real_3, ghost_3), - ] - - fig, axes = plt.subplots(1, 3, figsize=(18, 6)) - last_sc = None - for ax, (label, peaks, real, ghost) in zip(axes, cases): - result = run_experiment(label, peaks, real, ghost) - last_sc = plot_experiment(ax, label, *result, - real_zones=real, ghost_zones=ghost) - - fig.colorbar(last_sc, ax=axes, label="brightness (adapt image)", - fraction=0.025, pad=0.02) - out = os.path.join(os.path.dirname(__file__), "ghost_peak.png") - fig.savefig(out, dpi=130, bbox_inches="tight") - print(f"\nfigure written: {out}") - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/files/ghost_peak_findings.md b/files/ghost_peak_findings.md deleted file mode 100644 index bdbc8fe7f..000000000 --- a/files/ghost_peak_findings.md +++ /dev/null @@ -1,179 +0,0 @@ -# Ghost-Peak Findings — Phase 1.5 - -Empirical results from `files/ghost_peak_experiment.py`. Companion to -`files/cdf_audit.md`. Issue -[#322](https://github.com/PyAutoLabs/PyAutoArray/issues/322). - -## Hypothesis - -A separable per-axis CDF cannot adapt to multi-modal source planes (multiple -bright regions) without creating ghost peaks at the cross-product of marginal -modes. - -## Result — **CONFIRMED** - -| Configuration | mean(real) | mean(ghost) | ghost / real | -|---------------|-----------:|------------:|-------------:| -| K = 1 (control, one peak at origin) | 1006 px | — | n/a | -| K = 2 (diagonal: peaks at ±(0.5, 0.5)) | 226 px | 224 px | **0.99** | -| K = 3 (offset: 3 non-axis-aligned peaks) | 109 px | 105 px | **0.96** | - -Visual proof: `files/ghost_peak.png`. The K=2 panel shows a 2×2 grid of dense -mesh zones at all four (±0.5, ±0.5) corners — two coincide with the bright -source peaks (green circles), two land on empty background (red ×). The -warped mesh wastes **roughly half the pixel budget** on ghosts. - -The K=3 panel shows the corresponding 3×3 grid of dense zones, of which 3 are -real and 6 are ghosts — **two-thirds of the adaptive budget is wasted**. - -## Why it happens (mechanism) - -Per-axis CDF construction (`rectangular_spline.py:create_transforms_spline`) -sorts traced points independently along y and x. The two marginal CDFs -steepen at the y- and x-projections of every peak. When the inverse CDF is -evaluated on a uniform unit-square grid, both axes concentrate edges at the -projections — and the outer product places dense pixels at **all K² -combinations** of (y-projection × x-projection), only K of which are real -source peaks. - -This is a structural property of any separable transform; no choice of -`weight_power`, `weight_floor`, or `spline_deg` can avoid it. Composing -multiple physics signals (magnification + brightness + residuals + caustic) -into the same `mesh_weight_map` slot — what the original Phase 3 plan was -building — inherits the failure because all signals still funnel through the -same separable transform. - -## Implication for the project - -The multi-component framework shipped in -`autoarray/inversion/mesh/interpolator/density_components.py` (Phase 2: -`compose_density` + `uniform_density_component`) is still useful for -combining physics signals on **single-mode** sources, but it does **not** -address the user's stated goal of adapting to multiple source-plane regions -simultaneously. That goal requires either: - -1. A non-separable transform (breaks axis-aligned-rectangle pixel shape), or -2. Multiple meshes (one per source region) stitched together (breaks - single-mesh assumption), or -3. A rotation hack that handles the K=2 case (most common in practice). - -## Three concrete paths forward, smallest-first - -### Path A — PCA rotation hack (smallest change, partial fix) - -For any multi-modal brightness map, PCA on the brightness-weighted point -cloud gives a principal axis. Rotate the source plane into a frame where -that axis is x (or y), run the existing separable CDF, rotate the warped -mesh back. - -- **K = 2 (binary peaks)**: peaks become collinear on the principal axis → - marginal-y collapses (one mode), marginal-x has both modes → outer product - gives **only 2 dense zones, both real, zero ghosts**. -- **K = 3 colinear**: same — perfect fix. -- **K = 3 non-colinear**: PCA finds the dominant axis but the off-axis peak - still creates ghosts. Reduces ghost density but doesn't eliminate them. -- **K ≥ 4**: PCA helps proportionally less. - -**Effort**: small. A new `RectangularRotatedAdaptImage` subclass that wraps -`RectangularSplineAdaptImage` with pre/post rotation steps. - -**JAX**: trivial (rotation is a single matrix multiply, JAX-friendly). - -**Honest limitation**: helps the most common case (binary / colinear) but -isn't a complete fix. - -### Path B — Pre-segment + multi-sub-mesh (moderate change, full fix) - -Detect K modes in the adapt image (Gaussian mixture, k-means on -brightness-weighted points, or quantile thresholding). Per detected mode, -construct a separate `RectangularSplineAdaptImage` covering that mode's -bounding box. Combine the K sub-meshes into a `MultiRectangularAdapt` -container that presents a single mesh interface to the inversion pipeline. - -- **Each sub-mesh** uses unchanged CDF code. -- **No ghosts** within each sub-mesh because each is single-mode. -- **Coverage**: gaps between sub-meshes either remain unmodelled (acceptable - for high-contrast multi-modal sources) or get a low-resolution fallback - sub-mesh. - -**Effort**: moderate. New container class, segmentation logic (probably -JAX-incompatible — would run as a numpy preprocessing step before the -JIT-compiled likelihood), per-sub-mesh inversion integration. - -**JAX**: the per-sub-mesh inversions stay JAX-compatible; the segmentation -step (clustering / GMM fit) is once per fit, not per likelihood eval, so it -can be numpy-only. - -**Honest limitation**: requires choosing K (or detecting it). Segmentation -quality determines reconstruction quality. - -### Path C — Knothe-Rosenblatt with non-axis-aligned cells (largest change, principled) - -Replace the separable transform with a Knothe-Rosenblatt 2D transport: -y = F₁⁻¹(u), x = F₂⁻¹(v | y). Conditional x-CDF depends on y, so x-edges vary -per y-stripe. Cells become quadrilaterals (still 4 corners, still -combinatorially-rectangular neighbours, but no longer axis-aligned). - -- **No ghosts** — the transform follows the full joint density. -- **Topology**: rectangular *connectivity* preserved. Cell *geometry* - becomes general quadrilateral. -- **Areas calculation** (Phase 1 audit §5) — needs Jacobian instead of - `outer(dy, dx)`. Still tractable but ≠ outer product. -- **Bilinear scatter**: needs bilinear-on-quadrilateral, a bigger change to - `adaptive_rectangular_mappings_weights_via_interpolation_from_spline`. -- **Spline machinery**: the conditional CDF requires a 2D function — a - tensor-product spline or one polynomial per y-bin. - -**Effort**: large. New CDF, new bilinear, new areas, new mesh class. The -test surface roughly doubles. - -**JAX**: feasible, but the conditional CDF needs careful design to avoid -dynamic shapes. - -**Honest limitation**: the "fixed rectangular topology" phrase in the prompt -might mean axis-aligned rectangles, in which case Path C violates the -constraint. Need user confirmation before committing. - -## What I'd revert / keep from the work so far - -- **Keep**: Phase 1 audit (`cdf_audit.md`), Phase 2 framework - (`density_components.py` + tests, currently 9/9 passing), this findings doc - and the experiment script. -- **Drop**: the four physical-density factory functions - (`magnification_density`, `brightness_density`, `residual_density`, - `caustic_proximity_density`) — already reverted out of - `density_components.py`. They were aimed at the wrong target. - -## Recommendation - -Pick Path A first. It's small, it handles the most common multi-peak case -(binary / colinear sources) cleanly, and it provides a real-world reference -point for measuring how much Path B/C would buy on harder cases. We can -re-extend Path A → B → C as the user's actual source complexity demands, -without backing-out work each step. - -## Path A empirical follow-up — **WORKS** - -`files/pca_rotation_experiment.py` reproduces the ghost-peak diagnostic with -a brightness-weighted PCA pre-rotation. Compares baseline vs rotated for the -same three test cases. - -| Case | ghost/real (baseline) | ghost/real (rotated) | real-peak density | -|------|----------------------:|---------------------:|------------------:| -| K = 2 diagonal | 0.99 | **0.00** | +117% (226 → 490.5) | -| K = 2 axis-aligned | 0.00 | 0.00 | unchanged (467 → 453) | -| K = 3 offset (triangle) | 0.96 | **0.06** | +46% (109 → 159.3) | - -Visual proof: `files/pca_rotation.png`. The K = 2 diagonal mesh rotates by -+45.75° to align with the peaks — the cyan mesh grid hugs the diagonal -band between the real peaks, both ghost zones land in sparse regions. For -K = 3 the PCA axis runs through two of the three peaks; the third peak -still receives heavy pixel density and the off-axis ghost zones are nearly -empty. Even my pessimistic K=3 prediction in §"Path A" above was wrong by -an order of magnitude — non-colinear K=3 also gets a clean fix. - -Practical implication: a `RectangularRotatedAdaptImage` subclass adding -brightness-weighted PCA rotation in front of the existing CDF would -materially fix the multi-modal-source case for all realistic K. The cost -is a 2×2 matrix multiply on the input grid plus a rotation-back step on -mesh-geometry outputs — JAX-trivial. diff --git a/files/pca_rotation.png b/files/pca_rotation.png deleted file mode 100644 index adc576ac6..000000000 Binary files a/files/pca_rotation.png and /dev/null differ diff --git a/files/pca_rotation_experiment.py b/files/pca_rotation_experiment.py deleted file mode 100644 index 2bec5a6d6..000000000 --- a/files/pca_rotation_experiment.py +++ /dev/null @@ -1,298 +0,0 @@ -""" -Path-A diagnostic: does a PCA-aligned rotation of the source plane eliminate -the ghost-peak failure observed in ``ghost_peak_experiment.py``? - -Hypothesis ----------- -The separable per-axis CDF concentrates pixels at the outer product of -marginal modes. If a multi-modal source's peaks are aligned with the y/x -axes, there are no off-axis "ghost" combinations — all marginal modes line -up on the real peaks. PCA on the brightness-weighted point cloud finds the -principal axis; rotating the source plane into that frame makes K=2 peaks -exactly axis-aligned (and helps for higher K). - -Method ------- -1. Reuse the ghost-peak setup: synthetic 2-peak Gaussian brightness on - uniformly sampled traced points. -2. Compute PCA on traced points weighted by brightness. Rotate everything - into the principal frame. -3. Run the same spline-CDF + outer-product mesh as the baseline. -4. Quantify pixel density at the real peaks and the (now potentially - collapsed) ghost positions, in the rotated frame. -5. Visualise rotated-frame warped mesh AND source-frame warped mesh - (rotated back). - -Run from the workspace root (or a task-bundle root) with its ``activate.sh`` -sourced:: - - source ./activate.sh - python PyAutoArray/files/pca_rotation_experiment.py -""" -import os -import sys - -import numpy as np -import matplotlib - -matplotlib.use("Agg") -import matplotlib.pyplot as plt - -from autoarray.inversion.mesh.interpolator.rectangular_spline import ( - SPLINE_CDF_DEFAULT_DEG, - create_transforms_spline, -) - - -# --------------------------------------------------------------------------- -# Synthetic source (identical to ghost_peak_experiment.py) -# --------------------------------------------------------------------------- - - -def sample_traced_points(n: int, seed: int = 0) -> np.ndarray: - rng = np.random.default_rng(seed) - return rng.uniform(-1.0, 1.0, size=(n, 2)) - - -def brightness_at(points: np.ndarray, peaks, sigma: float = 0.1) -> np.ndarray: - b = np.zeros(points.shape[0]) - for (py, px) in peaks: - dy = points[:, 0] - py - dx = points[:, 1] - px - b += np.exp(-(dy ** 2 + dx ** 2) / (2 * sigma ** 2)) - return b - - -# --------------------------------------------------------------------------- -# Brightness-weighted PCA → rotation matrix -# --------------------------------------------------------------------------- - - -def pca_rotation_matrix(points: np.ndarray, weights: np.ndarray) -> np.ndarray: - """Return a 2x2 rotation matrix R such that R @ (point - centroid) is in - the principal-axis frame: column 0 = direction of max brightness-weighted - variance, column 1 = perpendicular. - - The matrix is constructed from the eigenvectors of the weighted covariance, - so applying it to (y, x) row vectors (via ``points @ R.T``) rotates them - into the PCA frame. - """ - w = weights / weights.sum() - centroid = (w[:, None] * points).sum(axis=0) - centred = points - centroid - cov = (w[:, None, None] * centred[:, :, None] * centred[:, None, :]).sum(axis=0) - - eigvals, eigvecs = np.linalg.eigh(cov) - # eigh returns ascending eigenvalues; we want the largest first. - order = np.argsort(eigvals)[::-1] - eigvecs = eigvecs[:, order] - # Make rotation deterministic by enforcing a sign convention. - if eigvecs[0, 0] < 0: - eigvecs[:, 0] = -eigvecs[:, 0] - if np.linalg.det(eigvecs) < 0: - eigvecs[:, 1] = -eigvecs[:, 1] - # Rows of R are the new basis vectors (so R @ x rotates x into PCA frame). - R = eigvecs.T - return R, centroid - - -def rotate_points(points: np.ndarray, R: np.ndarray, centroid: np.ndarray) -> np.ndarray: - return (points - centroid) @ R.T - - -def unrotate_points(points: np.ndarray, R: np.ndarray, centroid: np.ndarray) -> np.ndarray: - return points @ R + centroid - - -# --------------------------------------------------------------------------- -# Warped mesh edges in whatever frame the inputs live in -# --------------------------------------------------------------------------- - - -def warped_mesh_edges(traced, mesh_weight_map, mesh_size, deg=SPLINE_CDF_DEFAULT_DEG): - mu = traced.mean(axis=0) - scale = traced.std(axis=0).min() - traced_scaled = (traced - mu) / scale - - _, inv_transform = create_transforms_spline( - traced_scaled, mesh_weight_map=mesh_weight_map, deg=deg, xp=np - ) - - edges_y_unit = np.linspace(1.0, 0.0, mesh_size + 1) - edges_x_unit = np.linspace(0.0, 1.0, mesh_size + 1) - grid = np.stack([edges_y_unit, edges_x_unit], axis=1) - warped = inv_transform(grid) * scale + mu - - return warped[:, 0], warped[:, 1] - - -def pixel_centres_from_edges(y_edges, x_edges): - yc = 0.5 * (y_edges[:-1] + y_edges[1:]) - xc = 0.5 * (x_edges[:-1] + x_edges[1:]) - yy, xx = np.meshgrid(yc, xc, indexing="ij") - return np.stack([yy.ravel(), xx.ravel()], axis=1) - - -def count_in_disc(centres, py, px, radius): - d2 = (centres[:, 0] - py) ** 2 + (centres[:, 1] - px) ** 2 - return int(np.sum(d2 < radius ** 2)) - - -# --------------------------------------------------------------------------- -# Experiment -# --------------------------------------------------------------------------- - - -def run_one(label, peaks, mesh_size=40, n_traced=2000, sigma=0.1, - deg=11, probe_radius=0.15): - traced = sample_traced_points(n_traced, seed=0) - brightness = brightness_at(traced, peaks, sigma=sigma) - weight_map = np.clip(brightness, 1e-12, None) - weight_map = weight_map / weight_map.sum() - - # ---- Baseline (no rotation) ---- - yE_b, xE_b = warped_mesh_edges(traced, weight_map, mesh_size, deg) - centres_b = pixel_centres_from_edges(yE_b, xE_b) - - real_zones = list(peaks) - ghost_zones = [(py, px) for py, _ in peaks for _, px in peaks - if (py, px) not in real_zones] - - real_b = [count_in_disc(centres_b, py, px, probe_radius) for py, px in real_zones] - ghost_b = [count_in_disc(centres_b, py, px, probe_radius) for py, px in ghost_zones] - - # ---- PCA-rotated ---- - R, centroid = pca_rotation_matrix(traced, brightness) - traced_rot = rotate_points(traced, R, centroid) - yE_r, xE_r = warped_mesh_edges(traced_rot, weight_map, mesh_size, deg) - centres_r_rotframe = pixel_centres_from_edges(yE_r, xE_r) - centres_r_sourceframe = unrotate_points(centres_r_rotframe, R, centroid) - - real_r = [count_in_disc(centres_r_sourceframe, py, px, probe_radius) - for py, px in real_zones] - ghost_r = [count_in_disc(centres_r_sourceframe, py, px, probe_radius) - for py, px in ghost_zones] - - mean_real_b = np.mean(real_b) if real_b else 0.0 - mean_ghost_b = np.mean(ghost_b) if ghost_b else 0.0 - mean_real_r = np.mean(real_r) if real_r else 0.0 - mean_ghost_r = np.mean(ghost_r) if ghost_r else 0.0 - - print(f"\n=== {label} ===") - angle_deg = np.degrees(np.arctan2(R[0, 1], R[0, 0])) - print(f" PCA principal-axis angle: {angle_deg:+.2f}°") - print(f" baseline: mean(real) = {mean_real_b:.1f}, " - f"mean(ghost) = {mean_ghost_b:.1f}, " - f"ghost/real = {mean_ghost_b / max(mean_real_b, 1e-12):.2f}") - print(f" rotated: mean(real) = {mean_real_r:.1f}, " - f"mean(ghost) = {mean_ghost_r:.1f}, " - f"ghost/real = {mean_ghost_r / max(mean_real_r, 1e-12):.2f}") - - return { - "label": label, - "traced": traced, - "brightness": brightness, - "peaks": peaks, - "real_zones": real_zones, - "ghost_zones": ghost_zones, - "yE_b": yE_b, "xE_b": xE_b, "centres_b": centres_b, - "R": R, "centroid": centroid, - "yE_r": yE_r, "xE_r": xE_r, - "centres_r_rotframe": centres_r_rotframe, - "centres_r_sourceframe": centres_r_sourceframe, - "ratios": (mean_ghost_b / max(mean_real_b, 1e-12), - mean_ghost_r / max(mean_real_r, 1e-12)), - } - - -def plot_overlay(ax, traced, brightness, centres, real_zones, ghost_zones, - grid_lines=None, title=""): - sc = ax.scatter(traced[:, 1], traced[:, 0], c=brightness, cmap="hot", - s=4, alpha=0.45, marker=".") - if grid_lines is not None: - ys_min, ys_max, xs_min, xs_max = grid_lines["bbox"] - for ye in grid_lines["y_edges"]: - ax.plot([xs_min, xs_max], [ye, ye], color="cyan", lw=0.25, alpha=0.5) - for xe in grid_lines["x_edges"]: - ax.plot([xe, xe], [ys_min, ys_max], color="cyan", lw=0.25, alpha=0.5) - ax.scatter(centres[:, 1], centres[:, 0], c="white", s=2, - edgecolors="black", linewidths=0.2) - for (py, px) in real_zones: - ax.plot(px, py, "o", color="lime", markersize=12, - markeredgecolor="black", markeredgewidth=1.5) - for (py, px) in ghost_zones: - ax.plot(px, py, "x", color="red", markersize=12, markeredgewidth=2) - ax.set_xlim(-1.05, 1.05) - ax.set_ylim(-1.05, 1.05) - ax.set_aspect("equal") - ax.set_title(title) - ax.set_xlabel("source x") - ax.set_ylabel("source y") - return sc - - -def main(): - # Three test cases reusing the ghost-peak setup - cases = [ - ("K=2 diagonal", [(-0.5, -0.5), (0.5, 0.5)]), - ("K=2 axis-aligned (control)", [(0.0, -0.5), (0.0, 0.5)]), - ("K=3 offset", [(-0.5, -0.5), (0.0, 0.5), (0.5, -0.2)]), - ] - - fig, axes = plt.subplots(len(cases), 2, figsize=(14, 5.0 * len(cases))) - last_sc = None - - for row, (label, peaks) in zip(axes, cases): - out = run_one(label, peaks) - - # Left: baseline (unrotated) - last_sc = plot_overlay( - row[0], out["traced"], out["brightness"], out["centres_b"], - out["real_zones"], out["ghost_zones"], - grid_lines={ - "y_edges": out["yE_b"], "x_edges": out["xE_b"], - "bbox": (out["yE_b"].min(), out["yE_b"].max(), - out["xE_b"].min(), out["xE_b"].max()), - }, - title=f"{label} — baseline\nghost/real = {out['ratios'][0]:.2f}", - ) - - # Right: PCA-rotated mesh, drawn in source-plane frame. - # The mesh edges are in the rotated frame, so we have to construct the - # full set of cell-corner points there and unrotate them all for the - # axis-aligned plot. For visualisation we just plot the unrotated - # pixel centres + the rotated bounding box, leaving the grid lines - # implicit (drawing rotated lines via plain plot() means many segments). - last_sc = plot_overlay( - row[1], out["traced"], out["brightness"], - out["centres_r_sourceframe"], - out["real_zones"], out["ghost_zones"], - grid_lines=None, - title=f"{label} — PCA-rotated\nghost/real = {out['ratios'][1]:.2f}", - ) - - # Draw the rotated grid as line segments transformed to source frame. - yE_r, xE_r = out["yE_r"], out["xE_r"] - R, centroid = out["R"], out["centroid"] - # Horizontal lines in rotated frame: y = yE_r[i], x in [xE_r.min, xE_r.max] - for ye in yE_r: - pts_rot = np.array([[ye, xE_r.min()], [ye, xE_r.max()]]) - pts_src = unrotate_points(pts_rot, R, centroid) - row[1].plot(pts_src[:, 1], pts_src[:, 0], color="cyan", - lw=0.25, alpha=0.5) - # Vertical lines in rotated frame: x = xE_r[j], y in [yE_r.min, yE_r.max] - for xe in xE_r: - pts_rot = np.array([[yE_r.min(), xe], [yE_r.max(), xe]]) - pts_src = unrotate_points(pts_rot, R, centroid) - row[1].plot(pts_src[:, 1], pts_src[:, 0], color="cyan", - lw=0.25, alpha=0.5) - - fig.colorbar(last_sc, ax=axes, label="brightness", - fraction=0.025, pad=0.02) - out_path = os.path.join(os.path.dirname(__file__), "pca_rotation.png") - fig.savefig(out_path, dpi=130, bbox_inches="tight") - print(f"\nfigure written: {out_path}") - - -if __name__ == "__main__": - sys.exit(main())