from __future__ import annotations
from dataclasses import dataclass, field
from itertools import product
from typing import Iterable, Optional
import numpy as np
import pandas as pd
from .models import MatrixPopulationModel
from .perturbations import Target
from .simulation import SimulationResult, simulate_dynamics
[docs]
@dataclass(frozen=True)
class PerturbationGrid:
"""Discrete perturbation space used for vulnerability calculations.
Parameters
----------
magnitudes : iterable of float
Proportional reductions applied to the selected demographic component.
Values must lie in the interval ``[0, 1]``.
durations : iterable of int
Number of consecutive projection intervals during which each
perturbation event is active.
periods : iterable of int
Number of projection intervals between the onset of consecutive
perturbation events.
Notes
-----
The package represents recurrence using ``period`` rather than frequency
because simulations are performed in discrete projection intervals. Grids can
also be constructed from event frequencies with :meth:`from_frequencies`.
"""
magnitudes: Iterable[float]
durations: Iterable[int]
periods: Iterable[int]
def __post_init__(self) -> None:
mags = np.asarray(list(self.magnitudes), dtype=float)
durations = np.asarray(list(self.durations), dtype=int)
periods = np.asarray(list(self.periods), dtype=int)
if mags.size == 0 or durations.size == 0 or periods.size == 0:
raise ValueError("magnitudes, durations and periods must be non-empty.")
if np.any((mags < 0) | (mags > 1)):
raise ValueError("magnitudes must lie in [0, 1].")
if np.any(durations < 0):
raise ValueError("durations must be non-negative integers.")
if np.any(periods <= 0):
raise ValueError("periods must be positive integers.")
object.__setattr__(self, "magnitudes", tuple(float(x) for x in mags))
object.__setattr__(self, "durations", tuple(int(x) for x in durations))
object.__setattr__(self, "periods", tuple(int(x) for x in periods))
[docs]
@classmethod
def from_frequencies(
cls,
*,
magnitudes: Iterable[float],
durations: Iterable[int],
frequencies: Iterable[float],
generation_time: float = 1.0,
rounding: str = "nearest",
) -> "PerturbationGrid":
"""Create a perturbation grid from event frequencies.
Parameters
----------
magnitudes : iterable of float
Proportional reductions applied to the selected demographic
component.
durations : iterable of int
Number of consecutive projection intervals during which each event
is active.
frequencies : iterable of float
Event frequencies, interpreted as the number of events per
``generation_time`` projection intervals.
generation_time : float, default=1.0
Number of projection intervals defining the reference time unit for
``frequencies``. For example, if frequencies are specified as events
per generation, this should be the generation time expressed in
projection intervals.
rounding : {"nearest", "floor", "ceil"}, default="nearest"
Rule used to convert non-integer recurrence periods into integer
projection intervals.
Returns
-------
PerturbationGrid
Grid with recurrence periods computed as
``period = generation_time / frequency``.
Notes
-----
Rounded duplicate periods are removed. The minimum period is one
projection interval.
"""
if generation_time <= 0:
raise ValueError("generation_time must be positive.")
periods = []
for f in frequencies:
f = float(f)
if f <= 0:
raise ValueError("frequencies must be positive.")
raw = generation_time / f
if rounding == "nearest":
p = int(round(raw))
elif rounding == "floor":
p = int(np.floor(raw))
elif rounding == "ceil":
p = int(np.ceil(raw))
else:
raise ValueError("rounding must be 'nearest', 'floor', or 'ceil'.")
periods.append(max(1, p))
return cls(magnitudes=magnitudes, durations=durations, periods=sorted(set(periods)))
[docs]
def scenarios(self, *, skip_infeasible: bool = True):
"""Iterate over perturbation scenarios in the grid.
Parameters
----------
skip_infeasible : bool, default=True
If True, omit regimes where ``duration > period``. If False, yield
these regimes with ``feasible=False``.
Yields
------
tuple
Tuple ``(magnitude, duration, period, feasible)``.
"""
for magnitude, duration, period in product(self.magnitudes, self.durations, self.periods):
feasible = period >= duration
if feasible or not skip_infeasible:
yield float(magnitude), int(duration), int(period), bool(feasible)
[docs]
@dataclass(frozen=True)
class GridResult:
"""Result of a perturbation-grid simulation.
Attributes
----------
table : pandas.DataFrame
Long-format table with one row per perturbation regime. Standard columns
include ``target``, ``magnitude``, ``duration``, ``period``,
``feasible``, ``population_reduction``, ``final_population`` and
``baseline_final_population``.
trajectories : dict
Dictionary mapping ``(magnitude, duration, period)`` tuples to
:class:`~demovuln.simulation.SimulationResult` objects. This dictionary
is populated only when ``return_trajectories=True`` is used in
:func:`run_grid`.
"""
table: pd.DataFrame
trajectories: dict[tuple[float, int, int], SimulationResult] = field(default_factory=dict)
@property
def vulnerability(self) -> float:
"""Integrated vulnerability, computed as mean population reduction."""
return compute_vulnerability(self.table)
[docs]
def compute_vulnerability(table: pd.DataFrame, *, column: str = "population_reduction") -> float:
"""Compute integrated vulnerability from grid-simulation output.
Parameters
----------
table : pandas.DataFrame
Long-format output table returned by :func:`run_grid`.
column : str, default="population_reduction"
Name of the column containing percent population reduction values.
Returns
-------
float
Mean population reduction across perturbation regimes. Missing values
are ignored, so infeasible regimes encoded as ``NaN`` do not contribute.
Raises
------
ValueError
If ``column`` is not present in ``table``.
"""
if column not in table:
raise ValueError(f"{column!r} not found in table.")
values = pd.to_numeric(table[column], errors="coerce")
return float(values.mean())
[docs]
def run_grid(
model: MatrixPopulationModel | np.ndarray,
*,
target: Target,
grid: PerturbationGrid,
t_max: int,
recovery_steps: int = 0,
start: int = 0,
initial_state: Optional[np.ndarray] = None,
normalize_by_lambda: bool = True,
survival_affects_fecundity: bool = True,
custom_mask: Optional[np.ndarray] = None,
return_trajectories: bool = False,
skip_infeasible: bool = True,
force_during_recovery: bool = False,
) -> GridResult:
"""Evaluate a perturbation grid and compute demographic vulnerability.
Parameters
----------
model : MatrixPopulationModel or array_like of shape (n_stages, n_stages)
Unperturbed matrix population model.
target : {"adult_survival", "juvenile_survival", "fecundity", "all", "custom"}
Demographic component affected by the perturbations.
grid : PerturbationGrid
Discrete perturbation space, defined by magnitudes, durations and
recurrence periods.
t_max : int
Number of projection intervals in the perturbation-forcing window.
recovery_steps : int, default=0
Number of additional unperturbed projection intervals after the forcing
window.
start : int, default=0
Projection interval at which the first perturbation event starts.
initial_state : array_like of shape (n_stages,), optional
Initial population vector. If omitted, the stable stage distribution of
the unperturbed model is used.
normalize_by_lambda : bool, default=True
If True, divide baseline and perturbed projection matrices by the
dominant eigenvalue of the unperturbed matrix.
survival_affects_fecundity : bool, default=True
If True, survival perturbations scale whole source-stage columns,
including fecundity entries.
custom_mask : array_like of bool, optional
Boolean matrix selecting entries to perturb when ``target="custom"``.
return_trajectories : bool, default=False
If True, store the full :class:`~demovuln.simulation.SimulationResult`
object for every feasible regime.
skip_infeasible : bool, default=True
If True, omit regimes where ``duration > period``. If False, retain them
in the output table with ``feasible=False`` and missing response values.
force_during_recovery : bool, default=False
If True, scheduled perturbation events continue during the recovery
window.
Returns
-------
GridResult
Object containing the full grid table and, optionally, individual
simulation trajectories.
Notes
-----
The integrated vulnerability metric is available as
``result.vulnerability`` and is computed as the mean percent population
reduction across feasible perturbation regimes.
"""
if not isinstance(model, MatrixPopulationModel):
model = MatrixPopulationModel(np.asarray(model, dtype=float))
rows: list[dict] = []
trajectories: dict[tuple[float, int, int], SimulationResult] = {}
for magnitude, duration, period, feasible in grid.scenarios(skip_infeasible=skip_infeasible):
if not feasible:
rows.append(
{
"target": target,
"magnitude": magnitude,
"duration": duration,
"period": period,
"frequency": np.nan,
"feasible": False,
"population_reduction": np.nan,
"final_population": np.nan,
"baseline_final_population": np.nan,
}
)
continue
sim = simulate_dynamics(
model,
target=target,
magnitude=magnitude,
duration=duration,
period=period,
t_max=t_max,
recovery_steps=recovery_steps,
start=start,
initial_state=initial_state,
normalize_by_lambda=normalize_by_lambda,
survival_affects_fecundity=survival_affects_fecundity,
custom_mask=custom_mask,
return_stage_vectors=return_trajectories,
force_during_recovery=force_during_recovery,
)
rows.append(
{
"target": target,
"magnitude": magnitude,
"duration": duration,
"period": period,
"frequency": 1.0 / period,
"feasible": True,
"population_reduction": sim.reduction,
"final_population": sim.final_population,
"baseline_final_population": sim.baseline_final_population,
}
)
if return_trajectories:
trajectories[(magnitude, duration, period)] = sim
return GridResult(table=pd.DataFrame(rows), trajectories=trajectories)