docs: update README and add benchmarks and output documentation

## What's PyAETNA?

PyAETNA is a flexible, configuration-driven pipeline for fitting

transit, phase‑curve light curves and radial velocities using a Bayesian

framework with high-performance JAX/NumPyro sampling.  It handles

instrument-specific systematics via polynomials, Fourier series, or

Gaussian processes, and produces high-quality diagnostic reports

with minimal user effort.

PyAETNA is a configuration-driven pipeline that fits transit and phase-curve light curves, radial velocities and

transit timing variations in a Bayesian framework, on a CPU or a GPU. It is built on JAX, NumPyro and BlackJAX; one

YAML file describes the whole analysis, and the same model runs unchanged on any device. Instrument systematics are

handled with polynomials, Fourier series or Gaussian processes, and a run ends with diagnostic reports ready to read.

### Features

- **JAX/NumPyro Backend**: JIT-compiled, GPU-accelerated MCMC sampling

- **Three samplers, chosen explicitly** (`production_sampler`):

  - `ns`: GPU-vectorised nested sampling, returning the posterior and ln Z (Bayes factors for model comparison)

  - `nuts`: the No-U-Turn Sampler (gradient-based)

  - `ess`, `aies`: ensemble slice sampling and the affine-invariant ensemble sampler (gradient-free)

- **Peak search before sampling**: optimisations from many prior draws find the peaks of the posterior; NUTS and

  the ensembles start from the best one with its curvature as scales, and several peaks (e.g. period aliases) are

  reported with one configuration per peak to run separately

- **Runs until converged**: NUTS and the ensembles can stop as soon as the slowest parameter has the requested number

  of independent draws (`target_ess`)

- **Named Priors & Global Tying**: tie parameters across different model components (e.g., sharing a stellar rotation period between multiple instrument and corresponding detrending schemes)

- **Comprehensive Light Curve Modeling**:

  - Primary transits, phase curves and secondary eclipses (via `jaxoplanet`)

  - **Transit Timing Variations (TTV)**: per-epoch free transit times with O−C posterior

  - Phase curves with hotspot offsets (Lambert sphere)

  - Ellipsoidal variations with tidal lag

  - Doppler beaming

- **Flexible Systematics Treatment**:

  - Polynomial detrending

  - Error-In-Variables linear detrending

  - Fourier series

  - Gaussian Processes (tinygp with O(N) and O(N³) kernels)

- **Robust Uncertainty Quantification**: 

  - White noise and red noise components

  - Correlated noise modeling

  - Confidence bands from posterior samples

- **Models**: transits, secondary eclipses and phase curves (with hotspot offset, ellipsoidal variation and Doppler

  beaming), the Rossiter–McLaughlin effect, Keplerian radial velocities, and free transit times (TTVs)

- **Systematics and noise**: polynomials, Fourier series, errors-in-variables detrending, Gaussian processes

  (including multi-dimensional ones shared between light curves and activity indicators), white and correlated noise

- **Joint fits**: several instruments and data types at once, with shared orbital parameters and independent

  systematics; parameters can be shared across components through named priors

- **Bayesian inference**: nested sampling, which returns the posterior and the Bayesian evidence for model

  comparison, and the NUTS and ensemble samplers for parameter estimation; you choose the sampler, and a peak search

  before sampling checks that the posterior has a single peak (or hands you one configuration per peak)

- **Reports and derived physics**: PDF reports with fits, residuals, corner and trace plots, periodograms and goodness

  of fit, and the physical quantities derived from the posterior (stellar and planetary properties, energy balance,

  albedos)

- **Reproducible and restartable**: checkpointed runs, portable result files, full provenance of every run

## Why Use PyAETNA?

- **Single YAML, complete analysis.** Define data paths, priors,

  systematics, MCMC strategy, and post-processing physics in one

  human-readable configuration file — no boilerplate Python scripts required.

- **GPU-accelerated end-to-end.** The entire forward model, likelihood, and

  gradient computation runs inside JIT-compiled JAX, delivering order-of-

  magnitude speedups on GPU/TPU hardware without code changes.

- **Transparent sampling.** You choose the sampler; a peak search before it finds where the posterior is, so NUTS

  and the ensembles start there without a warmup, and a posterior with several peaks is never sampled silently as

  if it had one.

- **Checkpointed, restartable runs.** The sampler state is saved

  periodically, so that long runs can survive crashes, and subsequent

  launches resume exactly where they left off.

- **Portable results.** Samples, diagnostics, the MAP, the evidence and the

  model components are written as plain `.npz` and `.json` files that open

  with NumPy alone and stay readable as PyAETNA changes.

- **Mixed-precision ephemeris firewall.** Transit mid-times are folded in

  float64 while the rest of the model evaluates in float32 for speed,

  eliminating the catastrophic cancellation that plagues naïve

  single-precision transit fits.

- **Rich systematics toolkit.** Polynomial trends, Fourier series, and

  scalable Gaussian Processes (O(N) quasiseparable or O(N³) dense kernels

  via tinygp) can be combined freely, scoped per-visit or per-instrument,

  and stacked on the same dataset.

- **Joint multi-instrument, multi-type fitting.** Photometric light curves

  and radial velocities from different instruments are modelled

  simultaneously with fully shared orbital parameters and independent

  per-instrument systematics.

- **High-quality and informative plots.** Diagnostic PDFs include phase-folded fits

  with residuals, best-fit models per visit, corner plots, trace plots,

  Lomb–Scargle periodograms, and an optional physics report with Cowan &

  Agol energy-balance, brightness temperatures, and passband-integrated

  albedos to clearly check the outcome of model exploration.

- **Physics post-processing in one step.** A dedicated `physics:` stage

  derives limb-darkening coefficients, orbital geometry, equilibrium

  temperatures, Bond albedos, and heat redistribution efficiencies —

  propagating full MCMC uncertainties via Monte Carlo resampling.

- **Model-atmosphere support.** BT-Settl and PHOENIX synthetic spectra are

  downloaded on the fly (with local caching) for passband-integrated

  brightness-temperature inversion, replacing the default blackbody

  approximation with no additional libraries to install.

- **One file, one command.** The data, priors, systematics, sampler and post-processing are set in a single readable

  configuration file, which an [online configurator](http://pyaetna-e92eeb.pages.ict.inaf.it/) can write for you; no

  scripting is needed.

- **Fast, on a laptop or a server.** The forward model and its derivatives are compiled with JAX and evaluated

  on many points at once, so a GPU gives up to 7 times the speed of the CPU of the same machine. On the benchmark

  problems PyAETNA was 1.6 to 45 times faster than allesfitter on the CPU, and faster than pyaneti with the same

  MCMC.

- **Trustworthy evidence.** Its nested sampler biases the Bayesian evidence 6 to 10 times less than dynesty on two of three test problems

  of known evidence, and the bias cancels in the Bayes factor between models.

- **Honest about its limits.** Nested sampling in PyAETNA needs 3 to 5 times more likelihood evaluations than

  juliet or PyORBIT (about 4 times slower on a CPU, close on a GPU), and keeps only one of several separate peaks;

  the peak search warns of the latter. The comparisons, with the cases where other codes are better, are in

  **[How PyAETNA compares with other codes](docs/benchmarks.md)**.

- **Agrees with the established codes.** On shared data sets, posteriors match those of juliet, allesfitter, pyaneti

  and PyORBIT within 0.7σ for almost all parameters.

- **Transparent and robust.** Numerical precision is handled automatically where it matters (transit times), runs

  resume after a crash, and results are plain `.npz` and `.json` files that open with NumPy alone.

## Quick Start

## Documentation

Comprehensive documentation is available:

PyAETNA is presented in the software paper Scandariato et al. (in preparation). Additional documentation is available at the following pages:

- **[Configuration Guide](docs/config/index.md)** - Complete reference for all YAML configuration fields

  - MCMC sampling strategies (NUTS, ensemble samplers AIES/ESS, multimodal detection)

  - Prior distributions (uniform, normal, truncnorm, lognormal, beta)

  - Detrending options (polynomials, Fourier, Gaussian Processes)

  - Best practices and troubleshooting

  A complete **online configurator** is available at http://pyaetna-e92eeb.pages.ict.inaf.it/. The INAF Pages server does not offer HTTPS yet; if your browser or network refuses plain HTTP, open `docs/config-generator.html` from a copy of the repository (cloned, or downloaded over HTTPS from the repository page): it is a single static page and works offline.

- **[How PyAETNA compares with other codes](docs/benchmarks.md)** - Speed, evidence reliability and agreement with juliet, allesfitter, pyaneti and PyORBIT, with the cases where other codes are better

- **[Installation & Environment Guide](docs/dependencies.md)** - Everything about setting up PyAETNA

  - Installation instructions (CPU and GPU)

  - Version compatibility and troubleshooting

  - Conda and Docker alternatives

- **Output & Diagnostics.** After a successful run, PyAETNA creates an output directory named `PyAETNAouts/` inside the specified `data.working_dir`. Every file is named `<config>_<kind>`, where `<config>` is the configuration file name without `.yaml`, so several configurations fitted to the same data keep separate outputs. The folder contains:

- **[Configuration Guide](docs/config/index.md)** - Complete reference for all YAML configuration fields

  - Samplers (nested sampling, NUTS, ensemble samplers), the peak search and convergence targets

  - Prior distributions (uniform, normal, truncnorm, lognormal, beta)

  - Detrending options (polynomials, Fourier, Gaussian Processes)

  - Best practices and troubleshooting

  - **Log (`<config>_run.log`)**: the record of the run: provenance and device, data and priors, the sampler's progress, the posterior summary with R-hat and effective sample sizes (or the Bayesian evidence), goodness of fit, derived quantities, warnings and the BibTeX entries to cite. **[Reading the run log](docs/outputs/run_log.md)** goes through it block by block and lists what to check first.

  - **Results (`<config>_results.json`)**: the run settings (sampler, device, precision, chains, draws, seed), the provenance (code revision, package versions), R-hat and ESS for every parameter, the fixed parameters, the MAP, the goodness of fit (χ², maximum likelihood, AIC, AICc, BIC), the nested-sampling evidence, the instruments and the configuration.

  - **Peaks (`<config>_peaks.json`)**: the peaks of the posterior found by the peak search before sampling, best first: energy, maximum log-likelihood, approximate (Laplace) evidence and weight, whether the peak is significant, position and widths of the parameters, and its Gaussian approximation (`z`, `cov`), from which a later run can start (`sampling.init_from_peaks`). With several significant peaks and `sampling.on_multiple_peaks: stop`, the run stops here and writes one configuration per peak, `<config>_peak<k>.yaml`, next to the original ([details](docs/config/sampling.md#the-peak-search-all-samplers)).

  - **Posterior samples (`<config>_samples.npz`)**: one array per parameter, shaped `(chain, draw)`, and the `potential_energy` of every draw when the sampler records it (NUTS, nested sampling).

  - **Model components (`<config>_components.npz`)**: the data (`bjd`, `y`, `dy`, `instrument`, `visit`) and, on every data point, the MAP model and its components (`map/total`, `map/planet`, `map/gp_mean`, …), the posterior median (`median/…`) and the 16th–84th percentile band (`q16_84/…`).

  - **Report (`<config>_report.pdf`)** and, with a `physics:` section, **physics report (`<config>_physics.pdf`)**: phase-folded fits, residuals, corner plots and derived quantities. **[Reading the PDF reports](docs/outputs/report.md)** shows an example of every page and how to read it.

  - **Checkpoint (`<config>_checkpoint.pckl`)**: the sampler state, used to resume a run and by the reports. It holds PyAETNA objects that a later version may not be able to open: keep the `.json` and `.npz` files as the record of a run.

- **[Online configurator](http://pyaetna-e92eeb.pages.ict.inaf.it/)** - Write or edit a configuration file in the browser

  - The INAF Pages server does not offer HTTPS yet; if your browser or network refuses plain HTTP, open `docs/config-generator.html` from a copy of the repository (cloned, or downloaded over HTTPS from the repository page): it is a single static page and works offline

  The `.json` and `.npz` files hold only numbers, strings and arrays, so they can be read with `json` and `numpy` (or any language) without PyAETNA; `PyAETNA.outputs` has loaders for them. To load, plot and rebuild the model from them, see the [Loading and plotting results](docs/load_binary_model.ipynb) notebook.

- **[Output & Diagnostics](docs/outputs/index.md)** - What a run writes and how to read it

  - The run log, the results and sample files, the peaks of the posterior and the PDF reports

  - Guides to [reading the run log](docs/outputs/run_log.md) and [the PDF reports](docs/outputs/report.md)

  - Loading the results without PyAETNA

Additional resources:

- Test files in `tests/`

If you use PyAETNA in your research, please add the following citation in your publication:

#TODO add citation to PyAETNA

#TODO add bibtex citation to PyAETNA

Please also cite the relevant works

listed in [docs/CITATIONS.md](docs/CITATIONS.md). For your convenience, at the end of every run the

[← Back to README](../README.md) · [Configuration reference](config/index.md) · [Installation](dependencies.md)

# How PyAETNA compares with other codes

This page summarises the benchmarks of the PyAETNA software paper: the same data, priors and sampler settings fitted

with PyAETNA and with other widely used codes, on the same machine, each code at its best number of CPU cores.

The scripts, data and results are in [`paper_benchmarks/`](../paper_benchmarks/README.md); every number below comes from

the `results/*.json` files there, and the paper tables are generated from them. The benchmarks were run with the version

circulated to the co-authors on 2026-10-05; PyAETNA's samplers have changed since (see the

[CHANGELOG](../CHANGELOG.md)) and the affected benchmarks are being repeated, so figures may be updated.

The comparison codes are juliet (dynesty), allesfitter (dynesty and emcee), pyaneti (its ensemble MCMC) and PyORBIT

(dynesty, with the ARoME Rossiter–McLaughlin model).

## Do the codes agree?

Yes. On the shared data sets (K2-140, Pi Men, WASP-18, TOI-216, K2-100, TOI-5398), the posteriors of PyAETNA agree

with those of the other codes within 0.7σ for almost all parameters. The larger differences are traced to the other

codes' choices (the coarse model grid of allesfitter, its stellar file for TOI-216, error weights, priors) and are

discussed in the paper. PyAETNA's closed-form Rossiter–McLaughlin anomaly agrees with ARoME to 0.27 m/s.

## Where PyAETNA is faster

Wall time on the same laptop (24 logical CPUs; the GPU is an RTX 5070 Ti), including start-up and compilation:

| Data set | Reference code | Reference | PyAETNA, CPU | PyAETNA, GPU |

|---|---|---|---|---|

| WASP-18 (phase curve) | allesfitter, MCMC | 43 min | 64 s (same MCMC), 37 s (nested sampling) | 32 s, 10 s |

| Pi Men, circular orbit | allesfitter, nested sampling | 1.0 h | 83 s | 32 s |

| Pi Men, eccentric orbit | allesfitter, nested sampling | 1.2 h | 121 s | 42 s |

| TOI-216 (transit times) | allesfitter, nested sampling | 1.5 h | 54 min | 32 min |

| K2-100 (multi-GP) | pyaneti, same MCMC | 175 s | 147 s | 100 s |

- Against **allesfitter**, PyAETNA is 1.6 to 45 times faster on the CPU, with a similar number of likelihood

  evaluations for the same sampler (WASP-18: 2.5 million each): the difference is the cost of one evaluation, which

  PyAETNA compiles with JAX.

- Against **pyaneti**, with the same ensemble MCMC (same walkers, steps and start), PyAETNA is faster on the CPU.

- A **GPU** makes PyAETNA up to 7 times faster than its own CPU run on the same machine, and never slower on these

  problems (K2-100 nested sampling is the closest: 422 s against 426 s). Compared with the reference codes on their

  CPUs, the GPU runs are 1.8 to 120 times faster, except against PyORBIT on TOI-5398.

## Where PyAETNA is slower

- **Nested sampling against juliet and PyORBIT** (K2-140: 475 s against 119 s on the CPU; TOI-5398: 52 s against

  14 s): PyAETNA needs 3 to 5 times more likelihood evaluations, because its slice sampler spends more evaluations per

  new point than dynesty's random walk. On the GPU the gap closes for K2-140 (64 s). The next section explains what the

  extra evaluations buy.

- **TOI-216 on the CPU** (54 min against allesfitter's 1.5 h) is a smaller gain than for the other data sets, because

  PyAETNA's nested sampler uses 19 million likelihood evaluations against allesfitter's 2.5 million; the GPU result is

  32 min.

- The **nautilus** sampler needs 5 to 37 times fewer evaluations than PyAETNA's nested sampler for the same evidence,

  but is slower in wall time on these models, because its neural networks train on the CPU. It is faster only for

  the K2-100 multi-GP likelihood (about 50 times fewer evaluations).

The complete list, with causes and planned improvements, is kept in the repository

(`localDocs/performance/PYAETNA_SLOWER_CASES.md`).

## Reliability of the Bayesian evidence

Nested sampling returns the evidence $Z$ that model comparison needs, but a sampler that takes too few steps per

new point reports an evidence that is biased high. On problems whose evidence is known exactly (24 runs of PyAETNA,

8 of dynesty, 500 live points), the bias of $\ln Z$ in nats is:

| Problem | PyAETNA | dynesty, slice | dynesty, random walk |

|---|---|---|---|

| Correlated 20-D Gaussian | +0.09 ± 0.06 | +0.70 ± 0.10 | +0.91 ± 0.15 |

| Gaussian cut by the prior bounds | +0.22 ± 0.11 | +1.54 ± 0.11 | +1.25 ± 0.17 |

| Curved posterior ($c = 0.3$) | +3.89 ± 0.07 | +5.66 ± 0.26 | +3.16 ± 1.82 |

On the first two problems, PyAETNA's bias is 6 to 10 times smaller than dynesty's, and in the Bayes factor between

two nested models, where the biases largely cancel, it stays below 0.25 nats. On the K2-140 problem, dynesty's random

walk makes the scatter of $\ln Z$ five times larger than its own reported error, while PyAETNA's scatter is close

to its error (1.0 to 1.5 times). This is the price of PyAETNA's extra evaluations.

**Limitation.** On a posterior with two well-separated peaks, PyAETNA's nested sampler keeps only one of them in every

run, while dynesty keeps both. This is why PyAETNA looks for the peaks before sampling and stops, with one

configuration per peak, if it finds several (see [Sampling](config/sampling.md#the-peak-search-all-samplers)).

## Which sampler to choose

On the benchmark problems, nested sampling was 4 to 600 times more efficient than NUTS in effective samples per

second, including a synthetic system with up to 100 free transit times (106 parameters), where NUTS had not converged

after 10 000 draws per chain. Ensemble samplers are efficient on few parameters and degrade above 20 to 30. The

recommendations, and how each sampler reaches convergence with and without the peak search, are in

[Sampling](config/sampling.md#choosing-a-sampler); the latter comparison is being finalised.

## Reproducing the results

Each benchmark folder contains the data, configurations and scripts, and `./rerun_all.sh` repeats the analysis. See the

[reproducibility policy](../paper_benchmarks/README.md).

[← Back to README](../../README.md) · [Reading the run log](run_log.md) · [Reading the PDF reports](report.md)

# Output & Diagnostics

After a successful run, PyAETNA creates an output directory named `PyAETNAouts/` inside the specified

`data.working_dir`. Every file is named `<config>_<kind>`, where `<config>` is the configuration file name without

`.yaml`, so several configurations fitted to the same data keep separate outputs. The folder contains:

- **Log (`<config>_run.log`)**: the record of the run: provenance and device, data and priors, the sampler's progress,

  the posterior summary with R-hat and effective sample sizes (or the Bayesian evidence), goodness of fit, derived

  quantities, warnings and the BibTeX entries to cite. **[Reading the run log](run_log.md)** goes through it block by

  block and lists what to check first.

- **Results (`<config>_results.json`)**: the run settings (sampler, device, precision, chains, draws, seed), the

  provenance (code revision, package versions), R-hat and ESS for every parameter, the fixed parameters, the MAP, the

  goodness of fit (χ², maximum likelihood, AIC, AICc, BIC), the nested-sampling evidence, the instruments and the

  configuration.

- **Peaks (`<config>_peaks.json`)**: the peaks of the posterior found by the peak search before sampling, best first:

  energy, maximum log-likelihood, approximate (Laplace) evidence and weight, whether the peak is significant, position

  and widths of the parameters, and its Gaussian approximation (`z`, `cov`), from which a later run can start

  (`sampling.init_from_peaks`). With several significant peaks and `sampling.on_multiple_peaks: stop`, the run stops

  here and writes one configuration per peak, `<config>_peak<k>.yaml`, next to the original

  ([details](../config/sampling.md#the-peak-search-all-samplers)).

- **Posterior samples (`<config>_samples.npz`)**: one array per parameter, shaped `(chain, draw)`, and the

  `potential_energy` of every draw when the sampler records it (NUTS, nested sampling).

- **Model components (`<config>_components.npz`)**: the data (`bjd`, `y`, `dy`, `instrument`, `visit`) and, on every

  data point, the MAP model and its components (`map/total`, `map/planet`, `map/gp_mean`, …), the posterior median

  (`median/…`) and the 16th–84th percentile band (`q16_84/…`).

- **Report (`<config>_report.pdf`)** and, with a `physics:` section, **physics report (`<config>_physics.pdf`)**:

  phase-folded fits, residuals, corner plots and derived quantities. **[Reading the PDF reports](report.md)** shows an

  example of every page and how to read it.

- **Checkpoint (`<config>_checkpoint.pckl`)**: the sampler state, used to resume a run and by the reports. It holds

  PyAETNA objects that a later version may not be able to open: keep the `.json` and `.npz` files as the record of a

  run.

The `.json` and `.npz` files hold only numbers, strings and arrays, so they can be read with `json` and `numpy` (or any

language) without PyAETNA; `PyAETNA.outputs` has loaders for them. To load, plot and rebuild the model from them, see

the [Loading and plotting results](../load_binary_model.ipynb) notebook.

---

**Last Updated:** 2026-10-06  

**PyAETNA Version:** 4.0.0 (unreleased)