diff --git a/pyflowbatt/analysis.py b/pyflowbatt/analysis.py index 885893a..2ffd7be 100644 --- a/pyflowbatt/analysis.py +++ b/pyflowbatt/analysis.py @@ -914,12 +914,17 @@ def analyse_sample( try: df = read_to_bdf(f) params, Z_fit = eis.analyse( - df, label=f.stem, model=config.eis_model, guess_init=has_ml_guess + df, + label=f.stem, + model=config.eis_model, + guess_init=has_ml_guess, + max_fit_error=config.eis_max_fit_error, ) - fig, _ax = eis.plot(df, Z_fit) + fig, _ax = eis.plot(df, Z_fit, max_fit_error=config.eis_max_fit_error) fig.savefig(results_dir / f"eis_{tag}.png") plt.close(fig) - eis_res[tag] = params + if params is not None: + eis_res[tag] = params df["Real Impedance Fit / ohm"] = np.real(Z_fit) df["Imaginary Impedance Fit / ohm"] = np.imag(Z_fit) df_save_bdf(df, results_dir / f"eis_{tag}.x", save_format=save_format) @@ -929,7 +934,7 @@ def analyse_sample( (results_dir / f"eis_{tag}.x").with_suffix(bdf_suffix) ) tracked_mpr_inputs[f"eis_{tag}"] = [Path(eis_file)] - for name, values in params.items(): + for name, values in (params or {}).items(): rows.append({"file": f.stem, "tag": tag, "name": name, **values}) except Exception as e: logger.error("Failed to fit %s: %s", f.stem, str(e)) # noqa: TRY400 diff --git a/pyflowbatt/config.py b/pyflowbatt/config.py index 0bfe8d2..0172625 100644 --- a/pyflowbatt/config.py +++ b/pyflowbatt/config.py @@ -277,6 +277,7 @@ class PyFlowBattConfig: [eis] model = "L0-R0-(R1,CPE1)-(R2,CPE2)" # any fasteis circuit string summary_params = ["R0.r", "R1.r"] # fit parameters reported in the summary + max_fit_error = 0.1 # fit-quality gate above which parameters aren't reported [cv] v_min = 0.4 # voltage window used to plot @@ -308,6 +309,7 @@ class PyFlowBattConfig: eis_tag_patterns: dict[str, list[str]] = field(default_factory=dict) eis_model: str = "L0-R0-(R1,CPE1)-(R2,CPE2)" # fasteis circuit string fitted to EIS eis_summary_params: list[str] = field(default_factory=lambda: ["R0.r", "R1.r", "R2.r"]) + eis_max_fit_error: float = 0.1 # RMS fit error relative to |Z| above which EIS isn't reported sample_name_pattern: str = r"^\d+_.+_.+$" sample_name: str | None = None # overrides pattern entirely lsv_threshold: int = 8 # numeric cutoff for pre/post when only one LSV file is found @@ -409,11 +411,19 @@ def to_resolved_toml(self) -> str: for toml_key, attr in _TECHNIQUE_KEYS.items() }, "eis_tags": dict(self.eis_tag_patterns), - "eis": {"model": self.eis_model, "summary_params": self.eis_summary_params}, + "eis": { + "model": self.eis_model, + "summary_params": self.eis_summary_params, + "max_fit_error": self.eis_max_fit_error, + }, "cv": {toml_key: getattr(self, attr) for toml_key, attr in _CV_KEYS.items()}, } section_defaults: dict[str, dict[str, object]] = { - "eis": {"model": defaults.eis_model, "summary_params": defaults.eis_summary_params}, + "eis": { + "model": defaults.eis_model, + "summary_params": defaults.eis_summary_params, + "max_fit_error": defaults.eis_max_fit_error, + }, "cv": {toml_key: getattr(defaults, attr) for toml_key, attr in _CV_KEYS.items()}, } for name, table in sections.items(): @@ -481,6 +491,11 @@ def _apply_toml(self, path: Path, data: dict, text: str) -> None: logger.warning( "Ignoring bad eis.summary_params in %s (expected list of strings)", path ) + if "max_fit_error" in eis: + if _is_number(eis["max_fit_error"]): + self.eis_max_fit_error = float(eis["max_fit_error"]) + else: + logger.warning("Ignoring bad eis.max_fit_error in %s (expected number)", path) for key in ("extensions", "extra_extensions", "include_files"): if key in data: @@ -619,9 +634,12 @@ def _apply_toml(self, path: Path, data: dict, text: str) -> None: ## See https://empaeconversion.github.io/fasteis/models/ ## summary_params lists the fit parameters reported in the summary for each EIS tag, ## named as in the EIS sheet, e.g. ["R0.r", "R1.r", "CPE1.alpha"]. +## max_fit_error is the RMS distance between fit and data relative to |Z|, e.g. 0.1 is +## 10% off on average. Fits worse than this are plotted but not reported. # [eis] # model = "L0-R0-(R1,CPE1)-(R2,CPE2)" -# summary_params = ["R0.r"] +# summary_params = ["R0.r", "R1.r", "R2.r"] +# max_fit_error = 0.1 ## --- CV analysis parameters --- # [cv] diff --git a/pyflowbatt/eis.py b/pyflowbatt/eis.py index ff8f966..1d4cad7 100644 --- a/pyflowbatt/eis.py +++ b/pyflowbatt/eis.py @@ -16,6 +16,7 @@ DEFAULT_MODEL = "L0-R0-(R1,CPE1)-(R2,CPE2)" MODELS_URL = "https://empaeconversion.github.io/fasteis/models/" +MAX_FIT_ERROR = 0.1 # fasteis unit strings -> units for summary row labels PRETTY_UNITS: dict[str, str] = { @@ -99,18 +100,33 @@ def check_model(model: str, summary_params: list[str]) -> tuple[dict[str, str], return units, True +def _impedance(df: pd.DataFrame) -> np.ndarray: + """Measured complex impedance in ohm.""" + return (df["Real Impedance / ohm"] + 1j * df["Imaginary Impedance / ohm"]).to_numpy() + + +def fit_error(df: pd.DataFrame, Z_fit: np.ndarray) -> float: + """RMS distance between fit and data relative to |Z|, e.g. 0.05 is 5% off on average.""" + Z = _impedance(df) + return float(np.sqrt(np.mean(np.abs((Z_fit - Z) / np.abs(Z)) ** 2))) + + def analyse( df: pd.DataFrame, label: str = "", model: str = DEFAULT_MODEL, *, guess_init: bool | None = None, -) -> tuple[dict[str, dict], np.ndarray]: - """Fit EIS to an equivalent circuit model, warning if the fit looks unreliable.""" + max_fit_error: float = MAX_FIT_ERROR, +) -> tuple[dict[str, dict] | None, np.ndarray]: + """Fit EIS to an equivalent circuit model, warning if the fit looks unreliable. + + Parameters are None if the fit error is above max_fit_error. + """ import fasteis # noqa: PLC0415 f = df["Frequency / Hz"] - Z = df["Real Impedance / ohm"] + 1j * df["Imaginary Impedance / ohm"] + Z = _impedance(df) circuit = fasteis.Circuit(model) res = circuit.fit(f, Z, guess_init=guess_init) Z_fit = res.circuit.impedance(f) @@ -123,23 +139,44 @@ def analyse( logger.warning("- EIS fit %s did not converge", label) if all(p["err"] is None for p in params.values()): logger.warning("- EIS fit %s has no standard errors, parameters are degenerate", label) + error = fit_error(df, Z_fit) + if error > max_fit_error: + logger.warning( + "- EIS fit %s error is very high (%.3f>%s), not reporting parameters", + label, + error, + max_fit_error, + ) + return None, Z_fit return params, Z_fit -def plot(df: pd.DataFrame, Z_fit: np.ndarray) -> tuple[Figure, Axes]: +def plot( + df: pd.DataFrame, Z_fit: np.ndarray, max_fit_error: float = MAX_FIT_ERROR +) -> tuple[Figure, Axes]: """Nyquist plot fit result.""" - Z = df["Real Impedance / ohm"] + 1j * df["Imaginary Impedance / ohm"] + Z = _impedance(df) + error = fit_error(df, Z_fit) + bad = error > max_fit_error fig, ax = plt.subplots() # Residual segments from each data point to its fitted point data_pts = np.column_stack([np.real(Z), -np.imag(Z)]) fit_pts = np.column_stack([np.real(Z_fit), -np.imag(Z_fit)]) segments = list(zip(data_pts, fit_pts, strict=True)) ax.plot(np.real(Z), -np.imag(Z), "k.-", label="Data") - ax.plot(np.real(Z_fit), -np.imag(Z_fit), "C0.-", label="Fit") + ax.plot( + np.real(Z_fit), + -np.imag(Z_fit), + "C0.--" if bad else "C0.-", + label="Fit (bad)" if bad else "Fit", + ) ax.add_collection( LineCollection(segments, colors="grey", linewidths=0.8, label="Residual", zorder=1) ) - ax.legend() + error_text = f"Fit error = {error:.3f}" + if bad: + error_text += f"\n(>{max_fit_error}, result ignored)" + ax.legend(title=error_text) ax.set_xlabel("Real Impedance / ohm") ax.set_ylabel("Imaginary Impedance / ohm") fig.tight_layout() diff --git a/tests/test_config.py b/tests/test_config.py index 7205ad6..80a2d12 100644 --- a/tests/test_config.py +++ b/tests/test_config.py @@ -1,6 +1,7 @@ """Tests for pyflowbatt.toml configuration and technique file classification.""" import json +import logging import re import zipfile from pathlib import Path @@ -2460,7 +2461,7 @@ def fake_analyse(_df: object, label: str, **_kwargs: object) -> tuple[dict, np.n } return params, np.array([1 + 1j]) - def fake_plot(_df: object, _z_fit: np.ndarray) -> tuple: + def fake_plot(_df: object, _z_fit: np.ndarray, **_kwargs: object) -> tuple: import matplotlib.pyplot as plt return plt.subplots() @@ -2498,6 +2499,108 @@ def fake_check_model(_model: str, _params: list[str]) -> tuple[dict[str, str], b assert np.isnan(df.loc["EIS R0.r pre-50%SOC / Ω", "Value"]) +def test_analyse_sample_leaves_out_rejected_eis_fit( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + """A fit rejected for its fit error is left out of the summary and EIS sheet but plotted.""" + import numpy as np + import pandas as pd + + from pyflowbatt import analysis as analysis_module + + sample = _write_two_eis_sample(tmp_path / "sample") + (sample / "pyflowbatt.toml").write_text("[eis]\nmax_fit_error = 0.05\n") + max_errors = [] + + def fake_analyse( + _df: object, label: str, max_fit_error: float, **_kwargs: object + ) -> tuple[dict | None, np.ndarray]: + max_errors.append(max_fit_error) + if label == "s_06_PEIS_A": + return None, np.array([1 + 1j]) + return {"R0.r": {"value": 1.0, "err": 0.1, "unit": "ohm"}}, np.array([1 + 1j]) + + def fake_plot(_df: object, _z_fit: np.ndarray, max_fit_error: float) -> tuple: + import matplotlib.pyplot as plt + + max_errors.append(max_fit_error) + return plt.subplots() + + monkeypatch.setattr(analysis_module, "read_to_bdf", lambda _f: pd.DataFrame({"x": [0]})) + monkeypatch.setattr(analysis_module, "df_save_bdf", lambda *_a, **_k: None) + monkeypatch.setattr(analysis_module.eis, "analyse", fake_analyse) + monkeypatch.setattr(analysis_module.eis, "plot", fake_plot) + monkeypatch.setattr(analysis_module.eis, "check_model", lambda *_a: ({"R0.r": "ohm"}, True)) + + config = analysis_module.PyFlowBattConfig.load(sample, home=tmp_path / "home") + analysis_module.analyse_sample(sample, save_format=None, config=config) + + assert max_errors == [0.05] * 4 + assert (sample / "results" / "eis_post.png").exists() + xlsx = sample / "results" / "summary.xlsx" + summary = pd.read_excel(xlsx, sheet_name="Summary", index_col=0) + assert summary.loc["EIS R0.r pre / Ω", "Value"] == 1.0 + assert np.isnan(summary.loc["EIS R0.r post / Ω", "Value"]) + assert list(pd.read_excel(xlsx, sheet_name="EIS")["tag"]) == ["pre"] + + +def test_eis_fit_error_relative_to_modulus() -> None: + """Fit error is the RMS of |Z_fit - Z| / |Z|, so it doesn't scale with point count.""" + import numpy as np + import pandas as pd + + from pyflowbatt import eis + + for n in (10, 100): + Z = np.full(n, 3 + 4j) + df = pd.DataFrame({"Real Impedance / ohm": Z.real, "Imaginary Impedance / ohm": Z.imag}) + assert eis.fit_error(df, Z) == 0 + assert eis.fit_error(df, Z + 0.5j) == pytest.approx(0.1) + + +def test_eis_analyse_rejects_bad_fit(caplog: pytest.LogCaptureFixture) -> None: + """A fit with error above max_fit_error warns and returns no parameters.""" + import fasteis + import numpy as np + import pandas as pd + + from pyflowbatt import eis + + f = np.logspace(5, -1, 40) + Z = fasteis.Circuit("R0-(R1,C1)").with_values([1.0, 10.0, 1e-2]).impedance(f) + df = pd.DataFrame( + {"Frequency / Hz": f, "Real Impedance / ohm": Z.real, "Imaginary Impedance / ohm": Z.imag} + ) + params, Z_fit = eis.analyse(df, model="R0-(R1,C1)", guess_init=False) + assert params is not None + assert eis.fit_error(df, Z_fit) < 1e-3 + + with caplog.at_level(logging.WARNING): + params, Z_fit = eis.analyse(df, label="flat", model="R0", guess_init=False) + assert params is None + assert len(Z_fit) == len(f) + assert "EIS fit flat error is very high" in caplog.text + + +def test_eis_max_fit_error_toml_override(tmp_path: Path) -> None: + """[eis] max_fit_error overrides the default, and a non-number is ignored.""" + from pyflowbatt import eis + + sample_dir = tmp_path / "sample" + sample_dir.mkdir() + home = tmp_path / "home" + toml = sample_dir / "pyflowbatt.toml" + assert PyFlowBattConfig.load(sample_dir, home=home).eis_max_fit_error == eis.MAX_FIT_ERROR + + toml.write_text("[eis]\nmax_fit_error = 0.3\n") + config = PyFlowBattConfig.load(sample_dir, home=home) + assert config.eis_max_fit_error == 0.3 + assert "max_fit_error = 0.3" in config.to_resolved_toml() + + toml.write_text('[eis]\nmax_fit_error = "low"\n') + assert PyFlowBattConfig.load(sample_dir, home=home).eis_max_fit_error == eis.MAX_FIT_ERROR + + def test_custom_eis_tag_descriptions() -> None: """Custom EIS tags get descriptions naming the tag, standard ones keep theirs.""" from pyflowbatt import analysis as analysis_module