diff --git a/doc/_basic_usage.md b/doc/_basic_usage.md index 361d6bb..02a8fea 100644 --- a/doc/_basic_usage.md +++ b/doc/_basic_usage.md @@ -27,6 +27,59 @@ for i, row in circuits.iterrows(): ae.visualization.print_inference_results(circuits) ``` +## Posterior quality diagnostics + +After Bayesian inference, AutoEIS can score each posterior parameter +distribution for Gaussian-like behavior. The diagnostic combines Q-Q plot +linearity, skewness, and KDE mode count, then labels each parameter as +`Normal`, `Multi-Peak`, `Left-Skewed`, `Right-Skewed`, `Touch-Bound`, or +`Invalid`. + +```python +result = results[0] + +# Return one row per inferred parameter. +quality = result.assess_posterior_quality(threshold=79) +print(quality[["parameter", "score", "shape", "status"]]) + +# The same API is also available as a module-level function. +quality = ae.diagnostics.assess_posterior_quality(result.mcmc, threshold=79) + +# Optional: save KDE + histogram plots for each parameter. +quality = result.assess_posterior_quality( + threshold=79, + save_plots=True, + save_dir="posterior_quality", +) + +# Add posterior scores back to a circuits dataframe and rank models. +circuits = ae.diagnostics.add_posterior_quality(circuits, results, threshold=79) +``` + +## ECM topology vectors + +ECM topology utilities are available from `ae.ecm`. They operate on the same +AutoEIS circuit strings, but produce normalized topology trees and compact +vectors for tracking or comparing model structure. + +```python +circuit = "R1-[P2,R3]-C4" + +# Simplify mergeable R/C/L branches without changing the existing parser API. +simplified = ae.ecm.simplify_ecm("R1-R2-[C1-C2,R3]-R4") + +# Convert ECM topology to an integer vector and back. +ecm_vector = ae.ecm.ecm_to_vector(circuit) +roundtrip = ae.ecm.vector_to_ecm(ecm_vector) + +# Convert numeric EIS values with matching series/parallel syntax. +eis_vector = ae.ecm.eis_to_vector("10-[2.5,30]-1e-6") +values = ae.ecm.extract_eis_values(eis_vector) + +# Add vectors to a circuits dataframe. +circuits = ae.ecm.add_ecm_vectors(circuits) +``` + :::{seealso} While the above example should work out of the box, it is recommended to use AutoEIS in a step by step fashion to have more control over the analysis process. Furthermore, you'll learn more about the inner workings of AutoEIS this way. An example notebook that demonstrates how to use AutoEIS in more details can be found [here](https://github.com/AUTODIAL/AutoEIS/blob/develop/examples/autoeis_demo.ipynb). ::: diff --git a/doc/modules.rst b/doc/modules.rst index 8e0dd4c..0e56085 100644 --- a/doc/modules.rst +++ b/doc/modules.rst @@ -3,7 +3,9 @@ API Reference ############# -AutoEIS contains five modules: ``io`` for reading and writing data, ``core`` for finding the best fit and performing the analysis, ``visualization`` for plotting, ``utils`` for helper functions used in ``core``, and ``parser`` for parsing circuit strings. +AutoEIS contains modules for reading and writing data, finding the best fit +and performing the analysis, posterior diagnostics, ECM topology tools, +plotting, helper utilities, and parsing circuit strings. .. automodule:: autoeis @@ -13,6 +15,8 @@ AutoEIS contains five modules: ``io`` for reading and writing data, ``core`` for autoeis.io autoeis.core + autoeis.diagnostics + autoeis.ecm autoeis.visualization autoeis.utils autoeis.parser diff --git a/pyproject.toml b/pyproject.toml index 543dd38..578ee5f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,6 +46,7 @@ dependencies = [ "pyparsing>=3", "rich", "scikit-learn>=1.4", + "scipy", "seaborn", "tqdm", "numexpr", diff --git a/src/autoeis/__init__.py b/src/autoeis/__init__.py index 8bf6e82..00c393d 100644 --- a/src/autoeis/__init__.py +++ b/src/autoeis/__init__.py @@ -26,7 +26,16 @@ def _setup_logger(): config = _Settings() _setup_logger() -from . import core, io, metrics, parser, utils, visualization # noqa: F401, E402 +from . import ( # noqa: F401, E402 + core, + diagnostics, + ecm, + io, + metrics, + parser, + utils, + visualization, +) from .core import * # noqa: E402 from .version import __equivalent_circuits_jl_version__, __version__ # noqa: F401, E402 from .visualization import rich_print # noqa: F401, E402 diff --git a/src/autoeis/diagnostics.py b/src/autoeis/diagnostics.py new file mode 100644 index 0000000..8462494 --- /dev/null +++ b/src/autoeis/diagnostics.py @@ -0,0 +1,552 @@ +""" +Posterior distribution diagnostics for Bayesian inference results. + +.. currentmodule:: autoeis.diagnostics + +.. autosummary:: + :toctree: generated/ + + posterior_gaussianity_score + assess_posterior_sample + assess_posterior_quality + plot_posterior_quality + +""" + +from __future__ import annotations + +import logging +from collections.abc import Mapping, Sequence +from dataclasses import dataclass +from pathlib import Path +from typing import Any + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from scipy.stats import gaussian_kde, probplot + +log = logging.getLogger(__name__) + +__all__ = [ + "PosteriorQuality", + "posterior_gaussianity_score", + "assess_posterior_sample", + "assess_posterior_quality", + "assess_inference_results", + "add_posterior_quality", + "plot_posterior_quality", +] + +_AUXILIARY_SAMPLE_NAMES = {"lp__", "log_likelihood"} + + +@dataclass(frozen=True) +class PosteriorQuality: + """Quality assessment for one posterior parameter distribution.""" + + parameter: str + score: float + r2: float + skewness: float + modes: int + shape: str + status: str + threshold: float + peak_locations: tuple[float, ...] = () + peak_densities: tuple[float, ...] = () + + @property + def good(self) -> bool: + """Whether the posterior distribution passes the score threshold.""" + return self.status == "Good" + + def to_dict(self) -> dict[str, Any]: + """Return a serializable dictionary representation.""" + return { + "parameter": self.parameter, + "score": self.score, + "r2": self.r2, + "skewness": self.skewness, + "modes": self.modes, + "shape": self.shape, + "status": self.status, + "good": self.good, + "threshold": self.threshold, + "peak_locations": self.peak_locations, + "peak_densities": self.peak_densities, + } + + +def _as_sample_array(sample: Sequence[float] | np.ndarray) -> np.ndarray: + """Convert a posterior sample to a clean flat finite array.""" + sample = np.asarray(sample, dtype=float).ravel() + sample = sample[np.isfinite(sample)] + if sample.size < 3: + raise ValueError("Posterior sample must contain at least three finite values.") + return sample + + +def _skewness(sample: np.ndarray) -> float: + """Calculate skewness using central moments.""" + centered = sample - np.mean(sample) + m2 = np.mean(centered**2) + if m2 <= 0: + return 0.0 + m3 = np.mean(centered**3) + return float(m3 / (m2**1.5)) + + +def _qq_r2(sample: np.ndarray) -> float: + """Return R^2 of the normal Q-Q plot linear fit.""" + (theoretical, ordered), _ = probplot(sample, dist="norm") + if np.allclose(ordered, ordered[0]): + return 1.0 + + slope, intercept = np.polyfit(theoretical, ordered, deg=1) + fitted = slope * theoretical + intercept + ss_res = np.sum((ordered - fitted) ** 2) + ss_tot = np.sum((ordered - np.mean(ordered)) ** 2) + if ss_tot <= 0: + return 1.0 + return float(np.clip(1 - ss_res / ss_tot, 0, 1)) + + +def _detect_peaks( + x: np.ndarray, + y: np.ndarray, + *, + min_height_fraction: float = 0.4, + threshold_scale: float = 6e-5, + min_distance_fraction: float = 0.4, +) -> list[int]: + """Detect KDE peaks by consensus thresholding and non-maximum suppression.""" + if x.size < 3 or y.size < 3 or not np.any(y): + return [] + + thresholds = [ + threshold_scale * (float(np.max(y)) - float(np.min(y))), + threshold_scale * float(np.std(y)), + threshold_scale * float(np.subtract(*np.percentile(y, [75, 25]))), + ] + + candidates_by_method = [] + for threshold in thresholds: + peaks = [] + for idx in range(1, len(y) - 1): + is_peak = y[idx] > y[idx - 1] and y[idx] > y[idx + 1] + is_prominent = y[idx] - min(y[idx - 1], y[idx + 1]) > threshold + if is_peak and is_prominent: + peaks.append(idx) + candidates_by_method.append(set(peaks)) + + consensus = [] + for idx in sorted(set.union(*candidates_by_method) if candidates_by_method else set()): + votes = sum(idx in candidates for candidates in candidates_by_method) + if votes >= 2: + consensus.append(idx) + + height_threshold = min_height_fraction * float(np.max(y)) + candidates = [idx for idx in consensus if y[idx] >= height_threshold] + candidates.sort(key=lambda idx: y[idx], reverse=True) + + dx = float(np.median(np.diff(np.sort(x)))) + if dx <= 0: + return sorted(candidates) + + bandwidth = 1.06 * float(np.std(x)) * len(x) ** (-1 / 5) + min_distance = max(1, int(round((min_distance_fraction * bandwidth) / dx))) + + kept = [] + for idx in candidates: + if all(abs(idx - other) >= min_distance for other in kept): + kept.append(idx) + + return sorted(kept) + + +def _estimate_modes( + sample: np.ndarray, n_grid: int +) -> tuple[int, np.ndarray, np.ndarray, list[int]]: + """Estimate posterior modes from a KDE curve.""" + if np.allclose(sample, sample[0]): + x = np.array([sample[0], sample[0]]) + y = np.array([1.0, 1.0]) + return 1, x, y, [0] + + x = np.linspace(float(np.min(sample)), float(np.max(sample)), n_grid) + try: + y = gaussian_kde(sample)(x) + except np.linalg.LinAlgError: + log.debug("KDE failed for posterior sample; treating distribution as invalid.") + return 0, x, np.zeros_like(x), [] + + peaks = _detect_peaks(x, y) + return len(peaks), x, y, peaks + + +def posterior_gaussianity_score(r2: float, skewness: float, modes: int) -> float: + """Compute the PART I Gaussianity score from Q-Q fit, skewness, and modes.""" + if not np.isfinite(r2) or not np.isfinite(skewness) or modes < 1: + return 0.0 + score = r2 * 100 + score -= abs(skewness) * 5 + score -= modes * 9 + return float(max(0, round(score, 2))) + + +def _classify_distribution(skewness: float, modes: int, skew_threshold: float) -> str: + """Classify a posterior shape from mode count and skewness.""" + if not np.isfinite(skewness) or modes < 1: + return "Invalid" + if modes > 1: + return "Multi-Peak" + if abs(skewness) >= 2: + return "Touch-Bound" + if skewness <= -skew_threshold: + return "Left-Skewed" + if skewness >= skew_threshold: + return "Right-Skewed" + return "Normal" + + +def assess_posterior_sample( + sample: Sequence[float] | np.ndarray, + parameter: str = "parameter", + *, + threshold: float = 79, + skew_threshold: float = 0.5, + n_grid: int = 300, +) -> PosteriorQuality: + """Assess one posterior parameter distribution. + + Parameters + ---------- + sample + Posterior samples for one parameter. + parameter + Parameter name used in returned diagnostics and plots. + threshold + Minimum Gaussianity score required for ``status == "Good"``. + skew_threshold + Absolute skewness threshold used to label skewed distributions. + n_grid + Number of KDE grid points used for mode detection. + + Returns + ------- + PosteriorQuality + Diagnostic metrics and classification for the sample. + """ + sample = _as_sample_array(sample) + + r2 = _qq_r2(sample) + skew = _skewness(sample) + modes, x, y, peaks = _estimate_modes(sample, n_grid=n_grid) + score = posterior_gaussianity_score(r2, skew, modes) + shape = _classify_distribution(skew, modes, skew_threshold) + status = "Good" if score >= threshold else "Bad" + + return PosteriorQuality( + parameter=parameter, + score=score, + r2=round(r2, 6), + skewness=round(skew, 6), + modes=modes, + shape=shape, + status=status, + threshold=threshold, + peak_locations=tuple(float(x[idx]) for idx in peaks), + peak_densities=tuple(float(y[idx]) for idx in peaks), + ) + + +def _extract_samples(obj: Any) -> Mapping[str, np.ndarray]: + """Extract sample mapping from MCMC, InferenceResult, or mapping-like input.""" + if isinstance(obj, Mapping): + return obj + if hasattr(obj, "samples"): + return obj.samples + if hasattr(obj, "get_samples"): + return obj.get_samples() + if hasattr(obj, "mcmc") and hasattr(obj.mcmc, "get_samples"): + return obj.mcmc.get_samples() + raise TypeError( + "Expected an InferenceResult, a numpyro MCMC object, or a mapping of samples." + ) + + +def assess_posterior_quality( + samples: Any, + *, + threshold: float = 79, + var_names: Sequence[str] | None = None, + skip_auxiliary: bool = True, + save_plots: bool = False, + save_dir: str | Path | None = None, + n_grid: int = 300, +) -> pd.DataFrame: + """Assess all posterior parameter distributions in an inference result. + + Parameters + ---------- + samples + An AutoEIS ``InferenceResult``, NumPyro ``MCMC`` object, or mapping from + parameter names to posterior sample arrays. + threshold + Minimum Gaussianity score required for a parameter to pass. + var_names + Optional subset of parameter names to assess. + skip_auxiliary + Skip sampler bookkeeping variables such as ``lp__`` and ``log_likelihood``. + save_plots + If True, save one KDE/histogram diagnostic plot per parameter. + save_dir + Directory for saved plots. Defaults to ``posterior_quality`` when + ``save_plots=True``. + n_grid + Number of KDE grid points used for mode detection. + + Returns + ------- + pandas.DataFrame + One row per parameter with score, shape, and pass/fail status. + """ + sample_mapping = _extract_samples(samples) + requested = set(var_names) if var_names is not None else None + plot_dir = Path(save_dir or "posterior_quality") if save_plots else None + + results = [] + for name, sample in sample_mapping.items(): + if skip_auxiliary and name in _AUXILIARY_SAMPLE_NAMES: + continue + if requested is not None and name not in requested: + continue + + result = assess_posterior_sample( + sample, + str(name), + threshold=threshold, + n_grid=n_grid, + ) + results.append(result.to_dict()) + + if plot_dir is not None: + fig, _ = plot_posterior_quality(sample, result=result, save_dir=plot_dir) + plt.close(fig) + + columns = [ + "parameter", + "score", + "r2", + "skewness", + "modes", + "shape", + "status", + "good", + "threshold", + "peak_locations", + "peak_densities", + ] + return pd.DataFrame(results, columns=columns) + + +def assess_inference_results( + results: Sequence[Any], + *, + threshold: float = 79, + save_plots: bool = False, + save_dir: str | Path | None = None, + n_grid: int = 300, +) -> pd.DataFrame: + """Summarize posterior quality for a sequence of inference results.""" + rows = [] + for idx, result in enumerate(results): + circuit = getattr(result, "circuit", None) + converged = bool(getattr(result, "converged", True)) + if converged: + result_save_dir = None + if save_plots: + result_save_dir = Path(save_dir or "posterior_quality") / f"inference_{idx}" + quality = assess_posterior_quality( + result, + threshold=threshold, + save_plots=save_plots, + save_dir=result_save_dir, + n_grid=n_grid, + ) + bad = quality.loc[~quality["good"], "parameter"].to_list() + posterior_score = float(quality["score"].sum()) + else: + quality = pd.DataFrame() + bad = ["Not Converged"] + posterior_score = 0.0 + + rows.append( + { + "circuit": circuit, + "InferenceResult": result, + "Posterior Score": posterior_score, + "Bad Parameters": ", ".join(bad) if bad else "None", + "n_bad_parameters": len(bad), + "posterior_quality": quality, + } + ) + + return pd.DataFrame(rows) + + +def add_posterior_quality( + circuits: pd.DataFrame, + results: Sequence[Any] | None = None, + *, + freq: np.ndarray | None = None, + Z: np.ndarray | None = None, + threshold: float = 79, + save_plots: bool = False, + save_dir: str | Path | None = None, + compute_metrics: bool = False, + sort: bool = True, + n_grid: int = 300, +) -> pd.DataFrame: + """Add posterior quality scores to a circuits dataframe. + + Parameters + ---------- + circuits + Circuits dataframe, usually the dataframe passed to Bayesian inference. + results + Optional sequence of AutoEIS ``InferenceResult`` objects. If omitted, + ``circuits`` must already contain an ``InferenceResult`` column. + freq, Z + Frequency and impedance arrays. Required only when ``compute_metrics=True``. + threshold + Minimum per-parameter Gaussianity score required to pass. + save_plots + If True, save one KDE/histogram diagnostic plot per parameter. + save_dir + Directory for saved plots. + compute_metrics + If True, also call :func:`autoeis.core.compute_fitness_metrics`. + sort + Sort circuits without bad parameters first, then by posterior score. + n_grid + Number of KDE grid points used for mode detection. + + Returns + ------- + pandas.DataFrame + Copy of ``circuits`` with posterior quality columns added. + """ + circuits = circuits.copy() + if results is None: + if "InferenceResult" not in circuits: + raise ValueError("Provide results or a circuits dataframe with InferenceResult.") + results = circuits["InferenceResult"].to_list() + elif len(results) != len(circuits): + raise ValueError("results must have the same length as circuits.") + else: + circuits["InferenceResult"] = list(results) + + quality = assess_inference_results( + results, + threshold=threshold, + save_plots=save_plots, + save_dir=save_dir, + n_grid=n_grid, + ) + for column in ["Posterior Score", "Bad Parameters", "n_bad_parameters"]: + circuits[column] = quality[column].to_list() + circuits["Posterior Quality"] = quality["posterior_quality"].to_list() + + if compute_metrics: + if freq is None or Z is None: + raise ValueError("freq and Z are required when compute_metrics=True.") + from .core import compute_fitness_metrics + + circuits = compute_fitness_metrics(circuits, freq, Z) + + if sort: + has_bad = circuits["Bad Parameters"].str.strip().str.lower() != "none" + circuits = circuits.assign(_has_bad_parameters=has_bad) + circuits.sort_values( + by=["_has_bad_parameters", "Posterior Score"], + ascending=[True, False], + inplace=True, + ignore_index=True, + ) + circuits.drop(columns=["_has_bad_parameters"], inplace=True) + + circuits["Posterior Score Rank"] = np.arange(1, len(circuits) + 1) + return circuits + + +def plot_posterior_quality( + sample: Sequence[float] | np.ndarray, + *, + result: PosteriorQuality | None = None, + parameter: str | None = None, + threshold: float = 79, + n_grid: int = 300, + save_dir: str | Path | None = None, + ax: Sequence[plt.Axes] | None = None, +) -> tuple[plt.Figure, np.ndarray]: + """Plot posterior KDE, detected peaks, and histogram diagnostics.""" + sample = _as_sample_array(sample) + if result is None: + result = assess_posterior_sample( + sample, + parameter=parameter or "parameter", + threshold=threshold, + n_grid=n_grid, + ) + + modes, x, y, peaks = _estimate_modes(sample, n_grid=n_grid) + if ax is None: + fig, axes = plt.subplots(2, 1, figsize=(8, 6), layout="constrained") + else: + axes = np.asarray(ax) + fig = axes[0].figure + + axes[0].plot(x, y, label="KDE Distribution") + if peaks: + axes[0].scatter( + x[peaks], + y[peaks], + color="red", + marker="x", + s=90, + label=f"Peaks: {modes}", + ) + + summary = ( + f"Score: {result.score:.2f}\n" + f"R2: {result.r2:.4f}\n" + f"Skewness: {result.skewness:.4f}" + ) + axes[0].text( + 0.98, + 0.95, + summary, + transform=axes[0].transAxes, + va="top", + ha="right", + bbox={"facecolor": "white", "alpha": 0.75, "edgecolor": "0.5"}, + ) + axes[0].set_title(f"{result.parameter}: {result.shape} | {result.status}") + axes[0].set_xlabel("Parameter value") + axes[0].set_ylabel("Density") + axes[0].legend() + + axes[1].hist(sample, bins=30, alpha=0.75, color="tab:blue", edgecolor="black") + axes[1].set_title("Posterior sample histogram") + axes[1].set_xlabel("Parameter value") + axes[1].set_ylabel("Count") + + if save_dir is not None: + save_dir = Path(save_dir) + save_dir.mkdir(parents=True, exist_ok=True) + safe_name = "".join(c if c.isalnum() or c in "._-" else "_" for c in result.parameter) + fig.savefig(save_dir / f"{safe_name}-{result.score:.2f}.png", dpi=150) + + return fig, axes diff --git a/src/autoeis/ecm.py b/src/autoeis/ecm.py new file mode 100644 index 0000000..ad9f539 --- /dev/null +++ b/src/autoeis/ecm.py @@ -0,0 +1,914 @@ +""" +Utilities for ECM topology simplification and vectorization. + +.. currentmodule:: autoeis.ecm + +.. autosummary:: + :toctree: generated/ + + parse_circuit_tree + normalize_tree + simplify_ecm + compare_ecm_topology + deduplicate_ecms + ecm_to_vector + vector_to_ecm + parse_eis_value_tree + eis_to_vector + ecm_to_eis_string + +""" + +from __future__ import annotations + +import json +import math +import re +from collections import OrderedDict, defaultdict +from collections.abc import Iterable, Mapping, Sequence +from pathlib import Path +from typing import Any + +import numpy as np +import pandas as pd + +COMPONENT_RE = re.compile(r"[RCLPW]\d+", re.IGNORECASE) +VALUE_RE = re.compile(r"[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?") + +ECM_START_CODES = {"series": 0, "parallel": 1} +ECM_LEAF_CODES = {"R": 2, "C": 3, "L": 4, "P": 5, "W": 6} +ECM_END_CODES = {"series": -2, "parallel": -3} + +ECM_CODE_TO_NODE = { + **{v: k for k, v in ECM_START_CODES.items()}, + **{v: k for k, v in ECM_LEAF_CODES.items()}, +} +ECM_END_CODE_TO_NODE = {v: k for k, v in ECM_END_CODES.items()} + +EIS_START_CODES = {"series": 0, "parallel": 1} +EIS_VALUE_CODE = 2 +EIS_END_CODES = {"series": -2, "parallel": -3} + +CircuitTree = str | tuple[str, list["CircuitTree"]] +ValueTree = float | tuple[str, list["ValueTree"]] + +__all__ = [ + "CircuitTree", + "ValueTree", + "split_top_level", + "parse_circuit_tree", + "tree_to_circuit", + "normalize_tree", + "simplify_tree", + "simplify_ecm", + "compare_ecm_topology", + "deduplicate_ecms", + "topology_to_weighted_vector", + "cosine_similarity", + "structure_similarity", + "analyze_circuit_structure", + "compare_circuit_structures", + "ecm_to_vector", + "vector_to_ecm", + "parse_eis_value_tree", + "eis_to_vector", + "extract_eis_values", + "replace_eis_values", + "simplify_ecm_eis_pair", + "ecm_to_eis_string", + "geometric_mean_frequency", + "compress_parameter_values", + "order_parameters_by_topology", + "select_best_circuit_row", + "extract_vectorizable_params", + "add_ecm_vectors", + "ECMVectorTracker", + "EISVectorTracker", +] + + +def split_top_level(text: str, sep: str) -> list[str]: + """Split a circuit string on a separator only at bracket depth zero.""" + parts: list[str] = [] + buf: list[str] = [] + depth = 0 + for char in text: + if char == "[": + depth += 1 + elif char == "]": + depth -= 1 + if depth < 0: + raise ValueError(f"Unbalanced bracket in {text!r}.") + + if char == sep and depth == 0: + parts.append("".join(buf)) + buf = [] + else: + buf.append(char) + + if depth != 0: + raise ValueError(f"Unbalanced bracket in {text!r}.") + parts.append("".join(buf)) + return parts + + +def _component_token(token: str, preserve_labels: bool) -> str: + token = token.strip() + match = COMPONENT_RE.fullmatch(token) + if match: + value = match.group(0) + return value.upper() if preserve_labels else value[0].upper() + + type_match = re.fullmatch(r"[RCLPW]", token, flags=re.IGNORECASE) + if type_match: + return token.upper() + + raise ValueError(f"Invalid ECM component token: {token!r}.") + + +def parse_circuit_tree(circuit: str, *, preserve_labels: bool = True) -> CircuitTree: + """Parse an AutoEIS circuit string into a nested series/parallel tree.""" + text = circuit.strip().replace(" ", "") + if not text: + raise ValueError("Circuit string is empty.") + + series_parts = split_top_level(text, "-") + if len(series_parts) > 1: + return ( + "series", + [ + parse_circuit_tree(part, preserve_labels=preserve_labels) + for part in series_parts + ], + ) + + if text.startswith("[") and text.endswith("]"): + inner = text[1:-1] + parallel_parts = split_top_level(inner, ",") + if len(parallel_parts) < 2: + raise ValueError(f"Parallel group must contain at least two children: {text!r}.") + return ( + "parallel", + [ + parse_circuit_tree(part, preserve_labels=preserve_labels) + for part in parallel_parts + ], + ) + + return _component_token(text, preserve_labels) + + +def _component_type(component: str) -> str: + return component[0].upper() + + +def _tree_sort_key(node: CircuitTree | ValueTree) -> str: + return tree_to_circuit(node) if isinstance(node, tuple) else str(node) + + +def normalize_tree( + tree: CircuitTree, + *, + sort: bool = True, + flatten: bool = True, +) -> CircuitTree: + """Normalize a circuit tree by flattening equal groups and optionally sorting.""" + if isinstance(tree, str): + return tree + + label, children = tree + normalized_children = [ + normalize_tree(child, sort=sort, flatten=flatten) for child in children + ] + + flat: list[CircuitTree] = [] + for child in normalized_children: + if flatten and isinstance(child, tuple) and child[0] == label: + flat.extend(child[1]) + else: + flat.append(child) + + if len(flat) == 1: + return flat[0] + if sort: + flat = sorted(flat, key=_tree_sort_key) + return (label, flat) + + +def simplify_tree(tree: CircuitTree, *, preserve_order: bool = True) -> CircuitTree | None: + """Simplify an ECM tree by merging same-type R/C/L leaves within each group.""" + if isinstance(tree, str): + return tree + + label, children = tree + simplified = [simplify_tree(child, preserve_order=preserve_order) for child in children] + simplified = [child for child in simplified if child is not None] + if not simplified: + return None + if len(simplified) == 1: + return simplified[0] + + flattened: list[CircuitTree] = [] + for child in simplified: + if isinstance(child, tuple) and child[0] == label: + flattened.extend(child[1]) + else: + flattened.append(child) + + representatives: dict[str, str] = {} + output: list[CircuitTree] = [] + for child in flattened: + if isinstance(child, str) and _component_type(child) in {"R", "C", "L"}: + ctype = _component_type(child) + if ctype not in representatives: + representatives[ctype] = child + output.append(child) + else: + output.append(child) + + if not preserve_order: + output = sorted(output, key=_tree_sort_key) + if len(output) == 1: + return output[0] + return (label, output) + + +def tree_to_circuit(tree: CircuitTree | ValueTree) -> str: + """Convert a circuit or numeric value tree back to a string.""" + if isinstance(tree, tuple): + label, children = tree + sep = "-" if label == "series" else "," + body = sep.join(tree_to_circuit(child) for child in children) + return body if label == "series" else f"[{body}]" + if isinstance(tree, float): + return f"{tree:.12g}" + return str(tree) + + +def simplify_ecm(circuit: str, *, preserve_labels: bool = True) -> str: + """Simplify an AutoEIS ECM string while preserving the first label per merged type.""" + tree = parse_circuit_tree(circuit, preserve_labels=preserve_labels) + simplified = simplify_tree(tree) + if simplified is None: + return "" + return tree_to_circuit(simplified) + + +def _topology_tree(circuit: str) -> CircuitTree: + return normalize_tree(parse_circuit_tree(circuit, preserve_labels=False)) + + +def compare_ecm_topology(circuit1: str, circuit2: str, *, simplify: bool = True) -> bool: + """Return True when two ECM strings have equivalent normalized topology.""" + tree1 = parse_circuit_tree(circuit1, preserve_labels=False) + tree2 = parse_circuit_tree(circuit2, preserve_labels=False) + if simplify: + tree1 = simplify_tree(tree1) or tree1 + tree2 = simplify_tree(tree2) or tree2 + return normalize_tree(tree1) == normalize_tree(tree2) + + +def deduplicate_ecms(circuits: Iterable[str], *, simplify: bool = True) -> list[str]: + """Remove duplicate ECM topologies while preserving first occurrences.""" + unique: list[str] = [] + for circuit in circuits: + is_duplicate = any( + compare_ecm_topology(circuit, other, simplify=simplify) for other in unique + ) + if not is_duplicate: + unique.append(circuit) + return unique + + +def topology_to_weighted_vector(tree: CircuitTree | str) -> dict[str, float]: + """Convert an ECM topology tree to a sparse weighted path vector.""" + if isinstance(tree, str): + tree = parse_circuit_tree(tree, preserve_labels=False) + + weights: dict[str, float] = defaultdict(float) + + def visit(node: CircuitTree, path: str = "root", depth: int = 0) -> None: + if isinstance(node, str): + weights[f"{path}:{_component_type(node)}"] += 3 / (1 + depth) + return + + label, children = node + weights[f"{path}:{label}"] += 6 / (1 + depth) + for idx, child in enumerate(children): + visit(child, f"{path}-{label}[{idx}]", depth + 1) + + visit(tree) + return dict(weights) + + +def cosine_similarity(vector1: Mapping[str, float], vector2: Mapping[str, float]) -> float: + """Compute cosine similarity between two sparse vectors.""" + keys = set(vector1) | set(vector2) + dot = sum(vector1.get(key, 0.0) * vector2.get(key, 0.0) for key in keys) + norm1 = math.sqrt(sum(value * value for value in vector1.values())) + norm2 = math.sqrt(sum(value * value for value in vector2.values())) + return dot / (norm1 * norm2) if norm1 and norm2 else 0.0 + + +def structure_similarity(circuit1: str | CircuitTree, circuit2: str | CircuitTree) -> float: + """Compute weighted topology similarity between two ECMs.""" + return cosine_similarity( + topology_to_weighted_vector(circuit1), + topology_to_weighted_vector(circuit2), + ) + + +def analyze_circuit_structure( + tree: CircuitTree | str, + circuit_name: str = "Circuit", +) -> dict[str, Any]: + """Return depth, element, and connection summaries for an ECM topology.""" + if isinstance(tree, str): + tree = parse_circuit_tree(tree, preserve_labels=True) + + info: dict[str, Any] = { + "circuit_name": circuit_name, + "depth_info": {}, + "max_depth": 0, + "total_elements": 0, + "element_types": {}, + "connection_patterns": [], + } + + def visit(node: CircuitTree, depth: int = 0) -> None: + info["max_depth"] = max(info["max_depth"], depth) + info["depth_info"].setdefault( + depth, + {"connections": [], "elements": [], "element_count": 0}, + ) + + if isinstance(node, str): + info["depth_info"][depth]["elements"].append(node) + info["depth_info"][depth]["element_count"] += 1 + info["total_elements"] += 1 + ctype = _component_type(node) + info["element_types"][ctype] = info["element_types"].get(ctype, 0) + 1 + return + + label, children = node + info["depth_info"][depth]["connections"].append(label) + info["connection_patterns"].append(f"Depth{depth}:{label}({len(children)})") + for child in children: + visit(child, depth + 1) + + visit(tree) + return info + + +def compare_circuit_structures( + structure1: Mapping[str, Any], + structure2: Mapping[str, Any], +) -> dict[str, Any]: + """Compare two structure summaries returned by :func:`analyze_circuit_structure`.""" + comparison: dict[str, Any] = { + "depth_comparison": {}, + "element_comparison": {}, + "summary": { + "max_depth_diff": structure1["max_depth"] - structure2["max_depth"], + "total_elements_diff": structure1["total_elements"] - structure2["total_elements"], + }, + } + + element_types = set(structure1["element_types"]) | set(structure2["element_types"]) + for ctype in element_types: + count1 = structure1["element_types"].get(ctype, 0) + count2 = structure2["element_types"].get(ctype, 0) + comparison["element_comparison"][ctype] = { + "structure1": count1, + "structure2": count2, + "diff": count1 - count2, + } + + max_depth = max(structure1["max_depth"], structure2["max_depth"]) + empty = {"connections": [], "elements": [], "element_count": 0} + for depth in range(max_depth + 1): + depth1 = structure1["depth_info"].get(depth, empty) + depth2 = structure2["depth_info"].get(depth, empty) + comparison["depth_comparison"][depth] = { + "connections_1": depth1["connections"], + "connections_2": depth2["connections"], + "elements_1": depth1["elements"], + "elements_2": depth2["elements"], + "element_count_diff": depth1["element_count"] - depth2["element_count"], + "connection_match": set(depth1["connections"]) == set(depth2["connections"]), + "element_match": set(depth1["elements"]) == set(depth2["elements"]), + } + return comparison + + +def ecm_to_vector(circuit: str | CircuitTree, *, simplify: bool = True) -> list[int]: + """Convert an ECM topology to a bracketed integer vector.""" + tree = ( + parse_circuit_tree(circuit, preserve_labels=False) + if isinstance(circuit, str) + else circuit + ) + if simplify: + tree = simplify_tree(tree, preserve_order=True) or tree + + vector: list[int] = [] + + def visit(node: CircuitTree) -> None: + if isinstance(node, tuple): + label, children = node + vector.append(ECM_START_CODES[label]) + for child in children: + visit(child) + vector.append(ECM_END_CODES[label]) + else: + vector.append(ECM_LEAF_CODES[_component_type(node)]) + + visit(tree) + return vector + + +def vector_to_ecm(vector: Sequence[int], *, start_index: int = 1) -> str: + """Reconstruct a labeled ECM string from an integer topology vector.""" + stack: list[tuple[str, list[CircuitTree]]] = [] + roots: list[CircuitTree] = [] + + for code in vector: + code = int(code) + if code in ECM_START_CODES.values(): + stack.append((ECM_CODE_TO_NODE[code], [])) + elif code in ECM_END_CODE_TO_NODE: + expected = ECM_END_CODE_TO_NODE[code] + if not stack or stack[-1][0] != expected: + raise ValueError(f"Unbalanced ECM vector near end code {code}.") + label, children = stack.pop() + node: CircuitTree = (label, children) + if stack: + stack[-1][1].append(node) + else: + roots.append(node) + elif code in ECM_LEAF_CODES.values(): + leaf = ECM_CODE_TO_NODE[code] + if stack: + stack[-1][1].append(leaf) + else: + roots.append(leaf) + else: + raise ValueError(f"Unknown ECM vector code: {code}.") + + if stack: + raise ValueError("Unbalanced ECM vector: some groups were not closed.") + if len(roots) != 1: + raise ValueError("ECM vector must decode to exactly one root node.") + + counters = {key: start_index for key in ECM_LEAF_CODES} + + def label_tree(node: CircuitTree) -> CircuitTree: + if isinstance(node, str): + ctype = _component_type(node) + label = f"{ctype}{counters[ctype]}" + counters[ctype] += 1 + return label + label, children = node + return (label, [label_tree(child) for child in children]) + + return tree_to_circuit(label_tree(roots[0])) + + +def parse_eis_value_tree(expr: str) -> ValueTree: + """Parse a numeric EIS expression such as ``10-[20,30]`` into a value tree.""" + text = expr.replace(" ", "") + idx = 0 + + def parse_number() -> float: + nonlocal idx + match = VALUE_RE.match(text[idx:]) + if not match: + raise ValueError(f"Expected numeric value at position {idx} in {text!r}.") + idx += len(match.group(0)) + return float(match.group(0)) + + def parse_node() -> ValueTree: + nonlocal idx + if idx < len(text) and text[idx] == "[": + return parse_parallel() + return parse_number() + + def parse_series() -> ValueTree: + nonlocal idx + children = [parse_node()] + while idx < len(text) and text[idx] == "-": + idx += 1 + children.append(parse_node()) + return ("series", children) if len(children) > 1 else children[0] + + def parse_parallel() -> ValueTree: + nonlocal idx + idx += 1 + children = [parse_series()] + while idx < len(text) and text[idx] == ",": + idx += 1 + children.append(parse_series()) + if idx >= len(text) or text[idx] != "]": + raise ValueError(f"Unclosed parallel group in {text!r}.") + idx += 1 + return ("parallel", children) + + tree = parse_series() + if idx != len(text): + raise ValueError(f"Unexpected trailing text in EIS value expression: {text[idx:]!r}.") + return tree + + +def eis_to_vector(eis: str | ValueTree) -> list[int | float]: + """Convert a numeric EIS value tree to a vector with value markers.""" + tree = parse_eis_value_tree(eis) if isinstance(eis, str) else eis + vector: list[int | float] = [] + + def visit(node: ValueTree) -> None: + if isinstance(node, tuple): + label, children = node + vector.append(EIS_START_CODES[label]) + for child in children: + visit(child) + vector.append(EIS_END_CODES[label]) + else: + vector.extend([EIS_VALUE_CODE, float(node)]) + + visit(tree) + return vector + + +def extract_eis_values(eis_vector: Sequence[int | float]) -> list[float]: + """Extract numeric leaves from an EIS vector.""" + values: list[float] = [] + idx = 0 + while idx < len(eis_vector): + if int(eis_vector[idx]) == EIS_VALUE_CODE and idx + 1 < len(eis_vector): + values.append(float(eis_vector[idx + 1])) + idx += 2 + else: + idx += 1 + return values + + +def replace_eis_values( + eis_vector: Sequence[int | float], + values: Sequence[float], +) -> list[int | float]: + """Replace numeric leaves in an EIS vector, preserving its topology codes.""" + output = list(eis_vector) + value_idx = 0 + idx = 0 + while idx < len(output) - 1: + if int(output[idx]) == EIS_VALUE_CODE: + if value_idx < len(values): + output[idx + 1] = float(values[value_idx]) + value_idx += 1 + idx += 2 + else: + idx += 1 + return output + + +def _is_number(value: Any) -> bool: + try: + return bool(np.isfinite(float(value))) + except (TypeError, ValueError): + return False + + +def simplify_ecm_eis_pair( + ecm: str | CircuitTree, + eis: str | ValueTree, +) -> tuple[CircuitTree, ValueTree]: + """Simplify matching ECM and numeric EIS trees together for R/C/L runs.""" + ecm_tree = parse_circuit_tree(ecm, preserve_labels=True) if isinstance(ecm, str) else ecm + eis_tree = parse_eis_value_tree(eis) if isinstance(eis, str) else eis + + def combine(label: str, ctype: str, values: list[float]) -> float | None: + if label == "series": + if ctype in {"R", "L"}: + return float(sum(values)) + if ctype == "C": + denom = sum(1 / val for val in values if val != 0) + return float("inf") if denom == 0 else float(1 / denom) + if label == "parallel": + if ctype in {"R", "L"}: + denom = sum(1 / val for val in values if val != 0) + return float("inf") if denom == 0 else float(1 / denom) + if ctype == "C": + return float(sum(values)) + return None + + def simplify_pair( + ecm_node: CircuitTree, + eis_node: ValueTree, + ) -> tuple[CircuitTree, ValueTree, str | None]: + if isinstance(ecm_node, str): + simple_eis = float(eis_node) if _is_number(eis_node) else eis_node + return ecm_node, simple_eis, _component_type(ecm_node) + + label, ecm_children = ecm_node + if isinstance(eis_node, tuple) and eis_node[0] == label: + eis_children = list(eis_node[1]) + else: + eis_children = [eis_node] * len(ecm_children) + + items = [] + for child_ecm, child_eis in zip(ecm_children, eis_children): + simple_ecm, simple_eis, ctype = simplify_pair(child_ecm, child_eis) + items.append( + { + "ecm": simple_ecm, + "eis": simple_eis, + "type": ctype, + "number": float(simple_eis) + if ctype in {"R", "C", "L"} and _is_number(simple_eis) + else None, + } + ) + + merged = [] + idx = 0 + while idx < len(items): + item = items[idx] + ctype = item["type"] + if ctype in {"R", "C", "L"} and item["number"] is not None: + run = [item] + idx += 1 + while ( + idx < len(items) + and items[idx]["type"] == ctype + and items[idx]["number"] is not None + ): + run.append(items[idx]) + idx += 1 + if len(run) > 1: + value = combine(label, ctype, [entry["number"] for entry in run]) + if value is not None and np.isfinite(value): + merged.append({"ecm": run[0]["ecm"], "eis": value, "type": ctype}) + continue + merged.extend(run) + else: + merged.append(item) + idx += 1 + + ecm_out = [entry["ecm"] for entry in merged] + eis_out = [entry["eis"] for entry in merged] + if len(ecm_out) == 1: + ctype = _component_type(ecm_out[0]) if isinstance(ecm_out[0], str) else None + return ecm_out[0], eis_out[0], ctype + return (label, ecm_out), (label, eis_out), None + + simple_ecm, simple_eis, _ = simplify_pair(ecm_tree, eis_tree) + return simple_ecm, simple_eis + + +def ecm_to_eis_string( + ecm: str, + parameter_values: Mapping[str, float], + *, + value_fmt: str = ".6g", +) -> str: + """Replace ECM component labels with numeric values from a parameter mapping.""" + + def replace(match: re.Match[str]) -> str: + component = match.group(0) + if component not in parameter_values: + return component + return format(float(parameter_values[component]), value_fmt) + + return COMPONENT_RE.sub(replace, ecm) + + +def geometric_mean_frequency(freq: Sequence[float]) -> float: + """Return the geometric mean of positive finite frequencies.""" + values = np.asarray(freq, dtype=float) + values = values[np.isfinite(values) & (values > 0)] + if values.size == 0: + raise ValueError("No positive finite frequencies provided.") + return float(np.exp(np.mean(np.log(values)))) + + +def compress_parameter_values( + parameter_values: Mapping[str, float], + f_ref_hz: float, +) -> dict[str, float]: + """Compress parameter labels to one scalar per component for vector workflows.""" + output: dict[str, float] = {} + omega = 2 * math.pi * float(f_ref_hz) + + grouped: dict[str, dict[str, float]] = {} + for name, value in parameter_values.items(): + if re.fullmatch(r"[PW]\d+[wn]", name, flags=re.IGNORECASE): + grouped.setdefault(name[:-1], {})[name[-1].lower()] = float(value) + elif re.fullmatch(r"[RCL]\d+", name, flags=re.IGNORECASE): + output[name.upper()] = float(value) + + for base, values in grouped.items(): + coeff = values.get("w", np.nan) + exponent = values.get("n", np.nan) + base = base.upper() + if not np.isfinite(coeff) or not np.isfinite(exponent): + output[base] = np.nan + elif base.startswith("P"): + output[base] = float(coeff * omega ** (exponent - 1)) + elif base.startswith("W"): + output[base] = float(coeff / omega**exponent) + + return output + + +def order_parameters_by_topology( + topology: str, + scalar_map: Mapping[str, float], +) -> tuple[OrderedDict[str, float], list[float]]: + """Order scalar component values according to left-to-right ECM topology.""" + seen = set() + ordered_items = [] + for component in COMPONENT_RE.findall(topology): + component = component.upper() + if component in scalar_map and component not in seen: + ordered_items.append((component, scalar_map[component])) + seen.add(component) + + for key, value in scalar_map.items(): + key_upper = key.upper() + if key_upper not in seen: + ordered_items.append((key, value)) + + ordered = OrderedDict(ordered_items) + return ordered, [value for _, value in ordered_items] + + +def _bad_parameters_count(value: Any) -> int: + if value is None: + return 0 + if isinstance(value, (list, tuple, set)): + return len(value) + text = str(value).strip() + if text == "" or text.lower() == "none" or text == "[]": + return 0 + if "," in text: + return len([part for part in text.split(",") if part.strip()]) + return 1 + + +def _is_finite_scalar(value: Any) -> bool: + try: + return bool(np.isfinite(float(value))) + except (TypeError, ValueError): + return False + + +def select_best_circuit_row(circuits: pd.DataFrame) -> pd.Series: + """Select the best ranked circuit row using posterior quality and fit metrics.""" + df = circuits.copy() + if "WAIC (real)" in df: + df["WAIC (sum)"] = df["WAIC (real)"] + df["WAIC (imag)"] + elif "WAIC (mag)" in df: + df["WAIC (sum)"] = df["WAIC (mag)"] + df["WAIC (phase)"] + else: + df["WAIC (sum)"] = 0.0 + + bad_parameters = ( + df["Bad Parameters"] + if "Bad Parameters" in df + else pd.Series(["None"] * len(df), index=df.index) + ) + df["_bad_params_count"] = bad_parameters.apply(_bad_parameters_count) + if (df["_bad_params_count"] == 0).any(): + df = df[df["_bad_params_count"] == 0].copy() + + converged = ( + df["converged"] if "converged" in df else pd.Series([True] * len(df), index=df.index) + ) + divergences = ( + df["divergences"] if "divergences" in df else pd.Series([0] * len(df), index=df.index) + ) + r2_real = ( + df["R^2 (ravg)"] if "R^2 (ravg)" in df else pd.Series([0] * len(df), index=df.index) + ) + r2_imag = ( + df["R^2 (iavg)"] if "R^2 (iavg)" in df else pd.Series([0] * len(df), index=df.index) + ) + mape_real = ( + df["MAPE (ravg)"] if "MAPE (ravg)" in df else pd.Series([0] * len(df), index=df.index) + ) + mape_imag = ( + df["MAPE (iavg)"] if "MAPE (iavg)" in df else pd.Series([0] * len(df), index=df.index) + ) + n_params = ( + df["n_params"] if "n_params" in df else pd.Series([0] * len(df), index=df.index) + ) + + df["_not_converged"] = ~converged.astype(bool) + df["_divergences"] = divergences.fillna(0).astype(float) + df["_r2_neg"] = -(r2_real.fillna(0) + r2_imag.fillna(0)) + df["_mape"] = mape_real.fillna(0) + mape_imag.fillna(0) + df["_n_params"] = n_params.fillna(0).astype(float) + + sort_keys = [ + "Posterior Score Rank" if "Posterior Score Rank" in df else "WAIC (sum)", + "WAIC (sum)", + "_not_converged", + "_divergences", + "_bad_params_count", + "_r2_neg", + "_mape", + "_n_params", + ] + df = df.sort_values(sort_keys, ascending=[True] * len(sort_keys), ignore_index=True) + return df.iloc[0] + + +def extract_vectorizable_params( + row: pd.Series, + *, + use_sanitized_for_gaussian: bool = False, + raw_fallback: bool = True, +) -> dict[str, float]: + """Extract one scalar value per posterior parameter from an annotated row.""" + details = row.get("Posterior Details", {}) or {} + means = row.get("Posterior Means", {}) or {} + sanitized = row.get("Sanitized Params", {}) or {} + + output: dict[str, float] = {} + for parameter in sorted(set(details) | set(means) | set(sanitized)): + info = details.get(parameter, {}) + use_sanitized = use_sanitized_for_gaussian or info.get("is_gaussian") is False + value = ( + sanitized.get(parameter, np.nan) + if use_sanitized + else means.get(parameter, np.nan) + ) + if use_sanitized and raw_fallback and not _is_finite_scalar(value): + value = means.get(parameter, np.nan) + output[parameter] = float(value) if value is not None else np.nan + return output + + +def add_ecm_vectors( + circuits: pd.DataFrame, + *, + circuit_column: str = "circuitstring", + vector_column: str = "ECM Vector", + simplify: bool = True, +) -> pd.DataFrame: + """Return a copy of a circuits dataframe with an ECM vector column added.""" + df = circuits.copy() + df[vector_column] = df[circuit_column].apply( + lambda circuit: ecm_to_vector(circuit, simplify=simplify) + ) + return df + + +class ECMVectorTracker: + """Small JSON-backed tracker for ECM topology vectors.""" + + def __init__(self, path: str | Path): + self.path = Path(path) + self.vectors: list[list[int]] = self._load() + + def _load(self) -> list[list[int]]: + if not self.path.exists(): + return [] + return json.loads(self.path.read_text()) + + def save(self) -> None: + self.path.write_text(json.dumps(self.vectors, indent=2)) + + def add(self, circuit: str, *, simplify: bool = True) -> list[int]: + vector = ecm_to_vector(circuit, simplify=simplify) + self.vectors.append(vector) + self.save() + return vector + + def latest_similarity(self) -> float | None: + if len(self.vectors) < 2: + return None + previous = vector_to_ecm(self.vectors[-2]) + current = vector_to_ecm(self.vectors[-1]) + return structure_similarity(previous, current) + + +class EISVectorTracker: + """Small JSON-backed tracker for numeric EIS vectors.""" + + def __init__(self, path: str | Path): + self.path = Path(path) + self.vectors: list[list[int | float]] = self._load() + + def _load(self) -> list[list[int | float]]: + if not self.path.exists(): + return [] + return json.loads(self.path.read_text()) + + def save(self) -> None: + self.path.write_text(json.dumps(self.vectors, indent=2)) + + def add(self, eis: str | ValueTree) -> list[int | float]: + vector = eis_to_vector(eis) + self.vectors.append(vector) + self.save() + return vector diff --git a/src/autoeis/utils.py b/src/autoeis/utils.py index 5403e0f..58f0b59 100644 --- a/src/autoeis/utils.py +++ b/src/autoeis/utils.py @@ -1310,6 +1310,16 @@ def print_summary(self): """Prints a summary of the inference results.""" self.mcmc.print_summary() + def assess_posterior_quality(self, **kwargs): + """Assess posterior parameter distributions from this inference result. + + See :func:`autoeis.diagnostics.assess_posterior_quality` for keyword + arguments such as ``threshold``, ``var_names``, and ``save_plots``. + """ + from .diagnostics import assess_posterior_quality + + return assess_posterior_quality(self, **kwargs) + class ImpedanceData: """Container for impedance data.""" diff --git a/tests/test_diagnostics.py b/tests/test_diagnostics.py new file mode 100644 index 0000000..3f157c5 --- /dev/null +++ b/tests/test_diagnostics.py @@ -0,0 +1,106 @@ +import numpy as np +import pandas as pd + +import autoeis as ae + + +def test_assess_posterior_sample_classifies_normal_distribution(): + rng = np.random.default_rng(0) + sample = rng.normal(loc=10, scale=2, size=2000) + + result = ae.diagnostics.assess_posterior_sample(sample, "R1", threshold=75) + + assert result.parameter == "R1" + assert result.shape == "Normal" + assert result.status == "Good" + assert result.r2 > 0.99 + assert result.modes == 1 + + +def test_assess_posterior_sample_detects_multipeak_distribution(): + rng = np.random.default_rng(1) + sample = np.concatenate( + [ + rng.normal(loc=-3, scale=0.25, size=1000), + rng.normal(loc=3, scale=0.25, size=1000), + ] + ) + + result = ae.diagnostics.assess_posterior_sample(sample, "P2w", threshold=75) + + assert result.shape == "Multi-Peak" + assert result.status == "Bad" + assert result.modes == 2 + + +def test_assess_posterior_quality_accepts_sample_mapping_and_skips_auxiliary(): + rng = np.random.default_rng(2) + samples = { + "R1": rng.normal(size=500), + "P2n": rng.lognormal(mean=0, sigma=0.8, size=500), + "lp__": rng.normal(size=500), + } + + diagnostics = ae.diagnostics.assess_posterior_quality(samples, threshold=75) + + assert list(diagnostics["parameter"]) == ["R1", "P2n"] + assert "lp__" not in diagnostics["parameter"].to_list() + assert set( + [ + "parameter", + "score", + "r2", + "skewness", + "modes", + "shape", + "status", + "good", + ] + ).issubset(diagnostics.columns) + + +def test_plot_posterior_quality_saves_png(tmp_path): + rng = np.random.default_rng(3) + sample = rng.normal(size=500) + + fig, _ = ae.diagnostics.plot_posterior_quality( + sample, parameter="sigma.mag", save_dir=tmp_path + ) + + try: + saved = list(tmp_path.glob("sigma.mag-*.png")) + assert len(saved) == 1 + finally: + fig.clf() + + +def test_add_posterior_quality_updates_and_ranks_circuits_dataframe(): + class FakeInferenceResult: + def __init__(self, circuit, samples): + self.circuit = circuit + self.samples = samples + self.converged = True + + rng = np.random.default_rng(4) + results = [ + FakeInferenceResult("R1", {"R1": rng.normal(size=1000)}), + FakeInferenceResult( + "R1-[R2,P3]", + { + "R1": np.concatenate( + [ + rng.normal(loc=-3, scale=0.2, size=500), + rng.normal(loc=3, scale=0.2, size=500), + ] + ) + }, + ), + ] + circuits = pd.DataFrame({"circuitstring": ["R1", "R1-[R2,P3]"]}) + + ranked = ae.diagnostics.add_posterior_quality(circuits, results, threshold=75) + + assert "Posterior Score" in ranked + assert "Posterior Quality" in ranked + assert list(ranked["Posterior Score Rank"]) == [1, 2] + assert ranked.iloc[0]["Bad Parameters"] == "None" diff --git a/tests/test_ecm.py b/tests/test_ecm.py new file mode 100644 index 0000000..1c9a101 --- /dev/null +++ b/tests/test_ecm.py @@ -0,0 +1,77 @@ +import numpy as np +import pandas as pd + +import autoeis as ae + + +def test_parse_normalize_and_compare_ecm_topology(): + tree = ae.ecm.parse_circuit_tree("R1-[P2,R3]", preserve_labels=False) + + assert tree == ("series", ["R", ("parallel", ["P", "R"])]) + assert ae.ecm.compare_ecm_topology("R1-C2-L3", "L9-R8-C7") + assert not ae.ecm.compare_ecm_topology("R1-C2-L3", "R1-C2-P3") + + +def test_simplify_ecm_preserves_representative_labels(): + simplified = ae.ecm.simplify_ecm("R1-R2-[C1-C2,R3]-R4") + + assert simplified == "R1-[C1,R3]" + + +def test_deduplicate_ecms_preserves_first_occurrence(): + circuits = ["R1-C2", "C9-R8", "R1-[P2,R3]", "R4-[R5,P6]"] + + unique = ae.ecm.deduplicate_ecms(circuits) + + assert unique == ["R1-C2", "R1-[P2,R3]"] + + +def test_ecm_vectorization_roundtrip(): + vector = ae.ecm.ecm_to_vector("R1-[P2,R3]", simplify=False) + + assert vector == [0, 2, 1, 5, 2, -3, -2] + assert ae.ecm.vector_to_ecm(vector) == "R1-[P1,R2]" + + +def test_eis_vectorization_and_value_replacement(): + vector = ae.ecm.eis_to_vector("10-[2.5,30]-1e-6") + + assert vector == [0, 2, 10.0, 1, 2, 2.5, 2, 30.0, -3, 2, 1e-6, -2] + assert ae.ecm.extract_eis_values(vector) == [10.0, 2.5, 30.0, 1e-6] + replaced = ae.ecm.replace_eis_values(vector, [1, 2, 3, 4]) + assert ae.ecm.extract_eis_values(replaced) == [1.0, 2.0, 3.0, 4.0] + + +def test_simplify_matching_ecm_and_eis_pair(): + ecm_tree, eis_tree = ae.ecm.simplify_ecm_eis_pair("R1-R2-[C1,C2]", "10-20-[2,3]") + + assert ae.ecm.tree_to_circuit(ecm_tree) == "R1-C1" + assert ae.ecm.tree_to_circuit(eis_tree) == "30-5" + + +def test_ecm_to_eis_string_and_parameter_scalar_ordering(): + parameter_values = {"R1": 10, "P2": 0.5, "R3": 30} + + eis_string = ae.ecm.ecm_to_eis_string("R1-[P2,R3]", parameter_values) + + assert eis_string == "10-[0.5,30]" + + +def test_compress_and_order_parameter_values(): + f_ref = 1 / (2 * np.pi) + params = {"R1": 10, "P2w": 4, "P2n": 1, "C3": 1e-6} + + compressed = ae.ecm.compress_parameter_values(params, f_ref) + ordered, vector = ae.ecm.order_parameters_by_topology("R1-[P2,C3]", compressed) + + assert compressed["P2"] == 4 + assert list(ordered) == ["R1", "P2", "C3"] + assert vector == [10.0, 4.0, 1e-6] + + +def test_add_ecm_vectors_to_dataframe(): + circuits = pd.DataFrame({"circuitstring": ["R1-[P2,R3]"]}) + + out = ae.ecm.add_ecm_vectors(circuits, simplify=False) + + assert out["ECM Vector"].iloc[0] == [0, 2, 1, 5, 2, -3, -2]