Analysis

The aimmd.analysis subpackage provides the numerical post-processing used for adaptive binning and for turning sampled paths into rates and free-energy profiles: bin construction and merging, binomial confidence intervals, committor fields by relaxation, rate-estimate extraction from logs, and path-lineage utilities.

Numerical analysis helpers for AIMMD.

This module collects standalone routines used by the trainer and by analysis workflows. Functions operate on NumPy arrays and on PathEnsemble-like objects (AIMMD’s path-sampling convention), without depending on high-level AIMMD classes.

Overview Bins and grid helpers

marginal outer bins at -inf and/or +inf).

plus and minus infinity.

ones moving away from the transition state.

Extremes used to define bins

from free excursions started from each terminal state.

transition paths by sampling values adjacent to the endpoints.

Simple statistics

a two-sided confidence interval using Beta quantiles.

Rate-estimate extraction

logs and returns the time series of forward/backward rate estimates.

The parser is tailored to AIMMD train*.log files. It scans the file in order and recognizes a three-line pattern:

  • a line containing "k12 estimate": extracts the forward estimate

following "estimate:" (up to the next "[" token),

  • the next relevant line: extracts the backward estimate in the same way,

  • the next relevant line: extracts the time coordinate as the number

preceding the token "frames".

The returned arrays (t, k12, k21) therefore reflect the sequence of estimates printed by the training logger and the corresponding cumulative simulated time/frames at which they were reported.

2D committor solver

by relaxation (finite-difference discretization) with coarse-to-fine refinement.

Track and plot paths lineage

worker.log file produced by AIMMD for each shooting chain.

graphics.

Dependencies and expectations

  • PathEnsemble conventions: pathensemble.types() returns per-path type strings and pathensemble[...] supports boolean mask selection.

  • Type-pattern matching is delegated to aimmd.pathensemble.utils.match_patterns().

Notes

The bin construction functions assume that the projected reaction-coordinate values order terminal states consistently (as in AIMMD training). Bounds are made robust by explicit cutoffs and by falling back from free-excursion statistics to transition statistics when required.

aimmd.analysis.utils.find_extremes_with_free_simulations(pathensemble, states='ARB', source='values')[source]

Estimate bin extremes using free excursions.

This helper is used by compute_bins() to obtain left/right limits for a bin range using only free excursions originating from each terminal state.

It extracts two subsets:

  • free excursions originating from a and returning to a (pattern f'{a}{r}.{a}')

  • free excursions originating from b and returning to b (pattern f'{b}{r}.{b}')

From these, it computes:

  • e1: sorted per-path maxima (descending) from a-side free excursions

  • e2: sorted per-path minima (ascending) from b-side free excursions

The index used is n_transitions + 1, where n_transitions counts paths matching either direction of the requested transition.

Parameters:

pathensemble – Ensemble supporting:

  • types() returning per-path type strings

  • boolean-mask selection pathensemble[mask]

  • per-path extrema selectors max and min.

Parameters:

states – Triplet of state labels (a, r, b):

  • a: left terminal state

  • r: reactive region label

  • b: right terminal state

Parameters:

source (str, default='values') – Name of the per-frame stream used when computing per-path extrema.

Returns:

(begin, end) suggested extremes for binning.

Return type:

(float, float)

Notes

If no free excursions are present on a side, the corresponding list is replaced with [-inf] or [+inf] to keep the selection logic stable.

aimmd.analysis.utils.find_extremes_with_transitions(pathensemble, states='ARB', source='values')[source]

Estimate bin extremes using transition trajectories.

This helper is used by compute_bins() when free excursions do not give reasonable bounds (or when explicitly requested via find_extremes_with).

It identifies transition paths in either direction (A→B or B→A) using type-pattern matching and then collects representative boundary-adjacent values from those paths:

  • For paths with type prefix equal to states (A→R→B):

  • begin uses _position(1, source) (near the start, excluding boundary)

  • end uses _position(-2, source) (near the end, excluding boundary)

  • For the reverse direction:

  • begin/end are swapped accordingly.

The returned extremes are the medians over the collected values.

Parameters:

pathensemble – Ensemble supporting:

  • types() returning per-path type strings

  • integer indexing returning Path objects

  • Path method _position(i, source).

Parameters:

  • states (str, default='ARB') – Triplet (a, r, b) defining which transitions are considered.

  • source (str, default='values') – Stream name passed to _position.

Returns:

(begin, end) suggested extremes for binning.

Return type:

(float, float)

Notes

If no transition trajectories are present, returns (0., 0.).

aimmd.analysis.utils.compute_bins(pathensemble, nbins, cutoff_max=20.0, cutoff_min=0.5, find_extremes_with='free', source='values', states='ARB', terminal_bin_extension='all')[source]

Construct 1D bin boundaries for projection/analysis.

This function returns a 1D array of bin boundaries of length nbins + 1. The finite (non plus or minus inf) part of the grid spans the interval [begin, end], where (begin, end) are inferred from the ensemble using one of the extreme estimators:

The returned boundaries are then optionally given marginal outer bins by replacing the first and/or last boundary with -inf and/or +inf. The finite interior boundaries are still uniformly spaced between begin and end.

Parameters:

  • pathensemble (PathEnsemble-like) – Ensemble used to infer the finite interval.

  • nbins – Number of bins (including marginal bins if requested).

The returned array has length nbins + 1. If nbins <= 0, returns an empty array. :type nbins: int :param cutoff_max: Hard clip for the finite interval. If the free-based extremes exceed this range, transition-based extremes are used and clipped. :type cutoff_max: float, default=20. :param cutoff_min: Inner clip to keep the finite interval away from zero: begin <= -cutoff_min and end >= +cutoff_min when possible. :type cutoff_min: float, default=0.5 :param find_extremes_with: Preferred heuristic for obtaining (begin, end). Even when ‘free’ is requested, transition-based extremes are used if the free-based bounds fall outside the configured cutoffs. :type find_extremes_with: {‘free’, ‘transitions’}, default=’free’ :param source: Passed to the extreme estimators. :type source: str, default=’values’ :param states: Triplet (a, r, b) used by the extreme estimators. :type states: str, default=’ARB’ :param terminal_bin_extension: Which marginal bins to include by replacing outer boundaries:

  • ‘all’ : replace both ends (first boundary = -inf, last boundary = +inf)

  • string containing a : replace the left boundary with -inf

  • string containing b : replace the right boundary with +inf

The replacement changes the outermost bin to be unbounded. The finite interior bins still span [begin, end]. :type terminal_bin_extension: {‘all’} or str, default=’all’

Returns:

1D array of bin boundaries. If marginal bins are enabled, the first and/or

last element is ±inf. Otherwise all boundaries are finite. :rtype: numpy.ndarray

Notes

Implementation detail:

  • The number of finite boundaries generated by np.linspace is nbins + 1 - left - right where left/right indicate marginal bins.

  • If marginal bins are requested, -inf and/or +inf are then inserted as the first/last boundary, replacing the corresponding finite boundary.

aimmd.analysis.utils.bin_centers(bins)[source]

Return bin centers, supporting ±inf marginal bins.

Parameters:

bins (array-like) – 1D array of bin boundaries of length >= 2. Bin i spans [bins[i], bins[i+1]).

Returns:

Array of length len(bins) - 1 with one center per bin.

Return type:

numpy.ndarray

Notes

For bins with infinite boundaries, finite extrapolated centers are produced using linear extrapolation from the nearest finite bins.

aimmd.analysis.utils.merge_empty_bins(bins, keepers, *histograms, center=0)[source]

Merge low-occupancy bins into the closest occupied bins moving away from the reference center. Any additional histograms are merged accordingly by summation.

An empty bin is merged only if there is an occupied bin further outward on the same side of center. Thus, bins left of or at center merge leftward, while bins right of center merge rightward. Empty edge bins with no occupied bin further outward are left unchanged.

Parameters:

bins – One-dimensional array of bin boundaries with shape (k + 1,).

Bin i spans [bins[i], bins[i + 1]). :type bins: np.ndarray :param keepers: Indices or boolean mask identifying bins that are occupied and therefore kept as merge targets. :type keepers: np.ndarray :param *histograms: One or more one-dimensional histograms of shape (len(bins) - 1,) defined on the original bins. Each histogram is merged to match the returned merged_bins by summing the values of all original bins contributing to each merged bin. If one of the histograms is None, just leave unchanged. :type *histograms: np.ndarray :param center: Reference value defining the outward merge direction. Bins whose centers are less than or equal to center merge toward smaller values; bins whose centers are greater than center merge toward larger values. :type center: float, default=0

Returns:

  • merged_bins (np.ndarray) – Updated bin boundaries after merging.

  • merged_bin_counts (np.ndarray) – Number of original bins contributing to each merged bin.

  • *merged_histograms (np.ndarray) – The input histograms after merging by summation, returned in the

same order as provided.

aimmd.analysis.utils.merge_marginal_bins(bins, *values, min_values=3)[source]

Merge low-occupancy marginal bins.

This helper is used to stabilize analysis when the outermost bins are too sparse. It merges bins from the left edge and from the right edge such that:

  1. all bins from the second to the second-to-last are a contiguous subset of the original bins;

  2. the above bins contain at least min_values counts for all the datasets provided;

  3. given the above conditions, the number of bins is maximum.

The occupancy criterion is computed from the per-dataset histograms: np.histogram(v, bins=bins)[0] and then taking the minimum over datasets.

Parameters:

  • bins (np.ndarray) – 1D array of bin boundaries of shape (k+1,). Bin i spans [bins[i], bins[i+1]).

  • *values – One or more 1D arrays of raw samples. Each provided dataset contributes

its histogram to the “minimum occupancy” decision. :type *values: np.ndarray :param min_values: Minimum required histogram count in each marginal bin. :type min_values: int, default=3

Returns:

  • merged_bins (np.ndarray) – Updated bin boundaries after merging.

  • merged_bin_counts (np.ndarray) – Number of original bins merged into each new bin. This can be used to

scale selections consistently when downstream code depends on the original binning density.

Notes

The function always preserves the first and last boundary of bins.

aimmd.analysis.utils.binomial_mean_and_confidence_interval(r1, r2, alpha=0.95)[source]

Binomial mean and two-sided confidence interval.

Parameters:

r1 – Two outcome counts. The code interprets:

  • n = r1 + r2

  • k = r2

and computes the confidence interval for p = k/n. :type r1: int :param r2: Two outcome counts. The code interprets:

  • n = r1 + r2

  • k = r2

and computes the confidence interval for p = k/n. :type r2: int :param alpha: Confidence level. :type alpha: float, default=0.95

Returns:

(p, lower, upper), where p = k/n and (lower, upper) is a two-sided

beta-based confidence interval. :rtype: (float, float, float)

Notes

This uses the standard Beta quantile construction:

  • lower = Beta(a/2; k, n-k+1)

  • upper = Beta(1-a/2; k+1, n-k)

with special-case handling for k=0 and k=n.

aimmd.analysis.utils.extract_rate_estimates_from_log_file(fname)[source]

Parse rate estimates from an AIMMD training log.

This helper reads AIMMD train*.log files produced during training and extracts the time series of rate-constant estimates printed in the log.

Parsed quantities The function returns three arrays:

  • t : cumulative simulated time (as printed by the logger)

  • k12 : rate constant estimate for the state 1 to state 2 transition

  • k21 : rate constant estimate for the state 2 to state 1 transition

Parameters:

fname (str or path-like) – Path to a training log file (e.g., trainARB.log).

Returns:

(t, k12, k21) as NumPy arrays of dtype float.

Return type:

(np.ndarray, np.ndarray, np.ndarray)

Notes

This parser assumes the log format where:

  • a line containing 'k12 estimate' holds the k12 value after 'estimate:'

  • the next relevant line holds the k21 value in the same format

  • the subsequent relevant line begins with the number of frames followed by the token 'frames' (used as the time coordinate)

The parsing logic is intentionally minimal and matches the existing logger output exactly.

aimmd.analysis.utils.extract_bias_reweighted_rate_estimates_from_log_file(fname)[source]

Parse bias-reweighted rate estimates from an AIMMD training log.

Complements extract_rate_estimates_from_log_file() for biased runs (record_bias=True). The bias-reweighted lines appear after the raw estimate block in each training round:

k12 estimate: 1.234e-03 [1/dt] k21 estimate: 5.678e-04 [1/dt] 12500 frames (excluded margins) Bias check passed … k12 bias-reweighted: 1.200e-03 [1/dt] k21 bias-reweighted: 5.500e-04 [1/dt]

The time coordinate t is taken from the frames line (same as for the raw estimates), so the returned arrays align index-by-index with those returned by extract_rate_estimates_from_log_file().

Parameters:

fname (str or path-like) – Path to a training log file (e.g. trainARB.log).

Returns:

(t, k12_rw, k21_rw) as NumPy arrays. Empty arrays if the log

contains no bias-reweighted lines (e.g. record_bias=False run). :rtype: (np.ndarray, np.ndarray, np.ndarray)

aimmd.analysis.utils.solve_committor_by_relaxation(X, Y, Fx, Fy, A, B, P0, progress=[5, 4, 2, 1])[source]

Compute committor in 2D with a relaxation method.

This routine solves the steady-state committor equation on a 2D grid using a finite-difference discretization and iterative relaxation.

The method enforces:

  • P=0 on region A

  • P=1 on region B

  • reflecting (zero-gradient) boundary conditions at the grid edges

Parameters:

  • X (np.ndarray) – 2D arrays defining the grid coordinates. Shapes must match.

  • Y (np.ndarray) – 2D arrays defining the grid coordinates. Shapes must match.

  • Fx (np.ndarray) – 2D arrays of drift/force components on the same grid.

  • Fy (np.ndarray) – 2D arrays of drift/force components on the same grid.

  • A (np.ndarray of bool) – Boolean masks on the grid indicating the A and B regions.

  • B (np.ndarray of bool) – Boolean masks on the grid indicating the A and B regions.

  • P0 – Initial guess for the committor field. Updated in-place by coarse-to-fine

passes. :type P0: np.ndarray :param progress: Coarsening schedule. For each split, the solver runs on the subgrid [::split, ::split] and, if split>1, interpolates back to the full grid before continuing at higher resolution. :type progress: list[int], default=[5, 4, 2, 1]

Returns:

Final committor estimate on the full grid.

Return type:

np.ndarray

Notes

  • This function uses tqdm(progress) but does not import tqdm. The caller must ensure that tqdm is available in scope.

  • The update rule clamps P to [0,1] each iteration and enforces boundary conditions after each sweep.

  • When refining, interpolation uses scipy.interpolate.RegularGridInterpolator.

aimmd.analysis.utils.find_path_lineages(*shooting_chains, verbose=False)[source]

Reconstructs path lineages by parsing the specific worker.log for each chain.

This function attaches a ._previous attribute to path objects. It dynamically identifies the correct log file for every path by looking at its directory.

Parameters:

  • *shooting_chains (list of PathChain objects) – One or more chains containing path objects to be linked.

  • verbose (bool, default=False) – If True, displays a progress bar for the parsing process.

aimmd.analysis.utils.plot_path_lineages(*shooting_chains, out='path_tree.pdf', states='ARB', source='values', vmin=-20, vmax=20, fields=None, show=False)[source]

Visualizes the genealogical history of AIMMD path sampling results ensembles using an independently-packed grid layout.

Each lineage is treated as an autonomous vertical unit, allowing columns to stack tightly without being constrained by the length of chains in adjacent tiles. Labels are centered relative to the data axis of each column.