from __future__ import annotations
import contextlib
import io
import json
import re
from datetime import datetime
from pathlib import Path
from typing import Any, Dict, List
import numpy as np
try:
import msgspec
except ImportError: # pragma: no cover - exercised by minimal installations
msgspec = None
from pycurves_lib.data.modified_bases import parent_base_name
from pycurves_lib.topology.base_annotations import (
annotate_context,
base_pair_observed_geometry_tag,
render_section_m,
)
from pycurves_lib.io.curves_visualization_payload import VisualizationPayloadMixin
def _to_jsonable(value: Any):
if isinstance(value, np.ndarray):
return _to_jsonable(value.tolist())
if isinstance(value, np.generic):
return _to_jsonable(value.item())
if isinstance(value, float):
return value if np.isfinite(value) else None
if isinstance(value, dict):
return {str(k): _to_jsonable(v) for k, v in value.items()}
if isinstance(value, (list, tuple)):
return [_to_jsonable(v) for v in value]
return value
def _json_encode_hook(value: Any):
if isinstance(value, np.ndarray):
return value.tolist()
if isinstance(value, np.generic):
return value.item()
raise TypeError(f"Unsupported JSON value: {type(value).__name__}")
def _json_dumps(value: Any) -> str:
"""Return compact JSON, using msgspec directly when it is installed."""
if msgspec is not None:
return msgspec.json.encode(value, enc_hook=_json_encode_hook).decode("utf-8") + "\n"
return json.dumps(
_to_jsonable(value),
separators=(",", ":"),
allow_nan=False,
) + "\n"
def _all_finite(values) -> bool:
return bool(np.all(np.isfinite(np.asarray(values, dtype=float))))
STEP_PARAMETERS = ("shift", "slide", "rise", "tilt", "roll", "twist")
AXIS_PARAMETERS = ("xdisp", "ydisp", "inclin", "tip")
BASE_BASE_PARAMETERS = ("shear", "stretch", "stagger", "buckle", "propel", "opening")
SUGAR_PUCKERS = (
"C3'-endo", "C4'-exo", "O1'-endo", "C1'-exo", "C2'-endo",
"C3'-exo", "C4'-endo", "O1'-exo", "C1'-endo", "C2'-exo",
)
def _parameter_record(values, parameter_names, **metadata):
"""Build one structured parameter row, or None if the values are invalid."""
if values is None:
return None
values = np.asarray(values, dtype=float)
if not _all_finite(values[:len(parameter_names)]):
return None
row = dict(metadata)
for idx, name in enumerate(parameter_names):
row[name] = values[idx]
return row
def _nullable_parameter_record(values, parameter_names, **metadata):
"""Build a structured row, preserving invalid/missing values as JSON null."""
row = dict(metadata)
if values is None:
values = [None] * len(parameter_names)
else:
values = np.asarray(values, dtype=float)
for idx, name in enumerate(parameter_names):
try:
value = float(values[idx])
except (TypeError, ValueError, IndexError):
value = np.nan
row[name] = value if np.isfinite(value) else None
return row
def _json_number(value):
"""Return a finite JSON number, otherwise None."""
try:
number = float(value)
except (TypeError, ValueError):
return None
return number if np.isfinite(number) else None
def _backbone_number(value):
"""Return a finite backbone value, otherwise None."""
return _json_number(value)
class CurvesOutputFormatter(VisualizationPayloadMixin):
"""Render pyCurves results as Curves-style text or structured JSON/Pandas."""
def __init__(
self,
runner,
visualization: bool = False,
fast_groove: bool = False,
):
self.runner = runner
self.include_visualization = visualization
self.fast_groove = bool(fast_groove and not visualization)
self._annotation_cache = None
def get_dataframes(self):
"""Return a dictionary of Pandas DataFrames for all computed parameters."""
try:
import pandas as pd
except ImportError as exc:
raise ImportError(
"Pandas or one of its dependencies is not installed. "
"Run 'pip install pandas python-dateutil' to use get_dataframes()."
) from exc
self._require_results()
dataframes = self._build_dataframes()
dfs = {}
for k, v in dataframes.items():
if k == "groove" and isinstance(v, dict) and "data" in v:
flat_list = []
for lvl, l_data in v["data"].items():
for sub_lvl, sl_data in l_data["sub_levels"].items():
row = {
"atom_defining_backbone": v.get("atom_defining_backbone", ""),
"total_levels": v.get("levels", 0),
"total_sub_levels": v.get("sub_levels", 0),
"level": int(lvl),
"base_pair": l_data.get("base_pair", ""),
"sub_level": int(sub_lvl),
}
row.update(sl_data)
flat_list.append(row)
dfs[k] = pd.DataFrame(flat_list)
else:
dfs[k] = pd.DataFrame(v)
return dfs
[docs]
def render_curves_text(self, header_name: str = "pyCurves") -> str:
self._require_results()
annotations = self._annotations()
outaxe_text = self._capture(self.runner.calc.outaxe).rstrip()
records = self._build_dataframes()
pieces = [
self._header(header_name=header_name),
self.runner.analysis_log.rstrip(),
outaxe_text,
self._render_section_l_local_base_base(records).rstrip(),
]
pieces.append(render_section_m(annotations).rstrip())
return "\n\n".join(piece for piece in pieces if piece) + "\n"
def _build_json_payload(self) -> Dict[str, Any]:
self._require_results()
dataframes = self._build_dataframes()
payload = {
"program": "pyCurves",
"format": "pycurves-slim-v2",
"frame_convention": self._frame_convention_payload(),
"analysis_options": self._analysis_options_payload(),
"inputs": {
"inpfile": self.runner.inpfile,
"pdbfile": self.runner.pdbfile,
"generated_inpfiles": self.runner.generated_inpfiles,
},
"sequence": self._sequence_payload(),
"dataframes": dataframes,
}
if self.include_visualization:
payload["rebuild_frame_mappings"] = self._rebuild_frame_mappings_payload()
payload["visualization"] = self._visualization_payload(dataframes)
return payload
def _rebuild_frame_mappings_payload(self) -> Dict[str, Any]:
"""Return relative interaction/standard frame maps for exact rebuilding."""
ctx = self.runner.ctx
calc = self.runner.calc
standard_frames = np.asarray(ctx.params.frames, dtype=float)
interaction_frames = np.asarray(
getattr(ctx.params, "shape_frames", standard_frames),
dtype=float,
)
annotations = {}
for annotation in self._base_pair_rows(self._annotations()):
level = annotation.get("level")
strand_1 = int(annotation.get("strand_1", 1))
strand_2 = int(annotation.get("strand_2", 2))
if level is None or 1 not in (strand_1, strand_2):
continue
partner = (strand_2 if strand_1 == 1 else strand_1) - 1
key = (partner, int(level))
if key not in annotations or annotation.get("reference_pair", True):
annotations[key] = annotation
variants = getattr(ctx, "step_pair_frame_variants", {}) or {}
reverse_z = getattr(ctx, "local_step_reverse_z", {}) or {}
normal_signs = getattr(ctx, "pair_normal_signs", {}) or {}
normal_modes = getattr(ctx, "pair_normal_branch_modes", {}) or {}
levels = []
for partner in range(1, ctx.nst):
for level in range(1, ctx.nux + 1):
if not (calc._has_level(0, level) and calc._has_level(partner, level)):
continue
annotation = annotations.get((partner, level), {})
geometry = annotation.get("contact_geometry") or {}
lw_family = base_pair_observed_geometry_tag(annotation) or "cWW"
members = []
for strand, role in ((0, "primary"), (partner, "partner")):
standard = standard_frames[strand, level]
interaction = interaction_frames[strand, level]
if not (_all_finite(standard) and _all_finite(interaction)):
continue
rotation = interaction[:3] @ standard[:3].T
translation = interaction[:3] @ (standard[3] - interaction[3])
annotation_side = next(
(
side for side in (1, 2)
if int(annotation.get(f"strand_{side}", side)) == strand + 1
),
1 if strand == 0 else 2,
)
edge = str(
geometry.get(f"edge_{annotation_side}")
or annotation.get(f"edge_{annotation_side}")
or (lw_family[annotation_side] if len(lw_family) > annotation_side else "W")
).upper()
members.append({
"strand": strand + 1,
"role": role,
"edge": edge,
"interaction_from_standard": {
"rotation": rotation,
"translation": translation,
},
})
key = (0, partner, level)
row = {
"level": level,
"partner_strand": partner + 1,
"lw_family": lw_family,
"frame_mode": str(
geometry.get("frame_mode")
or annotation.get("frame_mode")
or "fitted_base"
),
"pair_frame_branch": {
"step_from_interaction_axis_signs": list(
variants.get(key, (1, 1, 1))
),
"normal_sign": int(normal_signs.get(key, 1)),
"normal_branch_mode": str(normal_modes.get(key, "fitted_base")),
},
"members": members,
}
next_level = level + 1
step_key = (0, partner, level, next_level)
if step_key in reverse_z:
row["outgoing_step_branch"] = {
"next_level": next_level,
"reverse_z": bool(reverse_z[step_key]),
}
levels.append(row)
return {
"schema": "pycurves-rebuild-frame-mappings-v1",
"scope": "standard-base-and-c1-prime frames",
"reference_frame_convention": getattr(ctx.cfg, "frame_convention", "standard"),
"transform_definition": {
"name": "interaction_from_standard",
"equation": "x_interaction = rotation @ x_standard + translation",
"frame_storage": "global = axes.T @ local + origin",
"rotation_layout": "row-major",
"translation_unit": "angstrom",
},
"pair_branch_definition": {
"equation": "step_axes = diag(axis_signs) @ interaction_pair_axes",
"reverse_z_correction": (
"Before rebuilding, negate Rise and Twist when "
"outgoing_step_branch.reverse_z is true"
),
},
"levels": levels,
}
def _frame_convention_payload(self) -> Dict[str, Any]:
cfg = self.runner.ctx.cfg
name = getattr(cfg, "frame_convention", "standard")
library = getattr(self.runner.ctx, "base_reference_library", None)
return {
"name": name,
"compatible_with": ["Curves+", "3DNA", "x3dna"] if name == "standard" else ["Curves 5.3"],
"reference_source": getattr(library, "source", "setup.f hardcoded bref"),
"axis_convention": getattr(cfg, "axis_convention", "global"),
}
def _analysis_options_payload(self) -> Dict[str, Any]:
cfg = self.runner.ctx.cfg
return {
"continuous_strands": bool(getattr(self.runner, "continuous_strands", False)),
"duplex_only": bool(getattr(self.runner, "duplex_only", False)),
"comb": bool(getattr(cfg, "comb", False)),
"groove": bool(getattr(cfg, "grv", False)),
"ends": bool(getattr(cfg, "ends", False)),
"mini": bool(getattr(cfg, "mini", False)),
"axis_convention": getattr(cfg, "axis_convention", "global"),
"axis_weighting": bool(getattr(cfg, "axis_weighting", False)),
"altloc": getattr(self.runner.ctx.molecule, "altloc_selection", "first"),
"available_altlocs": list(getattr(self.runner.ctx.molecule, "available_altlocs", ())),
}
def _sequence_payload(self) -> Dict[str, Any]:
ctx = self.runner.ctx
strands = []
for strand in range(ctx.nst):
residues = []
sequence = []
for level in range(1, ctx.nux + 1):
if ctx.ni_map[strand, level - 1] <= 0:
continue
subunit = int(ctx.ni_map[strand, level - 1])
atom_idx = int(ctx.molecule.subunit_boundaries[subunit - 1])
residue_name = str(ctx.molecule.residue_names[atom_idx]).strip()
chain_id = str(ctx.molecule.chain_ids[atom_idx]).strip() if ctx.molecule.chain_ids is not None else ""
parent_base = parent_base_name(residue_name)
sequence.append(str(parent_base or "N")[:1])
residues.append({
"level": level,
"chain_id": chain_id,
"residue_name": residue_name,
"parent_base": parent_base,
"residue_id": int(ctx.molecule.residue_ids[atom_idx]),
"insertion_code": str(ctx.molecule.insertion_codes[atom_idx]).strip() if ctx.molecule.insertion_codes is not None else "",
})
strands.append({
"strand": strand + 1,
"sequence": "".join(sequence),
"residues": residues,
})
return {
"n_strands": ctx.nst,
"n_levels": ctx.nux,
"strands": strands,
}
def _primary_sequence_levels(self) -> List[int]:
"""Return actual strand-1 residue levels, excluding Curves gap levels."""
ctx = self.runner.ctx
calc = self.runner.calc
return [level for level in range(1, ctx.nux + 1) if calc._has_level(0, level)]
def _step_label_between(self, strand: int, first_level: int, second_level: int) -> str:
calc = self.runner.calc
first = calc._residue_label(strand, first_level)
second = calc._residue_label(strand, second_level)
if first is None or second is None:
return "-"
first_base, _, first_id = first
second_base, _, second_id = second
return f"{first_base}{first_id:3d}/{second_base}{second_id:3d}"
def _build_dataframes(self) -> Dict[str, Any]:
"""Build the slim dataframe dictionary used by JSON, CSV, and viewers."""
ctx = self.runner.ctx
calc = self.runner.calc
annotations = self._annotations()
curvesplus_axis = self._uses_curvesplus_axis()
records = {}
# 1. Global Base-Axis
if not curvesplus_axis:
base_axis = []
for strand in range(ctx.nst):
for level in range(1, ctx.nux + 1):
if not calc._has_level(strand, level): continue
label = calc._residue_unit_label(strand, level)
if label is None: continue
res_name, res_id = label
p = ctx.params.helical[strand, level]
if not _all_finite([p[0], p[1], p[3], p[4]]):
continue
base_axis.append({
"strand": strand + 1,
"level": level,
"residue_name": res_name,
"residue_id": res_id,
"xdisp": p[0],
"ydisp": p[1],
"inclin": p[3],
"tip": p[4]
})
records["global_base_axis"] = base_axis
primary_levels = self._primary_sequence_levels()
# 2. Global/Curves+ Base Pair-Axis and global Base-Base
bp_axis = []
base_base = []
if ctx.cfg.comb and ctx.nst > 1:
# Curves+ builds its BP-axis from the strand-1/2 mean frame only.
# Legacy Curves reports a separate strand-1/partner series.
axis_partners = range(1, 2) if curvesplus_axis else range(1, ctx.nst)
for partner_strand in axis_partners:
base_base_by_level = (
dict(calc._global_base_base_values_by_level(partner_strand))
if not curvesplus_axis
else {}
)
for sequence_index, level in enumerate(primary_levels, start=1):
has_pair = calc._has_level(0, level) and calc._has_level(partner_strand, level)
duplex = calc._duplex_id(0, partner_strand, level)
axis_row = _nullable_parameter_record(
calc._global_base_pair_axis_values(partner_strand, level) if has_pair else None,
AXIS_PARAMETERS,
partner_strand=partner_strand + 1,
sequence_index=sequence_index,
level=level,
duplex=duplex,
)
axis_row["axis_weight"] = calc._axis_pair_weight(partner_strand, level)
bp_axis.append(axis_row)
if curvesplus_axis:
continue
base_base.append(_nullable_parameter_record(
base_base_by_level.get(level) if has_pair else None,
BASE_BASE_PARAMETERS,
partner_strand=partner_strand + 1,
sequence_index=sequence_index,
level=level,
duplex=duplex,
))
if curvesplus_axis:
records["curvesplus_base_pair_axis"] = bp_axis
else:
records["global_base_pair_axis"] = bp_axis
records["global_base_base"] = base_base
if ctx.cfg.comb and ctx.nst > 1:
local_base_base = []
for partner_strand in range(1, ctx.nst):
for sequence_index, level in enumerate(primary_levels, start=1):
has_pair = calc._has_level(0, level) and calc._has_level(partner_strand, level)
values = calc._local_base_base_values(partner_strand, level) if has_pair else None
local_base_base.append(_nullable_parameter_record(
values,
BASE_BASE_PARAMETERS,
partner_strand=partner_strand + 1,
sequence_index=sequence_index,
level=level,
duplex=calc._duplex_id(0, partner_strand, level),
))
records["local_base_base"] = local_base_base
# 3. Global Inter-Base
if not curvesplus_axis:
inter_base = []
for strand in range(ctx.nst):
for level in range(1, ctx.nux):
if not (calc._has_level(strand, level) and calc._has_level(strand, level+1)): continue
step = calc._step_label(strand, level+1)
row = _parameter_record(
ctx.params.inter_base[strand, level+1],
STEP_PARAMETERS,
strand=strand + 1,
level=level,
step=step,
)
if row is not None:
inter_base.append(row)
records["global_inter_base"] = inter_base
# 4. Local Inter-Base
local_inter_base = []
for strand in range(ctx.nst):
for level in range(1, ctx.nux):
if not (calc._has_level(strand, level) and calc._has_level(strand, level+1)): continue
step = calc._step_label(strand, level+1)
row = _parameter_record(
calc.local_inter_base[level+1, :, strand],
STEP_PARAMETERS,
strand=strand + 1,
level=level,
step=step,
)
if row is not None:
local_inter_base.append(row)
records["local_inter_base"] = local_inter_base
# 5. Inter-Base Pair
if ctx.cfg.comb and ctx.nst > 1:
global_inter_bp = []
local_inter_bp = []
pair_step_partners = range(1, 2) if curvesplus_axis else range(1, ctx.nst)
for partner_strand in pair_step_partners:
for sequence_index, (level, next_level) in enumerate(zip(primary_levels, primary_levels[1:]), start=1):
adjacent = next_level == level + 1
has_step = (
adjacent
and calc._has_level(0, level)
and calc._has_level(partner_strand, level)
and calc._has_level(0, next_level)
and calc._has_level(partner_strand, next_level)
)
duplex = f"{calc._duplex_id(0, partner_strand, level)}/{calc._duplex_id(0, partner_strand, next_level)}"
step = self._step_label_between(0, level, next_level)
if not curvesplus_axis:
global_row = _nullable_parameter_record(
calc._global_inter_base_pair_values(partner_strand, next_level) if has_step else None,
STEP_PARAMETERS,
partner_strand=partner_strand + 1,
sequence_index=sequence_index,
level=level,
next_level=next_level,
step=step,
duplex=duplex,
)
global_inter_bp.append(global_row)
local_values = None
if has_step:
local_values = (
calc._global_inter_base_pair_values(partner_strand, next_level)
if curvesplus_axis
else calc.local_inter_base_pair[next_level, :, partner_strand]
)
local_row = _nullable_parameter_record(
local_values,
STEP_PARAMETERS,
partner_strand=partner_strand + 1,
sequence_index=sequence_index,
level=level,
next_level=next_level,
step=step,
duplex=duplex,
)
local_inter_bp.append(local_row)
if not curvesplus_axis:
records["global_inter_base_pair"] = global_inter_bp
records["local_inter_base_pair"] = local_inter_bp
# 6. Backbone Parameters
backbone = []
for strand in range(ctx.nst):
for level in range(1, ctx.nux + 1):
if not calc._has_level(strand, level):
missing = [
"c1_c2", "c2_c3", "phase", "amplitude", "c1_prime", "c2_prime",
"c3_prime", "chi", "gamma", "delta", "epsilon", "zeta", "alpha", "beta",
]
backbone.append({
"strand": strand + 1,
"level": level,
"residue_name": None,
"residue_id": None,
"chain_id": None,
"insertion_code": None,
"gap": True,
"valid": False,
"status": "gap",
"warnings": ["topology_gap", "missing_parameter_values"],
"missing_parameters": missing,
"pucker": None,
**{name: None for name in missing},
})
continue
subunit = abs(int(ctx.ni_map[strand, level - 1]))
atom_idx = int(ctx.molecule.subunit_boundaries[subunit - 1])
res_name = str(ctx.molecule.residue_names[atom_idx]).strip()
res_id = int(ctx.molecule.residue_ids[atom_idx])
chain_id = str(ctx.molecule.chain_ids[atom_idx]).strip() if ctx.molecule.chain_ids is not None else ""
insertion_code = str(ctx.molecule.insertion_codes[atom_idx]).strip() if ctx.molecule.insertion_codes is not None else ""
tor = ctx.backbone.torsions[strand, level]
pucker = ctx.backbone.sugar_pucker[strand, level]
phase = _backbone_number(pucker[1])
pucker_label = None
if phase is not None:
pucker_index = int((phase % 360.0) / 36.0)
pucker_index = max(0, min(len(SUGAR_PUCKERS) - 1, pucker_index))
pucker_label = SUGAR_PUCKERS[pucker_index]
row = {
"strand": strand + 1,
"level": level,
"residue_name": res_name,
"residue_id": res_id,
"chain_id": chain_id,
"insertion_code": insertion_code,
"gap": False,
"c1_c2": _backbone_number(tor[4]),
"c2_c3": _backbone_number(tor[5]),
"phase": phase,
"amplitude": _backbone_number(pucker[0]),
"pucker": pucker_label,
"c1_prime": _backbone_number(tor[0]),
"c2_prime": _backbone_number(tor[1]),
"c3_prime": _backbone_number(tor[2]),
"chi": _backbone_number(tor[12]),
"gamma": _backbone_number(tor[8]),
"delta": _backbone_number(tor[9]),
"epsilon": _backbone_number(tor[6]),
"zeta": _backbone_number(tor[7]),
"alpha": _backbone_number(tor[10]),
"beta": _backbone_number(tor[11]),
}
parameter_names = (
"c1_c2", "c2_c3", "phase", "amplitude", "c1_prime", "c2_prime",
"c3_prime", "chi", "gamma", "delta", "epsilon", "zeta", "alpha", "beta",
)
missing_parameters = [name for name in parameter_names if row[name] is None]
warnings = []
if bool(ctx.backbone.flag[strand, level]):
warnings.append("incomplete_sugar_geometry")
if missing_parameters:
warnings.append("missing_parameter_values")
row.update({
"valid": not warnings,
"status": "complete" if not warnings else "partial",
"warnings": warnings,
"missing_parameters": missing_parameters,
})
backbone.append(row)
records["backbone"] = backbone
if not curvesplus_axis:
records["global_axis_curvature"] = self._axis_curvature_records()
bending_rows, bending_summary = self._axis_bending_records()
records["global_axis_bending"] = bending_rows
records["global_axis_bending_summary"] = bending_summary
# 7. Groove Parameters
if ctx.cfg.comb and ctx.nst > 1 and getattr(ctx.cfg, "grv", False):
if not hasattr(calc, "groove_params"):
if self.fast_groove:
self._capture(lambda: calc.groove(fast=True))
else:
self._capture(calc.groove)
if hasattr(calc, "groove_params") and calc.groove_params:
records["groove"] = self._groove_records(calc.groove_params)
records["base_pair_observations"] = self._base_pair_observation_records(
annotations
)
records["annotations"] = self._slim_annotation_records(annotations)
return records
def _uses_curvesplus_axis(self) -> bool:
cfg = self.runner.ctx.cfg
axis = str(getattr(cfg, "axis_convention", "global")).lower().replace("-", "_")
return (
str(getattr(cfg, "frame_convention", "standard")).lower() == "standard"
and axis in {"local", "curvesplus", "curves_plus", "curves+", "canal"}
)
def _render_section_l_local_base_base(self, records: Dict[str, Any]) -> str:
rows = records.get("local_base_base", [])
if not rows:
return ""
grouped: Dict[int, List[Dict[str, Any]]] = {}
for row in rows:
grouped.setdefault(int(row.get("partner_strand", 0)), []).append(row)
lines = [
" --------------------------------",
" |L| Local Base-Base Parameters |",
" --------------------------------",
]
for partner_strand in sorted(grouped):
lines.extend([
"",
f" Strand 1 with strand {partner_strand} ...",
"",
" Duplex Shear Stretch Stagger Buckle Propel Opening",
" (Sx) (Sy) (Sz) (kappa) (omega) (sigma)",
])
for row in sorted(grouped[partner_strand], key=lambda item: int(item.get("level", 0))):
if not _all_finite([row.get(name) for name in BASE_BASE_PARAMETERS]):
continue
lines.append(
f" {int(row.get('level', 0)):3d}) {str(row.get('duplex', '')):<9s} "
f"{float(row.get('shear', 0.0)):8.2f} {float(row.get('stretch', 0.0)):8.2f} "
f"{float(row.get('stagger', 0.0)):8.2f} {float(row.get('buckle', 0.0)):8.2f} "
f"{float(row.get('propel', 0.0)):8.2f} {float(row.get('opening', 0.0)):8.2f}"
)
return "\n".join(lines)
def _axis_curvature_records(self) -> List[Dict[str, Any]]:
ctx = self.runner.ctx
calc = self.runner.calc
rows = []
section_strands = [0] if ctx.cfg.comb else list(range(ctx.nst))
for strand in section_strands:
_, _, start, stop = calc._axis_bounds(strand)
for level in range(start + 1, stop + 1):
values = calc.vkin[level, :, strand]
if not _all_finite(values[:7]):
continue
step = calc._step_label(strand, level) if calc._has_level(strand, level - 1) and calc._has_level(strand, level) else "-"
rows.append({
"strand": strand + 1,
"level": level - 1,
"next_level": level,
"step": step,
"ax": values[0],
"ay": values[1],
"ainc": values[2],
"atip": values[3],
"adis": values[4],
"angle": values[5],
"path": values[6],
"dc": int(calc.dcod[level, strand]) if hasattr(calc, "dcod") else 0,
})
return rows
def _axis_bending_records(self):
ctx = self.runner.ctx
calc = self.runner.calc
rows = []
summary = []
section_strands = [0] if ctx.cfg.comb else list(range(ctx.nst))
for strand in section_strands:
_, _, start, stop = calc._axis_bounds(strand)
axis_points = ctx.params.ox if ctx.cfg.comb else self.runner.opt.hho[:, :, strand]
path_length = float(sum(calc.vkin[level, 6, strand] for level in range(start + 1, stop + 1)))
end_to_end = float(np.linalg.norm(axis_points[stop] - axis_points[start]))
shortening = 100.0 * (1.0 - end_to_end / path_length) if path_length > 0.0 else 0.0
start_point = axis_points[start]
end_to_end_vector = axis_points[stop] - start_point
end_to_end_unit = end_to_end_vector / (np.linalg.norm(end_to_end_vector) + 1e-12)
for level in range(start, stop + 1):
rel = axis_points[level] - start_point
projection = float(np.dot(rel, end_to_end_unit))
offset = float(np.linalg.norm(rel - end_to_end_unit * projection))
local_direction = 0.0 if level in {start, stop} else float(calc.bend[level, strand])
label = calc._residue_label(strand, level)
if label is None:
residue_name = "-"
residue_id = None
else:
residue_name, _, residue_id = label
rows.append({
"strand": strand + 1,
"level": level,
"residue_name": residue_name,
"residue_id": residue_id,
"offset": _json_number(offset),
"local_direction": _json_number(local_direction),
})
summary.append({
"strand": strand + 1,
"path_length": _json_number(path_length),
"end_to_end": _json_number(end_to_end),
"shortening_percent": _json_number(shortening),
"overall_bend_uu": _json_number(calc.bend[start, strand]),
"overall_bend_pp": _json_number(calc.bend[stop, strand]),
})
return rows, summary
@staticmethod
def _base_pair_rows(
annotations: Dict[str, List[Dict[str, Any]]],
) -> List[Dict[str, Any]]:
return annotations.get(
"frame_base_pair_observations",
annotations.get("base_pair_annotations", []),
)
def _base_pair_observation_records(
self,
annotations: Dict[str, List[Dict[str, Any]]],
) -> List[Dict[str, Any]]:
"""Return one compact record for each reference or newly observed pair."""
rows: List[Dict[str, Any]] = []
for row in self._base_pair_rows(annotations):
rows.append({
"pair_id": row.get("pair_id", ""),
"reference_pair": bool(row.get("reference_pair", True)),
"level": row.get("level"),
"residue_1": row.get("residue_1"),
"residue_2": row.get("residue_2"),
"pair_status": row.get("pair_status", "uncertain"),
"pairing_mode": row.get("pairing_mode", ""),
"observed_lw_family": base_pair_observed_geometry_tag(row),
"reference_lw_family": row.get("reference_lw_family", ""),
"candidate_mode": row.get("candidate_mode", ""),
"classification_status": row.get(
"classification_status", "unassigned"
),
"helical_context": row.get("helical_context", ""),
"diagnostic_flags": list(row.get("diagnostic_flags") or []),
})
return rows
def _slim_annotation_records(self, annotations: Dict[str, List[Dict[str, Any]]]) -> List[Dict[str, Any]]:
rows: List[Dict[str, Any]] = []
for row in annotations.get("modified_base_annotations", []):
rows.append({
"annotation_type": "modified_base",
"severity": "info",
"level": row.get("level"),
"location": f"strand {row.get('strand')} level {row.get('level')}",
"code": "modified_base",
"message": f"{row.get('residue_name', '')} fitted with {row.get('parent_base', '')} parent-base template.".strip(),
"residue": row.get("residue_name"),
"parent_base": row.get("parent_base"),
})
for row in annotations.get("backbone_links", []):
relation = row.get("topology_relation", "")
if relation == "sequential":
continue
rows.append({
"annotation_type": "backbone_link",
"severity": "warn" if row.get("topology_discontinuity") else "info",
"level": row.get("level"),
"location": f"strand {row.get('strand')} level {row.get('level')}",
"code": relation,
"message": (
f"O3'-P link {row.get('source_residue_id')} -> "
f"{row.get('target_residue_id')} follows molecular connectivity."
),
"source_subunit": row.get("source_subunit"),
"target_subunit": row.get("target_subunit"),
"bond_source": row.get("bond_source"),
"distance": row.get("distance"),
})
return rows
def _groove_records(self, groove_params: Dict[str, Any]) -> Dict[str, Any]:
groove = _to_jsonable(groove_params)
if self.include_visualization:
return groove
for level_data in groove.get("data", {}).values():
for sub_level_data in level_data.get("sub_levels", {}).values():
sub_level_data.pop("geometry", None)
return groove
def _annotations(self) -> Dict[str, List[Dict[str, Any]]]:
if self._annotation_cache is None:
ctx = self.runner.ctx
existing = getattr(ctx, "annotations", {})
frame_observations = existing.get("frame_base_pair_observations")
self._annotation_cache = (
existing
if "base_pair_annotations" in existing
else annotate_context(ctx)
)
if frame_observations is not None:
self._annotation_cache["frame_base_pair_observations"] = frame_observations
from pycurves_lib.core.parameter_conventions import (
_annotate_pair_normal_branches,
)
_annotate_pair_normal_branches(
ctx,
getattr(ctx, "pair_normal_signs", {}) or {},
getattr(ctx, "pair_normal_branch_modes", {}) or {},
getattr(ctx, "pair_glycosidic_details", {}) or {},
)
return self._annotation_cache
def _header(self, header_name: str) -> str:
cfg = self.runner.ctx.cfg
namelist = self._read_namelist()
today = datetime.now().strftime("%d %b %y")
file_name = self.runner.pdbfile or namelist.get("file", "")
lis_name = namelist.get("lis", Path(self.runner.inpfile or "").stem)
pdb_name = namelist.get("pdb", "")
title = f"****** {header_name.upper():^17s} *****"
return "\n".join([
"",
" *********************************** **************",
f" {title:<35s} * {today:>9s} *",
" *********************************** **************",
"",
"",
f" FILE : {str(file_name):<32s} LIS : {str(lis_name):<32s}",
f" dna : {namelist.get('dna', ''):<32s} axin : {namelist.get('axin', ''):<32s}",
f" axout: {namelist.get('axout', ''):<32s} daf : {namelist.get('daf', ''):<32s}",
f" PDB : {str(pdb_name):<32s}",
"",
f" acc : {cfg.acc:8.3f} wid : {cfg.wid:8.3f}",
"",
(
f" maxn : {cfg.maxn:5d} ior : {cfg.ior:5d} ibond: {cfg.ibond:5d} "
f"splin: {cfg.spline:5d} break: {cfg.break_lvl:5d}"
),
f" nleve: {getattr(cfg, 'nlevel', 3):5d} nbac : {getattr(cfg, 'nbac', 7):5d}",
"",
(
f" ends : {self._tf(cfg.ends):>5s} supp : {self._tf(cfg.supp):>5s} "
f"COMB : {self._tf(cfg.comb):>5s} dinu : {self._tf(cfg.dinu):>5s} mini : {self._tf(cfg.mini):>5s}"
),
(
f" rest : {self._tf(cfg.rest):>5s} line : {self._tf(cfg.line):>5s} "
f"zaxe : {self._tf(cfg.zaxe):>5s} FIT : {self._tf(cfg.fit):>5s} test : {self._tf(cfg.test):>5s}"
),
(
f" GRV : {self._tf(cfg.grv):>5s} old : {self._tf(cfg.old):>5s} "
f"axonl: {self._tf(cfg.axonly):>5s}"
),
])
def _read_namelist(self) -> Dict[str, str]:
if not self.runner.inpfile:
return {}
text = Path(self.runner.inpfile).read_text(encoding="utf-8", errors="ignore")
values = {}
for key, value in re.findall(r"([A-Za-z][A-Za-z0-9_]*)\s*=\s*([^,\s&]+)", text):
values[key.lower()] = value.strip().strip("'\"")
return values
@staticmethod
def _capture(func) -> str:
buf = io.StringIO()
with contextlib.redirect_stdout(buf):
func()
return buf.getvalue()
@staticmethod
def _tf(value: bool) -> str:
return "T" if value else "F"
def _require_results(self) -> None:
if self.runner.ctx is None or self.runner.opt is None or self.runner.calc is None:
raise RuntimeError("Call analyze() before rendering output.")