Total emission

The basin total emitted mass from a TMTPP posterior, as posterior samples split into detected and inferred-missing mass, with the mass-rate over the observed windows. M_total = Λ E_f[Q] D needs a mean event duration D in the reciprocal of the flux's time unit, so the mass is in kg. The rate is M_total / T_eff, with T_eff the length of the union of the observed windows, the same windows Λ integrates over. Tier V.

totals

Total-emission estimator (implementation in plumax.population._src.totals).

TotalEmission dataclass

TotalEmission(samples_kg: Float[Array, ' member'], detected_kg: Float[Array, ' member'], missing_kg: Float[Array, ' member'], t_effective_s: float)

Posterior samples of a total emission, split into detected and missing mass.

samples_kg = detected_kg + missing_kg draw by draw.

Attributes:

Name	Type	Description
samples_kg	Float[Array, ' member']	M_total, the un-thinned total mass [kg].
detected_kg	Float[Array, ' member']	Λ E_f[Q P_d(Q)] D, the mass the model expects to be detected [kg].
missing_kg	Float[Array, ' member']	Λ E_f[Q (1 − P_d(Q))] D, the inferred missing mass [kg].
t_effective_s	float	T_eff, the length of the union of the observed windows [s].

rate_tg_per_yr

rate_tg_per_yr() -> Float[Array, ' member']

Mean emission rate while observed, M_total / T_eff [Tg yr⁻¹].

Returns:

Type	Description
Float[Array, ' member']	Posterior samples of the rate, over a Julian year.

Source code in src/plumax/population/_src/totals.py

def rate_tg_per_yr(self) -> Float[Array, " member"]:
    """Mean emission rate while observed, M_total / T_eff [Tg yr⁻¹].

    Returns:
        Posterior samples of the rate, over a Julian year.
    """
    # M_total / T_eff [kg s⁻¹]
    rate_kg_per_s = self.samples_kg / self.t_effective_s
    return rate_kg_per_s * SECONDS_PER_YEAR / KG_PER_TG

total_emission

total_emission(posterior: TMTPPPosterior | Mapping[str, TMTPPPosterior], *, horizon: tuple[float, float] | Sequence[tuple[float, float]], event_duration: float, log_pod: Callable[[Array], Array] | None, q_min: float, q_max: float, n_nodes: int = 2001, nodes_per_day: int = NODES_PER_DAY, pods: Sequence[InstrumentPOD] | None = None, tod_origin: float = 0.0) -> TotalEmission

Posterior samples of the basin total, with the inferred missing mass.

M_total⁽ᵉ⁾ = Λ⁽ᵉ⁾ E_f⁽ᵉ⁾[Q] D, M_missing⁽ᵉ⁾ = Λ⁽ᵉ⁾ E_f⁽ᵉ⁾[Q (1 − P_d(Q))] D.

E_f[Q] is the mark's truncated_mean, the partial first moment over [q_min, q_max], not a mean conditional on that range. See the module docstring for the mass and mass-rate convention.

Algorithm

1. Merge horizon into disjoint observed windows; T_eff is their length.
2. For each draw, Λ over those windows, from the kernel's closed form or Simpson's rule.
3. For each draw, E_f[Q] in closed form and E_f[Q (1 − P_d)] by log-space Simpson over [q_min, q_max].
4. Detected mass is the difference, floored at zero against round-off.
5. Sum strata draw by draw: their posteriors are independent, so any pairing of their draws is a draw of the sum.

Parameters:

Name	Type	Description	Default
posterior	TMTPPPosterior | Mapping[str, TMTPPPosterior]	A TMTPP posterior, or one per stratum (the design's v1 spatial shortcut: facility classes that sum to the basin).	required
horizon	tuple[float, float] | Sequence[tuple[float, float]]	The observed window (t₀, t₁) or a sequence of them [days]; they may overlap.	required
event_duration	float	Mean event duration D > 0, in the reciprocal of Q's time unit.	required
log_pod	Callable[[Array], Array] | None	log P_d on rates in the marks' units, as :func:~plumax.population.fitting.tmtpp_iw_model takes it, or None for certain detection. A fixed union is log_pod_union(pods).	required
q_min	float	Lower limit of the mark integrals, > 0.	required
q_max	float	Upper limit of the mark integrals, finite and > q_min.	required
n_nodes	int	Simpson nodes of the missing-mass integral.	2001
nodes_per_day	int	Quadrature density of Λ for kernels without a closed form. Lower it for long windows when the kernel is smooth.	NODES_PER_DAY
pods	Sequence[InstrumentPOD] | None	The instruments of a fit with a calibration prior; requires log_pod=None. Each draw then uses its own union P_d^∪, from its pod.<name>.* sites. Every instrument counts as looking throughout horizon. With tod_bins, the missing mass sums over the hours: Σ_b Λ_b E_f[Q (1 − P_d(Q; b))] D.	None
tod_origin	float	Local hour of day at t = 0 [hours], for tod_bins.	0.0

Returns:

Name	Type	Description
The	TotalEmission	class:TotalEmission.

Raises:

Type	Description
ValueError	For an empty window, invalid mark limits or duration, no strata, strata with different numbers of draws, or both log_pod and pods.

Examples:

Two draws of a constant rate, where a tenth of the events are detected:

>>> import math
>>> import jax.numpy as jnp
>>> from plumax.population.fitting import TMTPPPosterior
>>> from plumax.population.marks import LognormalMark
>>> samples = {
...     "lambda_0": jnp.array([1.0, 2.0]),
...     "mark.mu": jnp.zeros(2),
...     "mark.sigma": jnp.ones(2),
... }
>>> posterior = TMTPPPosterior("constant", LognormalMark(), samples)
>>> total = total_emission(
...     posterior,
...     horizon=(0.0, 10.0),
...     event_duration=2.0,
...     log_pod=lambda q: math.log(0.1) + jnp.zeros_like(q),
...     q_min=1e-6,
...     q_max=1e6,
... )
>>> # M_total = λ T e^{1/2} D
>>> expected = jnp.array([1.0, 2.0]) * 10.0 * math.exp(0.5) * 2.0
>>> bool(jnp.allclose(total.samples_kg, expected))
True
>>> bool(jnp.allclose(total.detected_kg, 0.1 * expected, rtol=1e-4))
True
>>> total.t_effective_s
864000.0

Source code in src/plumax/population/_src/totals.py

def total_emission(
    posterior: TMTPPPosterior | Mapping[str, TMTPPPosterior],
    *,
    horizon: tuple[float, float] | Sequence[tuple[float, float]],
    event_duration: float,
    log_pod: Callable[[Array], Array] | None,
    q_min: float,
    q_max: float,
    n_nodes: int = 2001,
    nodes_per_day: int = NODES_PER_DAY,
    pods: Sequence[InstrumentPOD] | None = None,
    tod_origin: float = 0.0,
) -> TotalEmission:
    """Posterior samples of the basin total, with the inferred missing mass.

    M_total⁽ᵉ⁾ = Λ⁽ᵉ⁾ E_f⁽ᵉ⁾[Q] D,  M_missing⁽ᵉ⁾ = Λ⁽ᵉ⁾ E_f⁽ᵉ⁾[Q (1 − P_d(Q))] D.

    E_f[Q] is the mark's ``truncated_mean``, the partial first moment over
    [q_min, q_max], not a mean conditional on that range. See the module docstring for
    the mass and mass-rate convention.

    Algorithm:
        1. Merge ``horizon`` into disjoint observed windows; T_eff is their length.
        2. For each draw, Λ over those windows, from the kernel's closed form or
           Simpson's rule.
        3. For each draw, E_f[Q] in closed form and E_f[Q (1 − P_d)] by log-space
           Simpson over [q_min, q_max].
        4. Detected mass is the difference, floored at zero against round-off.
        5. Sum strata draw by draw: their posteriors are independent, so any pairing
           of their draws is a draw of the sum.

    Args:
        posterior: A TMTPP posterior, or one per stratum (the design's v1 spatial
            shortcut: facility classes that sum to the basin).
        horizon: The observed window ``(t₀, t₁)`` or a sequence of them [days]; they
            may overlap.
        event_duration: Mean event duration D > 0, in the reciprocal of Q's time unit.
        log_pod: log P_d on rates in the marks' units, as
            :func:`~plumax.population.fitting.tmtpp_iw_model` takes it, or ``None``
            for certain detection. A fixed union is ``log_pod_union(pods)``.
        q_min: Lower limit of the mark integrals, > 0.
        q_max: Upper limit of the mark integrals, finite and > ``q_min``.
        n_nodes: Simpson nodes of the missing-mass integral.
        nodes_per_day: Quadrature density of Λ for kernels without a closed form.
            Lower it for long windows when the kernel is smooth.
        pods: The instruments of a fit with a calibration prior; requires
            ``log_pod=None``. Each draw then uses its own union P_d^∪, from its
            ``pod.<name>.*`` sites. Every instrument counts as looking throughout
            ``horizon``. With ``tod_bins``, the missing mass sums over the hours:
            Σ_b Λ_b E_f[Q (1 − P_d(Q; b))] D.
        tod_origin: Local hour of day at t = 0 [hours], for ``tod_bins``.

    Returns:
        The :class:`TotalEmission`.

    Raises:
        ValueError: For an empty window, invalid mark limits or duration, no strata,
            strata with different numbers of draws, or both ``log_pod`` and ``pods``.

    Examples:
        Two draws of a constant rate, where a tenth of the events are detected:

        >>> import math
        >>> import jax.numpy as jnp
        >>> from plumax.population.fitting import TMTPPPosterior
        >>> from plumax.population.marks import LognormalMark
        >>> samples = {
        ...     "lambda_0": jnp.array([1.0, 2.0]),
        ...     "mark.mu": jnp.zeros(2),
        ...     "mark.sigma": jnp.ones(2),
        ... }
        >>> posterior = TMTPPPosterior("constant", LognormalMark(), samples)
        >>> total = total_emission(
        ...     posterior,
        ...     horizon=(0.0, 10.0),
        ...     event_duration=2.0,
        ...     log_pod=lambda q: math.log(0.1) + jnp.zeros_like(q),
        ...     q_min=1e-6,
        ...     q_max=1e6,
        ... )
        >>> # M_total = λ T e^{1/2} D
        >>> expected = jnp.array([1.0, 2.0]) * 10.0 * math.exp(0.5) * 2.0
        >>> bool(jnp.allclose(total.samples_kg, expected))
        True
        >>> bool(jnp.allclose(total.detected_kg, 0.1 * expected, rtol=1e-4))
        True
        >>> total.t_effective_s
        864000.0
    """
    if log_pod is not None and pods is not None:
        raise ValueError("total_emission: pass log_pod or pods, not both.")
    strata = (
        [posterior]
        if isinstance(posterior, TMTPPPosterior)
        else list(posterior.values())
    )
    if not strata:
        raise ValueError("total_emission: no strata given.")
    if len({stratum.num_samples for stratum in strata}) != 1:
        raise ValueError(
            "total_emission: strata need the same number of draws to sum, got "
            f"{[stratum.num_samples for stratum in strata]}."
        )
    if not 0.0 < q_min < q_max < float("inf"):
        raise ValueError(
            f"total_emission: need 0 < q_min < q_max < inf, got ({q_min}, {q_max})."
        )
    if not event_duration > 0.0:
        raise ValueError(
            f"total_emission: event_duration must be > 0, got {event_duration}."
        )
    windows = _merge_windows(horizon)
    parts = [
        _stratum_emission(
            stratum,
            windows=windows,
            event_duration=event_duration,
            log_pod=log_pod,
            pods=pods,
            q_min=q_min,
            q_max=q_max,
            n_nodes=n_nodes,
            nodes_per_day=nodes_per_day,
            tod_origin=tod_origin,
        )
        for stratum in strata
    ]
    # M_total = Σ_strata M_total,stratum
    samples_kg = sum(total for total, _ in parts)
    missing_kg = sum(missing for _, missing in parts)
    # T_eff = |∪_w w|
    effective_days = sum(stop - start for start, stop in windows)
    return TotalEmission(
        samples_kg=samples_kg,
        detected_kg=jnp.maximum(samples_kg - missing_kg, 0.0),
        missing_kg=missing_kg,
        t_effective_s=effective_days * SECONDS_PER_DAY,
    )