feat: peak-search starts spread over the CPU host devices; ESS walkers stay...

      <input type="number" id="m-wpfp" data-key="walkers_per_free_parameter" value="3" min="2" step="1">

      <div class="help">Coefficient multiplied by n_free_params to determine the number of ensemble walkers. Must be &ge; 2. The actual walker count is the first even number &ge; coefficient &times; n_free_params (at least 3&nbsp;&times;&nbsp;n_free_params when ESS uses <code>KDEMove</code>). Typical values: 3&ndash;5.</div>

    </div>

    <div class="fg">

      <label><span class="key-name">num_walkers</span></label>

      <input type="number" id="m-nw" data-key="num_walkers" data-emit="set" min="4" step="2" placeholder="(from walkers_per_free_parameter)">

      <div class="help">Optional: the exact number of walkers (even), instead of <code>walkers_per_free_parameter</code>, e.g. to match the ensemble of another code such as emcee.</div>

    </div>

    <div class="fg">

      <label><span class="key-name">num_samples_if_ensemble</span></label>

      <input type="number" id="m-nse" data-key="num_samples_if_ensemble" value="5000" min="1">

| Run length | Its own stopping rule (`dlogz_if_ns`) | `num_samples_if_nuts`, or until converged (`target_ess`) | `num_samples_if_ensemble`, or until converged (`target_ess`) |

**Rule of thumb.** If you need Bayesian evidence to compare models (circular

vs eccentric, one vs two planets, with or without a GP), use **`ns`**. It was

also the more efficient route against NUTS in every benchmark of the software paper. For parameter estimation of a

model with few parameters, an ensemble started at the peak is often the fastest. NUTS is the gradient-based

alternative. Whatever the sampler, a peak search runs first (below): it finds separate peaks, such as period aliases,

which no sampler of PyAETNA weights correctly within one run.

vs eccentric, one vs two planets, with or without a GP), use **`ns`**. It was also the fastest route to a converged

posterior in every benchmark of the software paper ([comparison](../benchmarks.md#pyaetnas-samplers-compared)), and it

does not depend on the peak search. **NUTS**, started from the peak, converged on every benchmark, up to 106

parameters, in 1 to 7 times the time of nested sampling: a good independent check. The **ensembles** converged only up

to about 20 parameters, and never from the prior: keep them for reproducing an `emcee`-like analysis. Whatever the

sampler, a peak search runs first (below): it finds separate peaks, such as period aliases, which no sampler of PyAETNA

weights correctly within one run. It is a probe, not a guarantee: with few starts it can miss a peak, so repeat it (or

raise `peak_search_starts`) when a posterior is known to be complex.

![The sampling branches](img/branches.png)

| Sampler | Pros | Cons | Use it when |

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

| **`ns`** | $\ln Z$ for model comparison; independent posterior draws, no burn-in to judge; starts from the prior, so it cannot get stuck near a bad starting point; clear stopping rule (`dlogz_if_ns`); vectorised on GPU, parallel on CPU cores; checkpointed | Cost grows with the information gained over the prior (so very wide priors cost more) and roughly with the square of the number of free parameters; needs proper priors (all PyAETNA priors are); of separate peaks all but one are lost, and $\ln Z$ is overestimated by several nats on strongly curved posteriors | You need Bayes factors; recommended default |

| **`nuts`** | Gradient-based moves; starts at the peak with the parameter scales of its curvature; can run until converged; checkpointed | No $\ln Z$; cannot move between separate peaks; long trajectories on narrow, correlated or curved posteriors | Parameter estimation, many parameters; a check of a nested-sampling result |

| **`ess`** / **`aies`** | Gradient-free; robust to rough likelihood surfaces; very efficient from the peak at low dimension; can run until converged; checkpointed | No $\ln Z$; walkers cannot move between separate peaks; mixing degrades above ~20–30 parameters; the burn-in is discarded in the reports | Parameter estimation with few parameters; reproducing an `emcee`-like analysis |

| **`nuts`** | Gradient-based moves; starts at the peak with the parameter scales of its curvature; converged from the peak on every benchmark, up to 106 parameters; can run until converged; checkpointed | No $\ln Z$; cannot move between separate peaks; needs the peak search (from the prior it failed on three of five benchmarks); long trajectories on narrow, correlated or curved posteriors | Parameter estimation, many parameters; a check of a nested-sampling result |

| **`ess`** / **`aies`** | Gradient-free; robust to rough likelihood surfaces; reproduce `emcee`-like analyses; can run until converged; checkpointed | No $\ln Z$; walkers cannot move between separate peaks; slow from about 20 parameters on, even from the peak, and never converged from the prior in the paper's benchmarks; the burn-in is discarded in the reports | Reproducing an `emcee`-like analysis; parameter estimation with few parameters |

> 💡 **Terminology.** PyAETNA's `ess` is **Ensemble Slice Sampling**

> (Karamanis & Beutler 2020; Karamanis, Beutler & Peacock 2021), *not*

posterior for each peak. The log lists the peaks, and `<config>_peaks.json` saves them (with the Gaussian

approximations). The log also gives the time of one likelihood and one gradient evaluation.

**Cost.** On the paper's benchmarks the search took 16 s to about 3 min on an H100 GPU, set more by the cost of the

model (a GP over the transits is the most expensive) than by the number of parameters. On CPUs the starts are spread

over the host devices that PyAETNA makes of the cores (32 of them for 32 starts on the 60-core server), and the search

took 28 to 95 s there, 2 to 3 times the GPU time. It

is worth its cost before NUTS and the ensembles (2–25 % of a NUTS run; without it NUTS failed on three of the five

benchmarks); before nested sampling it only tells you, before any sampling, whether the posterior has separate peaks.

NUTS chains and ensemble walkers then start from draws of the Gaussian approximation of the best peak, so they need

no warmup: NUTS uses the curvature as its first mass matrix, and the ensemble starts already spread like the

posterior.

## Ensemble samplers (`production_sampler: ess` or `aies`)

- **Walkers**: `walkers_per_free_parameter` × $n_{\text{free\_params}}$, the first even number, at least 4 (at least

  $3\times n_{\text{free\_params}}$ when ESS uses `KDEMove`). They start from draws of the Gaussian approximation of

  the best peak.

  $3\times n_{\text{free\_params}}$ when ESS uses `KDEMove`), or exactly `num_walkers` (an even number), e.g. to

  match another code's ensemble. They start from draws of the Gaussian approximation of the best peak.

- **Steps**: `num_samples_if_ensemble` per walker, or fewer with `target_ess`. The first steps are a burn-in, short

  when the walkers start at the peak: discard it in the reports with `reports.discard_extra_warmup` (with

  `target_ess`, the log gives the burn-in, 2 autocorrelation times).

| `sampling.dense_mass` | no | `false` | bool | Dense mass matrix |

| `sampling.adapt_mass_matrix_if_nuts` | no | `true` | bool | Adapt the mass matrix during the adaptation (it starts from the curvature of the peak) |

| **Ensemble samplers** | | | | |

| `sampling.walkers_per_free_parameter` | when `ess`, `aies` | — | float ($\ge 2$) | Walkers per free parameter (rounded up to an even number, at least 4; at least $3\times n_{\text{params}}$ with `KDEMove`) |

| `sampling.num_walkers` | no | — | int (even, $\ge 4$) | Exact number of walkers, instead of `walkers_per_free_parameter` |

| `sampling.walkers_per_free_parameter` | when `ess`, `aies` (unless `num_walkers`) | — | float ($\ge 2$) | Walkers per free parameter (rounded up to an even number, at least 4; at least $3\times n_{\text{params}}$ with `KDEMove`) |

| `sampling.num_samples_if_ensemble` | when `ess`, `aies` | — | int ($> 0$) | Steps per walker (the maximum with `target_ess`) |

| `sampling.ess_moves` | no | `{DifferentialMove: 0.7, KDEMove: 0.3}` | dict | ESS move weights |

| `sampling.aies_moves` | no | `{DEMove: 0.5, StretchMove: 0.5}` | dict | AIES move weights |

data:

  file_path: test_spectral_mixture.dat

  working_dir: /dati/software/pyLib/PyAETNA/examples

  working_dir: ./

data_type:

- LC

detrending:

            "target_ess",

            "use_float64",

            "walkers_per_free_parameter",

            "num_walkers",

        }

    ),

    "star": frozenset(

    }

    if extra:

        out.update(extra)

    atomic_write(save_path, lambda f: dill.dump(out, f, protocol=4))

    # protocol 5 writes contiguous arrays without an in-memory copy (a 9 GB transit-time array of a long run)

    atomic_write(save_path, lambda f: dill.dump(out, f, protocol=5))

# Top-level sections that may change between a run and its extension.