Deborah.Miriam.Reweighting

Deborah.Miriam.Reweighting.ReweightingSolver — Type

mutable struct ReweightingSolver{T}

Mutable solver object for computing reweighting factors over a collection of lattice-QCD ensembles.

This struct encapsulates all state required to iteratively solve for reweighting weights across an EnsembleArray{T}, including convergence control parameters and the flattened global weight vector. It is designed to be initialized once per reweighting task and then updated in-place during the solve loop.

Type Parameters

T: Numeric type for floating-point parameters (e.g., Float64)

Constructor

ReweightingSolver(
    ens::Ensemble.EnsembleArray{T},
    maxiter::Int = 500,
    eps::T = 1e-12
) -> ReweightingSolver{T}

Construct a ReweightingSolver{T} for a given ensemble array.

Constructor Arguments

ens: Target EnsembleArray{T} to operate on
maxiter: Maximum number of iterations allowed in the solver
eps: Convergence threshold for the iterative update

Constructor Behavior

Computes nconf_all as the total number of configurations across all ensembles contained in ens
Initializes the global reweighting vector w to a length-nconf_all vector of ones

Constructor Returns

A fully initialized ReweightingSolver{T} instance ready for iteration.

Fields

ens::Ensemble.EnsembleArray{T}: Ensemble array on which multi-ensemble reweighting is performed
maxiter::Int: Maximum number of iterations for solving
eps::T: Convergence threshold
w::Vector{T}: Current flattened reweighting vector (length nconf_all)
nconf_all::Int: Total number of configurations across all ensembles

Notes

The weight vector w is stored in a flattened form that concatenates configurations from all ensembles in ens in their internal order.
This struct is intentionally mutable, as w and internal solver state are updated iteratively during the reweighting procedure.

source

Deborah.Miriam.Reweighting.calc_f! — Method

calc_f!(
    rw::ReweightingSolver{T}, 
    info_file::String, 
    tag::String, 
    jobid::Union{Nothing, String}=nothing
) -> Nothing where T

Solve for the free energy differences ($f_i$) between ensembles using Newton-Raphson iteration, with analytic Jacobian and result logging.

Arguments

rw::ReweightingSolver{T}: Solver instance with ensemble data and parameters.
info_file::String: Path to the TOML file for logging results.
tag::String: Tag label for the TOML solver block.
jobid::Union{Nothing,String}: Optional job identifier for logging context.

Behavior

Initializes $f$ with maximum action differences.
Runs Newton iteration with analytic Jacobian.
Updates ensemble fields ens[i].f in place (except fixed ensemble).
Logs solver details and residuals to TOML.

Returns

Nothing (modifies solver state and logs status in place).

Notes

This function provides the core solver for the MBAR self-consistency equations. For each target ensemble $i = 1, \dots, n_{\mathrm{ens}} - 1$, we define the denominator

\[D_{j;m}^{(i)}(f) = \sum_{k=1}^{n_{\mathrm{ens}}} N_k \; \exp\!\left[ \Delta S_{j;m}(\kappa_i) - \Delta S_{j;m}(\kappa_k) - f_i + f_k \right] \, ,\]

where $j$ labels the source ensemble and $m$ indexes configurations within that ensemble. Using this quantity, the normalization sum for ensemble $i$ is given by

\[T_i(f) = \sum_{j=1}^{n_{\mathrm{ens}}} \sum_{m=1}^{N_j} \frac{1}{D_{j;m}^{(i)}(f)} \, .\]

The MBAR residual equations are then written in the compact form

\[F_i(f) = \log T_i(f) = 0, \qquad i = 1, \dots, n_{\mathrm{ens}} - 1 \, .\]

These nonlinear equations enforce the self-consistency condition $T_i(f) = 1$ for each target ensemble, corresponding to statistical equilibrium among all ensembles in the reweighting framework. The solution determines the optimal free-energy offsets $(f_1, \dots, f_{n_{\mathrm{ens}}-1})$, with one ensemble fixed to set the overall normalization.

Workflow

Initialization Calls set_init_guess! to provide robust starting values
\[f_i^{(0)} = \max_{1 \le j \le n_{\mathrm{ens}}} \left( \max_{1 \le m \le N_j} \Delta S_{j;m}(\kappa_i) \right) \,.\]
ensuring stability of the first Newton step.
Residual and Jacobian Uses fdf! to evaluate the MBAR residual vector $F(f)$ and its analytic Jacobian $J(f) = \dfrac{\partial F}{\partial f}$.
Nonlinear Solve Calls NLsolve.nlsolve to iteratively update
\[f \;\mapsto\; f - J^{-1}(f)\,F(f).\]
Normalization One ensemble's free energy (typically the last) is held fixed so that only $n_{\mathrm{ens}}-1$ unknowns are solved.
Logging Solver diagnostics (convergence, iterations, residual norms, solutions) are appended to the specified TOML file under a tagged section.

source

Deborah.Miriam.Reweighting.calc_w! — Method

calc_w!(
    rw::ReweightingSolver{T}, 
    paramT::Ensemble.Params{T}
) -> Nothing

Compute normalized reweighting weights $w_{j;m}$ for all configurations, given a target parameter set $\theta_T$.

Arguments

rw::ReweightingSolver{T}: The solver containing ensemble data and offsets $f_k$.
paramT::Ensemble.Params{T}: Target parameter set for$\kappa_T$.

Returns

Nothing (fills rw.w in place with weights of length rw.nconf_all).

Notes

Each configuration $\mathcal{U}_{j;m}$ with ensemble $j$, gauge configuration index $m$ is assigned a weight that incorporates action differences and the previously determined free energy offsets $f_k$. These weights allow expectation values at $\kappa_T$ to be estimated from configurations generated at multiple ensembles.

Formula

For a configuration $\mathcal{U}_{j;m}$, the reweighting weight is given by

\[ w_{j;m}(\kappa_T) = \frac{\exp \!\left[ - \Delta S_{j;m}(\kappa_T) \right]}{ \displaystyle\sum_{k=1}^{n_{\mathrm{ens}}} N_k \; \exp \!\left[ - \Delta S_{j;m}(\kappa_k) + f_k - \mathcal{X} \right] } \,,\]

where

$N_k$: number of configurations in ensemble $k$,
$f_k$: free energy offset for ensemble $k$.
$\mathcal{X}$: stabilizer of the exponentials for numerical safety (see below).
$\Delta S_{j;m}(\kappa_k) \equiv S(\kappa_k; \mathcal{U}_{j;m}) - S(\kappa_j;\mathcal{U}_{j;m})$: action difference evaluated on configuration $m$ from ensemble $j$ under parameter $\kappa_k$ where we note that

\[ \exp \left[ -\,\Delta S_{j;m}(\kappa_k) \right] = \left[ \left( \frac{\det M(\kappa_k)}{\det M(\kappa_j)} \right)^{N_{\text{f}}} \right]_{j;m} \,.\]

Numerical Stability

Direct evaluation of exponentials can cause overflow or underflow. To stabilize, the implementation uses the log-sum-exp trick: first determine the global maximum exponent across all configurations

\[\mathcal{X} = \max_{1 \le i,j \le n_{\mathrm{ens}}} \left( \max_{1 \le m \le N_j} \left\{ \Delta S_{i;m}(\kappa_T) - \Delta S_{i;m}(\kappa_j) + f_j \right\} \right) \,, \]

and then subtract this value before exponentiation:

\[\sum_{j=k}^{n_{\mathrm{ens}}} N_k \exp \left[ \Delta S_{j;m}(\kappa_T) - \Delta S_{j;m}(\kappa_k) + f_j - \mathcal{X} \right].\]

Since every term in the denominator is shifted by the same constant ($-\mathcal{X}$), the final weights

\[w_{j;m}(\kappa_T) = \frac{1}{\displaystyle{\sum_k N_k \exp[\cdots - \mathcal{X}]}}\]

remain unchanged. Only the intermediate exponential evaluations are affected, which makes the computation numerically safe.

Usage in Observables

The target expectation value is then estimated as

\[ \left\langle \Omega(\kappa_T;\mathcal{U}) \right\rangle = \frac{\displaystyle \sum_{j=1}^{n_{\mathrm{ens}}} \sum_{m=1}^{N_j} w_{j;m}(\kappa_T) \; \Omega(\kappa_T;\mathcal{U}_{j;m})} {\displaystyle \sum_{j=1}^{n_{\mathrm{ens}}} \sum_{m=1}^{N_j} w_{j;m}(\kappa_T)} \, .\]

source

Deborah.Miriam.Reweighting.fdf! — Method

fdf!(
    F::Vector{T}, 
    f::Vector{T}, 
    rw::ReweightingSolver{T}
) -> (Vector{T}, Matrix{T})

Compute the residual vector F and Jacobian matrix J for the reweighting equations in a system of n ensembles, using the current guess f for the free energy offsets.

Arguments

F::Vector{T}: Output residual vector (size $n_{\mathrm{ens}} - 1$)
f::Vector{T}: Current guess for free energy offsets (size $n_{\mathrm{ens}} - 1$)
rw::ReweightingSolver{T}: Reweighting solver containing ensemble data

Returns

F::Vector{T}: Updated residual vector
J::Matrix{T}: Computed Jacobian matrix (size $(n_{\mathrm{ens}} - 1) \times (n_{\mathrm{ens}} - 1)$)

Notes

This is part of a Newton-Raphson solver for determining optimal reweighting factors.

The equations implemented here are mathematically identical to the MBAR (Multistate Bennett Acceptance Ratio) self-consistency equations. MBAR provides the statistically optimal (minimum-variance, asymptotically unbiased) estimator for the relative free energies of multiple ensembles, given finite sampling from each.

In the Lattice QCD literature, the same formalism is more often referred to simply as multi-ensemble reweighting or Ferrenberg-Swendsen reweighting rather than by the name MBAR, but the underlying equations coincide.

Problem Setup

We have $n_{ens}$ ensembles indexed by $k=1,\dots,n_{\mathrm{ens}}$. Each ensemble has:

Number of configurations: $N_k = \texttt{ens[k].nconf}$,
Parameter (e.g. coupling, potential parameter): $\theta_k = \texttt{ens[k].param}$. Here we mainly perform reweighting in $\kappa$. Hence, we write $\kappa_k$ for clarity hereafter.

For each configuration $\mathcal{U}_{j;m}$ drawn from ensemble $j$, (where $m$ denotes the gauge configuration index), we can evaluate the reduced action (or potential) under another parameter $\kappa_k$:

\[\Delta S_{j;m}(\kappa_k) \equiv S(\kappa_k; \mathcal{U}_{j;m}) - S(\kappa_j;\mathcal{U}_{j;m})\]

We fix one ensemble (say $k=n_{\mathrm{ens}}$, the last index) as the reference and solve for the free energy offsets

\[\mathbf{f} = (f_1, f_2, \dots, f_{n_{\mathrm{ens}} - 1}).\]

Residual Equations

For each target ensemble $i = 1,\dots,n_{\mathrm{ens}} - 1$, define

\[D_{j;m}^{(i)}(f) = \sum_{k=1}^{n_{\mathrm{ens}}} N_k \exp\!\left[ E_{j;m}^{(i,k)} \right].\]

where

\[E_{j;m}^{(i,k)} = \Delta S_{j;m}(\kappa_i) - \Delta S_{j;m}(\kappa_k) - f_i + f_k \,.\]

Then define the normalization sum

\[T_i(f) = \sum_{j=1}^{n_{\mathrm{ens}}}\; \sum_{m=1}^{N_j} \frac{1}{D_{j;m}^{(i)}(f)}.\]

The residual equations are

\[F_i(f) = \log T_i(f), \qquad i=1,\dots,n_{\mathrm{ens}} - 1.\]

The system of equations to be solved is

\[F_i(f) = 0, \qquad i=1,\dots,n_{\mathrm{ens}} - 1,\]

which enforces

\[T_i(f) = 1.\]

This is precisely the MBAR self-consistency condition for the free energies.

Jacobian Matrix

We compute the Jacobian entries

\[J_{i \ell}(f) = \frac{\partial F_i}{\partial f_\ell}.\]

From the definition,

\[\frac{\partial F_i}{\partial f_\ell} = \frac{1}{T_i(f)} \sum_{j,m} \left( -\frac{1}{(D_{j;m}^{(i)}(f))^2} \frac{\partial D_{j;m}^{(i)}}{\partial f_\ell} \right).\]

The derivative of the denominator is

\[\frac{\partial D_{j;m}^{(i)}}{\partial f_\ell} = N_\ell \; \exp\!\left[ E_{j;m}^{(i,\ell)} \right] - \delta_{i\ell}\,D_{j;m}^{(i)}(f).\]

where

\[E_{j;m}^{(i,\ell)} = \Delta S_{(j \to i);m} - \Delta S_{(j \to \ell);m} - f_i + f_\ell.\]

Substituting back,

\[J_{i \ell}(f) = \frac{1}{T_i(f)}\sum_{j,m} \left( -\frac{N_\ell \; \exp\!\left[E_{j;m}^{(i,\ell)}\right]}{\left(D_{j;m}^{(i)}(f)\right)^2} \right) + \delta_{i\ell} \,.\]

Thus:

Off-diagonal entries ($i \neq \ell$) come only from the exponential ratio term,
Diagonal entries ($i=\ell$) acquire an additional $+1$.

Summary

Unknowns: free energy offsets $f_1,\dots,f_{n_{\mathrm{ens}} - 1}$,
Equations: $F_i(f)=0$, enforcing MBAR consistency,
Residual vector:
\[F_i(f) = \log\left(\sum_{j,m} \frac{1}{D_{j;m}^{(i)}(f)}\right),\]
Jacobian matrix:
\[J_{i\ell}(f) = \frac{1}{T_i(f)} \sum_{j,m} \left(-\frac{N_\ell \; \exp\!\left[E_{j;m}^{(i,\ell)}\right]}{(D_{j;m}^{(i)}(f))^2}\right) + \delta_{i\ell}.\]

This system is solved iteratively using Newton-Raphson updates:

\[f \;\mapsto\; f - J^{-1}(f)\,F(f).\]

source

Deborah.Miriam.Reweighting.set_init_guess! — Method

set_init_guess!(
    rw::ReweightingSolver{T}
) -> Nothing where T

Initialize the free energy estimates f in the given ReweightingSolver.

Arguments

rw::ReweightingSolver{T}: The solver object containing all ensembles

Returns

Nothing (modifies ensemble fields in-place)

Notes

For each ensemble $i$, the initial value $f_i$ is set to the maximum action difference Deborah.Miriam.EnsembleUtils.dS across all configurations $j$ from all ensembles, evaluated under the parameter $\theta_i$, or $\kappa_i$ more specifically:

\[f_i^{(0)} = \max_{1 \le j \le n_{\mathrm{ens}}} \left( \max_{1 \le m \le N_j} \Delta S_{j;m}(\kappa_i) \right) \,, \qquad \Delta S_{j;m}(\kappa_i) \equiv S(\kappa_i; \mathcal{U}_{j;m}) - S(\kappa_j;\mathcal{U}_{j;m}) \,.\]

This provides a robust starting point for the Newton iteration used in multi-ensemble reweighting (MBAR-type equations). Choosing the maximum value ensures numerical stability by avoiding underestimation of the exponential weights.

source