Ensemble comparison

Tools to statistically compare distributions of a feature (or a 2-column joint feature) across systems and/or replicas. Each compute_* function is also exposed as a Dataset method of the same name.

These tools group by Dataset.unique_replica_col rather than the raw replica column, so comparisons stay correct even when replica labels (e.g. plain numbers 1, 2, 3…) are reused identically across systems – see Dataset.

Jensen-Shannon divergence

Compares two distributions via their Jensen-Shannon divergence, a symmetric, bounded ([0, ln(2)]) measure well suited to answering “are these two systems sampling different states?” – with bootstrap confidence intervals and an automatic, human-readable verdict.

RepLikCompare.EnsembleComparison.compute_jsd.compute_jsd(df, sim_name_col, replica_col, feature_col, group_by='replica', is_angle=False, angle_unit='auto', bins='auto', iters=200, frac=1.0, replace=True, seed=42, confidence_level=0.95, ratio_bootstrap=5000, ratio_seed=0, outlier_z_threshold=1.2, plot_heatmap=False, heatmap_kwargs=None, show_values=True, verbose=True, explain=True, explain_kwargs=None)[source]

Compare N systems x M replicas on an already-computed feature (1 or 2 columns), via the Jensen-Shannon Divergence (JSD), with bootstrap confidence intervals. No hypothesis test / p-value: only bootstrap means, percentile CIs, and descriptive magnitude/verdict heuristics.

Parameters:

df (pandas.DataFrame) – “Long” dataframe (one row per frame), with the sim_name_col / replica_col columns and the feature column(s).
sim_name_col (str) – Name of the column identifying the system (e.g. “A”, “B”).
replica_col (str) – Name of the column identifying the replica, unique per system. Only used when group_by=”replica” (the default); ignored (but still required) when group_by=”system”.
feature_col (str or (str, str)) – Name of the column to compare (1D JSD), or a tuple of 2 column names for a joint 2D JSD (e.g. (“PC1”, “PC2”)).
group_by ({"replica", "system"}, optional) – “replica” (default): the rigorous mode. Groups by replica_col, so intra-system (replica-to-replica) variability is computed and used as the noise floor against which the inter-system JSD is judged (see summarize_systems – produces verdict, effect_magnitude, and the full replica-level detail tables). “system”: a quicker, pooled mode. Groups by sim_name_col directly – all replicas of a system are fused into one series before histogramming. Produces only a direct inter-system JSD (see summarize_systems_pooled); no intra-system spread, no verdict, no replica-level detail tables (replica_pairs/ replica_outliers/system_overall are omitted). Useful for a fast overview, but more optimistic than “replica” mode: it ignores that some replicas may diverge from each other.
is_angle (bool or (bool, bool), optional) – True if the column is an angle/dihedral (circular treatment: fixed JSD range [-pi, pi], after conversion to radians if needed). Tuple (bool, bool) to specify each column independently in 2D. Default False.
angle_unit (str or (str, str), optional) – “auto” (automatic degrees/radians detection from value amplitude), “deg” or “rad” to force. Default “auto”.
bins (str or int or (int, int), optional) – “auto”/”sqrt”/”sturges”, a fixed integer, or (bins_x, bins_y) in 2D. Default “auto”.
iters (int, optional) – Number of bootstrap draws. Default 200.
frac (float, optional) – Fraction of each group resampled at each draw. Default 1.0.
replace (bool, optional) – Resample with replacement. Default True.
seed (int, optional) – Seed for the main bootstrap. Default 42.
confidence_level (float, optional) – Confidence level for percentile CIs. Default 0.95.
ratio_bootstrap (int, optional) – Number of draws for the inter/intra ratio CI. Only used in group_by=”replica” mode. Default 5000.
ratio_seed (int, optional) – Seed for the ratio bootstrap. Default 0.
outlier_z_threshold (float, optional) – deviation_score threshold above which a replica is flagged “worth_checking”. Only used in group_by=”replica” mode. Default 1.2.
plot_heatmap (bool, optional) – If True, plot the JSD score heatmap (added to the return dict under the “heatmap_ax” key). Default False.
heatmap_kwargs (dict, optional) – Extra arguments passed to plot_jsd_heatmap (e.g. triangular=”lower”, context=”talk”, ax=…).
show_values (bool, optional) – If True, display the mean JSD value and bootstrap confidence interval in each heatmap cell. If False, hide the annotations for a more compact heatmap.
verbose (bool, optional) – Show progress (detected angle unit, resolved bins, bootstrap progress bar if tqdm is installed). Default True.
explain (bool, optional) – If True (default), automatically call explain_jsd(tables) at the end and print a formatted, human-readable interpretation (verdicts / magnitudes / outlier replicas / table legend) right below the log output and the heatmap. Set to False for a silent call that only returns the dict (e.g. in a loop over many features). Set to True and use explain_kwargs={“top_outliers”: N} to only show the N most extreme outlier replicas.
explain_kwargs (dict, optional) – Extra arguments passed to explain_jsd when explain=True (currently just top_outliers).

Returns:

tables (dict) –

“system_summary” : system-vs-system comparison. In group_by=”replica” mode: intra/inter/verdict/effect size (see summarize_systems). In group_by=”system” mode: pooled inter-system JSD only (see summarize_systems_pooled).
”replica_pairs” : replica-vs-replica comparison within each system (group_by=”replica” only).
”replica_outliers” : per-replica heuristic (deviation_score/ assessment) (group_by=”replica” only).
”system_overall” : aggregated verdict across features, only if >= 2 systems (group_by=”replica” only).
”heatmap_ax” : heatmap Axes (only if plot_heatmap=True).
”__description__” : description of each table above.

Wasserstein distance

Compares two distributions via the (first-order) Wasserstein/ earth-mover’s distance – unbounded, expressed in the feature’s native units, and computed directly on the raw samples rather than a histogram. Useful as a complement to the JSD when the magnitude of the shift between distributions matters, not just whether it’s detectable.

RepLikCompare.EnsembleComparison.compute_wasserstein.compute_wasserstein(df, sim_name_col, replica_col, feature_col, group_by='replica', is_angle=False, angle_unit='auto', iters=200, frac=1.0, replace=True, seed=42, confidence_level=0.95, ratio_bootstrap=5000, ratio_seed=0, outlier_z_threshold=1.2, plot_heatmap=False, heatmap_kwargs=None, verbose=True, explain=True, explain_kwargs=None)[source]

Compare N systems x M replicas on an already-computed feature (1 column – 2D Wasserstein is not implemented, use compute_jsd for a 2-column feature), via the (1st order) Wasserstein distance, with bootstrap confidence intervals. No hypothesis test / p-value.

Unlike the JSD, the Wasserstein distance is NOT bounded: it is expressed in the feature’s own physical units, computed exactly from the raw samples (no histogram/binning step at all – there is no “number of bins” to report or tune for this metric), and has no universal absolute-magnitude scale (only the inter/intra ratio is reported, no “jsd_magnitude”-style column).

Parameters:

df (pandas.DataFrame) – “Long” dataframe (one row per frame), with the sim_name_col / replica_col columns and the feature column.
sim_name_col (str) – Name of the column identifying the system (e.g. “A”, “B”).
replica_col (str) – Name of the column identifying the replica, unique per system. Only used when group_by=”replica” (the default); ignored (but still required) when group_by=”system”.
feature_col (str) – Name of the column to compare.
group_by ({"replica", "system"}, optional) – “replica” (default): the rigorous mode. Groups by replica_col, so intra-system (replica-to-replica) variability is computed and used as the noise floor against which the inter-system distance is judged (see summarize_systems – produces verdict, effect_magnitude, and the full replica-level detail tables). “system”: a quicker, pooled mode. Groups by sim_name_col directly – all replicas of a system are fused into one series before computing the distance. Produces only a direct inter-system distance (see summarize_systems_pooled); no intra-system spread, no verdict, no replica-level detail tables (replica_pairs/replica_outliers/system_overall are omitted). Useful for a fast overview, but more optimistic than “replica” mode: it ignores that some replicas may diverge from each other.
is_angle (bool, optional) – True if the column is an angle/dihedral (circular approximation of the Wasserstein distance via a rotation grid search). Default False.
angle_unit (str, optional) – “auto”/”deg”/”rad”. Default “auto”.
iters (int, optional) – Number of bootstrap draws. Default 200.
frac (float, optional) – Fraction of each group resampled at each draw. Default 1.0.
replace (bool, optional) – Resample with replacement. Default True.
seed (int, optional) – Seed for the main bootstrap. Default 42.
confidence_level (float, optional) – Confidence level for percentile CIs. Default 0.95.
ratio_bootstrap (int, optional) – Number of draws for the inter/intra ratio CI. Only used in group_by=”replica” mode. Default 5000.
ratio_seed (int, optional) – Seed for the ratio bootstrap. Default 0.
outlier_z_threshold (float, optional) – deviation_score threshold above which a replica is flagged “worth_checking”. Only used in group_by=”replica” mode. Default 1.2.
plot_heatmap (bool, optional) – If True, plot the Wasserstein score heatmap (added to the return dict under the “heatmap_ax” key). Default False.
heatmap_kwargs (dict, optional) – Extra arguments passed to plot_wasserstein_heatmap.
verbose (bool, optional) – Show progress. Default True.
explain (bool, optional) – If True (default), automatically call explain_wasserstein(tables) at the end and print a formatted, human-readable interpretation (verdicts / outlier replicas / table legend) right below the log output and the heatmap. Set to False for a silent call that only returns the dict (e.g. in a loop over many features). Use explain_kwargs={“top_outliers”: N} to only show the N most extreme outlier replicas.
explain_kwargs (dict, optional) – Extra arguments passed to explain_wasserstein when explain=True (currently just top_outliers).

Returns:

tables (dict) –

“system_summary” : system-vs-system comparison. In group_by=”replica” mode: intra/inter/verdict/effect ratio (see summarize_systems). In group_by=”system” mode: pooled inter-system distance only (see summarize_systems_pooled).
”replica_pairs” : replica-vs-replica comparison within each system (group_by=”replica” only).
”replica_outliers” : per-replica heuristic (deviation_score/ assessment) (group_by=”replica” only).
”system_overall” : aggregated verdict across features, only if >= 2 systems (group_by=”replica” only).
”heatmap_ax” : heatmap Axes (only if plot_heatmap=True).
”__description__” : description of each table above.

Convergence diagnostics

Checks whether a single trajectory has been run long enough for a given feature to look converged – a per-trajectory diagnostic, answering “has this run converged?” rather than “do these two runs agree?” (which is what compute_jsd/compute_wasserstein are for).

RepLikCompare.EnsembleComparison.compute_convergence.compute_convergence(df, feature_col, time_col, sim_name_col=None, replica_col=None, n_blocks=10, fractions=(0.25, 0.5, 0.75, 1.0), bins='auto', plot=False, **kwargs)[source]

Assess whether feature_col has converged over the course of the trajectory. See the module docstring for the three diagnostics computed (cumulative stats at fixed fractions, first-half vs second-half divergence, and block-averaging + growing-window JSD).

Parameters:

df (pandas.DataFrame) – Long-format dataframe, one row per frame.
feature_col (str) – Scalar feature to track (e.g. RMSD, radius of gyration, a collective variable). Must already be numeric – for circular quantities, precompute a suitable projection first.
time_col (str) – Column giving frame order/time; df is sorted by this column (ascending) within each group before slicing into fractions/ blocks/halves.
sim_name_col (str or None, default None) – If given, convergence is computed independently per (sim_name, replica) group. If both are None, the whole df is treated as a single trajectory.
replica_col (str or None, default None) – If given, convergence is computed independently per (sim_name, replica) group. If both are None, the whole df is treated as a single trajectory.
n_blocks (int, default 10) – Number of contiguous blocks used for the finer-grained block-averaging / growing-window-JSD-to-full diagnostic (#3). Trajectories shorter than n_blocks frames are skipped (see block_table[“note”]).
fractions (tuple of float, default (0.25, 0.5, 0.75, 1.0)) – Cumulative fractions of the trajectory (ordered by time_col) at which to report mean/std of feature_col, e.g. the default reports stats over frames [0:25%], [0:50%], [0:75%], [0:100%]. Values must be in (0, 1], increasing.
bins (int or "auto", default "auto") – Histogram binning used for both the first-half-vs-second-half divergence and the growing-window-JSD-to-full diagnostic, forwarded to numpy.histogram. A shared bin edge set (from the full trajectory’s range) is used throughout a group, so all histograms stay directly comparable.
plot (bool, default False) – If True (and matplotlib is installed), also return a three-panel figure per group: cumulative mean +/- std, first-half vs second-half histogram overlay, and growing-window JSD-to-full decay.
**kwargs – Currently unused; reserved for future diagnostics (e.g. autocorrelation time estimation).

Returns:

dict with keys –

“cumulative_table”pandas.DataFrame: One row per (group, fraction): sim_name, replica (if applicable), fraction, n_frames, mean, std. This is diagnostic #1 – e.g. compare the mean/std columns across increasing fraction values to check stability of RMSD/Rg.
”half_split_table”pandas.DataFrame: One row per group: sim_name, replica (if applicable), n_frames_first_half, n_frames_second_half, mean_first_half, mean_second_half, std_first_half, std_second_half, jsd_half_split, wasserstein_half_split. This is diagnostic #2. Small jsd_half_split / wasserstein_half_split (relative to the feature’s own scale) indicate the two halves sample the same distribution.
”block_table”pandas.DataFrame: One row per (group, block): block_index, block_start_time, block_end_time, mean, std, n_frames. Diagnostic #3 (block-averaging part).
”jsd_table”pandas.DataFrame: One row per (group, block boundary k): n_blocks_included, jsd_to_full. Diagnostic #3 (growing-window part).
”figures”dict[group_key -> (fig, axes)], only present if: plot=True.

Ensemble comparison

Jensen-Shannon divergence

Wasserstein distance

Convergence diagnostics

Notebook example