import numpy as np
import networkx as nx
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize, LinearSegmentedColormap, ListedColormap
from matplotlib.patches import FancyArrowPatch
from adjustText import adjust_text
# [2026-10-24 added: the direction rule shared with the ranked lists]
from twinfer.inference.ranked_edges import pair_directions
__version__ = "2026-07-16-one-reciprocal-edge-curved-v6"
# ============================================================
# 1. Colormap
# ============================================================
def make_reds_blues_colormap(vmin=-0.05, vmax=0.18):
"""
Negative weights: red
Zero: white
Positive weights: blue
"""
vmin = float(vmin)
vmax = float(vmax)
if np.isclose(vmin, vmax):
if np.isclose(vmin, 0.0):
vmin = -1.0
vmax = 1.0
elif vmin > 0:
vmin = 0.0
else:
vmax = 0.0
zero_position = np.clip(
(0.0 - vmin) / (vmax - vmin),
0.0,
1.0,
)
n_total = 256
n_reds = int(round(zero_position * n_total))
n_blues = n_total - n_reds
denominator = max(
abs(vmin),
abs(vmax),
1e-12,
)
red_intensity = abs(vmin) / denominator
blue_intensity = abs(vmax) / denominator
low_clip = 0.50
color_parts = []
if n_reds > 0:
reds = plt.cm.Reds(
np.linspace(
low_clip
+ (1.0 - low_clip) * red_intensity,
low_clip,
n_reds,
)
)
color_parts.append(reds)
whites = np.ones((1, 4)) # pure white at 0
color_parts.append(whites)
if n_blues > 0:
blues = plt.cm.Blues(
np.linspace(
low_clip,
low_clip
+ (1.0 - low_clip) * blue_intensity,
n_blues,
)
)
color_parts.append(blues)
colors = np.vstack(color_parts)
# ListedColormap uses `colors` directly as the lookup table (no
# resampling), so the single pure-white row placed at index n_reds
# stays exactly white. LinearSegmentedColormap.from_list would instead
# resample this list down to its own N-entry LUT, and since the white
# row is a single entry out of ~257, it almost always got interpolated
# away -- 0 rendered as a red/blue blend, never pure white.
return ListedColormap(colors, name="RedsBlues")
# ============================================================
# 2. Angular helpers
# ============================================================
def _wrap_angle(angle):
"""Map an angle to [-pi, pi]."""
return np.arctan2(
np.sin(angle),
np.cos(angle),
)
def _circular_angle_distance(angle_a, angle_b):
"""Smallest absolute angular difference."""
return abs(
_wrap_angle(angle_a - angle_b)
)
def _circular_mean(angles):
"""Circular mean of a sequence of angles."""
angles = np.asarray(
angles,
dtype=float,
)
return np.arctan2(
np.mean(np.sin(angles)),
np.mean(np.cos(angles)),
)
def _centered_slots(count):
"""
Examples
--------
1 -> [0]
2 -> [-0.5, +0.5]
3 -> [-1, 0, +1]
"""
if count <= 0:
return np.empty(
0,
dtype=float,
)
return (
np.arange(count, dtype=float)
- (count - 1) / 2.0
)
# ============================================================
# 3. Group crowded terminals
# ============================================================
def _group_terminals(
terminals,
angle_tolerance,
):
"""
Group terminals whose natural contact directions are close.
terminals:
[((source, target), angle), ...]
"""
if not terminals:
return []
proximity_graph = nx.Graph()
proximity_graph.add_nodes_from(
range(len(terminals))
)
for first in range(len(terminals)):
for second in range(
first + 1,
len(terminals),
):
angle_difference = (
_circular_angle_distance(
terminals[first][1],
terminals[second][1],
)
)
if (
angle_difference
<= float(angle_tolerance)
):
proximity_graph.add_edge(
first,
second,
)
groups = []
for component in nx.connected_components(
proximity_graph
):
groups.append(
[
terminals[index]
for index in component
]
)
return groups
def _assign_desired_contact_angles(
terminals,
angle_tolerance,
terminal_spread,
):
"""
Assign evenly spaced contact angles to crowded terminals.
"""
desired_angles = {}
crowded_edges = set()
groups = _group_terminals(
terminals,
angle_tolerance=angle_tolerance,
)
for group in groups:
if len(group) == 1:
edge, natural_angle = group[0]
desired_angles[edge] = float(
natural_angle
)
continue
natural_angles = np.array(
[
natural_angle
for _, natural_angle in group
],
dtype=float,
)
mean_angle = _circular_mean(
natural_angles
)
ordered_group = sorted(
group,
key=lambda item: _wrap_angle(
item[1] - mean_angle
),
)
slots = _centered_slots(
len(ordered_group)
)
for (
edge,
natural_angle,
), slot in zip(
ordered_group,
slots,
):
desired_angle = (
mean_angle
+ slot * float(terminal_spread)
)
desired_angles[edge] = float(
_wrap_angle(desired_angle)
)
crowded_edges.add(edge)
return desired_angles, crowded_edges
# ============================================================
# 4. Straight-line terminal shifts
# ============================================================
def _build_edge_shift_fractions(
graph,
positions,
terminal_angle_tolerance,
terminal_spread,
reciprocal_separation,
max_shift_fraction,
):
"""
Compute a perpendicular whole-line shift for each edge.
The shift is represented as a fraction of node radius.
"""
source_desired_angles = {}
target_desired_angles = {}
source_crowded = set()
target_crowded = set()
# --------------------------------------------------------
# Assign desired source and target contact angles
# --------------------------------------------------------
for node in graph.nodes():
node_position = np.asarray(
positions[node],
dtype=float,
)
outgoing_terminals = []
for _, target in graph.out_edges(node):
target_position = np.asarray(
positions[target],
dtype=float,
)
vector = (
target_position
- node_position
)
natural_source_angle = np.arctan2(
vector[1],
vector[0],
)
outgoing_terminals.append(
(
(node, target),
natural_source_angle,
)
)
(
node_source_angles,
node_source_crowded,
) = _assign_desired_contact_angles(
outgoing_terminals,
angle_tolerance=(
terminal_angle_tolerance
),
terminal_spread=terminal_spread,
)
source_desired_angles.update(
node_source_angles
)
source_crowded.update(
node_source_crowded
)
incoming_terminals = []
for source, _ in graph.in_edges(node):
source_position = np.asarray(
positions[source],
dtype=float,
)
vector = (
source_position
- node_position
)
natural_target_angle = np.arctan2(
vector[1],
vector[0],
)
incoming_terminals.append(
(
(source, node),
natural_target_angle,
)
)
(
node_target_angles,
node_target_crowded,
) = _assign_desired_contact_angles(
incoming_terminals,
angle_tolerance=(
terminal_angle_tolerance
),
terminal_spread=terminal_spread,
)
target_desired_angles.update(
node_target_angles
)
target_crowded.update(
node_target_crowded
)
# --------------------------------------------------------
# Convert contact-angle requests into line shifts
# --------------------------------------------------------
shift_fractions = {}
for source, target in graph.edges():
edge = (
source,
target,
)
source_position = np.asarray(
positions[source],
dtype=float,
)
target_position = np.asarray(
positions[target],
dtype=float,
)
direction = (
target_position
- source_position
)
natural_source_angle = np.arctan2(
direction[1],
direction[0],
)
natural_target_angle = _wrap_angle(
natural_source_angle + np.pi
)
source_shift = 0.0
target_shift = 0.0
if edge in source_crowded:
desired_source_angle = (
source_desired_angles[edge]
)
source_delta = _wrap_angle(
desired_source_angle
- natural_source_angle
)
source_shift = np.sin(
source_delta
)
if edge in target_crowded:
desired_target_angle = (
target_desired_angles[edge]
)
target_delta = _wrap_angle(
desired_target_angle
- natural_target_angle
)
# Target orientation is reversed relative to source.
target_shift = -np.sin(
target_delta
)
# Target terminal receives priority because arrowheads
# and repression bars are drawn at the target.
if edge in target_crowded:
shift_fraction = target_shift
if (
edge in source_crowded
and not np.isclose(
source_shift,
0.0,
)
and not np.isclose(
target_shift,
0.0,
)
and np.sign(source_shift)
== np.sign(target_shift)
):
shift_fraction += (
0.15 * source_shift
)
elif edge in source_crowded:
shift_fraction = source_shift
else:
shift_fraction = 0.0
shift_fractions[edge] = float(
shift_fraction
)
# --------------------------------------------------------
# Reciprocal pairs still receive slight endpoint separation.
#
# This is necessary because curving only changes the middle of
# the edge; without endpoint separation, arrowheads/bars could
# still occupy exactly the same point.
# --------------------------------------------------------
processed_pairs = set()
for source, target in graph.edges():
if not graph.has_edge(
target,
source,
):
continue
pair_key = frozenset(
(source, target)
)
if pair_key in processed_pairs:
continue
processed_pairs.add(
pair_key
)
forward_edge = (
source,
target,
)
reverse_edge = (
target,
source,
)
forward_shift = shift_fractions.get(
forward_edge,
0.0,
)
reverse_shift = shift_fractions.get(
reverse_edge,
0.0,
)
# Reverse edges have opposite perpendicular vectors.
# Therefore physical separation is proportional to:
#
# abs(forward_shift + reverse_shift)
current_separation = abs(
forward_shift + reverse_shift
)
required_separation = float(
reciprocal_separation
)
if (
current_separation
< required_separation
):
pair_sum = (
forward_shift
+ reverse_shift
)
direction_sign = (
1.0
if pair_sum >= 0.0
else -1.0
)
additional_shift = (
required_separation
- current_separation
) / 2.0
forward_shift += (
direction_sign
* additional_shift
)
reverse_shift += (
direction_sign
* additional_shift
)
shift_fractions[forward_edge] = float(
np.clip(
forward_shift,
-max_shift_fraction,
max_shift_fraction,
)
)
shift_fractions[reverse_edge] = float(
np.clip(
reverse_shift,
-max_shift_fraction,
max_shift_fraction,
)
)
for edge in shift_fractions:
shift_fractions[edge] = float(
np.clip(
shift_fractions[edge],
-max_shift_fraction,
max_shift_fraction,
)
)
return shift_fractions
# ============================================================
# 5. Radial staggering of crowded target symbols
# ============================================================
def _build_terminal_depth_slots(
graph,
positions,
terminal_angle_tolerance,
):
"""
Give crowded incoming symbols different radial distances.
This helps when large arrowheads or repression bars are located
close together on the same side of a node.
"""
depth_slots = {}
for target in graph.nodes():
target_position = np.asarray(
positions[target],
dtype=float,
)
incoming_terminals = []
for source, _ in graph.in_edges(target):
source_position = np.asarray(
positions[source],
dtype=float,
)
vector = (
source_position
- target_position
)
angle = np.arctan2(
vector[1],
vector[0],
)
incoming_terminals.append(
(
(source, target),
angle,
)
)
groups = _group_terminals(
incoming_terminals,
angle_tolerance=(
terminal_angle_tolerance
),
)
for group in groups:
if len(group) == 1:
edge, angle = group[0]
depth_slots[edge] = 0.0
continue
group_angles = np.array(
[
angle
for _, angle in group
],
dtype=float,
)
mean_angle = _circular_mean(
group_angles
)
# The terminal nearest the center of the sector stays
# closest to the node. Others move progressively back.
ordered_group = sorted(
group,
key=lambda item: abs(
_wrap_angle(
item[1] - mean_angle
)
),
)
for depth_level, (
edge,
angle,
) in enumerate(ordered_group):
depth_slots[edge] = float(
depth_level
)
return depth_slots
# ============================================================
# 6. Shifted straight-line geometry
# ============================================================
def _shifted_circle_contacts(
source_center,
target_center,
node_radius,
shift_fraction,
):
"""
Shift the entire source-target line and intersect it with both
node circles.
"""
source_center = np.asarray(
source_center,
dtype=float,
)
target_center = np.asarray(
target_center,
dtype=float,
)
center_difference = (
target_center
- source_center
)
center_distance = float(
np.linalg.norm(
center_difference
)
)
if center_distance <= 1e-12:
raise ValueError(
"Source and target positions must differ."
)
unit = (
center_difference
/ center_distance
)
perpendicular = np.array(
[
-unit[1],
unit[0],
]
)
shift_fraction = float(
np.clip(
shift_fraction,
-0.999,
0.999,
)
)
shift_distance = (
shift_fraction
* float(node_radius)
)
radial_fraction = np.sqrt(
max(
1.0
- shift_fraction ** 2,
0.0,
)
)
radial_distance = (
radial_fraction
* float(node_radius)
)
start = (
source_center
+ perpendicular * shift_distance
+ unit * radial_distance
)
end = (
target_center
+ perpendicular * shift_distance
- unit * radial_distance
)
return (
start,
end,
unit,
perpendicular,
)
# ============================================================
# 7. Curve tangent
# ============================================================
def _arc3_end_tangent(
start,
end,
radius,
):
"""
Approximate the tangent of Matplotlib's arc3 curve at the end.
This is used to orient a curved repression bar perpendicular
to the curve rather than perpendicular to the center-to-center
line.
"""
start = np.asarray(
start,
dtype=float,
)
end = np.asarray(
end,
dtype=float,
)
difference = (
end - start
)
control_point = (
0.5 * (start + end)
+ float(radius)
* np.array(
[
difference[1],
-difference[0],
]
)
)
tangent = (
end - control_point
)
tangent_norm = float(
np.linalg.norm(tangent)
)
if tangent_norm <= 1e-12:
tangent = difference
tangent_norm = float(
np.linalg.norm(tangent)
)
if tangent_norm <= 1e-12:
return np.array(
[1.0, 0.0]
)
return (
tangent / tangent_norm
)
# ============================================================
# 8. Main plotting function
# ============================================================
[docs]
def plot_grn(
matrix,
node_labels=None,
figsize=None,
node_color="#1a1a1a",
font_size=None,
label_font_size=None,
font_color="black",
title=None,
ax=None,
fixed_pos=None,
show_colorbar=True,
vmin=None,
vmax=None,
*,
terminal_spread=0.38,
terminal_angle_tolerance=1.10,
reciprocal_separation=0.30,
reciprocal_curve=0.18,
terminal_depth_step=0.35,
max_shift_fraction=0.92,
edge_halo=True,
edge_styles=None,
undirected_edges=None,
colorbar_label="Edge Weight",
# [2026-10-24 added: edge_colors]
edge_colors=None,
# [2026-10-24 added: edge_signs]
edge_signs=None,
):
"""
Plot a signed, weighted directed GRN.
Matrix convention
-----------------
matrix[i, j] represents:
node i -> node j
Positive value:
activation arrow
Negative value:
repression bar
Reciprocal-pair rule
--------------------
If both u -> v and v -> u exist:
lower-index source -> higher-index target:
straight
higher-index source -> lower-index target:
curved
Parameters
----------
terminal_spread
Angular spacing between neighboring crowded terminals.
terminal_angle_tolerance
Natural terminal directions within this angular distance are
treated as crowded.
reciprocal_separation
Small endpoint separation retained for reciprocal edges.
reciprocal_curve
Magnitude of the curve used for one edge in each reciprocal
pair.
terminal_depth_step
Radial staggering of crowded target symbols, in node-radius
units.
edge_styles
Optional dict {(source, target): Matplotlib line style}, e.g.
{(0, 1): "--"} draws the edge 0 -> 1 dashed. Edges not in the
dict are solid. Arrowheads and repression bars stay solid.
undirected_edges
Optional collection of (source, target) edges drawn as a plain
line without arrowhead or repression bar (an association with
no direction). Give each such pair once: matrix[source, target]
must be non-zero and matrix[target, source] zero.
colorbar_label
Label of the colour bar (when show_colorbar is True).
edge_colors
Optional dict {(source, target): Matplotlib colour}, e.g.
{(0, 1): "tab:blue"}: the edge 0 -> 1 (line, arrowhead or bar)
is drawn in that colour instead of the colour of its weight.
Edges not in the dict are coloured by weight.
edge_signs
Optional dict {(source, target): +1 or -1}: arrowhead (+1) or
repression bar (-1) of that edge, instead of the sign of its
weight; the colour still follows the weight. Edges not in the
dict take the sign of their weight.
Notes
-----
`network_from_scenarios` builds the matrix, edge_styles,
undirected_edges and edge_colors from the output of
twinfer.inference.covariance_decomposition_z_scores.infer_with_twinfer.
"""
matrix = np.asarray(
matrix,
dtype=float,
)
if matrix.ndim != 2:
raise ValueError(
f"Matrix is not 2D: shape={matrix.shape}"
)
if matrix.shape[0] != matrix.shape[1]:
raise ValueError(
f"Matrix is not square: shape={matrix.shape}"
)
if not np.all(np.isfinite(matrix)):
raise ValueError(
"Matrix contains NaN or infinite values."
)
n_nodes = int(
matrix.shape[0]
)
if n_nodes == 0:
raise ValueError(
"Matrix must contain at least one node."
)
if figsize is None:
side_length = max(
6.0,
min(
20.0,
np.sqrt(n_nodes) * 2.5,
),
)
figsize = (
side_length,
side_length,
)
node_radius = max(
0.012,
min(
0.05,
0.15 / np.sqrt(n_nodes),
),
)
if font_size is None:
font_size = max(
6,
min(
12,
int(
35
/ np.log1p(n_nodes)
),
),
)
if label_font_size is None:
label_font_size = max(
font_size + 2,
int(
round(
font_size * 1.30
)
),
)
if node_labels is None:
node_labels = [
f"$g_{{{node + 1}}}$"
for node in range(n_nodes)
]
if len(node_labels) != n_nodes:
raise ValueError(
"node_labels must contain one label per node."
)
# --------------------------------------------------------
# Build graph
# --------------------------------------------------------
graph = nx.DiGraph()
graph.add_nodes_from(
range(n_nodes)
)
for source in range(n_nodes):
for target in range(n_nodes):
if (
source != target
and matrix[source, target] != 0
):
graph.add_edge(
source,
target,
weight=float(
matrix[source, target]
),
)
# --------------------------------------------------------
# Layout
# --------------------------------------------------------
if fixed_pos is not None:
missing_nodes = [
node
for node in graph.nodes()
if node not in fixed_pos
]
if missing_nodes:
raise ValueError(
"fixed_pos is missing nodes: "
f"{missing_nodes}"
)
positions = {
node: np.asarray(
fixed_pos[node],
dtype=float,
)
for node in graph.nodes()
}
else:
# Clockwise circular fallback.
angles = np.linspace(
3.0 * np.pi / 4.0,
3.0 * np.pi / 4.0
- 2.0 * np.pi,
n_nodes,
endpoint=False,
)
positions = {
node: np.array(
[
0.5
+ 0.36 * np.cos(angle),
0.5
+ 0.36 * np.sin(angle),
]
)
for node, angle in enumerate(
angles
)
}
# --------------------------------------------------------
# Figure
# --------------------------------------------------------
if ax is None:
fig, ax = plt.subplots(
figsize=figsize
)
own_figure = True
else:
fig = ax.get_figure()
own_figure = False
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
ax.set_aspect("equal")
ax.axis("off")
if title:
ax.set_title(
title,
fontsize=font_size + 2,
fontweight="bold",
pad=10,
)
figure_scale = (
min(figsize) / 8.0
)
head_scale = figure_scale
# --------------------------------------------------------
# Edge colors
# --------------------------------------------------------
edge_weights = np.array(
[
data["weight"]
for _, _, data
in graph.edges(data=True)
],
dtype=float,
)
if edge_weights.size == 0:
minimum_weight = (
-1.0
if vmin is None
else float(vmin)
)
maximum_weight = (
1.0
if vmax is None
else float(vmax)
)
else:
minimum_weight = (
min(
0.0,
float(edge_weights.min()),
)
if vmin is None
else float(vmin)
)
maximum_weight = (
max(
0.0,
float(edge_weights.max()),
)
if vmax is None
else float(vmax)
)
if np.isclose(
minimum_weight,
maximum_weight,
):
if np.isclose(
minimum_weight,
0.0,
):
minimum_weight = -1.0
maximum_weight = 1.0
elif minimum_weight > 0:
minimum_weight = 0.0
else:
maximum_weight = 0.0
colormap = make_reds_blues_colormap(
minimum_weight,
maximum_weight,
)
color_norm = Normalize(
vmin=minimum_weight,
vmax=maximum_weight,
)
# --------------------------------------------------------
# Terminal positions
# --------------------------------------------------------
edge_shift_fractions = (
_build_edge_shift_fractions(
graph=graph,
positions=positions,
terminal_angle_tolerance=(
terminal_angle_tolerance
),
terminal_spread=terminal_spread,
reciprocal_separation=(
reciprocal_separation
),
max_shift_fraction=(
max_shift_fraction
),
)
)
terminal_depth_slots = (
_build_terminal_depth_slots(
graph=graph,
positions=positions,
terminal_angle_tolerance=(
terminal_angle_tolerance
),
)
)
line_width = (
2.0 * figure_scale
)
terminal_gap = (
node_radius * 0.40
)
edge_styles = dict(edge_styles or {})
edge_colors = dict(edge_colors or {})
edge_signs = dict(edge_signs or {})
undirected_edges = {
tuple(edge)
for edge in (undirected_edges or ())
}
# --------------------------------------------------------
# Draw edges
# --------------------------------------------------------
for source, target, data in graph.edges(
data=True
):
edge = (
source,
target,
)
weight = float(
data["weight"]
)
# [2026-10-23 replaced: undirected edges get no repression bar]
# repression = (
# weight < 0
# )
undirected = (
edge in undirected_edges
)
# [2026-10-24 replaced: edge_signs sets arrow / bar apart from the colour]
# repression = (
# weight < 0
# and not undirected
# )
repression = (
edge_signs.get(edge, weight) < 0
and not undirected
)
line_style = edge_styles.get(
edge,
"solid",
)
# [2026-10-24 replaced: edge_colors overrides the colour of the weight]
# color = colormap(
# color_norm(weight)
# )
color = edge_colors.get(
edge,
colormap(
color_norm(weight)
),
)
shift_fraction = (
edge_shift_fractions.get(
edge,
0.0,
)
)
(
start,
end,
unit,
perpendicular,
) = _shifted_circle_contacts(
source_center=positions[source],
target_center=positions[target],
node_radius=node_radius,
shift_fraction=shift_fraction,
)
terminal_depth_level = (
terminal_depth_slots.get(
edge,
0.0,
)
)
terminal_depth = (
terminal_depth_level
* float(terminal_depth_step)
* node_radius
)
# ----------------------------------------------------
# Reciprocal curve rule
# ----------------------------------------------------
is_reciprocal = graph.has_edge(
target,
source,
)
curve_this_edge = (
is_reciprocal
and source > target
)
if curve_this_edge:
# Bow the curved edge farther toward its current shifted
# side rather than bending it back toward the straight edge.
if np.isclose(
shift_fraction,
0.0,
):
curve_sign = 1.0
else:
# For Matplotlib arc3, positive rad bows toward
# -perpendicular, hence the negative sign.
curve_sign = -np.sign(
shift_fraction
)
curve_radius = (
curve_sign
* abs(float(reciprocal_curve))
)
else:
curve_radius = 0.0
connection_style = (
f"arc3,rad={curve_radius}"
)
# ----------------------------------------------------
# Repression
# ----------------------------------------------------
if repression:
bar_end = (
end
- unit
* (
terminal_gap
+ terminal_depth
)
)
if curve_this_edge:
if edge_halo:
ax.annotate(
"",
xy=bar_end,
xytext=start,
arrowprops=dict(
arrowstyle="-",
color="white",
lw=line_width + 2.4,
connectionstyle=(
connection_style
),
shrinkA=0,
shrinkB=0,
),
zorder=1.3,
)
ax.annotate(
"",
xy=bar_end,
xytext=start,
arrowprops=dict(
arrowstyle="-",
color=color,
lw=line_width,
linestyle=line_style,
connectionstyle=(
connection_style
),
shrinkA=0,
shrinkB=0,
),
zorder=2.0,
)
terminal_tangent = (
_arc3_end_tangent(
start,
bar_end,
curve_radius,
)
)
bar_perpendicular = np.array(
[
-terminal_tangent[1],
terminal_tangent[0],
]
)
else:
if edge_halo:
ax.plot(
[
start[0],
bar_end[0],
],
[
start[1],
bar_end[1],
],
color="white",
lw=line_width + 2.4,
solid_capstyle="round",
zorder=1.3,
)
ax.plot(
[
start[0],
bar_end[0],
],
[
start[1],
bar_end[1],
],
color=color,
lw=line_width,
linestyle=line_style,
solid_capstyle="round",
zorder=2.0,
)
bar_perpendicular = (
perpendicular
)
# Original larger repression bar.
bar_half_length = (
node_radius
* 0.50
* head_scale
)
point_1 = (
bar_end
- bar_perpendicular
* bar_half_length
)
point_2 = (
bar_end
+ bar_perpendicular
* bar_half_length
)
if edge_halo:
ax.plot(
[
point_1[0],
point_2[0],
],
[
point_1[1],
point_2[1],
],
color="white",
lw=(
line_width * 1.50
+ 2.4
),
solid_capstyle="butt",
zorder=2.2,
)
ax.plot(
[
point_1[0],
point_2[0],
],
[
point_1[1],
point_2[1],
],
color=color,
lw=line_width * 1.50,
solid_capstyle="butt",
zorder=2.6,
)
# ----------------------------------------------------
# Activation
# ----------------------------------------------------
else:
# Original larger arrowhead.
head_length = (
node_radius
* 1.00
* head_scale
)
head_width = (
node_radius
* 0.70
* head_scale
)
arrow_style = (
f"-|>,"
f"head_length="
f"{head_length / node_radius:.2f},"
f"head_width="
f"{head_width / node_radius:.2f}"
)
arrow_end = (
end
- unit * terminal_depth
)
# Undirected edge: plain line, symmetric gaps at both nodes.
if undirected:
arrow_style = "-"
arrow_end = (
end
- unit
* 1.50
* terminal_gap
)
arrow_start = (
start
+ unit
* 1.50
* terminal_gap
)
if edge_halo:
ax.annotate(
"",
xy=arrow_end,
xytext=arrow_start,
arrowprops=dict(
arrowstyle=arrow_style,
color="white",
lw=line_width + 2.4,
connectionstyle=(
connection_style
),
shrinkA=0,
shrinkB=0,
),
zorder=1.3,
)
ax.annotate(
"",
xy=arrow_end,
xytext=arrow_start,
arrowprops=dict(
arrowstyle=arrow_style,
color=color,
lw=line_width,
linestyle=line_style,
connectionstyle=(
connection_style
),
shrinkA=0,
shrinkB=0,
),
zorder=2.0,
)
# --------------------------------------------------------
# Draw nodes
# --------------------------------------------------------
for node in graph.nodes():
ax.add_patch(
plt.Circle(
positions[node],
node_radius * 1.06,
facecolor="white",
edgecolor="none",
zorder=3.5,
)
)
ax.add_patch(
plt.Circle(
positions[node],
node_radius,
facecolor=node_color,
edgecolor="white",
linewidth=1.3,
zorder=4.0,
)
)
# --------------------------------------------------------
# Labels
# --------------------------------------------------------
centroid = np.mean(
[
positions[node]
for node in graph.nodes()
],
axis=0,
)
label_texts = []
for node in graph.nodes():
x, y = positions[node]
outward = (
np.array([x, y])
- centroid
)
outward_norm = float(
np.linalg.norm(outward)
)
if outward_norm > 1e-6:
label_direction = (
outward / outward_norm
)
else:
label_direction = np.array(
[0.0, 1.0]
)
label_x = (
x
+ label_direction[0]
* node_radius
* 2.0
)
label_y = (
y
+ label_direction[1]
* node_radius
* 2.0
)
label_text = ax.text(
label_x,
label_y,
node_labels[node],
ha="center",
va="center",
fontsize=label_font_size,
color=font_color,
fontweight="bold",
zorder=6,
)
label_texts.append(
label_text
)
if fixed_pos is None:
node_x_values = [
positions[node][0]
for node in graph.nodes()
]
node_y_values = [
positions[node][1]
for node in graph.nodes()
]
adjust_text(
label_texts,
x=node_x_values,
y=node_y_values,
ax=ax,
force_points=(0.8, 0.8),
force_text=(0.2, 0.2),
ensure_inside_axes=True,
arrowprops=dict(
arrowstyle="-",
color="#bdbdc4",
lw=0.4,
alpha=0.4,
),
)
# --------------------------------------------------------
# Colorbar
# --------------------------------------------------------
if show_colorbar:
scalar_mappable = (
plt.cm.ScalarMappable(
cmap=colormap,
norm=color_norm,
)
)
scalar_mappable.set_array([])
fig.colorbar(
scalar_mappable,
ax=ax,
orientation="vertical",
fraction=0.03,
pad=0.02,
# [2026-10-23 replaced: label is now an argument]
# label="Edge Weight",
label=colorbar_label,
)
if own_figure:
fig.tight_layout()
return fig, ax
[docs]
def network_from_scenarios(
result,
genes=None,
time="tb",
# [2026-10-24 replaced: edge colour by z_rho_Delta (twin-difference correlation, regulation)]
# weight="z_rho",
# [2026-10-24 replaced: default None = z_dag of the drawn direction with two time points, else z_rho_Delta]
# weight="z_rho_Delta",
weight=None,
solid_scenarios=(2, 3),
dashed_scenarios=(4,),
direction_alpha=0.05,
# [2026-10-24 added: scenario_colors]
scenario_colors=None,
# [2026-10-24 added: z_dag_threshold, two_way_scenarios (pair_directions)]
z_dag_threshold=1.96,
two_way_scenarios=(2,),
):
"""
Turn the output of `infer_with_twinfer` into the inputs of `plot_grn`.
Edge rules
----------
- inference scenario in `solid_scenarios` (default 2 regulation,
3 regulation + heterogeneity): solid edge.
- inference scenario in `dashed_scenarios` (default 4 heterogeneity,
a shared state, not regulation): dashed edge without direction.
- any other scenario (1 nothing, or undefined): no edge.
- direction of a solid edge (only when the table has several time
points), `twinfer.inference.ranked_edges.pair_directions`, the same
rule as the z_dag ranked list:
scenario in `two_way_scenarios` (default 2, regulation): each
direction on its own, |z_dag| > z_dag_threshold (z_dag_xy for
gene_x -> gene_y, z_dag_yx for gene_y -> gene_x); both can be
drawn (mutual regulation, two curved edges);
other solid scenarios (3, regulation + a shared state, which
correlates both directions): one direction, the sign of
delta = |rho_dag_xy| - |rho_dag_yx|, when p_delta < direction_alpha.
Otherwise the edge is drawn without direction.
- edge colour = `weight`. Default (None): "z_dag" when the table has
two time points (the z_dag of the drawn direction; an edge without
direction takes the larger of |z_dag_xy| and |z_dag_yx|; a pair
without cross-time values, see cross_correlation_pairs of
infer_with_twinfer, is coloured by z_rho_Delta), otherwise
"z_rho_Delta". Any other column (without _ta / _tb) is used as is.
- activation arrow vs repression bar of a directed edge = the sign of
its z_dag (the cross-time correlation in the drawn direction),
passed to plot_grn as edge_signs, whatever the colour column.
- `scenario_colors` (optional dict {scenario: colour}, e.g.
{3: "tab:blue", 4: "grey"}): edges of these scenarios are drawn
in that colour instead of the colour of their weight.
Parameters
----------
result
DataFrame returned by infer_with_twinfer, one row per gene pair.
With several time points, keep one (t_a, t_b) pair of rows.
genes
Node order (default: genes in order of appearance in result).
time
"ta" or "tb": which twin group's scenario and weight to use
when the table has several time points (default the later, tb).
scenario_colors
Optional dict {inference scenario: Matplotlib colour}; see the
edge rules.
Returns
-------
plot_kwargs
dict with matrix, node_labels, edge_styles, undirected_edges,
edge_signs (and edge_colors when scenario_colors is given):
call ``plot_grn(**plot_kwargs)`` (other plot_grn options can be
added, e.g. ``plot_grn(**plot_kwargs, colorbar_label="z_rho")``).
edges
DataFrame, one row per drawn edge: source, target, weight (the
colour value), inference_scenario, line_style, directed, z_dag
(of the drawn direction; NaN without cross-time values).
"""
suffix = (
""
if "inference_scenario" in result.columns
else f"_{time}"
)
scenario_column = f"inference_scenario{suffix}"
# [2026-10-24 replaced: weight "z_dag" (default with two time points); z_rho_Delta is its fallback column]
# weight_column = f"{weight}{suffix}"
#
# for column in (scenario_column, weight_column):
has_cross = {"z_dag_xy", "z_dag_yx", "delta", "p_delta"} <= set(result.columns)
if weight is None:
weight = "z_dag" if has_cross else "z_rho_Delta"
if weight == "z_dag" and not has_cross:
raise ValueError("weight z_dag needs a result with two time points (columns z_dag_xy, z_dag_yx, delta, p_delta)")
weight_column = f"z_rho_Delta{suffix}" if weight == "z_dag" else f"{weight}{suffix}"
for column in (scenario_column, weight_column):
if column not in result.columns:
raise ValueError(
f"result has no column {column}"
)
if result.duplicated(["gene_x", "gene_y"]).any():
raise ValueError(
"result has several rows for a gene pair: "
"keep one time pair (t_a, t_b)"
)
if genes is None:
genes = list(
dict.fromkeys(
list(result["gene_x"])
+ list(result["gene_y"])
)
)
index = {
gene: position
for position, gene in enumerate(genes)
}
# [2026-10-24 replaced: has_cross above]
# has_direction = {"delta", "p_delta"} <= set(result.columns)
matrix = np.zeros(
(len(genes), len(genes))
)
edge_styles = {}
edge_colors = {}
# [2026-10-24 added: edge_signs, rho_dag of the drawn direction]
edge_signs = {}
undirected_edges = set()
edges = []
for row in result.itertuples(index=False):
row = row._asdict()
scenario = row[scenario_column]
if pd.isna(scenario):
continue
scenario = int(scenario)
if scenario in solid_scenarios:
line_style = "solid"
elif scenario in dashed_scenarios:
line_style = "--"
else:
continue
if row["gene_x"] not in index or row["gene_y"] not in index:
continue
# [2026-10-24 replaced: directions from pair_directions (both directions in scenario 2), colour z_dag of each direction]
# source = index[row["gene_x"]]
# target = index[row["gene_y"]]
#
# directed = (
# line_style == "solid"
# and has_direction
# and row["p_delta"] < direction_alpha
# )
#
# if directed and row["delta"] < 0:
# source, target = target, source
#
# matrix[source, target] = row[weight_column]
# ...
if has_cross:
directions = pair_directions(row["z_dag_xy"], row["z_dag_yx"], row["delta"], row["p_delta"], scenario,
direction_alpha, z_dag_threshold, two_way_scenarios)
else:
directions = [(True, np.nan, False)]
if line_style != "solid": # dashed (a shared state): one edge without direction
directions = [(directions[0][0], directions[0][1], False)] if len(directions) == 1 else \
[max(((f, z, False) for f, z, _ in directions), key=lambda d: abs(d[1]))]
for forward, z_dag, directed in directions:
source = index[row["gene_x"] if forward else row["gene_y"]]
target = index[row["gene_y"] if forward else row["gene_x"]]
value = z_dag if weight == "z_dag" and not pd.isna(z_dag) else row[weight_column]
matrix[source, target] = value
edge_styles[(source, target)] = line_style
if scenario_colors and scenario in scenario_colors:
edge_colors[(source, target)] = scenario_colors[scenario]
if directed:
edge_signs[(source, target)] = 1 if z_dag >= 0 else -1
else:
undirected_edges.add(
(source, target)
)
edges.append(
dict(
source=genes[source],
target=genes[target],
weight=value,
inference_scenario=scenario,
line_style=line_style,
directed=directed,
z_dag=z_dag,
)
)
plot_kwargs = dict(
matrix=matrix,
node_labels=list(genes),
edge_styles=edge_styles,
undirected_edges=undirected_edges,
edge_signs=edge_signs,
)
if scenario_colors:
plot_kwargs["edge_colors"] = edge_colors
edges = pd.DataFrame(
edges,
columns=[
"source",
"target",
"weight",
"inference_scenario",
"line_style",
"directed",
"z_dag",
],
)
return plot_kwargs, edges
# ============================================================
# Ported from TwINFER_function_scripts/correlation_analysis_helpers.py
# (unchanged) — plot_grn above is not a drop-in replacement (different
# signature: adjacency matrix + labels vs. correlation matrix + gene
# list + edges); disposition vs. plot_grn is an open REORG_CHECKLIST item.
# ============================================================
def shrink_arrow_endpoints(x1, y1, x2, y2, shrink_source=0.2, shrink_target=0.2, lateral_offset=0.0, scaled_shrink_val = 0.25):
dx, dy = x2 - x1, y2 - y1
dist = np.hypot(dx, dy)
if dist == 0:
return (x1, y1), (x2, y2)
scaled_shrink = scaled_shrink_val * dist # 25% of distance, adjust as needed
shrink_source = scaled_shrink
shrink_target = scaled_shrink*0.8
ux, uy = dx / dist, dy / dist
orth_x, orth_y = -uy, ux
x1_off = x1 + shrink_source * ux + lateral_offset * orth_x
y1_off = y1 + shrink_source * uy + lateral_offset * orth_y
x2_off = x2 - shrink_target * ux + lateral_offset * orth_x
y2_off = y2 - shrink_target * uy + lateral_offset * orth_y
return (x1_off, y1_off), (x2_off, y2_off)
def flat_t_head_arrow(start, end, color='red', linewidth=2, rad=0.2, abs_weight=0.5, ax=None):
"""Draw repression arrow with T-head scaled by abs_weight."""
arrow = FancyArrowPatch(
start, end,
connectionstyle=f"arc3,rad={rad}",
arrowstyle='-',
color=color,
linewidth=linewidth,
zorder=1
)
ax.add_patch(arrow)
x1, y1 = start
x2, y2 = end
dx, dy = x2 - x1, y2 - y1
dist = np.hypot(dx, dy)
if dist == 0:
return
ux, uy = dx / dist, dy / dist
tx, ty = x2, y2
px, py = -uy, ux
# Scale T-head size with abs_weight (capped)
t_len = 0.03 + 0.07 * min(abs_weight, 1.0)
t_xs = [tx - px * t_len, tx + px * t_len]
t_ys = [ty - py * t_len, ty + py * t_len]
ax.plot(t_xs, t_ys, color=color, linewidth=linewidth, solid_capstyle='round', zorder=2)
def polygon_layout(gene_list, radius=1.0):
n = len(gene_list)
angles = np.linspace(0, 2 * np.pi, n, endpoint=False)
return {
gene_list[i]: (radius * np.cos(angle), radius * np.sin(angle))
for i, angle in enumerate(angles)
}
def plot_network(correlation_matrix, gene_list, edges, title=None):
DG = nx.DiGraph()
for gene in gene_list:
DG.add_node(gene)
# Add all gene-gene interactions (even if directionless)
# Add all gene-gene interactions (only if both exist in correlation_matrix)
for g1 in gene_list:
for g2 in gene_list:
if g1 == g2:
continue
if g1 in correlation_matrix.index and g2 in correlation_matrix.columns:
if (g1, g2) in edges:
val = correlation_matrix.loc[g1, g2]
if pd.notna(val):
DG.add_edge(g1, g2, weight=val)
# Define node positions and colors
pos = polygon_layout(gene_list, radius=max(2, len(gene_list) / 2)) if len(gene_list) > 2 else {
gene_list[0]: (0, 0),
gene_list[1]: (1, 0)
}
fig, ax = plt.subplots(figsize=(9, 9))
node_colors = np.array([
correlation_matrix.loc[g, g] if g in correlation_matrix.index and g in correlation_matrix.columns else 0
for g in gene_list
])
v = 1
norm = Normalize(vmin=-v, vmax=v)
cmap = make_reds_blues_colormap()
node_rgba = cmap(norm(node_colors))
node_rgba[:, -1] = 0.8
nx.draw_networkx_nodes(DG, pos, ax=ax, node_color=node_rgba, node_size=6000, edgecolors='black', linewidths=1.5)
labels = {node: f"$g{int(node.split('_')[-1])}$" for node in DG.nodes()}
nx.draw_networkx_labels(DG, pos, labels=labels, font_size=16, ax=ax)
# Draw edges
for u, v in DG.edges():
raw_weight = correlation_matrix.loc[u, v]
if pd.isnull(raw_weight) or raw_weight == 0:
continue
x1, y1 = pos[u]
x2, y2 = pos[v]
offset = -0.15 if DG.has_edge(v, u) else 0.0
rad = 0.2 if DG.has_edge(v, u) else 0.0
start, end = shrink_arrow_endpoints(x1, y1, x2, y2, lateral_offset=offset)
color = cmap(norm(raw_weight))
abs_weight = abs(raw_weight)
if raw_weight > 0:
head_length = min(10 + abs_weight * 10, 20)
head_width = min(10 + abs_weight * 10, 20)
arrow_style = f'->,head_length={head_length},head_width={head_width}'
# Positive correlation → activation (arrow)
arrow = FancyArrowPatch(
start, end,
connectionstyle=f'arc3,rad={rad}',
arrowstyle=arrow_style,
mutation_scale=1,
color=color,
linewidth=min(abs_weight * 10, 5.0),
zorder=1
)
ax.add_patch(arrow)
elif raw_weight < 0:
# Negative correlation → repression (flat T-bar)
flat_t_head_arrow(
start, end,
color=color,
linewidth=min(abs_weight * 10, 5.0),
rad=rad,
abs_weight=abs_weight,
ax=ax
)
ax.set_title(title or "Inferred GRN", fontsize=16, fontweight='bold')
ax.axis('off')
if len(pos) > 0:
x_vals, y_vals = zip(*pos.values())
x_range = max(x_vals) - min(x_vals)
y_range = max(y_vals) - min(y_vals)
ax.set_xlim(min(x_vals) - 0.3 * x_range, max(x_vals) + 0.3 * x_range)
ax.set_ylim(min(y_vals) - 0.3 * y_range, max(y_vals) + 0.3 * y_range)
if len(gene_list) == 2:
ax.set_xlim(-0.5, len(gene_list) - 0.5)
ax.set_ylim(-0.5, 0.5)
ax.set_aspect('equal')
plt.subplots_adjust(left=0.05, right=0.95, top=0.9, bottom=0.1)
else:
plt.subplots_adjust(left=0.02, right=0.98, top=0.95, bottom=0.05)
plt.tight_layout()
plt.show()