Coverage for src/gwtransport/logremoval.py: 96%
46 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""
2Log Removal Calculations for First-Order Decay Processes.
4This module provides utilities to calculate log removal values from first-order decay
5processes, including pathogen inactivation and radioactive decay. The module supports
6basic log removal calculations and parallel flow arrangements where multiple flow paths
7operate simultaneously.
9First-Order Decay Model
10-----------------------
11The log removal from any first-order decay process is:
13 Log Removal = log10_decay_rate * residence_time
15where ``log10_decay_rate`` has units [log10/day] and ``residence_time`` has units [days].
16This is equivalent to exponential decay ``C_out/C_in = 10^(-mu * t)``, where mu is the
17log10 decay rate and t is residence time. The natural-log decay rate constant lambda [1/day]
18is related to mu by ``lambda = mu * ln(10)``.
20This model applies to any process that follows first-order kinetics:
22- **Pathogen inactivation**: viruses, bacteria, and protozoa lose infectivity over time
23- **Radioactive decay**: isotopes used for groundwater dating (tritium, CFC, SF6)
24- **Chemical degradation**: first-order breakdown of contaminants
26Pathogen Removal in Bank Filtration
27------------------------------------
28For pathogen removal during soil passage, total removal consists of two distinct mechanisms
29(Schijven and Hassanizadeh, 2000):
311. **Inactivation (time-dependent)**: Pathogens lose infectivity over time through biological
32 decay. This follows first-order kinetics and is modeled by this module as
33 ``LR_decay = log10_decay_rate * residence_time``. The inactivation rate depends strongly
34 on temperature and pathogen type.
362. **Attachment (geometry-dependent)**: Pathogens are physically removed by adsorption to soil
37 grains and straining. This depends on aquifer geometry, distance, soil properties, and pH,
38 and is NOT modeled by this module. Users should add this component separately based on
39 site-specific data.
41Total log removal = LR_decay (this module) + LR_attachment (user-specified).
43At the Castricum dune recharge site, Schijven et al. (1999) found that attachment contributed
44approximately 97% of total MS2 removal, with inactivation contributing only 3%. Inactivation
45rates for common model viruses at 10 degrees C are typically 0.02-0.11 log10/day (Schijven and
46Hassanizadeh, 2000, Table 7).
48Gamma-distribution parameter notation
49-------------------------------------
50Several functions are parameterized by a gamma distribution. The parameter prefix marks
51*which* physical quantity is gamma-distributed, because two distinct quantities appear here:
53- ``rt_alpha`` / ``rt_beta`` / ``rt_loc`` (or the equivalent ``rt_mean`` / ``rt_std`` /
54 ``rt_loc``) parameterize the gamma distribution of the **residence time** (used by
55 :func:`gamma_pdf`, :func:`gamma_cdf`, :func:`gamma_mean`).
56- ``apv_alpha`` / ``apv_beta`` / ``apv_loc`` (or the equivalent ``apv_mean`` / ``apv_std`` /
57 ``apv_loc``) parameterize the gamma distribution of the **aquifer pore volume** (used by
58 :func:`gamma_find_flow_for_target_mean`).
60These prefixes are intentional and load-bearing: residence time and pore volume are different
61quantities, so a bare ``alpha`` / ``beta`` / ``loc`` would be ambiguous in this module. Both the
62shape/scale and the mean/std pairs are validated through :func:`gwtransport.gamma.parse_parameters`,
63so invalid parameters (e.g. a negative shape) raise ``ValueError`` rather than silently returning
64an unphysical result.
66Available functions:
68- :func:`residence_time_to_log_removal` - Multiply residence times [days] by a log10 decay rate
69 [log10/day] to get log removal, returning an array of the same shape as the input. Negative
70 residence times or a negative rate give negative log removal; the caller interprets the sign.
72- :func:`decay_rate_to_log10_decay_rate` - Convert a natural-log first-order decay rate constant
73 lambda [1/day] to a log10 decay rate mu [log10/day] via ``mu = lambda / ln(10)``.
75- :func:`log10_decay_rate_to_decay_rate` - Convert a log10 decay rate mu [log10/day] to a
76 natural-log first-order decay rate constant lambda [1/day] via ``lambda = mu * ln(10)``.
78- :func:`parallel_mean` - Combine the log removals of parallel flow paths into the single log
79 removal of their blended outflow, ``-log10(sum(F_i * 10^(-LR_i)))``. Flow fractions default to
80 an equal split and must sum to 1.0 along the reduction axis; ``axis=None`` reduces over the
81 flattened input and returns a scalar.
83- :func:`gamma_pdf` - Probability density of log removal when the residence time is (shifted)
84 gamma-distributed: a gamma density with shape ``rt_alpha``, scale ``log10_decay_rate * rt_beta``,
85 and location ``log10_decay_rate * rt_loc``, evaluated at the requested log removal values.
87- :func:`gamma_cdf` - Cumulative distribution of log removal for the same shifted-gamma residence
88 time, i.e. the fraction of the extracted water whose log removal does not exceed ``r``.
90- :func:`gamma_mean` - Effective (flow-weighted, concentration-mixing) mean log removal for
91 gamma-distributed residence time,
92 ``log10_decay_rate * rt_loc + rt_alpha * log10(1 + rt_beta * log10_decay_rate * ln(10))``. This
93 is the continuous counterpart of :func:`parallel_mean` and lies below the arithmetic mean
94 ``log10_decay_rate * (rt_alpha * rt_beta + rt_loc)``, because short-residence-time paths
95 dominate the mixed outflow.
97- :func:`gamma_find_flow_for_target_mean` - Invert :func:`gamma_mean` for the flow rate that
98 attains a target effective mean log removal, given a gamma-distributed **aquifer pore volume**
99 (hence the ``apv_`` prefix; dividing the pore volume by the flow gives the residence time).
100 Closed form when ``apv_loc == 0``, otherwise a bracketed :func:`scipy.optimize.brentq` root in
101 ``1 / flow``. Requires a positive ``target_mean`` and a positive ``log10_decay_rate``.
103This file is part of gwtransport which is released under AGPL-3.0 license.
104See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
105"""
107import numpy as np
108import numpy.typing as npt
109from scipy import optimize, stats
111from gwtransport.gamma import parse_parameters
114def residence_time_to_log_removal(
115 *, residence_times: npt.ArrayLike, log10_decay_rate: float
116) -> npt.NDArray[np.floating]:
117 """
118 Compute log removal given residence times and a log10 decay rate.
120 ``Log Removal = log10_decay_rate * residence_time``, equivalent to the exponential
121 decay ``C_out/C_in = 10^(-log10_decay_rate * residence_time)``.
123 Parameters
124 ----------
125 residence_times : array-like
126 Residence times (days), of any shape. Negative values produce negative log
127 removal (mathematical amplification); the caller interprets the sign.
128 log10_decay_rate : float
129 Log10 decay rate mu (log10/day). Negative values correspond to first-order
130 production rather than decay.
132 Returns
133 -------
134 log_removals : ndarray
135 Log removal values, same shape as ``residence_times``. Log 1 is a 90%
136 reduction, log 2 a 99% reduction, log 3 a 99.9% reduction.
138 See Also
139 --------
140 decay_rate_to_log10_decay_rate : Convert natural-log decay rate to log10 decay rate
141 gamma_mean : Effective mean log removal for gamma-distributed residence times
142 gwtransport.residence_time.full : Compute residence times from flow and pore volume
143 :ref:`concept-residence-time` : Time in aquifer determines pathogen contact time
145 Examples
146 --------
147 >>> from gwtransport.logremoval import residence_time_to_log_removal
148 >>> residence_time_to_log_removal(
149 ... residence_times=[10.0, 20.0, 50.0], log10_decay_rate=0.2
150 ... ) # doctest: +NORMALIZE_WHITESPACE
151 array([ 2., 4., 10.])
152 """
153 return log10_decay_rate * np.asarray(residence_times, dtype=float)
156def decay_rate_to_log10_decay_rate(decay_rate: float) -> float:
157 """
158 Convert a natural-log decay rate constant to a log10 decay rate: mu = lambda / ln(10).
160 Parameters
161 ----------
162 decay_rate : float
163 Natural-log first-order decay rate constant lambda (1/day), e.g. ``np.log(2) / half_life``.
165 Returns
166 -------
167 log10_decay_rate : float
168 Log10 decay rate mu (log10/day).
170 See Also
171 --------
172 log10_decay_rate_to_decay_rate : Inverse conversion
174 Examples
175 --------
176 >>> import numpy as np
177 >>> from gwtransport.logremoval import decay_rate_to_log10_decay_rate
178 >>> decay_rate_to_log10_decay_rate(np.log(2) / 30) # doctest: +ELLIPSIS
179 np.float64(0.01003...)
180 """
181 return decay_rate / np.log(10)
184def log10_decay_rate_to_decay_rate(log10_decay_rate: float) -> float:
185 """
186 Convert a log10 decay rate to a natural-log decay rate constant: lambda = mu * ln(10).
188 Parameters
189 ----------
190 log10_decay_rate : float
191 Log10 decay rate mu (log10/day).
193 Returns
194 -------
195 decay_rate : float
196 Natural-log first-order decay rate constant lambda (1/day).
198 See Also
199 --------
200 decay_rate_to_log10_decay_rate : Inverse conversion
202 Examples
203 --------
204 >>> from gwtransport.logremoval import log10_decay_rate_to_decay_rate
205 >>> log10_decay_rate_to_decay_rate(0.2) # doctest: +ELLIPSIS
206 np.float64(0.4605...)
207 """
208 return log10_decay_rate * np.log(10)
211def parallel_mean(
212 *, log_removals: npt.ArrayLike, flow_fractions: npt.ArrayLike | None = None, axis: int | None = None
213) -> np.floating | npt.NDArray[np.floating]:
214 """
215 Calculate the weighted average log removal for a system with parallel flows.
217 This function computes the overall log removal efficiency of a parallel
218 filtration system. If flow_fractions is not provided, it assumes equal
219 distribution of flow across all paths.
221 The calculation uses the formula:
223 Total Log Removal = -log10(sum(F_i * 10^(-LR_i)))
225 Where:
226 - F_i = fraction of flow through system i (decimal, sum to 1.0)
227 - LR_i = log removal of system i
229 Parameters
230 ----------
231 log_removals : array-like
232 Array of log removal values for each parallel flow.
233 Each value represents the log10 reduction of pathogens.
234 For multi-dimensional arrays, the parallel mean is computed along
235 the specified axis.
237 flow_fractions : array-like, optional
238 Array of flow fractions for each parallel flow.
239 Must sum to 1.0 along the specified axis and have compatible shape
240 with log_removals. If None, equal flow distribution is assumed
241 (default is None).
243 axis : int, optional
244 Axis along which to compute the parallel mean for multi-dimensional
245 arrays. If None, the reduction matches the way ``np.mean`` / ``np.sum``
246 treat ``axis=None``: the parallel mean is computed over the flattened
247 input (default is None).
249 Returns
250 -------
251 np.floating or ndarray
252 The combined log removal value for the parallel system. Returns a
253 scalar when axis=None, otherwise an array with the specified axis
254 removed.
256 Raises
257 ------
258 ValueError
259 If ``flow_fractions`` does not sum to 1.0 along the specified axis.
261 See Also
262 --------
263 residence_time_to_log_removal : Compute log removal from residence times
265 Notes
266 -----
267 Log removal is a logarithmic measure of pathogen reduction:
269 - Log 1 = 90% reduction
270 - Log 2 = 99% reduction
271 - Log 3 = 99.9% reduction
273 For parallel flows, the combined removal is typically less effective
274 than the best individual removal but better than the worst.
275 For systems in series, log removals would be summed directly.
277 Examples
278 --------
279 >>> import numpy as np
280 >>> from gwtransport.logremoval import parallel_mean
281 >>> # Three parallel streams with equal flow and log removals of 3, 4, and 5
282 >>> log_removals = np.array([3, 4, 5])
283 >>> parallel_mean(log_removals=log_removals)
284 np.float64(3.431798275933005)
286 >>> # Two parallel streams with weighted flow
287 >>> log_removals = np.array([3, 5])
288 >>> flow_fractions = np.array([0.7, 0.3])
289 >>> parallel_mean(log_removals=log_removals, flow_fractions=flow_fractions)
290 np.float64(3.153044674980176)
292 >>> # Multi-dimensional array: parallel mean along axis 1
293 >>> log_removals_2d = np.array([[3, 4, 5], [2, 3, 4]])
294 >>> parallel_mean(log_removals=log_removals_2d, axis=1)
295 array([3.43179828, 2.43179828])
296 """
297 log_removals = np.asarray(log_removals, dtype=float)
298 decimal_reductions = 10.0 ** (-log_removals)
299 if flow_fractions is None:
300 return -np.log10(np.mean(decimal_reductions, axis=axis))
301 flow_fractions = np.asarray(flow_fractions, dtype=float)
302 if not np.all(np.isclose(np.sum(flow_fractions, axis=axis), 1.0)):
303 msg = "flow_fractions must sum to 1.0 (along the specified axis)"
304 raise ValueError(msg)
305 return -np.log10(np.sum(flow_fractions * decimal_reductions, axis=axis))
308def gamma_pdf(
309 *,
310 r: npt.ArrayLike,
311 rt_alpha: float | None = None,
312 rt_beta: float | None = None,
313 rt_loc: float = 0.0,
314 rt_mean: float | None = None,
315 rt_std: float | None = None,
316 log10_decay_rate: float,
317) -> npt.NDArray[np.floating]:
318 """
319 Compute the PDF of log removal given (shifted) gamma-distributed residence time.
321 With residence time ``T = T0 + rt_loc`` where ``T0 ~ Gamma(rt_alpha, rt_beta)``,
322 the log removal ``R = mu * T`` follows a shifted gamma distribution with shape
323 ``rt_alpha``, scale ``mu * rt_beta``, and location ``mu * rt_loc``.
325 The residence-time distribution is specified with either ``(rt_alpha, rt_beta)`` or
326 ``(rt_mean, rt_std)`` (optionally shifted by ``rt_loc``); both are routed through
327 :func:`gwtransport.gamma.parse_parameters`.
329 Parameters
330 ----------
331 r : array-like
332 Log removal values at which to compute the PDF.
333 rt_alpha : float, optional
334 Shape parameter of the gamma distribution for residence time. Must be positive.
335 rt_beta : float, optional
336 Scale parameter of the gamma distribution for residence time (days). Must be positive.
337 rt_loc : float, optional
338 Location (minimum residence time, days) of the residence time distribution.
339 Must be non-negative. Default is ``0.0``.
340 rt_mean : float, optional
341 Mean residence time (days). Alternative to ``rt_alpha``; supply with ``rt_std``.
342 Must be strictly greater than ``rt_loc``.
343 rt_std : float, optional
344 Standard deviation of the residence time (days). Alternative to ``rt_beta``;
345 supply with ``rt_mean``. Must be positive.
346 log10_decay_rate : float
347 Log10 decay rate mu (log10/day). Relates residence time to
348 log removal via R = mu * T.
350 Returns
351 -------
352 pdf : ndarray
353 PDF values corresponding to the input r values.
355 Raises
356 ------
357 ValueError
358 If parameter validation in :func:`gwtransport.gamma.parse_parameters` fails
359 (e.g. ``rt_loc`` negative, non-positive shape/scale, or neither/both
360 parameter pairs supplied).
362 See Also
363 --------
364 gamma_cdf : Cumulative distribution function of log removal
365 gamma_mean : Mean of the log removal distribution
366 """
367 rt_alpha, rt_beta, rt_loc = parse_parameters(mean=rt_mean, std=rt_std, loc=rt_loc, alpha=rt_alpha, beta=rt_beta)
368 return stats.gamma.pdf(r, a=rt_alpha, loc=log10_decay_rate * rt_loc, scale=log10_decay_rate * rt_beta)
371def gamma_cdf(
372 *,
373 r: npt.ArrayLike,
374 rt_alpha: float | None = None,
375 rt_beta: float | None = None,
376 rt_loc: float = 0.0,
377 rt_mean: float | None = None,
378 rt_std: float | None = None,
379 log10_decay_rate: float,
380) -> npt.NDArray[np.floating]:
381 """
382 Compute the CDF of log removal given (shifted) gamma-distributed residence time.
384 With residence time ``T = T0 + rt_loc`` where ``T0 ~ Gamma(rt_alpha, rt_beta)``,
385 the CDF is ``P(R <= r) = P(mu*(T0 + rt_loc) <= r) =
386 P(T0 <= (r - mu*rt_loc)/mu)`` which is the CDF of a shifted gamma distribution
387 with location ``mu * rt_loc``.
389 The residence-time distribution is specified with either ``(rt_alpha, rt_beta)`` or
390 ``(rt_mean, rt_std)`` (optionally shifted by ``rt_loc``); both are routed through
391 :func:`gwtransport.gamma.parse_parameters`.
393 Parameters
394 ----------
395 r : array-like
396 Log removal values at which to compute the CDF.
397 rt_alpha : float, optional
398 Shape parameter of the gamma distribution for residence time. Must be positive.
399 rt_beta : float, optional
400 Scale parameter of the gamma distribution for residence time (days). Must be positive.
401 rt_loc : float, optional
402 Location (minimum residence time, days) of the residence time distribution.
403 Must be non-negative. Default is ``0.0``.
404 rt_mean : float, optional
405 Mean residence time (days). Alternative to ``rt_alpha``; supply with ``rt_std``.
406 Must be strictly greater than ``rt_loc``.
407 rt_std : float, optional
408 Standard deviation of the residence time (days). Alternative to ``rt_beta``;
409 supply with ``rt_mean``. Must be positive.
410 log10_decay_rate : float
411 Log10 decay rate mu (log10/day). Relates residence time to
412 log removal via R = mu * T.
414 Returns
415 -------
416 cdf : ndarray
417 CDF values corresponding to the input r values.
419 Raises
420 ------
421 ValueError
422 If parameter validation in :func:`gwtransport.gamma.parse_parameters` fails
423 (e.g. ``rt_loc`` negative, non-positive shape/scale, or neither/both
424 parameter pairs supplied).
426 See Also
427 --------
428 gamma_pdf : Probability density function of log removal
429 gamma_mean : Mean of the log removal distribution
430 """
431 rt_alpha, rt_beta, rt_loc = parse_parameters(mean=rt_mean, std=rt_std, loc=rt_loc, alpha=rt_alpha, beta=rt_beta)
432 return stats.gamma.cdf(r, a=rt_alpha, loc=log10_decay_rate * rt_loc, scale=log10_decay_rate * rt_beta)
435def gamma_mean(
436 *,
437 rt_alpha: float | None = None,
438 rt_beta: float | None = None,
439 rt_loc: float = 0.0,
440 rt_mean: float | None = None,
441 rt_std: float | None = None,
442 log10_decay_rate: float,
443) -> float:
444 """
445 Compute the effective (parallel) mean log removal for (shifted) gamma-distributed residence time.
447 When water travels through multiple flow paths with gamma-distributed
448 residence times, the effective log removal is determined by mixing the
449 output concentrations (not by averaging individual log removals). For a
450 shifted gamma distribution ``T = T0 + rt_loc`` with ``T0 ~ Gamma(alpha, beta)``,
451 factoring the moment generating function gives:
453 LR_eff = -log10(E[10^(-mu*T)])
454 = -log10(10^(-mu*rt_loc) * E[10^(-mu*T0)])
455 = mu * rt_loc + alpha * log10(1 + beta * mu * ln(10))
457 The ``rt_loc`` term shifts the whole log-removal distribution by a constant
458 ``mu * rt_loc``; the alpha/beta term is unchanged. This is always less than
459 the arithmetic mean ``mu * (alpha * beta + rt_loc)`` because short residence
460 time paths contribute disproportionately to the output concentration.
462 The residence-time distribution is specified with either ``(rt_alpha, rt_beta)`` or
463 ``(rt_mean, rt_std)`` (optionally shifted by ``rt_loc``); both are routed through
464 :func:`gwtransport.gamma.parse_parameters`.
466 Parameters
467 ----------
468 rt_alpha : float, optional
469 Shape parameter of the gamma distribution for residence time. Must be positive.
470 rt_beta : float, optional
471 Scale parameter of the gamma distribution for residence time (days). Must be positive.
472 rt_loc : float, optional
473 Location (minimum residence time, days) of the residence time distribution.
474 Must be non-negative. Default is ``0.0``.
475 rt_mean : float, optional
476 Mean residence time (days). Alternative to ``rt_alpha``; supply with ``rt_std``.
477 Must be strictly greater than ``rt_loc``.
478 rt_std : float, optional
479 Standard deviation of the residence time (days). Alternative to ``rt_beta``;
480 supply with ``rt_mean``. Must be positive.
481 log10_decay_rate : float
482 Log10 decay rate mu (log10/day).
484 Returns
485 -------
486 mean : float
487 Effective (parallel) mean log removal value.
489 Raises
490 ------
491 ValueError
492 If parameter validation in :func:`gwtransport.gamma.parse_parameters` fails
493 (e.g. ``rt_loc`` negative, non-positive shape/scale, or neither/both
494 parameter pairs supplied).
496 See Also
497 --------
498 gamma_find_flow_for_target_mean : Find flow for target mean log removal
499 parallel_mean : Discrete version of this calculation
500 gamma_pdf : PDF of the log removal distribution
501 gamma_cdf : CDF of the log removal distribution
502 :ref:`concept-pore-volume-distribution` : Why residence times are distributed
503 """
504 rt_alpha, rt_beta, rt_loc = parse_parameters(mean=rt_mean, std=rt_std, loc=rt_loc, alpha=rt_alpha, beta=rt_beta)
505 return log10_decay_rate * rt_loc + rt_alpha * np.log1p(rt_beta * log10_decay_rate * np.log(10)) / np.log(10)
508def gamma_find_flow_for_target_mean(
509 *,
510 target_mean: float,
511 apv_alpha: float | None = None,
512 apv_beta: float | None = None,
513 apv_loc: float = 0.0,
514 apv_mean: float | None = None,
515 apv_std: float | None = None,
516 log10_decay_rate: float,
517) -> float:
518 """
519 Find the flow rate that produces a target effective mean log removal.
521 Given a (shifted) gamma-distributed aquifer pore volume with parameters
522 ``(apv_alpha, apv_beta, apv_loc)``, the residence time distribution at flow
523 ``Q`` is a shifted gamma with shape ``apv_alpha``, scale ``apv_beta/Q``, and
524 location ``apv_loc/Q``. From :func:`gamma_mean`:
526 LR_eff = mu * apv_loc / Q + apv_alpha * log10(1 + (apv_beta/Q) * mu * ln(10))
528 For ``apv_loc == 0`` this is closed-form:
530 Q = apv_beta * mu * ln(10) / (10^(target_mean / apv_alpha) - 1)
532 For ``apv_loc > 0`` the equation is transcendental and solved numerically
533 with :func:`scipy.optimize.brentq` by bracketing the root in ``1/Q``.
535 The pore-volume distribution is specified with either ``(apv_alpha, apv_beta)`` or
536 ``(apv_mean, apv_std)`` (optionally shifted by ``apv_loc``); both are routed through
537 :func:`gwtransport.gamma.parse_parameters`.
539 Parameters
540 ----------
541 target_mean : float
542 Target effective mean log removal value. Must be positive.
543 apv_alpha : float, optional
544 Shape parameter of the gamma distribution for aquifer pore volume. Must be positive.
545 apv_beta : float, optional
546 Scale parameter of the gamma distribution for aquifer pore volume. Must be positive.
547 apv_loc : float, optional
548 Location (minimum aquifer pore volume) of the gamma distribution.
549 Must be non-negative. Default is ``0.0``.
550 apv_mean : float, optional
551 Mean aquifer pore volume. Alternative to ``apv_alpha``; supply with ``apv_std``.
552 Must be strictly greater than ``apv_loc``.
553 apv_std : float, optional
554 Standard deviation of the aquifer pore volume. Alternative to ``apv_beta``;
555 supply with ``apv_mean``. Must be positive.
556 log10_decay_rate : float
557 Log10 decay rate mu (log10/day).
559 Returns
560 -------
561 flow : float
562 Flow rate (same units as apv_beta per day) that produces the
563 target mean log removal.
565 Raises
566 ------
567 ValueError
568 If ``target_mean`` is not positive, if ``log10_decay_rate`` is not positive
569 (no decay can never produce a positive target log removal), or if parameter
570 validation in :func:`gwtransport.gamma.parse_parameters` fails (e.g. ``apv_loc``
571 negative, non-positive shape/scale, or neither/both parameter pairs supplied).
573 See Also
574 --------
575 gamma_mean : Compute effective mean log removal for given parameters
576 """
577 apv_alpha, apv_beta, apv_loc = parse_parameters(
578 mean=apv_mean, std=apv_std, loc=apv_loc, alpha=apv_alpha, beta=apv_beta
579 )
580 if target_mean <= 0:
581 msg = "target_mean must be positive"
582 raise ValueError(msg)
583 if log10_decay_rate <= 0:
584 # Without decay, the effective mean log removal is identically zero
585 # regardless of flow, so no finite flow can attain a positive target.
586 msg = "log10_decay_rate must be positive to attain a positive target_mean"
587 raise ValueError(msg)
589 ln10 = np.log(10)
590 flow_closed_form = apv_beta * log10_decay_rate * ln10 / (10 ** (target_mean / apv_alpha) - 1)
592 if apv_loc == 0.0:
593 return float(flow_closed_form)
595 # Solve target = mu*apv_loc*u + apv_alpha*log10(1 + apv_beta*mu*ln(10)*u) for u = 1/flow.
596 # Both terms are monotonically increasing in u, so f(u) - target is monotonic with a
597 # unique positive root. Bracket: at u = 1/flow_closed_form the alpha/beta term alone
598 # equals target_mean, so the full f overshoots by exactly mu*apv_loc*u_upper > 0.
599 u_upper = 1.0 / flow_closed_form
601 def residual(u: float) -> float:
602 return float(
603 log10_decay_rate * apv_loc * u
604 + apv_alpha * np.log1p(apv_beta * log10_decay_rate * ln10 * u) / ln10
605 - target_mean
606 )
608 u_root = optimize.brentq(residual, 0.0, u_upper)
609 return 1.0 / float(u_root) # type: ignore[arg-type]