From 4260e98d3fe7d3586415070d0e1c0dbbb73aed7a Mon Sep 17 00:00:00 2001 From: Defne Nihal Ertugrul Date: Thu, 9 Apr 2026 20:14:26 +0300 Subject: [PATCH 1/6] fix(constants): use math.pi instead of truncated 3.14159 for column cross-section MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit PlantGeometry.COLUMN_CROSS_SECTION_M2 was computed with the literal 3.14159, which truncates pi at the 6th significant figure. For the 0.1 m absorber column the resulting area is off by ~8e-9 m² — small in absolute terms but it propagates into every K_OGa and H_OG calculation that uses the default geometry, and the discrepancy is visible against any independent re-derivation. Use math.pi so the constant is exact to double precision. Co-Authored-By: Claude Opus 4.6 (1M context) --- src/meapy/constants.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/meapy/constants.py b/src/meapy/constants.py index bd52245..0681d11 100644 --- a/src/meapy/constants.py +++ b/src/meapy/constants.py @@ -15,6 +15,8 @@ from __future__ import annotations +import math + __all__ = [ "MEAProperties", "PlantGeometry", @@ -81,7 +83,7 @@ class PlantGeometry: COLUMN_HEIGHT_M: float = 6.2 # m COLUMN_DIAMETER_M: float = 0.1 # m (100 mm ID) - COLUMN_CROSS_SECTION_M2: float = 3.14159 * (0.1 / 2) ** 2 # π r² ≈ 7.854e-3 m² + COLUMN_CROSS_SECTION_M2: float = math.pi * (0.1 / 2) ** 2 # π r² ≈ 7.854e-3 m² SAMPLING_HEIGHTS_M: list[float] = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5] # ports 1–6 C100_PLATE_AREA_M2: float = 0.30 # m² (manufacturer spec) C200_PLATE_AREA_M2: float = 0.25 # m² From 8c3fccccf5bc59598a06fc78646c0ebacdad9dca Mon Sep 17 00:00:00 2001 From: Defne Nihal Ertugrul Date: Thu, 9 Apr 2026 20:32:38 +0300 Subject: [PATCH 2/6] refactor(constants): freeze class-level mutable defaults MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit PlantGeometry.SAMPLING_HEIGHTS_M, ExperimentLabels.LABELS, and ExperimentLabels.COOLING_WATER_FLOWRATES_KG_H were defined as plain list/dict attributes on a class with no instances. Because Python class attributes are shared, any caller that mutated these (e.g. SAMPLING_HEIGHTS_M.append(...) or COOLING_WATER_FLOWRATES_KG_H["A"] = 0) would silently corrupt the value for every other consumer in the same process. The risk is real for notebook workflows where the same kernel is reused across analyses. Switch the lists to tuples and wrap the dict in types.MappingProxyType so mutations raise TypeError immediately. The public read API is unchanged — indexing, iteration, len(), and `in` all behave identically. Co-Authored-By: Claude Opus 4.6 (1M context) --- src/meapy/constants.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/src/meapy/constants.py b/src/meapy/constants.py index 0681d11..454e9a3 100644 --- a/src/meapy/constants.py +++ b/src/meapy/constants.py @@ -16,6 +16,8 @@ from __future__ import annotations import math +from types import MappingProxyType +from typing import Mapping __all__ = [ "MEAProperties", @@ -84,7 +86,7 @@ class PlantGeometry: COLUMN_HEIGHT_M: float = 6.2 # m COLUMN_DIAMETER_M: float = 0.1 # m (100 mm ID) COLUMN_CROSS_SECTION_M2: float = math.pi * (0.1 / 2) ** 2 # π r² ≈ 7.854e-3 m² - SAMPLING_HEIGHTS_M: list[float] = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5] # ports 1–6 + SAMPLING_HEIGHTS_M: tuple[float, ...] = (0.5, 1.5, 2.5, 3.5, 4.5, 5.5) # ports 1–6 C100_PLATE_AREA_M2: float = 0.30 # m² (manufacturer spec) C200_PLATE_AREA_M2: float = 0.25 # m² PACKING_TYPE: str = "Sulzer MellapakPlus 252.Y" @@ -136,11 +138,13 @@ class ExperimentLabels: experiment in kg/h, as reported in Table 5.2.1 of the main report. """ - LABELS: list[str] = ["A", "B", "C", "D", "E"] - COOLING_WATER_FLOWRATES_KG_H: dict[str, float] = { - "A": 722.0, - "B": 719.0, - "C": 716.0, - "D": 541.0, - "E": 1615.0, - } + LABELS: tuple[str, ...] = ("A", "B", "C", "D", "E") + COOLING_WATER_FLOWRATES_KG_H: Mapping[str, float] = MappingProxyType( + { + "A": 722.0, + "B": 719.0, + "C": 716.0, + "D": 541.0, + "E": 1615.0, + } + ) From 04a9a84b5c6c73ecf56ec61c3445b1c7f5d8064e Mon Sep 17 00:00:00 2001 From: Defne Nihal Ertugrul Date: Thu, 9 Apr 2026 20:33:53 +0300 Subject: [PATCH 3/6] refactor(heat_transfer): warn when thermal efficiency exceeds unity MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The efficiency() docstring already states that η > 1 indicates measurement inconsistency (instrumentation drift, heat gain, or a swapped stream assignment), but the function never surfaced this — it returned the value silently and downstream U / ε / NTU calculations propagated the inconsistency. Emit a logger.warning with the offending duties and concrete remediation hints. Behaviour is unchanged for callers that don't configure logging (the package-level NullHandler swallows it), so this is non-breaking for existing scripts but immediately visible to anyone running with logging.basicConfig(level=logging.WARNING) or structured-log handlers. Co-Authored-By: Claude Opus 4.6 (1M context) --- src/meapy/heat_transfer.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/meapy/heat_transfer.py b/src/meapy/heat_transfer.py index 29fee7a..5f0147a 100644 --- a/src/meapy/heat_transfer.py +++ b/src/meapy/heat_transfer.py @@ -302,6 +302,16 @@ def efficiency(q_cold_w: float, q_hot_w: float) -> float: raise ValueError("q_hot_w must not be zero; efficiency is undefined.") eta = q_cold_w / abs(q_hot_w) logger.debug("efficiency η = %.4f", eta) + if eta > 1.0: + logger.warning( + "Thermal efficiency η = %.4f exceeds unity (Q_cold=%.2f W, |Q_hot|=%.2f W). " + "This indicates a measurement inconsistency: instrumentation drift, heat " + "gain from surroundings, or a swapped hot/cold stream assignment. Verify " + "transmitter calibration before trusting downstream U/ε/NTU values.", + eta, + q_cold_w, + abs(q_hot_w), + ) return eta From b9f319f56024ae22fba53ad3fe23c30cfee7e95f Mon Sep 17 00:00:00 2001 From: Defne Nihal Ertugrul Date: Thu, 9 Apr 2026 20:37:15 +0300 Subject: [PATCH 4/6] =?UTF-8?q?perf(utils):=20vectorise=20rolling=5Fmean?= =?UTF-8?q?=20with=20cumsum=20(O(n=C2=B7w)=20=E2=86=92=20O(n))?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The previous implementation called np.mean inside a Python for-loop, paying NumPy's per-call overhead once per element. cProfile against a 100k-sample input showed 1.2M function calls and 222 ms wall time — completely dominated by 100k np.mean dispatches. Replace the loop with the standard cumulative-sum trick: build a padded cumsum once, then index it twice to get every trailing-window sum in a single vectorised expression. The leading edge (where a full window is not yet available) is handled by counting actual window length per index, which preserves the previous behaviour exactly (verified by the existing 19 unit tests, all passing). Benchmark (Apple M-series, Python 3.13.5, NumPy 2.x, window=60, best-of-5): size before after speedup ----------- ------------ ----------- -------- 1,000 ~2.2 ms 0.013 ms ~170x 10,000 27.2 ms 0.078 ms 349x 100,000 222.5 ms 0.893 ms 249x 1,000,000 2,223.3 ms 11.69 ms 190x cProfile call count for size=100k drops from 1,200,006 to 11. Also adds benchmarks/bench_rolling_mean.py — a reproducible script with the exact baseline + cProfile dump so future regressions are easy to spot. Co-Authored-By: Claude Opus 4.6 (1M context) --- benchmarks/bench_rolling_mean.py | 67 ++++++++++++++++++++++++++++++++ src/meapy/utils.py | 19 ++++++--- 2 files changed, 81 insertions(+), 5 deletions(-) create mode 100644 benchmarks/bench_rolling_mean.py diff --git a/benchmarks/bench_rolling_mean.py b/benchmarks/bench_rolling_mean.py new file mode 100644 index 0000000..272f10b --- /dev/null +++ b/benchmarks/bench_rolling_mean.py @@ -0,0 +1,67 @@ +"""Micro-benchmark for meapy.utils.rolling_mean. + +Run with: + python benchmarks/bench_rolling_mean.py + +Profiles the rolling-mean computation across array sizes representative of +plant data: 1k samples (a single short steady-state run), 10k (a full hour +at 1 Hz), 100k (a multi-hour campaign), and 1M (a week of 1 Hz logging). + +Baseline numbers captured on Apple M-series, Python 3.13, NumPy 2.x. +See the commit message of the perf optimisation for before/after deltas. +""" + +from __future__ import annotations + +import cProfile +import pstats +import time +from io import StringIO + +import numpy as np + +from meapy.utils import rolling_mean + +SIZES = (1_000, 10_000, 100_000, 1_000_000) +WINDOW = 60 # 1-minute trailing window at 1 Hz +REPEATS = 5 + + +def time_one(size: int, window: int, repeats: int) -> float: + """Return best-of-N wall-clock time for a single rolling_mean call.""" + arr = np.random.default_rng(0).standard_normal(size) + best = float("inf") + for _ in range(repeats): + t0 = time.perf_counter() + rolling_mean(arr, window) + dt = time.perf_counter() - t0 + if dt < best: + best = dt + return best + + +def main() -> None: + print(f"rolling_mean benchmark — window={WINDOW}, best-of-{REPEATS}") + print("-" * 56) + print(f"{'size':>10} {'time (ms)':>12} {'throughput (M-elem/s)':>22}") + print("-" * 56) + for size in SIZES: + t = time_one(size, WINDOW, REPEATS) + thr = size / t / 1e6 + print(f"{size:>10} {t * 1e3:>12.3f} {thr:>22.2f}") + print("-" * 56) + + # cProfile dump for the largest size — useful for spotting hotspots. + arr = np.random.default_rng(0).standard_normal(100_000) + pr = cProfile.Profile() + pr.enable() + rolling_mean(arr, WINDOW) + pr.disable() + buf = StringIO() + pstats.Stats(pr, stream=buf).sort_stats("cumulative").print_stats(10) + print("\ncProfile (size=100_000, top 10 by cumulative):") + print(buf.getvalue()) + + +if __name__ == "__main__": + main() diff --git a/src/meapy/utils.py b/src/meapy/utils.py index 6e41b42..03cf0c8 100644 --- a/src/meapy/utils.py +++ b/src/meapy/utils.py @@ -211,8 +211,17 @@ def rolling_mean(arr: npt.ArrayLike, window: int) -> npt.NDArray[np.float64]: if window < 1: raise ValueError(f"window must be ≥ 1, got {window!r}.") a = np.asarray(arr, dtype=float).ravel() - result = np.empty_like(a) - for i in range(len(a)): - start = max(0, i - window + 1) - result[i] = np.mean(a[start : i + 1]) - return result + n = a.size + if n == 0: + return a + # Cumulative-sum trick: O(n) instead of the O(n·window) Python loop. + # csum[k] = sum(a[:k]); the trailing-window sum at index i is + # csum[i+1] - csum[max(0, i+1-window)], divided by the actual window + # length (which grows from 1 up to `window` over the leading edge). + csum = np.empty(n + 1, dtype=float) + csum[0] = 0.0 + np.cumsum(a, out=csum[1:]) + idx = np.arange(1, n + 1) + starts = np.maximum(0, idx - window) + counts = idx - starts # same as np.minimum(idx, window) + return (csum[idx] - csum[starts]) / counts From e13ebad32b11669b0f1c4ad3b332bee0fb7f1434 Mon Sep 17 00:00:00 2001 From: Defne Nihal Ertugrul Date: Thu, 9 Apr 2026 20:41:03 +0300 Subject: [PATCH 5/6] fix(mass_transfer): deprecate silently-ignored m_slope arg in ntu_og (#5) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit ntu_og(y_bottom, y_top, m_slope=...) accepted m_slope and the docstring promised "a rigorous calculation including the equilibrium back-pressure" when a non-zero value was passed. The function body never references m_slope at any point — every call returns the dilute-limit ln(Y_b / Y_t) regardless. Callers running the supposed "rigorous mode" were silently getting the dilute approximation, an optimistic NTU on real plant data. This is the v0.1.x semver-clean fix: * Non-zero m_slope now emits DeprecationWarning pointing at issue #5 and documenting that the parameter will be removed in v0.2.0. * The default m_slope=0.0 path stays warning-free, so existing scripts and the existing test suite are unaffected. * Two new tests pin the behaviour: - non-zero m_slope warns and returns the same value as m_slope=0 - default-call path is warning-free under warnings.simplefilter("error") A proper Colburn-equation entry point taking the absorption factor explicitly will land in v0.2.0 (tracked in #5). Fixes #5 Co-Authored-By: Claude Opus 4.6 (1M context) --- src/meapy/mass_transfer.py | 27 ++++++++++++++++++++++----- tests/unit/test_mass_transfer.py | 15 +++++++++++++++ 2 files changed, 37 insertions(+), 5 deletions(-) diff --git a/src/meapy/mass_transfer.py b/src/meapy/mass_transfer.py index 4144e76..7f1022f 100644 --- a/src/meapy/mass_transfer.py +++ b/src/meapy/mass_transfer.py @@ -29,6 +29,7 @@ import logging import math +import warnings from collections.abc import Sequence import numpy as np @@ -315,14 +316,17 @@ def ntu_og(y_bottom: float, y_top: float, m_slope: float = 0.0) -> float: NTU_OG ≈ ln(Y_bottom / Y_top) - For a rigorous calculation including the equilibrium back-pressure, - pass the Henry's law slope *m_slope*. - Args: y_bottom: CO₂ mole fraction entering the absorber (bottom). y_top: CO₂ mole fraction leaving the absorber (top). - m_slope: Henry's law slope m = y* / x (dimensionless). Set to 0 - (default) to neglect equilibrium back-pressure. + m_slope: **Deprecated and currently ignored.** In v0.1.x this + parameter was documented as enabling an equilibrium-corrected + (Colburn) calculation, but no such code path exists — the + function always returns the dilute-limit ``ln(Y_b / Y_t)``. + A non-zero value now raises ``DeprecationWarning``; the + parameter will be removed in v0.2.0. For a rigorous Colburn + calculation, use a forthcoming dedicated entry point that + takes the absorption factor explicitly (see issue #5). Returns: Dimensionless NTU_OG. @@ -336,6 +340,19 @@ def ntu_og(y_bottom: float, y_top: float, m_slope: float = 0.0) -> float: if y_top >= y_bottom: raise ValueError(f"y_top ({y_top}) must be < y_bottom ({y_bottom}) for absorption.") + if m_slope != 0.0: + warnings.warn( + "ntu_og(m_slope=...) has never been implemented and is silently " + "ignored — the function always returns the dilute-limit " + "ln(Y_bottom / Y_top). The parameter will be removed in meapy " + "v0.2.0; until then, calling with a non-zero m_slope produces " + "the same result as calling with m_slope=0. See " + "https://github.com/defnalk/meapy/issues/5 for the planned " + "Colburn-based replacement API.", + DeprecationWarning, + stacklevel=2, + ) + Y_b = mole_fraction_to_ratio(y_bottom) Y_t = mole_fraction_to_ratio(y_top) return math.log(Y_b / Y_t) diff --git a/tests/unit/test_mass_transfer.py b/tests/unit/test_mass_transfer.py index 6a13024..135b43d 100644 --- a/tests/unit/test_mass_transfer.py +++ b/tests/unit/test_mass_transfer.py @@ -3,6 +3,7 @@ from __future__ import annotations import math +import warnings import numpy as np import pytest @@ -193,6 +194,20 @@ def test_invalid_fractions_raise(self): with pytest.raises(ValueError): ntu_og(0.02, 0.14) # y_top ≥ y_bottom + def test_nonzero_m_slope_warns_and_returns_dilute_limit(self): + """m_slope is deprecated (issue #5): warn but return the same value as m_slope=0.""" + baseline = ntu_og(0.14, 0.02) + with pytest.warns(DeprecationWarning, match="m_slope"): + warned = ntu_og(0.14, 0.02, m_slope=0.85) + assert warned == pytest.approx(baseline, rel=1e-12) + + def test_zero_m_slope_does_not_warn(self): + """The default m_slope=0 path must remain warning-free.""" + with warnings.catch_warnings(): + warnings.simplefilter("error") # any warning becomes an error + ntu_og(0.14, 0.02) + ntu_og(0.14, 0.02, m_slope=0.0) + # --------------------------------------------------------------------------- # hog From 9c5270d99db352baf431c6d4ae334d9f7f2a8a0b Mon Sep 17 00:00:00 2001 From: Defne Nihal Ertugrul Date: Thu, 9 Apr 2026 20:52:39 +0300 Subject: [PATCH 6/6] fix(constants): import Mapping from collections.abc (ruff UP035) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit CI was failing on this branch with: UP035 Import from `collections.abc` instead: `Mapping` --> src/meapy/constants.py:20 Move the Mapping import. Behaviour unchanged — Mapping is the same ABC in either location since Python 3.9. --- src/meapy/constants.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/meapy/constants.py b/src/meapy/constants.py index 454e9a3..80de839 100644 --- a/src/meapy/constants.py +++ b/src/meapy/constants.py @@ -16,8 +16,8 @@ from __future__ import annotations import math +from collections.abc import Mapping from types import MappingProxyType -from typing import Mapping __all__ = [ "MEAProperties",