Individual-Based Model (IBM) API Reference

Super-individual tracking with bioenergetics, predation, reproduction, behavior, and spatial integration.

Base Classes

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.

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).
    """

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

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	Dimensionless energy reserve index (0-1 typical range).	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

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
        Dimensionless energy reserve index (0-1 typical range).
    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).
    """

    id: int
    n_represented: float
    weight: float
    length: float
    age: float
    energy_reserve: float
    patch_idx: int
    is_mature: bool
    sex: int

    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

Bioenergetics

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

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

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

    # Handle energy reserve: surplus increases reserve, deficit drains it
    if net_energy >= 0.0:
        # Surplus: store a fraction in energy reserve
        new_energy_reserve = energy_reserve + net_energy / params.energy_density
    else:
        # Deficit: drain energy reserve first
        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))

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

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)

Predation

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

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

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 []

    # 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()

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)

Reproduction

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

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

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

Behavior (Spatial Movement)

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.

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

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)

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

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
        )
        # Apply inertia bonus to current patch
        if p == current_patch:
            score += (1.0 - params.base_speed) * score
        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
    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
    )

EwE Integration

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.

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
    prey_dict = extract_prey_availability(QQ, group_idx, n_groups)
    prey_array = np.zeros(n_groups)
    for prey_idx, rate in prey_dict.items():
        if prey_idx < n_groups:
            prey_array[prey_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
    for prey_idx in range(len(result.consumption_by_prey)):
        consumed = result.consumption_by_prey[prey_idx]
        if consumed != 0.0 and prey_idx < len(deriv):
            deriv[prey_idx] -= consumed / dt

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

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]))

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

Smelt Implementation

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.

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)

    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 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)
        month = env_forcing.get("month", 6)
        zoo_peak_day = env_forcing.get("zoo_peak_day", 120.0)

        biomass_before = self.get_aggregate_biomass()

        # 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 1: Forage + Grow
        # ================================================================
        for ind in self.individuals:
            # Maximum consumption scales allometrically with weight
            # Using a simple Cmax = 0.1 * weight^0.7 * dt * 365 (per timestep)
            max_consumption = 0.1 * (ind.weight**0.7) * dt * 365.0

            # Scale max_consumption by number represented for total group take
            # But adaptive_forage works per-individual, consumption is per-individual
            allocation = adaptive_forage(
                prey_available=prey_dict,
                max_consumption=max_consumption,
                individual_length=ind.length,
                params=sp.foraging,
            )

            # Total consumption by this super-individual (per real individual)
            ind_consumption = sum(allocation.values())

            # Accumulate consumption in tonnes:
            # per-individual consumption * n_represented / 1e6
            for prey_idx, amount in allocation.items():
                if prey_idx < self.n_groups:
                    total_consumption[prey_idx] += amount * ind.n_represented / 1e6

            # Grow: update weight and energy reserve
            new_weight, new_energy = growth_step(
                weight=ind.weight,
                energy_reserve=ind.energy_reserve,
                consumption=ind_consumption,
                temperature=temperature,
                is_mature=ind.is_mature,
                dt=dt,
                params=sp.bioenerg,
            )

            ind.weight = new_weight
            ind.energy_reserve = new_energy
            # Update length from new weight
            ind.length = allometric_length(
                ind.weight, sp.bioenerg.a_length, sp.bioenerg.b_length
            )

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

        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

        # 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

        # Record consumption
        self._last_consumption = total_consumption

        # ================================================================
        # Phase 5: Spatial movement (only when spatial context provided)
        # ================================================================
        patch_biomass = None
        if spatial_context is not None:
            from pypath.ibm.behavior import calculate_movement_probabilities

            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)

        # 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,
        )

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)
    month = env_forcing.get("month", 6)
    zoo_peak_day = env_forcing.get("zoo_peak_day", 120.0)

    biomass_before = self.get_aggregate_biomass()

    # 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 1: Forage + Grow
    # ================================================================
    for ind in self.individuals:
        # Maximum consumption scales allometrically with weight
        # Using a simple Cmax = 0.1 * weight^0.7 * dt * 365 (per timestep)
        max_consumption = 0.1 * (ind.weight**0.7) * dt * 365.0

        # Scale max_consumption by number represented for total group take
        # But adaptive_forage works per-individual, consumption is per-individual
        allocation = adaptive_forage(
            prey_available=prey_dict,
            max_consumption=max_consumption,
            individual_length=ind.length,
            params=sp.foraging,
        )

        # Total consumption by this super-individual (per real individual)
        ind_consumption = sum(allocation.values())

        # Accumulate consumption in tonnes:
        # per-individual consumption * n_represented / 1e6
        for prey_idx, amount in allocation.items():
            if prey_idx < self.n_groups:
                total_consumption[prey_idx] += amount * ind.n_represented / 1e6

        # Grow: update weight and energy reserve
        new_weight, new_energy = growth_step(
            weight=ind.weight,
            energy_reserve=ind.energy_reserve,
            consumption=ind_consumption,
            temperature=temperature,
            is_mature=ind.is_mature,
            dt=dt,
            params=sp.bioenerg,
        )

        ind.weight = new_weight
        ind.energy_reserve = new_energy
        # Update length from new weight
        ind.length = allometric_length(
            ind.weight, sp.bioenerg.a_length, sp.bioenerg.b_length
        )

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

    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

    # 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

    # Record consumption
    self._last_consumption = total_consumption

    # ================================================================
    # Phase 5: Spatial movement (only when spatial context provided)
    # ================================================================
    patch_biomass = None
    if spatial_context is not None:
        from pypath.ibm.behavior import calculate_movement_probabilities

        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)

    # 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,
    )

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()

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)

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

    @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 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,
    )