Reject bad EIS fits

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

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]

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()

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