Ensemble reweighting (frame weights)

Overview

By default every analysis routine in SOURSOP treats all frames of a trajectory as contributing equally to an ensemble average. In many situations — enhanced-sampling simulations, Markov-state-model reweighting, maximum-entropy / experimentally-restrained reweighting, or simple importance reweighting — each frame should instead contribute according to a statistical weight.

Every SOURSOP function that returns an ensemble-average value therefore accepts an optional weights keyword: a per-frame vector that defines how much each frame contributes to the average. The default, weights=False, is an exact no-op and reproduces the original unweighted behaviour bit-for-bit.

Such a weight vector can be supplied from any external source (an MSM, an enhanced-sampling estimator, importance weights), or generated within SOURSOP by reweighting the ensemble against experimental data. Two complementary maximum-entropy reweighters are provided:

• BME / iBME — Bayesian Maximum Entropy, a tunable penalty balancing fit and diversity (see ssbme); and

• COPER / iCOPER — Convex Optimization for Ensemble Reweighting, a hard chi-squared constraint with no free regularisation parameter (see sscoper).

The weights returned by either satisfy the contract below and can be passed straight into any reweighting-capable method.

Reweighting in SOURSOP is deterministic. A weighted average is the closed-form expectation under the supplied weights — there is no stochastic resampling, so results are exactly reproducible and do not depend on a random seed.

The weight vector contract

A valid weights vector is a probability vector over frames:

• one entry per frame (len(weights) == n_frames);

• every entry in the closed interval [0, 1];

• all entries finite (no nan / inf);

• the entries sum to 1 within a small tolerance etol (default 1e-7).

These conditions are enforced centrally by soursop.ssutils.validate_weights(). If any condition fails an SSException is raised with a descriptive message, rather than silently producing a meaningless number.

Most reweighting-capable methods also expose an etol keyword so the sum-to-one tolerance can be tightened or relaxed.

Interaction with stride

When a method is called with both a frame stride and weights, the weight vector is first subsampled (weights[::stride]) and then renormalised so it still sums to 1 over the retained frames. A warning is emitted because per-stride reweighting is rarely what you want unless the weights were computed for exactly those frames.

The deterministic helpers

All weighting is funnelled through a small set of shared helpers in soursop.ssutils so that every method applies weights identically:

• validate_weights() — validate (and stride-normalise) a weight vector; the single source of truth.

• weighted_mean() — sum(w * x).

• weighted_rms() — sqrt(sum(w * x**2)) (the polymer-physics RMS order parameter).

• weighted_std() — the reliability-weighted population standard deviation. Frame weights are probability weights with no associated sample size, so no ddof correction is well-defined; the population estimator is the unambiguous, reproducible choice.

• weighted_corr() — the weighted Pearson correlation (used by the correlation observables).

Example

import numpy as np
from soursop.sstrajectory import SSTrajectory

TrajOb  = SSTrajectory('traj.xtc', 'start.pdb')
protein = TrajOb.proteinTrajectoryList[0]

n = protein.n_frames

# 1. default: every frame contributes equally
rg_per_frame = protein.get_radius_of_gyration()        # array, length n

# 2. uniform weights == the ordinary mean of the per-frame values
w_uniform = np.full(n, 1.0 / n)
rg_mean   = protein.get_radius_of_gyration(weights=w_uniform)
assert np.isclose(rg_mean, rg_per_frame.mean())

# 3. a one-hot weight isolates a single frame
w_one      = np.zeros(n); w_one[10] = 1.0
rg_frame10 = protein.get_radius_of_gyration(weights=w_one)
assert np.isclose(rg_frame10, rg_per_frame[10])

# 4. a real reweighting vector (e.g. from MSM / MaxEnt), here from BME
from soursop.ssbme import BME, ExperimentalObservable
obs  = [ExperimentalObservable(value=23.0, uncertainty=1.0, name="Rg")]
calc = protein.get_radius_of_gyration().reshape(-1, 1)   # (n, 1)
w    = BME(obs, calc).fit(theta=2.0, auto_theta=False).weights
dmap_w, _ = protein.get_distance_map(weights=w)
nu, A0    = protein.get_scaling_exponent(weights=w)[:2]

Behaviour at the extremes

• Uniform weights (1/n everywhere) reproduce the ordinary (unweighted) ensemble mean to within floating-point precision.

• One-hot weights (a single 1.0) return exactly that frame’s value.

• Invalid weights (wrong length, an element outside [0, 1], a non-finite element, or a sum that differs from 1 by more than etol) raise an SSException.

What can be reweighted

Any method that collapses the trajectory to an ensemble value accepts weights. This includes the global dimensions (get_radius_of_gyration, get_hydrodynamic_radius, get_asphericity, get_end_to_end_distance, get_gyration_tensor), the per-frame getters’ ensemble means, the distance/contact maps, get_Q, the polymer-scaling observables (get_internal_scaling, get_internal_scaling_RMS, get_scaling_exponent, get_local_to_global_correlation), the SASA summaries (get_all_SASA, get_site_accessibility, get_regional_SASA), get_angle_decay, get_local_collapse, get_end_to_end_vs_rg_correlation, the dihedral mutual information, and the SSTrajectory overall/inter-chain observables.

A small number of routines describe a distribution or a pairwise frame-vs-frame quantity rather than a single ensemble average — for example get_internal_scaling(mean_vals=False) and get_local_heterogeneity. A single per-frame probability vector is not well-defined for those, so they raise an SSException if weights is supplied (instead of silently returning a questionable number).

API reference

soursop.ssutils.validate_weights(weights, n_frames, stride=1, etol=1e-07)

Validate (and stride-normalise) a per-frame re-weighting vector.

This is the single, shared entry point used by every SOURSOP function that accepts a weights keyword. It enforces that the weights form a proper per-frame probability vector so that all downstream deterministic weighted averages (weighted_mean(), weighted_std(), weighted_rms(), weighted_corr()) are well defined and consistent.

The literal False (or None) is the “no weighting” sentinel and is passed straight through, so a default of weights=False is a strict no-op.

Validation (each failure raises SSException):

1. False / None -> returned unchanged (no-op).

2. castable to a 1-D numpy.float64 array.

3. all elements finite (no nan / inf).

4. exactly one weight per frame (len == n_frames).

5. every element in the closed interval [0, 1].

6. if stride > 1: the vector is subsampled (weights[::stride]) and renormalised to sum to 1 so the strided weighted average is still a proper expectation (this fixes the historical behaviour where strided weights silently no longer summed to 1).

7. |sum(weights) - 1| < etol (checked after any stride/renormalise).

Parameters:

• weights (array_like, False or None) – Per-frame weights, or the False/None no-op sentinel.

• n_frames (int) – Number of frames the (unstrided) weight vector must match.

• stride (int, optional) – Frame stride that will be applied to the trajectory. Default 1.

• etol (float, optional) – Tolerance on |sum(weights) - 1|. Default 1e-7.

Returns:

False if the input was False/None; otherwise a numpy.float64 array of length ceil(n_frames / stride) that sums to 1 within etol.

Return type:

numpy.ndarray or False

Raises:

SSException – If any of the validation conditions above fail.

Example

>>> import numpy as np
>>> from soursop.ssutils import validate_weights
>>> validate_weights(False, 10) is False
True
>>> w = validate_weights(np.full(10, 0.1), 10)
>>> float(w.sum())
1.0

soursop.ssutils.weighted_mean(x, weights, axis=0)

Deterministic frame-weighted mean (sum(w * x) along axis).

Thin wrapper around numpy.average(). Because weights is a validated probability vector (sum == 1) this is exactly the re-weighted ensemble expectation and is numerically identical to the numpy.average(x, axis, weights=weights) calls already used elsewhere in SOURSOP.

Parameters:

• x (array_like) – Per-frame data; the frame axis is axis.

• weights (numpy.ndarray) – Validated per-frame weights (see validate_weights()).

• axis (int, optional) – Axis to average over (the frame axis). Default 0.

Returns:

x averaged over axis with the given weights.

Return type:

numpy.ndarray or float

Example

>>> import numpy as np
>>> weighted_mean(np.array([1.0, 3.0]), np.array([0.25, 0.75]))
2.5

soursop.ssutils.weighted_rms(x, weights, axis=0)

Deterministic frame-weighted root-mean-square (sqrt(sum(w*x^2))).

The polymer-physics order parameter used by the internal-scaling and scaling-exponent routines is an RMS distance, so this is the weighted analogue of sqrt(mean(x**2)).

Parameters:

• x (array_like) – Per-frame data; the frame axis is axis.

• weights (numpy.ndarray) – Validated per-frame weights (see validate_weights()).

• axis (int, optional) – Axis to reduce over. Default 0.

Returns:

The weighted RMS of x over axis.

Return type:

numpy.ndarray or float

Example

>>> import numpy as np
>>> float(weighted_rms(np.array([3.0, 4.0]), np.array([0.5, 0.5])))
3.5355339059327378

soursop.ssutils.weighted_std(x, weights, axis=0)

Deterministic frame-weighted (population) standard deviation.

Uses the reliability-weighted population estimator sqrt(sum(w * (x - mean)**2)). Frame weights here are probability weights with no associated sample size, so there is no well-defined ddof (Bessel-style) correction; the population estimator is the unambiguous, reproducible choice and the de-facto standard for re-weighted molecular-dynamics ensembles.

Parameters:

• x (array_like) – Per-frame data; the frame axis is axis.

• weights (numpy.ndarray) – Validated per-frame weights (see validate_weights()).

• axis (int, optional) – Axis to reduce over. Default 0.

Returns:

The weighted population standard deviation of x over axis.

Return type:

numpy.ndarray or float

Example

>>> import numpy as np
>>> float(weighted_std(np.array([1.0, 1.0]), np.array([0.5, 0.5])))
0.0

soursop.ssutils.weighted_corr(a, b, weights)

Deterministic frame-weighted Pearson correlation between two vectors.

Computed from the weighted covariance matrix (numpy.cov(..., ddof=0, aweights=weights)), matching the ddof=0 / aweights convention already used by get_local_to_global_correlation.

Parameters:

• a (array_like) – Equal-length per-frame vectors.

• b (array_like) – Equal-length per-frame vectors.

• weights (numpy.ndarray) – Validated per-frame weights (see validate_weights()).

Returns:

The weighted Pearson correlation coefficient of a and b.

Return type:

float

Example

>>> import numpy as np
>>> w = np.full(4, 0.25)
>>> round(float(weighted_corr(np.array([1.,2,3,4]), np.array([2.,4,6,8]), w)), 6)
1.0

Shared reweighting primitives

These helpers in soursop.ssutils are shared by the BME (ssbme) and COPER (sscoper) reweighters so that both expose an identical interface.

class soursop.ssutils.ExperimentalObservable(value: float, uncertainty: float, constraint: str = 'equality', name: str | None = None, group: str | None = None)

Container for a single experimental observable.

Shared by soursop.ssbme and soursop.sscoper (both re-export this class), so the user-facing syntax is identical for BME, iBME, COPER and iCOPER.

Parameters:

• value (float) – The experimental value of the observable.

• uncertainty (float) – The experimental uncertainty (standard deviation). Must be positive.

• constraint (str, optional) –

  Type of constraint, one of:

  • "equality" (default): observable should match value within uncertainty.

  • "upper": observable should not exceed value (deviations below value are not penalized).

  • "lower": observable should not fall below value (deviations above value are not penalized).

• name (str, optional) – Optional human-readable name/description.

• group (str, optional) – Optional data-type label. Used by soursop.sscoper.COPER to impose a separate per-group chi-squared constraint (chi2_alpha <= limit for each group, as in Leung et al. 2016); ignored by BME / iBME. Observables without a group are pooled into a single default group.

Raises:

SSException – If uncertainty is not positive or constraint is invalid.

get_bounds() → Tuple[float | None, float | None]

Optimization bounds on the Lagrange multiplier for this observable.

Returns:

(None, None) for equality, (0.0, None) for upper, (None, 0.0) for lower.

Return type:

tuple

soursop.ssutils.relative_entropy(w0, w1)

Relative entropy (Kullback-Leibler divergence) of w1 from w0.

Parameters:

• w0 (numpy.ndarray) – Reference (prior) weights, normalized to sum to 1.

• w1 (numpy.ndarray) – Posterior weights, normalized to sum to 1.

Returns:

sum_i w1_i * log(w1_i / w0_i) over frames with non-negligible posterior weight.

Return type:

float

soursop.ssutils.weighted_linear_regression(x, y, sample_weight, fit_intercept=True)

Closed-form weighted least-squares regression of y on x.

A small numpy replacement for sklearn.linear_model.LinearRegression (SOURSOP does not depend on scikit-learn). Solves min_{a,b} sum_i s_i (y_i - (a x_i + b))^2. Used by the iterative scale/offset reweighters (iBME, iCOPER).

Parameters:

• x (numpy.ndarray) – Independent variable, shape (n,).

• y (numpy.ndarray) – Dependent variable, shape (n,).

• sample_weight (numpy.ndarray) – Per-sample weights, shape (n,).

• fit_intercept (bool, optional) – If True fit slope and intercept; if False force the intercept to zero (slope only). Default True.

Returns:

(slope, intercept). intercept is 0.0 when fit_intercept is False.

Return type:

tuple of float

soursop.ssutils.find_optimal_theta(x_values, y_values, method='perpendicular')

Select the L-curve knee from two paired metric arrays.

Generic knee-finder shared by BME’s theta_scan (chi-squared vs. relative entropy across theta) and COPER’s chi2_limit_scan (chi-squared vs. relative entropy across the chi-squared limit).

Parameters:

• x_values (numpy.ndarray) – First metric per scan point (e.g. final chi-squared).

• y_values (numpy.ndarray) – Second metric per scan point (e.g. relative entropy).

• method (str, optional) – "perpendicular" (default) or "curvature".

Returns:

(optimal_idx, method_name).

Return type:

tuple

Raises:

SSException – If method is unknown.

soursop.ssutils.validate_reweighting_inputs(observables, calculated_values, initial_weights)

Validate constructor arguments shared by the reweighter classes.

Used by soursop.ssbme.BME / iBME and soursop.sscoper.COPER / iCOPER.

Parameters:

• observables (list of ExperimentalObservable) – Experimental observables (non-empty).

• calculated_values (numpy.ndarray) – Per-frame calculated values, shape (n_frames, n_observables).

• initial_weights (numpy.ndarray or None) – Optional prior weights, one per frame.

Raises:

SSException – If any argument is malformed or dimensions are inconsistent.

soursop.ssutils.constraint_chi_squared(weights, calculated_values, observables, indices=None)

Constraint-aware reduced chi-squared for a weight vector.

equality observables always penalize deviations; upper / lower only penalize the disallowed side. Shared by BME (total chi-squared) and COPER (per-data-type chi-squared via indices).

Parameters:

• weights (numpy.ndarray) – Frame weights, shape (n_frames,).

• calculated_values (numpy.ndarray) – Per-frame calculated values, shape (n_frames, n_observables).

• observables (list of ExperimentalObservable) – The experimental observables (columns of calculated_values).

• indices (sequence of int, optional) – Restrict the chi-squared to this subset of observable columns (used for COPER per-group constraints). Defaults to all columns.

Returns:

Mean of (diff / sigma)^2 over the selected observables.

Return type:

float