Commit af86d9a5 authored by Gaetano Scandariato's avatar Gaetano Scandariato
Browse files

docs: benchmark page and README from the sampler comparison with the CPU peak search

parent 342cf497
Loading
Loading
Loading
Loading
+9 −10
Changes for README.md: 9 added lines, 10 removed lines.
Original line number Diff line number Diff line
@@ -48,17 +48,16 @@ handled with polynomials, Fourier series or Gaussian processes, and a run ends w
  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
  on many points at once, so the same configuration runs on a CPU or, faster, on a GPU. On the benchmark problems
  PyAETNA was much faster than allesfitter and faster than pyaneti with the same MCMC.
- **Trustworthy evidence.** On problems of known evidence, its nested sampler biases the Bayesian evidence several
  times less than dynesty, and the bias cancels in the Bayes factor between models.
- **Honest about its limits.** Its nested sampler needs more likelihood evaluations than juliet's or PyORBIT's, so it
  is slower than they are on a CPU, and it keeps only one of several separate peaks, which the peak search before
  sampling warns of. The numbers, 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.
- **Agrees with the established codes.** On shared data sets, its posteriors match those of juliet, allesfitter,
  pyaneti and PyORBIT.
- **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.

+100 −68
Changes for docs/benchmarks.md: 100 added lines, 68 removed lines.
Original line number Diff line number Diff line
[← Back to README](../README.md) · [Configuration reference](config/index.md) · [Installation](dependencies.md)

<!-- GENERATED by paper_benchmarks/make_docs.py from paper_benchmarks/*/results/*.json -- edit that script, not this page. -->

# 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.
This page summarises the benchmarks of the PyAETNA software paper: the same data and priors fitted with PyAETNA and
with other widely used codes, on the same machines, each code at its best number of CPU cores, and PyAETNA set up as
close as possible to the other code (the same sampler where the other code runs an ensemble MCMC, with the same
walkers, steps and starting points; no peak search before sampling, since none of the other codes optimises first).
The scripts, data and results are in [`paper_benchmarks/`](../paper_benchmarks/README.md), and every number below is
computed from its `results/*.json` files.

The comparison codes are juliet (dynesty), allesfitter (dynesty and emcee), pyaneti (its ensemble MCMC) and PyORBIT
(dynesty, with the ARoME Rossiter–McLaughlin model).
(dynesty, with the ARoME Rossiter–McLaughlin model). The machines are a laptop (24 logical CPUs, an RTX 5070 Ti GPU)
and a server (60 logical CPUs, an H100 GPU).

## 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
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 a fraction of their widths for almost all parameters. The few larger differences are
traced to the other codes' choices (model grids, stellar inputs, error weights, priors) and are discussed in the
paper, parameter by parameter.

## Run times

Wall times, including start-up and, for PyAETNA, compilation. In parentheses, the number of logical CPUs of the run:
the best of a scan for each code and problem. "NS": nested sampling; "MCMC": the reference code's ensemble sampler,
reproduced by PyAETNA. Last column: likelihood evaluations of the laptop CPU run.

| Data set | Code | Sampler | Laptop CPU | Laptop GPU | Server CPU | Server GPU | Evaluations |
|---|---|---|---|---|---|---|---|
| K2-140 | juliet | NS (dynesty, 500 live points) | 119 s (20) | – | 63 s (60) | – | 1.06 M |
| K2-140 | PyAETNA | NS (slice, 500 live points) | 8 min (24) | 56 s | 4 min (60) | 41 s | 5.71 M |
| Pi Men (circ) | allesfitter | NS (dynesty, 500 live points) | 62 min (24) | – | 21 min (60) | – | 3.16 M |
| Pi Men (circ) | PyAETNA | NS (slice, 500 live points) | 75 s (24) | 27 s | 67 s (60) | 22 s | 5.15 M |
| Pi Men (ecc) | allesfitter | NS (dynesty, 500 live points) | 73 min (24) | – | 23 min (60) | – | 3.66 M |
| Pi Men (ecc) | PyAETNA | NS (slice, 500 live points) | 117 s (24) | 32 s | 89 s (60) | 26 s | 6.03 M |
| WASP-18 | allesfitter | MCMC (emcee DEMove, 500 walkers, 5000 steps) | 43 min (24) | – | 18 min (60) | – | 2.50 M |
| WASP-18 | PyAETNA | MCMC (emcee DEMove, 500 walkers, 5000 steps) | 60 s (24) | 23 s | 45 s (60) | 27 s | 2.50 M |
| WASP-18 | PyAETNA | NS (slice, 500 live points) | 36 s (24) | 9 s | 32 s (60) | 10 s | 1.78 M |
| TOI-216 | allesfitter | NS (dynesty dynamic, 500 live points) | 88 min (24) | – | 34 min (60) | – | 2.48 M |
| TOI-216 | PyAETNA | NS (slice, 500 live points) | 53 min (12) | 29 min | 21 min (60) | 9 min | 19.11 M |
| K2-100 | pyaneti | MCMC (stretch move, 100 chains, 10000 steps) | 3 min (12) | – | 3 min (32) | – | 1.00 M |
| K2-100 | PyAETNA | MCMC (stretch move, 100 chains, 10000 steps) | 2 min (24) | 99 s | 89 s (60) | 58 s | 1.00 M |
| K2-100 | PyAETNA | NS (slice, 500 live points) | 7 min (24) | 7 min | 3 min (60) | 89 s | 4.43 M |
| TOI-5398 | PyORBIT | NS (dynesty, 500 live points) | 14 s (24) | – | 9 s (48) | – | 0.33 M |
| TOI-5398 | PyAETNA | NS (slice, 500 live points) | 49 s (24) | 13 s | 33 s (60) | 13 s | 0.92 M |

**Where PyAETNA is faster.**

- Against **allesfitter**, PyAETNA is 1.6 to 50 times faster on the same CPUs. With the same
  sampler (WASP-18, emcee's moves) both codes evaluate the likelihood the same number of times: the difference is the
  cost of one evaluation, which PyAETNA compiles with JAX.
- Against **pyaneti**, with the same ensemble MCMC, PyAETNA is 1.2 to 1.8 times faster.
- A **GPU** makes PyAETNA up to 8.5 times faster than the CPU of the same machine, and it is never
  slower (the smallest ratio, 1.02, is for K2-100 (NS, laptop)). Against the reference codes on their
  CPUs, the GPU runs are 1.1 to 137 times faster, except for TOI-5398 against PyORBIT.

**Where PyAETNA is slower.**

- **Nested sampling against juliet and PyORBIT**: PyAETNA is 3.5 to 4.0 times slower than juliet and
  3.4 to 3.8 times slower than PyORBIT on the same CPUs, because its slice sampler needs
  2.7 to 5.4 times more likelihood evaluations than dynesty's random walk. The next section shows 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`).
The complete list of the cases where PyAETNA is slower, with their causes and planned improvements, is kept in the
repository.

## 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:
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.
| 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 5.8 to 11 times smaller than dynesty's, and in the Bayes factor
between two nested models, where the biases largely cancel, it stays below 0.21 nats.

**Limitation.** On a posterior with two well-separated peaks, PyAETNA's nested sampler keeps only one of them, while
dynesty keeps both. This is why PyAETNA looks for the peaks before sampling and, by default, stops with one
configuration per peak if it finds several (see [Sampling](config/sampling.md#the-peak-search-all-samplers)).

## PyAETNA's samplers compared

Each of PyAETNA's samplers was run on five benchmarks until the slowest parameter had 1000
independent draws and the chains agreed (R-hat ≤ 1.01), with and without the peak search before sampling. Wall time to
that point; "not converged" when the run reached its maximum length or its time limit first.

| Data set | Machine | Peak search | Nested sampling | NUTS | ESS | AIES |
|---|---|---|---|---|---|---|
| WASP-18 | server GPU | yes | 26 s | 3 min | 107 s | 110 s |
| WASP-18 | server GPU | no | 22 s | 2 min | not converged (35 min) | not converged (16 min) |
| WASP-18 | server CPU | yes | 72 s | 4 min | 7 min | 3 min |
| WASP-18 | server CPU | no | 32 s | 4 min | not converged (120 min) | not converged (13 min) |
| K2-140 | server GPU | yes | 67 s | 7 min | 24 min | not converged (27 min) |
| K2-140 | server GPU | no | 29 s | not converged (2.8 h) | not converged (71 min) | not converged (27 min) |
| K2-140 | server CPU | yes | 5 min | 27 min | not converged (119 min) | not converged (48 min) |
| K2-140 | server CPU | no | 4 min | not converged (107 min) | not converged (119 min) | not converged (44 min) |
| K2-100 | server GPU | yes | 119 s | 2 min | 23 min | 9 min (on a secondary peak) |
| K2-100 | server GPU | no | 90 s | not converged (3.0 h) | not converged (2.7 h) | not converged (28 min) |
| K2-100 | server CPU | yes | 5 min | 7 min | 20 min | 10 min |
| K2-100 | server CPU | no | 3 min | not converged | not converged (120 min) | not converged (29 min) |
| TOI-216 | server GPU | yes | 12 min | 43 min | not converged (2.1 h) | not converged (28 min) |
| TOI-216 | server GPU | no | 9 min | not converged (3.0 h) | not converged (3.0 h) | not converged (25 min) |
| TTV system, 106 parameters | server GPU | yes | 12 min | 25 min | not converged (95 min) | not converged (34 min) |
| TTV system, 106 parameters | server GPU | no | 12 min | 26 min | not converged (2.3 h) | not converged (35 min) |

Where they converged, NUTS and the ensembles agree with nested sampling. One AIES run on K2-100 converged on a
secondary peak: its peak search missed the main one, and a convergence test, which compares chains with each other,
cannot see a peak that no chain has visited. Nested sampling, which starts from the prior, does not depend on the
search. The choice of sampler is discussed in [Sampling](config/sampling.md#choosing-a-sampler).

## 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).
Each benchmark folder contains the data, configurations and scripts, and `./rerun_all.sh` repeats the analysis. See
the [reproducibility policy](../paper_benchmarks/README.md).