Degree-class SIR on uncorrelated networks
The SIR epidemic on an uncorrelated network of power-law degrees, in heterogeneous mean field.
- Name:
heterogeneous_mean_field_sir - Field:
biology/epidemiology - Formalism: ordinary differential equations
- References: Bertin 2021; Moreno, Pastor-Satorras, and Vespignani 2002; Barrat, Barthelemy, and Vespignani 2008
Factory
@register_model(
field='biology/epidemiology',
tags=('network', 'mean field', 'ode', 'epidemic threshold', 'scale-free', 'final size'),
references=(
'Bertin 2021',
'Moreno, Pastor-Satorras, and Vespignani 2002',
'Barrat, Barthelemy, and Vespignani 2008',
),
defaults={'until': 100.0, 'dt': 0.5},
observables=lambda system: [
ObservableSpec('infected', _infected_share, group='network', reduce='sum'),
ObservableSpec('recovered', _recovered_share, group='network', reduce='sum'),
],
)
def heterogeneous_mean_field_sir(
exponent: float = 2.5,
k_min: int = 3,
k_max: int = 100,
*,
beta: float = 0.1,
mu: float = 1.0,
infected: float = 1e-4,
) -> System:
"""Return the SIR epidemic on an uncorrelated network of power-law degrees, in heterogeneous mean field.
The degrees run from `k_min` to `k_max` with `P(k) ~ k^-exponent`. An outbreak grows when
`beta / mu > <k> / (<k^2> - <k>)`, a threshold close to `<k> / <k^2>` on heterogeneous networks, which vanishes on
scale-free networks as `k_max` grows. Below it the epidemic dies out after infecting a share of the order of the
initial one, and above it a finite share of the population ends recovered. The hubs are infected first.
Args:
exponent:
The exponent of the degree distribution.
k_min:
The smallest degree.
k_max:
The largest degree.
beta:
The rate of transmission along one contact.
mu:
The recovery rate.
infected:
The initial infected density of every class.
Returns:
system:
The model, with the densities `network.s`, `network.i`, and `network.r` of the degree classes from `k_min`
to `k_max`, members of the group `network`, and the shares `infected` and `recovered` of the whole
population, the last one the final epidemic size at the end of the outbreak.
Examples:
>>> from skcomplex.models.biology.epidemiology import heterogeneous_mean_field_sir
>>> heterogeneous_mean_field_sir(k_min=2, k_max=5).elements[0].check().prepare().size_
4
"""
degrees = np.arange(k_min, k_max + 1, dtype=float)
weights = degrees**-exponent / np.sum(degrees**-exponent)
classes = GroupElement(
'network',
element=ImmunizingDegreeClassElement,
size=len(degrees),
values={'s': 1.0 - infected, 'i': infected, 'r': 0.0, 'beta': beta, 'mu': mu, 'k': degrees, 'P': weights},
)
interactions = [
CouplingInteraction(
'infection',
target='network',
read='i',
write='i',
function=_infection,
reduce='normalized-mean',
mode='group',
weights=_edge_weight,
include_self=True,
),
CouplingInteraction(
'contagion',
target='network',
read='s',
write='s',
function=_susceptible_loss,
reduce='normalized-mean',
mode='group',
weights=_edge_weight,
include_self=True,
),
TransferInteraction('recovery', source='network', target='network', state='i', into='r', amount=_recovery),
]
return System(elements=[classes], interactions=interactions)
Components
class ImmunizingDegreeClassElement(Element):
"""Susceptible, infected, and recovered densities of the nodes of one degree, in heterogeneous mean field.
A susceptible node of degree `k` is infected at rate `beta k theta` and an infected node recovers at rate `mu`, so
`d i_k / dt = -mu i_k + beta k s_k theta` and `d r_k / dt = mu i_k`. On an uncorrelated network, the probability
that a contact of a susceptible node is infected is `theta = sum_k (k - 1) P(k) i_k / <k>`, a total over the
degree classes, since an infected neighbour of degree `k` has at most `k - 1` other contacts. The classes are the
members of a group, and the infections and the recoveries are interactions of the group, which
`heterogeneous_mean_field_sir` adds.
Args:
label:
The label of the class or the group of classes.
visual:
How the element is drawn, or `None` for `default_visual_`.
s:
The initial susceptible density.
i:
The initial infected density.
r:
The initial recovered density.
beta:
The rate of transmission along one contact.
mu:
The recovery rate.
k:
The degree of the class.
P:
The probability `P(k)` of the degree.
Examples:
>>> from skcomplex.models.biology.epidemiology import ImmunizingDegreeClassElement
>>> GroupElement(
... 'network',
... element=ImmunizingDegreeClassElement,
... size=3,
... values={'k': [1.0, 2.0, 3.0], 'P': [0.5, 0.3, 0.2]},
... ).size
3
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
s: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.999,
i: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.001,
r: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
mu: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
P: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.s = s
self.i = i
self.r = r
self.beta = beta
self.mu = mu
self.k = k
self.P = P
class CouplingInteraction(NeighborInteraction):
"""Coupling of every member to other members, through a reduction over them of a function of both members.
Every member `i` gets the reduction over other members `j` of `function(namespace)`, where `namespace.member` reads
the states of `i`, `namespace.other` those of `j`, `namespace.edge` the member of the edge group that joins them, or
`None`, and the namespace also holds the parameters and the time `t`. The function is the difference `x_j - x_i` of
the state `read` by default. The members `j` are, by `mode`:
- `'edges'`: the neighbours of `i` along the edges of its lattice or network, or within the `neighborhood`.
- `'group'`: every other member of the group, or every member of the group `source` when it is given.
In the mode `'group'` the roles may also bind separate elements. A target that lists elements couples every element
to the other elements of the list, a single target element is a hub that reduces over every member of the group
`source`, and a source that is one element or a list of elements couples every member of the target to them. The
reduction over separate elements is written out term by term.
- `'next'`: the next member of the group in member order, the first after the last.
- `'hyperedges'`: the other members of every hyperedge around `i` of `size` members, or of every size, on a network
with hyperedges. The function then reads `namespace.member`, `namespace.others`, the other members of one
hyperedge in its order, and `namespace.members`, all its members, which reductions over them read, and the
reduction runs over the hyperedges around `i`. The default function is the sum over the other members of the
difference of `read`.
With `include_self`, the reduction also runs over `j = i`. Every term `f_ij` has a weight `w_ij`: on edges, the
state `weights` of the edge group, which a subclass binds with a second role, over the group, the entry `W_ij` of
the matrix `weights` or of the shaped parameter of the interaction that `weights` names, which a subclass declares,
and in every mode a callable `weights(namespace)` of the namespace of the term, such as an
indicator of the states of both members. The reductions are `'sum'`, `sum_j w_ij f_ij`, `'mean'`, the sum divided
by the number of members `j`, `'max'` and `'min'`, the largest and smallest `w_ij f_ij`, `'normalized-mean'`,
`sum_j w_ij f_ij / sum_j w_ij`, and `'circular-mean'`, the angle `arctan2(sum_j w_ij sin f_ij, sum_j w_ij cos
f_ij)`. The state `read` of `j` is read `delay` earlier, and the one of `i` shifted by the phase lag `lag`, when
they are not 0.
In the form `'rate'` the reduction adds to the rate of the state `write`, and in the form `'map'` it sets the next
value of `write` every step of discrete time to `combine(namespace, reduction)`, where `namespace.member` reads the
states of `i` and the namespace holds the parameters and the time, or to the reduction itself by default. The sum
and the mean of the rate on the edges of a lattice or network run once per edge along both directions, and every
other coupling runs once per member. On a `FieldElement`, a rate that is a constant times the difference of the
ends becomes that constant times the field's Laplacian times its squared spacing, and on a lattice the other end is
read across the field's boundary conditions, so beyond fixed edges it reads the fixed value.
In the form `'agent'` the target is a group of agents, and every agent runs the coupling in its own activation in
the phase `phase`. Every `RadiusSpec` of `radii` gives the vector `function(agent, other, model)` of every other
agent of the source group, or of the target group, within its radius, reduced by its sum or mean. Only the radii
of the highest priority that hold another agent count, and their reductions add up. The states `write` of the
agent form a vector `v`, which becomes `combine(agent, model, reduction, priority)`, or `v` plus the reduction by
default, where `reduction` and `priority` are `None` when no radius holds another agent. `clamp` then bounds the
angle the vector turns and its norm.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group of the members, or in the mode `'group'` the path of one element or the paths of
separate elements.
source:
The path of the group, the element, or the paths of the separate elements of the other members in the mode
`'group'`, or `None` for the members of the target.
read:
The name of the state that the default function, the delay, and the lag read.
write:
The name of the state that the coupling changes, or the names of the components of the vector of the agent
form.
function:
The term of one other member, a callable of the namespace, or `None` for the difference of `read`.
reduce:
`'sum'`, `'mean'`, `'max'`, `'min'`, `'normalized-mean'`, or `'circular-mean'`.
mode:
`'edges'`, `'group'`, `'next'`, or `'hyperedges'`.
size:
The number of members of the hyperedges of the mode `'hyperedges'`, or `None` for every size.
weights:
The name of the state of the edge group in the mode `'edges'`, a matrix with one row per member and one
column per other member or the name of a parameter of the interaction declared
`ParameterSpec(shape=('target', 'source'))` in the mode `'group'`, a callable of the namespace of a term in
any mode, or `None` for no weights.
include_self:
Whether the reduction also runs over the member itself.
lag:
The phase lag added to the state `read` of the member.
delay:
The delay of the state `read` of the other member, 0 for none.
form:
`'rate'`, `'map'`, or `'agent'`.
combine:
The next value in the form `'map'`, a callable of the namespace of the member and the reduction, or in the
form `'agent'` a callable of the agent, the model view, the reduction, and its priority, or `None` for the
default.
neighborhood:
The neighbourhood of the mode `'edges'`, such as `RadiusNeighborhood(1.0)` in a continuous space, or `None`
for the edges of the lattice or network.
layer:
The layer of the network whose edges the interaction runs on, or `None` for the graph of the network.
radii:
The terms of the agent form over the radii around an agent.
clamp:
The clamp of the turn and the norm of the vector of the agent form, or `None`.
phase:
The phase of the step the agent and the map forms run in, or `None` for the default phase.
Raises:
ValidationError:
A variant is unknown, the function or the combination is not callable, the weights, the source, the
reduction, or the neighbourhood do not suit the mode and the form, or the radii, the clamp, the phase, or
the written states do not suit the form, raised by `check`.
Examples:
>>> from types import SimpleNamespace
>>> from sage.all__sagemath_symbolics import sin, var
>>> from skcomplex.core.interactions import CouplingInteraction
>>> source, target = SimpleNamespace(theta=var('a')), SimpleNamespace(theta=var('b'))
>>> sine = lambda s: sin(s.other.theta - s.member.theta)
>>> coupling = CouplingInteraction('coupling', target='oscillators', read='theta', write='theta', function=sine)
>>> [law.value for law in coupling.check().prepare().run(SimpleNamespace(source=source, target=target)).laws_]
[sin(a - b), sin(-a + b)]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str | Sequence[str], RoleSpec()],
source: Annotated[str | Sequence[str] | None, RoleSpec(optional=True)] = None,
read: str = 'x',
write: str | Sequence[str] = 'x',
function: NamespaceFunction | None = None,
reduce: str = 'sum',
mode: str = 'edges',
size: int | None = None,
weights: Weights = None,
include_self: bool = False,
lag: Annotated[float, ParameterSpec()] = 0.0,
delay: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
form: str = 'rate',
combine: Combination | AgentCombination | None = None,
neighborhood: BaseNeighborhood | None = None,
layer: str | None = None,
radii: Sequence[RadiusSpec] = (),
clamp: ClampSpec | None = None,
phase: str | None = None,
) -> None:
super().__init__(label, layer=layer)
self.target = target
self.source = source
self.read = read
self.write = write
self.function = function
self.reduce = reduce
self.mode = mode
self.size = size
self.weights = weights
self.include_self = include_self
self.lag = lag
self.delay = delay
self.form = form
self.combine = combine
self.neighborhood = neighborhood
self.radii = radii
self.clamp = clamp
self.phase = phase
def _check_variants(self) -> None:
"""Raise an error for an unknown mode, reduction, or form, or a function or combination that is not callable."""
self._check_variant('mode', self.mode, MODES, (BaseSpec,))
for name, value, known in (
('reduce', self.reduce, REDUCTIONS),
('form', self.form, FORMS),
):
if value not in known:
error_msg = f'The {name} of coupling is {value!r}. Use one of {", ".join(known)}.'
raise ValidationError(error_msg)
for role, item in (('function', self.function), ('combine', self.combine)):
if item is not None and not callable(item):
error_msg = f'The {role} of coupling must be callable or None, but it is {item!r}.'
raise ValidationError(error_msg)
def _check_reduction(self) -> None:
"""Raise an error for a source outside the mode group, or a reduction that its mode cannot run."""
if self.source is not None and self.mode != 'group':
error_msg = f'The source {self.source!r} of coupling needs the mode group, but the mode is {self.mode!r}.'
raise ValidationError(error_msg)
if not isinstance(self.target, str) and self.mode != 'group':
error_msg = (
f'The coupling of the separate elements {self.target!r} needs the mode group, not {self.mode!r}.'
)
raise ValidationError(error_msg)
if self.mode == 'next' and self.reduce not in SUMS:
error_msg = f'The coupling to the next member reduces by sum or mean, not {self.reduce!r}.'
raise ValidationError(error_msg)
if self.reduce in EXTREMES and self.mode == 'group' and self._excludes_self() and isinstance(self.target, str):
error_msg = (
f'The {self.reduce} of coupling over its own group runs over every member. Pass include_self=True or '
'a separate source group.'
)
raise ValidationError(error_msg)
if self.reduce in EXTREMES and (self._has_matrix() or self._names_matrix()):
error_msg = f'The {self.reduce} of coupling takes no matrix of weights. Give the weights as a callable.'
raise ValidationError(error_msg)
def _check_structure(self) -> None:
"""Raise an error for weights or a neighbourhood that do not suit the mode and the form."""
if self.neighborhood is not None and (
self.mode != 'edges' or not isinstance(self.neighborhood, BaseNeighborhood)
):
error_msg = f'The neighbourhood {self.neighborhood!r} of coupling needs a neighbourhood in the mode edges.'
raise ValidationError(error_msg)
if isinstance(self.weights, str) and not self._is_per_edge() and not self._names_matrix():
error_msg = (
f'The weights {self.weights!r} of coupling name a state of the edges, which only the sum or mean rate '
'on the edges of a lattice or network reads, or a shaped parameter in the mode group. Give a matrix '
'in the mode group, or a callable.'
)
raise ValidationError(error_msg)
if self._names_matrix() and (self._excludes_self() or not isinstance(self.target, str) or self._is_separate()):
error_msg = (
f'The weights {self.weights!r} of coupling name a shaped parameter, whose matrix product runs over '
'every member of a group. Pass include_self=True or a separate source group, and bind groups.'
)
raise ValidationError(error_msg)
if self._has_matrix() and self.mode != 'group':
error_msg = 'A matrix of weights of coupling needs the mode group. Name an edge state on the edges.'
raise ValidationError(error_msg)
def _check_hyperedges(self) -> None:
"""Raise an error for a size outside the mode hyperedges, or a setting the mode hyperedges does not take."""
if self.size is not None and (
self.mode != 'hyperedges'
or isinstance(self.size, bool)
or not isinstance(self.size, int)
or self.size < MIN_SIZE
):
error_msg = (
f'The size {self.size!r} of coupling needs the mode hyperedges and a whole number of 2 or more, but '
f'the mode is {self.mode!r}.'
)
raise ValidationError(error_msg)
if self.mode == 'hyperedges' and (self.include_self or self._is_nonzero('delay') or self.form == 'agent'):
error_msg = (
'The coupling over hyperedges takes no include_self, no delay, and no agent form. Couple the member '
'to itself with another interaction.'
)
raise ValidationError(error_msg)
def _check_agent(self) -> None:
"""Raise an error for radii, a clamp, or several written states outside the agent form, or no radii."""
if self.form != 'agent':
if self.radii or self.clamp is not None or not isinstance(self.write, str):
error_msg = (
f'The coupling {self.label!r} takes radii, a clamp, and several written states in the agent form '
f'only, but its form is {self.form!r}.'
)
raise ValidationError(error_msg)
if self.phase is not None and self.form != 'map':
error_msg = (
f'The coupling {self.label!r} takes a phase in the agent and the map forms, but its form is '
f'{self.form!r}.'
)
raise ValidationError(error_msg)
return
if not self.radii or not isinstance(self.target, str):
error_msg = f'The agent form of coupling {self.label!r} needs one group of agents and at least one radius.'
raise ValidationError(error_msg)
for spec in self.radii:
spec.check()
if self.clamp is not None:
self.clamp.check()
def _check(self) -> None:
"""Raise an error for unknown variants, functions that cannot be called, or weights that do not suit the mode.
Raises:
ValidationError:
A variant is unknown, the function or the combination is not callable, or the weights, the source, the
reduction, or the neighbourhood do not suit the mode and the form.
"""
super()._check()
self._check_variants()
self._check_reduction()
self._check_structure()
self._check_hyperedges()
self._check_agent()
def _get_written(self) -> str:
"""Return the name of the one state that the symbolic forms write."""
return cast(str, self.write)
def _has_matrix(self) -> bool:
"""Return whether the weights are a matrix."""
return self.weights is not None and not isinstance(self.weights, str) and not callable(self.weights)
def _names_matrix(self) -> bool:
"""Return whether the weights name a shaped parameter of the interaction in the mode group."""
return isinstance(self.weights, str) and self.mode == 'group'
def _is_separate(self) -> bool:
"""Return whether the source lists separate elements."""
return self.source is not None and not isinstance(self.source, str)
def _excludes_self(self) -> bool:
"""Return whether a reduction over the own group leaves the member out."""
return self.source is None and not self.include_self
def _is_nonzero(self, name: str) -> bool:
"""Return whether a parameter, read under the name `_get_parameter` returns, has a value other than 0."""
return bool(getattr(self, self._get_parameter(name)) != 0)
def _is_per_edge(self) -> bool:
"""Return whether the coupling is a rate written once per edge of a lattice or network."""
return self.mode == 'edges' and self.form == 'rate' and self.neighborhood is None and self.reduce in SUMS
def _parameters(self, namespace: SimpleNamespace) -> dict[str, object]:
"""Return the parameters and the time of a namespace, without the ends and the edge."""
ends = {self._get_role(end) for end in (*ENDS, 'edge')}
return {name: value for name, value in vars(namespace).items() if name not in ends}
def _get_term(
self,
namespace: SimpleNamespace,
member: SimpleNamespace,
other: SimpleNamespace,
edge: SimpleNamespace | None = None,
) -> Target:
"""Return the term of one other member, with the delay and the lag applied to the read state."""
member_values, other_values = dict(vars(member)), dict(vars(other))
if self._is_nonzero('delay'):
other_values[self.read] = delayed(other_values[self.read], self._get_parameter_symbol('delay', namespace))
if self._is_nonzero('lag'):
member_values[self.read] = member_values[self.read] + self._get_parameter_symbol('lag', namespace)
if self.function is None:
return other_values[self.read] - member_values[self.read]
return self.function(
SimpleNamespace(
**self._parameters(namespace),
member=SimpleNamespace(**member_values),
other=SimpleNamespace(**other_values),
edge=edge,
),
)
def _get_weight(
self,
namespace: SimpleNamespace,
member: SimpleNamespace,
other: SimpleNamespace,
edge: SimpleNamespace | None = None,
) -> Target:
"""Return the weight of the term of one other member from callable weights, or 1 for other weights."""
if not callable(self.weights):
return SR(1)
return self.weights(SimpleNamespace(**self._parameters(namespace), member=member, other=other, edge=edge))
def _get_matrix(self, namespace: SimpleNamespace) -> npt.ArrayLike | scipy.sparse.sparray | Expression | None:
"""Return the matrix of weights, with a zero diagonal unless the member is included, or `None`.
Weights that name a shaped parameter give the symbol of that parameter.
"""
if self._names_matrix():
name = cast(str, self.weights)
if not hasattr(namespace, name):
error_msg = (
f'The weights {name!r} of coupling {self.label!r} name no parameter of the interaction. Declare '
f"it as `ParameterSpec(shape=('target', 'source'))`."
)
raise ValidationError(error_msg)
return cast(Expression, getattr(namespace, name))
if not self._has_matrix():
return None
weights = cast(npt.ArrayLike | scipy.sparse.sparray, self.weights)
if not self._excludes_self():
return weights
if scipy.sparse.issparse(weights):
matrix = scipy.sparse.csr_array(weights, dtype=float)
matrix.setdiag(0)
matrix.eliminate_zeros()
return matrix
matrix = np.array(weights, dtype=float)
np.fill_diagonal(matrix, 0.0)
return matrix
@staticmethod
def _map(value: Target, operation: Callable[[Expression], Expression]) -> Target:
"""Return an operation applied to an expression, or to every component of a vector."""
if isinstance(value, FreeModuleElement):
return value.parent()([operation(SR(item)) for item in value])
return operation(SR(value))
@staticmethod
def _shifted(member: SimpleNamespace) -> SimpleNamespace:
"""Return the namespace of the next member of the group in member order."""
return SimpleNamespace(
**{
name: CouplingInteraction._map(value, lambda item: shift_member(item, 1))
for name, value in vars(member).items()
if isinstance(value, Expression | FreeModuleElement)
},
)
def _get_reduced(self, total_of: Callable[[Part], Target], number: Target) -> Target:
"""Return the reduction from sums of parts of the weighted terms and the number of members, by reduction."""
if self.reduce == 'sum':
return total_of(lambda value, weight: weight * value)
if self.reduce == 'mean':
return total_of(lambda value, weight: weight * value) / number
if self.reduce == 'normalized-mean':
return total_of(lambda value, weight: weight * value) / total_of(lambda value, weight: weight)
sines = total_of(lambda value, weight: weight * sin(value))
return arctan2(sines, total_of(lambda value, weight: weight * cos(value)))
def _get_group_total(self, namespace: SimpleNamespace, part: Part) -> Target:
"""Return the sum of a part of the weighted term over the other members of the group."""
source, target = getattr(namespace, self._get_role('source')), self._get_member(namespace, 'target')
matrix = self._get_matrix(namespace)
value, weight = self._get_term(namespace, target, source), self._get_weight(namespace, target, source)
reduction = self._map(part(value, weight), lambda item: total(item, over=source, weights=matrix))
if self._excludes_self() and matrix is None:
reduction = reduction - part(
self._get_term(namespace, target, target),
self._get_weight(namespace, target, target),
)
return reduction
@staticmethod
def _get_extreme(values: list[Target], extreme: Callable[..., Expression]) -> Target:
"""Return the largest or smallest of several values, component by component for vectors."""
if len(values) == 1:
return values[0]
first = values[0]
if isinstance(first, FreeModuleElement):
columns = zip(*[list(value) for value in values], strict=True)
return first.parent()([extreme(*[SR(item) for item in column]) for column in columns])
return extreme(*[SR(value) for value in values])
@staticmethod
def _get_sum(pairs: list[tuple[Target, Target]], part: Part) -> Target:
"""Return the sum of a part of the weighted terms of separate elements."""
parts = [part(value, weight) for value, weight in pairs]
return sum(parts[1:], parts[0])
def _get_separate_reduction(self, namespace: SimpleNamespace) -> Target:
"""Return the reduction over the separate elements of `namespace.source`, written out term by term."""
if self._has_matrix():
error_msg = (
'A matrix of weights of coupling needs groups. Give the weights of separate elements as a callable.'
)
raise ValidationError(error_msg)
target = self._get_member(namespace, 'target')
others = [*getattr(namespace, self._get_role('source')), *([target] if self.include_self else [])]
pairs = [
(self._get_term(namespace, target, other), self._get_weight(namespace, target, other)) for other in others
]
if not pairs:
error_msg = f'The coupling {self.label!r} of separate elements has no other element. Bind two or more.'
raise ValidationError(error_msg)
if self.reduce in EXTREMES:
extreme = maximum if self.reduce == 'max' else minimum
return self._get_extreme([value * weight for value, weight in pairs], extreme)
return self._get_reduced(lambda part: self._get_sum(pairs, part), SR(len(pairs)))
def _get_group_reduction(self, namespace: SimpleNamespace) -> Target:
"""Return the reduction over every other member of the group, or every member of the source group."""
source, target = getattr(namespace, self._get_role('source')), self._get_member(namespace, 'target')
if isinstance(source, tuple):
return self._get_separate_reduction(namespace)
if self.reduce in EXTREMES:
extreme = maximum if self.reduce == 'max' else minimum
value = self._get_term(namespace, target, source) * self._get_weight(namespace, target, source)
return self._map(value, lambda item: extreme(item, over=source))
number = count(over=source) - (1 if self._excludes_self() else 0)
return self._get_reduced(lambda part: self._get_group_total(namespace, part), number)
def _get_next_reduction(self, namespace: SimpleNamespace) -> Target:
"""Return the reduction over the next member of the group."""
target = self._get_member(namespace, 'target')
shifted = self._shifted(target)
reduction = self._get_term(namespace, target, shifted) * self._get_weight(namespace, target, shifted)
if self.include_self:
own = self._get_term(namespace, target, target) * self._get_weight(namespace, target, target)
reduction = reduction + own
return reduction / 2 if self.reduce == 'mean' and self.include_self else reduction
def _get_neighbor_total(self, namespace: SimpleNamespace, part: Part) -> Target:
"""Return the sum of a part of the weighted term over the neighbours and the member itself."""
source, target = getattr(namespace, self._get_role('source')), self._get_member(namespace, 'target')
value, weight = self._get_term(namespace, target, source), self._get_weight(namespace, target, source)
reduction = self._map(part(value, weight), lambda item: neighbors(item, over=self.neighborhood))
return reduction + part(self._get_term(namespace, target, target), self._get_weight(namespace, target, target))
def _get_neighbor_reduction(self, namespace: SimpleNamespace) -> Target:
"""Return the reduction over the neighbours of the member."""
source, target = getattr(namespace, self._get_role('source')), self._get_member(namespace, 'target')
term = self._get_term(namespace, target, source)
weight = self._get_weight(namespace, target, source) if callable(self.weights) else None
if not self.include_self:
return self._map(
term,
lambda item: neighbors(item, over=self.neighborhood, reduce=self.reduce, weight=weight),
)
if self.reduce in EXTREMES:
extreme = maximum if self.reduce == 'max' else minimum
others = neighbors(term, over=self.neighborhood, reduce=self.reduce, weight=weight)
own = self._get_term(namespace, target, target) * self._get_weight(namespace, target, target)
return extreme(others, own)
number = neighbors(SR(1), over=self.neighborhood, reduce='count') + 1
return self._get_reduced(lambda part: self._get_neighbor_total(namespace, part), number)
def _get_hyperedge_term(self, namespace: SimpleNamespace, around: SimpleNamespace) -> Target:
"""Return the term of one hyperedge around a member, with the lag applied to the read state of the member."""
member = around.member
if self._is_nonzero('lag'):
lag = self._get_parameter_symbol('lag', namespace)
member = SimpleNamespace(**{**vars(member), self.read: getattr(member, self.read) + lag})
if self.function is None:
members = around.members
own = getattr(around.member, self.read)
others = total(getattr(members, self.read), over=members) - own
return others - (count(over=members) - 1) * getattr(member, self.read)
return self.function(SimpleNamespace(**self._parameters(namespace), **{**vars(around), 'member': member}))
def _get_hyperedge_reduction(self, namespace: SimpleNamespace) -> Target:
"""Return the reduction over the hyperedges around the member.
Args:
namespace:
The namespace of the states and parameters of the interaction.
Returns:
reduction:
The reduction over the hyperedges around the member.
Raises:
ValidationError:
The members live on no network with hyperedges.
"""
target = self._get_member(namespace, 'target')
others, members = getattr(target, 'others', None), getattr(target, 'members', None)
if others is None or members is None:
error_msg = (
f'The coupling {self.label!r} over hyperedges binds {self.target!r}, whose members live on no network '
'with hyperedges. Give the group a network space with hyperedges.'
)
raise ValidationError(error_msg)
if self.size is not None:
others = others[: self.size - 1]
around = SimpleNamespace(member=target, others=others, members=members)
term = self._get_hyperedge_term(namespace, around)
weight = None
if callable(self.weights):
weight = self.weights(SimpleNamespace(**self._parameters(namespace), **vars(around)))
return self._map(term, lambda item: hyperedges(item, reduce=self.reduce, size=self.size, weight=weight))
def _get_reduction(self, namespace: SimpleNamespace) -> Target:
"""Return the reduction that one member gets, by mode."""
if self.mode == 'hyperedges':
return self._get_hyperedge_reduction(namespace)
if self.mode == 'group':
return self._get_group_reduction(namespace)
if self.mode == 'next':
return self._get_next_reduction(namespace)
if self.mode == 'edges':
return self._get_neighbor_reduction(namespace)
raise self._get_variant_error(self.mode, 'mode', '_get_reduction')
def _get_edge_term(self, namespace: SimpleNamespace, member: SimpleNamespace, other: SimpleNamespace) -> Target:
"""Return the term one end of an edge gets from the other end, weighted and divided for a mean."""
edge = getattr(namespace, self._get_role('edge'), None)
term = self._get_term(namespace, member, other, edge)
if isinstance(self.weights, str):
term = term * getattr(edge, self.weights)
elif callable(self.weights):
term = term * self._get_weight(namespace, member, other, edge)
if self.reduce == 'mean':
term = term / (member.degree + int(self.include_self))
return term
def _get_edge_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the rates of both ends of an edge, with the rate of the member itself when it is included."""
source, target = getattr(namespace, self._get_role('source')), self._get_member(namespace, 'target')
laws: list[BaseLaw] = [
RateLaw(getattr(target, self._get_written()), self._get_edge_term(namespace, target, source)),
RateLaw(getattr(source, self._get_written()), self._get_edge_term(namespace, source, target)),
]
if self.include_self:
own = self._get_term(namespace, target, target) * self._get_weight(namespace, target, target)
divisor = target.degree + 1 if self.reduce == 'mean' else 1
laws.append(RateLaw(getattr(target, self._get_written()), own / divisor))
return laws
@staticmethod
def _get_bound(agent: AgentViewProtocol, bound: float | str) -> float:
"""Return a bound given as a number or as the name of a parameter of the agent."""
return float(cast(float, getattr(agent, bound))) if isinstance(bound, str) else float(bound)
def _get_agent_reduction(
self,
agent: AgentViewProtocol,
model: ModelViewProtocol,
) -> tuple[Array | None, int | None]:
"""Return the sum of the reduced terms of the radii of the highest priority that hold another agent."""
groups = set(model.groups(self._get_role('source' if self.source is not None else 'target')))
found: dict[float, list[AgentViewProtocol]] = {}
levels: dict[int, list[Array]] = {}
for spec in self.radii:
radius = self._get_bound(agent, spec.radius)
if radius not in found:
found[radius] = [other for other in agent.neighbors(radius=radius) if other.group in groups]
others = found[radius]
if not others:
continue
terms = np.array([np.asarray(spec.function(agent, other, model), dtype=float) for other in others])
levels.setdefault(spec.priority, []).append(
terms.mean(axis=0) if spec.reduce == 'mean' else terms.sum(axis=0),
)
if not levels:
return None, None
priority = max(levels)
return np.sum(levels[priority], axis=0), priority
@staticmethod
def _turn(current: Array, value: Array, max_angle: float) -> Array:
"""Return a vector turned from the current direction toward its own by at most an angle, keeping its norm."""
norm, size = float(np.linalg.norm(value)), float(np.linalg.norm(current))
if norm == 0 or size == 0:
return value
direction, facing = value / norm, current / size
angle = float(np.arccos(np.clip(facing @ direction, -1.0, 1.0)))
if angle <= max_angle:
return value
if angle == np.pi:
direction = np.eye(len(facing))[int(np.argmin(np.abs(facing)))]
normal = direction - (facing @ direction) * facing
normal = normal / np.linalg.norm(normal)
return norm * (np.cos(max_angle) * facing + np.sin(max_angle) * normal)
def _get_clamped(self, agent: AgentViewProtocol, current: Array, value: Array) -> Array:
"""Return a new vector with its turn and its norm clamped."""
clamp = self.clamp
if clamp is None:
return value
if clamp.max_angle is not None:
value = self._turn(current, value, self._get_bound(agent, clamp.max_angle))
norm = float(np.linalg.norm(value))
if norm == 0:
return value
low = 0.0 if clamp.min_norm is None else self._get_bound(agent, clamp.min_norm)
high = np.inf if clamp.max_norm is None else self._get_bound(agent, clamp.max_norm)
return value * min(max(norm, low), high) / norm
def _update_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Set the written states of one agent to the clamped combination of the reduction over its radii."""
if not agent.alive:
return
names = [self.write] if isinstance(self.write, str) else list(self.write)
current = np.array([float(cast(float, getattr(agent, name))) for name in names])
reduction, priority = self._get_agent_reduction(agent, model)
if self.combine is not None:
combine = cast(AgentCombination, self.combine)
value = np.asarray(combine(agent, model, reduction, priority), dtype=float)
else:
value = current if reduction is None else current + reduction
for name, item in zip(names, self._get_clamped(agent, current, value), strict=True):
setattr(agent, name, float(item))
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the rates of both ends of every edge, or the rate or next value of every member.
Args:
namespace:
The namespace of an edge, with the ends `source` and `target`.
Returns:
laws:
The laws of the coupling.
"""
if self.form == 'agent':
return [RuleLaw(self._update_member, phase=self.phase, members=self._get_role('target'))]
if self._is_per_edge():
return self._get_edge_laws(namespace)
target = self._get_member(namespace, 'target')
reduction = self._get_reduction(namespace)
if self.form == 'rate':
return [RateLaw(getattr(target, self._get_written()), reduction)]
state = getattr(target, self._get_written())
if self.combine is None:
return [UpdateLaw(state, reduction, phase=self.phase)]
member = SimpleNamespace(**self._parameters(namespace), member=target)
return [UpdateLaw(state, cast(Combination, self.combine)(member, reduction), phase=self.phase)]
class TransferInteraction(Interaction):
"""Transfer of amounts of states from the members of a source role to the members of a target role.
A leg moves an amount out of a state of the source member, or out of `Source()`, into a state of the target member,
or into `Sink()`. The transfer has a first leg, `state` into `into` by `amount`, and any further `legs`, each
`(state, into, amount)` or `(state, into, amount, flow)`. The amounts of all legs are read before any of them
moves, so a leg sees the states as they were before the transfer. A guard skips the transfer where it does not
hold.
The symbolic form reports the amount of every leg as a flow. A leg named by its `flow` is reported under that
name: under the element or group, as `<path>.<flow>`, when the source and the target are one element or group,
since the flow then moves a quantity within it like a flow of its own, and under the interaction, as
`<label>.<flow>`, otherwise, as it is for a subclass that names its roles otherwise. A leg without a name is
reported under the interaction as `<label>.flow_<index>`.
The form chooses how the transfer lowers:
- `'symbolic'` lowers to `TransferLaw`s between the source and the target of one element, or member `i` of both
groups. The amounts, the guard, and the rates are numbers or functions of the namespace, and the roles are read
as `namespace.source` and `namespace.target`. With a `rate`, the amount is the size of a jump at the rate, and
with a `reverse_rate` jumps also move the amounts back from the target to the source at that rate, so a leg from
`Source()` returns its amount to the environment. With an `offset`, the source is a group and every member
transfers along an edge of the group order, to the member `offset` places along it, which `namespace.target`
reads. Past the ends of a `'fixed'` boundary no member is there, and the edge moves nothing. A draw in an amount
is one draw per edge, which both of its ends read, so the amounts of a step move in parallel. A leg end
`'<role>.<state>'` names the state of the element or group another role binds, so one transfer moves amounts
between several holders, such as the payments of a circuit of sectors, and a subclass that declares parameters
reads them in its amounts. A subclass that names its roles otherwise returns their names from `_get_role`, or
the namespaces of the source and the target of every pair from `_get_holders`.
- `'agent'` lowers to an agent rule, whose amounts and guard are numbers or functions of the source agent, the
target agent, and the model view. The pairing `'members'` pairs member `i` of both groups, `'random'` pairs every
source member with a random target member, `'pairs'` draws disjoint random pairs, `'neighbors'` pairs every
source member with the target members in the cells or nodes around it, `'neighbor-pairs'` draws disjoint random
pairs of neighbours, every source member in random order taking a random unpaired target member around it, and
`'cell'` pairs the source and target members that stand in one cell or node. With a `fallback` group, the members
of that group are the targets of a step in which no target member lives. With the pairing `'site'` the target
role binds a field, and every source member moves its amounts into the site of the field at its cell or node, so
a negative amount takes from the site. The functions of this pairing take the cell of the source agent in place
of a target agent. With `remove`, the source or the target member of every pair that transfers is removed, as an
agent eaten whole. The named pairings other than `'members'` and `'site'` are shorthands of a `PartnerSpec`, which
pairs every source member with partners within a radius, by score, or with a fallback. A `settlement` receives
both partners and the amounts of the legs, read once, and changes the states in place of the legs, so one draw
credits both ends. A subclass that names its roles otherwise returns their names from `_get_role`.
Under the schedule `'random-sequential'` the pairs act in random order. Under `'member'` every member of the
`actor` role acts in its own activation, together with the rules of its element and the other member
interactions, in the order the system lists them: it transfers with its partners, the members of the other role
that the pairing `'random'`, `'neighbors'`, or `'cell'` gives it, or with `partners` of them drawn at random, and
with itself under the pairing `'members'` within one group.
A leg end of the agent form may name `'<reference>.<state>'`, the state of the member that a reference state of
the agent names, such as `'bank.deposits'` for the deposits of a customer's bank, declared
`StateSpec(reference='banks')`.
Args:
label:
The label of the interaction, unique among the parts of its system.
source:
The path of the element or group the amounts leave.
target:
The path of the element or group the amounts enter, the source again for a transfer within members.
state:
The name of the state of the source the first leg takes from, or `Source()` for the environment.
into:
The name of the state of the target the first leg adds to, `Sink()` for the environment, or `None` for
`state`.
amount:
The amount of the first leg, a number or a function.
flow:
The name the symbolic form reports the first leg under, or the empty string for `flow_0`.
pairing:
How members are paired, one of `'members'`, `'random'`, `'pairs'`, `'neighbors'`, `'neighbor-pairs'`,
`'cell'`, or `'site'`, the agent form alone taking the others than `'members'`, or, in the agent form, a
`PartnerSpec`, or a function `pairing(sources, targets, rng, capacity)` of the source and target agents
that returns their pairs, whose capacity is `partners`, or 1 without it. Under the schedule `'member'` the
spec or the function receives the actor alone as its first sequence and the actor's partners are the second
items of its pairs.
legs:
Further legs, each a tuple of a state or `Source()`, a state, `Sink()`, or `None`, an amount, and
optionally the name of its flow.
fallback:
The path of the group whose members the agent form pairs with instead of the target members in a step in
which no target member lives, or `None` for no fallback.
guard:
A function that returns whether a pair transfers, a condition of the namespace or a boolean of the agents,
or `None` for every pair.
rate:
The rate of jumps of the symbolic form, a number or a function of the namespace, or `None` for a
continuous flow.
reverse_rate:
The rate of the jumps of the symbolic form that move the amounts back from the target to the source, a
number or a function of the namespace, or `None` for no reverse jumps.
offset:
The number of places along the group order from every source member to its target member, or 0 to pair
the members of the roles.
boundary:
Past the ends of the group order, `'periodic'` closes a ring and `'fixed'` holds no member.
form:
The lowering, `'symbolic'` to transfer laws or `'agent'` to an agent rule.
schedule:
The order of the agent form, `'random-sequential'` over the pairs, or `'member'` over the actors.
actor:
The role whose members act under the schedule `'member'`, `'source'` or `'target'`.
partners:
The number of partners an actor draws at random under the schedule `'member'`, or `None` for all of them.
remove:
The member of a pair that transfers that is removed, `'source'`, `'target'`, or `None`.
phase:
The phase the rule of the agent form runs in, or `None` for the default phase.
settlement:
A function `settlement(source, target, amounts, model)` that settles a pair of agents from the amounts of
the legs, read once, in place of moving them, or `None` to move them.
Raises:
ValidationError:
The form, the pairing, the schedule, the actor, the removal, or the boundary is unknown, a leg is
malformed, two legs share a flow name, the symbolic form has an agent pairing, a fallback, or a settlement,
the agent form a rate, an offset, or a flow name, a reverse rate has no rate, an offset is not a whole
number or pairs two groups, a leg end of the symbolic form names no role, the settlement cannot be called
or pairs with a site, or the schedule `'member'` has a pairing it cannot draw, raised by `check`.
Examples:
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import TransferInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> from typing import Annotated
>>> class TankElement(Element):
... def __init__(self, label, visual=None, *, level: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label, visual=visual)
... self.level = level
>>> flow = TransferInteraction(
... 'flow', source='upper', target='lower', state='level', amount=lambda s: s.source.level / 2
... )
>>> system = System(elements=[TankElement('upper', level=4.0), TankElement('lower')], interactions=[flow])
>>> trajectory = Simulation(system).check().prepare().run(until=1.0).trajectory_
>>> round(float(trajectory.values('upper.level')[-1] + trajectory.values('lower.level')[-1]), 6)
4.0
"""
def __init__(
self,
label: str,
*,
source: Annotated[str, RoleSpec()],
target: Annotated[str, RoleSpec()],
state: str | Source,
into: str | Sink | None = None,
amount: Amount,
flow: str = '',
pairing: str | Pairing = 'members',
legs: Sequence[Leg] = (),
fallback: Annotated[str | None, RoleSpec()] = None,
guard: Callable[..., object] | None = None,
rate: Amount | None = None,
reverse_rate: Amount | None = None,
offset: int = 0,
boundary: str = 'periodic',
form: str = 'symbolic',
schedule: str = 'random-sequential',
actor: str = 'source',
partners: int | None = None,
remove: str | None = None,
phase: str | None = None,
settlement: Settlement | None = None,
) -> None:
super().__init__(label)
self.source = source
self.target = target
self.state = state
self.into = into
self.amount = amount
self.flow = flow
self.pairing = pairing
self.legs = legs
self.fallback = fallback
self.guard = guard
self.rate = rate
self.reverse_rate = reverse_rate
self.offset = offset
self.boundary = boundary
self.form = form
self.schedule = schedule
self.actor = actor
self.partners = partners
self.remove = remove
self.phase = phase
self.settlement = settlement
def _get_legs(self) -> list[NamedLeg]:
"""Return every leg with the state it adds to resolved and the name of its flow, empty when it has none."""
legs: list[NamedLeg] = []
for state, into, amount, *flow in ((self.state, self.into, self.amount, self.flow), *self.legs):
name = flow[0] if flow else ''
legs.append((state, state if into is None and isinstance(state, str) else into, amount, name))
return legs
def _check_leg(self, state: object, into: object, amount: object, flow: object) -> None:
"""Raise an error for a leg with unusable ends, amount, or flow name, or a symbolic end that names no role."""
if not isinstance(state, str | Source) or not isinstance(into, str | Sink):
error_msg = (
f'A leg of {self.label!r} moves {state!r} into {into!r}. Give a state name or `Source()` it leaves, '
'and a state name or `Sink()` it enters.'
)
raise ValidationError(error_msg)
if isinstance(state, Source) and isinstance(into, Sink):
error_msg = f'A leg of {self.label!r} moves from `Source()` into `Sink()`. Name a state of one end.'
raise ValidationError(error_msg)
if not callable(amount) and not isinstance(amount, int | float):
error_msg = f'The amount {amount!r} of {self.label!r} is neither a number nor a function.'
raise ValidationError(error_msg)
if not isinstance(flow, str):
error_msg = f'A leg of {self.label!r} names its flow {flow!r}. Give a string, or the empty string.'
raise ValidationError(error_msg)
roles = [name for name, spec in self._get_specs().items() if isinstance(spec, RoleSpec)]
referenced = [end for end in (state, into) if isinstance(end, str) and '.' in end]
unknown = [end for end in referenced if end.partition('.')[0] not in roles]
if unknown and self.form == 'symbolic':
error_msg = (
f'{self.label!r} moves {unknown[0]!r}, but its symbolic form names the states of roles, one of '
f'{", ".join(roles)}. Name a role, or use `form="agent"` for a referenced member.'
)
raise ValidationError(error_msg)
def _check_flows(self) -> None:
"""Raise an error for flow names that two legs share, or that the agent form, which reports no flows, takes."""
names = [leg[3] for leg in self._get_legs() if leg[3]]
if len(set(names)) < len(names):
error_msg = f'Two legs of {self.label!r} share a flow name among {names!r}. Give every flow its own name.'
raise ValidationError(error_msg)
if names and self.form == 'agent':
error_msg = f'The agent form of {self.label!r} reports no flows, but it names {names!r}. Remove the names.'
raise ValidationError(error_msg)
def _check_space(self, space: BaseSpace | None) -> None:
"""Raise an error for a pairing of cells, nodes, or sites in a continuous space."""
if isinstance(self.pairing, str) and self.pairing in GRID_PAIRINGS:
alternative = "a radius pairing, `PartnerSpec(candidates=<radius>)`, or pairing='random'"
self._check_grid('pairing', self.pairing, space, alternative)
self._check_partners(self.pairing, space)
def _check(self) -> None:
"""Raise an error for an unknown form or pairing, a malformed leg, or an option the form does not take.
Raises:
ValidationError:
The form, the pairing, the schedule, the actor, the removal, or the boundary is unknown, a leg is
malformed, two legs share a flow name, the symbolic form has an agent pairing, a fallback, or a leg end
that names no role, the agent form has a rate, an offset, or a flow name, a reverse rate has no rate,
an offset is not a whole number or pairs two groups, or the schedule `'member'` has a pairing it
cannot draw.
"""
super()._check()
if self.form not in FORMS:
error_msg = f'Unknown form {self.form!r} of {self.label!r}. Use one of {", ".join(FORMS)}.'
raise ValidationError(error_msg)
if self.form == 'symbolic' and self.pairing != 'members':
error_msg = f'The symbolic form of {self.label!r} pairs members. Use `form="agent"` for {self.pairing!r}.'
raise ValidationError(error_msg)
self._check_variant('pairing', self.pairing, PAIRINGS, (PartnerSpec,), functions=True)
if isinstance(self.pairing, PartnerSpec):
self.pairing.check()
self._check_settlement()
if self.form == 'symbolic' and self.fallback is not None:
error_msg = f'The symbolic form of {self.label!r} pairs members. Use `form="agent"` for a fallback group.'
raise ValidationError(error_msg)
for leg in self.legs:
if not isinstance(leg, tuple) or len(leg) not in LEG_LENGTHS:
error_msg = f'A leg of {self.label!r} is {leg!r}. Give a tuple `(state, into, amount[, flow])`.'
raise ValidationError(error_msg)
for state, into, amount, flow in self._get_legs():
self._check_leg(state, into, amount, flow)
self._check_flows()
if self.rate is not None and self.form == 'agent':
error_msg = f'{self.label!r} has a rate, which only the symbolic form takes. Give `rate=None`.'
raise ValidationError(error_msg)
if self.guard is not None and not callable(self.guard):
error_msg = f'The guard of {self.label!r} is {self.guard!r}. Give a function, or None.'
raise ValidationError(error_msg)
self._check_edges()
self._check_schedule()
def _check_settlement(self) -> None:
"""Raise an error for a settlement that cannot be called, or that the form or the pairing does not take."""
if self.settlement is None:
return
if not callable(self.settlement):
error_msg = f'The settlement of {self.label!r} is {self.settlement!r}. Give a function, or None.'
raise ValidationError(error_msg)
if self.form != 'agent' or self.pairing == 'site':
error_msg = f'{self.label!r} settles a pair of agents. Use `form="agent"` and a pairing other than `site`.'
raise ValidationError(error_msg)
def _check_edges(self) -> None:
"""Raise an error for an unusable reverse rate, offset, or boundary."""
if self.reverse_rate is not None and self.rate is None:
error_msg = f'{self.label!r} has a reverse rate but no rate. Give the rate of the forward jumps.'
raise ValidationError(error_msg)
reverse = self.reverse_rate
if reverse is not None and not callable(reverse) and not isinstance(reverse, int | float):
error_msg = f'The reverse rate {reverse!r} of {self.label!r} is neither a number nor a function.'
raise ValidationError(error_msg)
if isinstance(self.offset, bool) or not isinstance(self.offset, int):
error_msg = f'The offset of {self.label!r} is {self.offset!r}. Give a whole number of members.'
raise ValidationError(error_msg)
if self.boundary not in BOUNDARIES:
error_msg = f'Unknown boundary {self.boundary!r} of {self.label!r}. Use one of {", ".join(BOUNDARIES)}.'
raise ValidationError(error_msg)
if self.offset and (self.form != 'symbolic' or self.source != self.target):
error_msg = (
f'{self.label!r} transfers along the order of one group by an offset. Give the group as both source '
'and target, and `form="symbolic"`.'
)
raise ValidationError(error_msg)
def _check_schedule(self) -> None:
"""Raise an error for an unknown schedule, actor, or removal, or partners the schedule cannot draw."""
for name, value, variants in (
('schedule', self.schedule, SCHEDULES),
('actor', self.actor, ACTORS),
('removal', self.remove, (*REMOVALS, None)),
):
if value not in variants:
error_msg = f'Unknown {name} {value!r} of {self.label!r}. Use one of {", ".join(map(str, variants))}.'
raise ValidationError(error_msg)
own = self.pairing == 'members' and self.source == self.target
if self.schedule == 'member' and self.pairing in {'members', 'pairs', 'neighbor-pairs'} and not own:
error_msg = (
f'The actors of {self.label!r} draw partners by `random`, `neighbors`, `cell`, or `site`, or pair '
'with themselves by `members` within one group.'
)
raise ValidationError(error_msg)
if self.pairing == 'site' and self.actor != 'source':
error_msg = f'The sources of {self.label!r} move amounts into a field. Give `actor="source"`.'
raise ValidationError(error_msg)
if self.schedule == 'member' and self.form != 'agent':
error_msg = f'The schedule `member` of {self.label!r} runs agent rules. Use `form="agent"`.'
raise ValidationError(error_msg)
if isinstance(self.pairing, PartnerSpec) and self.schedule == 'member' and self.pairing.capacity is not None:
error_msg = (
f'The partners of {self.label!r} have a capacity of one step, but under the schedule `member` every '
'actor draws alone. Use the schedule `random-sequential`, or no capacity.'
)
raise ValidationError(error_msg)
if isinstance(self.pairing, PartnerSpec) and self.partners is not None:
error_msg = f'{self.label!r} draws partners by a `PartnerSpec`. Give its capacity, not `partners`.'
raise ValidationError(error_msg)
if self.partners is not None and (isinstance(self.partners, bool) or not self.partners >= 1):
error_msg = f'The partners of {self.label!r} are {self.partners!r}. Give a positive whole number, or None.'
raise ValidationError(error_msg)
def _get_owner(self) -> str:
"""Return the path of the element or group the source and target roles both bind, or the empty string."""
roles = self._get_specs()
bound = all(isinstance(roles.get(self._get_role(role)), RoleSpec) for role in ('source', 'target'))
return self.source if bound and self.source == self.target else ''
def _get_holders(self, namespace: SimpleNamespace) -> list[tuple[SimpleNamespace, SimpleNamespace]]:
"""Return the namespaces of the source and the target of every pair the symbolic form transfers between."""
return [(getattr(namespace, self._get_role('source')), getattr(namespace, self._get_role('target')))]
def _shift(self, holder: SimpleNamespace) -> SimpleNamespace:
"""Return the namespace of a group with every symbol read at the member the offset names."""
if not self.offset:
return holder
return SimpleNamespace(
**{
name: (
shift_member(value, self.offset, boundary=self.boundary) if isinstance(value, Expression) else value
)
for name, value in vars(holder).items()
},
)
@staticmethod
def _get_holder(namespace: SimpleNamespace, end: str, own: SimpleNamespace) -> object:
"""Return the symbol of a leg end, the state of the role it names or the state of its own side."""
role, _, state = end.rpartition('.')
return getattr(getattr(namespace, role) if role else own, state)
def _get_law(
self,
namespace: SimpleNamespace,
leg: NamedLeg,
rate: Amount | None,
reverse: bool,
) -> TransferLaw:
"""Return the transfer law of one leg of a pair, forward or back, reported under its flow name."""
state, into, amount, flow = leg
owner = self._get_owner() if flow else ''
held = (
state
if isinstance(state, Source)
else self._get_holder(namespace, state, self._get_member(namespace, 'source'))
)
sink = (
into
if isinstance(into, Sink)
else self._get_holder(namespace, str(into), self._get_member(namespace, 'target'))
)
value = amount(namespace) if callable(amount) else amount
if reverse:
held, sink = (
Source(sink.label) if isinstance(sink, Sink) else sink,
Sink(held.label) if isinstance(held, Source) else held,
)
flow = f'{flow}_reverse' if flow else ''
if rate is None:
if self.guard is not None:
value = where(self.guard(namespace), value, 0)
return TransferLaw(held, sink, value, label=flow, owner=owner)
jumps = rate(namespace) if callable(rate) else rate
if self.guard is not None:
jumps = where(self.guard(namespace), jumps, 0)
return TransferLaw(held, sink, value, jumps, label=flow, owner=owner)
def _get_transfer_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the transfer laws of every leg of every pair, forward and back, of the symbolic form."""
laws: list[BaseLaw] = []
for source, target in self._get_holders(namespace):
paired = SimpleNamespace(
**{**vars(namespace), self._get_role('source'): source, self._get_role('target'): self._shift(target)},
)
laws += [self._get_law(paired, leg, self.rate, False) for leg in self._get_legs()]
if self.reverse_rate is not None:
laws += [self._get_law(paired, leg, self.reverse_rate, True) for leg in self._get_legs()]
return laws
def _get_pairing(self) -> MatchLaw:
"""Return the law of the pairs of the pairing, whose capacity for a function is `partners`, or 1 without it."""
pairing = self.pairing
spec = NAMED_PARTNERS[pairing] if isinstance(pairing, str) else pairing
return MatchLaw(self.source, self.target, spec, capacity=self.partners or 1)
def _get_pairs(self, model: ModelViewProtocol) -> list[tuple[AgentViewProtocol, AgentViewProtocol]]:
"""Return the pairs of source and target agents of one step in random order."""
sources = model.agents(self._get_role('source'))
targets = sources if self.source == self.target else model.agents(self._get_role('target'))
if self.fallback is not None and not any(target.alive for target in targets):
targets = model.agents(self._get_role('fallback'))
pairs: list[tuple[AgentViewProtocol, AgentViewProtocol]] = []
if self.pairing == 'members':
pairs = list(zip(sources, targets, strict=False))
else:
pairs = cast(
list[tuple[AgentViewProtocol, AgentViewProtocol]],
self._get_pairing().get_pairs(sources, targets, model.rng),
)
return [pairs[int(index)] for index in model.rng.permutation(len(pairs))]
def _evaluate(
self,
amount: Amount,
source: AgentViewProtocol,
target: object,
model: ModelViewProtocol,
) -> float:
"""Return the amount of a leg between two agents as a number."""
value = amount(source, target, model) if callable(amount) else amount
if not isinstance(value, int | float | np.floating | np.integer):
error_msg = f'An amount of {self.label!r} between agents is {value!r}. Return a number.'
raise ValidationError(error_msg)
return float(value)
def _transfer_pair(self, source: AgentViewProtocol, target: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Transfer between two living agents, and remove the member of the pair that a removal names."""
if source.alive and target.alive and self.transfer(source, target, model) and self.remove is not None:
(source if self.remove == 'source' else target).remove()
def _deposit(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Move the amounts of every leg from an agent into the site of the target field at its cell or node."""
cell = agent.location_key
if not agent.alive or (self.guard is not None and not self.guard(agent, cell, model)):
return
legs = self._get_legs()
amounts = [self._evaluate(amount, agent, cell, model) for _, _, amount, _ in legs]
field = self._get_role('target')
for (state, into, _, _), amount in zip(legs, amounts, strict=True):
if isinstance(state, str):
setattr(agent, state, getattr(agent, state) - amount)
if isinstance(into, str):
model.deposit(f'{field}.{into}', cell, amount)
def _get_partners(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> list[AgentViewProtocol]:
"""Return the partners of an actor by the pairing, or the number of them drawn at random."""
if self.pairing == 'members':
return [agent]
role = self._get_role('source' if self.actor == 'target' else 'target')
if callable(self.pairing) or isinstance(self.pairing, PartnerSpec):
others = [other for other in model.agents(role) if other is not agent]
pairs = self._get_pairing().get_pairs([agent], others, model.rng)
return [cast(AgentViewProtocol, partner) for _, partner in pairs]
if self.pairing == 'cell':
groups = set(model.groups(role))
others = [other for other in agent.agents_here() if other.group in groups]
else:
others = [other for other in model.agents(role) if other is not agent]
if self.pairing == 'neighbors':
cells = set(agent.neighbor_cells())
others = [other for other in others if other.location_key in cells]
elif self.pairing == 'random' and others:
others = [others[int(model.rng.integers(len(others)))]]
if self.partners is None or len(others) <= self.partners:
return others
return [others[int(index)] for index in model.rng.permutation(len(others))[: self.partners]]
def _transfer_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Transfer between one actor and its partners, in the activation of the actor."""
if self.pairing == 'site':
self._deposit(agent, model)
return
for partner in self._get_partners(agent, model):
if agent.alive:
pair = (partner, agent) if self.actor == 'target' else (agent, partner)
self._transfer_pair(*pair, model)
def _transfer_pairs(self, model: ModelViewProtocol) -> None:
"""Transfer between every pair of agents of one step, or into the sites of the source agents."""
if self.pairing == 'site':
sources = model.agents(self._get_role('source'))
for index in model.rng.permutation(len(sources)):
self._deposit(sources[int(index)], model)
return
for source, target in self._get_pairs(model):
self._transfer_pair(source, target, model)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the transfer laws of the symbolic form, or the rule of the agent form.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
A transfer law per leg, or the rule that transfers between the pairs of agents.
"""
if self.form == 'agent' and self.schedule == 'member':
return [RuleLaw(self._transfer_member, phase=self.phase, members=self._get_role(self.actor))]
if self.form == 'agent':
return [RuleLaw(self._transfer_pairs, phase=self.phase)]
return self._get_transfer_laws(namespace)
def transfer(self, source: AgentViewProtocol, target: AgentViewProtocol, model: ModelViewProtocol) -> bool:
"""Move the amounts of every leg from one source agent to one target agent, when the guard holds.
Args:
source:
The agent the amounts leave.
target:
The agent the amounts enter.
model:
The model view of the rule.
Returns:
transferred:
Whether the guard held, so the amounts moved.
Raises:
ValidationError:
An amount function returns something other than a number.
SimulationError:
A leg end names the state of a referenced member, and the reference state names no member.
"""
if self.guard is not None and not self.guard(source, target, model):
return False
legs = self._get_legs()
amounts = [self._evaluate(amount, source, target, model) for _, _, amount, _ in legs]
if self.settlement is not None:
self.settlement(source, target, amounts, model)
return True
ends = [
(
self._get_end(source, state) if isinstance(state, str) else None,
self._get_end(target, into) if isinstance(into, str) else None,
)
for state, into, _, _ in legs
]
for (taken, given), amount in zip(ends, amounts, strict=True):
if taken is not None:
holder, name = taken
setattr(holder, name, getattr(holder, name) - amount)
if given is not None:
holder, name = given
setattr(holder, name, getattr(holder, name) + amount)
return True
def _get_end(self, agent: AgentViewProtocol, name: str) -> tuple[AgentViewProtocol, str]:
"""Return the agent that holds the state a leg end names and the name of the state.
A name `'<reference>.<state>'` names the state of the member that the reference state of the agent names.
"""
reference, _, state = name.rpartition('.')
if not reference:
return agent, name
holder = agent.referenced(reference)
if holder is None:
error_msg = (
f'A leg of {self.label!r} moves {name!r}, but the state {reference!r} of agent {agent.id} names no '
'member.'
)
raise SimulationError(error_msg)
return holder, state
Equations
network.s[i]: diff(s_network(t), t) == -(_reduced_1*beta_network*k_network*s_network(t) - _reduced_2*beta_network*k_network*s_network(t))/_reduced_0
network.i[i]: diff(i_network(t), t) == -mu_network*i_network(t) + (_reduced_1*beta_network*k_network*s_network(t) - _reduced_2*beta_network*k_network*s_network(t))/_reduced_0
network.r[i]: diff(r_network(t), t) == mu_network*i_network(t)