Coverage for src/gwtransport/residence_time.py: 95%

1""" 

2Residence Time Calculations for Retarded Compound Transport. 

3 

4This module provides functions to compute residence times for compounds traveling through 

5aquifer systems, accounting for flow variability, pore volume, and retardation due to 

6physical or chemical interactions with the aquifer matrix. Residence time represents the 

7duration a compound spends traveling from infiltration to extraction points, depending on 

8flow rate (higher flow yields shorter residence time), pore volume (larger volume yields 

9longer residence time), and retardation factor (interaction with matrix yields longer 

10residence time). 

11 

12The residence times are resolved per pore volume (:func:`full`), collapsed over a discrete 

13equally-weighted pore-volume distribution (:func:`mean`), or in closed form over a (shifted) gamma 

14aquifer pore-volume distribution with no discretization (:func:`gamma`). 

15 

16Spin-up period 

17-------------- 

18The spin-up **region** is determined entirely by the supplied flow record (``tedges``, which 

19fixes the cumulative throughflow volume ``V`` from ``0`` at the record start to ``V_end`` at the 

20record end) together with the retarded pore volume ``retardation_factor * V_p`` -- it is not a 

21length you set. A residence time for an output time needs the corresponding parcel to stay inside 

22the flow record: 

23 

24* ``direction='extraction_to_infiltration'`` looks **back** to the infiltration event, so the 

25 spin-up sits at the **start** of the output record: the residence time of a pore volume ``V_p`` 

26 needs ``V(t) >= retardation_factor * V_p`` (the extracted water was infiltrated before the record 

27 began otherwise). 

28* ``direction='infiltration_to_extraction'`` looks **forward** to the extraction event, so the 

29 spin-up sits at the **end** of the output record: it needs ``V_end - V(t) >= retardation_factor * 

30 V_p`` (the infiltrated water is extracted after the record ends otherwise). 

31 

32The spin-up therefore lengthens with both the pore volume and the retardation factor, and is 

33longest for the largest pore volumes of a distribution. 

34 

35What happens in that region is governed by a ``spinup`` policy, following the package convention 

36(see :mod:`gwtransport.advection`); :func:`full`, :func:`mean` and 

37:func:`gamma` all share the contract ``spinup={'constant'} | None | float in 

38[0, 1]`` and the default is ``"constant"`` everywhere: 

39 

40* ``"constant"`` (default) **warm-starts** by extrapolating the boundary flow (flow held constant 

41 at its first/last value), so no in-record output is ``NaN``. 

42* ``None`` is strict (no extrapolation), marking a pore volume ``NaN`` for any output bin its parcel 

43 leaves the record within. Where the pore-volume axis is collapsed -- :func:`mean` over a 

44 discrete set, :func:`gamma` over the continuum -- the bin mean then **renormalizes** 

45 over the covered streamtubes / sub-mass, emitted wherever any coverage remains. 

46* a ``float`` covered-fraction threshold is the strict mode with a minimum coverage gate: the 

47 renormalized mean is emitted only where the covered streamtube fraction / sub-mass fraction is at 

48 least ``spinup`` (``0.0`` matches ``None``; larger values demand more coverage). For the 

49 per-pore-volume :func:`full` there is no axis to collapse, so the ``float`` behaves 

50 exactly like ``None``. 

51 

52Output bins lying wholly outside ``tedges`` are ``NaN`` under every policy. 

53 

54The :func:`fraction_explained_full` / :func:`fraction_explained_mean` / 

55:func:`fraction_explained_gamma` diagnostics report, per output bin, the advective fraction of the 

56pore-volume distribution that is out of spin-up (``1.0`` = advectively fully informed, ``0.0`` = 

57entirely in spin-up) and are the way to locate the spin-up region when the means warm-start over it. 

58They are **purely advective** -- molecular diffusion and microdispersion spread each bin over a 

59range of infiltration times that is not captured there, so no bin is fully informed once dispersion 

60is present (for that dispersive informed fraction use the captured kernel mass of the diffusion 

61coefficient matrix). 

62 

63Available functions: 

64 

65- :func:`full` - Flow-weighted mean residence time [days] over each output bin, resolved per pore volume: 

66 the full ``(n_pore_volumes, n_output_bins)`` array, without collapsing the pore-volume axis. The bin 

67 average is uniform in cumulative throughflow volume, and ``direction`` selects whether the time is the 

68 look-back to infiltration (``extraction_to_infiltration``) or the look-forward to extraction 

69 (``infiltration_to_extraction``). 

70 

71- :func:`mean` - Mean residence time [days] per output bin for a discrete aquifer pore-volume 

72 distribution: the :func:`full` array collapsed by averaging over the equally-weighted streamtubes that 

73 are valid in each bin, shape ``(n_output_bins,)``. 

74 

75- :func:`gamma` - Mean residence time [days] per output bin for a continuous (shifted) gamma aquifer 

76 pore-volume distribution, parameterized by either ``(mean, std, loc)`` or ``(alpha, beta, loc)``. The 

77 expectation is taken in closed form from regularized incomplete-gamma partial moments, so there is no 

78 pore-volume discretization and no accuracy/cost knob; shape ``(n_output_bins,)``. 

79 

80- :func:`fraction_explained_full` - Advective coverage per pore volume: the flow-weighted fraction of each 

81 output bin, in ``[0, 1]``, whose retarded parcel lies inside the supplied flow record. Returned as the 

82 full ``(n_pore_volumes, n_output_bins)`` array, mirroring :func:`full`. 

83 

84- :func:`fraction_explained_mean` - Equally-weighted mean of :func:`fraction_explained_full` over the 

85 discrete streamtubes, shape ``(n_output_bins,)``. 

86 

87- :func:`fraction_explained_gamma` - Closed-form expectation of the advective in-record indicator over a 

88 (shifted) gamma pore-volume distribution, shape ``(n_output_bins,)`` -- the continuum analogue of 

89 :func:`fraction_explained_mean`, evaluated from the antiderivative of the shifted-gamma CDF. 

90 

91- :func:`freundlich_retardation` - Concentration-dependent retardation factors 

92 ``R = 1 + (rho_b / theta) * k_f * (1 / n) * C ** (1 / n - 1)`` from the Freundlich isotherm 

93 ``s = k_f * C ** (1 / n)``, one per concentration entry, for use as the ``retardation_factor`` input of 

94 the transport functions. 

95 

96This file is part of gwtransport which is released under AGPL-3.0 license. 

97See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details. 

98""" 

99 

100import numpy as np 

101import numpy.typing as npt 

102import pandas as pd 

103from scipy.stats import gamma as gamma_dist 

104 

105from gwtransport._time import tedges_to_days 

106from gwtransport.gamma import parse_parameters 

107from gwtransport.utils import cumulative_flow_volume, linear_interpolate 

108 

109# Relative slack on the covered-fraction spin-up gate. The covered sub-mass ``den`` equals its 

110# fully-covered reference only up to summation-reassociation ulps, so an exact ``den >= threshold * 

111# reference`` comparison spuriously rejects fully-covered bins at the strictest threshold (1.0). The 

112# slack is far above that float noise (band widths reach ~1e4 pieces, ~1e-12 relative) yet negligible 

113# as a physical coverage tolerance. 

114_SPINUP_GATE_RELTOL = 1e-9 

115 

116 

117def _resolve_spinup(spinup: str | float | None) -> tuple[bool, float]: 

118 """Normalize the residence-time ``spinup`` policy to ``(extrapolate, threshold)``. 

119 

120 ``'constant'`` -> ``(True, 0.0)``, ``None`` -> ``(False, 0.0)``, and a ``float`` in ``[0, 1]`` -> 

121 ``(False, float)``. See the module docstring (``Spin-up period``) for the public contract. 

122 

123 Returns 

124 ------- 

125 extrapolate : bool 

126 Whether to warm-start by extrapolating the boundary flow. 

127 threshold : float 

128 Covered-fraction gate applied where the pore-volume axis is collapsed (ignored by the 

129 per-pore-volume :func:`full`, which is strict per bin). 

130 

131 Raises 

132 ------ 

133 ValueError 

134 If ``spinup`` is not ``'constant'``, ``None``, or a float in ``[0, 1]``. 

135 """ 

136 if spinup == "constant": 

137 return True, 0.0 

138 if spinup is None: 

139 return False, 0.0 

140 if isinstance(spinup, int | float) and not isinstance(spinup, bool) and 0.0 <= spinup <= 1.0: 

141 return False, float(spinup) 

142 msg = "spinup should be 'constant', None, or a float in [0, 1]" 

143 raise ValueError(msg) 

144 

145 

146def _phi_setup( 

147 flow: npt.NDArray[np.floating], 

148 flow_cum: npt.NDArray[np.floating], 

149 tedges_days: npt.NDArray[np.floating], 

150 *, 

151 extrapolate: bool, 

152 pad: float, 

153) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating], npt.NDArray[np.floating], npt.NDArray[np.floating]]: 

154 """Build the antiderivative ``phi(x) = int_0^x T(w) dw`` of the cumulative-volume -> time map ``T``. 

155 

156 ``T`` is piecewise-linear (knots at the cumulative-volume edges), so ``phi`` is piecewise-quadratic 

157 with the same knots. With ``extrapolate`` the map is extended one anchor past each end of the record 

158 at the boundary flow rate (padded by ``pad``) for the ``"constant"`` warm-start; otherwise it is the 

159 raw record. Shared by :func:`full` and :func:`gamma`. 

160 

161 Returns 

162 ------- 

163 phi_v : ndarray 

164 Cumulative-volume knots of the map (extended by ``pad`` at each end when ``extrapolate``). 

165 phi_t : ndarray 

166 Time knots aligned with ``phi_v``. 

167 phi_knot : ndarray 

168 ``phi`` evaluated at each volume knot. 

169 phi_rate : ndarray 

170 Per-segment ``dV/dt`` (flow) rate between consecutive knots. 

171 """ 

172 # pad == 0 (no spin-up reach, e.g. retardation_factor or all pore volumes 0) carries no 

173 # extrapolation and would only add zero-width boundary segments (0/0 -> NaN phi_rate), so skip it. 

174 if extrapolate and pad > 0.0 and np.any(flow > 0.0): 

175 # The boundary extrapolation slope is 1/Q, anchored on the nearest strictly-positive flow so 

176 # a zero-flow boundary bin (whose flow_cum step is only the strictly-monotone ulp bump) does 

177 # not give a 1/0 extrapolation. 

178 positive_flow = flow[flow > 0.0] 

179 phi_v = np.concatenate([[flow_cum[0] - pad], flow_cum, [flow_cum[-1] + pad]]) 

180 phi_t = np.concatenate([ 

181 [tedges_days[0] - pad * (1.0 / positive_flow[0])], 

182 tedges_days, 

183 [tedges_days[-1] + pad * (1.0 / positive_flow[-1])], 

184 ]) 

185 else: 

186 phi_v, phi_t = flow_cum, tedges_days 

187 phi_dv = phi_v[1:] - phi_v[:-1] 

188 phi_rate = phi_dv / (phi_t[1:] - phi_t[:-1]) 

189 phi_knot = np.concatenate([[0.0], np.cumsum(phi_t[:-1] * phi_dv + phi_dv**2 / (2 * phi_rate))]) 

190 return phi_v, phi_t, phi_knot, phi_rate 

191 

192 

193def _eval_phi( 

194 x: npt.NDArray[np.floating], 

195 phi_v: npt.NDArray[np.floating], 

196 phi_t: npt.NDArray[np.floating], 

197 phi_knot: npt.NDArray[np.floating], 

198 phi_rate: npt.NDArray[np.floating], 

199 *, 

200 strict_nan: bool = False, 

201) -> npt.NDArray[np.floating]: 

202 """Evaluate the piecewise-quadratic antiderivative ``phi`` from :func:`_phi_setup` at ``x``. 

203 

204 ``x`` is clipped to the map range ``[phi_v[0], phi_v[-1]]`` (the warm-start extrapolation lives in 

205 that range when padded). With ``strict_nan`` any ``x`` outside the range returns ``NaN`` instead -- 

206 used by :func:`full` so an output bin whose parcel leaves the record is ``NaN``. 

207 

208 Returns 

209 ------- 

210 ndarray 

211 ``phi`` evaluated at ``x``, with the same shape as ``x``. 

212 """ 

213 x = np.asarray(x, dtype=float) 

214 xc = np.clip(x, phi_v[0], phi_v[-1]) 

215 j = np.clip(np.searchsorted(phi_v, xc, side="right") - 1, 0, len(phi_rate) - 1) 

216 dv = xc - phi_v[j] 

217 out = phi_knot[j] + phi_t[j] * dv + dv * dv / (2 * phi_rate[j]) 

218 if strict_nan: 

219 out = np.where((x < phi_v[0]) | (x > phi_v[-1]), np.nan, out) 

220 return out 

221 

222 

223def full( 

224 *, 

225 flow: npt.ArrayLike, 

226 tedges: pd.DatetimeIndex | np.ndarray, 

227 cout_tedges: pd.DatetimeIndex | np.ndarray, 

228 aquifer_pore_volumes: npt.ArrayLike, 

229 direction: str = "extraction_to_infiltration", 

230 retardation_factor: float = 1.0, 

231 spinup: str | float | None = "constant", 

232) -> npt.NDArray[np.floating]: 

233 r""" 

234 Compute the mean residence time over output bins, per pore volume. 

235 

236 The flow-weighted mean residence time is computed over each output interval 

237 ``[cout_tedges[i], cout_tedges[i + 1])`` and returned as the full 

238 ``(n_pore_volumes, n_output_bins)`` array -- one row per entry in 

239 ``aquifer_pore_volumes``, without collapsing the pore-volume axis. The average is uniform in 

240 cumulative throughflow volume, matching the package's bin-edge convention. 

241 

242 Parameters 

243 ---------- 

244 flow : array-like 

245 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one. 

246 tedges : pandas.DatetimeIndex 

247 Time edges for the flow data, as datetime64 objects, defining the flow intervals. 

248 cout_tedges : pandas.DatetimeIndex 

249 Output time edges as datetime64 objects; ``n + 1`` edges define ``n`` output bins. 

250 aquifer_pore_volumes : float or array-like 

251 Pore volume(s) of the aquifer [m³]. A single value or an array of pore volumes 

252 representing different flow paths. 

253 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional 

254 Direction of the flow calculation: 

255 

256 * 'extraction_to_infiltration': 

257 Extraction to infiltration modeling - how many days ago was the extracted water infiltrated. 

258 * 'infiltration_to_extraction': 

259 Infiltration to extraction modeling - how many days until the infiltrated water is extracted. 

260 

261 Default is 'extraction_to_infiltration'. 

262 retardation_factor : float, optional 

263 Retardation factor of the compound in the aquifer [dimensionless]. A value greater 

264 than 1.0 indicates the compound moves slower than water. Default is 1.0. 

265 spinup : {'constant'}, None, or float in [0, 1], optional 

266 How to treat the spin-up zone, where a pore volume's retarded look-back/forward parcel 

267 leaves the flow record. Matches the package convention (see :mod:`gwtransport.advection`). 

268 

269 * ``'constant'`` (default): warm-start -- extrapolate the cumulative-volume-to-time map past 

270 the record at the boundary flow rates (flow held constant at its first/last value), so 

271 the residence time stays finite. No left-edge (extraction) or right-edge (infiltration) 

272 spin-up ``NaN``. 

273 * ``None`` or a ``float`` in ``[0, 1]``: strict -- a pore volume whose parcel leaves the 

274 record at any point within an output bin is ``NaN`` for that bin (all-or-nothing per bin), 

275 with no extrapolation. This function returns the full per-pore-volume array, so there is no 

276 pore-volume axis to collapse; the ``float`` covered-fraction threshold therefore behaves 

277 identically to ``None`` here and only takes effect once the axis is collapsed in 

278 :func:`mean` / :func:`gamma`. 

279 

280 Output bins lying wholly outside ``tedges`` are ``NaN`` under either policy. 

281 

282 Returns 

283 ------- 

284 numpy.ndarray 

285 Mean residence time [days], shape ``(n_pore_volumes, n_output_bins)``. The first 

286 dimension corresponds to the pore volumes and the second to the ``cout_tedges`` bins. 

287 Negative or ``NaN`` ``flow`` makes the cumulative-volume map non-monotone or undefined; the 

288 whole array is returned as ``NaN`` (the function refuses rather than raising). 

289 

290 Raises 

291 ------ 

292 ValueError 

293 If ``tedges`` does not have exactly one more element than ``flow``. If 

294 ``direction`` is not ``'extraction_to_infiltration'`` or 

295 ``'infiltration_to_extraction'``. If ``spinup`` is not ``'constant'``, ``None``, or a float 

296 in ``[0, 1]``. 

297 

298 See Also 

299 -------- 

300 fraction_explained_full : Advective fraction of each output bin explained, per pore volume 

301 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction 

302 :ref:`concept-transport-equation` : Flow-weighted averaging convention 

303 

304 Notes 

305 ----- 

306 With the default ``spinup='constant'`` the spin-up zone is warm-started by extrapolating the 

307 boundary flow, so no in-record bin is ``NaN``; use :func:`fraction_explained_mean` (or 

308 ``spinup=None``) to locate the spin-up region. See the module docstring (``Spin-up period``) 

309 for the full rule. 

310 

311 The single-streamtube residence time :math:`\tau(V) = \mathrm{sign}\,[T(V + \mathrm{sign}\,R V_p) 

312 - T(V)]` is piecewise-linear in cumulative throughflow volume :math:`V` (:math:`T` is the 

313 volume :math:`\to` time map, :math:`\mathrm{sign} = -1` for ``extraction_to_infiltration`` and 

314 :math:`+1` for ``infiltration_to_extraction``). Its flow-weighted bin average is therefore a 

315 closed-form difference of the antiderivative :math:`\Phi(x) = \int_0^x T(w)\,dw` (piecewise- 

316 quadratic), evaluated at four points per pore volume and output bin: 

317 

318 .. math:: 

319 

320 \bar\tau 

321 = \frac{1}{\Delta V}\int_{V_\mathrm{lo}}^{V_\mathrm{hi}} \tau(V)\,dV 

322 = \frac{\mathrm{sign}}{\Delta V}\bigl[ 

323 \Phi(V_\mathrm{hi} + \mathrm{sign}\,R V_p) - \Phi(V_\mathrm{lo} + \mathrm{sign}\,R V_p) 

324 - \Phi(V_\mathrm{hi}) + \Phi(V_\mathrm{lo})\bigr], 

325 

326 where :math:`V` is cumulative throughflow volume (:math:`dV = Q\,dt`). This avoids materialising a 

327 per-streamtube integration grid, so memory and time scale as the output size 

328 :math:`O(n_\mathrm{pore\ volumes}\cdot n_\mathrm{bins})`. A zero-throughflow output bin 

329 (:math:`\Delta V \to 0`) has a fixed volume while output time advances, so it degenerates to the 

330 pointwise residence time at the bin's time midpoint. 

331 

332 Examples 

333 -------- 

334 >>> import pandas as pd 

335 >>> import numpy as np 

336 >>> from gwtransport.residence_time import full 

337 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-01-10", freq="D") 

338 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # Constant flow of 100 m³/day 

339 >>> mean_times = full( 

340 ... flow=flow_values, 

341 ... tedges=flow_dates, 

342 ... cout_tedges=flow_dates, 

343 ... aquifer_pore_volumes=200.0, 

344 ... direction="extraction_to_infiltration", 

345 ... ) 

346 >>> # 200 m³ / 100 m³/day = 2 days residence time; the default constant warm-start 

347 >>> # extrapolates the boundary flow, so the left-edge spin-up bins are also 2 days 

348 >>> print(mean_times) # doctest: +NORMALIZE_WHITESPACE 

349 [[2. 2. 2. 2. 2. 2. 2. 2. 2.]] 

350 """ 

351 if direction not in {"extraction_to_infiltration", "infiltration_to_extraction"}: 

352 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'" 

353 raise ValueError(msg) 

354 extrapolate, _ = _resolve_spinup(spinup) 

355 

356 aquifer_pore_volumes = np.atleast_1d(aquifer_pore_volumes) 

357 tedges = pd.DatetimeIndex(tedges) 

358 cout_tedges = pd.DatetimeIndex(cout_tedges) 

359 flow = np.asarray(flow, dtype=float) 

360 n_pv = len(aquifer_pore_volumes) 

361 n_out = len(cout_tedges) - 1 

362 

363 if len(tedges) != len(flow) + 1: 

364 msg = "tedges must have one more element than flow" 

365 raise ValueError(msg) 

366 if np.any(flow < 0) or np.any(np.isnan(flow)): 

367 return np.full((n_pv, n_out), np.nan) 

368 

369 tedges_days = tedges_to_days(tedges) 

370 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0]) 

371 # Plateaus in flow_cum from Q = 0 bins make the V -> t inversion multi-valued; bump duplicates by 

372 # the smallest representable amount so the inverse map is single-valued. 

373 flow_cum = cumulative_flow_volume(flow, np.diff(tedges_days), strictly_monotone=True) 

374 

375 # Sign convention: sign = -1 for extraction_to_infiltration, +1 for infiltration_to_extraction; 

376 # the look-back/forward parcel sits at volume V + shift and tau(V) = sign * (T(V + shift) - T(V)) is 

377 # piecewise-linear in V (T is the volume -> time map). Its flow-weighted bin average is a closed- 

378 # form difference of an antiderivative phi with phi' = T (any additive constant cancels in the 

379 # difference below; over the extrapolated map T is the extended map), so no per-streamtube 

380 # integration grid is built (memory/time O(n_pore_volumes * n_bins), not O(n_pore_volumes^2 * n_flow)). 

381 sign = -1.0 if direction == "extraction_to_infiltration" else 1.0 

382 shift = sign * retardation_factor * aquifer_pore_volumes # (n_pv,) 

383 

384 # phi over the cumulative-volume -> time map. With spinup="constant" the map is extrapolated past 

385 # the record at the boundary flow (padded by the largest reach R * max(V_p)) so phi warm-starts the 

386 # spin-up; otherwise phi is NaN outside the record so a bin whose parcel leaves it becomes NaN. 

387 pad = retardation_factor * float(aquifer_pore_volumes.max()) if aquifer_pore_volumes.size else 0.0 

388 phi_v, phi_t, phi_knot, phi_rate = _phi_setup(flow, flow_cum, tedges_days, extrapolate=extrapolate, pad=pad) 

389 # The map is only actually extended when there is a positive boundary flow to extrapolate; with 

390 # all-zero flow (or spinup=None) it stays the raw record, so out-of-record look-backs are NaN. 

391 strict = not (extrapolate and bool(np.any(flow > 0.0))) 

392 

393 vol_out = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days, left=np.nan, right=np.nan) 

394 v_lo = vol_out[:-1] 

395 v_hi = vol_out[1:] 

396 dvol = v_hi - v_lo 

397 bins_within = np.isfinite(v_lo) & np.isfinite(v_hi) 

398 

399 phi_base = _eval_phi(vol_out, phi_v, phi_t, phi_knot, phi_rate, strict_nan=strict) # (n_out + 1,) 

400 phi_shift = _eval_phi(vol_out[None, :] + shift[:, None], phi_v, phi_t, phi_knot, phi_rate, strict_nan=strict) 

401 

402 tol = 1e6 * np.spacing(float(np.max(np.abs(flow_cum)))) if flow_cum.size else 0.0 

403 with np.errstate(divide="ignore", invalid="ignore"): 

404 result = ( 

405 sign * ((phi_shift[:, 1:] - phi_shift[:, :-1]) - (phi_base[1:] - phi_base[:-1])[None, :]) / dvol[None, :] 

406 ) 

407 # Zero-throughflow output bin (dvol -> 0): the volume is fixed while output time advances, so the 

408 # flow-weighted average degenerates to the pointwise tau at the bin time midpoint, sign*(T(V_lo + 

409 # shift) - t_mid). A direct ratio there would catastrophically cancel. Only in-record zero-flow 

410 # bins reach it, so skip the interp entirely when none are present. 

411 degenerate = bins_within & (dvol <= tol) 

412 if np.any(degenerate): 

413 t_mid = 0.5 * (cout_tedges_days[:-1] + cout_tedges_days[1:]) 

414 nan_outside = np.nan if strict else None 

415 with np.errstate(divide="ignore", invalid="ignore"): 

416 t_lookback = linear_interpolate( 

417 x_ref=phi_v, y_ref=phi_t, x_query=v_lo[None, :] + shift[:, None], left=nan_outside, right=nan_outside 

418 ) 

419 point = sign * (t_lookback - t_mid[None, :]) 

420 result = np.where(dvol[None, :] > tol, result, point) 

421 return np.where(bins_within[None, :], result, np.nan) 

422 

423 

424def mean( 

425 *, 

426 flow: npt.ArrayLike, 

427 tedges: pd.DatetimeIndex | np.ndarray, 

428 cout_tedges: pd.DatetimeIndex | np.ndarray, 

429 aquifer_pore_volumes: npt.ArrayLike, 

430 direction: str = "extraction_to_infiltration", 

431 retardation_factor: float = 1.0, 

432 spinup: str | float | None = "constant", 

433) -> npt.NDArray[np.floating]: 

434 r""" 

435 Compute the mean residence time over output bins for a discrete APVD. 

436 

437 The mean is taken over a **discrete** set of equally-weighted aquifer pore volumes -- one 

438 streamtube per entry in ``aquifer_pore_volumes``. Each streamtube's flow-weighted bin average 

439 is computed with :func:`full` and the pore-volume axis is then collapsed to a 

440 single per-output-bin series by averaging over the streamtubes that are valid in each bin. For 

441 a continuous (shifted) gamma pore-volume distribution evaluated in closed form, use 

442 :func:`gamma`. 

443 

444 The mean is over the valid streamtubes, 

445 

446 .. math:: 

447 

448 \bar\tau_b = \frac{1}{|V_b|}\sum_{i \in V_b} \tau_{i,b}, 

449 \qquad V_b = \{\, i : \tau_{i,b}\ \mathrm{finite} \,\}. 

450 

451 With the default ``spinup='constant'`` every streamtube is finite within the flow record 

452 (the boundary flow is extrapolated), so this is simply the mean over all pore volumes; with 

453 ``spinup=None`` it renormalizes over the streamtubes that have broken through. 

454 

455 Parameters 

456 ---------- 

457 flow : array-like 

458 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one. 

459 tedges : pandas.DatetimeIndex 

460 Time edges for the flow data, as datetime64 objects, defining the flow intervals. 

461 cout_tedges : pandas.DatetimeIndex 

462 Output time edges as datetime64 objects; ``n + 1`` edges define ``n`` output bins. 

463 aquifer_pore_volumes : array-like 

464 Discrete pore volumes [m³], one per (equally-weighted) streamtube. A single value 

465 collapses to the per-streamtube mean of :func:`full`. 

466 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional 

467 Direction of the flow calculation: 

468 * 'extraction_to_infiltration': how many days ago was the extracted water infiltrated 

469 * 'infiltration_to_extraction': how many days until the infiltrated water is extracted 

470 Default is 'extraction_to_infiltration'. 

471 retardation_factor : float, optional 

472 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0. 

473 spinup : {'constant'}, None, or float in [0, 1], optional 

474 Spin-up policy, sharing the contract of :func:`gamma`. ``'constant'`` 

475 (default) warm-starts by extrapolating the boundary flow so no in-record bin is ``NaN``; 

476 ``None`` leaves spin-up streamtubes ``NaN`` and the mean renormalizes over those that have 

477 broken through (emitted wherever at least one streamtube is valid). A ``float`` in 

478 ``[0, 1]`` is the covered-fraction threshold: the renormalized mean is emitted only where 

479 the fraction of valid streamtubes is at least ``spinup`` (``0.0`` matches ``None``; ``1.0`` 

480 demands every streamtube; larger values demand more streamtubes to have broken through). Use 

481 :func:`fraction_explained_mean` to 

482 locate the spin-up region. 

483 

484 Returns 

485 ------- 

486 numpy.ndarray 

487 Mean residence time [days], shape ``(n_output_bins,)``. Output bins with no valid 

488 streamtube (outside the flow record, or -- with ``spinup=None`` -- fully in the spin-up 

489 zone) are NaN; with a ``float`` ``spinup`` so are bins whose valid-streamtube fraction is 

490 below the threshold. Negative or ``NaN`` ``flow`` makes the cumulative-volume map non-monotone 

491 or undefined; the whole series is returned as ``NaN`` (the function refuses rather than raising). 

492 

493 See Also 

494 -------- 

495 gamma : Exact closed-form mean for a continuous (shifted) gamma APVD 

496 full : Per-pore-volume mean residence time over output bins 

497 fraction_explained_mean : Advective fraction of each output bin explained by the record 

498 gwtransport.gamma.bins : Discretize a gamma APVD into pore-volume bins 

499 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction 

500 

501 Notes 

502 ----- 

503 With ``spinup=None`` the spin-up is **all-or-nothing per streamtube**: a streamtube whose 

504 look-back/forward parcel leaves the flow record part-way through an output bin has a ``NaN`` bin 

505 average (inherited from :func:`full`) and is dropped from that bin's mean 

506 entirely, rather than contributing its partially-covered share; the bin is ``NaN`` only once 

507 every streamtube is in spin-up. In that mode the discrete mean differs from 

508 :func:`gamma`, which renormalizes over the covered sub-mass exactly. See the 

509 module docstring (``Spin-up period``) for the full rule. 

510 

511 Examples 

512 -------- 

513 >>> import pandas as pd 

514 >>> import numpy as np 

515 >>> from gwtransport.residence_time import mean 

516 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-02-10", freq="D") 

517 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # 100 m³/day 

518 >>> tau_bar = mean( 

519 ... flow=flow_values, 

520 ... tedges=flow_dates, 

521 ... cout_tedges=flow_dates, 

522 ... aquifer_pore_volumes=[400.0, 600.0], # two equally-weighted streamtubes 

523 ... ) 

524 >>> # Deep in the record: mean pore volume 500 / 100 m³/day = 5 days 

525 >>> float(np.round(tau_bar[-1], 6)) 

526 5.0 

527 """ 

528 _, threshold = _resolve_spinup(spinup) 

529 rt = full( 

530 flow=flow, 

531 tedges=tedges, 

532 cout_tedges=cout_tedges, 

533 aquifer_pore_volumes=aquifer_pore_volumes, 

534 direction=direction, 

535 retardation_factor=retardation_factor, 

536 spinup=spinup, 

537 ) 

538 

539 # Mean over the streamtubes that are valid (non-NaN) in each output bin; bins with no valid 

540 # streamtube reduce to 0/0 and are NaN. With spinup='constant' every in-record streamtube is 

541 # finite, so this is the plain mean; otherwise it renormalizes over the broken-through set and 

542 # a float covered-fraction threshold further NaNs bins where too few streamtubes have arrived. 

543 n_streamtubes = rt.shape[0] 

544 valid_count = np.isfinite(rt).sum(axis=0) 

545 with np.errstate(invalid="ignore"): 

546 bin_mean = np.nansum(rt, axis=0) / valid_count 

547 if threshold > 0.0: 

548 bin_mean = np.where(valid_count >= threshold * n_streamtubes, bin_mean, np.nan) 

549 return bin_mean 

550 

551 

552def gamma( 

553 *, 

554 flow: npt.ArrayLike, 

555 tedges: pd.DatetimeIndex | np.ndarray, 

556 cout_tedges: pd.DatetimeIndex | np.ndarray, 

557 mean: float | None = None, 

558 std: float | None = None, 

559 loc: float = 0.0, 

560 alpha: float | None = None, 

561 beta: float | None = None, 

562 direction: str = "extraction_to_infiltration", 

563 retardation_factor: float = 1.0, 

564 spinup: str | float | None = "constant", 

565 _max_tile_elements: int = 1_000_000, 

566) -> npt.NDArray[np.floating]: 

567 r""" 

568 Compute the mean residence time over output bins for a (shifted) gamma APVD. 

569 

570 The expectation over a (shifted) gamma aquifer pore-volume distribution (APVD), 

571 parameterized by either ``(mean, std, loc)`` or ``(alpha, beta, loc)``, is taken in closed 

572 form -- no pore-volume binning and no ``n_bins`` accuracy/cost knob. The bin mean is 

573 flow-weighted (uniform in cumulative volume), matching the bin-edge convention of the 

574 package, and a single per-output-bin series is returned. 

575 

576 The single-streamtube residence time is piecewise-linear in the pore volume :math:`V_p`, so 

577 its per-bin time integral :math:`G_b(V_p) = \int_{\mathrm{bin}} \tau\,dV` is piecewise- 

578 quadratic in :math:`V_p` and the covered length :math:`L_b(V_p)` piecewise-linear. The bin 

579 mean is the ratio of two closed-form integrals against the gamma density -- its zeroth, 

580 first and second partial moments (regularized incomplete gamma) -- formed once after 

581 integrating. The ``spinup`` policy sets what happens where part of the APVD lacks flow 

582 history: ``'constant'`` (default) extrapolates the boundary flow over the full distribution 

583 (the package default warm-start), while a ``float`` threshold renormalizes the mean over the 

584 covered sub-mass (``0.0`` reproduces the exact covered-sub-mass conditional mean). 

585 

586 Parameters 

587 ---------- 

588 flow : array-like 

589 Flow rate of water in the aquifer [m³/day]. Length matches ``tedges`` minus one. 

590 tedges : pandas.DatetimeIndex 

591 Time edges for the flow data, as datetime64 objects, defining the flow intervals. 

592 cout_tedges : pandas.DatetimeIndex 

593 Output time edges as datetime64 objects; ``n + 1`` edges define ``n`` output bins. 

594 mean : float, optional 

595 Mean of the gamma APVD [m³]. Must be strictly greater than ``loc``. Provide either 

596 ``(mean, std)`` or ``(alpha, beta)``. 

597 std : float, optional 

598 Standard deviation of the gamma APVD [m³]. Must be positive. 

599 loc : float, optional 

600 Location (lower bound of support) of the gamma APVD [m³]; a guaranteed minimum pore 

601 volume. Must satisfy ``0 <= loc < mean``. Default is 0.0. 

602 alpha : float, optional 

603 Shape parameter of the gamma APVD (must be > 0). 

604 beta : float, optional 

605 Scale parameter of the gamma APVD (must be > 0). 

606 direction : {'extraction_to_infiltration', 'infiltration_to_extraction'}, optional 

607 Direction of the flow calculation: 

608 * 'extraction_to_infiltration': how many days ago was the extracted water infiltrated 

609 * 'infiltration_to_extraction': how many days until the infiltrated water is extracted 

610 Default is 'extraction_to_infiltration'. 

611 retardation_factor : float, optional 

612 Retardation factor of the compound in the aquifer [dimensionless]. Default is 1.0. 

613 spinup : {'constant'}, None, or float in [0, 1], optional 

614 How to treat the spin-up zone, where part of the gamma APVD lacks flow history. Matches 

615 the package convention (see :mod:`gwtransport.advection`). 

616 

617 * ``'constant'`` (default): warm-start -- extrapolate the cumulative-volume-to-time map past 

618 the record at the boundary flow rates (flow held constant at its first/last value) and 

619 integrate the full distribution, so no in-record bin is ``NaN``. 

620 * ``None`` or a ``float`` in ``[0, 1]``: renormalize the mean over the covered sub-mass, 

621 emitting a bin only where the covered fraction of the distribution is at least the 

622 threshold. ``None`` and ``0.0`` both give the exact covered-sub-mass conditional mean 

623 (emit whenever any sub-mass is covered); larger values demand a larger covered fraction, 

624 and ``1.0`` requires the full distribution to be covered. 

625 

626 Output bins lying wholly outside ``tedges`` are ``NaN`` under either policy. 

627 

628 Returns 

629 ------- 

630 numpy.ndarray 

631 APVD-mean residence time [days], shape ``(n_output_bins,)``. Output bins outside the flow 

632 record are NaN; with a ``float`` ``spinup`` so are bins whose covered fraction is below the 

633 threshold. Negative or ``NaN`` ``flow`` makes the cumulative-volume map non-monotone or 

634 undefined; the whole series is returned as ``NaN`` (the function refuses rather than raising). 

635 

636 Raises 

637 ------ 

638 ValueError 

639 If ``tedges`` does not have exactly one more element than ``flow``. If ``direction`` 

640 is not ``'extraction_to_infiltration'`` or ``'infiltration_to_extraction'``. If ``spinup`` 

641 is not ``'constant'``, ``None``, or a float in ``[0, 1]``. Gamma parameter validation is 

642 delegated to :func:`gwtransport.gamma.parse_parameters`. 

643 

644 See Also 

645 -------- 

646 mean : Equally-weighted mean for a discrete set of pore volumes 

647 full : Per-pore-volume mean residence time over output bins 

648 fraction_explained_mean : Advective fraction of each output bin explained by the record 

649 gwtransport.gamma.bins : Discretize a gamma APVD into pore-volume bins 

650 :ref:`concept-residence-time` : Time in aquifer between infiltration and extraction 

651 :ref:`concept-gamma-distribution` : Two-parameter pore volume model 

652 

653 Notes 

654 ----- 

655 With the default ``spinup='constant'`` the spin-up is warm-started exactly as in 

656 :func:`mean` (constant-boundary-flow extrapolation) -- the same warm-start policy applies across 

657 the whole distribution (``mean`` discretizes it, ``gamma`` integrates it in closed form). With 

658 ``spinup=0.0`` the spin-up is instead handled by exact covered-sub-mass renormalization: each 

659 output bin integrates over only the pore-volume sub-range with sufficient flow history. See the 

660 module docstring (``Spin-up period``) for the full rule. 

661 

662 The closed form uses regularized incomplete-gamma CDFs. For an extremely narrow APVD 

663 (``alpha = (mean / std) ** 2`` above ~1e7, i.e. ``std / mean`` below ~3e-4) SciPy's incomplete 

664 gamma loses precision in its far tail and the result degrades to ~1e-5 relative; such a 

665 distribution is effectively a single pore volume, so :func:`mean` with one bin is the exact 

666 alternative there. 

667 

668 Examples 

669 -------- 

670 >>> import pandas as pd 

671 >>> import numpy as np 

672 >>> from gwtransport.residence_time import gamma 

673 >>> flow_dates = pd.date_range(start="2023-01-01", end="2023-02-10", freq="D") 

674 >>> flow_values = np.full(len(flow_dates) - 1, 100.0) # 100 m³/day 

675 >>> tau_bar = gamma( 

676 ... flow=flow_values, 

677 ... tedges=flow_dates, 

678 ... cout_tedges=flow_dates, 

679 ... mean=500.0, 

680 ... std=100.0, 

681 ... direction="extraction_to_infiltration", 

682 ... ) 

683 >>> # Deep in the record the mean residence time approaches mean / flow = 5 days 

684 >>> float(np.round(tau_bar[-1], 6)) 

685 5.0 

686 """ 

687 if direction not in {"extraction_to_infiltration", "infiltration_to_extraction"}: 

688 msg = "direction should be 'extraction_to_infiltration' or 'infiltration_to_extraction'" 

689 raise ValueError(msg) 

690 extrapolate, spinup_threshold = _resolve_spinup(spinup) 

691 

692 alpha, beta, loc = parse_parameters(mean=mean, std=std, loc=loc, alpha=alpha, beta=beta) 

693 

694 flow = np.asarray(flow, dtype=float) 

695 tedges = pd.DatetimeIndex(tedges) 

696 cout_tedges = pd.DatetimeIndex(cout_tedges) 

697 n_out = len(cout_tedges) - 1 

698 

699 if len(tedges) != len(flow) + 1: 

700 msg = "tedges must have one more element than flow" 

701 raise ValueError(msg) 

702 if np.any(flow < 0) or np.any(np.isnan(flow)): 

703 return np.full(n_out, np.nan) 

704 

705 sign = -1.0 if direction == "extraction_to_infiltration" else 1.0 

706 r = retardation_factor 

707 

708 tedges_days = tedges_to_days(tedges) 

709 flow_cum = cumulative_flow_volume(flow, np.diff(tedges_days), strictly_monotone=True) 

710 v_end = flow_cum[-1] 

711 n_edges = len(flow_cum) 

712 

713 # Finite support: drop the gamma tails (mass ~1e-13, far below the discretization error this 

714 # closed form replaces). Restricting the integral to [support_lo, support_hi] is what keeps 

715 # the per-bin flow-edge band -- and the gamma CDF evaluations -- bounded. 

716 tail = 1e-13 

717 support_lo = float(loc + gamma_dist.ppf(tail, alpha, scale=beta)) 

718 support_hi = float(loc + gamma_dist.ppf(1.0 - tail, alpha, scale=beta)) 

719 

720 # phi(v) = int_0^v T(w) dw with T the inverse cumulative-volume map (piecewise-linear); phi is 

721 # piecewise-quadratic with knots at the cumulative-volume edges, so the per-bin time integral of 

722 # tau is a difference of phi at the look-back/forward limits. With spinup="constant" the map is 

723 # extended past the record at the boundary flow rates (one anchor each end, padded by the largest 

724 # reach r*support_hi) so phi extrapolates the spin-up; with a float spinup it stays clipped to 

725 # [0, v_end] (the covered sub-mass only). 

726 phi_v, phi_t, phi_knot, phi_rate = _phi_setup( 

727 flow, flow_cum, tedges_days, extrapolate=extrapolate, pad=r * support_hi 

728 ) 

729 

730 cout_tedges_days = tedges_to_days(cout_tedges, ref=tedges[0]) 

731 vol_out = linear_interpolate(x_ref=tedges_days, y_ref=flow_cum, x_query=cout_tedges_days, left=np.nan, right=np.nan) 

732 good_all = np.isfinite(vol_out[:-1]) & np.isfinite(vol_out[1:]) & (vol_out[1:] > vol_out[:-1]) 

733 if not np.any(good_all): 

734 return np.full(n_out, np.nan) 

735 v_lo_all = np.where(good_all, vol_out[:-1], 0.0) 

736 v_hi_all = np.where(good_all, vol_out[1:], 1.0) 

737 

738 # Fixed band of flow edges each bin's look-back/forward sweep can cross over the supported 

739 # pore volumes. band_width is the global maximum so every tile shares one column layout; a 

740 # spurious column clips to the support and merely splits a quadratic piece without changing 

741 # its integral, so the band is an exact superset. 

742 a_min = v_lo_all - r * support_hi if sign < 0 else v_lo_all + r * support_lo 

743 a_max = v_hi_all - r * support_lo if sign < 0 else v_hi_all + r * support_hi 

744 jlo_all = np.clip(np.searchsorted(flow_cum, a_min, "left") - 1, 0, n_edges - 1) 

745 jhi_all = np.clip(np.searchsorted(flow_cum, a_max, "right"), 0, n_edges - 1) 

746 band_width = int((jhi_all - jlo_all).max()) + 1 

747 band = np.arange(band_width) 

748 

749 # Tile over output bins to bound peak memory. Each bin carries ~3*band_width pieces through a 

750 # (3 nodes, 3 phi-args) stack of 9 * n_pieces elements; size the tile to that element budget. 

751 n_pieces = 3 * band_width + 3 

752 tile = max(1, _max_tile_elements // (9 * n_pieces)) 

753 

754 result = np.full(n_out, np.nan) 

755 for t0 in range(0, n_out, tile): 

756 t1 = min(t0 + tile, n_out) 

757 nt = t1 - t0 

758 good = good_all[t0:t1] 

759 v_lo = v_lo_all[t0:t1] 

760 v_hi = v_hi_all[t0:t1] 

761 cols = np.clip(jlo_all[t0:t1, None] + band[None, :], 0, n_edges - 1) 

762 fcb = flow_cum[cols] 

763 vlo = v_lo[:, None] 

764 vhi = v_hi[:, None] 

765 

766 # Direction-specific breakpoint families: the V_p values where a phi argument (a clipped 

767 # look-back/forward limit) crosses a flow edge or a validity bound. Clip to the support 

768 # and sort to form the integration pieces (zero-width pieces contribute nothing). 

769 if sign < 0: 

770 cand = np.concatenate([(vhi - fcb) / r, (vlo - fcb) / r, fcb / r, vlo / r, vhi / r], axis=1) 

771 else: 

772 cand = np.concatenate( 

773 [(fcb - vlo) / r, (fcb - vhi) / r, (v_end - fcb) / r, (v_end - vlo) / r, (v_end - vhi) / r], axis=1 

774 ) 

775 np.clip(cand, support_lo, support_hi, out=cand) 

776 cand.sort(axis=1) 

777 edges = np.concatenate([np.full((nt, 1), support_lo), cand, np.full((nt, 1), support_hi)], axis=1) 

778 lo = edges[:, :-1] 

779 hi = edges[:, 1:] 

780 mid = 0.5 * (lo + hi) 

781 half = 0.5 * (hi - lo) 

782 nodes = np.stack([lo, mid, hi], axis=0) # (3, nt, n_pieces) 

783 

784 # G(V_p) at the three quadrature nodes (piece lo/mid/hi) as a difference of phi. The constant 

785 # term phi(v_hi)/phi(v_lo) does not depend on V_p, so evaluate it once per bin and batch only 

786 # the three V_p-dependent phi arguments. With spinup="constant" the full bin is integrated 

787 # against the extrapolated phi (length is the full bin width); with a float spinup the limits 

788 # clamp to the streamtube's covered sub-interval and the covered length renormalizes. 

789 if sign < 0: 

790 a_hi = vhi[None] - r * nodes 

791 if extrapolate: 

792 v_start = np.broadcast_to(vlo[None], a_hi.shape) 

793 a_lo = vlo[None] - r * nodes 

794 length = np.broadcast_to(vhi[None] - vlo[None], a_hi.shape) 

795 else: 

796 v_start = np.maximum(vlo[None], r * nodes) 

797 a_lo = np.maximum(vlo[None] - r * nodes, 0.0) 

798 length = np.maximum(vhi[None] - v_start, 0.0) 

799 phi_const = _eval_phi(v_hi, phi_v, phi_t, phi_knot, phi_rate) 

800 phi_stack = _eval_phi(np.stack([v_start, a_hi, a_lo]), phi_v, phi_t, phi_knot, phi_rate) 

801 g = phi_const[None, :, None] - phi_stack[0] - phi_stack[1] + phi_stack[2] 

802 else: 

803 a_lo = vlo[None] + r * nodes 

804 if extrapolate: 

805 v_stop = np.broadcast_to(vhi[None], a_lo.shape) 

806 a_hi = vhi[None] + r * nodes 

807 length = np.broadcast_to(vhi[None] - vlo[None], a_lo.shape) 

808 else: 

809 v_stop = np.minimum(vhi[None], v_end - r * nodes) 

810 a_hi = np.minimum(vhi[None] + r * nodes, v_end) 

811 length = np.maximum(v_stop - vlo[None], 0.0) 

812 phi_const = _eval_phi(v_lo, phi_v, phi_t, phi_knot, phi_rate) 

813 phi_stack = _eval_phi(np.stack([a_hi, a_lo, v_stop]), phi_v, phi_t, phi_knot, phi_rate) 

814 g = phi_stack[0] - phi_stack[1] - phi_stack[2] + phi_const[None, :, None] 

815 g = np.where(length > 0, g, 0.0) 

816 length = np.where(length > 0, length, 0.0) 

817 g_lo, g_mid, g_hi = g 

818 l_lo, l_mid, l_hi = length 

819 

820 # Gamma partial moments over each piece: one CDF per shape on the piece edges, then diff. 

821 # M1/M2 follow from the shifted-gamma partial-moment identities. 

822 cdf_edges = edges - loc 

823 f0 = gamma_dist.cdf(cdf_edges, alpha, scale=beta) 

824 f1 = gamma_dist.cdf(cdf_edges, alpha + 1, scale=beta) 

825 f2 = gamma_dist.cdf(cdf_edges, alpha + 2, scale=beta) 

826 m0 = np.diff(f0, axis=1) 

827 d1 = np.diff(f1, axis=1) 

828 m1 = alpha * beta * d1 + loc * m0 

829 d2 = np.diff(f2, axis=1) 

830 m2 = alpha * (alpha + 1) * beta**2 * d2 + 2 * loc * alpha * beta * d1 + loc * loc * m0 

831 

832 # Piece-centred gamma moments mu_k = int (x - mid)^k f over each piece. They are formed by 

833 # shifting the raw partial moments (m1 - mid m0, m2 - 2 mid m1 + mid^2 m0), which for R > 1 

834 # subtracts large near-equal terms (|mid| up to support_hi, with m1/m2 carrying the matching 

835 # mid powers) and catastrophically cancels on near-degenerate pieces. They are, however, 

836 # bounded purely by the piece geometry -- |x - mid| <= half over the piece -- so clip each to 

837 # its exact range. This is not a strict no-op on well-conditioned pieces (a raw value may sit 

838 # a rounding step outside the tight bound; clamping it back is a negligible correction). It 

839 # earns its place on the near-degenerate pieces, where the cancellation noise -- which 

840 # otherwise pairs with the 1/half^2 second-difference below to manufacture a spurious 

841 # contribution -- is forced back into the physically valid band. 

842 b0 = np.maximum(m0, 0.0) 

843 mu1 = np.clip(m1 - mid * m0, -half * b0, half * b0) 

844 mu2 = np.clip(m2 - 2 * mid * m1 + mid**2 * m0, 0.0, half**2 * b0) 

845 

846 # G is quadratic per piece, L linear; reconstruct each from its three nodes (Lagrange- 

847 # exact) centred at mid and contract with the moments: int (a(x-mid)^2 + b(x-mid) + c) f 

848 # = a mu2 + b mu1 + c mu0. Zero-width pieces (duplicate breakpoints) divide by half == 0; 

849 # their masked result is discarded by np.where. 

850 safe = half > 0 

851 with np.errstate(divide="ignore", invalid="ignore"): 

852 g_a = np.where(safe, (g_lo - 2 * g_mid + g_hi) / (2 * half**2), 0.0) 

853 g_b = np.where(safe, (g_hi - g_lo) / (2 * half), 0.0) 

854 l_b = np.where(safe, (l_hi - l_lo) / (2 * half), 0.0) 

855 int_g = g_a * mu2 + g_b * mu1 + g_mid * m0 

856 int_l = l_b * mu1 + l_mid * m0 

857 

858 num = int_g.sum(axis=1) 

859 den = int_l.sum(axis=1) 

860 covered = good & (den > 0) 

861 if spinup_threshold > 0.0: 

862 # float spinup: emit only where the covered fraction of the APVD reaches the threshold. 

863 # den is the covered sub-mass; (v_hi - v_lo) * m0.sum is the fully-covered reference (equal 

864 # up to reassociation ulps when fully covered). The relative slack keeps the strictest 

865 # threshold (1.0) from rejecting fully-covered bins on that float noise. 

866 covered &= den >= (spinup_threshold - _SPINUP_GATE_RELTOL) * (v_hi - v_lo) * m0.sum(axis=1) 

867 tile_result = np.full(nt, np.nan)