Numerical Solvers

The object-oriented nature of Species (Populations), Parameters, and the Model Domain allows flexibility and modularity while solving of populations over time. Numerical solvers may be developed to address unique logic associated with projects and may be swapped or combined while simulating populations.

All solvers include calls to the Model Domain pre-solve check functions:

• Error Checking
• Species Inheritance

Available, documented solvers include:

• Discrete Explicit - solvers.discrete_explicit

Error Checking

Error checking is performed using the following:

solvers.error_check(domain)

Ensure there are no data conflicts within a domain prior to calling a solver. Current checks include:

• Make sure at least one Species is included
• Make sure age (stage) groups do not have overlapping ages
• Make sure species defined in inter-species relationships are included

Parameters:domain (Domain) – Domain instance

Species Inheritance

Population and Carrying Capacity data may be applied to a Species-level instance, although Sex and AgeGroup-level species may exist in the model Domain. As such, these data are cascaded to children classes in advance of solving, which reduces complexity of numerical solvers and their ability to recognize data inheritance. Inheritance typically involves:

1. Check whether children of Species and Sex-level classes exist in the model domain
2. Collect data related to parent classes
3. Divide parent data (\(param\)) by the number of children (\(n\)), \(\frac{param}{n}\)

solvers.inherit(domain, start_time, end_time)

Divide population and carrying capacity among children of species if they exist, and apply max age truncation to all children species.

Parameters:

• domain – Domain instance
• start_time (int) – A simulation start time that coincides with that of the specified solver time range
• end_time (int) – A simulation stop time that coincides with that of the specified solver time range

Discrete Explicit

This solver is designed to solve species populations on a discrete time interval using explicit functions. At each time step between the given start time and end time, species populations from the previous time step are collected and combined with parameters from the current time step to calculate and save the populations at the current time step. As such, populations must be provided at the first time step, and the first time step is not solved in the discrete_explicit solver.

Calculation graphs are created for each species, sex, age (stage) group, at each time using the dask library, which enables:

• Disk input/output to be optimized while reading HDF5 data
• Mathematical parameter optimization
• Parallelism where possible
• Memory-management to ensure RAM usage is predictable and constrained

At each time step (starting with the second time step), populations are solved using the following methods

Note

When describing calculations, note that they are performed independently at each element in the model domain

The two primary drivers during simulations are the carrying capacity and populations. At each time step, the total populations of each sex (gender) if they are included in the domain are collected from the previous time step:

\[p_{o(t)}=\sum_{i=1}^{n}p_{a(t-1)}\]

where \(p_o\) is the initial population of each sex (gender) used in calculations, \(p_a\) is the population associated with each age and, and \(t\) is time.

If the total_density keyword argument is specified, the total population of the entire species is also calcualted, which will be used to calculate density relationships (as opposed to doing so independently for each sex or age (stage) group).

Each species, sex (gender) and age (stage) group are then iterated and populations of each are iterated and solved independently. While solving the population, parameters at the current time step are first collected for the species.

If parameters are dependent on inter-species relationships, the density of the contributing species \(\rho\) is calculated using the population \(p_o\) from the previous time step, and the carrying capacity \(k\) at the current time step:

\[\rho=\frac{p_o}{k}\]

Circularity may exist when calculating inter-species relationships (i.e. two or more species have derived \(k\) based on one another), which is why derived carrying capacity is performed with a coefficient. In such cases, the carrying capacity data are first collected for all interdependent species, and are subsequently used to calculate derived values, rather than using derived values to calculate other derived values.

Note

Stochastic variability included in parameters is applied after values are derived from inter-species relationships

Parameters collected at the beginning of the time step include:

• Mortality for each mortality driver
• The sum of Carrying Capacity - \(\sum_{i=1}^{n}k_{i(t)}\)
• Total populations for the species, males, and females (if this level of stratification exists in the domain)
• Total populations of reproducing species
• Total populations of species that contribute to density (specified in the Species (Populations) keyword arg contributes_to_density)
• Density
• The sum of Fecundity - \(\sum_{i=1}^{n}f_{i(t)}\), where \(f\) is the fecundity type

Fecundity Method

Fecundity values are scaled by two coefficients, which yields the effective fecundity:

1. A coefficient (\(y\)) derived from a lookup table (methods explained here.) based on density (\(x\)).
2. A coefficient derived from density dependency on fecundity. This value is calculated using the density_fecundity_threshold, density_fecundity_max, and fecundity_reduction_rate keyword arguments in the Fecundity class.

Density dependency on fecundity (coefficient 2 above) is calculated linearly, as follows:

\[min\{1,\frac{\rho-l}{max\{0,u-l\}}\}\cdot \Delta\]

where:

\(u\) is density_fecundity_max, the upper threshold

\(l\) is density_fecundity_threshold, the lower threshold

\(\Delta\) is fecundity_reduction_rate, the rate of fecundity reduction at the upper threshold

Prior to other calculations, total offspring (\(b\)) at the current time step are calculated using the total population of the reproducing group (ex. females) and the effective fecundity rate \(f_e\),

\[b(t)=p_o(t)f_e(t)\]

recalling that \(p_o\) is the total population from \(t-1\).

Offspring are allocated to either males, females, or neither (depending on model parameterization) of the minimum age in the domain at \(t+1\).

Following the calcualtion of offspring, all ages in the group (stage) are iterated and mortality, dispersal, and \(t+1\) propagation are completed.

Mortality Method

Mortality is applied following the calculation of offspring in four steps while solving the model. The first step includes the calculation of mortality as a result of all mortality driver parameters. Mortality as a result of each driver, \(m\), is calculated using \(p_o\):

\[m_{i(t)}=p_{o(t)}q_{i(t)}\]

where \(i\) is the mortality driver with the rate \(q\).

In the case where the sum of all mortality driver rates exceed 1, they are scaled proportionally:

\[m_{i(t)}=\frac{q_{i(t)}}{\sum_{i=1}^{n}q_{i(t)}}\cdot p_{o(t)}\]

Following the drivers, the remaining mortality calculations are a result of Implicit mortality types. These include:

• Old age; and
• Density-dependent mortality.

Populations are first moved using any provided Dispersal methods prior to applying density-dependent mortality. Density-dependent mortality is calculated using the density_threshold and density_scale keyword arguments provided to the Species (Populations) constructor:

\[\frac{max\{0,\rho-l\}}{1-l}\cdot \Delta\]

where,

\(\rho\) is the population density

\(l\) is density_threshold, and is the lower limit of density dependent mortality

\(\Delta\) is density_scale, and is the rate of density dependent mortality at a density of 1.

Old age mortality is only applied if the live_past_max keyword argument in the Species (Populations) instance is False. In this case, populations of the maximum specified age in the domain will not propagate to \(t+1\), and will be tracked as mortality as a result of old age.

Prior to propagating a population to \(t+1\), the minimum_viable_population keyword argument of the species is checked to determine whether a minimum viable population calculation should take place. This form of mortality is applied last at the current time step.

Dispersal Method

Dispersal is applied prior to implementing density-dependent mortality to give populations a chance to move before they succumb to the effects of density. Dispersal is applied using any number of the Dispersal methods inherent to species, in the order that they were applied to the Species (Populations) object.

Conversion & Immigration/Emigration

Two processes that take place during solving are conversion from one species to another, and immigration/emigration. The latter are applied by adding populations to the domain at time slices other than the first time in the simulation. Positive values will apply (immigrate) populations to the domain at a given time, while negative values will enforce a population loss from the system (emigration).

The Discrete Explicit Class

class solvers.discrete_explicit(domain, start_time, end_time, **kwargs)

Solve the domain from a start time to an end time on a discrete time step associated with the Domain.

The total populations from the previous year are used to calculate derived parameters (density, fecundity, inter-species relationships, etc.).

Parameters used during execution are stored in the domain under the params key, while magnitude related to mortality and fecundity are stored under the flux key at each time step.

To run the simulation, use .execute()

Parameters:

• domain (Domain) – Domain instance populated with at least one species
• start_time (int) – Start time of the simulation
• end_time (int) – End time of the simulation

Keyword Arguments:  

total_density (bool) –

Use the total species population for density calculations (as opposed to populations of sexes or groups) (Default: True)

calculate_parameters(species, sex, group, time)

Collect all required parameters at a time step

carrying_capacity_total(species, datasets)

Calculate the total carrying capacity for a species with the provided datasets, usually a result of a Domain.get_carrying_capacity() call. This updates self.carrying_capacity_arrays to avoid excessive or infinite recursion depth when collecting carrying capacity that may be circular.

Parameters:

• datasets – list of tuples with instance - h5py.Dataset pairs
• dependent (list) – Tracks the species used in species-dependent carrying capacity calculations

Returns:

dask array of total carrying capacity

collect_parameter(param, data)

Collect data associated with a Fecundity, CarryingCapacity, or Mortality instance. If a species is linked, data will be None and density of the linked species will be calculated.

Parameters:

• param – CarryingCapacity, Fecundity, or Mortality instance
• data – h5py.Dataset key

Returns:

A dask array with the collected data

execute()

Run the simulation

get_counter(species, group, sex, age, time)

Collect the consecutive number of years a non-zero population has existed

Parameters:

• species (str) – Species name key
• group (str) – Group name key
• sex (str) – Sex name key
• age (int) – Age of the population
• time (int) – Current model time

Returns:

Integer count

interspecies_carrying_capacity(parent_species, param)

Because inter-species dependencies may be nested or circular, the carrying capacity of all connected species is calculated herein and added to the self.carrying_capacity_arrays object. Future calls of any of the connected species will have data available, avoiding redundant calculations.

Parameters:

• parent_species (Species) – One species that is tied to others
• param (CarryingCapacity) – Input Carrying Capacity

Returns:

coefficient array for the given carrying capacity parameter

population_total(datasets, honor_density=True, honor_reproduction=False, honor_global=False)

Calculate the total population for a species with the provided datasets, usually a result of a Domain.all_population() call. This updates self.population_arrays to avoid repetitive IO

Parameters:

• species – Species key
• datasets (list) – population dataset keys, retrieved from popdyn.Domain.all_population()
• honor_density (bool) – Use the density contribution attribute from a species
• honor_reproduction (bool) – Use the density contribution to fecundity attribute from a species

Returns:

Dask array of total population

prepare_domain(start_time, end_time)

Error check and apply inheritance to domain

propagate(species, sex, group, time)

Create a graph of calculations to propagate populations and record parameters. A dict of key pointers - dask array pairs are created and returned

Note

• If age gaps exist, they will be filled by empty groups and will propagate
• Currently, female:male ratio uses the total species population even if densities are calculated at the age group level

Parameters:

• species (str) – Species key (Species.name_key)
• sex (str) – Sex ('male', 'female', or None)
• group (str) – Group key (AgeGroup.group_key)
• time (int) – Time slice

totals(all_species, time)

Collect all population and carrying capacity datasets, and calculate species totals if necessary. Dynamic calculations of carrying capacity are made by calling collect_parameter(), which may also update the carrying_capacity or population array dicts.

Parameters:

• all_species (list) – All species in the domain
• time (int) – Time slice to collect totals