from typing import Literal
from warnings import warn
from pyiron_snippets.deprecate import deprecate
import matplotlib.pyplot as plt
import matplotlib.patheffects as patheffects
import seaborn as sns
import numpy as np
import pandas as pd
import shapely
from shapely.ops import polylabel
from matplotlib.colors import to_rgba
from .calculate import calc_phase_diagram, cluster, cluster_T_c, _join_phase_unit, _apply_series
from .features import Locus
import landau.poly as poly
__all__ = [
"plot_phase_diagram",
"plot_mu_phase_diagram",
"plot_1d_mu_phase_diagram",
"plot_1d_T_phase_diagram",
]
def _text_with_outline(ax, x, y, s, *, outline_width=3, **kwargs):
"""Draw text with a solid white outline so it stays legible over any fill.
A small reusable wrapper around :meth:`matplotlib.axes.Axes.text` that
strokes the glyphs with a white outline (via matplotlib path effects)
instead of drawing an opaque box behind them, keeping a label readable on
top of coloured regions or tielines. Extra keyword arguments are forwarded
to ``ax.text``.
Returns the created :class:`matplotlib.text.Text`.
"""
kwargs.setdefault(
"path_effects",
[patheffects.withStroke(linewidth=outline_width, foreground="white")],
)
return ax.text(x, y, s, **kwargs)
def _shapely_polygon(coords):
"""Valid shapely polygon from an (N, 2) coordinate array, or ``None``.
A self-intersecting outline (as the TSP-based poly methods can produce) is
repaired with :func:`shapely.make_valid`; if the repair splits it into
several pieces, the largest one is kept. Degenerate (empty or zero-area)
input gives ``None``.
"""
poly = shapely.Polygon(coords)
if not poly.is_valid:
poly = shapely.make_valid(poly)
if isinstance(poly, shapely.MultiPolygon):
poly = max(poly.geoms, key=lambda g: g.area)
if not isinstance(poly, shapely.Polygon) or poly.is_empty or poly.area == 0:
return None
return poly
def _largest_inscribed_circle_center(polygon_xy, ax):
"""Centre of the largest circle inscribable in a polygon, in data units.
Uses shapely's pole of inaccessibility (:func:`shapely.ops.polylabel`),
which lies inside even concave or crescent-shaped regions. Phase-diagram
axes are strongly anisotropic (``c`` spans ~1, ``T`` spans hundreds of
kelvin), so the coordinates are normalised by the axis data-ranges before
the search and mapped back afterwards; this yields the visually – rather
than numerically – largest circle.
Returns ``(x, y)`` in data coordinates, or ``None`` for a degenerate
polygon.
"""
x0, x1 = ax.get_xlim()
y0, y1 = ax.get_ylim()
sx = (x1 - x0) or 1.0
sy = (y1 - y0) or 1.0
coords = np.asarray(polygon_xy, dtype=float)
poly = _shapely_polygon(np.column_stack([coords[:, 0] / sx, coords[:, 1] / sy]))
if poly is None:
return None
point = polylabel(poly, tolerance=1e-3)
return point.x * sx, point.y * sy
def _label_fits(poly_px, text, renderer):
"""Whether ``text``'s rendered bounding box lies inside a pixel-space polygon."""
return poly_px.contains(shapely.box(*text.get_window_extent(renderer).extents))
def _group_overlapping_intervals(intervals, gap=0.0):
"""Group 1-d intervals that overlap, transitively, closer than ``gap``.
Returns lists of indices into ``intervals``; within a group every interval
overlaps (or comes within ``gap`` of) another one, so the members have to be
laid out together along the other axis.
"""
order = sorted(range(len(intervals)), key=lambda i: intervals[i][0])
groups, hi = [], None
for i in order:
lo_i, hi_i = intervals[i]
if hi is None or lo_i > hi + gap:
groups.append([])
hi = hi_i
else:
hi = max(hi, hi_i)
groups[-1].append(i)
return groups
def _add_inline_polygon_labels(ax, polys):
"""Label each phase polygon in place instead of drawing a legend box.
Every polygon is annotated with its phase name (with trailing apostrophes
for repeated stability regions, matching :func:`plot_polygons`), with a
white outline. The text is black: the polygon fill already carries the
phase colour, and black with a white stroke stays legible even over the
pale pastel fills.
Placement tries three positions in turn, keeping the first whose rendered
bounding box fits inside the polygon:
1. horizontal at the centre of the largest inscribed circle,
2. rotated by 90° at the same point — for tall, narrow regions,
3. still rotated, but moved horizontally off the polygon — for line phases
too thin to hold any label. The label sits just right of the polygon
unless that would leave the axes (a terminal line phase at the right
edge), in which case it sits to the left; it is clamped into the axes
both ways, vertically too.
Offset labels are no longer anchored inside their own polygon, so two close
line phases can land on top of each other; a final pass fans labels with
overlapping horizontal extents apart vertically (via :func:`_spread_labels`,
within the axes), so each one slides along its line instead of covering its
neighbour.
Args:
ax: matplotlib Axes the polygons were drawn on.
polys: Series of matplotlib Polygons indexed as in :func:`get_polygons`.
"""
renderer = _get_renderer(ax.figure)
axbb = ax.get_window_extent(renderer)
pad = 0.004 * axbb.width # clearance between an offset label box and its polygon
moved = [] # (text, cx, cy, width, height) of off-polygon labels, in pixels
for key, p in polys.items():
if isinstance(key, tuple):
phase, rep = key
else:
phase, rep = key, 0
center = _largest_inscribed_circle_center(p.get_xy(), ax)
if center is None:
continue
text = _text_with_outline(
ax, center[0], center[1], phase + "'" * rep,
ha="center", va="center", fontsize="small", fontweight="bold",
color="black",
)
poly_px = _shapely_polygon(ax.transData.transform(p.get_xy()))
if poly_px is None or _label_fits(poly_px, text, renderer):
continue
text.set_rotation(90)
if _label_fits(poly_px, text, renderer):
continue
# Too thin even for a rotated label (a line phase): move it beside the
# polygon, keeping the rotation.
bbox = text.get_window_extent(renderer)
half_w = bbox.width / 2
minx, _, maxx, _ = poly_px.bounds
cx = maxx + pad + half_w
if cx + half_w > axbb.x1:
cx = minx - pad - half_w
cx = min(max(cx, axbb.x0 + half_w), axbb.x1 - half_w)
cy = ax.transData.transform(center)[1]
moved.append((text, cx, cy, bbox.width, bbox.height))
# Vertical overlap pass: spread offset labels whose horizontal extents
# collide. _spread_labels also clamps every stack — singletons included —
# into the axes vertically.
inv = ax.transData.inverted()
intervals = [(cx - w / 2, cx + w / 2) for _, cx, _, w, _ in moved]
for group in _group_overlapping_intervals(intervals, gap=pad):
spread = _spread_labels(
[moved[i][2] for i in group],
[moved[i][4] for i in group],
axbb.y0, axbb.y1, gap=pad,
)
for i, cy in zip(group, spread):
text, cx = moved[i][0], moved[i][1]
text.set_position(inv.transform((cx, cy)))
def cluster_phase(df, distance_threshold=0.5): # 0.5 hand-tuned
"""Cluster the stable, single phase regions.
When a (e.g solid solution) phase has multiple disconnected regions of stability, the make_poly and
make_concave_poly functions give wrong results, because they draw a single polygon.
Instead this function adds two new columns `phase_unit` and `phase_id` and the latter will always refer to only a
single connected stability region. `phase_unit` enumerates disconnected regions of one phase.
Args:
df: DataFrame with columns 'phase', 'T', 'c'.
distance_threshold: Passed to :func:`~landau.calculate.cluster_T_c`. Lower values
(e.g. 0.1) split more aggressively and are needed when two stable segments of
the same phase are close in concentration space.
"""
df["phase_unit"] = _apply_series(
df.groupby("phase", group_keys=False),
lambda g: cluster_T_c(g, distance_threshold=distance_threshold),
"phase_unit",
)
df["phase_id"] = _join_phase_unit(df["phase"], df["phase_unit"])
return df
def get_polygons(
df,
poly_method: Literal["concave", "segments", "fasttsp", "tsp", "segment-fasttsp", "segment-tsp"] | poly.AbstractPolyMethod | None = None,
variables: list[str] | None = None,
distance_threshold: float = 0.5, # hand-tuned
**kwargs,
):
"""Turn the stable phase regions in df into polygons.
Args:
df (pandas.DataFrame):
Input data containing columns for the variables and 'phase', 'stable'.
poly_method (str or poly.AbstractPolyMethod, optional):
The method to use for polygon construction.
variables (list of str, optional):
The columns in df to use as coordinates for the polygons. Defaults to ["c", "T"].
distance_threshold (float, optional):
Passed to :func:`cluster_phase`. Lower values split disconnected stable regions
more aggressively. Default is 0.5.
**kwargs:
Passed to poly.handle_poly_method.
Returns:
pandas.Series:
The constructed polygons, indexed by phase and phase_unit.
"""
if variables is None:
variables = ["c", "T"]
df = df.query("stable").copy()
df = cluster_phase(df, distance_threshold=distance_threshold)
if (df.phase_unit == -1).any():
warn("Clustering of phase points failed for some points, dropping them.")
df = df.query("phase_unit>=0")
poly_method = poly.handle_poly_method(poly_method, **kwargs)
return poly_method.apply(df, variables=variables)
def plot_polygons(polys, color_map, ax=None):
"""Plot the given polygons to a matplotlib axis.
Args:
polys (pandas.Series):
The polygons to plot, as returned by get_polygons.
color_map (dict):
Mapping from phase names to colors.
ax (matplotlib.axes.Axes, optional):
The axis to plot on. If None, plt.gca() is used.
"""
if ax is None:
ax = plt.gca()
for i, (phase, p) in enumerate(polys.items()):
with np.errstate(divide="ignore"):
p.zorder = 1 / p.get_extents().size.prod()
if isinstance(phase, tuple):
phase, rep = phase
else:
rep = 0
p.set_color(color_map[phase])
p.set_edgecolor("k")
p.set_label(phase + "'" * rep)
ax.add_patch(p)
def _plot_triplepoints(df, ax=None, variables=None):
"""Mark the three-phase invariants of a phase diagram.
Triple points are tagged :attr:`~landau.features.Locus.TRIPLE` in the
``locus`` column of a refined :func:`~landau.calculate.calc_phase_diagram`
frame; the three coexisting phases share one ``(mu, T)``.
The mark depends on the axes:
* In a concentration-temperature diagram (``variables[0] == "c"``) a triple
point is an isothermal line joining the three coexisting compositions, so
one horizontal line is drawn across the concentration span of each
``(mu, T)`` group.
* In a chemical-potential-temperature diagram (``variables[0] == "mu"``) the
three phases collapse onto a single ``(mu, T)`` point, so a black marker is
drawn there.
Args:
df (pandas.DataFrame):
Phase-diagram frame carrying a ``locus`` column. Unrefined frames
have no triple points and draw nothing.
ax (matplotlib.axes.Axes, optional):
The axis to plot on.
variables (list[str], optional):
The ``[x, y]`` axis variables; defaults to ``["c", "T"]``.
"""
if variables is None:
variables = ["c", "T"]
if ax is None:
ax = plt.gca()
if "locus" not in df.columns:
return
triple = df[df["locus"] == Locus.TRIPLE]
if variables[0] == "c":
for (_mu, T), grp in triple.groupby(["mu", "T"], sort=False)[["c"]]:
ax.hlines(T, grp["c"].min(), grp["c"].max(), color="k", zorder=-2, alpha=0.5, lw=2)
elif variables[0] == "mu":
for (mu, T), _grp in triple.groupby(["mu", "T"], sort=False):
ax.plot(mu, T, marker="o", color="k", zorder=3)
def _set_axis_for(axis_var: str, df_stable, element: str | None, ax) -> None:
"""Configure the x-axis of a phase diagram based on the axis variable.
Args:
axis_var: The x-axis variable name, either ``"c"`` (concentration) or ``"mu"``
(chemical potential).
df_stable: Stable-phase rows of the diagram DataFrame (used only for ``"mu"``
to determine finite x-limits).
element: Optional element symbol used to build the axis label.
ax: The :class:`matplotlib.axes.Axes` to configure.
Raises:
ValueError: If ``axis_var`` is not ``"c"`` or ``"mu"``.
"""
if axis_var == "c":
ax.set_xlim(0, 1)
if element is not None:
ax.set_xlabel(rf"$c_\mathrm{{{element}}}$")
else:
ax.set_xlabel("$c$")
elif axis_var == "mu":
mus = df_stable["mu"].unique()
mus = mus[np.isfinite(mus)]
if len(mus) > 0:
ax.set_xlim(mus.min(), mus.max())
if element is not None:
ax.set_xlabel(rf"$\Delta\mu_\mathrm{{{element}}}$ [eV]")
else:
ax.set_xlabel(r"$\Delta\mu$ [eV]")
else:
raise ValueError(
f"Unknown coordinate system: variables[0]={axis_var!r}. Expected 'c' or 'mu'."
)
def _plot_phase_diagram(
df,
alpha=0.1,
element=None,
min_c_width=1e-2,
color_override: dict[str, str] = {},
triplepoints=False,
poly_method: Literal["concave", "segments", "fasttsp", "tsp", "segment-fasttsp", "segment-tsp"] | poly.AbstractPolyMethod | None = None,
variables: list[str] | None = None,
inline_legend=True,
legend=True,
ax=None,
):
if variables is None:
variables = ["c", "T"]
if ax is None:
ax = plt.gca()
df_stable = df.query("stable")
color_map = get_phase_colors(df_stable.phase.unique(), color_override)
polys = get_polygons(df, poly_method=poly_method, variables=variables, min_c_width=min_c_width, alpha=alpha)
plot_polygons(polys, color_map, ax=ax)
if triplepoints:
_plot_triplepoints(df, ax=ax, variables=variables)
_set_axis_for(variables[0], df_stable, element, ax)
ax.set_ylim(df_stable["T"].min(), df_stable["T"].max())
# Inline labels need the final axis limits to place each label at the centre
# of its polygon, so this runs after the limits above are set.
if legend:
if inline_legend:
_add_inline_polygon_labels(ax, polys)
else:
ax.legend(ncols=2)
ax.set_ylabel("$T$ [K]")
[docs]
@deprecate(
alpha="Pass a poly method from landau.poly to poly_method",
min_c_width="Pass a poly method from landau.poly to poly_method",
tielines="Use triplepoints instead",
)
def plot_phase_diagram(
df,
alpha=0.1,
element=None,
min_c_width=1e-2,
color_override: dict[str, str] = {},
triplepoints=False,
poly_method: Literal["concave", "segments", "fasttsp", "tsp", "segment-fasttsp", "segment-tsp"] | poly.AbstractPolyMethod | None = None,
variables: list[str] | None = None,
inline_legend=True,
legend=True,
ax=None,
tielines=None,
):
if tielines is not None:
triplepoints = tielines
return _plot_phase_diagram(
df,
alpha=alpha,
element=element,
min_c_width=min_c_width,
color_override=color_override,
triplepoints=triplepoints,
poly_method=poly_method,
variables=variables,
inline_legend=inline_legend,
legend=legend,
ax=ax,
)
def get_phase_colors(phase_names, override: dict[str, str] | None = None):
if override is None:
override = {}
# the default map
color_map = dict(zip(phase_names, sns.palettes.SEABORN_PALETTES["pastel"]))
# disregard overriden phases that are not present
override = {p: c for p, c in override.items() if p in color_map}
# if the override uses the same colors as the default map, multiple phases
# would be mapped to the same color; so instead let's update the color map of phases that would
# use the same color as a phase in the override to use the default colors of the overriden phases
# instead
duplicates_map = {c: color_map[o] for o, c in override.items()}
diff = {k: duplicates_map[c] for k, c in color_map.items() if c in duplicates_map}
color_map.update(diff | override)
return color_map
[docs]
@deprecate(alpha="Pass a poly method from landau.poly to poly_method")
def plot_mu_phase_diagram(
df,
alpha=0.1,
element=None,
color_override: dict[str, str] = {},
triplepoints=False,
poly_method: Literal["concave", "segments", "fasttsp", "tsp", "segment-fasttsp", "segment-tsp"] | poly.AbstractPolyMethod | None = None,
inline_legend=True,
legend=True,
ax=None,
):
return _plot_phase_diagram(
df,
alpha=alpha,
element=element,
color_override=color_override,
triplepoints=triplepoints,
poly_method=poly_method,
variables=["mu", "T"],
inline_legend=inline_legend,
legend=legend,
ax=ax,
)
def _assign_segment_ids(df: pd.DataFrame, scan_col: str) -> pd.Series:
"""Assign contiguous-segment IDs for a 1d scan along ``scan_col``.
Ordering the points by ``scan_col``, a phase begins a new segment whenever its
``stable`` flag flips. This is threshold-free and independent of grid density
or of how derived quantities (concentration, ...) vary along the cut, so two
disjoint metastable branches of one phase – e.g. a middle phase that is
unstable both below and above its stable window – are never joined into one
line. Pass the result as ``units=`` to seaborn so it draws each segment
separately.
The cut axis is a parameter, so the helper extends unchanged to generalised
1d diagrams along arbitrary T / mu / ... cuts: pass whichever column orders
the points along the cut.
Assumes every phase is sampled at every scan point – as produced by
:func:`~landau.calculate.calc_phase_diagram` with ``keep_unstable=True`` – so
a stability flip is the only way a phase's run along the cut can break.
Args:
df: DataFrame with 'phase', 'stable', and ``scan_col`` columns.
scan_col: Column ordering the points along the cut ('mu', 'T', ...).
Returns:
String Series aligned with ``df.index``, one unique value per segment.
"""
ordered = df.sort_values(scan_col)
seg = ordered.groupby("phase", group_keys=False).apply(
lambda g: (g["stable"] != g["stable"].shift()).cumsum().to_frame("_seg"),
include_groups=False,
)["_seg"]
return (
ordered["phase"].astype(str)
+ "_"
+ ordered["stable"].astype(int).astype(str)
+ "_"
+ seg.astype(str)
).reindex(df.index)
def _bridge_unstable_segments(df: pd.DataFrame, scan_col: str) -> pd.DataFrame:
"""Extend each unstable branch up to the exact transition point.
A refined transition (``border``) row sits exactly at a stability flip but is
marked ``stable``, so it anchors the solid line of both coexisting phases.
The phase that turns metastable across the flip only resumes at the next
sampled point, leaving a gap between the transition and the start of its
dashed branch – the wider the gap, the coarser the grid. For every such
border row this duplicates it as an unstable point tagged with the adjacent
unstable segment's ``_seg_id`` (the metastable side of the flip), so the
dashed branch is drawn right up to the transition.
Operates on the ``_seg_id``-tagged frame and returns it with the bridge rows
appended; the originals are untouched, so the solid line still reaches the
same point. A no-op when no ``border`` column is present.
Args:
df: DataFrame with 'phase', 'stable', '_seg_id', ``scan_col`` and
optionally 'border' columns.
scan_col: Column ordering the points along the cut ('mu', 'T', ...).
Returns:
``df`` with one duplicated unstable row per (border, adjacent-unstable)
pair appended, reindexed.
"""
if "border" not in df.columns:
return df
src_idx, seg_ids = [], []
for _, g in df.groupby("phase", sort=False):
g = g.sort_values(scan_col)
idx = g.index.to_numpy()
stable = g["stable"].to_numpy()
border = g["border"].to_numpy()
seg = g["_seg_id"].to_numpy()
for i in range(len(g)):
if not (border[i] and stable[i]):
continue
# The flip is stable on one side, metastable on the other; bridge to
# whichever neighbour(s) along the cut are unstable.
for j in (i - 1, i + 1):
if 0 <= j < len(g) and not stable[j]:
src_idx.append(idx[i])
seg_ids.append(seg[j])
if not src_idx:
return df
bridges = df.loc[src_idx].copy()
bridges["stable"] = False
bridges["_seg_id"] = seg_ids
return pd.concat([df, bridges], ignore_index=True)
def _subtract_reference_phase(df, scan_col, reference_phase):
"""Subtract reference phase's phi from all phases along scan_col.
The reference phase must be sampled at every grid point of the scan. Otherwise
the ``np.interp`` below silently fabricates the missing reference values --
clamping to a constant beyond the reference's range, or bridging an interior
gap with a chord -- and corrupts every phase's curve. A frame keeping only
stable rows triggers this: the reference is then present only within its own
stability window(s), so recompute the diagram with ``keep_unstable=True``.
"""
if reference_phase not in df["phase"].values:
raise ValueError(f"reference_phase {reference_phase!r} not found in data")
ref = df.loc[df["phase"] == reference_phase, [scan_col, "phi"]].sort_values(scan_col)
# Refined border rows are excluded: the reference legitimately has no row there
# (refinement adds rows only for the two transitioning phases), and np.interp
# interpolates them from the dense grid neighbours on either side.
grid = df.loc[~df["border"], scan_col] if "border" in df.columns else df[scan_col]
grid = np.unique(grid)
missing = np.setdiff1d(grid, ref[scan_col].to_numpy())
if missing.size:
raise ValueError(
f"reference_phase {reference_phase!r} is not sampled across the full "
f"{scan_col} range ({missing.size} of {grid.size} grid points have no "
f"reference row); recompute the diagram with keep_unstable=True."
)
df["phi"] = df["phi"] - np.interp(df[scan_col], ref[scan_col], ref["phi"])
return df
def _bold_math(label: str) -> str:
"""Wrap mathtext segments of ``label`` in ``\\mathbf`` so subscripts render bold.
matplotlib's ``fontweight="bold"`` only affects the regular-text portions of a
string; content inside ``$...$`` keeps the regular math weight. Wrapping each
math segment in ``\\mathbf{...}`` bolds it to match while leaving arbitrary
LaTeX inside renderable. A malformed string (odd number of ``$``) is returned
unchanged.
"""
parts = label.split("$")
if len(parts) % 2 == 0: # unbalanced '$' -> leave untouched
return label
for i in range(1, len(parts), 2):
if parts[i]:
parts[i] = r"\mathbf{" + parts[i] + "}"
return "$".join(parts)
def _get_renderer(fig):
"""Return a renderer for *fig*, drawing the canvas first so text can be measured.
``Text.get_window_extent`` needs a renderer to report the rendered pixel size of a
label. ``Figure.canvas.get_renderer`` exists on the Agg backend; other backends
expose it via the private ``Figure._get_renderer`` (matplotlib >= 3.6).
"""
fig.canvas.draw()
if hasattr(fig.canvas, "get_renderer"):
return fig.canvas.get_renderer()
return fig._get_renderer()
def _spread_labels(centers, heights, lo, hi, gap=0.0):
"""Nudge label centers apart so their vertical extents never overlap.
Each label *i* occupies ``[center - heights[i]/2, center + heights[i]/2]``. The
returned centers keep the input order's stacking, sit within ``[lo, hi]`` where the
stack fits, and move as little as possible from the requested ``centers``. Works in
whatever 1-d coordinate the caller passes (display pixels are convenient because a
rendered text height is constant there regardless of the data scale).
Args:
centers: Desired center coordinate of each label.
heights: Full extent of each label in the same units as ``centers``.
lo, hi: Lower / upper bounds the stack should stay within.
gap: Extra clearance to keep between adjacent labels.
Returns:
List of adjusted centers, one per input label, in the input order.
"""
n = len(centers)
if n == 0:
return []
order = sorted(range(n), key=lambda i: centers[i])
adj = list(centers)
# Upward pass: push each label up just enough to clear the one below it.
prev_top = lo
for i in order:
half = heights[i] / 2
c = max(centers[i], prev_top + half)
adj[i] = c
prev_top = c + half + gap
# Downward pass: pull labels down to respect the top bound while staying separated.
next_bot = hi
for i in reversed(order):
half = heights[i] / 2
c = min(adj[i], next_bot - half)
adj[i] = c
next_bot = c - half - gap
return adj
def _phase_visible_in_band(phi, lo, hi):
"""Whether a line with potentials ``phi`` (in scan order) shows inside ``[lo, hi]``.
Visible if any sampled point falls within the band, or if two consecutive points
straddle it (the connecting segment crosses the whole window). Used to drop the
label of a phase whose line an applied ``ylim`` pushes entirely out of view.
"""
phi = np.asarray(phi, dtype=float)
phi = phi[np.isfinite(phi)]
if phi.size == 0:
return False
if np.any((phi >= lo) & (phi <= hi)):
return True
below, above = phi < lo, phi > hi
crosses = (below[:-1] & above[1:]) | (above[:-1] & below[1:])
return bool(np.any(crosses))
def _place_side_labels(ax, df, scan_col, phase_colors, top_texts=()):
"""Label every phase at the end of its line, reserving room adaptively.
A phase is labelled at its right-hand line end by default. The horizontal space
reserved for the stack is derived from the widest *rendered* label (measured in
pixels), not from any string-length heuristic, so the axis is widened by exactly as
much as the labels need. Within a stack the labels are spread vertically so they do
not overlap, clamped to the current y-limits.
If the current y-limit would clip a label off the top (its right-end value lies above
the visible window), that phase is instead labelled at its left-hand line end on a
mirrored left-hand stack, which is reserved and laid out by the same rules. A phase
whose line the y-limit pushes entirely out of view is not labelled at all.
``top_texts`` are the bold top stable-phase labels; both stacks are capped just below
their rendered bottom so they never intrude into that band.
"""
fig = ax.figure
x_min, x_max = df[scan_col].min(), df[scan_col].max()
span0 = x_max - x_min
# Settle autoscaled limits (and obtain a renderer) before deciding which side each
# label belongs to or measuring any text.
renderer = _get_renderer(fig)
lo_d, hi_d = sorted(ax.get_ylim())
axbb = ax.get_window_extent(renderer)
# Ceiling for both stacks: just below the bold top-label band so they don't collide.
hi_px = axbb.y1
if len(top_texts):
top_bottom = min(t.get_window_extent(renderer).y0 for t in top_texts)
hi_px = max(min(hi_px, top_bottom - 0.01 * axbb.height), axbb.y0)
# Split phases: those visible at the right end label on the right; those whose
# right end is above the window label at their left end instead. A phase whose
# whole line is pushed out of view by the y-limit gets no label at all.
right, left = [], []
for phase, group in df.groupby("phase"):
g = group.sort_values(scan_col)
phi = g["phi"].to_numpy()
if not _phase_visible_in_band(phi, lo_d, hi_d):
continue
right_y = g["phi"].iloc[-1]
if right_y > hi_d:
left.append((phase, g["phi"].iloc[0]))
else:
right.append((phase, right_y))
def make_texts(entries, ha):
return [
ax.text(
x_min, y, phase, transform=ax.transData,
ha=ha, va="center", fontsize="small",
color=phase_colors.get(phase, "black"), clip_on=True,
)
for phase, y in entries
]
right_texts = make_texts(right, "left")
left_texts = make_texts(left, "right")
def measure(texts):
if not texts:
return 0.0, []
ext = [t.get_window_extent(renderer) for t in texts]
return max(e.width for e in ext) / axbb.width, [e.height for e in ext]
gap_frac, margin_frac = 0.02, 0.01
f_right, h_right = measure(right_texts)
f_left, h_left = measure(left_texts)
reserve_r = (gap_frac + f_right + margin_frac) if right_texts else 0.0
reserve_l = (gap_frac + f_left + margin_frac) if left_texts else 0.0
# Widen the axis so the line data occupies the middle (1 - reserve_l - reserve_r) of
# it and each stack sits in its reserved strip.
denom = max(1.0 - reserve_l - reserve_r, 0.2)
total_span = span0 / denom
x0 = x_min - reserve_l * total_span
ax.set_xlim(x0, x0 + total_span)
def place(texts, entries, heights, anchor_frac):
if not texts:
return
anchor_x = x0 + anchor_frac * total_span
# Spread the true (unclamped) line-terminal pixels: _spread_labels bounds the
# result to the axes box itself, so the labels stay in the window while keeping
# the terminals' vertical order. Pre-clamping would tie targets that fall
# outside the window and collapse that order.
target_px = [ax.transData.transform((anchor_x, y))[1] for _, y in entries]
placed_px = _spread_labels(target_px, heights, axbb.y0, hi_px)
inv = ax.transData.inverted()
for t, py in zip(texts, placed_px):
y_d = inv.transform((axbb.x0, py))[1]
t.set_position((anchor_x, y_d))
place(right_texts, right, h_right, 1.0 - margin_frac - f_right)
place(left_texts, left, h_left, f_left + margin_frac)
def _place_transition_labels(ax, positions, labels, *, side, **text_kw):
"""Label transition lines with vertical text, spread to avoid overlaps.
Each transition at ``positions[i]`` is annotated with ``labels[i]`` as rotated
(vertical) text just inside the bottom of the axes. Each label is first offset
to one side of its line, then the whole row is fanned apart along x by the same
:func:`_spread_labels` routine that stacks the side labels vertically -- this
removes every mutual overlap while preserving the labels' left-to-right order.
A final best-effort pass shifts any label that still straddles a dotted
transition line clear of it, but only within the slack its neighbours leave, so
it never re-introduces an overlap; lines packed closer than a label is wide
cannot all be avoided and are left as they fall.
``side`` ('left' or 'right') biases the initial offset of each label relative
to its own line.
Returns the created :class:`~matplotlib.text.Text` artists.
"""
if len(positions) == 0:
return []
order = np.argsort(positions)
positions = [positions[i] for i in order]
labels = [labels[i] for i in order]
renderer = _get_renderer(ax.figure)
axbb = ax.get_window_extent(renderer)
# y in axes fraction (blended transform) so a zoomed ylim can't fling the text
# out of view; only the x position is spread.
xform = ax.get_xaxis_transform()
texts = [
_text_with_outline(
ax, x, 0.02, s, transform=xform,
rotation="vertical", ha="center", va="bottom", fontsize="small", zorder=100, **text_kw,
)
for x, s in zip(positions, labels)
]
widths = [t.get_window_extent(renderer).width for t in texts]
# Spreading is done in display pixels (a rendered width is constant there), then
# mapped back to the data x the blended transform expects.
lines_px = [ax.transData.transform((x, 0.0))[0] for x in positions]
pad = 0.004 * axbb.width # clearance kept between a label box and a line
half = [w / 2 for w in widths]
offset = -1 if side == "left" else 1
# Offset each label to one side of its line, then fan them apart (order- and
# bound-preserving) so no two boxes overlap.
centers = [lines_px[i] + offset * (half[i] + pad) for i in range(len(texts))]
centers = _spread_labels(centers, widths, axbb.x0, axbb.x1, gap=pad)
# Best-effort: shift any label that still covers a transition line off it,
# clamped to the slack between its neighbours (and the axes) so the no-overlap
# guarantee from the spread survives.
def covered(c, i):
return sum(c - half[i] - pad < L < c + half[i] + pad for L in lines_px)
for i, c in enumerate(centers):
lo = axbb.x0 + half[i] if i == 0 else centers[i - 1] + half[i - 1] + half[i] + pad
hi = axbb.x1 - half[i] if i == len(centers) - 1 else centers[i + 1] - half[i + 1] - half[i] - pad
if lo > hi or covered(c, i) == 0:
continue
# Candidate slots sit the box just left or right of each line; clamp to the
# neighbour slack and keep the reachable one that covers the fewest lines.
slots = [L + s * (half[i] + pad) for L in lines_px for s in (-1, 1)]
cands = [c] + [min(max(s, lo), hi) for s in slots]
centers[i] = min(cands, key=lambda cc: (covered(cc, i), abs(cc - c)))
inv = ax.transData.inverted()
for t, c in zip(texts, centers):
t.set_position((inv.transform((c, 0.0))[0], 0.02))
return texts
def _add_1d_phase_legend(ax, df, scan_col, top_labels=True, side_labels=True, ylim=None):
"""Annotate a 1d phase diagram with inline phase labels.
top_labels
Ticks on the top spine mark each transition boundary (positions where
``border == True``) and the stable phase name is placed near the top of
the axis, centered between adjacent boundaries.
side_labels
The default seaborn legend is removed and, when unstable lines are
present, every phase is annotated by name at the right end of its final
line segment inside the axis.
Args:
ax: matplotlib Axes with a seaborn lineplot already rendered.
df: DataFrame with 'phase', 'stable', ``scan_col``, 'phi', and
(if available) 'border' columns.
scan_col: Scan-axis column ('mu' or 'T').
top_labels: If True, add the top-spine ticks and stable-phase labels.
side_labels: If True, remove the default seaborn legend and add the
right-end side labels.
ylim: If given, applied via ``ax.set_ylim`` before the labels are placed,
so the side labels are clamped to (and spread within) this window.
A scalar is treated as ``(None, ylim)`` — an upper bound only.
"""
if ylim is not None:
if np.isscalar(ylim):
ylim = (None, ylim)
ax.set_ylim(ylim)
# Extract phase → color from the seaborn legend (used by both label sets).
phase_colors = {}
legend = ax.get_legend()
if legend is not None:
white = to_rgba("w")
gray2 = to_rgba(".2")
for handle, text in zip(legend.legend_handles, legend.texts):
try:
c = handle.get_color()
except AttributeError:
continue
if to_rgba(c) in (white, gray2):
continue
phase_colors[text.get_text()] = c
# The side labels replace the default seaborn legend.
if side_labels:
legend.remove()
# Store so tests (or user code) can still access the color map.
ax._landau_phase_colors = phase_colors
top_texts = []
if top_labels and "border" in df.columns:
# Interior transition positions only.
x_min = df[scan_col].min()
x_max = df[scan_col].max()
transitions = sorted(
t for t in df.loc[df["border"], scan_col].unique()
if x_min < t < x_max
)
boundaries = [x_min] + transitions + [x_max]
# Top-spine ticks at transition boundaries (no tick labels; just marks).
ax2 = ax.secondary_xaxis("top")
ax2.set_ticks(transitions)
ax2.set_xticklabels([])
ax2.tick_params(direction="in", length=6)
# Phase-name labels near the top of the axis, centered between boundaries.
# get_xaxis_transform(): x in data coords, y in axes [0,1] fraction.
xform = ax.get_xaxis_transform()
for lo, hi in zip(boundaries[:-1], boundaries[1:]):
mid = (lo + hi) / 2
# Strict inequalities exclude the border rows themselves, which can have
# two phases marked stable simultaneously (the transition point).
mask = (df[scan_col] > lo) & (df[scan_col] < hi) & df["stable"]
stable_phases = df.loc[mask, "phase"].unique()
if len(stable_phases) != 1:
# A stable window narrower than the grid spacing has no sample
# strictly inside its two borders; the phase marked stable on
# both enclosing rows is the one stable throughout.
at_lo = set(df.loc[(df[scan_col] == lo) & df["stable"], "phase"])
at_hi = set(df.loc[(df[scan_col] == hi) & df["stable"], "phase"])
stable_phases = sorted(at_lo & at_hi)
if len(stable_phases) != 1:
raise RuntimeError(
f"expected exactly one stable phase in [{lo}, {hi}], "
f"got {list(stable_phases)}"
)
phase = stable_phases[0]
# White outline keeps the bold label legible on top of tielines.
top_texts.append(_text_with_outline(
ax, mid, 0.97, _bold_math(phase),
transform=xform,
ha="center", va="top", fontsize="small", fontweight="bold",
color=phase_colors.get(phase, "black"),
))
if side_labels and not df["stable"].all():
_place_side_labels(ax, df, scan_col, phase_colors, top_texts)
[docs]
def plot_1d_mu_phase_diagram(
df,
ax=None,
show=True,
mark_transitions=True,
reference_phase=None,
top_labels=True,
side_labels=True,
ylim=None):
"""
Plot a one dimensional isothermal phase diagram of the semi-grandcanonical
potential as function of the chemical potential difference.
Args:
df (pandas.DataFrame):
Input data containing columns for chemical potential difference ('mu'),
semi-grandcanonical potential ('phi'), phase name ('phase'), stability
('stable'), and optionally a 'border' column indicating phase transition.
ax (matplotlib.axes.Axes, optional):
Existing matplotlib Axes to plot on. If None, a new figure and axes are created.
mark_transitions (bool, optional):
If True, all transition temperatures are marked on the plot. Defaults to True.
reference_phase (str, optional):
If given, subtract this phase's potential from all other phases before
plotting so that the reference phase lies at zero throughout.
top_labels (bool, optional):
If True, label the stable phase of each segment near the top of the
axis. Defaults to True.
side_labels (bool, optional):
If True, remove the default seaborn legend and label every phase at
the right end of its line instead. Defaults to True.
ylim (tuple or float, optional):
If given, applied like :func:`matplotlib.pyplot.ylim`. A scalar is
treated as ``(None, ylim)`` (upper bound only). It also bounds the
side-label stack: a label whose line end is above the window is moved
to a mirrored stack on the left, and a phase whose line is pushed
entirely out of view is not labelled.
Returns:
matplotlib.axes.Axes:
The Axes object with the phase diagram plot.
"""
if len(df['T'].unique()) > 1:
raise ValueError("data contains more than one temperature!")
if ax is None:
fig, ax = plt.subplots()
df = df.sort_values("mu").copy()
if reference_phase is not None:
df = _subtract_reference_phase(df, "mu", reference_phase)
df["_seg_id"] = _assign_segment_ids(df, scan_col="mu")
sns.lineplot(
data=_bridge_unstable_segments(df, scan_col="mu"),
x='mu', y='phi',
hue='phase', hue_order=sorted(df.phase.unique()),
style='stable', style_order=[True, False],
units='_seg_id', estimator=None, errorbar=None,
ax=ax,
)
_add_1d_phase_legend(ax, df, scan_col="mu", top_labels=top_labels, side_labels=side_labels, ylim=ylim)
if 'border' not in df.columns:
return ax
if mark_transitions:
# The marker dot stays in data coords and is simply clipped if the crossing
# lies outside the window; the labels are placed (and spread) by
# _place_transition_labels.
positions, labels = [], []
for mt, dd in df.query("mu.min()<mu<mu.max() and border").groupby("mu"):
ft = dd['phi'].iloc[0]
ax.axvline(mt, color='k', linestyle='dotted', alpha=.5)
ax.scatter(mt, ft, marker='o', c='k', zorder=10)
positions.append(mt)
labels.append(rf"$\Delta\mu = {mt:.03f}\,\mathrm{{eV}}$")
_place_transition_labels(ax, positions, labels, side="left")
ax.set_xlabel("Chemical Potential Difference [eV]")
ylabel = "Semi-grandcanonical Potential [eV/atom]"
if reference_phase is not None:
ylabel = f"Semi-grandcanonical Potential\nrelative to {reference_phase} [eV/atom]"
ax.set_ylabel(ylabel)
return ax
[docs]
def plot_1d_T_phase_diagram(
df,
ax=None,
mark_transitions=True,
show=True,
reference_phase=None,
top_labels=True,
side_labels=True,
ylim=None,
):
"""
Plots a one-dimensional equipotential phase diagram as a function of temperature.
Args:
df (pandas.DataFrame):
Input data containing columns for temperature ('T'), semi-grandcanonical
potential ('phi'), phase name ('phase'), and optionally a 'border' column
indicating phase transition.
ax (matplotlib.axes.Axes, optional):
Existing matplotlib Axes to plot on. If None, a new figure and axes are created.
mark_transitions (bool, optional):
If True, all transition temperatures are marked on the plot. Defaults to True.
reference_phase (str, optional):
If given, subtract this phase's potential from all other phases before
plotting so that the reference phase lies at zero throughout.
top_labels (bool, optional):
If True, label the stable phase of each segment near the top of the
axis. Defaults to True.
side_labels (bool, optional):
If True, remove the default seaborn legend and label every phase at
the right end of its line instead. Defaults to True.
ylim (tuple or float, optional):
If given, applied like :func:`matplotlib.pyplot.ylim`. A scalar is
treated as ``(None, ylim)`` (upper bound only). It also bounds the
side-label stack: a label whose line end is above the window is moved
to a mirrored stack on the left, and a phase whose line is pushed
entirely out of view is not labelled.
Returns:
matplotlib.axes.Axes:
The Axes object with the phase diagram plot.
"""
if len(df.mu.unique()) > 1:
raise ValueError("Data contains more than one chemical potential!")
if ax is None:
fig, ax = plt.subplots()
df = df.copy()
if reference_phase is not None:
df = _subtract_reference_phase(df, "T", reference_phase)
df["_seg_id"] = _assign_segment_ids(df, scan_col="T")
sns.lineplot(
data=_bridge_unstable_segments(df, scan_col="T"),
x='T', y='phi',
hue='phase', hue_order=sorted(df.phase.unique()),
style='stable', style_order=[True, False],
units='_seg_id', estimator=None, errorbar=None,
ax=ax,
)
_add_1d_phase_legend(ax, df, scan_col="T", top_labels=top_labels, side_labels=side_labels, ylim=ylim)
if 'border' not in df.columns:
return ax
if mark_transitions:
# The marker dot stays in data coords and is simply clipped if the crossing
# lies outside the window; the labels are placed (and spread) by
# _place_transition_labels.
positions, labels = [], []
for Tt, dd in df.query("T.min()<T<T.max() and border").groupby("T"):
ft = dd['phi'].iloc[0]
ax.axvline(Tt, color='k', linestyle='dotted', alpha=.5)
ax.scatter(Tt, ft, marker='o', c='k', zorder=10)
positions.append(Tt)
labels.append(rf"$T = {Tt:.0f}\,\mathrm{{K}}$")
_place_transition_labels(ax, positions, labels, side="right")
ax.set_xlabel("Temperature [K]")
ylabel = "Semi-grandcanonical potential [eV/atom]"
if reference_phase is not None:
ylabel = f"Semi-grandcanonical potential\nrelative to {reference_phase} [eV/atom]"
ax.set_ylabel(ylabel)
return ax
# ---------------------------------------------------------------------------
# Excess free energy plot
# ---------------------------------------------------------------------------
def _curve_obstacles(ax):
"""Drawn curves and scatter markers of *ax* as one shapely geometry in pixels.
Phase curves become :class:`shapely.LineString`\\ s and scatter dots (line-phase
markers, hull vertices) become :class:`shapely.Point`\\ s buffered by their
marker radius, all in display coordinates so a label's pixel bounding box can be
tested against them directly. The horizontal reference line (a blended-transform
``axhline``) is skipped so labels are not pushed off the zero line. Returns the
unioned geometry, or ``None`` when nothing is drawn yet.
"""
geoms = []
for line in ax.lines:
if line.get_transform() is not ax.transData: # skip refline / blended artists
continue
disp = line.get_transform().transform(line.get_xydata())
disp = disp[np.isfinite(disp).all(axis=1)]
if len(disp) >= 2:
geoms.append(shapely.LineString(disp))
dpi = ax.figure.dpi
for coll in ax.collections:
offsets = np.asarray(coll.get_offsets(), dtype=float)
if offsets.size == 0:
continue
disp = coll.get_offset_transform().transform(offsets)
sizes = coll.get_sizes()
for i, (px, py) in enumerate(disp):
if not (np.isfinite(px) and np.isfinite(py)):
continue
s = sizes[i % len(sizes)] if len(sizes) else 36.0
radius = max((np.sqrt(s) / 2.0) * dpi / 72.0, 1.0) # points^2 -> pixel radius
geoms.append(shapely.Point(px, py).buffer(radius))
return shapely.union_all(geoms) if geoms else None
def _add_inline_curve_labels(ax, entries):
"""Label curves/points just off the data line instead of via a legend box.
Each entry ``(label, x, y, color, side)`` is anchored at ``(x, y)`` on the
curve (or dot) and drawn one label-height clear of it rather than on top of it:
``side="above"`` lifts the label into the open (convex) side of a free-energy
curve, ``side="below"`` drops it under a lower-hull line-phase dot. Labels are
white-outlined in the phase colour with mathtext subscripts bolded by
:func:`_bold_math`, then nudged apart along ``y`` in display pixels by
:func:`_spread_labels` so they never overlap, while ``x`` stays anchored to the
curve. A final pass tests each label's pixel box (via shapely) against the
drawn curves, scatter markers (:func:`_curve_obstacles`) and the
already-placed labels, pushing it further in its ``side`` direction until it no
longer overlaps any of them or it reaches the axes edge. Mirrors the inline
labelling used by the 2d (:func:`_add_inline_polygon_labels`) and 1d
(:func:`_place_side_labels`) plots.
Args:
ax: matplotlib Axes the curves were drawn on.
entries: iterable of ``(label, x, y, color, side)`` anchors, one per phase,
with ``side`` one of ``"above"`` / ``"below"``.
Returns the created :class:`~matplotlib.text.Text` artists.
"""
entries = list(entries)
if not entries:
return []
renderer = _get_renderer(ax.figure)
axbb = ax.get_window_extent(renderer)
texts = [
_text_with_outline(
ax, x, y, _bold_math(label),
ha="center", va="center", fontsize="small", fontweight="bold",
color=color, zorder=10,
)
for label, x, y, color, _side in entries
]
heights = [t.get_window_extent(renderer).height for t in texts]
pad = 0.01 * axbb.height # clearance kept between a label box and the line
target_px = [
ax.transData.transform((x, y))[1] + (h / 2 + pad) * (1 if side == "above" else -1)
for (_label, x, y, _color, side), h in zip(entries, heights)
]
placed_px = _spread_labels(target_px, heights, axbb.y0, axbb.y1)
inv = ax.transData.inverted()
for t, (_label, x, *_rest), py in zip(texts, entries, placed_px):
t.set_position((x, inv.transform((axbb.x0, py))[1]))
# Final pass: push any label still overlapping a curve, a marker or an earlier
# label further out until its pixel box clears them. The box only translates in
# y, so candidates are tested by shifting its bounds rather than re-measuring text.
obstacles = _curve_obstacles(ax)
step = 2.0 # px
clearance = 0.012 * axbb.height # small gap kept between a label box and a line
for t, (_label, _x, _y, _color, side) in zip(texts, entries):
e = t.get_window_extent(renderer)
half = e.height / 2
base_cy = (e.y0 + e.y1) / 2
lo_c, hi_c = axbb.y0 + half, axbb.y1 - half
# Label box at vertical centre ``cy``, inflated by ``clearance`` so a label
# stops a hair short of a curve, marker or neighbour rather than flush against
# it. Only y moves, so the x bounds stay fixed.
def _clear_box(cy, e=e):
return shapely.box(e.x0 - clearance, cy - half - clearance,
e.x1 + clearance, cy + half + clearance)
if obstacles is not None and _clear_box(base_cy).intersects(obstacles):
# Scan the preferred side first; if it is blocked to the axes edge (e.g. a
# dot pinned against the floor) try the opposite side.
sign = 1 if side == "above" else -1
for direction in (sign, -sign):
cy = base_cy + direction * step
while lo_c <= cy <= hi_c and _clear_box(cy).intersects(obstacles):
cy += direction * step
if lo_c <= cy <= hi_c:
t.set_position((t.get_position()[0], inv.transform((axbb.x0, cy))[1]))
break
e = t.get_window_extent(renderer)
box = shapely.box(e.x0, e.y0, e.x1, e.y1)
obstacles = box if obstacles is None else shapely.union_all([obstacles, box])
return texts
def plot_excess_free_energy(
df,
col_wrap=3,
height=3.0,
aspect=1.3,
color_override=None,
convex_hull=True,
inline_legend=True,
):
"""Plot excess free energy vs concentration for competing phases.
Takes a pre-computed DataFrame from ``calc_phase_diagram(..., keep_unstable=True)``
and delegates to seaborn ``relplot``.
When ``convex_hull=True``, stable solution phases render as solid curves,
metastable/unstable regions as faded lines (same colour, alpha=0.4), and
the common-tangent construction as black dotted segments — one segment per
coexistence region (grouped by ``mu``) — with black vertex markers.
Line phases always render as a single coloured scatter dot.
When ``convex_hull=False``, all solution phases render as plain solid curves
regardless of stability.
Args:
df: DataFrame from ``calc_phase_diagram(..., keep_unstable=True)`` with
columns ``c``, ``f_excess``, ``phase``, ``T``, ``stable``, and
optionally ``border`` and ``mu``.
col_wrap: Maximum subplot columns per row.
height: Height of each facet in inches (multiple temperatures only).
aspect: Width-to-height ratio of each facet (multiple temperatures only).
color_override: Optional ``dict[name -> color]`` overriding phase colours.
convex_hull: If True, distinguish stable (solid curves) from metastable
(faded lines) and overlay the common-tangent segments in black.
If False, all solution phases render as plain solid curves.
inline_legend: If True (default), drop the figure legend and label each
phase just off its line -- solution curves above the centre of their
largest continuous region, line phases below their dot -- white-outlined
and colour-matched, with overlapping labels spread apart vertically. If
False, keep the figure legend box on the right.
Returns:
For a single temperature, the current :class:`matplotlib.axes.Axes`
holding the plain lineplot; no figure is allocated, so pre-create one
to control its size. For multiple temperatures, a
:class:`seaborn.FacetGrid` with one column per temperature (figure via
``.fig``, axes via ``.axes``).
"""
import matplotlib.lines as mlines
df = df.copy()
if df.empty:
raise ValueError("df is empty.")
if "border" not in df.columns:
df["border"] = False
temperatures = sorted(df["T"].unique())
col_wrap = min(col_wrap, len(temperatures))
# Line phases have a fixed concentration — detect by zero range across all rows.
phase_names = list(df["phase"].unique())
c_range = df.groupby("phase")["c"].apply(lambda s: s.max() - s.min())
line_phase_names = set(c_range[c_range < 1e-9].index)
muted_colors = sns.color_palette("muted")
palette = {name: muted_colors[i % len(muted_colors)] for i, name in enumerate(phase_names)}
if color_override:
palette.update({k: v for k, v in color_override.items() if k in palette})
df_sol = df[~df["phase"].isin(line_phase_names)].copy()
df_lp = df[df["phase"].isin(line_phase_names)].copy()
sol_palette = {k: v for k, v in palette.items() if k not in line_phase_names}
sol_hue_order = [n for n in phase_names if n not in line_phase_names]
if convex_hull:
base_data = df_sol[df_sol["stable"]].copy()
if not base_data.empty:
base_data = cluster_phase(base_data, distance_threshold=0.1)
units_col = "phase_id"
else:
units_col = None
else:
base_data = df_sol
units_col = None
# A single temperature does not need a facet grid; draw a plain lineplot onto
# the current axes (the caller controls figure allocation) and return that axes.
# Multiple temperatures fan out into a FacetGrid.
single = len(temperatures) == 1
if single:
single_ax = plt.gca()
if not base_data.empty:
sns.lineplot(
data=base_data,
x="c",
y="f_excess",
hue="phase",
hue_order=sol_hue_order,
units=units_col,
palette=sol_palette,
estimator=None,
errorbar=None,
linewidth=2.5,
ax=single_ax,
)
g = None
axes = [single_ax]
elif base_data.empty:
# No solution-phase rows (e.g. all phases are line phases, or all solution rows are
# unstable). sns.relplot cannot create facets from an empty DataFrame, so build
# the grid structure directly and skip the line-drawing step.
g = sns.FacetGrid(
data=df[["T"]].drop_duplicates(),
col="T",
col_order=temperatures,
col_wrap=col_wrap,
height=height,
aspect=aspect,
)
axes = list(g.axes.flat)
else:
g = sns.relplot(
data=base_data,
x="c",
y="f_excess",
hue="phase",
hue_order=sol_hue_order,
units=units_col,
col="T",
col_wrap=col_wrap,
palette=sol_palette,
height=height,
aspect=aspect,
kind="line",
estimator=None,
errorbar=None,
linewidth=2.5,
)
axes = list(g.axes.flat)
facet_entries = [] # (ax, [(label, x, y, color), ...]) for inline labelling
for ax, T_val in zip(axes, temperatures):
sub_all = df[df["T"] == T_val]
sub_sol = df_sol[df_sol["T"] == T_val]
sub_lp = df_lp[df_lp["T"] == T_val]
if convex_hull:
# Metastable solution phases: faded lines, same colour as stable.
# T is constant within this facet, so cluster_T_c reduces to c-only
# clustering, correctly splitting disjoint metastable c-ranges.
unstable = sub_sol[~sub_sol["stable"]]
for pname, grp in unstable.groupby("phase"):
seg_ids = cluster_T_c(grp, distance_threshold=0.1)
color = sol_palette.get(pname, "gray")
for seg_id in seg_ids.unique():
seg = grp.loc[seg_ids == seg_id].sort_values("c")
ax.plot(
seg["c"].values, seg["f_excess"].values,
color=color, alpha=0.4, lw=2.5, zorder=2,
)
# Line phases: single colored dot at fixed concentration.
for _, row in sub_lp.drop_duplicates("phase").iterrows():
ax.scatter(
[row["c"]], [row["f_excess"]],
color=palette.get(row["phase"], "k"),
zorder=5, s=80,
)
if convex_hull:
# Common-tangent lines: one dotted segment per coexistence region.
# Rows sharing the same mu value belong to one two-phase equilibrium,
# so grouping by mu gives one tangent line per coexistence pair.
bd_all = sub_all[sub_all["border"]]
for _mu, grp in bd_all.groupby("mu"):
grp_sorted = grp.drop_duplicates(subset=["c", "f_excess"]).sort_values("c")
if len(grp_sorted) >= 2:
ax.plot(
grp_sorted["c"].values, grp_sorted["f_excess"].values,
ls="dotted", color="k", zorder=3, lw=1.5,
)
# Hull vertex markers: deduplicated across all mu values.
bd_unique = bd_all.drop_duplicates(subset=["c", "f_excess"]).sort_values("c")
if not bd_unique.empty:
ax.scatter(
bd_unique["c"].values, bd_unique["f_excess"].values,
color="k", s=25, zorder=7,
)
if inline_legend:
# Label each solution phase above its largest continuous region: a
# free-energy curve is convex, so the space above it is open whether the
# phase is stable (lower-hull) or metastable (upper arc). Anchoring on
# the largest continuous c-region keeps the label on the prominent stretch
# rather than a sliver (a phase stable only in the dilute corners would
# otherwise pull its label into a corner). Within that arc the x-anchor
# is fanned out by the phase's rank instead of always sitting at the arc
# centre, so phases whose arcs all span the range don't pile their labels
# in the middle. Line phases sit on the lower hull, so their labels go
# below the dot. T is constant within a facet, so cluster_T_c splits
# regions by c alone. Placement is deferred until refline/limits settle.
arcs = []
for pname in sol_hue_order:
pg = sub_sol[sub_sol["phase"] == pname].dropna(subset=["f_excess"])
if pg.empty:
continue
region = pg
if len(region) > 1:
seg = cluster_T_c(region, distance_threshold=0.1)
extent = region.groupby(seg)["c"].agg(lambda c: c.max() - c.min())
region = region.loc[seg == extent.idxmax()]
region = region.sort_values("c")
arcs.append((
pname, region["c"].to_numpy(), region["f_excess"].to_numpy(),
sol_palette.get(pname, "k"),
))
# Spread the anchors left-to-right: order the arcs by their centre, then
# place each phase's label at the rank-derived fraction of its own arc.
# The fraction is mapped into the arc's inner 60% so a label never lands
# on the curve's edge where it turns up.
arcs.sort(key=lambda a: (0.5 * (a[1][0] + a[1][-1]), a[0]))
n = len(arcs)
margin = 0.2 # keep anchors clear of each arc's turning-up ends
entries = []
for i, (pname, cs, fs, color) in enumerate(arcs):
frac = margin + (1 - 2 * margin) * (i + 0.5) / n
x = cs[0] + frac * (cs[-1] - cs[0])
y = np.interp(x, cs, fs)
entries.append((pname, x, y, color, "above"))
for _, row in sub_lp.drop_duplicates("phase").iterrows():
entries.append((row["phase"], row["c"], row["f_excess"], palette.get(row["phase"], "k"), "below"))
facet_entries.append((ax, entries))
# The single-temperature path draws onto one axes whose legend lives on the axes
# itself; the FacetGrid path keeps its legend on ``g._legend``. Resolve both here.
legend = g._legend if g is not None else axes[0].get_legend()
figure = g.figure if g is not None else axes[0].figure
if inline_legend:
# Inline labels replace the legend box; drop seaborn's legend.
if legend is not None:
legend.remove()
# Add line phases to the figure legend.
elif not df_lp.empty:
lp_handles = [
mlines.Line2D(
[0], [0], marker="o", color="w",
markerfacecolor=palette.get(n, "k"), markersize=8, label=n,
)
for n in sorted(line_phase_names)
]
if legend is not None:
existing_handles = list(legend.legend_handles)
existing_labels = [t.get_text() for t in legend.texts]
legend.remove()
else:
existing_handles = []
existing_labels = []
figure.legend(
existing_handles + lp_handles,
existing_labels + [h.get_label() for h in lp_handles],
title="phase",
bbox_to_anchor=(1.02, 0.5),
loc="center left",
borderaxespad=0,
)
if g is not None:
g.refline(y=0)
g.set_titles("T = {col_name:.0f} K")
g.set(xlabel="Concentration", ylabel="Free Energy of Formation")
else:
ax = axes[0]
ax.axhline(0, color=".5", linestyle="--")
ax.set(xlabel="Concentration", ylabel="Free Energy of Formation")
# Place inline labels last so they see the final axis limits (refline/relplot).
for ax, entries in facet_entries:
_add_inline_curve_labels(ax, entries)
return g if g is not None else axes[0]