Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 67 additions & 0 deletions benchmarks/bench_rolling_mean.py
Original file line number Diff line number Diff line change
@@ -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()
26 changes: 16 additions & 10 deletions src/meapy/constants.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,10 @@

from __future__ import annotations

import math
from collections.abc import Mapping
from types import MappingProxyType

__all__ = [
"MEAProperties",
"PlantGeometry",
Expand Down Expand Up @@ -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"
Expand Down Expand Up @@ -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,
}
)
10 changes: 10 additions & 0 deletions src/meapy/heat_transfer.py
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down
27 changes: 22 additions & 5 deletions src/meapy/mass_transfer.py
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@

import logging
import math
import warnings
from collections.abc import Sequence

import numpy as np
Expand Down Expand Up @@ -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.
Expand All @@ -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)
Expand Down
19 changes: 14 additions & 5 deletions src/meapy/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
15 changes: 15 additions & 0 deletions tests/unit/test_mass_transfer.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
from __future__ import annotations

import math
import warnings

import numpy as np
import pytest
Expand Down Expand Up @@ -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
Expand Down
Loading