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/constants.py b/src/meapy/constants.py index bd52245..80de839 100644 --- a/src/meapy/constants.py +++ b/src/meapy/constants.py @@ -15,6 +15,10 @@ from __future__ import annotations +import math +from collections.abc import Mapping +from types import MappingProxyType + __all__ = [ "MEAProperties", "PlantGeometry", @@ -81,8 +85,8 @@ 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² - SAMPLING_HEIGHTS_M: list[float] = [0.5, 1.5, 2.5, 3.5, 4.5, 5.5] # ports 1–6 + COLUMN_CROSS_SECTION_M2: float = math.pi * (0.1 / 2) ** 2 # π r² ≈ 7.854e-3 m² + 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" @@ -134,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, + } + ) 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 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/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 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