Individual-Based Model (IBM) API Reference

The IBM module tracks populations of super-individuals with explicit body size, age, energy reserves, and spatial location. Each super-individual represents a cohort of real organisms. The IBM integrates into Ecosim's derivative loop: IBM groups override the standard ODE derivative for their functional group while the rest of the food web runs normally.

See the IBM Parameterization Guide for practical tuning advice.

Base Classes

Core data structures: SuperIndividual, IBMGroup (abstract base), IBMStepResult, and SpatialContext.

pypath.ibm.base

Base data structures for the IBM (Individual-Based Model) module.

Provides the foundational dataclasses and abstract base class used by all IBM group implementations in PyPath. These structures define how individual organisms (super-individuals) are represented and how IBM groups interface with the Ecosim population dynamics engine.

Classes:

Name	Description
SuperIndividual	Represents a cohort of biologically identical organisms.
IBMStepResult	Return type for a single IBM integration step.
IBMGroup	Abstract base class that all concrete IBM implementations must subclass.

LifeStage

Bases: IntEnum

Life stage codes used by IBM super-individuals.

Source code in pypath/ibm/base.py

class LifeStage(IntEnum):
    """Life stage codes used by IBM super-individuals."""

    EGG = 0
    YOLK_SAC = 1
    LARVA = 2
    JUVENILE = 3
    ADULT = 4

SuperIndividual dataclass

A super-individual representing a cohort of identical organisms.

Each super-individual tracks the state of n_represented biologically identical fish (or other organisms). This is the fundamental unit of the IBM module.

Parameters:

Name	Type	Description	Default
id	int	Unique identifier for this super-individual.	required
n_represented	float	Number of real individuals this super-individual represents.	required
weight	float	Individual body weight (grams).	required
length	float	Individual body length (cm).	required
age	float	Age (years).	required
energy_reserve	float	Energy reserve (grams-equivalent, initialized as weight * 0.1).	required
patch_idx	int	Spatial patch index where this super-individual currently resides.	required
is_mature	bool	Whether this super-individual has reached sexual maturity.	required
sex	int	Sex code (0 = female, 1 = male, or other coding as needed).	required
life_stage	int	Life stage code (0=egg, 1=yolk_sac, 2=larva, 3=juvenile, 4=adult).	4
degree_days	float	Accumulated thermal development in degree-days (°C·day), egg stage only.	0.0
starvation_days	float	Consecutive days without sufficient feeding.	0.0
yolk_energy_kj	float	Yolk energy in kJ (only meaningful for life_stage 0-1).	0.0

Source code in pypath/ibm/base.py

@dataclass
class SuperIndividual:
    """A super-individual representing a cohort of identical organisms.

    Each super-individual tracks the state of *n_represented* biologically
    identical fish (or other organisms). This is the fundamental unit of
    the IBM module.

    Parameters
    ----------
    id : int
        Unique identifier for this super-individual.
    n_represented : float
        Number of real individuals this super-individual represents.
    weight : float
        Individual body weight (grams).
    length : float
        Individual body length (cm).
    age : float
        Age (years).
    energy_reserve : float
        Energy reserve (grams-equivalent, initialized as weight * 0.1).
    patch_idx : int
        Spatial patch index where this super-individual currently resides.
    is_mature : bool
        Whether this super-individual has reached sexual maturity.
    sex : int
        Sex code (0 = female, 1 = male, or other coding as needed).
    life_stage : int
        Life stage code (0=egg, 1=yolk_sac, 2=larva, 3=juvenile, 4=adult).
    degree_days : float
        Accumulated thermal development in degree-days (°C·day), egg stage only.
    starvation_days : float
        Consecutive days without sufficient feeding.
    yolk_energy_kj : float
        Yolk energy in kJ (only meaningful for life_stage 0-1).
    """

    id: int
    n_represented: float
    weight: float
    length: float
    age: float
    energy_reserve: float
    patch_idx: int
    is_mature: bool
    sex: int
    life_stage: int = 4  # 0=egg, 1=yolk_sac, 2=larva, 3=juvenile, 4=adult
    degree_days: float = 0.0  # accumulated thermal development (°C·day), egg stage only
    starvation_days: float = 0.0  # consecutive days without sufficient feeding
    yolk_energy_kj: float = (
        0.0  # yolk energy in kJ (only meaningful for life_stage 0-1)
    )

    def total_biomass_tonnes(self) -> float:
        """Return total biomass represented by this super-individual in tonnes.

        Computed as ``n_represented * weight / 1e6``, assuming weight is in
        grams and 1 tonne = 1e6 grams.

        Returns
        -------
        float
            Total biomass in metric tonnes.
        """
        return self.n_represented * self.weight / 1e6

total_biomass_tonnes

total_biomass_tonnes() -> float

Return total biomass represented by this super-individual in tonnes.

Computed as n_represented * weight / 1e6, assuming weight is in grams and 1 tonne = 1e6 grams.

Returns:

Type	Description
float	Total biomass in metric tonnes.

Source code in pypath/ibm/base.py

def total_biomass_tonnes(self) -> float:
    """Return total biomass represented by this super-individual in tonnes.

    Computed as ``n_represented * weight / 1e6``, assuming weight is in
    grams and 1 tonne = 1e6 grams.

    Returns
    -------
    float
        Total biomass in metric tonnes.
    """
    return self.n_represented * self.weight / 1e6

IBMStepResult dataclass

Result of a single IBM integration step.

Returned by :meth:IBMGroup.compute_step to communicate the IBM state back to the Ecosim integration loop.

Parameters:

Name	Type	Description	Default
biomass	float	Total group biomass (tonnes) after this step.	required
production	float	Net production during this step (tonnes).	required
consumption_by_prey	ndarray	1-D array of shape (n_groups,) giving the biomass consumed from each prey group during this step.	required
mortality_count	float	Number of individuals that died during this step.	required
recruitment_count	float	Number of new individuals recruited during this step.	required

Source code in pypath/ibm/base.py

@dataclass
class IBMStepResult:
    """Result of a single IBM integration step.

    Returned by :meth:`IBMGroup.compute_step` to communicate the IBM
    state back to the Ecosim integration loop.

    Parameters
    ----------
    biomass : float
        Total group biomass (tonnes) after this step.
    production : float
        Net production during this step (tonnes).
    consumption_by_prey : np.ndarray
        1-D array of shape ``(n_groups,)`` giving the biomass consumed
        from each prey group during this step.
    mortality_count : float
        Number of individuals that died during this step.
    recruitment_count : float
        Number of new individuals recruited during this step.
    """

    biomass: float
    production: float
    consumption_by_prey: np.ndarray
    mortality_count: float
    recruitment_count: float
    patch_biomass: Optional[np.ndarray] = None

SpatialContext dataclass

Spatial data passed to IBM groups during Ecospace simulations.

When an IBM group is part of a spatial simulation, this context provides the patch-level environmental information needed for movement decisions.

Parameters:

Name	Type	Description	Default
adjacency	Any	Sparse adjacency matrix of shape (n_patches, n_patches). Typically scipy.sparse.csr_matrix.	required
habitat_quality	ndarray	Per-patch habitat quality for this group, shape (n_patches,).	required
food_density	ndarray	Per-patch total prey biomass, shape (n_patches,).	required
predator_density	ndarray	Per-patch total predator biomass, shape (n_patches,).	required
n_patches	int	Number of spatial patches.	required

Source code in pypath/ibm/base.py

@dataclass
class SpatialContext:
    """Spatial data passed to IBM groups during Ecospace simulations.

    When an IBM group is part of a spatial simulation, this context provides
    the patch-level environmental information needed for movement decisions.

    Parameters
    ----------
    adjacency : Any
        Sparse adjacency matrix of shape ``(n_patches, n_patches)``.
        Typically ``scipy.sparse.csr_matrix``.
    habitat_quality : np.ndarray
        Per-patch habitat quality for this group, shape ``(n_patches,)``.
    food_density : np.ndarray
        Per-patch total prey biomass, shape ``(n_patches,)``.
    predator_density : np.ndarray
        Per-patch total predator biomass, shape ``(n_patches,)``.
    n_patches : int
        Number of spatial patches.
    """

    adjacency: Any  # scipy.sparse.csr_matrix (avoid hard import)
    habitat_quality: np.ndarray
    food_density: np.ndarray
    predator_density: np.ndarray
    n_patches: int

IBMGroup

Bases: ABC

Abstract base class for IBM group implementations.

Every IBM-managed functional group must subclass IBMGroup and implement all four abstract methods. The Ecosim integration loop calls these methods to delegate dynamics to the IBM engine while keeping the rest of the food-web coupled.

Parameters:

Name	Type	Description	Default
group_index	int	One-based index of this group in the Ecopath/Ecosim model (0 is reserved for the "Outside" placeholder).	required
n_groups	int	Total number of functional groups in the model (used to size consumption arrays).	required

Attributes:

Name	Type	Description
group_index	int	One-based index of this group.
n_groups	int	Total number of groups in the model.
individuals	List[SuperIndividual]	Population of super-individuals managed by this group.

Source code in pypath/ibm/base.py

class IBMGroup(ABC):
    """Abstract base class for IBM group implementations.

    Every IBM-managed functional group must subclass ``IBMGroup`` and
    implement all four abstract methods. The Ecosim integration loop
    calls these methods to delegate dynamics to the IBM engine while
    keeping the rest of the food-web coupled.

    Parameters
    ----------
    group_index : int
        One-based index of this group in the Ecopath/Ecosim model
        (0 is reserved for the "Outside" placeholder).
    n_groups : int
        Total number of functional groups in the model (used to size
        consumption arrays).

    Attributes
    ----------
    group_index : int
        One-based index of this group.
    n_groups : int
        Total number of groups in the model.
    individuals : List[SuperIndividual]
        Population of super-individuals managed by this group.
    """

    def __init__(self, group_index: int, n_groups: int) -> None:
        self.group_index: int = group_index
        self.n_groups: int = n_groups
        self.individuals: List[SuperIndividual] = []

    @abstractmethod
    def compute_step(
        self,
        prey_available: np.ndarray,
        predation_pressure: float,
        env_forcing: Dict[str, Any],
        dt: float,
        spatial_context: Optional["SpatialContext"] = None,
    ) -> IBMStepResult:
        """Advance the IBM population by one time step.

        Parameters
        ----------
        prey_available : np.ndarray
            1-D array of shape ``(n_groups,)`` giving the available biomass
            of each prey group.
        predation_pressure : float
            Total predation mortality pressure on this group from other
            predators in the Ecosim food web.
        env_forcing : Dict[str, Any]
            Dictionary of environmental forcing values (temperature, etc.).
        dt : float
            Time step size (years).
        spatial_context : SpatialContext, optional
            Spatial patch data for Ecospace simulations. When ``None``
            (default), no spatial movement is performed.

        Returns
        -------
        IBMStepResult
            Aggregated results of this time step.
        """

    @abstractmethod
    def get_aggregate_biomass(self) -> float:
        """Return the total biomass (tonnes) of all super-individuals.

        Returns
        -------
        float
            Sum of biomass across all super-individuals.
        """

    @abstractmethod
    def get_consumption_by_prey(self) -> np.ndarray:
        """Return the consumption vector by prey group.

        Returns
        -------
        np.ndarray
            1-D array of shape ``(n_groups,)`` with biomass consumed from
            each prey group.
        """

    @abstractmethod
    def initialize_from_ecosim(
        self,
        biomass: float,
        params: Dict[str, Any],
        n_super_individuals: int = 500,
    ) -> None:
        """Initialize the IBM population from Ecosim equilibrium state.

        Parameters
        ----------
        biomass : float
            Initial total biomass (tonnes) from Ecopath.
        params : Dict[str, Any]
            Species-specific biological parameters.
        n_super_individuals : int, optional
            Number of super-individuals to create (default 500).
        """

compute_step abstractmethod

compute_step(prey_available: ndarray, predation_pressure: float, env_forcing: Dict[str, Any], dt: float, spatial_context: Optional['SpatialContext'] = None) -> IBMStepResult

Advance the IBM population by one time step.

Parameters:

Name	Type	Description	Default
prey_available	ndarray	1-D array of shape (n_groups,) giving the available biomass of each prey group.	required
predation_pressure	float	Total predation mortality pressure on this group from other predators in the Ecosim food web.	required
env_forcing	Dict[str, Any]	Dictionary of environmental forcing values (temperature, etc.).	required
dt	float	Time step size (years).	required
spatial_context	SpatialContext	Spatial patch data for Ecospace simulations. When None (default), no spatial movement is performed.	None

Returns:

Type	Description
IBMStepResult	Aggregated results of this time step.

Source code in pypath/ibm/base.py

@abstractmethod
def compute_step(
    self,
    prey_available: np.ndarray,
    predation_pressure: float,
    env_forcing: Dict[str, Any],
    dt: float,
    spatial_context: Optional["SpatialContext"] = None,
) -> IBMStepResult:
    """Advance the IBM population by one time step.

    Parameters
    ----------
    prey_available : np.ndarray
        1-D array of shape ``(n_groups,)`` giving the available biomass
        of each prey group.
    predation_pressure : float
        Total predation mortality pressure on this group from other
        predators in the Ecosim food web.
    env_forcing : Dict[str, Any]
        Dictionary of environmental forcing values (temperature, etc.).
    dt : float
        Time step size (years).
    spatial_context : SpatialContext, optional
        Spatial patch data for Ecospace simulations. When ``None``
        (default), no spatial movement is performed.

    Returns
    -------
    IBMStepResult
        Aggregated results of this time step.
    """

get_aggregate_biomass abstractmethod

get_aggregate_biomass() -> float

Return the total biomass (tonnes) of all super-individuals.

Returns:

Type	Description
float	Sum of biomass across all super-individuals.

Source code in pypath/ibm/base.py

@abstractmethod
def get_aggregate_biomass(self) -> float:
    """Return the total biomass (tonnes) of all super-individuals.

    Returns
    -------
    float
        Sum of biomass across all super-individuals.
    """

get_consumption_by_prey abstractmethod

get_consumption_by_prey() -> np.ndarray

Return the consumption vector by prey group.

Returns:

Type	Description
ndarray	1-D array of shape (n_groups,) with biomass consumed from each prey group.

Source code in pypath/ibm/base.py

@abstractmethod
def get_consumption_by_prey(self) -> np.ndarray:
    """Return the consumption vector by prey group.

    Returns
    -------
    np.ndarray
        1-D array of shape ``(n_groups,)`` with biomass consumed from
        each prey group.
    """

initialize_from_ecosim abstractmethod

initialize_from_ecosim(biomass: float, params: Dict[str, Any], n_super_individuals: int = 500) -> None

Initialize the IBM population from Ecosim equilibrium state.

Parameters:

Name	Type	Description	Default
biomass	float	Initial total biomass (tonnes) from Ecopath.	required
params	Dict[str, Any]	Species-specific biological parameters.	required
n_super_individuals	int	Number of super-individuals to create (default 500).	500

Source code in pypath/ibm/base.py

@abstractmethod
def initialize_from_ecosim(
    self,
    biomass: float,
    params: Dict[str, Any],
    n_super_individuals: int = 500,
) -> None:
    """Initialize the IBM population from Ecosim equilibrium state.

    Parameters
    ----------
    biomass : float
        Initial total biomass (tonnes) from Ecopath.
    params : Dict[str, Any]
        Species-specific biological parameters.
    n_super_individuals : int, optional
        Number of super-individuals to create (default 500).
    """

Bioenergetics

Wisconsin bioenergetics model: Q10 temperature scaling, allometric metabolism, assimilation, and growth. Drives individual weight and energy reserve changes each timestep.

Key functions:

• growth_step() — single-individual energy budget and weight update
• growth_step_batch() — vectorized batch processing for all individuals (uses NumPy arrays, significantly faster for large populations)
• q10_temperature_factor() — Q10 temperature scaling
• metabolism() — weight-dependent metabolic rate with temperature correction
• allometric_length() — weight-to-length conversion
• BioenergParams — dataclass with all bioenergetics parameters

pypath.ibm.bioenergetics

Wisconsin bioenergetics model for the IBM integration.

Implements the energy-budget approach used to drive individual fish growth in the Individual-Based Model. Each super-individual's weight and energy reserve are updated each timestep based on consumption, metabolism, specific dynamic action (SDA), and (for mature fish) reproduction costs.

The core equation follows the Wisconsin model framework:

net_energy = assimilated_consumption - metabolism - SDA - reproduction_cost

Temperature dependence is modelled with a Q10 formulation, and allometric scaling converts weight to length.

Functions:

Name	Description
q10_temperature_factor	Compute Q10 temperature scaling factor.
allometric_length	Convert weight to length using an allometric power law.
metabolism	Compute standard metabolic rate.
assimilation	Compute assimilated consumption.
growth_step	Advance weight and energy reserve by one timestep.

Classes:

Name	Description
BioenergParams	Dataclass holding all bioenergetics parameters for a species.

BioenergParams dataclass

Parameters for the Wisconsin bioenergetics model.

Holds species-specific constants that govern metabolism, assimilation, growth, and reproduction in the IBM.

Parameters:

Name	Type	Description	Default
ra	float	Metabolic rate intercept (g O2 / g fish / day at reference temperature).	required
rb	float	Metabolic rate weight exponent (typically negative, indicating per-gram metabolic rate decreases with body size).	required
q10	float	Q10 temperature coefficient -- the factor by which metabolic rate increases for every 10 degree C rise in temperature.	required
t_ref	float	Reference temperature (degrees C) at which ra was measured.	required
sda_fraction	float	Specific dynamic action fraction (0-1). The proportion of consumption allocated to the energetic cost of digestion.	required
unassimilated_fraction	float	Fraction of consumption that is not assimilated (0-1), i.e. lost as faeces and excretion.	required
a_length	float	Allometric coefficient for the weight-to-length conversion (length = a_length * weight ** b_length).	required
b_length	float	Allometric exponent for the weight-to-length conversion.	required
energy_density	float	Energy content per gram of fish tissue (kJ/g). Used to convert net energy (kJ) to weight change (g). Default is 5.0.	5.0
reproduction_fraction	float	Fraction of net surplus energy allocated to reproduction when the individual is mature. Default is 0.3.	0.3

Source code in pypath/ibm/bioenergetics.py

@dataclass
class BioenergParams:
    """Parameters for the Wisconsin bioenergetics model.

    Holds species-specific constants that govern metabolism, assimilation,
    growth, and reproduction in the IBM.

    Parameters
    ----------
    ra : float
        Metabolic rate intercept (g O2 / g fish / day at reference temperature).
    rb : float
        Metabolic rate weight exponent (typically negative, indicating
        per-gram metabolic rate decreases with body size).
    q10 : float
        Q10 temperature coefficient -- the factor by which metabolic rate
        increases for every 10 degree C rise in temperature.
    t_ref : float
        Reference temperature (degrees C) at which ``ra`` was measured.
    sda_fraction : float
        Specific dynamic action fraction (0-1). The proportion of
        consumption allocated to the energetic cost of digestion.
    unassimilated_fraction : float
        Fraction of consumption that is not assimilated (0-1),
        i.e. lost as faeces and excretion.
    a_length : float
        Allometric coefficient for the weight-to-length conversion
        (length = a_length * weight ** b_length).
    b_length : float
        Allometric exponent for the weight-to-length conversion.
    energy_density : float
        Energy content per gram of fish tissue (kJ/g). Used to convert
        net energy (kJ) to weight change (g). Default is 5.0.
    reproduction_fraction : float
        Fraction of net surplus energy allocated to reproduction when
        the individual is mature. Default is 0.3.
    """

    ra: float
    rb: float
    q10: float
    t_ref: float
    sda_fraction: float
    unassimilated_fraction: float
    a_length: float
    b_length: float
    energy_density: float = 5.0
    reproduction_fraction: float = 0.3

    def __post_init__(self):
        assert self.energy_density > 0, (
            f"energy_density must be > 0, got {self.energy_density}"
        )

q10_temperature_factor

q10_temperature_factor(temp: float, t_ref: float, q10: float) -> float

Compute the Q10 temperature scaling factor.

The Q10 model scales a rate measured at t_ref to a new temperature temp using the formula:

factor = q10 ** ((temp - t_ref) / 10.0)

Parameters:

Name	Type	Description	Default
temp	float	Current water temperature (degrees C).	required
t_ref	float	Reference temperature (degrees C).	required
q10	float	Q10 coefficient (dimensionless, typically 1.5 -- 3.0).	required

Returns:

Type	Description
float	Multiplicative scaling factor (1.0 at t_ref).

Source code in pypath/ibm/bioenergetics.py

def q10_temperature_factor(temp: float, t_ref: float, q10: float) -> float:
    """Compute the Q10 temperature scaling factor.

    The Q10 model scales a rate measured at ``t_ref`` to a new temperature
    ``temp`` using the formula:

        factor = q10 ** ((temp - t_ref) / 10.0)

    Parameters
    ----------
    temp : float
        Current water temperature (degrees C).
    t_ref : float
        Reference temperature (degrees C).
    q10 : float
        Q10 coefficient (dimensionless, typically 1.5 -- 3.0).

    Returns
    -------
    float
        Multiplicative scaling factor (1.0 at ``t_ref``).
    """
    return q10 ** ((temp - t_ref) / 10.0)

allometric_length

allometric_length(weight: float, a: float, b: float) -> float

Convert body weight to body length using an allometric power law.

Computes length = a * weight ** b. Returns 0.0 if weight is zero or negative.

Parameters:

Name	Type	Description	Default
weight	float	Individual body weight (grams).	required
a	float	Allometric coefficient.	required
b	float	Allometric exponent.	required

Returns:

Type	Description
float	Body length (cm), or 0.0 if weight <= 0.

Source code in pypath/ibm/bioenergetics.py

def allometric_length(weight: float, a: float, b: float) -> float:
    """Convert body weight to body length using an allometric power law.

    Computes ``length = a * weight ** b``. Returns 0.0 if weight is zero
    or negative.

    Parameters
    ----------
    weight : float
        Individual body weight (grams).
    a : float
        Allometric coefficient.
    b : float
        Allometric exponent.

    Returns
    -------
    float
        Body length (cm), or 0.0 if weight <= 0.
    """
    if weight <= 0.0:
        return 0.0
    return a * weight**b

metabolism

metabolism(weight: float, temperature: float, params: BioenergParams) -> float

Compute the standard metabolic rate for an individual.

Uses the allometric form with Q10 temperature dependence:

rate = ra * weight^rb * q10_factor

The result is in the same units as ra (g O2 / g fish / day) and represents the per-gram daily metabolic cost scaled for temperature.

Parameters:

Name	Type	Description	Default
weight	float	Individual body weight (grams).	required
temperature	float	Current water temperature (degrees C).	required
params	BioenergParams	Bioenergetics parameters.	required

Returns:

Type	Description
float	Metabolic rate (g O2 / g fish / day, temperature-adjusted).

Source code in pypath/ibm/bioenergetics.py

def metabolism(weight: float, temperature: float, params: BioenergParams) -> float:
    """Compute the standard metabolic rate for an individual.

    Uses the allometric form with Q10 temperature dependence:

        rate = ra * weight^rb * q10_factor

    The result is in the same units as ``ra`` (g O2 / g fish / day) and
    represents the per-gram daily metabolic cost scaled for temperature.

    Parameters
    ----------
    weight : float
        Individual body weight (grams).
    temperature : float
        Current water temperature (degrees C).
    params : BioenergParams
        Bioenergetics parameters.

    Returns
    -------
    float
        Metabolic rate (g O2 / g fish / day, temperature-adjusted).
    """
    q10_factor = q10_temperature_factor(temperature, params.t_ref, params.q10)
    return params.ra * (weight**params.rb) * q10_factor

assimilation

assimilation(consumption: float, params: BioenergParams) -> float

Compute the assimilated portion of consumption.

Removes the unassimilated fraction (faeces + excretion) from total consumption:

assimilated = consumption * (1 - unassimilated_fraction)

Parameters:

Name	Type	Description	Default
consumption	float	Total consumption (energy or mass units).	required
params	BioenergParams	Bioenergetics parameters.	required

Returns:

Type	Description
float	Assimilated consumption.

Source code in pypath/ibm/bioenergetics.py

def assimilation(consumption: float, params: BioenergParams) -> float:
    """Compute the assimilated portion of consumption.

    Removes the unassimilated fraction (faeces + excretion) from
    total consumption:

        assimilated = consumption * (1 - unassimilated_fraction)

    Parameters
    ----------
    consumption : float
        Total consumption (energy or mass units).
    params : BioenergParams
        Bioenergetics parameters.

    Returns
    -------
    float
        Assimilated consumption.
    """
    return consumption * (1.0 - params.unassimilated_fraction)

growth_step

growth_step(weight: float, energy_reserve: float, consumption: float, temperature: float, is_mature: bool, dt: float, params: BioenergParams) -> tuple[float, float]

Advance an individual's weight and energy reserve by one timestep.

Implements the Wisconsin bioenergetics budget for a single integration step. The energy budget is:

assim = assimilation(consumption)
sda   = consumption * sda_fraction
met   = metabolism(weight, temperature) * dt * 365
net   = assim - met - sda
reproduction_cost = net * reproduction_fraction  (if mature and net > 0)
weight_change = net / energy_density

The dt * 365 factor converts the daily metabolic rate (ra) to the appropriate fraction of a year represented by dt.

Surplus energy is added to the energy reserve; deficits drain the reserve first before reducing body weight. Weight is clamped to a minimum of 0.1 grams.

Parameters:

Name	Type	Description	Default
weight	float	Current individual body weight (grams).	required
energy_reserve	float	Current energy reserve (dimensionless index).	required
consumption	float	Total consumption during this timestep.	required
temperature	float	Current water temperature (degrees C).	required
is_mature	bool	Whether the individual has reached sexual maturity.	required
dt	float	Timestep size (fraction of a year).	required
params	BioenergParams	Bioenergetics parameters.	required

Returns:

Type	Description
tuple[float, float]	(new_weight, new_energy_reserve) after this timestep.

Source code in pypath/ibm/bioenergetics.py

def growth_step(
    weight: float,
    energy_reserve: float,
    consumption: float,
    temperature: float,
    is_mature: bool,
    dt: float,
    params: BioenergParams,
) -> tuple[float, float]:
    """Advance an individual's weight and energy reserve by one timestep.

    Implements the Wisconsin bioenergetics budget for a single integration
    step. The energy budget is:

        assim = assimilation(consumption)
        sda   = consumption * sda_fraction
        met   = metabolism(weight, temperature) * dt * 365
        net   = assim - met - sda
        reproduction_cost = net * reproduction_fraction  (if mature and net > 0)
        weight_change = net / energy_density

    The ``dt * 365`` factor converts the daily metabolic rate (``ra``) to
    the appropriate fraction of a year represented by ``dt``.

    Surplus energy is added to the energy reserve; deficits drain the
    reserve first before reducing body weight. Weight is clamped to a
    minimum of 0.1 grams.

    Parameters
    ----------
    weight : float
        Current individual body weight (grams).
    energy_reserve : float
        Current energy reserve (dimensionless index).
    consumption : float
        Total consumption during this timestep.
    temperature : float
        Current water temperature (degrees C).
    is_mature : bool
        Whether the individual has reached sexual maturity.
    dt : float
        Timestep size (fraction of a year).
    params : BioenergParams
        Bioenergetics parameters.

    Returns
    -------
    tuple[float, float]
        ``(new_weight, new_energy_reserve)`` after this timestep.
    """
    # Assimilated energy
    assim = assimilation(consumption, params)

    # Specific dynamic action (cost of digestion)
    sda = consumption * params.sda_fraction

    # Metabolic cost: ra is daily, dt is in years, so multiply by dt * 365
    met = metabolism(weight, temperature, params) * dt * 365.0

    # Net energy balance
    net_energy = assim - met - sda

    # Reproduction cost for mature fish with positive surplus
    repro_cost = 0.0
    if is_mature and net_energy > 0.0:
        repro_cost = net_energy * params.reproduction_fraction
        net_energy -= repro_cost

    # Convert net energy to weight change
    weight_change = net_energy / params.energy_density

    # Update weight and energy reserve
    new_weight = weight + weight_change

    # Update energy reserve: surplus increases it, deficit drains it
    new_energy_reserve = energy_reserve + net_energy / params.energy_density

    # Enforce minimum weight
    new_weight = max(new_weight, 0.1)

    # Energy reserve cannot go below zero
    new_energy_reserve = max(new_energy_reserve, 0.0)

    return (float(new_weight), float(new_energy_reserve))

growth_step_batch

growth_step_batch(weights: ndarray, energy_reserves: ndarray, consumptions: ndarray, temperature: 'float | np.ndarray', is_mature: ndarray, dt: float, params: BioenergParams) -> tuple[np.ndarray, np.ndarray]

Vectorized growth step for all individuals at once.

Parameters:

Name	Type	Description	Default
weights	ndarray	Body weights (grams), shape (n,).	required
energy_reserves	ndarray	Energy reserves, shape (n,).	required
consumptions	ndarray	Total consumption per individual, shape (n,).	required
temperature	float or ndarray	Water temperature (degrees C). Scalar or per-individual array.	required
is_mature	ndarray	Boolean array, shape (n,).	required
dt	float	Timestep (fraction of a year).	required
params	BioenergParams	Bioenergetics parameters.	required

Returns:

Type	Description
tuple[ndarray, ndarray]	(new_weights, new_energy_reserves).

Source code in pypath/ibm/bioenergetics.py

def growth_step_batch(
    weights: np.ndarray,
    energy_reserves: np.ndarray,
    consumptions: np.ndarray,
    temperature: "float | np.ndarray",
    is_mature: np.ndarray,
    dt: float,
    params: BioenergParams,
) -> tuple[np.ndarray, np.ndarray]:
    """Vectorized growth step for all individuals at once.

    Parameters
    ----------
    weights : np.ndarray
        Body weights (grams), shape ``(n,)``.
    energy_reserves : np.ndarray
        Energy reserves, shape ``(n,)``.
    consumptions : np.ndarray
        Total consumption per individual, shape ``(n,)``.
    temperature : float or np.ndarray
        Water temperature (degrees C). Scalar or per-individual array.
    is_mature : np.ndarray
        Boolean array, shape ``(n,)``.
    dt : float
        Timestep (fraction of a year).
    params : BioenergParams
        Bioenergetics parameters.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        ``(new_weights, new_energy_reserves)``.
    """
    assim = consumptions * (1 - params.unassimilated_fraction)
    sda = consumptions * params.sda_fraction
    q10_factor = q10_temperature_factor(temperature, params.t_ref, params.q10)
    met = params.ra * (weights**params.rb) * q10_factor * dt * 365.0
    net_energy = assim - met - sda

    # Reproduction cost for mature fish with positive surplus
    repro_cost = np.where(
        np.asarray(is_mature, dtype=bool) & (net_energy > 0),
        net_energy * params.reproduction_fraction,
        0.0,
    )
    net_energy = net_energy - repro_cost

    weight_change = net_energy / params.energy_density
    new_weights = np.maximum(weights + weight_change, 0.1)
    new_energy_reserves = np.maximum(energy_reserves + weight_change, 0.0)

    return new_weights, new_energy_reserves

thornton_lessem

thornton_lessem(temp: float, CQ: float, CTO: float, CTM: float, CTL: float, CK1: float, CK4: float) -> float

Thornton-Lessem temperature function (Fish Bioenergetics 3.0/4.0).

Produces a dome-shaped temperature response curve widely used in fish bioenergetics models (Thornton & Lessem 1978). The result f(T) ranges from ~CK1 at CQ (cold extreme) to ~0.98 near CTO-CTM (optimum) to ~CK4 at CTL (hot extreme).

Parameters:

Name	Type	Description	Default
temp	float	Current water temperature (degrees C).	required
CQ	float	Lower temperature where rate = CK1 fraction of maximum (~T_min).	required
CTO	float	Temperature where ascending limb reaches 0.98 of max (~T_opt).	required
CTM	float	Temperature where descending limb is still 0.98 of max (> CTO).	required
CTL	float	Upper temperature where rate = CK4 fraction of maximum (~T_max).	required
CK1	float	Small fraction of max at CQ (typically 0.01-0.05).	required
CK4	float	Small fraction of max at CTL (typically 0.01-0.05).	required

Returns:

Type	Description
float	Temperature scaling factor in [0, ~0.98], dome-shaped.

Source code in pypath/ibm/bioenergetics.py

def thornton_lessem(
    temp: float,
    CQ: float,
    CTO: float,
    CTM: float,
    CTL: float,
    CK1: float,
    CK4: float,
) -> float:
    """Thornton-Lessem temperature function (Fish Bioenergetics 3.0/4.0).

    Produces a dome-shaped temperature response curve widely used in fish
    bioenergetics models (Thornton & Lessem 1978).  The result f(T) ranges
    from ~CK1 at CQ (cold extreme) to ~0.98 near CTO-CTM (optimum) to
    ~CK4 at CTL (hot extreme).

    Parameters
    ----------
    temp : float
        Current water temperature (degrees C).
    CQ : float
        Lower temperature where rate = CK1 fraction of maximum (~T_min).
    CTO : float
        Temperature where ascending limb reaches 0.98 of max (~T_opt).
    CTM : float
        Temperature where descending limb is still 0.98 of max (> CTO).
    CTL : float
        Upper temperature where rate = CK4 fraction of maximum (~T_max).
    CK1 : float
        Small fraction of max at CQ (typically 0.01-0.05).
    CK4 : float
        Small fraction of max at CTL (typically 0.01-0.05).

    Returns
    -------
    float
        Temperature scaling factor in [0, ~0.98], dome-shaped.
    """
    # Parameter validation
    if not (0 < CK1 < 1):
        raise ValueError(f"CK1 must be in (0, 1), got {CK1}")
    if not (0 < CK4 < 1):
        raise ValueError(f"CK4 must be in (0, 1), got {CK4}")
    if CTO <= CQ:
        raise ValueError(f"CTO must be > CQ, got CTO={CTO}, CQ={CQ}")
    if CTL <= CTM:
        raise ValueError(f"CTL must be > CTM, got CTL={CTL}, CTM={CTM}")

    if temp < CQ or temp > CTL:
        return 0.0
    G1 = (1.0 / (CTO - CQ)) * math.log(0.98 * (1.0 - CK1) / (CK1 * 0.02))
    G2 = (1.0 / (CTL - CTM)) * math.log(0.98 * (1.0 - CK4) / (CK4 * 0.02))
    L1 = math.exp(min(G1 * (temp - CQ), 700))
    L2 = math.exp(min(G2 * (CTL - temp), 700))
    K_A = (CK1 * L1) / (1.0 + CK1 * (L1 - 1.0))
    K_B = (CK4 * L2) / (1.0 + CK4 * (L2 - 1.0))
    result = K_A * K_B
    if not math.isfinite(result):
        return 0.0
    return max(0.0, result)

oxygen_scalar

oxygen_scalar(o2: float, pcrit: float) -> float

Compute oxygen limitation scalar.

Returns 1.0 when dissolved oxygen is at or above the critical threshold (pcrit), and scales linearly to 0.0 as oxygen approaches zero.

Parameters:

Name	Type	Description	Default
o2	float	Dissolved oxygen concentration (mg/L).	required
pcrit	float	Critical oxygen threshold (mg/L) below which limitation begins.	required

Returns:

Type	Description
float	Oxygen scalar in [0, 1].

Source code in pypath/ibm/bioenergetics.py

def oxygen_scalar(o2: float, pcrit: float) -> float:
    """Compute oxygen limitation scalar.

    Returns 1.0 when dissolved oxygen is at or above the critical threshold
    (``pcrit``), and scales linearly to 0.0 as oxygen approaches zero.

    Parameters
    ----------
    o2 : float
        Dissolved oxygen concentration (mg/L).
    pcrit : float
        Critical oxygen threshold (mg/L) below which limitation begins.

    Returns
    -------
    float
        Oxygen scalar in [0, 1].
    """
    if o2 >= pcrit:
        return 1.0
    return max(0.0, o2 / pcrit)

growth_step_batch_ontogenetic

growth_step_batch_ontogenetic(weights: ndarray, energy_reserves: ndarray, consumptions: ndarray, temperature: 'float | np.ndarray', is_mature: ndarray, dt: float, bioenerg_params: BioenergParams, larval_params: 'LarvalParams') -> tuple[np.ndarray, np.ndarray]

Vectorized growth step with ontogenetic (size-dependent) interpolation.

Extends the Wisconsin bioenergetics model with sigmoid-interpolated activity multiplier, assimilation efficiency, and basal metabolism split (Rs + Ra). At large (adult) body sizes the results converge to those of :func:growth_step_batch when larval_params.rs_a * (1 + am_max) == ra.

Parameters:

Name	Type	Description	Default
weights	ndarray	Body weights (grams), shape (n,).	required
energy_reserves	ndarray	Energy reserves, shape (n,).	required
consumptions	ndarray	Total consumption per individual, shape (n,).	required
temperature	float or ndarray	Water temperature (degrees C). Scalar or per-individual array.	required
is_mature	ndarray	Boolean array, shape (n,).	required
dt	float	Timestep (fraction of a year).	required
bioenerg_params	BioenergParams	Standard bioenergetics parameters.	required
larval_params	LarvalParams	Ontogenetic interpolation parameters.	required

Returns:

Type	Description
tuple[ndarray, ndarray]	(new_weights, new_energy_reserves).

Source code in pypath/ibm/bioenergetics.py

def growth_step_batch_ontogenetic(
    weights: np.ndarray,
    energy_reserves: np.ndarray,
    consumptions: np.ndarray,
    temperature: "float | np.ndarray",
    is_mature: np.ndarray,
    dt: float,
    bioenerg_params: BioenergParams,
    larval_params: "LarvalParams",
) -> tuple[np.ndarray, np.ndarray]:
    """Vectorized growth step with ontogenetic (size-dependent) interpolation.

    Extends the Wisconsin bioenergetics model with sigmoid-interpolated
    activity multiplier, assimilation efficiency, and basal metabolism split
    (Rs + Ra).  At large (adult) body sizes the results converge to those of
    :func:`growth_step_batch` when ``larval_params.rs_a * (1 + am_max) == ra``.

    Parameters
    ----------
    weights : np.ndarray
        Body weights (grams), shape ``(n,)``.
    energy_reserves : np.ndarray
        Energy reserves, shape ``(n,)``.
    consumptions : np.ndarray
        Total consumption per individual, shape ``(n,)``.
    temperature : float or np.ndarray
        Water temperature (degrees C). Scalar or per-individual array.
    is_mature : np.ndarray
        Boolean array, shape ``(n,)``.
    dt : float
        Timestep (fraction of a year).
    bioenerg_params : BioenergParams
        Standard bioenergetics parameters.
    larval_params : LarvalParams
        Ontogenetic interpolation parameters.

    Returns
    -------
    tuple[np.ndarray, np.ndarray]
        ``(new_weights, new_energy_reserves)``.
    """
    lp = larval_params
    bp = bioenerg_params

    # --- Sigmoid helper (vectorized) ---
    def _sigmoid(x: np.ndarray) -> np.ndarray:
        x = np.clip(x, -30, 30)
        return 1.0 / (1.0 + np.exp(-x))

    # --- Activity multiplier: size-dependent ---
    am = lp.am_min + (lp.am_max - lp.am_min) * _sigmoid(
        (weights - lp.w_activity_mid) / lp.w_activity_scale
    )

    # --- Metabolism: Rs * (1 + am) ---
    q10_factor = q10_temperature_factor(temperature, bp.t_ref, bp.q10)
    # rs_a * w^rb gives per-gram basal rate; multiply by (1 + am) for total met rate
    met = lp.rs_a * (weights**bp.rb) * q10_factor * (1.0 + am) * dt * 365.0

    # --- Assimilation efficiency: size-dependent ---
    ae = lp.ae_min + (lp.ae_max - lp.ae_min) * _sigmoid(
        (weights - lp.w_ae_mid) / lp.w_ae_scale
    )
    assim = consumptions * ae

    # --- SDA ---
    sda = consumptions * bp.sda_fraction

    # --- Net energy ---
    net_energy = assim - met - sda

    # --- Reproduction cost for mature fish with positive surplus ---
    repro_cost = np.where(
        np.asarray(is_mature, dtype=bool) & (net_energy > 0),
        net_energy * bp.reproduction_fraction,
        0.0,
    )
    net_energy = net_energy - repro_cost

    # --- Weight and energy reserve update ---
    weight_change = net_energy / bp.energy_density
    # Use a lower minimum weight (0.0001g) than the adult model (0.1g)
    # to accommodate larvae that naturally weigh < 0.1g.
    new_weights = np.maximum(weights + weight_change, 0.0001)
    new_energy_reserves = np.maximum(energy_reserves + weight_change, 0.0)

    return new_weights, new_energy_reserves

Predation

Size-structured predation mortality using a log-normal selectivity curve. Distributes Ecosim group-level mortality across super-individuals based on body length.

pypath.ibm.predation

Size-structured predation module for the IBM.

Distributes Ecosim group-level predation mortality across IBM super-individuals based on body size using a log-normal selectivity curve. Larger or smaller fish relative to the optimal prey length experience proportionally less predation pressure, reflecting size-dependent vulnerability to predators.

Functions:

Name	Description
size_selectivity	Log-normal selectivity based on prey body length.
distribute_mortality	Allocate group-level mortality across super-individuals.
apply_predation_mortality	Apply predation mortality and return surviving individuals.

Classes:

Name	Description
PredationParams	Dataclass holding size-selectivity parameters.

PredationParams dataclass

Parameters for size-structured predation selectivity.

Parameters:

Name	Type	Description	Default
optimal_prey_length	float	Prey body length (cm) at which predation selectivity is maximised.	required
selectivity_sd	float	Standard deviation of the log-normal selectivity curve (in log-space units). Larger values yield a flatter curve.	required

Source code in pypath/ibm/predation.py

@dataclass
class PredationParams:
    """Parameters for size-structured predation selectivity.

    Parameters
    ----------
    optimal_prey_length : float
        Prey body length (cm) at which predation selectivity is maximised.
    selectivity_sd : float
        Standard deviation of the log-normal selectivity curve (in log-space
        units).  Larger values yield a flatter curve.
    """

    optimal_prey_length: float
    selectivity_sd: float

size_selectivity

size_selectivity(length: float, params: PredationParams) -> float

Compute log-normal size selectivity for a given prey length.

The selectivity peaks at 1.0 when length equals params.optimal_prey_length and decays symmetrically in log-space as the prey deviates from the optimal size.

Parameters:

Name	Type	Description	Default
length	float	Body length of the prey (cm).	required
params	PredationParams	Predation selectivity parameters.	required

Returns:

Type	Description
float	Selectivity value in [0.0, 1.0]. Returns 0.0 if length <= 0.

Source code in pypath/ibm/predation.py

def size_selectivity(length: float, params: PredationParams) -> float:
    """Compute log-normal size selectivity for a given prey length.

    The selectivity peaks at 1.0 when *length* equals
    ``params.optimal_prey_length`` and decays symmetrically in log-space
    as the prey deviates from the optimal size.

    Parameters
    ----------
    length : float
        Body length of the prey (cm).
    params : PredationParams
        Predation selectivity parameters.

    Returns
    -------
    float
        Selectivity value in [0.0, 1.0].  Returns 0.0 if *length* <= 0.
    """
    if length <= 0.0:
        return 0.0
    z = math.log(length / params.optimal_prey_length) / params.selectivity_sd
    return math.exp(-0.5 * z * z)

distribute_mortality

distribute_mortality(individuals: List[SuperIndividual], total_mortality_rate: float, dt: float, params: PredationParams) -> List[float]

Distribute group-level mortality across super-individuals by size.

Deaths are allocated proportionally to each individual's selectivity-weighted abundance (n_represented * selectivity).

Parameters:

Name	Type	Description	Default
individuals	List[SuperIndividual]	Current population of super-individuals.	required
total_mortality_rate	float	Annual mortality rate for the functional group (yr^-1).	required
dt	float	Time-step size (years).	required
params	PredationParams	Size-selectivity parameters.	required

Returns:

Type	Description
List[float]	Number of deaths for each super-individual. Each entry is capped at the individual's n_represented.

Source code in pypath/ibm/predation.py

def distribute_mortality(
    individuals: List[SuperIndividual],
    total_mortality_rate: float,
    dt: float,
    params: PredationParams,
) -> List[float]:
    """Distribute group-level mortality across super-individuals by size.

    Deaths are allocated proportionally to each individual's
    selectivity-weighted abundance (``n_represented * selectivity``).

    Parameters
    ----------
    individuals : List[SuperIndividual]
        Current population of super-individuals.
    total_mortality_rate : float
        Annual mortality rate for the functional group (yr^-1).
    dt : float
        Time-step size (years).
    params : PredationParams
        Size-selectivity parameters.

    Returns
    -------
    List[float]
        Number of deaths for each super-individual.  Each entry is
        capped at the individual's ``n_represented``.
    """
    if not individuals:
        return []

    if params.optimal_prey_length <= 0:
        raise ValueError(
            f"optimal_prey_length must be > 0, got {params.optimal_prey_length}"
        )
    if params.selectivity_sd <= 0:
        raise ValueError(f"selectivity_sd must be > 0, got {params.selectivity_sd}")

    # Extract attributes into NumPy arrays for vectorised computation
    n_repr = np.array([ind.n_represented for ind in individuals])
    lengths = np.array([ind.length for ind in individuals])

    total_n = n_repr.sum()
    total_deaths = total_n * total_mortality_rate * dt

    # Vectorised log-normal selectivity: sel = exp(-0.5 * z^2), 0 for length <= 0
    positive = lengths > 0.0
    z = np.zeros(len(individuals))
    z[positive] = (
        np.log(lengths[positive] / params.optimal_prey_length) / params.selectivity_sd
    )
    sel = np.where(positive, np.exp(-0.5 * z * z), 0.0)

    # Selectivity-weighted abundance
    weighted = n_repr * sel
    total_weighted = weighted.sum()

    if total_weighted == 0.0:
        return [0.0] * len(individuals)

    # Distribute deaths proportionally, cap at n_represented
    deaths_arr = np.minimum(total_deaths * weighted / total_weighted, n_repr)

    return deaths_arr.tolist()

apply_predation_mortality

apply_predation_mortality(individuals: List[SuperIndividual], total_mortality_rate: float, dt: float, params: PredationParams) -> List[SuperIndividual]

Apply predation mortality and return surviving super-individuals.

Creates shallow copies of the input individuals, reduces their n_represented by the allocated deaths, and removes any individual whose n_represented drops to zero or below. The original individuals list is not modified.

Parameters:

Name	Type	Description	Default
individuals	List[SuperIndividual]	Current population of super-individuals (not modified).	required
total_mortality_rate	float	Annual mortality rate for the functional group (yr^-1).	required
dt	float	Time-step size (years).	required
params	PredationParams	Size-selectivity parameters.	required

Returns:

Type	Description
List[SuperIndividual]	Surviving super-individuals with updated n_represented.

Source code in pypath/ibm/predation.py

def apply_predation_mortality(
    individuals: List[SuperIndividual],
    total_mortality_rate: float,
    dt: float,
    params: PredationParams,
) -> List[SuperIndividual]:
    """Apply predation mortality and return surviving super-individuals.

    Creates shallow copies of the input individuals, reduces their
    ``n_represented`` by the allocated deaths, and removes any
    individual whose ``n_represented`` drops to zero or below.
    The original *individuals* list is **not** modified.

    Parameters
    ----------
    individuals : List[SuperIndividual]
        Current population of super-individuals (not modified).
    total_mortality_rate : float
        Annual mortality rate for the functional group (yr^-1).
    dt : float
        Time-step size (years).
    params : PredationParams
        Size-selectivity parameters.

    Returns
    -------
    List[SuperIndividual]
        Surviving super-individuals with updated ``n_represented``.
    """
    if not individuals:
        return []

    deaths = distribute_mortality(individuals, total_mortality_rate, dt, params)

    survivors: List[SuperIndividual] = []
    for ind, d in zip(individuals, deaths):
        survivor = copy.copy(ind)
        survivor.n_represented = ind.n_represented - d
        if survivor.n_represented > 0.0:
            survivors.append(survivor)

    return survivors

Reproduction

Stochastic spawning with weight-dependent fecundity and Cushing match/mismatch larval survival. Mature females produce eggs; surviving larvae become new recruit super-individuals.

pypath.ibm.reproduction

Stochastic reproduction module for the IBM.

Implements spawning, fecundity calculation, and the Cushing match/mismatch hypothesis for larval survival in Baltic smelt. Mature females produce eggs proportional to body weight, and larval survival depends on the temporal overlap between hatching and zooplankton peak abundance.

Functions:

Name	Description
calculate_fecundity	Weight-dependent egg production per female.
larval_survival_probability	Gaussian match/mismatch survival for larvae.
spawn	Determine total egg production for a super-individual.
create_recruits	Create new super-individual recruits from surviving larvae.

Classes:

Name	Description
ReproductionParams	Dataclass holding reproduction and larval survival parameters.

ReproductionParams dataclass

Parameters for stochastic reproduction and larval survival.

Parameters:

Name	Type	Description	Default
fecundity_coefficient	float	Coefficient in the fecundity-weight relationship (eggs = coefficient * weight ^ exponent).	required
fecundity_exponent	float	Exponent in the fecundity-weight relationship.	required
larval_base_survival	float	Base survival probability at perfect match (0-1).	required
zooplankton_match_window	float	Width of the Gaussian match/mismatch window (days).	required
maturity_energy_threshold	float	Minimum energy reserve (kJ) required for spawning.	required
spawning_temp_threshold	float	Minimum water temperature (C) for spawning to occur.	required
larval_duration_days	int	Duration of the larval phase (days).	required
recruit_weight	float	Body weight (g) of a newly recruited individual.	required
recruit_length	float	Body length (cm) of a newly recruited individual.	required

Source code in pypath/ibm/reproduction.py

@dataclass
class ReproductionParams:
    """Parameters for stochastic reproduction and larval survival.

    Parameters
    ----------
    fecundity_coefficient : float
        Coefficient in the fecundity-weight relationship
        (eggs = coefficient * weight ^ exponent).
    fecundity_exponent : float
        Exponent in the fecundity-weight relationship.
    larval_base_survival : float
        Base survival probability at perfect match (0-1).
    zooplankton_match_window : float
        Width of the Gaussian match/mismatch window (days).
    maturity_energy_threshold : float
        Minimum energy reserve (kJ) required for spawning.
    spawning_temp_threshold : float
        Minimum water temperature (C) for spawning to occur.
    larval_duration_days : int
        Duration of the larval phase (days).
    recruit_weight : float
        Body weight (g) of a newly recruited individual.
    recruit_length : float
        Body length (cm) of a newly recruited individual.
    """

    fecundity_coefficient: float
    fecundity_exponent: float
    larval_base_survival: float
    zooplankton_match_window: float
    maturity_energy_threshold: float
    spawning_temp_threshold: float
    larval_duration_days: int
    recruit_weight: float
    recruit_length: float

calculate_fecundity

calculate_fecundity(weight: float, params: ReproductionParams) -> float

Calculate egg production for a single female of a given weight.

Uses a power-law relationship: eggs = coefficient * weight ^ exponent.

Parameters:

Name	Type	Description	Default
weight	float	Individual body weight (g).	required
params	ReproductionParams	Reproduction parameters containing fecundity coefficient and exponent.	required

Returns:

Type	Description
float	Number of eggs produced. Returns 0.0 if weight <= 0.

Source code in pypath/ibm/reproduction.py

def calculate_fecundity(weight: float, params: ReproductionParams) -> float:
    """Calculate egg production for a single female of a given weight.

    Uses a power-law relationship: ``eggs = coefficient * weight ^ exponent``.

    Parameters
    ----------
    weight : float
        Individual body weight (g).
    params : ReproductionParams
        Reproduction parameters containing fecundity coefficient and exponent.

    Returns
    -------
    float
        Number of eggs produced.  Returns 0.0 if *weight* <= 0.
    """
    if weight <= 0.0:
        return 0.0
    return params.fecundity_coefficient * weight**params.fecundity_exponent

larval_survival_probability

larval_survival_probability(spawn_day: float, zoo_peak_day: float, params: ReproductionParams) -> float

Compute larval survival probability using the Cushing match/mismatch hypothesis.

Survival follows a Gaussian function of the temporal mismatch between spawning and the zooplankton peak:

survival = base_survival * exp(-0.5 * (mismatch / match_window)^2)

Parameters:

Name	Type	Description	Default
spawn_day	float	Day of year when spawning occurs.	required
zoo_peak_day	float	Day of year of the zooplankton abundance peak.	required
params	ReproductionParams	Reproduction parameters containing base survival and match window.	required

Returns:

Type	Description
float	Survival probability in (0, base_survival].

Source code in pypath/ibm/reproduction.py

def larval_survival_probability(
    spawn_day: float,
    zoo_peak_day: float,
    params: ReproductionParams,
) -> float:
    """Compute larval survival probability using the Cushing match/mismatch hypothesis.

    Survival follows a Gaussian function of the temporal mismatch between
    spawning and the zooplankton peak:

        survival = base_survival * exp(-0.5 * (mismatch / match_window)^2)

    Parameters
    ----------
    spawn_day : float
        Day of year when spawning occurs.
    zoo_peak_day : float
        Day of year of the zooplankton abundance peak.
    params : ReproductionParams
        Reproduction parameters containing base survival and match window.

    Returns
    -------
    float
        Survival probability in (0, base_survival].
    """
    mismatch = abs(spawn_day - zoo_peak_day)
    z = mismatch / params.zooplankton_match_window
    return params.larval_base_survival * math.exp(-0.5 * z * z)

spawn

spawn(individual: SuperIndividual, temperature: float, params: ReproductionParams) -> float

Determine total egg production for a super-individual.

Only mature females (is_mature=True, sex=0) spawn, and only when temperature and energy reserves meet the required thresholds.

Parameters:

Name	Type	Description	Default
individual	SuperIndividual	The super-individual attempting to spawn.	required
temperature	float	Current water temperature (C).	required
params	ReproductionParams	Reproduction parameters.	required

Returns:

Type	Description
float	Total number of eggs produced (n_represented * fecundity_per_female). Returns 0.0 if spawning conditions are not met.

Source code in pypath/ibm/reproduction.py

def spawn(
    individual: SuperIndividual,
    temperature: float,
    params: ReproductionParams,
) -> float:
    """Determine total egg production for a super-individual.

    Only mature females (``is_mature=True``, ``sex=0``) spawn, and only
    when temperature and energy reserves meet the required thresholds.

    Parameters
    ----------
    individual : SuperIndividual
        The super-individual attempting to spawn.
    temperature : float
        Current water temperature (C).
    params : ReproductionParams
        Reproduction parameters.

    Returns
    -------
    float
        Total number of eggs produced (n_represented * fecundity_per_female).
        Returns 0.0 if spawning conditions are not met.
    """
    # Only mature females spawn
    if not individual.is_mature:
        return 0.0
    if individual.sex != 0:
        return 0.0

    # Check environmental and physiological thresholds
    if temperature < params.spawning_temp_threshold:
        return 0.0
    if individual.energy_reserve < params.maturity_energy_threshold:
        return 0.0

    fecundity_per_female = calculate_fecundity(individual.weight, params)
    return individual.n_represented * fecundity_per_female

create_recruits

create_recruits(total_eggs: float, spawn_day: float, zoo_peak_day: float, patch_idx: int, next_id: int, params: ReproductionParams, n_super_individuals: int = 1) -> List[SuperIndividual]

Create new super-individual recruits from surviving larvae.

Calculates the number of surviving larvae using :func:larval_survival_probability, then distributes the survivors evenly across n_super_individuals new :class:SuperIndividual objects.

Parameters:

Name	Type	Description	Default
total_eggs	float	Total number of eggs produced.	required
spawn_day	float	Day of year when spawning occurred.	required
zoo_peak_day	float	Day of year of the zooplankton abundance peak.	required
patch_idx	int	Spatial patch index where recruits are placed.	required
next_id	int	Starting ID for the new super-individuals.	required
params	ReproductionParams	Reproduction parameters.	required
n_super_individuals	int	Number of super-individuals to create (default 1).	1

Returns:

Type	Description
List[SuperIndividual]	New recruit super-individuals. Empty list if total survivors < 1.

Source code in pypath/ibm/reproduction.py

def create_recruits(
    total_eggs: float,
    spawn_day: float,
    zoo_peak_day: float,
    patch_idx: int,
    next_id: int,
    params: ReproductionParams,
    n_super_individuals: int = 1,
) -> List[SuperIndividual]:
    """Create new super-individual recruits from surviving larvae.

    Calculates the number of surviving larvae using
    :func:`larval_survival_probability`, then distributes the survivors
    evenly across *n_super_individuals* new :class:`SuperIndividual`
    objects.

    Parameters
    ----------
    total_eggs : float
        Total number of eggs produced.
    spawn_day : float
        Day of year when spawning occurred.
    zoo_peak_day : float
        Day of year of the zooplankton abundance peak.
    patch_idx : int
        Spatial patch index where recruits are placed.
    next_id : int
        Starting ID for the new super-individuals.
    params : ReproductionParams
        Reproduction parameters.
    n_super_individuals : int, optional
        Number of super-individuals to create (default 1).

    Returns
    -------
    List[SuperIndividual]
        New recruit super-individuals.  Empty list if total survivors < 1.
    """
    survival = larval_survival_probability(spawn_day, zoo_peak_day, params)
    n_survivors = total_eggs * survival

    if n_survivors < 1.0:
        return []

    n_per_si = n_survivors / n_super_individuals

    recruits: List[SuperIndividual] = []
    for i in range(n_super_individuals):
        recruit = SuperIndividual(
            id=next_id + i,
            n_represented=n_per_si,
            weight=params.recruit_weight,
            length=params.recruit_length,
            age=0.0,
            energy_reserve=params.recruit_weight * 0.1,
            patch_idx=patch_idx,
            is_mature=False,
            sex=random.choice([0, 1]),
        )
        recruits.append(recruit)

    return recruits

Behavior (Spatial Movement & Adaptive Foraging)

Score-based stochastic patch selection over a sparse adjacency graph, plus profitability-based adaptive prey allocation.

pypath.ibm.behavior

Behavior module for the IBM: spatial movement and adaptive foraging.

Handles two key behaviors for IBM super-individuals:

1. Spatial movement between ECOSPACE patches, using a weighted score of habitat quality, food density, and predator avoidance to compute movement probabilities across a sparse adjacency graph.

2. Adaptive foraging (prey selection), where super-individuals allocate consumption proportionally to prey profitability while respecting availability constraints.

Functions:

Name	Description
calculate_movement_probabilities	Compute per-patch movement probabilities for a super-individual.
move_individual	Move a super-individual to a new patch based on movement probabilities.
should_migrate	Determine whether migration conditions are met.
adaptive_forage	Allocate consumption across prey groups by profitability.

Classes:

Name	Description
MovementParams	Dataclass holding spatial movement parameters.
ForagingParams	Dataclass holding adaptive foraging parameters.

MovementParams dataclass

Parameters controlling spatial movement between patches.

Parameters:

Name	Type	Description	Default
base_speed	float	Base movement probability scaling (0-1). Low values increase the inertia bonus, making the individual more likely to stay.	required
habitat_weight	float	Weight for habitat quality in the movement score (0-1).	required
food_weight	float	Weight for food density in the movement score (0-1).	required
predator_weight	float	Weight for predator avoidance in the movement score (0-1).	required
migration_temp_threshold	float	Temperature (degrees C) above which spring migration can occur.	required
migration_months	tuple	Months (1-12) during which migration can occur.	(3, 4, 5)

Source code in pypath/ibm/behavior.py

@dataclass
class MovementParams:
    """Parameters controlling spatial movement between patches.

    Parameters
    ----------
    base_speed : float
        Base movement probability scaling (0-1).  Low values increase
        the inertia bonus, making the individual more likely to stay.
    habitat_weight : float
        Weight for habitat quality in the movement score (0-1).
    food_weight : float
        Weight for food density in the movement score (0-1).
    predator_weight : float
        Weight for predator avoidance in the movement score (0-1).
    migration_temp_threshold : float
        Temperature (degrees C) above which spring migration can occur.
    migration_months : tuple
        Months (1-12) during which migration can occur.
    """

    base_speed: float
    habitat_weight: float
    food_weight: float
    predator_weight: float
    migration_temp_threshold: float
    migration_months: Tuple[int, ...] = (3, 4, 5)

ForagingParams dataclass

Parameters controlling adaptive prey selection.

Parameters:

Name	Type	Description	Default
energy_content	ndarray	Energy per gram of each prey group (kJ/g). Shape (n_prey,).	required
handling_time	ndarray	Handling time per gram of each prey group. Shape (n_prey,).	required

Source code in pypath/ibm/behavior.py

@dataclass
class ForagingParams:
    """Parameters controlling adaptive prey selection.

    Parameters
    ----------
    energy_content : np.ndarray
        Energy per gram of each prey group (kJ/g).  Shape ``(n_prey,)``.
    handling_time : np.ndarray
        Handling time per gram of each prey group.  Shape ``(n_prey,)``.
    """

    energy_content: np.ndarray
    handling_time: np.ndarray

calculate_movement_probabilities

calculate_movement_probabilities(current_patch: int, adjacency: csr_matrix, habitat_quality: ndarray, food_density: ndarray, predator_density: ndarray, params: MovementParams) -> np.ndarray

Compute per-patch movement probabilities for a super-individual.

For each reachable patch (current patch plus its neighbors in the sparse adjacency matrix), a weighted score is computed as:

score = habitat_weight * habitat_quality
      + food_weight * food_density
      + predator_weight * (1 / (1 + predator_density))

The current patch receives an inertia bonus proportional to (1 - base_speed), making the individual more sedentary when base_speed is low.

Scores are normalized to probabilities summing to 1.0. If all scores are zero, the individual stays in its current patch with probability 1.0.

Parameters:

Name	Type	Description	Default
current_patch	int	Index of the patch currently occupied.	required
adjacency	csr_matrix	Sparse adjacency matrix (n_patches x n_patches).	required
habitat_quality	ndarray	Per-patch habitat quality, shape (n_patches,).	required
food_density	ndarray	Per-patch food density, shape (n_patches,).	required
predator_density	ndarray	Per-patch predator density, shape (n_patches,).	required
params	MovementParams	Movement parameters.	required

Returns:

Type	Description
ndarray	Probability array of shape (n_patches,) summing to 1.0.

Source code in pypath/ibm/behavior.py

def calculate_movement_probabilities(
    current_patch: int,
    adjacency: sp.csr_matrix,
    habitat_quality: np.ndarray,
    food_density: np.ndarray,
    predator_density: np.ndarray,
    params: MovementParams,
) -> np.ndarray:
    """Compute per-patch movement probabilities for a super-individual.

    For each reachable patch (current patch plus its neighbors in the
    sparse adjacency matrix), a weighted score is computed as:

        score = habitat_weight * habitat_quality
              + food_weight * food_density
              + predator_weight * (1 / (1 + predator_density))

    The current patch receives an inertia bonus proportional to
    ``(1 - base_speed)``, making the individual more sedentary when
    ``base_speed`` is low.

    Scores are normalized to probabilities summing to 1.0.  If all
    scores are zero, the individual stays in its current patch with
    probability 1.0.

    Parameters
    ----------
    current_patch : int
        Index of the patch currently occupied.
    adjacency : scipy.sparse.csr_matrix
        Sparse adjacency matrix (n_patches x n_patches).
    habitat_quality : np.ndarray
        Per-patch habitat quality, shape ``(n_patches,)``.
    food_density : np.ndarray
        Per-patch food density, shape ``(n_patches,)``.
    predator_density : np.ndarray
        Per-patch predator density, shape ``(n_patches,)``.
    params : MovementParams
        Movement parameters.

    Returns
    -------
    np.ndarray
        Probability array of shape ``(n_patches,)`` summing to 1.0.
    """
    n_patches = adjacency.shape[0]
    probs = np.zeros(n_patches, dtype=np.float64)

    # Identify reachable patches: self + neighbors
    row = adjacency.getrow(current_patch)
    neighbor_indices = row.indices.tolist()
    reachable = set(neighbor_indices)
    reachable.add(current_patch)

    # Compute scores for reachable patches
    for p in reachable:
        pred_avoidance = 1.0 / (1.0 + predator_density[p])
        score = (
            params.habitat_weight * habitat_quality[p]
            + params.food_weight * food_density[p]
            + params.predator_weight * pred_avoidance
        )
        probs[p] = score

    total = probs.sum()
    if total == 0.0:
        # All zero: stay in current patch
        probs[current_patch] = 1.0
        return probs

    probs /= total

    # Apply additive inertia bonus after normalization so zero-quality patches
    # still retain some probability of staying.  Then renormalize.
    inertia = 1.0 - params.base_speed
    if inertia > 0.0:
        probs[current_patch] = max(probs[current_patch], 0.0) + inertia
        probs /= probs.sum()

    return probs

move_individual

move_individual(individual: SuperIndividual, adjacency: csr_matrix, habitat_quality: ndarray, food_density: ndarray, predator_density: ndarray, params: MovementParams, rng: Optional[Generator] = None) -> SuperIndividual

Move a super-individual to a new patch based on movement probabilities.

Computes movement probabilities via :func:calculate_movement_probabilities and uses np.random.choice (or the supplied RNG) to select a destination patch. Returns a copy of the individual with updated patch_idx; the original is not modified.

Parameters:

Name	Type	Description	Default
individual	SuperIndividual	The super-individual to move.	required
adjacency	csr_matrix	Sparse adjacency matrix.	required
habitat_quality	ndarray	Per-patch habitat quality.	required
food_density	ndarray	Per-patch food density.	required
predator_density	ndarray	Per-patch predator density.	required
params	MovementParams	Movement parameters.	required
rng	Generator	Random number generator. If None, a default is created.	None

Returns:

Type	Description
SuperIndividual	A copy with potentially updated patch_idx.

Source code in pypath/ibm/behavior.py

def move_individual(
    individual: SuperIndividual,
    adjacency: sp.csr_matrix,
    habitat_quality: np.ndarray,
    food_density: np.ndarray,
    predator_density: np.ndarray,
    params: MovementParams,
    rng: Optional[np.random.Generator] = None,
) -> SuperIndividual:
    """Move a super-individual to a new patch based on movement probabilities.

    Computes movement probabilities via :func:`calculate_movement_probabilities`
    and uses ``np.random.choice`` (or the supplied RNG) to select a destination
    patch.  Returns a **copy** of the individual with updated ``patch_idx``; the
    original is not modified.

    Parameters
    ----------
    individual : SuperIndividual
        The super-individual to move.
    adjacency : scipy.sparse.csr_matrix
        Sparse adjacency matrix.
    habitat_quality : np.ndarray
        Per-patch habitat quality.
    food_density : np.ndarray
        Per-patch food density.
    predator_density : np.ndarray
        Per-patch predator density.
    params : MovementParams
        Movement parameters.
    rng : np.random.Generator, optional
        Random number generator.  If ``None``, a default is created.

    Returns
    -------
    SuperIndividual
        A copy with potentially updated ``patch_idx``.
    """
    if rng is None:
        rng = np.random.default_rng()

    probs = calculate_movement_probabilities(
        current_patch=individual.patch_idx,
        adjacency=adjacency,
        habitat_quality=habitat_quality,
        food_density=food_density,
        predator_density=predator_density,
        params=params,
    )

    n_patches = len(probs)
    new_patch = rng.choice(n_patches, p=probs)

    result = copy.copy(individual)
    result.patch_idx = int(new_patch)
    return result

should_migrate

should_migrate(temperature: float, month: int, params: MovementParams) -> bool

Determine whether migration conditions are met.

Migration occurs when the temperature is strictly above the threshold and the current month is in the configured migration months.

Parameters:

Name	Type	Description	Default
temperature	float	Current water temperature (degrees C).	required
month	int	Current month (1-12).	required
params	MovementParams	Movement parameters containing threshold and migration months.	required

Returns:

Type	Description
bool	True if migration conditions are met.

Source code in pypath/ibm/behavior.py

def should_migrate(
    temperature: float,
    month: int,
    params: MovementParams,
) -> bool:
    """Determine whether migration conditions are met.

    Migration occurs when the temperature is **strictly above** the
    threshold and the current month is in the configured migration months.

    Parameters
    ----------
    temperature : float
        Current water temperature (degrees C).
    month : int
        Current month (1-12).
    params : MovementParams
        Movement parameters containing threshold and migration months.

    Returns
    -------
    bool
        True if migration conditions are met.
    """
    return (
        temperature > params.migration_temp_threshold
        and month in params.migration_months
    )

adaptive_forage

adaptive_forage(prey_available: Dict[int, float], max_consumption: float, individual_length: float, params: ForagingParams) -> Dict[int, float]

Allocate consumption across prey groups by profitability.

Profitability of each prey group is defined as:

profitability = (energy_content / handling_time) * availability

Consumption is allocated proportionally to profitability, subject to availability constraints (cannot eat more than available for any group). When a group is availability-constrained, the surplus is redistributed to the remaining groups in a second pass.

Parameters:

Name	Type	Description	Default
prey_available	Dict[int, float]	Available biomass per prey group index.	required
max_consumption	float	Maximum total consumption for this individual.	required
individual_length	float	Body length of the individual (cm). Currently unused but reserved for future size-dependent foraging selectivity.	required
params	ForagingParams	Foraging parameters.	required

Returns:

Type	Description
Dict[int, float]	Consumption by prey group index.

Source code in pypath/ibm/behavior.py

def adaptive_forage(
    prey_available: Dict[int, float],
    max_consumption: float,
    individual_length: float,
    params: ForagingParams,
) -> Dict[int, float]:
    """Allocate consumption across prey groups by profitability.

    Profitability of each prey group is defined as:

        profitability = (energy_content / handling_time) * availability

    Consumption is allocated proportionally to profitability, subject to
    availability constraints (cannot eat more than available for any group).
    When a group is availability-constrained, the surplus is redistributed
    to the remaining groups in a second pass.

    Parameters
    ----------
    prey_available : Dict[int, float]
        Available biomass per prey group index.
    max_consumption : float
        Maximum total consumption for this individual.
    individual_length : float
        Body length of the individual (cm).  Currently unused but
        reserved for future size-dependent foraging selectivity.
    params : ForagingParams
        Foraging parameters.

    Returns
    -------
    Dict[int, float]
        Consumption by prey group index.
    """
    if not prey_available or max_consumption <= 0.0:
        return {k: 0.0 for k in prey_available}

    # Compute profitability for each available prey group
    profitabilities: Dict[int, float] = {}
    for group_idx, available in prey_available.items():
        if available <= 0.0:
            profitabilities[group_idx] = 0.0
            continue
        ec = params.energy_content[group_idx]
        ht = params.handling_time[group_idx]
        if ht <= 0.0:
            profitabilities[group_idx] = 0.0
            continue
        profitabilities[group_idx] = (ec / ht) * available

    total_profitability = sum(profitabilities.values())

    if total_profitability <= 0.0:
        return {k: 0.0 for k in prey_available}

    # Iteratively allocate consumption, respecting availability limits.
    # When a group hits its cap, redistribute the surplus.
    consumption: Dict[int, float] = {k: 0.0 for k in prey_available}
    remaining_consumption = max_consumption
    active_groups = set(prey_available.keys())

    for _ in range(len(prey_available) + 1):
        if remaining_consumption <= 0.0 or not active_groups:
            break

        # Profitability among active groups
        active_prof = {g: profitabilities[g] for g in active_groups}
        total_active = sum(active_prof.values())
        if total_active <= 0.0:
            break

        # Proportional allocation
        capped = False
        for g in list(active_groups):
            share = remaining_consumption * active_prof[g] / total_active
            available = prey_available[g] - consumption[g]
            if share >= available:
                # Capped by availability
                consumption[g] += available
                remaining_consumption -= available
                active_groups.discard(g)
                capped = True
            else:
                consumption[g] += share

        if not capped:
            # All allocations fit within availability
            remaining_consumption = 0.0
            break

    return consumption

EwE Integration

Bridge between the IBM and Ecosim/Ecospace derivative systems. apply_ibm_to_derivative() overrides the Ecosim derivative for IBM groups and subtracts IBM consumption from prey derivatives.

pypath.ibm.integration

IBM-Ecosim derivative override integration.

Provides the bridge functions that allow IBM-managed functional groups to override the standard Ecosim derivative calculation. When an IBM group is active, these functions extract the relevant food-web data from the Ecosim consumption matrix, delegate the dynamics to the IBM engine, and write the results back into the derivative vector.

Functions:

Name	Description
extract_prey_availability	Extract non-zero prey consumption rates for a given predator.
extract_predation_pressure	Sum total predation on a given prey from all living predators.
check_ibm_mass_balance	Validate that an IBM step result is physically plausible.
apply_ibm_to_derivative	Override the Ecosim derivative for an IBM group in-place.

extract_prey_availability

extract_prey_availability(QQ: ndarray, predator_idx: int, n_groups: int) -> Dict[int, float]

Extract non-zero prey consumption rates for a predator.

Reads column predator_idx of the consumption matrix QQ and returns a dictionary mapping each prey index to its consumption rate, excluding entries that are exactly zero.

Parameters:

Name	Type	Description	Default
QQ	ndarray	Consumption matrix of shape (n_groups+1, n_groups+1) where QQ[prey, predator] is the consumption rate.	required
predator_idx	int	1-based index of the predator group.	required
n_groups	int	Total number of functional groups (excluding the 0-index padding).	required

Returns:

Type	Description
Dict[int, float]	Mapping {prey_idx: consumption_rate} for all non-zero entries.

Source code in pypath/ibm/integration.py

def extract_prey_availability(
    QQ: np.ndarray, predator_idx: int, n_groups: int
) -> Dict[int, float]:
    """Extract non-zero prey consumption rates for a predator.

    Reads column ``predator_idx`` of the consumption matrix ``QQ`` and
    returns a dictionary mapping each prey index to its consumption rate,
    excluding entries that are exactly zero.

    Parameters
    ----------
    QQ : np.ndarray
        Consumption matrix of shape ``(n_groups+1, n_groups+1)`` where
        ``QQ[prey, predator]`` is the consumption rate.
    predator_idx : int
        1-based index of the predator group.
    n_groups : int
        Total number of functional groups (excluding the 0-index padding).

    Returns
    -------
    Dict[int, float]
        Mapping ``{prey_idx: consumption_rate}`` for all non-zero entries.
    """
    result: Dict[int, float] = {}
    for prey in range(1, n_groups + 1):
        rate = float(QQ[prey, predator_idx])
        if rate != 0.0:
            result[prey] = rate
    return result

extract_predation_pressure

extract_predation_pressure(QQ: ndarray, prey_idx: int, n_living: int) -> float

Sum total predation on a prey group from all living predators.

Parameters:

Name	Type	Description	Default
QQ	ndarray	Consumption matrix of shape (n_groups+1, n_groups+1).	required
prey_idx	int	1-based index of the prey group.	required
n_living	int	Number of living groups in the model.	required

Returns:

Type	Description
float	Sum of QQ[prey_idx, 1:n_living+1].

Source code in pypath/ibm/integration.py

def extract_predation_pressure(QQ: np.ndarray, prey_idx: int, n_living: int) -> float:
    """Sum total predation on a prey group from all living predators.

    Parameters
    ----------
    QQ : np.ndarray
        Consumption matrix of shape ``(n_groups+1, n_groups+1)``.
    prey_idx : int
        1-based index of the prey group.
    n_living : int
        Number of living groups in the model.

    Returns
    -------
    float
        Sum of ``QQ[prey_idx, 1:n_living+1]``.
    """
    return float(np.sum(QQ[prey_idx, 1 : n_living + 1]))

check_ibm_mass_balance

check_ibm_mass_balance(result: 'IBMStepResult', tolerance: float = 0.05) -> Tuple[bool, float]

Validate that an IBM step result is physically plausible.

Checks that biomass is non-negative and that no consumption entry is negative. The relative error returned is the magnitude of the largest violation (or 0.0 if none).

Parameters:

Name	Type	Description	Default
result	IBMStepResult	The result from an IBM integration step.	required
tolerance	float	Maximum acceptable relative error (default 0.05).	0.05

Returns:

Type	Description
Tuple[bool, float]	(is_balanced, relative_error) where is_balanced is True when all checks pass within the tolerance.

Source code in pypath/ibm/integration.py

def check_ibm_mass_balance(
    result: "IBMStepResult", tolerance: float = 0.05
) -> Tuple[bool, float]:
    """Validate that an IBM step result is physically plausible.

    Checks that biomass is non-negative and that no consumption entry
    is negative. The relative error returned is the magnitude of the
    largest violation (or 0.0 if none).

    Parameters
    ----------
    result : IBMStepResult
        The result from an IBM integration step.
    tolerance : float, optional
        Maximum acceptable relative error (default 0.05).

    Returns
    -------
    Tuple[bool, float]
        ``(is_balanced, relative_error)`` where ``is_balanced`` is True
        when all checks pass within the tolerance.
    """
    error = 0.0

    if result.biomass < 0.0:
        error = max(error, abs(result.biomass))
        return False, error

    if np.any(result.consumption_by_prey < 0.0):
        min_consumption = float(np.min(result.consumption_by_prey))
        error = max(error, abs(min_consumption))
        return False, error

    return True, error

apply_ibm_to_derivative

apply_ibm_to_derivative(deriv: ndarray, QQ: ndarray, BB: ndarray, ibm_group: 'IBMGroup', forcing: dict, dt: float, spatial_context: Optional['SpatialContext'] = None) -> None

Override the Ecosim derivative for an IBM-managed group in-place.

Calls the IBM group's compute_step method with prey availability and predation pressure extracted from the current Ecosim state, then writes the resulting biomass change into the derivative vector and subtracts IBM consumption from prey derivatives.

Parameters:

Name	Type	Description	Default
deriv	ndarray	Derivative vector (modified in-place). Shape (n_groups+1,).	required
QQ	ndarray	Consumption matrix (n_groups+1, n_groups+1).	required
BB	ndarray	Current biomass vector (n_groups+1,).	required
ibm_group	IBMGroup	The IBM group instance that will compute the step.	required
forcing	dict	Environmental forcing dictionary passed through to compute_step.	required
dt	float	Time step size in years.	required
spatial_context	SpatialContext or None	Spatial context for multi-patch IBM movement. When provided, it is forwarded to ibm_group.compute_step() so the IBM engine can distribute individuals across patches. Default is None.	None

Source code in pypath/ibm/integration.py

def apply_ibm_to_derivative(
    deriv: np.ndarray,
    QQ: np.ndarray,
    BB: np.ndarray,
    ibm_group: "IBMGroup",
    forcing: dict,
    dt: float,
    spatial_context: Optional["SpatialContext"] = None,
) -> None:
    """Override the Ecosim derivative for an IBM-managed group in-place.

    Calls the IBM group's ``compute_step`` method with prey availability
    and predation pressure extracted from the current Ecosim state, then
    writes the resulting biomass change into the derivative vector and
    subtracts IBM consumption from prey derivatives.

    Parameters
    ----------
    deriv : np.ndarray
        Derivative vector (modified in-place). Shape ``(n_groups+1,)``.
    QQ : np.ndarray
        Consumption matrix ``(n_groups+1, n_groups+1)``.
    BB : np.ndarray
        Current biomass vector ``(n_groups+1,)``.
    ibm_group : IBMGroup
        The IBM group instance that will compute the step.
    forcing : dict
        Environmental forcing dictionary passed through to ``compute_step``.
    dt : float
        Time step size in years.
    spatial_context : SpatialContext or None, optional
        Spatial context for multi-patch IBM movement. When provided, it is
        forwarded to ``ibm_group.compute_step()`` so the IBM engine can
        distribute individuals across patches. Default is ``None``.
    """
    group_idx = ibm_group.group_index
    n_groups = ibm_group.n_groups

    # Extract prey availability as a dict, then convert to array.
    # extract_prey_availability returns 1-based keys; convert to 0-based indices.
    prey_dict = extract_prey_availability(QQ, group_idx, n_groups)
    prey_array = np.zeros(n_groups)
    for prey_idx, rate in prey_dict.items():
        zero_idx = prey_idx - 1
        if 0 <= zero_idx < n_groups:
            prey_array[zero_idx] = rate

    # Extract predation pressure from all living predators
    # Use n_groups as upper bound for n_living (safe default)
    predation = extract_predation_pressure(QQ, group_idx, n_groups)

    # Delegate dynamics to the IBM engine
    result = ibm_group.compute_step(
        prey_available=prey_array,
        predation_pressure=predation,
        env_forcing=forcing,
        dt=dt,
        spatial_context=spatial_context,
    )

    # Override derivative: dB/dt = (new_biomass - current_biomass) / dt
    deriv[group_idx] = (result.biomass - BB[group_idx]) / dt

    # Subtract IBM consumption from prey derivatives.
    # consumption_by_prey is 0-based; deriv is 1-based (index 0 = "Outside").
    for prey_idx in range(len(result.consumption_by_prey)):
        consumed = result.consumption_by_prey[prey_idx]
        if consumed != 0.0 and (prey_idx + 1) < len(deriv):
            deriv[prey_idx + 1] -= consumed / dt

Smelt Implementation

Concrete IBMGroup implementation for Baltic smelt (Osmerus eperlanus), combining all sub-modules into a 5-phase compute_step(). Use as a template for building custom IBM species.

pypath.ibm.smelt

SmeltIBM concrete implementation for Baltic smelt (Osmerus eperlanus).

Orchestrates all IBM behavior modules -- bioenergetics, predation, foraging, reproduction, and growth -- into a single cohesive IBM group that can be injected into the Ecosim derivative loop.

The SmeltIBM is initialized from Ecopath equilibrium biomass and creates an age-structured population of super-individuals using Von Bertalanffy growth curves. Each time step, it runs up to five phases:

1. Forage + Grow: adaptive foraging followed by Wisconsin bioenergetics.
2. Reproduce: mature females spawn; surviving larvae become recruits.
3. Predation mortality: size-structured mortality from Ecosim predators.
4. Bookkeeping: add recruits, age individuals, remove senescent fish.
5. Spatial movement (optional): move individuals between patches.

Classes:

Name	Description
SmeltParams	Composite parameter dataclass combining all IBM sub-module parameters.
SmeltIBM	Concrete IBMGroup implementation for Baltic smelt.

SmeltParams dataclass

Composite parameters for the SmeltIBM.

Combines all sub-module parameter sets plus species-specific Von Bertalanffy growth parameters for Baltic smelt.

Parameters:

Name	Type	Description	Default
bioenerg	BioenergParams	Wisconsin bioenergetics model parameters.	required
predation	PredationParams	Size-structured predation selectivity parameters.	required
foraging	ForagingParams	Adaptive prey-selection parameters.	required
movement	MovementParams	Spatial movement parameters.	required
reproduction	ReproductionParams	Spawning and larval survival parameters.	required
vbgf_k_mean	float	Mean Von Bertalanffy growth coefficient K (yr^-1).	0.3
vbgf_k_sd	float	Standard deviation of K across individuals.	0.05
vbgf_linf_mean	float	Mean asymptotic body length Linf (cm).	25.0
vbgf_linf_sd	float	Standard deviation of Linf across individuals.	3.0
max_age	float	Maximum age (years) before natural senescence removal.	10.0

Source code in pypath/ibm/smelt.py

@dataclass
class SmeltParams:
    """Composite parameters for the SmeltIBM.

    Combines all sub-module parameter sets plus species-specific Von
    Bertalanffy growth parameters for Baltic smelt.

    Parameters
    ----------
    bioenerg : BioenergParams
        Wisconsin bioenergetics model parameters.
    predation : PredationParams
        Size-structured predation selectivity parameters.
    foraging : ForagingParams
        Adaptive prey-selection parameters.
    movement : MovementParams
        Spatial movement parameters.
    reproduction : ReproductionParams
        Spawning and larval survival parameters.
    vbgf_k_mean : float
        Mean Von Bertalanffy growth coefficient K (yr^-1).
    vbgf_k_sd : float
        Standard deviation of K across individuals.
    vbgf_linf_mean : float
        Mean asymptotic body length Linf (cm).
    vbgf_linf_sd : float
        Standard deviation of Linf across individuals.
    max_age : float
        Maximum age (years) before natural senescence removal.
    """

    bioenerg: BioenergParams
    predation: PredationParams
    foraging: ForagingParams
    movement: MovementParams
    reproduction: ReproductionParams
    vbgf_k_mean: float = 0.3
    vbgf_k_sd: float = 0.05
    vbgf_linf_mean: float = 25.0
    vbgf_linf_sd: float = 3.0
    max_age: float = 10.0
    max_super_individuals: int = 2000
    egg: Optional[EggParams] = None
    yolk_sac: Optional[YolkSacParams] = None
    larval: Optional[LarvalParams] = None
    oxygen: Optional[OxygenParams] = None
    zones: Optional[ZoneParams] = None

    @classmethod
    def baltic_defaults(cls) -> "SmeltParams":
        """Return default parameters for Baltic smelt (Osmerus eperlanus).

        These are literature-based defaults suitable for the Baltic Sea
        ecosystem.  All sub-parameter sets are populated with species-
        specific values.

        Returns
        -------
        SmeltParams
            Fully populated parameter object.
        """
        # Number of functional groups for foraging arrays (sized generously;
        # actual model may be smaller, but arrays are indexed by group index).
        n_groups_default = 20

        bioenerg = BioenergParams(
            ra=0.0033,  # g O2/g/day at reference temp
            rb=-0.227,  # metabolic weight exponent
            q10=2.1,  # Q10 temperature coefficient
            t_ref=10.0,  # reference temperature (C)
            sda_fraction=0.172,  # specific dynamic action
            unassimilated_fraction=0.27,  # unassimilated fraction
            a_length=0.55,  # allometric length coefficient
            b_length=0.333,  # allometric length exponent (cube root)
            energy_density=5.0,  # kJ/g tissue
            reproduction_fraction=0.3,
        )

        predation = PredationParams(
            optimal_prey_length=10.0,  # cm
            selectivity_sd=0.5,
        )

        foraging = ForagingParams(
            energy_content=np.full(n_groups_default, 4.0),  # kJ/g
            handling_time=np.full(n_groups_default, 1.0),  # time/g
        )

        movement = MovementParams(
            base_speed=0.3,
            habitat_weight=0.4,
            food_weight=0.4,
            predator_weight=0.2,
            migration_temp_threshold=4.0,  # C
            migration_months=(3, 4, 5),
        )

        reproduction = ReproductionParams(
            fecundity_coefficient=200.0,  # eggs per g^exponent
            fecundity_exponent=1.2,
            larval_base_survival=0.01,
            zooplankton_match_window=15.0,  # days
            maturity_energy_threshold=0.5,
            spawning_temp_threshold=4.0,  # C
            larval_duration_days=30,
            recruit_weight=0.5,  # g
            recruit_length=3.0,  # cm
        )

        return cls(
            bioenerg=bioenerg,
            predation=predation,
            foraging=foraging,
            movement=movement,
            reproduction=reproduction,
            vbgf_k_mean=0.3,
            vbgf_k_sd=0.05,
            vbgf_linf_mean=25.0,
            vbgf_linf_sd=3.0,
            max_age=10.0,
        )

    @classmethod
    def baltic_defaults_els(cls) -> "SmeltParams":
        """Return Baltic smelt defaults with early life stages enabled.

        Creates a standard ``baltic_defaults()`` parameter set and adds
        default EggParams, YolkSacParams, LarvalParams, and OxygenParams
        for mechanistic egg-to-larva modeling.

        Returns
        -------
        SmeltParams
            Fully populated parameter object with ELS enabled.
        """
        base = cls.baltic_defaults()
        base.egg = EggParams()
        base.yolk_sac = YolkSacParams()
        base.larval = LarvalParams()
        base.oxygen = OxygenParams()
        return base

    @classmethod
    def baltic_defaults_zonal(cls) -> "SmeltParams":
        """Return Baltic smelt defaults with ELS and zonal model enabled.

        Creates a ``baltic_defaults_els()`` parameter set and adds default
        ``ZoneParams`` for the 3-zone Curonian Lagoon spatial model.

        Returns
        -------
        SmeltParams
            Fully populated parameter object with ELS + zones enabled.
        """
        base = cls.baltic_defaults_els()
        base.zones = ZoneParams()
        return base

baltic_defaults classmethod

baltic_defaults() -> 'SmeltParams'

Return default parameters for Baltic smelt (Osmerus eperlanus).

These are literature-based defaults suitable for the Baltic Sea ecosystem. All sub-parameter sets are populated with species- specific values.

Returns:

Type	Description
SmeltParams	Fully populated parameter object.

Source code in pypath/ibm/smelt.py

@classmethod
def baltic_defaults(cls) -> "SmeltParams":
    """Return default parameters for Baltic smelt (Osmerus eperlanus).

    These are literature-based defaults suitable for the Baltic Sea
    ecosystem.  All sub-parameter sets are populated with species-
    specific values.

    Returns
    -------
    SmeltParams
        Fully populated parameter object.
    """
    # Number of functional groups for foraging arrays (sized generously;
    # actual model may be smaller, but arrays are indexed by group index).
    n_groups_default = 20

    bioenerg = BioenergParams(
        ra=0.0033,  # g O2/g/day at reference temp
        rb=-0.227,  # metabolic weight exponent
        q10=2.1,  # Q10 temperature coefficient
        t_ref=10.0,  # reference temperature (C)
        sda_fraction=0.172,  # specific dynamic action
        unassimilated_fraction=0.27,  # unassimilated fraction
        a_length=0.55,  # allometric length coefficient
        b_length=0.333,  # allometric length exponent (cube root)
        energy_density=5.0,  # kJ/g tissue
        reproduction_fraction=0.3,
    )

    predation = PredationParams(
        optimal_prey_length=10.0,  # cm
        selectivity_sd=0.5,
    )

    foraging = ForagingParams(
        energy_content=np.full(n_groups_default, 4.0),  # kJ/g
        handling_time=np.full(n_groups_default, 1.0),  # time/g
    )

    movement = MovementParams(
        base_speed=0.3,
        habitat_weight=0.4,
        food_weight=0.4,
        predator_weight=0.2,
        migration_temp_threshold=4.0,  # C
        migration_months=(3, 4, 5),
    )

    reproduction = ReproductionParams(
        fecundity_coefficient=200.0,  # eggs per g^exponent
        fecundity_exponent=1.2,
        larval_base_survival=0.01,
        zooplankton_match_window=15.0,  # days
        maturity_energy_threshold=0.5,
        spawning_temp_threshold=4.0,  # C
        larval_duration_days=30,
        recruit_weight=0.5,  # g
        recruit_length=3.0,  # cm
    )

    return cls(
        bioenerg=bioenerg,
        predation=predation,
        foraging=foraging,
        movement=movement,
        reproduction=reproduction,
        vbgf_k_mean=0.3,
        vbgf_k_sd=0.05,
        vbgf_linf_mean=25.0,
        vbgf_linf_sd=3.0,
        max_age=10.0,
    )

baltic_defaults_els classmethod

baltic_defaults_els() -> 'SmeltParams'

Return Baltic smelt defaults with early life stages enabled.

Creates a standard baltic_defaults() parameter set and adds default EggParams, YolkSacParams, LarvalParams, and OxygenParams for mechanistic egg-to-larva modeling.

Returns:

Type	Description
SmeltParams	Fully populated parameter object with ELS enabled.

Source code in pypath/ibm/smelt.py

@classmethod
def baltic_defaults_els(cls) -> "SmeltParams":
    """Return Baltic smelt defaults with early life stages enabled.

    Creates a standard ``baltic_defaults()`` parameter set and adds
    default EggParams, YolkSacParams, LarvalParams, and OxygenParams
    for mechanistic egg-to-larva modeling.

    Returns
    -------
    SmeltParams
        Fully populated parameter object with ELS enabled.
    """
    base = cls.baltic_defaults()
    base.egg = EggParams()
    base.yolk_sac = YolkSacParams()
    base.larval = LarvalParams()
    base.oxygen = OxygenParams()
    return base

baltic_defaults_zonal classmethod

baltic_defaults_zonal() -> 'SmeltParams'

Return Baltic smelt defaults with ELS and zonal model enabled.

Creates a baltic_defaults_els() parameter set and adds default ZoneParams for the 3-zone Curonian Lagoon spatial model.

Returns:

Type	Description
SmeltParams	Fully populated parameter object with ELS + zones enabled.

Source code in pypath/ibm/smelt.py

@classmethod
def baltic_defaults_zonal(cls) -> "SmeltParams":
    """Return Baltic smelt defaults with ELS and zonal model enabled.

    Creates a ``baltic_defaults_els()`` parameter set and adds default
    ``ZoneParams`` for the 3-zone Curonian Lagoon spatial model.

    Returns
    -------
    SmeltParams
        Fully populated parameter object with ELS + zones enabled.
    """
    base = cls.baltic_defaults_els()
    base.zones = ZoneParams()
    return base

SmeltIBM

Bases: IBMGroup

Concrete IBM group implementation for Baltic smelt.

Orchestrates bioenergetics, predation, foraging, and reproduction modules for an age-structured population of super-individuals representing Osmerus eperlanus.

Parameters:

Name	Type	Description	Default
group_index	int	One-based index of this group in the Ecopath/Ecosim model (0 is reserved for the "Outside" placeholder).	required
n_groups	int	Total number of functional groups in the model.	required
params	SmeltParams	Species-specific parameters for all IBM sub-modules.	required

Source code in pypath/ibm/smelt.py

class SmeltIBM(IBMGroup):
    """Concrete IBM group implementation for Baltic smelt.

    Orchestrates bioenergetics, predation, foraging, and reproduction
    modules for an age-structured population of super-individuals
    representing Osmerus eperlanus.

    Parameters
    ----------
    group_index : int
        One-based index of this group in the Ecopath/Ecosim model
        (0 is reserved for the "Outside" placeholder).
    n_groups : int
        Total number of functional groups in the model.
    params : SmeltParams
        Species-specific parameters for all IBM sub-modules.
    """

    def __init__(
        self,
        group_index: int,
        n_groups: int,
        params: SmeltParams,
    ) -> None:
        super().__init__(group_index, n_groups)
        self.params = params
        self._last_consumption: np.ndarray = np.zeros(n_groups)
        self._next_id: int = 0
        self._rng: np.random.Generator = np.random.default_rng(42)
        self._o2_warning_issued: bool = False
        self._temp_warning_issued: bool = False

    def initialize_from_ecosim(
        self,
        biomass: float,
        params: Dict[str, Any],
        n_super_individuals: int = 500,
    ) -> None:
        """Initialize an age-structured population from Ecosim equilibrium.

        Creates *n_super_individuals* super-individuals with ages
        distributed from 0.5 to max_age.  Lengths are computed from the
        Von Bertalanffy growth function, weights from inverse allometry,
        and ``n_represented`` from an exponential survival curve scaled
        so that total biomass matches the input.

        Parameters
        ----------
        biomass : float
            Initial total biomass (tonnes) from Ecopath.
        params : Dict[str, Any]
            Additional species-specific parameters (currently unused;
            all parameters come from ``self.params``).
        n_super_individuals : int, optional
            Number of super-individuals to create (default 500).
        """
        sp = self.params

        # Distribute ages uniformly from 0.5 to max_age
        ages = np.linspace(0.5, sp.max_age, n_super_individuals)

        # Draw individual VBGF parameters with slight variation
        k_vals = self._rng.normal(sp.vbgf_k_mean, sp.vbgf_k_sd, n_super_individuals)
        k_vals = np.clip(k_vals, 0.05, None)  # K must be positive

        linf_vals = self._rng.normal(
            sp.vbgf_linf_mean, sp.vbgf_linf_sd, n_super_individuals
        )
        linf_vals = np.clip(linf_vals, 5.0, None)  # Linf must be positive

        # Compute lengths from Von Bertalanffy: L = Linf * (1 - exp(-K * age))
        lengths = linf_vals * (1.0 - np.exp(-k_vals * ages))
        lengths = np.clip(lengths, 0.1, None)

        # Compute weights from inverse allometry: weight = (length / a)^(1/b)
        a = sp.bioenerg.a_length
        b = sp.bioenerg.b_length
        weights = (lengths / a) ** (1.0 / b)
        weights = np.clip(weights, 0.1, None)

        # Exponential survival curve: relative abundance decreases with age
        # Use a natural mortality rate of ~0.5/yr as a reasonable default
        natural_mortality = 0.5
        survival_weights = np.exp(-natural_mortality * ages)

        # Compute n_represented to match target biomass.
        # Each individual contributes: n_represented_i * weight_i / 1e6 tonnes.
        # We want sum_i(n_represented_i * weight_i / 1e6) = biomass.
        # Set n_represented_i proportional to survival_weights and solve for
        # the scaling constant.
        # total_biomass = scale * sum(survival_weights * weights) / 1e6
        raw_biomass_contribution = survival_weights * weights
        total_raw = np.sum(raw_biomass_contribution)

        if total_raw <= 0.0:
            logger.warning(
                "Cannot initialize SmeltIBM: total raw biomass contribution is zero."
            )
            return

        # scale factor so that sum(scale * survival_weights * weight / 1e6) = biomass
        scale = biomass * 1e6 / total_raw
        n_represented = survival_weights * scale

        # Maturity: assume fish mature at age >= 2 years
        maturity_age = 2.0

        sexes = self._rng.integers(0, 2, size=n_super_individuals)

        individuals: List[SuperIndividual] = []
        for i in range(n_super_individuals):
            ind = SuperIndividual(
                id=i,
                n_represented=float(n_represented[i]),
                weight=float(weights[i]),
                length=float(lengths[i]),
                age=float(ages[i]),
                energy_reserve=float(weights[i]) * 0.1,  # initial reserve
                patch_idx=0,
                is_mature=bool(ages[i] >= maturity_age),
                sex=int(sexes[i]),
            )
            individuals.append(ind)

        self.individuals = individuals
        self._next_id = n_super_individuals
        self._last_consumption = np.zeros(self.n_groups)

    def get_aggregate_biomass(self) -> float:
        """Return total biomass (tonnes) across all super-individuals.

        Returns
        -------
        float
            Sum of ``total_biomass_tonnes()`` for each individual.
        """
        return sum(ind.total_biomass_tonnes() for ind in self.individuals)

    def get_consumption_by_prey(self) -> np.ndarray:
        """Return the consumption vector from the last time step.

        Returns
        -------
        np.ndarray
            1-D array of shape ``(n_groups,)`` with biomass consumed
            from each prey group.
        """
        return self._last_consumption.copy()

    def _aggregate_by_patch(self, n_patches: int) -> np.ndarray:
        """Aggregate individual biomass by spatial patch.

        Parameters
        ----------
        n_patches : int
            Number of spatial patches.

        Returns
        -------
        np.ndarray
            1-D array of shape ``(n_patches,)`` with total biomass per patch.
        """
        if not self.individuals:
            return np.zeros(n_patches)
        patches = np.array([ind.patch_idx for ind in self.individuals])
        biomasses = np.array([ind.total_biomass_tonnes() for ind in self.individuals])
        valid = (patches >= 0) & (patches < n_patches)
        result = np.zeros(n_patches)
        np.add.at(result, patches[valid], biomasses[valid])
        return result

    def _consolidate_population(self) -> None:
        """Merge similar super-individuals when population exceeds the cap.

        Groups individuals by ``(life_stage, patch_idx)`` and merges the
        two most similar (by weight) within each group until the total count
        is at or below ``self.params.max_super_individuals``.  Biomass is
        conserved exactly: the merged individual inherits the weighted-mean
        weight and recomputes length from allometry.
        """
        cap = self.params.max_super_individuals
        if len(self.individuals) <= cap:
            return

        sp = self.params

        # Group by (life_stage, patch_idx)
        groups: Dict[tuple, List[SuperIndividual]] = defaultdict(list)
        for ind in self.individuals:
            groups[(ind.life_stage, ind.patch_idx, ind.sex)].append(ind)

        # Within each group, sort by weight and merge closest pairs
        while sum(len(v) for v in groups.values()) > cap:
            # Find the largest group to merge within
            largest_key = max(groups, key=lambda k: len(groups[k]))
            grp = groups[largest_key]
            if len(grp) < 2:
                break
            grp.sort(key=lambda x: x.weight)

            # Merge the two closest by weight (adjacent after sort)
            best_idx = 0
            best_diff = float("inf")
            for j in range(len(grp) - 1):
                diff = abs(grp[j + 1].weight - grp[j].weight)
                if diff < best_diff:
                    best_diff = diff
                    best_idx = j

            a = grp[best_idx]
            b = grp[best_idx + 1]

            # Save allometric params BEFORE mutation
            a_n = a.n_represented
            b_n = b.n_represented
            total_n = a_n + b_n

            if total_n > 0:
                # Weighted-mean weight (conserves biomass exactly)
                merged_weight = (a_n * a.weight + b_n * b.weight) / total_n
                merged_age = (a_n * a.age + b_n * b.age) / total_n
                merged_energy = (
                    a_n * a.energy_reserve + b_n * b.energy_reserve
                ) / total_n
                merged_dd = (a_n * a.degree_days + b_n * b.degree_days) / total_n
                merged_yolk = (
                    a_n * a.yolk_energy_kj + b_n * b.yolk_energy_kj
                ) / total_n

                a.n_represented = total_n
                a.weight = merged_weight
                a.age = merged_age
                a.energy_reserve = merged_energy
                a.degree_days = merged_dd
                a.yolk_energy_kj = merged_yolk
                a.starvation_days = (
                    a_n * a.starvation_days + b_n * b.starvation_days
                ) / total_n
                a.is_mature = a.is_mature or b.is_mature

                # Recompute length from allometry
                if a.life_stage < 3 and sp.larval is not None:
                    a.length = (
                        sp.larval.a_length_larval
                        * max(a.weight, 0.0) ** sp.larval.b_length_larval
                    )
                else:
                    a.length = (
                        sp.bioenerg.a_length
                        * max(a.weight, 0.0) ** sp.bioenerg.b_length
                    )

            # Remove b from the group
            grp.pop(best_idx + 1)

        # Rebuild self.individuals from groups
        self.individuals = []
        for grp in groups.values():
            self.individuals.extend(grp)

    def _resolve_forcing(self, env_forcing: dict, patch_idx: int) -> dict:
        """Resolve environment for a specific zone/patch.

        When ``env_forcing`` contains a ``'zone_forcing'`` dict keyed by
        zone index, the zone-specific values override the global defaults.
        If no zone override exists for the given *patch_idx*, the original
        ``env_forcing`` is returned unchanged.

        Parameters
        ----------
        env_forcing : dict
            Global environmental forcing dict, optionally containing a
            ``'zone_forcing'`` key mapping zone indices to override dicts.
        patch_idx : int
            Zone or patch index to resolve for.

        Returns
        -------
        dict
            Resolved forcing dict with zone-specific overrides applied.
        """
        zone_forcing = env_forcing.get("zone_forcing")
        if zone_forcing and patch_idx in zone_forcing:
            resolved = dict(env_forcing)  # shallow copy
            resolved.update(zone_forcing[patch_idx])
            return resolved
        return env_forcing

    def compute_step(
        self,
        prey_available: np.ndarray,
        predation_pressure: float,
        env_forcing: Dict[str, Any],
        dt: float,
        spatial_context: Optional[SpatialContext] = None,
    ) -> IBMStepResult:
        """Advance the SmeltIBM population by one time step.

        Executes up to five phases:

        1. **Forage + Grow**: For each individual, compute adaptive foraging
           allocation, then update weight and energy via bioenergetics.
        2. **Reproduce**: Mature females spawn; surviving larvae create recruits.
        3. **Predation mortality**: Apply size-structured predation.
        4. **Bookkeeping**: Age individuals, add recruits, remove senescent fish.
        5. **Spatial movement** (optional): Move individuals between patches
           when a spatial context is provided.

        Parameters
        ----------
        prey_available : np.ndarray
            1-D array of shape ``(n_groups,)`` giving available biomass
            per prey group.
        predation_pressure : float
            Total predation mortality rate on this group (yr^-1).
        env_forcing : Dict[str, Any]
            Environmental forcing with keys like ``"temperature"``,
            ``"month"``, ``"zoo_peak_day"``.
        dt : float
            Time step size (fraction of a year).
        spatial_context : SpatialContext, optional
            Spatial patch data for Ecospace simulations. When ``None``
            (default), no spatial movement is performed.

        Returns
        -------
        IBMStepResult
            Aggregated results of this time step.
        """
        sp = self.params
        temperature = env_forcing.get("temperature", 10.0)
        if "temperature" not in env_forcing and not self._temp_warning_issued:
            logger.debug(
                "SmeltIBM: 'temperature' missing from env_forcing, defaulting to 10.0"
            )
            self._temp_warning_issued = True
        month = env_forcing.get("month", 6)
        zoo_peak_day = env_forcing.get("zoo_peak_day", 120.0)

        biomass_before = self.get_aggregate_biomass()

        # Pre-compute zone forcing cache to avoid redundant _resolve_forcing calls
        _zone_cache: Dict[int, dict] = {}

        def _get_zone_forcing(zone_idx: int) -> dict:
            if zone_idx not in _zone_cache:
                _zone_cache[zone_idx] = self._resolve_forcing(env_forcing, zone_idx)
            return _zone_cache[zone_idx]

        # Convert prey_available ndarray to dict for adaptive_forage
        prey_dict: Dict[int, float] = {}
        for idx in range(len(prey_available)):
            if prey_available[idx] > 0.0:
                prey_dict[idx] = float(prey_available[idx])

        # Accumulate consumption across all individuals
        total_consumption = np.zeros(self.n_groups)

        # ================================================================
        # Phase 1a: Egg development (degree-day accumulation + hatching)
        # ================================================================
        if sp.egg is not None:
            dt_days = dt * 365.0
            if (
                "dissolved_oxygen" not in env_forcing
                and "o2" not in env_forcing
                and not self._o2_warning_issued
            ):
                logger.debug(
                    "SmeltIBM: 'dissolved_oxygen'/'o2' missing from env_forcing, "
                    "defaulting to %.1f (no limitation)",
                    _DEFAULT_O2,
                )
                self._o2_warning_issued = True
            active_individuals: List[SuperIndividual] = []

            for ind in self.individuals:
                if ind.life_stage == 0:
                    # Resolve zone-specific forcing for this individual
                    ind_env = _get_zone_forcing(ind.patch_idx)
                    ind_temp = ind_env.get("temperature", temperature)
                    ind_o2 = ind_env.get(
                        "dissolved_oxygen", ind_env.get("o2", _DEFAULT_O2)
                    )
                    # Accumulate degree-days
                    ind.degree_days = accumulate_degree_days(
                        ind.degree_days, ind_temp, sp.egg.t_zero, dt_days
                    )
                    # Apply egg mortality
                    hypoxia_rate = (
                        sp.oxygen.hypoxia_mortality_rate
                        if sp.oxygen is not None
                        else 0.5
                    )
                    ind.n_represented = apply_egg_mortality(
                        n_represented=ind.n_represented,
                        background_rate=sp.egg.background_mortality_rate,
                        dt_days=dt_days,
                        o2=ind_o2,
                        o2_lethal=sp.egg.o2_lethal,
                        degree_days=ind.degree_days,
                        dd_mortality=sp.egg.dd_mortality,
                        hypoxia_mortality_rate=hypoxia_rate,
                    )
                    if ind.n_represented < 1.0:
                        continue  # egg cohort dead
                    # Check hatching
                    if check_hatching(ind.degree_days, sp.egg.dd_hatch):
                        ind.life_stage = 1
                        if sp.yolk_sac is not None:
                            ind.yolk_energy_kj = sp.yolk_sac.initial_yolk_kj
                    active_individuals.append(ind)
                elif ind.life_stage == 1 and sp.yolk_sac is not None:
                    # Phase 1b: Yolk-sac depletion + first feeding
                    # Resolve zone-specific forcing for this individual
                    ind_env = _get_zone_forcing(ind.patch_idx)
                    ind_temp = ind_env.get("temperature", temperature)
                    ind_o2 = ind_env.get(
                        "dissolved_oxygen", ind_env.get("o2", _DEFAULT_O2)
                    )
                    yolk_rate = compute_yolk_depletion(
                        weight=ind.weight,
                        temperature=ind_temp,
                        rs_a_larval=sp.larval.rs_a_larval if sp.larval else 0.12,
                        rs_b=sp.bioenerg.rb,
                        q10=sp.bioenerg.q10,
                        t_ref=sp.bioenerg.t_ref,
                        oxycal=sp.yolk_sac.oxycal_kj_per_g_o2,
                        dt_days=dt_days,
                    )
                    # Oxygen stress: accelerate yolk depletion under hypoxia
                    if sp.oxygen is not None and ind_o2 < sp.oxygen.pcrit_yolk_sac:
                        stress_factor = 1.0 + 0.5 * (
                            1.0 - ind_o2 / sp.oxygen.pcrit_yolk_sac
                        )
                        yolk_rate *= stress_factor
                    ind.yolk_energy_kj = max(0.0, ind.yolk_energy_kj - yolk_rate)

                    # Oxygen lethal mortality for yolk-sac larvae
                    if (
                        sp.oxygen is not None
                        and sp.oxygen.o2_lethal_yolk_sac > 0.0
                        and ind_o2 < sp.oxygen.o2_lethal_yolk_sac
                    ):
                        o2_mort_rate = (
                            sp.yolk_sac.background_mortality_rate
                            + sp.oxygen.hypoxia_mortality_rate
                            * (1.0 - ind_o2 / sp.oxygen.o2_lethal_yolk_sac)
                        )
                        ind.n_represented *= np.exp(-o2_mort_rate * dt_days)
                        if ind.n_represented < 1.0:
                            ind.n_represented = 0.0
                            continue  # cohort dead

                    # Determine zoo_density from zone-resolved env or prey_available
                    zoo_density = ind_env.get("zoo_density", None)
                    if zoo_density is None and sp.larval is not None:
                        pidx = sp.larval.zooplankton_prey_idx
                        if pidx < len(prey_available):
                            zoo_density = (
                                prey_available[pidx] * sp.larval.zoo_conversion_factor
                            )
                        else:
                            zoo_density = 0.0
                    elif zoo_density is None:
                        zoo_density = 0.0

                    status = check_first_feeding(
                        yolk_energy_kj=ind.yolk_energy_kj,
                        threshold_kj=sp.yolk_sac.first_feeding_threshold_kj,
                        zoo_density=zoo_density,
                        minimum_prey=sp.yolk_sac.minimum_prey_density,
                        starvation_days=ind.starvation_days,
                        pnr=sp.yolk_sac.point_of_no_return,
                    )

                    if status == "feed":
                        ind.life_stage = 2
                        ind.energy_reserve = ind.weight * 0.1
                        ind.starvation_days = 0.0
                        active_individuals.append(ind)
                    elif status == "dead":
                        ind.n_represented = 0.0
                        # do not add to active — effectively removed
                    elif status == "starving":
                        ind.starvation_days += dt_days
                        active_individuals.append(ind)
                    else:
                        # "yolk_sac" — still on yolk
                        active_individuals.append(ind)
                else:
                    active_individuals.append(ind)

            self.individuals = active_individuals

        # Filter: feeding individuals (life_stage >= 2) go through Phase 1;
        # eggs (0) and yolk-sac larvae (1) are excluded.
        if sp.egg is not None:
            feeding = [ind for ind in self.individuals if ind.life_stage >= 2]
            early = [ind for ind in self.individuals if ind.life_stage < 2]
        else:
            feeding = self.individuals
            early = []

        # ================================================================
        # Phase 1: Forage + Grow (with ontogenetic blending for larvae)
        # ================================================================
        n_ind = len(feeding)
        ind_consumptions = np.empty(n_ind)

        dt_days = dt * 365.0

        if n_ind > 0:
            weights = np.array([ind.weight for ind in feeding])

            # Sigmoid helper for blending
            def _sigmoid(x):
                return 1.0 / (1.0 + np.exp(-np.clip(x, -30, 30)))

            # Foraging loop with consumption blending
            for i, ind in enumerate(feeding):
                w = ind.weight

                # Resolve zone-specific forcing for this individual
                ind_env = _get_zone_forcing(ind.patch_idx)
                ind_temp = ind_env.get("temperature", temperature)

                # Determine zoo_density for larval consumption (per-zone)
                zoo_density = ind_env.get("zoo_density", None)
                if zoo_density is None and sp.larval is not None:
                    pidx = sp.larval.zooplankton_prey_idx
                    if pidx < len(prey_available):
                        zoo_density = (
                            prey_available[pidx] * sp.larval.zoo_conversion_factor
                        )
                    else:
                        zoo_density = 0.0
                elif zoo_density is None:
                    zoo_density = 0.0

                # Oxygen limitation on consumption (per-zone)
                o2_feeding = ind_env.get(
                    "dissolved_oxygen", ind_env.get("o2", _DEFAULT_O2)
                )

                # Compute per-individual oxygen scalar based on life stage
                o2_scale = 1.0
                if sp.oxygen is not None:
                    if ind.life_stage == 2:
                        pcrit = sp.oxygen.pcrit_larva
                    elif ind.life_stage == 3:
                        pcrit = sp.oxygen.pcrit_juvenile
                    else:
                        pcrit = sp.oxygen.pcrit_adult
                    o2_scale = oxygen_scalar(o2_feeding, pcrit)

                if sp.larval is not None:
                    alpha = float(
                        _sigmoid(
                            (w - sp.larval.w_forage_mid) / sp.larval.w_forage_scale
                        )
                    )
                else:
                    alpha = 1.0  # pure adaptive foraging

                # --- Larval consumption (concentration-dependent) ---
                c_larval_vec = np.zeros(self.n_groups)
                if alpha < 0.99 and sp.larval is not None:
                    lp = sp.larval
                    f_temp = thornton_lessem(
                        ind_temp,
                        CQ=lp.cmax_CQ,
                        CTO=lp.cmax_CTO,
                        CTM=lp.cmax_CTM,
                        CTL=lp.cmax_CTL,
                        CK1=lp.cmax_CK1,
                        CK4=lp.cmax_CK4,
                    )
                    cmax = lp.cmax_c_a * (w**lp.cmax_c_b) * f_temp * dt_days
                    cmax *= o2_scale  # oxygen limitation
                    denom = zoo_density + lp.k_half_zoo
                    c_larval_scalar = cmax * (zoo_density / denom) if denom > 0 else 0.0
                    pidx = lp.zooplankton_prey_idx
                    if pidx < self.n_groups:
                        c_larval_vec[pidx] = c_larval_scalar

                # --- Adaptive foraging consumption ---
                c_adaptive_vec = np.zeros(self.n_groups)
                if alpha > 0.01:
                    # Pre-existing heuristic for the adaptive foraging path;
                    # differs from Thornton-Lessem Cmax intentionally to give
                    # a simpler, weight-only upper bound on adult consumption.
                    max_cons = 0.1 * (w**0.7) * dt_days
                    max_cons *= o2_scale  # oxygen limitation
                    allocation = adaptive_forage(
                        prey_available=prey_dict,
                        max_consumption=max_cons,
                        individual_length=ind.length,
                        params=sp.foraging,
                    )
                    for prey_idx, amount in allocation.items():
                        if prey_idx < self.n_groups:
                            c_adaptive_vec[prey_idx] = amount

                # --- Blend ---
                c_total_vec = (1.0 - alpha) * c_larval_vec + alpha * c_adaptive_vec

                ind_consumptions[i] = c_total_vec.sum()
                for prey_idx in range(self.n_groups):
                    if c_total_vec[prey_idx] > 0:
                        total_consumption[prey_idx] += (
                            c_total_vec[prey_idx] * ind.n_represented / 1e6
                        )

            # Batch growth step (vectorized bioenergetics)
            energy_reserves = np.array([ind.energy_reserve for ind in feeding])
            is_mature = np.array([ind.is_mature for ind in feeding])

            # Per-individual temperature from zone-specific forcing
            temps = np.array(
                [
                    _get_zone_forcing(ind.patch_idx).get("temperature", temperature)
                    for ind in feeding
                ]
            )

            if sp.larval is not None:
                new_weights, new_energies = growth_step_batch_ontogenetic(
                    weights=weights,
                    energy_reserves=energy_reserves,
                    consumptions=ind_consumptions,
                    temperature=temps,
                    is_mature=is_mature,
                    dt=dt,
                    bioenerg_params=sp.bioenerg,
                    larval_params=sp.larval,
                )
            else:
                new_weights, new_energies = growth_step_batch(
                    weights=weights,
                    energy_reserves=energy_reserves,
                    consumptions=ind_consumptions,
                    temperature=temps,
                    is_mature=is_mature,
                    dt=dt,
                    params=sp.bioenerg,
                )

            # Batch allometric length (use larval allometry based on life_stage)
            if sp.larval is not None:
                life_stages = np.array([ind.life_stage for ind in feeding])
                new_lengths = np.where(
                    life_stages < 3,
                    sp.larval.a_length_larval
                    * np.maximum(new_weights, 1e-6) ** sp.larval.b_length_larval,
                    sp.bioenerg.a_length
                    * np.maximum(new_weights, 1e-6) ** sp.bioenerg.b_length,
                )
            else:
                new_lengths = (
                    sp.bioenerg.a_length
                    * np.maximum(new_weights, 0.0) ** sp.bioenerg.b_length
                )

            # Write back to individuals
            for i, ind in enumerate(feeding):
                ind.weight = float(new_weights[i])
                ind.energy_reserve = float(new_energies[i])
                ind.length = float(new_lengths[i])

        # Apply oxygen-dependent lethal mortality to feeding larvae
        if sp.oxygen is not None and n_ind > 0:
            dt_days_mort = dt * 365.0
            for ind in feeding:
                if ind.life_stage == 2:
                    ind_env_o2 = _get_zone_forcing(ind.patch_idx)
                    o2_val = ind_env_o2.get(
                        "dissolved_oxygen", ind_env_o2.get("o2", _DEFAULT_O2)
                    )
                    if (
                        sp.oxygen.o2_lethal_larva > 0.0
                        and o2_val < sp.oxygen.o2_lethal_larva
                    ):
                        lp_bg = (
                            sp.larval.background_mortality_rate if sp.larval else 0.01
                        )
                        o2_mort = lp_bg + sp.oxygen.hypoxia_mortality_rate * (
                            1.0 - o2_val / sp.oxygen.o2_lethal_larva
                        )
                        ind.n_represented *= np.exp(-o2_mort * dt_days_mort)

            # Remove dead individuals
            feeding = [ind for ind in feeding if ind.n_represented >= 1.0]

        # Recombine early life stages with feeding individuals
        self.individuals = early + feeding

        # ================================================================
        # Phase 2: Reproduce
        # ================================================================
        recruits: List[SuperIndividual] = []
        # Approximate spawn day from month
        spawn_day = month * 30.0

        if sp.egg is not None:
            # ELS mode: collect eggs by zone and create egg cohorts
            eggs_by_zone: Dict[int, float] = defaultdict(float)
            for ind in self.individuals:
                eggs_from = spawn(ind, temperature, sp.reproduction)
                if eggs_from > 0:
                    eggs_by_zone[ind.patch_idx] += eggs_from
            for zone_idx, zone_eggs in eggs_by_zone.items():
                if zone_eggs > 0:
                    n_cohorts = min(
                        sp.egg.max_egg_cohorts,
                        max(1, int(zone_eggs / sp.egg.eggs_per_cohort)),
                    )
                    per_cohort = zone_eggs / n_cohorts
                    for _i in range(n_cohorts):
                        egg_si = SuperIndividual(
                            id=self._next_id,
                            n_represented=per_cohort,
                            weight=sp.egg.egg_weight,
                            length=sp.egg.egg_length_cm,
                            age=0.0,
                            energy_reserve=0.0,
                            patch_idx=zone_idx,
                            is_mature=False,
                            sex=0,
                            life_stage=0,
                            degree_days=0.0,
                        )
                        recruits.append(egg_si)
                        self._next_id += 1
        else:
            # Legacy mode: create recruits directly via create_recruits()
            for ind in self.individuals:
                total_eggs = spawn(ind, temperature, sp.reproduction)
                if total_eggs > 0.0:
                    new_recruits = create_recruits(
                        total_eggs=total_eggs,
                        spawn_day=spawn_day,
                        zoo_peak_day=zoo_peak_day,
                        patch_idx=ind.patch_idx,
                        next_id=self._next_id,
                        params=sp.reproduction,
                        n_super_individuals=1,
                    )
                    self._next_id += len(new_recruits)
                    recruits.extend(new_recruits)

        # ================================================================
        # Phase 3: Apply predation mortality
        # ================================================================
        n_before_predation = sum(ind.n_represented for ind in self.individuals)

        survivors = apply_predation_mortality(
            individuals=self.individuals,
            total_mortality_rate=predation_pressure,
            dt=dt,
            params=sp.predation,
        )

        n_after_predation = sum(ind.n_represented for ind in survivors)
        mortality_count = max(0.0, n_before_predation - n_after_predation)

        # ================================================================
        # Phase 4: Add recruits, age individuals, remove senescent fish
        # ================================================================
        # Age all survivors
        for ind in survivors:
            ind.age += dt

        # Juvenile transition: advance life_stage 2→3 when length threshold met
        if sp.larval is not None:
            juv_len = sp.larval.juvenile_length_cm
            for ind in survivors:
                if ind.life_stage == 2 and ind.length >= juv_len:
                    ind.life_stage = 3

        # Remove fish that exceeded max_age
        survivors = [ind for ind in survivors if ind.age <= sp.max_age]

        # Add recruits
        survivors.extend(recruits)

        # Update population
        self.individuals = survivors

        # Consolidate if over cap
        if len(self.individuals) > sp.max_super_individuals:
            self._consolidate_population()

        # Record consumption
        self._last_consumption = total_consumption

        # ================================================================
        # Phase 5: Spatial movement
        # ================================================================
        # Ontogenetic habitat constraints: allowed zones per life stage
        _ALLOWED_ZONES = {
            0: {0},  # Eggs: river only (sessile)
            1: {0, 1},  # Yolk-sac: river, lagoon
            2: {0, 1},  # Larvae: river, lagoon
            3: {1, 2},  # Juvenile: lagoon, coastal
            4: {0, 1, 2},  # Adult: all zones
        }

        patch_biomass = None
        if spatial_context is not None:
            for ind in self.individuals:
                probs = calculate_movement_probabilities(
                    current_patch=ind.patch_idx,
                    adjacency=spatial_context.adjacency,
                    habitat_quality=spatial_context.habitat_quality,
                    food_density=spatial_context.food_density,
                    predator_density=spatial_context.predator_density,
                    params=sp.movement,
                )
                ind.patch_idx = int(self._rng.choice(len(probs), p=probs))

            patch_biomass = self._aggregate_by_patch(spatial_context.n_patches)

        elif sp.zones is not None:
            # Standalone zonal mode (3-zone Curonian Lagoon model)
            conn = sp.zones.connectivity
            n_zones = len(conn)

            for ind in self.individuals:
                # Eggs: sessile — no movement
                if ind.life_stage == 0:
                    continue

                # Spawning migration: mature adults migrate to river
                if ind.life_stage == 4 and ind.is_mature:
                    ind_env = _get_zone_forcing(ind.patch_idx)
                    ind_temp = ind_env.get("temperature", temperature)
                    if should_migrate(ind_temp, month, sp.movement):
                        ind.patch_idx = 0  # migrate to spawning zone
                        continue

                # Passive drift for yolk-sac and larvae
                if ind.life_stage in (1, 2):
                    drift_rate = sp.zones.base_drift_rate
                    # Larvae have reduced drift as they grow
                    if ind.life_stage == 2:
                        drift_rate *= max(0.1, 1.0 - ind.length / 2.0)
                    if self._rng.random() < drift_rate:
                        probs = conn[ind.patch_idx]
                        new_zone = int(self._rng.choice(n_zones, p=probs))
                        allowed = _ALLOWED_ZONES.get(ind.life_stage, {0, 1, 2})
                        if new_zone in allowed:
                            ind.patch_idx = new_zone
                    continue

                # Active movement for juveniles and non-migrating adults
                if ind.life_stage >= 3:
                    probs = conn[ind.patch_idx].copy().astype(float)
                    # O2 behavioral avoidance: weight connectivity by oxygen quality
                    if sp.oxygen is not None:
                        pcrit = (
                            sp.oxygen.pcrit_juvenile
                            if ind.life_stage == 3
                            else sp.oxygen.pcrit_adult
                        )
                        for z in range(n_zones):
                            z_env = _get_zone_forcing(z)
                            z_o2 = z_env.get(
                                "dissolved_oxygen", z_env.get("o2", _DEFAULT_O2)
                            )
                            probs[z] *= oxygen_scalar(
                                z_o2, pcrit
                            ) * sp.oxygen.oxygen_avoidance_weight + (
                                1.0 - sp.oxygen.oxygen_avoidance_weight
                            )
                        # Renormalize
                        p_sum = probs.sum()
                        if p_sum > 0:
                            probs /= p_sum
                        else:
                            probs = conn[ind.patch_idx].copy().astype(float)
                    new_zone = int(self._rng.choice(n_zones, p=probs))
                    allowed = _ALLOWED_ZONES.get(ind.life_stage, {0, 1, 2})
                    if new_zone in allowed:
                        ind.patch_idx = new_zone

            patch_biomass = self._aggregate_by_patch(n_zones)

        # Compute results
        biomass_after = self.get_aggregate_biomass()
        production = biomass_after - biomass_before
        recruitment_count = sum(r.n_represented for r in recruits)

        return IBMStepResult(
            biomass=biomass_after,
            production=production,
            consumption_by_prey=total_consumption,
            mortality_count=mortality_count,
            recruitment_count=recruitment_count,
            patch_biomass=patch_biomass,
        )

initialize_from_ecosim

initialize_from_ecosim(biomass: float, params: Dict[str, Any], n_super_individuals: int = 500) -> None

Initialize an age-structured population from Ecosim equilibrium.

Creates n_super_individuals super-individuals with ages distributed from 0.5 to max_age. Lengths are computed from the Von Bertalanffy growth function, weights from inverse allometry, and n_represented from an exponential survival curve scaled so that total biomass matches the input.

Parameters:

Name	Type	Description	Default
biomass	float	Initial total biomass (tonnes) from Ecopath.	required
params	Dict[str, Any]	Additional species-specific parameters (currently unused; all parameters come from self.params).	required
n_super_individuals	int	Number of super-individuals to create (default 500).	500

Source code in pypath/ibm/smelt.py

def initialize_from_ecosim(
    self,
    biomass: float,
    params: Dict[str, Any],
    n_super_individuals: int = 500,
) -> None:
    """Initialize an age-structured population from Ecosim equilibrium.

    Creates *n_super_individuals* super-individuals with ages
    distributed from 0.5 to max_age.  Lengths are computed from the
    Von Bertalanffy growth function, weights from inverse allometry,
    and ``n_represented`` from an exponential survival curve scaled
    so that total biomass matches the input.

    Parameters
    ----------
    biomass : float
        Initial total biomass (tonnes) from Ecopath.
    params : Dict[str, Any]
        Additional species-specific parameters (currently unused;
        all parameters come from ``self.params``).
    n_super_individuals : int, optional
        Number of super-individuals to create (default 500).
    """
    sp = self.params

    # Distribute ages uniformly from 0.5 to max_age
    ages = np.linspace(0.5, sp.max_age, n_super_individuals)

    # Draw individual VBGF parameters with slight variation
    k_vals = self._rng.normal(sp.vbgf_k_mean, sp.vbgf_k_sd, n_super_individuals)
    k_vals = np.clip(k_vals, 0.05, None)  # K must be positive

    linf_vals = self._rng.normal(
        sp.vbgf_linf_mean, sp.vbgf_linf_sd, n_super_individuals
    )
    linf_vals = np.clip(linf_vals, 5.0, None)  # Linf must be positive

    # Compute lengths from Von Bertalanffy: L = Linf * (1 - exp(-K * age))
    lengths = linf_vals * (1.0 - np.exp(-k_vals * ages))
    lengths = np.clip(lengths, 0.1, None)

    # Compute weights from inverse allometry: weight = (length / a)^(1/b)
    a = sp.bioenerg.a_length
    b = sp.bioenerg.b_length
    weights = (lengths / a) ** (1.0 / b)
    weights = np.clip(weights, 0.1, None)

    # Exponential survival curve: relative abundance decreases with age
    # Use a natural mortality rate of ~0.5/yr as a reasonable default
    natural_mortality = 0.5
    survival_weights = np.exp(-natural_mortality * ages)

    # Compute n_represented to match target biomass.
    # Each individual contributes: n_represented_i * weight_i / 1e6 tonnes.
    # We want sum_i(n_represented_i * weight_i / 1e6) = biomass.
    # Set n_represented_i proportional to survival_weights and solve for
    # the scaling constant.
    # total_biomass = scale * sum(survival_weights * weights) / 1e6
    raw_biomass_contribution = survival_weights * weights
    total_raw = np.sum(raw_biomass_contribution)

    if total_raw <= 0.0:
        logger.warning(
            "Cannot initialize SmeltIBM: total raw biomass contribution is zero."
        )
        return

    # scale factor so that sum(scale * survival_weights * weight / 1e6) = biomass
    scale = biomass * 1e6 / total_raw
    n_represented = survival_weights * scale

    # Maturity: assume fish mature at age >= 2 years
    maturity_age = 2.0

    sexes = self._rng.integers(0, 2, size=n_super_individuals)

    individuals: List[SuperIndividual] = []
    for i in range(n_super_individuals):
        ind = SuperIndividual(
            id=i,
            n_represented=float(n_represented[i]),
            weight=float(weights[i]),
            length=float(lengths[i]),
            age=float(ages[i]),
            energy_reserve=float(weights[i]) * 0.1,  # initial reserve
            patch_idx=0,
            is_mature=bool(ages[i] >= maturity_age),
            sex=int(sexes[i]),
        )
        individuals.append(ind)

    self.individuals = individuals
    self._next_id = n_super_individuals
    self._last_consumption = np.zeros(self.n_groups)

get_aggregate_biomass

get_aggregate_biomass() -> float

Return total biomass (tonnes) across all super-individuals.

Returns:

Type	Description
float	Sum of total_biomass_tonnes() for each individual.

Source code in pypath/ibm/smelt.py

def get_aggregate_biomass(self) -> float:
    """Return total biomass (tonnes) across all super-individuals.

    Returns
    -------
    float
        Sum of ``total_biomass_tonnes()`` for each individual.
    """
    return sum(ind.total_biomass_tonnes() for ind in self.individuals)

get_consumption_by_prey

get_consumption_by_prey() -> np.ndarray

Return the consumption vector from the last time step.

Returns:

Type	Description
ndarray	1-D array of shape (n_groups,) with biomass consumed from each prey group.

Source code in pypath/ibm/smelt.py

def get_consumption_by_prey(self) -> np.ndarray:
    """Return the consumption vector from the last time step.

    Returns
    -------
    np.ndarray
        1-D array of shape ``(n_groups,)`` with biomass consumed
        from each prey group.
    """
    return self._last_consumption.copy()

compute_step

compute_step(prey_available: ndarray, predation_pressure: float, env_forcing: Dict[str, Any], dt: float, spatial_context: Optional[SpatialContext] = None) -> IBMStepResult

Advance the SmeltIBM population by one time step.

Executes up to five phases:

1. Forage + Grow: For each individual, compute adaptive foraging allocation, then update weight and energy via bioenergetics.
2. Reproduce: Mature females spawn; surviving larvae create recruits.
3. Predation mortality: Apply size-structured predation.
4. Bookkeeping: Age individuals, add recruits, remove senescent fish.
5. Spatial movement (optional): Move individuals between patches when a spatial context is provided.

Parameters:

Name	Type	Description	Default
prey_available	ndarray	1-D array of shape (n_groups,) giving available biomass per prey group.	required
predation_pressure	float	Total predation mortality rate on this group (yr^-1).	required
env_forcing	Dict[str, Any]	Environmental forcing with keys like "temperature", "month", "zoo_peak_day".	required
dt	float	Time step size (fraction of a year).	required
spatial_context	SpatialContext	Spatial patch data for Ecospace simulations. When None (default), no spatial movement is performed.	None

Returns:

Type	Description
IBMStepResult	Aggregated results of this time step.

Source code in pypath/ibm/smelt.py

def compute_step(
    self,
    prey_available: np.ndarray,
    predation_pressure: float,
    env_forcing: Dict[str, Any],
    dt: float,
    spatial_context: Optional[SpatialContext] = None,
) -> IBMStepResult:
    """Advance the SmeltIBM population by one time step.

    Executes up to five phases:

    1. **Forage + Grow**: For each individual, compute adaptive foraging
       allocation, then update weight and energy via bioenergetics.
    2. **Reproduce**: Mature females spawn; surviving larvae create recruits.
    3. **Predation mortality**: Apply size-structured predation.
    4. **Bookkeeping**: Age individuals, add recruits, remove senescent fish.
    5. **Spatial movement** (optional): Move individuals between patches
       when a spatial context is provided.

    Parameters
    ----------
    prey_available : np.ndarray
        1-D array of shape ``(n_groups,)`` giving available biomass
        per prey group.
    predation_pressure : float
        Total predation mortality rate on this group (yr^-1).
    env_forcing : Dict[str, Any]
        Environmental forcing with keys like ``"temperature"``,
        ``"month"``, ``"zoo_peak_day"``.
    dt : float
        Time step size (fraction of a year).
    spatial_context : SpatialContext, optional
        Spatial patch data for Ecospace simulations. When ``None``
        (default), no spatial movement is performed.

    Returns
    -------
    IBMStepResult
        Aggregated results of this time step.
    """
    sp = self.params
    temperature = env_forcing.get("temperature", 10.0)
    if "temperature" not in env_forcing and not self._temp_warning_issued:
        logger.debug(
            "SmeltIBM: 'temperature' missing from env_forcing, defaulting to 10.0"
        )
        self._temp_warning_issued = True
    month = env_forcing.get("month", 6)
    zoo_peak_day = env_forcing.get("zoo_peak_day", 120.0)

    biomass_before = self.get_aggregate_biomass()

    # Pre-compute zone forcing cache to avoid redundant _resolve_forcing calls
    _zone_cache: Dict[int, dict] = {}

    def _get_zone_forcing(zone_idx: int) -> dict:
        if zone_idx not in _zone_cache:
            _zone_cache[zone_idx] = self._resolve_forcing(env_forcing, zone_idx)
        return _zone_cache[zone_idx]

    # Convert prey_available ndarray to dict for adaptive_forage
    prey_dict: Dict[int, float] = {}
    for idx in range(len(prey_available)):
        if prey_available[idx] > 0.0:
            prey_dict[idx] = float(prey_available[idx])

    # Accumulate consumption across all individuals
    total_consumption = np.zeros(self.n_groups)

    # ================================================================
    # Phase 1a: Egg development (degree-day accumulation + hatching)
    # ================================================================
    if sp.egg is not None:
        dt_days = dt * 365.0
        if (
            "dissolved_oxygen" not in env_forcing
            and "o2" not in env_forcing
            and not self._o2_warning_issued
        ):
            logger.debug(
                "SmeltIBM: 'dissolved_oxygen'/'o2' missing from env_forcing, "
                "defaulting to %.1f (no limitation)",
                _DEFAULT_O2,
            )
            self._o2_warning_issued = True
        active_individuals: List[SuperIndividual] = []

        for ind in self.individuals:
            if ind.life_stage == 0:
                # Resolve zone-specific forcing for this individual
                ind_env = _get_zone_forcing(ind.patch_idx)
                ind_temp = ind_env.get("temperature", temperature)
                ind_o2 = ind_env.get(
                    "dissolved_oxygen", ind_env.get("o2", _DEFAULT_O2)
                )
                # Accumulate degree-days
                ind.degree_days = accumulate_degree_days(
                    ind.degree_days, ind_temp, sp.egg.t_zero, dt_days
                )
                # Apply egg mortality
                hypoxia_rate = (
                    sp.oxygen.hypoxia_mortality_rate
                    if sp.oxygen is not None
                    else 0.5
                )
                ind.n_represented = apply_egg_mortality(
                    n_represented=ind.n_represented,
                    background_rate=sp.egg.background_mortality_rate,
                    dt_days=dt_days,
                    o2=ind_o2,
                    o2_lethal=sp.egg.o2_lethal,
                    degree_days=ind.degree_days,
                    dd_mortality=sp.egg.dd_mortality,
                    hypoxia_mortality_rate=hypoxia_rate,
                )
                if ind.n_represented < 1.0:
                    continue  # egg cohort dead
                # Check hatching
                if check_hatching(ind.degree_days, sp.egg.dd_hatch):
                    ind.life_stage = 1
                    if sp.yolk_sac is not None:
                        ind.yolk_energy_kj = sp.yolk_sac.initial_yolk_kj
                active_individuals.append(ind)
            elif ind.life_stage == 1 and sp.yolk_sac is not None:
                # Phase 1b: Yolk-sac depletion + first feeding
                # Resolve zone-specific forcing for this individual
                ind_env = _get_zone_forcing(ind.patch_idx)
                ind_temp = ind_env.get("temperature", temperature)
                ind_o2 = ind_env.get(
                    "dissolved_oxygen", ind_env.get("o2", _DEFAULT_O2)
                )
                yolk_rate = compute_yolk_depletion(
                    weight=ind.weight,
                    temperature=ind_temp,
                    rs_a_larval=sp.larval.rs_a_larval if sp.larval else 0.12,
                    rs_b=sp.bioenerg.rb,
                    q10=sp.bioenerg.q10,
                    t_ref=sp.bioenerg.t_ref,
                    oxycal=sp.yolk_sac.oxycal_kj_per_g_o2,
                    dt_days=dt_days,
                )
                # Oxygen stress: accelerate yolk depletion under hypoxia
                if sp.oxygen is not None and ind_o2 < sp.oxygen.pcrit_yolk_sac:
                    stress_factor = 1.0 + 0.5 * (
                        1.0 - ind_o2 / sp.oxygen.pcrit_yolk_sac
                    )
                    yolk_rate *= stress_factor
                ind.yolk_energy_kj = max(0.0, ind.yolk_energy_kj - yolk_rate)

                # Oxygen lethal mortality for yolk-sac larvae
                if (
                    sp.oxygen is not None
                    and sp.oxygen.o2_lethal_yolk_sac > 0.0
                    and ind_o2 < sp.oxygen.o2_lethal_yolk_sac
                ):
                    o2_mort_rate = (
                        sp.yolk_sac.background_mortality_rate
                        + sp.oxygen.hypoxia_mortality_rate
                        * (1.0 - ind_o2 / sp.oxygen.o2_lethal_yolk_sac)
                    )
                    ind.n_represented *= np.exp(-o2_mort_rate * dt_days)
                    if ind.n_represented < 1.0:
                        ind.n_represented = 0.0
                        continue  # cohort dead

                # Determine zoo_density from zone-resolved env or prey_available
                zoo_density = ind_env.get("zoo_density", None)
                if zoo_density is None and sp.larval is not None:
                    pidx = sp.larval.zooplankton_prey_idx
                    if pidx < len(prey_available):
                        zoo_density = (
                            prey_available[pidx] * sp.larval.zoo_conversion_factor
                        )
                    else:
                        zoo_density = 0.0
                elif zoo_density is None:
                    zoo_density = 0.0

                status = check_first_feeding(
                    yolk_energy_kj=ind.yolk_energy_kj,
                    threshold_kj=sp.yolk_sac.first_feeding_threshold_kj,
                    zoo_density=zoo_density,
                    minimum_prey=sp.yolk_sac.minimum_prey_density,
                    starvation_days=ind.starvation_days,
                    pnr=sp.yolk_sac.point_of_no_return,
                )

                if status == "feed":
                    ind.life_stage = 2
                    ind.energy_reserve = ind.weight * 0.1
                    ind.starvation_days = 0.0
                    active_individuals.append(ind)
                elif status == "dead":
                    ind.n_represented = 0.0
                    # do not add to active — effectively removed
                elif status == "starving":
                    ind.starvation_days += dt_days
                    active_individuals.append(ind)
                else:
                    # "yolk_sac" — still on yolk
                    active_individuals.append(ind)
            else:
                active_individuals.append(ind)

        self.individuals = active_individuals

    # Filter: feeding individuals (life_stage >= 2) go through Phase 1;
    # eggs (0) and yolk-sac larvae (1) are excluded.
    if sp.egg is not None:
        feeding = [ind for ind in self.individuals if ind.life_stage >= 2]
        early = [ind for ind in self.individuals if ind.life_stage < 2]
    else:
        feeding = self.individuals
        early = []

    # ================================================================
    # Phase 1: Forage + Grow (with ontogenetic blending for larvae)
    # ================================================================
    n_ind = len(feeding)
    ind_consumptions = np.empty(n_ind)

    dt_days = dt * 365.0

    if n_ind > 0:
        weights = np.array([ind.weight for ind in feeding])

        # Sigmoid helper for blending
        def _sigmoid(x):
            return 1.0 / (1.0 + np.exp(-np.clip(x, -30, 30)))

        # Foraging loop with consumption blending
        for i, ind in enumerate(feeding):
            w = ind.weight

            # Resolve zone-specific forcing for this individual
            ind_env = _get_zone_forcing(ind.patch_idx)
            ind_temp = ind_env.get("temperature", temperature)

            # Determine zoo_density for larval consumption (per-zone)
            zoo_density = ind_env.get("zoo_density", None)
            if zoo_density is None and sp.larval is not None:
                pidx = sp.larval.zooplankton_prey_idx
                if pidx < len(prey_available):
                    zoo_density = (
                        prey_available[pidx] * sp.larval.zoo_conversion_factor
                    )
                else:
                    zoo_density = 0.0
            elif zoo_density is None:
                zoo_density = 0.0

            # Oxygen limitation on consumption (per-zone)
            o2_feeding = ind_env.get(
                "dissolved_oxygen", ind_env.get("o2", _DEFAULT_O2)
            )

            # Compute per-individual oxygen scalar based on life stage
            o2_scale = 1.0
            if sp.oxygen is not None:
                if ind.life_stage == 2:
                    pcrit = sp.oxygen.pcrit_larva
                elif ind.life_stage == 3:
                    pcrit = sp.oxygen.pcrit_juvenile
                else:
                    pcrit = sp.oxygen.pcrit_adult
                o2_scale = oxygen_scalar(o2_feeding, pcrit)

            if sp.larval is not None:
                alpha = float(
                    _sigmoid(
                        (w - sp.larval.w_forage_mid) / sp.larval.w_forage_scale
                    )
                )
            else:
                alpha = 1.0  # pure adaptive foraging

            # --- Larval consumption (concentration-dependent) ---
            c_larval_vec = np.zeros(self.n_groups)
            if alpha < 0.99 and sp.larval is not None:
                lp = sp.larval
                f_temp = thornton_lessem(
                    ind_temp,
                    CQ=lp.cmax_CQ,
                    CTO=lp.cmax_CTO,
                    CTM=lp.cmax_CTM,
                    CTL=lp.cmax_CTL,
                    CK1=lp.cmax_CK1,
                    CK4=lp.cmax_CK4,
                )
                cmax = lp.cmax_c_a * (w**lp.cmax_c_b) * f_temp * dt_days
                cmax *= o2_scale  # oxygen limitation
                denom = zoo_density + lp.k_half_zoo
                c_larval_scalar = cmax * (zoo_density / denom) if denom > 0 else 0.0
                pidx = lp.zooplankton_prey_idx
                if pidx < self.n_groups:
                    c_larval_vec[pidx] = c_larval_scalar

            # --- Adaptive foraging consumption ---
            c_adaptive_vec = np.zeros(self.n_groups)
            if alpha > 0.01:
                # Pre-existing heuristic for the adaptive foraging path;
                # differs from Thornton-Lessem Cmax intentionally to give
                # a simpler, weight-only upper bound on adult consumption.
                max_cons = 0.1 * (w**0.7) * dt_days
                max_cons *= o2_scale  # oxygen limitation
                allocation = adaptive_forage(
                    prey_available=prey_dict,
                    max_consumption=max_cons,
                    individual_length=ind.length,
                    params=sp.foraging,
                )
                for prey_idx, amount in allocation.items():
                    if prey_idx < self.n_groups:
                        c_adaptive_vec[prey_idx] = amount

            # --- Blend ---
            c_total_vec = (1.0 - alpha) * c_larval_vec + alpha * c_adaptive_vec

            ind_consumptions[i] = c_total_vec.sum()
            for prey_idx in range(self.n_groups):
                if c_total_vec[prey_idx] > 0:
                    total_consumption[prey_idx] += (
                        c_total_vec[prey_idx] * ind.n_represented / 1e6
                    )

        # Batch growth step (vectorized bioenergetics)
        energy_reserves = np.array([ind.energy_reserve for ind in feeding])
        is_mature = np.array([ind.is_mature for ind in feeding])

        # Per-individual temperature from zone-specific forcing
        temps = np.array(
            [
                _get_zone_forcing(ind.patch_idx).get("temperature", temperature)
                for ind in feeding
            ]
        )

        if sp.larval is not None:
            new_weights, new_energies = growth_step_batch_ontogenetic(
                weights=weights,
                energy_reserves=energy_reserves,
                consumptions=ind_consumptions,
                temperature=temps,
                is_mature=is_mature,
                dt=dt,
                bioenerg_params=sp.bioenerg,
                larval_params=sp.larval,
            )
        else:
            new_weights, new_energies = growth_step_batch(
                weights=weights,
                energy_reserves=energy_reserves,
                consumptions=ind_consumptions,
                temperature=temps,
                is_mature=is_mature,
                dt=dt,
                params=sp.bioenerg,
            )

        # Batch allometric length (use larval allometry based on life_stage)
        if sp.larval is not None:
            life_stages = np.array([ind.life_stage for ind in feeding])
            new_lengths = np.where(
                life_stages < 3,
                sp.larval.a_length_larval
                * np.maximum(new_weights, 1e-6) ** sp.larval.b_length_larval,
                sp.bioenerg.a_length
                * np.maximum(new_weights, 1e-6) ** sp.bioenerg.b_length,
            )
        else:
            new_lengths = (
                sp.bioenerg.a_length
                * np.maximum(new_weights, 0.0) ** sp.bioenerg.b_length
            )

        # Write back to individuals
        for i, ind in enumerate(feeding):
            ind.weight = float(new_weights[i])
            ind.energy_reserve = float(new_energies[i])
            ind.length = float(new_lengths[i])

    # Apply oxygen-dependent lethal mortality to feeding larvae
    if sp.oxygen is not None and n_ind > 0:
        dt_days_mort = dt * 365.0
        for ind in feeding:
            if ind.life_stage == 2:
                ind_env_o2 = _get_zone_forcing(ind.patch_idx)
                o2_val = ind_env_o2.get(
                    "dissolved_oxygen", ind_env_o2.get("o2", _DEFAULT_O2)
                )
                if (
                    sp.oxygen.o2_lethal_larva > 0.0
                    and o2_val < sp.oxygen.o2_lethal_larva
                ):
                    lp_bg = (
                        sp.larval.background_mortality_rate if sp.larval else 0.01
                    )
                    o2_mort = lp_bg + sp.oxygen.hypoxia_mortality_rate * (
                        1.0 - o2_val / sp.oxygen.o2_lethal_larva
                    )
                    ind.n_represented *= np.exp(-o2_mort * dt_days_mort)

        # Remove dead individuals
        feeding = [ind for ind in feeding if ind.n_represented >= 1.0]

    # Recombine early life stages with feeding individuals
    self.individuals = early + feeding

    # ================================================================
    # Phase 2: Reproduce
    # ================================================================
    recruits: List[SuperIndividual] = []
    # Approximate spawn day from month
    spawn_day = month * 30.0

    if sp.egg is not None:
        # ELS mode: collect eggs by zone and create egg cohorts
        eggs_by_zone: Dict[int, float] = defaultdict(float)
        for ind in self.individuals:
            eggs_from = spawn(ind, temperature, sp.reproduction)
            if eggs_from > 0:
                eggs_by_zone[ind.patch_idx] += eggs_from
        for zone_idx, zone_eggs in eggs_by_zone.items():
            if zone_eggs > 0:
                n_cohorts = min(
                    sp.egg.max_egg_cohorts,
                    max(1, int(zone_eggs / sp.egg.eggs_per_cohort)),
                )
                per_cohort = zone_eggs / n_cohorts
                for _i in range(n_cohorts):
                    egg_si = SuperIndividual(
                        id=self._next_id,
                        n_represented=per_cohort,
                        weight=sp.egg.egg_weight,
                        length=sp.egg.egg_length_cm,
                        age=0.0,
                        energy_reserve=0.0,
                        patch_idx=zone_idx,
                        is_mature=False,
                        sex=0,
                        life_stage=0,
                        degree_days=0.0,
                    )
                    recruits.append(egg_si)
                    self._next_id += 1
    else:
        # Legacy mode: create recruits directly via create_recruits()
        for ind in self.individuals:
            total_eggs = spawn(ind, temperature, sp.reproduction)
            if total_eggs > 0.0:
                new_recruits = create_recruits(
                    total_eggs=total_eggs,
                    spawn_day=spawn_day,
                    zoo_peak_day=zoo_peak_day,
                    patch_idx=ind.patch_idx,
                    next_id=self._next_id,
                    params=sp.reproduction,
                    n_super_individuals=1,
                )
                self._next_id += len(new_recruits)
                recruits.extend(new_recruits)

    # ================================================================
    # Phase 3: Apply predation mortality
    # ================================================================
    n_before_predation = sum(ind.n_represented for ind in self.individuals)

    survivors = apply_predation_mortality(
        individuals=self.individuals,
        total_mortality_rate=predation_pressure,
        dt=dt,
        params=sp.predation,
    )

    n_after_predation = sum(ind.n_represented for ind in survivors)
    mortality_count = max(0.0, n_before_predation - n_after_predation)

    # ================================================================
    # Phase 4: Add recruits, age individuals, remove senescent fish
    # ================================================================
    # Age all survivors
    for ind in survivors:
        ind.age += dt

    # Juvenile transition: advance life_stage 2→3 when length threshold met
    if sp.larval is not None:
        juv_len = sp.larval.juvenile_length_cm
        for ind in survivors:
            if ind.life_stage == 2 and ind.length >= juv_len:
                ind.life_stage = 3

    # Remove fish that exceeded max_age
    survivors = [ind for ind in survivors if ind.age <= sp.max_age]

    # Add recruits
    survivors.extend(recruits)

    # Update population
    self.individuals = survivors

    # Consolidate if over cap
    if len(self.individuals) > sp.max_super_individuals:
        self._consolidate_population()

    # Record consumption
    self._last_consumption = total_consumption

    # ================================================================
    # Phase 5: Spatial movement
    # ================================================================
    # Ontogenetic habitat constraints: allowed zones per life stage
    _ALLOWED_ZONES = {
        0: {0},  # Eggs: river only (sessile)
        1: {0, 1},  # Yolk-sac: river, lagoon
        2: {0, 1},  # Larvae: river, lagoon
        3: {1, 2},  # Juvenile: lagoon, coastal
        4: {0, 1, 2},  # Adult: all zones
    }

    patch_biomass = None
    if spatial_context is not None:
        for ind in self.individuals:
            probs = calculate_movement_probabilities(
                current_patch=ind.patch_idx,
                adjacency=spatial_context.adjacency,
                habitat_quality=spatial_context.habitat_quality,
                food_density=spatial_context.food_density,
                predator_density=spatial_context.predator_density,
                params=sp.movement,
            )
            ind.patch_idx = int(self._rng.choice(len(probs), p=probs))

        patch_biomass = self._aggregate_by_patch(spatial_context.n_patches)

    elif sp.zones is not None:
        # Standalone zonal mode (3-zone Curonian Lagoon model)
        conn = sp.zones.connectivity
        n_zones = len(conn)

        for ind in self.individuals:
            # Eggs: sessile — no movement
            if ind.life_stage == 0:
                continue

            # Spawning migration: mature adults migrate to river
            if ind.life_stage == 4 and ind.is_mature:
                ind_env = _get_zone_forcing(ind.patch_idx)
                ind_temp = ind_env.get("temperature", temperature)
                if should_migrate(ind_temp, month, sp.movement):
                    ind.patch_idx = 0  # migrate to spawning zone
                    continue

            # Passive drift for yolk-sac and larvae
            if ind.life_stage in (1, 2):
                drift_rate = sp.zones.base_drift_rate
                # Larvae have reduced drift as they grow
                if ind.life_stage == 2:
                    drift_rate *= max(0.1, 1.0 - ind.length / 2.0)
                if self._rng.random() < drift_rate:
                    probs = conn[ind.patch_idx]
                    new_zone = int(self._rng.choice(n_zones, p=probs))
                    allowed = _ALLOWED_ZONES.get(ind.life_stage, {0, 1, 2})
                    if new_zone in allowed:
                        ind.patch_idx = new_zone
                continue

            # Active movement for juveniles and non-migrating adults
            if ind.life_stage >= 3:
                probs = conn[ind.patch_idx].copy().astype(float)
                # O2 behavioral avoidance: weight connectivity by oxygen quality
                if sp.oxygen is not None:
                    pcrit = (
                        sp.oxygen.pcrit_juvenile
                        if ind.life_stage == 3
                        else sp.oxygen.pcrit_adult
                    )
                    for z in range(n_zones):
                        z_env = _get_zone_forcing(z)
                        z_o2 = z_env.get(
                            "dissolved_oxygen", z_env.get("o2", _DEFAULT_O2)
                        )
                        probs[z] *= oxygen_scalar(
                            z_o2, pcrit
                        ) * sp.oxygen.oxygen_avoidance_weight + (
                            1.0 - sp.oxygen.oxygen_avoidance_weight
                        )
                    # Renormalize
                    p_sum = probs.sum()
                    if p_sum > 0:
                        probs /= p_sum
                    else:
                        probs = conn[ind.patch_idx].copy().astype(float)
                new_zone = int(self._rng.choice(n_zones, p=probs))
                allowed = _ALLOWED_ZONES.get(ind.life_stage, {0, 1, 2})
                if new_zone in allowed:
                    ind.patch_idx = new_zone

        patch_biomass = self._aggregate_by_patch(n_zones)

    # Compute results
    biomass_after = self.get_aggregate_biomass()
    production = biomass_after - biomass_before
    recruitment_count = sum(r.n_represented for r in recruits)

    return IBMStepResult(
        biomass=biomass_after,
        production=production,
        consumption_by_prey=total_consumption,
        mortality_count=mortality_count,
        recruitment_count=recruitment_count,
        patch_biomass=patch_biomass,
    )