diff --git a/llms.txt b/llms.txt index 4f66f35..35935ee 100644 --- a/llms.txt +++ b/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. diff --git a/scripts/demonstrations/README.md b/scripts/demonstrations/README.md new file mode 100644 index 0000000..c49c60a --- /dev/null +++ b/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). diff --git a/scripts/demonstrations/demo_1_calibrate_1d.py b/scripts/demonstrations/demo_1_calibrate_1d.py new file mode 100644 index 0000000..85cf071 --- /dev/null +++ b/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.") diff --git a/scripts/demonstrations/demo_2_calibrate_imaging_ci.py b/scripts/demonstrations/demo_2_calibrate_imaging_ci.py new file mode 100644 index 0000000..3d4063b --- /dev/null +++ b/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.") diff --git a/scripts/demonstrations/demo_3_correct_1d.py b/scripts/demonstrations/demo_3_correct_1d.py new file mode 100644 index 0000000..bec4c57 --- /dev/null +++ b/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.") diff --git a/scripts/demonstrations/demo_4_aggregate_results.py b/scripts/demonstrations/demo_4_aggregate_results.py new file mode 100644 index 0000000..f5509d2 --- /dev/null +++ b/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.")