Adaptive SIS with rewiring
The SIS epidemic on a network whose susceptible nodes cut their edges to infected neighbours.
- Name:
adaptive_sis - Field:
biology/epidemiology - Formalism: jump process
- References: Gross, Dommar D'Lima, and Blasius 2006; Thurner et al. 2018; Pastor-Satorras et al. 2015
Factory
@register_model(
field='biology/epidemiology',
tags=('network', 'adaptive network', 'rewiring', 'jumps', 'epidemic threshold'),
references=("Gross, Dommar D'Lima, and Blasius 2006", 'Thurner et al. 2018', 'Pastor-Satorras et al. 2015'),
defaults={'until': 30.0, 'dt': 0.25},
)
def adaptive_sis(
n: int = 200,
k: float = 6.0,
*,
beta: float = 0.4,
gamma: float = 1.0,
w: float = 1.0,
infected: float = 0.2,
seed: int | None = 1,
) -> System:
"""Return the SIS epidemic on a network whose susceptible nodes cut their edges to infected neighbours.
Every edge between a susceptible and an infected node is rewired at rate `w`: the susceptible end keeps it and
joins it to another susceptible node. Rewiring isolates the infected nodes and raises the epidemic threshold above
its value on the fixed network, while the susceptible nodes cluster.
Args:
n:
The number of nodes.
k:
The mean degree of the initial Erdős-Rényi network.
beta:
The infection rate per edge to an infected neighbour.
gamma:
The recovery rate.
w:
The rewiring rate of an edge between a susceptible and an infected node.
infected:
The initial share of infected nodes.
seed:
The seed of the network and of the initial infections.
Returns:
system:
The model, with the states `hosts.i` of the variable group `hosts` on the nodes and the variable group
`contacts` on the edges.
Examples:
>>> from skcomplex.models.biology.epidemiology import adaptive_sis
>>> [component.label for component in adaptive_sis(n=20).elements]
['hosts', 'contacts']
"""
space = _random_network(n, k, seed)
initial = (np.random.default_rng(seed).random(n) < infected).astype(float)
hosts = GroupElement(
'hosts',
element=HostElement,
space=space,
on='sites',
variable=True,
capacity=n,
values={'i': initial, 'w': w},
visual=VisualSpec('sites', {'states': ['i'], 'low': 0, 'high': 1}),
)
edges = space.graph.number_of_edges()
contacts = GroupElement(
'contacts',
element=Element,
space=space,
on='edges',
variable=True,
capacity=edges,
visual=HIDDEN,
)
spreading = SpreadingInteraction('spreading', target='hosts', state='i', beta=beta, gamma=gamma)
return System(
elements=[hosts, contacts],
interactions=[spreading, _get_rewiring('source'), _get_rewiring('target')],
)
Components
class HostElement(Element):
"""Node of a network that is susceptible at `i = 0` or infected at `i = 1`, in continuous time.
A `SpreadingInteraction` infects and recovers it, and a `TopologyInteraction` rewires its edges at its rate `w`.
Args:
label:
The label of the node or of the group of nodes.
visual:
How the element is drawn, or `None` for `default_visual_`.
i:
Whether the node is infected.
w:
The rate at which the node, while susceptible, rewires an edge to an infected neighbour.
Examples:
>>> from skcomplex.models.biology.epidemiology import HostElement
>>> HostElement('node').w
1.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
i: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
w: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.i = i
self.w = w
class Element(BaseElement):
"""Base class of entities that carry state, such as particles, populations, agents, or firms.
A subclass takes its states, its parameters, and `visual` as arguments of its own `__init__`, annotating each state
with a `StateSpec` and each parameter with a `ParameterSpec` inside `Annotated`, and returns its intrinsic dynamics
from `_get_laws`, which returns no laws by default. A subclass drawn in its own way gives `visual` a default in its
signature.
An element whose laws have the agent trait, such as rules, alone or in a group, is an agent: under discrete time it
runs on the agent stepper, whose rules receive every member as an agent with its states as attributes. Its states
then hold any value, not only numbers. A member is removed only when its states tagged with a quantity are settled
to zero, and its states and parameters cannot take the names of the agent view, such as `location` or `alive`, or
of the attributes Mesa reads, such as `position` or `model`.
Args:
label:
The label of the element, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
Attributes:
default_visual_:
The visual of the element when it sets none.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.specs import StateSpec
>>> from skcomplex.specs import VisualSpec
>>> from sage.all__sagemath_symbolics import oo
>>> from sage.sets.real_set import RealSet
>>> class PopulationElement(Element):
... def __init__(
... self,
... label,
... n: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
... visual=VisualSpec('line'),
... ):
... super().__init__(label, visual=visual)
... self.n = n
>>> PopulationElement('rabbits', n=40)
PopulationElement(label='rabbits', n=40)
"""
@classmethod
def _get_view(cls) -> type | None:
"""Return the class of the views of the element's agents, or `None` for the engine's plain view.
An agent element names a subclass of the engine's `AgentView` that annotates its states and parameters, such as
`class WordView(AgentView): count: float`, so its rules read them with checked types.
"""
return None
def _get_laws(self, namespace: SimpleNamespace) -> Sequence[BasePart]:
"""Return the laws of the element for a namespace of symbols, none by default.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the element.
"""
return []
class SpreadingInteraction(TransitionInteraction):
"""Spreading of the active value of a discrete state along the edges of a network or the bonds of a lattice.
An inactive member becomes active at the rate `beta` times the number of its active neighbours, or times their
fraction of its neighbours. An active member recovers at the rate `gamma`, to the inactive value, or to the removed
value when there is one, and a removed member returns to the inactive value at the rate `omega`. Every move is a
transition at a rate, so the state changes by random events in continuous time.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group on the nodes or sites.
state:
The name of the discrete state of the members.
beta:
The rate of spreading from one active neighbour, or from all neighbours when they are all active.
gamma:
The rate of recovery of an active member.
omega:
The rate at which a removed member returns to the inactive value.
inactive:
The value of the state of an inactive member.
active:
The value of the state of an active member.
removed:
The value of the state of a recovered member, or `None` for members that recover to the inactive value.
form:
`'number'` for a rate proportional to the number of active neighbours, or `'fraction'` for one
proportional to their fraction of the neighbours.
layer:
The layer of the network whose edges the spreading runs on, or `None` for the graph of the network or the
neighbourhood of the lattice.
Raises:
ValidationError:
The form is unknown, or the values of the state are not distinct.
Examples:
>>> from skcomplex.models.biology.epidemiology import SpreadingInteraction
>>> spreading = SpreadingInteraction('spreading', target='sites', beta=2.0, gamma=1.0, form='fraction')
>>> spreading.check().form
'fraction'
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'state',
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
omega: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
inactive: float = 0.0,
active: float = 1.0,
removed: float | None = None,
form: str = 'number',
layer: str | None = None,
) -> None:
super().__init__(label, target=target, state=state)
self.beta = beta
self.gamma = gamma
self.omega = omega
self.inactive = inactive
self.active = active
self.removed = removed
self.form = form
self.layer = layer
def _check_transition(self) -> None:
"""Raise an error for an unknown form or values of the state that are not distinct.
Raises:
ValidationError:
The form is unknown, or the values of the state are not distinct.
"""
if self.form not in FORMS:
error_msg = f'The form of spreading is {self.form!r}. Use one of {", ".join(FORMS)}.'
raise ValidationError(error_msg)
values = [self.inactive, self.active, *([] if self.removed is None else [self.removed])]
if len(set(values)) < len(values):
error_msg = (
f'The inactive, active, and removed values of spreading are {values}. Give each a distinct value.'
)
raise ValidationError(error_msg)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the spreading from the active neighbours, the recovery, and the return of every member.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The jumps of the state.
"""
member = self._get_member(namespace, 'target')
value = getattr(member, self.state)
over = None if self.layer is None else AdjacencyNeighborhood(layer=self.layer)
reduce = 'mean' if self.form == 'fraction' else 'sum'
held = neighbors(where(getattr(member.other, self.state) == self.active, 1, 0), over=over, reduce=reduce)
recovered = self.inactive if self.removed is None else self.removed
laws: list[BaseLaw] = [
self._get_jump(value, self.inactive, self.active, self._get_parameter_symbol('beta', namespace) * held),
self._get_jump(value, self.active, recovered, self._get_parameter_symbol('gamma', namespace)),
]
if self.removed is not None:
laws.append(
self._get_jump(value, self.removed, self.inactive, self._get_parameter_symbol('omega', namespace)),
)
return laws
class TopologyInteraction(NeighborInteraction):
"""Change of the topology of a network: nodes and edges that are added or removed, and edges that are rewired.
The action runs for every node of the group a role binds, `'add-node'`, `'remove-node'`, `'replace-node'`, which
removes the node and adds a new one, and `'add-edge'`, or for every edge of its network, `'remove-edge'` and
`'rewire'`. Without a target, `'add-node'` runs once per step, adding one node from outside the network to the one
network whose groups are variable, as in growing networks. It fires at a `rate`, with a `probability` per step, or
with probability one, and only where the `condition` holds. The rate, the probability, and the condition are numbers
or functions of the namespace of an edge, whose ends are `namespace.source` and `namespace.target`, and an action of
the nodes reads the node as the source. A new node joins `edges` other nodes, a node adding edges joins `edges`
other nodes, and a rewiring keeps the end `keep` of its edge. New ends are drawn with probability proportional to
`toward(node)`, a function of the namespace of a candidate node, or uniformly. A new node or edge takes the
`values`, each a value or a function of the random number generator that draws it. The groups on the nodes and the
edges of the network are variable. With `hyperedge`, a new node instead forms a hyperedge of that size with
`hyperedge - 1` existing nodes and joins each of them by an edge: the first two are the ends of an edge drawn
uniformly, so a size of three attaches a triangle along an edge, and further members are drawn by `toward`. The
new hyperedge is a new member of the variable group on the hyperedges of the network.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group on the nodes of the network, or `None` for an addition of one node per step.
action:
The change, `'add-node'`, `'remove-node'`, `'replace-node'`, `'add-edge'`, `'remove-edge'`, or
`'rewire'`.
rate:
The rate of the change, a number or a function, or `None` for a probability per step.
probability:
The probability of the change per step, a number or a function, or `None` for probability one.
condition:
A function of the namespace that returns where the change may fire, or `None` everywhere.
toward:
A function of the namespace of a candidate node that returns its weight as a new end, or `None` for equal
weights.
edges:
The number of edges a new node or a node adding edges joins.
keep:
The end of an edge a rewiring keeps, `'source'` or `'target'`.
values:
The initial states of a new node or edge by name, each a value or a function of the random number
generator that draws it, which default to the defaults of the element type.
layer:
The layer of the network whose edges change, or `None` for the graph of the network.
hyperedge:
The size of the hyperedge every new node forms with existing nodes, or `None` for none.
Raises:
ValidationError:
The action or the kept end is unknown, both a rate and a probability are given, the number of edges is
not a whole number, an action other than `'add-node'` has no target, or a hyperedge is given to an action
other than `'add-node'` or has fewer than two members, raised by `check`.
Examples:
>>> import networkx as nx
>>> from skcomplex.core.elements import Element, GroupElement
>>> from skcomplex.core.interactions import TopologyInteraction
>>> from skcomplex.core.spaces import NetworkSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> nodes = GroupElement('nodes', element=Element, on='sites', variable=True, capacity=20)
>>> links = GroupElement('links', element=Element, on='edges', variable=True, capacity=20)
>>> expiry = TopologyInteraction('expiry', target='nodes', action='remove-edge')
>>> system = System(elements=[nodes, links], interactions=[expiry], space=NetworkSpace(nx.cycle_graph(5)))
>>> trajectory = Simulation(system, random_state=0).check().prepare().run(until=1).trajectory_
>>> len(trajectory.network(0).sources), len(trajectory.network(1).sources)
(5, 0)
A network that grows by triangles from one triangle:
>>> import xgi
>>> from skcomplex.core.spaces import DiscreteTime
>>> triangles = GroupElement('triangles', element=Element, on='hyperedges', variable=True, capacity=8)
>>> growth = TopologyInteraction('growth', action='add-node', hyperedge=3)
>>> space = NetworkSpace(nx.complete_graph(3), hyperedges=[[0, 1, 2]])
>>> system = System(
... elements=[nodes, links, triangles], interactions=[growth], space=space, time=DiscreteTime()
... )
>>> trajectory = Simulation(system, random_state=0).check().prepare().run(until=5).trajectory_
>>> trajectory.network(-1).hyperedges.shape
(6, 3)
"""
def __init__(
self,
label: str,
*,
target: Annotated[str | None, RoleSpec(optional=True)] = None,
action: str = 'add-edge',
rate: Chance | None = None,
probability: Chance | None = None,
condition: Callable[[SimpleNamespace], object] | None = None,
toward: Callable[[SimpleNamespace], object] | None = None,
edges: int = 1,
keep: str = 'source',
values: Mapping[str, object] | None = None,
layer: str | None = None,
hyperedge: int | None = None,
) -> None:
super().__init__(label, layer=layer)
self.target = target
self.action = action
self.rate = rate
self.probability = probability
self.condition = condition
self.toward = toward
self.edges = edges
self.keep = keep
self.values = values
self.hyperedge = hyperedge
def _check(self) -> None:
"""Raise an error for an unknown action or kept end, a rate next to a probability, or a bad number of edges.
Raises:
ValidationError:
The action or the kept end is unknown, both a rate and a probability are given, a condition or a
weight cannot be called, the number of edges is not a whole number of at least 0, an action other
than `'add-node'` has no target, or a hyperedge is given to an action other than `'add-node'`.
"""
super()._check()
self._check_variant('action', self.action, ACTIONS, (BaseSpec,))
if self.rate is not None and self.probability is not None:
error_msg = f'{self.label!r} fires at a rate and with a probability. Give one of them.'
raise ValidationError(error_msg)
if self.target is None and self.action != 'add-node':
error_msg = f'{self.label!r} runs {self.action!r} without a target. Give the group on the nodes.'
raise ValidationError(error_msg)
if self.hyperedge is not None and self.action != 'add-node':
error_msg = f'{self.label!r} runs {self.action!r} with a hyperedge. Only new nodes form hyperedges.'
raise ValidationError(error_msg)
for name in ('condition', 'toward'):
value = getattr(self, name)
if value is not None and not callable(value):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace, or None.'
raise ValidationError(error_msg)
for request in self._get_requests():
request.check()
def _get_requests(self) -> list[TopologyRequestSpec]:
"""Return the requests of the action, two for a replacement of a node."""
if self.action not in ACTIONS:
raise self._get_variant_error(self.action, 'action', '_get_requests')
values = tuple((self.values or {}).items())
kinds = ['add-node', 'remove-node'] if self.action == 'replace-node' else [self.action]
return [
TopologyRequestSpec(
kind,
self.edges,
self.toward,
self.keep,
self.layer,
values,
self.hyperedge if kind == 'add-node' else None,
)
for kind in kinds
]
@staticmethod
def _evaluate(value: Chance, namespace: SimpleNamespace) -> object:
"""Return a number, or the value of a function of the namespace."""
return value(namespace) if callable(value) else value
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the requests of the action, fired at the rate or with the probability where the condition holds.
Args:
namespace:
The namespace of an edge, with `namespace.source` and `namespace.target`.
Returns:
laws:
A jump law per request at a rate, or else an update law per request with its probability.
"""
chance: object = 1
if self.rate is not None:
chance = self._evaluate(self.rate, namespace)
elif self.probability is not None:
chance = self._evaluate(self.probability, namespace)
if self.condition is not None:
chance = where(self.condition(namespace), chance, 0)
if self.rate is not None:
return [JumpLaw(chance, {request: 1}) for request in self._get_requests()]
return [UpdateLaw(request, chance) for request in self._get_requests()]