Heterogeneous mean-field SIS
The SIS epidemic on a network with a power-law degree distribution, in heterogeneous mean field.
- Name:
heterogeneous_mean_field_sis - Field:
biology/epidemiology - Formalism: ordinary differential equations
- References: Pastor-Satorras and Vespignani 2001; Pastor-Satorras et al. 2015
Factory
@register_model(
field='biology/epidemiology',
tags=('network', 'mean field', 'ode', 'epidemic threshold', 'scale-free'),
references=('Pastor-Satorras and Vespignani 2001', 'Pastor-Satorras et al. 2015'),
defaults={'until': 50.0, 'dt': 0.5},
)
def heterogeneous_mean_field_sis(
exponent: float = 2.5,
k_min: int = 3,
k_max: int = 30,
*,
lam: float = 0.2,
rho: float = 0.01,
) -> System:
"""Return the SIS epidemic on a network with a power-law degree distribution, in heterogeneous mean field.
The degrees run from `k_min` to `k_max` with `P(k) ~ k^-exponent`. The infection persists above the epidemic
threshold `lam_c = <k> / <k^2>`, which vanishes on scale-free networks as `k_max` grows, and dies out below it.
Args:
exponent:
The exponent of the degree distribution.
k_min:
The smallest degree.
k_max:
The largest degree.
lam:
The spreading rate, the ratio of the infection and recovery rates.
rho:
The initial infected density of every class.
Returns:
system:
The model, with the densities `network.rho` of the degree classes from `k_min` to `k_max`, members of the
group `network`.
Examples:
>>> from skcomplex.models.biology.epidemiology import heterogeneous_mean_field_sis
>>> heterogeneous_mean_field_sis(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=DegreeClassElement,
size=len(degrees),
values={'rho': rho, 'lam': lam, 'k': degrees, 'P': weights},
)
interactions = [
RateInteraction('recovery', target='network', state='rho', rate=_recovery),
CouplingInteraction(
'infection',
target='network',
read='rho',
write='rho',
function=_infection,
reduce='normalized-mean',
mode='group',
weights=_edge_weight,
include_self=True,
),
]
return System(elements=[classes], interactions=interactions)
Components
class DegreeClassElement(Element):
"""Infected density of the nodes of one degree of a network, in heterogeneous mean field.
A node of degree `k` recovers at rate 1 and is infected at rate `lam k Theta`, where
`Theta = sum_k k P(k) rho_k / <k>` is the probability that an edge points to an infected node, a total over the
degree classes, so `d rho_k / dt = -rho_k + lam k (1 - rho_k) Theta`. The classes are the members of a group, and
the recovery and the infection are interactions of the group, which `heterogeneous_mean_field_sis` adds.
Args:
label:
The label of the class or the group of classes.
visual:
How the element is drawn, or `None` for `default_visual_`.
rho:
The initial infected density.
lam:
The spreading rate, the ratio of the infection and recovery rates.
k:
The degree of the class.
P:
The probability `P(k)` of the degree.
Examples:
>>> from skcomplex.core.elements import GroupElement
>>> from skcomplex.models.biology.epidemiology import DegreeClassElement
>>> GroupElement(
... 'network', element=DegreeClassElement, size=2, values={'k': [1.0, 2.0], 'P': [0.5, 0.5]}
... ).check().prepare().size_
2
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
rho: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
lam: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
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.rho = rho
self.lam = lam
self.k = k
self.P = P
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
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)]
Equations
network.rho[i]: diff(rho_network(t), t) == -(k_network*lam_network*rho_network(t) - k_network*lam_network)*_reduced_1/_reduced_0 - rho_network(t)