Fix circular-coordinate handling and histogram support in NMI

mdpath/mdpath.py

@@ -21,7 +21,7 @@
 import pickle
 
 from mdpath.src.structure import StructureCalculations, DihedralAngles
-from mdpath.src.mutual_information import NMICalculator
+from mdpath.src.mutual_information import NMICalculator, binning_diagnostic
 from mdpath.src.graph import GraphBuilder
 from mdpath.src.cluster import PatwayClustering
 from mdpath.src.visualization import MDPathVisualize
@@ -216,6 +216,26 @@ def main():
     parser.add_argument("-null-seed", type=int, default=0)
     parser.add_argument("-ddof", type=int, choices=(0, 1), default=0)
     parser.add_argument("-null-outdir", default="null_analysis")
+    parser.add_argument(
+        "-hist-range",
+        dest="hist_range",
+        nargs="+",
+        default=["-180", "180"],
+        help="Fixed histogram support for the dihedral-movement histograms, in "
+        "degrees, as two floats (default: -180 180). Fixing the bin edges is what "
+        "makes NMI comparable across replicas, blocks and bootstrap samples. Pass "
+        '"none" to restore the pre-fix behaviour, where each input derived its own '
+        "support from its own extrema.",
+    )
+    parser.add_argument(
+        "-diagnose-bins",
+        dest="diagnose_bins",
+        action="store_true",
+        help="Write binning_diagnostic.csv: per-residue histogram occupancy under "
+        "the chosen support and branch-cut exposure (frames with |dphi| > 180 "
+        "before wrapping). Use it to check whether the fixed support is too wide "
+        "for the observed dihedral movements.",
+    )
     parser.add_argument(
         "-experimental-c-corr",
         dest="experimental_c_corr",
@@ -225,6 +245,15 @@ def main():
     )
 
     args = parser.parse_args()
+    # -hist-range: two floats, or the literal "none" to restore data-derived binning.
+    if len(args.hist_range) == 1 and str(args.hist_range[0]).lower() == "none":
+        hist_range = None
+    elif len(args.hist_range) == 2:
+        hist_range = (float(args.hist_range[0]), float(args.hist_range[1]))
+        if hist_range[0] >= hist_range[1]:
+            raise ValueError("-hist-range needs an increasing (low, high) pair")
+    else:
+        raise ValueError('-hist-range takes two floats or the literal "none"')
     if not args.topology or not args.trajectory:
         print("Both trajectory and topology files are required!")
         exit()
@@ -333,8 +362,16 @@ def main():
             )
     print("\033[1mTrajectory is processed and ready for analysis.\033[0m")
 
+    if args.diagnose_bins:
+        binning_diagnostic(
+            df_all_residues,
+            num_bins=35,
+            hist_range=hist_range,
+            csv_path="binning_diagnostic.csv",
+        )
+
     # Calculate the mutual information and build the graph
-    nmi_calc = NMICalculator(df_all_residues, invert=invert)
+    nmi_calc = NMICalculator(df_all_residues, invert=invert, hist_range=hist_range)
     nmi_calc.entropy_df.to_csv("entropy_df.csv", index=False)
     nmi_calc.nmi_df.to_csv("nmi_df.csv", index=False)
     graph_builder = GraphBuilder(
@@ -354,6 +391,7 @@ def main():
             per_replica_dfs,
             pairs=graph_pairs,
             ddof=args.ddof,
+            hist_range=hist_range,
         )
         edge_conf_df = edge_conf_calc.build_edge_confidence_df()
         edge_conf_df_graph = EdgeConfidenceCalculator.filter_to_graph_edges(
@@ -401,6 +439,7 @@ def main():
                 seed=args.null_seed,
                 ddof=args.ddof,
                 n_jobs=num_parallel_processes,
+                hist_range=hist_range,
             )
             null_outputs = null_calc.run(args.null_outdir, args.null_block_ns)
             status = null_outputs["edge_null_confidence"].corr_status.value_counts(
@@ -471,6 +510,7 @@ def main():
             topology,
             structure_calc.last_res_num,
             graphdist,
+            hist_range=hist_range,
         )
         file_name = "path_confidence_intervals.txt"
         bootstrap_analysis.bootstrap_write(file_name)

mdpath/src/bootstrap.py

@@ -37,6 +37,8 @@ class BootstrapAnalysis:
 
         num_bins (int): Number of bins to group dihedral angle movements into for NMI calculation. Defaults to 35.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` histogram support forwarded to the NMI calculation. Each resample otherwise derives its own support, so path-confidence values inherit that inconsistency. Defaults to ``(-180.0, 180.0)``.
+
         common_counts (np.ndarray): Array with the counts of common paths between the original sample and bootstrap samples.
 
         path_confidence_intervals (dict): Dictionary with the confidence intervals for each path.
@@ -53,6 +55,7 @@ def __init__(
         last_residue: int,
         graphdist: int,
         num_bins: int = 35,
+        hist_range=(-180.0, 180.0),
     ) -> None:
         self.df_all_residues = df_all_residues
         self.df_distant_residues = df_distant_residues
@@ -63,6 +66,7 @@ def __init__(
         self.last_residue = last_residue
         self.graphdist = graphdist
         self.num_bins = num_bins
+        self.hist_range = hist_range
         self.common_counts, self.path_confidence_intervals = self.bootstrap_analysis()
 
     def create_bootstrap_sample(self, df_dihedral: pd.DataFrame) -> tuple:
@@ -100,7 +104,9 @@ def process_bootstrap_sample(
             bootstrap_pathways (list): List of paths within the bootstrap sample.
         """
         bootstrap_sample = self.create_bootstrap_sample(self.df_all_residues)
-        nmi_calculator = NMICalculator(bootstrap_sample, num_bins=self.num_bins)
+        nmi_calculator = NMICalculator(
+            bootstrap_sample, num_bins=self.num_bins, hist_range=self.hist_range
+        )
         bootstrap_mi_diff = nmi_calculator.nmi_df
         graph = GraphBuilder(
             self.pdb, self.last_residue, bootstrap_mi_diff, self.graphdist

mdpath/src/confidence.py

@@ -11,8 +11,18 @@
 These scores measure *consistency of the estimated edge magnitude* across
 replicas; they are not calibrated probabilities, significance thresholds, or
 statements about effect size (an edge with consistently negligible mutual
-information can still score highly). The finite-sampling block reference and its
-experimental corrected score live in :mod:`mdpath.src.null_model`.
+information can still score highly).
+
+.. important::
+
+   Per-replica values are only comparable to each other under a **fixed histogram
+   support**. With ``hist_range=None`` every replica derives its own bin edges from
+   its own extrema, so a single wider excursion in one replica changes that
+   replica's binning and the difference propagates directly into ``MI_std`` ->
+   ``MI_cv`` -> ``confidence``. The default fixed support exists for this reason.
+
+The finite-sampling block reference and its experimental corrected score live in
+:mod:`mdpath.src.null_model`.
 
 Classes
 --------
@@ -43,6 +53,12 @@ class EdgeConfidenceCalculator:
 
         num_bins (int): Number of histogram bins forwarded to the mutual information calculator.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` histogram support forwarded to the
+            mutual information calculator. Per-replica values are only comparable under a
+            fixed support: with ``None`` each replica derives its own bin edges from its own
+            extrema, and that difference flows straight into ``MI_std`` -> ``MI_cv`` ->
+            ``confidence``. ``None`` is retained only to reproduce pre-fix numbers.
+
         pairs: Optional restriction of residue pairs forwarded to the mutual information calculator; ``None`` evaluates all pairs.
 
         ddof (int): Delta degrees of freedom (0 or 1) used for the standard deviation and the finite-replica floor.
@@ -56,6 +72,7 @@ def __init__(
         num_bins: int = 35,
         pairs=None,
         ddof: int = 0,
+        hist_range=(-180.0, 180.0),
     ) -> None:
         if len(per_replica_dfs) < 2:
             raise ValueError("EdgeConfidenceCalculator requires at least two replicas.")
@@ -64,6 +81,7 @@ def __init__(
         self.per_replica_dfs = per_replica_dfs
         self.num_bins = num_bins
         self.pairs, self.ddof = pairs, ddof
+        self.hist_range = hist_range
         self.per_replica_nmi_dfs = self._compute_per_replica_nmi()
 
     def _compute_per_replica_nmi(self) -> List[pd.DataFrame]:
@@ -78,6 +96,7 @@ def _compute_per_replica_nmi(self) -> List[pd.DataFrame]:
                 num_bins=self.num_bins,
                 pairs=self.pairs,
                 show_progress=False,
+                hist_range=self.hist_range,
             ).nmi_df
             for frame in self.per_replica_dfs
         ]

mdpath/src/mutual_information.py

@@ -9,6 +9,11 @@
 --------
 
 :class:`NMICalculator`
+
+Functions
+---------
+
+:func:`binning_diagnostic`
 """
 
 import pandas as pd
@@ -29,6 +34,14 @@ class NMICalculator:
 
         num_bins (int): Number of bins to use for histogram calculations. Default is 35.
 
+        hist_range (tuple | None): Fixed ``(low, high)`` support for the
+            dihedral-movement histograms, in degrees. Default ``(-180, 180)``
+            fixes the bin edges so that values computed from different
+            replicas, blocks or bootstrap samples share one binning and are
+            directly comparable. Passing ``None`` restores the previous
+            behaviour, in which each input derived its own support from its
+            own extrema.
+
         pairs (list | None): Optional list of residue pairs restricting which pairs
             are evaluated (each pair is any two residue labels, e.g. ``42`` or
             ``"Res 42"``). When ``None`` every pair is computed as before. This only
@@ -48,12 +61,14 @@ def __init__(
         invert=False,
         pairs=None,
         show_progress: bool = True,
+        hist_range=(-180.0, 180.0),
     ) -> None:
         self.df_all_residues = df_all_residues
         self.num_bins = num_bins
         self.invert = invert
         self.pairs = pairs
         self.show_progress = show_progress
+        self.hist_range = hist_range
         self.nmi_df, self.entropy_df = self.NMI_calcs()
 
     def _column_pairs(self, columns):
@@ -88,7 +103,11 @@ def NMI_calcs(self):
         histograms = {}
         entropys = {}
         for col in columns:
-            hist, _ = np.histogram(self.df_all_residues[col], bins=self.num_bins)
+            hist, _ = np.histogram(
+                self.df_all_residues[col],
+                bins=self.num_bins,
+                range=self.hist_range,
+            )
             histograms[col] = hist
             entropys[col] = entropy(hist)
 
@@ -103,6 +122,14 @@ def NMI_calcs(self):
                 self.df_all_residues[col1],
                 self.df_all_residues[col2],
                 bins=self.num_bins,
+                # np.histogram wants (lo, hi); np.histogram2d wants
+                # [(lo, hi), (lo, hi)] and raises otherwise. Both accept None,
+                # so the None branch is a true no-op.
+                range=(
+                    None
+                    if self.hist_range is None
+                    else [self.hist_range, self.hist_range]
+                ),
             )
             mi = mutual_info_score(
                 histograms[col1], histograms[col2], contingency=hist_joint
@@ -121,3 +148,77 @@ def NMI_calcs(self):
             max_nmi_diff = nmi_df["MI Difference"].max()
             nmi_df["MI Difference"] = max_nmi_diff - nmi_df["MI Difference"]
         return nmi_df, entropy_df
+
+
+def binning_diagnostic(
+    df_all_residues: pd.DataFrame,
+    num_bins: int = 35,
+    hist_range=(-180.0, 180.0),
+    csv_path: str = "binning_diagnostic.csv",
+) -> pd.DataFrame:
+    """Report histogram occupancy and branch-cut exposure per residue.
+
+    A fixed support makes mutual information comparable across replicas, blocks and
+    bootstrap samples, but it costs resolution when the dihedral movements are narrow:
+    bins outside the sampled range sit empty. This quantifies that trade-off so the
+    width of the support can be chosen from data rather than guessed.
+
+    Args:
+        df_all_residues (pd.DataFrame): Dihedral angle movements, one column per residue.
+
+        num_bins (int): Number of histogram bins. Default is 35.
+
+        hist_range (tuple | None): Fixed ``(low, high)`` support to evaluate. ``None``
+            reports the data-derived binning instead.
+
+        csv_path (str): Where the per-residue table is written.
+
+    Returns:
+        pd.DataFrame: One row per residue with ``residue``, ``occupied_bins``,
+        ``occupied_fraction``, ``movement_min``, ``movement_max`` and, when the
+        upstream calculation recorded it, ``branch_cut_crossings`` (the number of
+        frames whose *unwrapped* difference exceeded 180 degrees).
+    """
+    crossings = df_all_residues.attrs.get("branch_cut_crossings", {})
+    records = []
+    for column in df_all_residues.columns:
+        values = np.asarray(df_all_residues[column], dtype=float)
+        values = values[np.isfinite(values)]
+        counts, _ = np.histogram(values, bins=num_bins, range=hist_range)
+        record = {
+            "residue": column,
+            "occupied_bins": int(np.count_nonzero(counts)),
+            "occupied_fraction": float(np.count_nonzero(counts) / num_bins),
+            "movement_min": float(values.min()) if values.size else np.nan,
+            "movement_max": float(values.max()) if values.size else np.nan,
+        }
+        try:
+            record["branch_cut_crossings"] = int(crossings[residue_id(column)])
+        except (KeyError, ValueError):
+            record["branch_cut_crossings"] = np.nan
+        records.append(record)
+
+    frame = pd.DataFrame(records)
+    if csv_path:
+        frame.to_csv(csv_path, index=False)
+
+    occupied = frame["occupied_bins"].to_numpy()
+    low, high = np.percentile(occupied, [5, 95]) if occupied.size else (np.nan, np.nan)
+    sparse = int(np.count_nonzero(occupied < 10))
+    print(
+        f"Binning diagnostic ({num_bins} bins, support "
+        f"{'data-derived' if hist_range is None else f'{hist_range[0]:g}..{hist_range[1]:g}'}): "
+        f"occupied bins median {np.median(occupied):.0f}, 5th/95th pct {low:.0f}/{high:.0f}; "
+        f"{sparse}/{len(occupied)} residues ({sparse / max(1, len(occupied)):.1%}) occupy <10 bins."
+    )
+    exposure = frame["branch_cut_crossings"]
+    if exposure.notna().any():
+        affected = int((exposure.fillna(0) > 0).sum())
+        print(
+            f"  branch-cut exposure: {affected}/{len(exposure)} residues "
+            f"({affected / max(1, len(exposure)):.1%}) had at least one frame with "
+            f"|dphi| > 180 before wrapping; {int(exposure.fillna(0).sum())} frames total."
+        )
+    if csv_path:
+        print(f"  per-residue values written to {csv_path}")
+    return frame