Neutral birth-death with hopping between boxes
Organisms that divide, die, and hop between patches on a periodic lattice, and gather in clusters.
- Name:
neutral_clustering - Field:
biology/ecology - Formalism: jump process
- References: Bertin 2021; Houchmandzadeh 2002; Houchmandzadeh 2009
Factory
@register_model(
field='biology/ecology',
tags=('neutral theory', 'clustering', 'jump process', 'lattice', 'birth-death', 'diffusion'),
references=('Bertin 2021', 'Houchmandzadeh 2002', 'Houchmandzadeh 2009'),
defaults={'until': 50.0, 'dt': 0.5, 'backend': NumbaBackend()},
observables=lambda system: [
ObservableSpec('mean', group='patches', variable='n', reduce='mean'),
ObservableSpec('variance', group='patches', variable='n', reduce='var'),
],
)
def neutral_clustering(
shape: int | Sequence[int] = 64,
*,
n: float = 10.0,
alpha: float = 1.0,
beta: float = 0.5,
) -> System:
"""Return organisms that divide, die, and hop between patches on a periodic lattice, and gather in clusters.
Births balance deaths, so the mean number of organisms per patch stays `n` while its variance grows. Without hops it
grows as `2 alpha n t`, so most patches empty while a few fill. Hops at the rate `beta` slow the growth to
`t^(1/2)` in one dimension and `ln t` in two, so the clusters persist in low dimension although the mean density
obeys a diffusion equation.
Args:
shape:
The number of patches of a ring, or the shape of a periodic lattice of patches.
n:
The initial number of organisms in every patch.
alpha:
The rate at which an organism divides, equal to the rate at which it dies.
beta:
The rate at which an organism hops to one neighbouring patch.
Returns:
system:
The model, with the numbers `patches.n` of organisms on the lattice `patches`, and their mean `mean` and
variance `variance` over the patches.
Examples:
>>> from skcomplex.models.biology.ecology import neutral_clustering
>>> neutral_clustering(8).elements[0].check().prepare().size_
8
"""
lattice = (shape,) if isinstance(shape, int) else tuple(shape)
patches = LatticeElement(
'patches',
element=HabitatPatchElement,
space=LatticeSpace(lattice),
values={'n': np.full(lattice, float(n)), 'alpha': alpha},
)
members = [f'patches[{index}]' for index in range(int(np.prod(lattice)))]
interactions: list[BaseInteraction] = [
JumpInteraction('division', target='patches', changes={'n': 1}, rate=lambda s: s.target.alpha * s.target.n),
JumpInteraction('death', target='patches', changes={'n': -1}, rate=lambda s: s.target.alpha * s.target.n),
]
if beta:
interactions.append(HoppingInteraction('hopping', patches=members, links=_neighbor_pairs(lattice), beta=beta))
return System(
elements=[patches],
interactions=interactions,
)
Components
class HabitatPatchElement(CompartmentElement):
"""Patch of space, a compartment of `n` organisms that divide and die at the same rate `alpha`.
The divisions and deaths are jumps of `neutral_clustering`, so the number of organisms is a symmetric random walk,
and an empty patch stays empty without arrivals from its neighbours.
Args:
label:
The label of the patch or the lattice of patches.
visual:
How the element is drawn, or `None` for `default_visual_`.
n:
The initial number of organisms.
alpha:
The rate at which an organism divides, equal to the rate at which it dies.
Examples:
>>> from skcomplex.models.biology.ecology import HabitatPatchElement
>>> HabitatPatchElement('patch', n=5.0).alpha
1.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = VisualSpec('image', {'low': 0.0}),
*,
n: Annotated[float, StateSpec(domain=NN)] = 10.0,
alpha: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual, n=n)
self.alpha = alpha
class JumpInteraction(Interaction):
"""Random event that changes states of the target by signed amounts, at a rate, by a probability, or after a wait.
With a `rate`, the event fires at the times of a Poisson process of that propensity in continuous time. With a
`probability`, it fires in a step of discrete time with that probability. With a `waiting` time, usually a draw from
a distribution with nonnegative support, it fires after a time drawn anew at every firing, as in renewal processes,
and every member of a bound group fires on its own clock. When the event fires it adds the `changes` to the states
and, with a waiting time, sets the states of `reset` to their new values. Every function reads `namespace.target`,
the states and parameters of the target, and the time `namespace.t`.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states the event changes.
changes:
A mapping from the name of a state of the target to the signed amount the event adds to it.
rate:
The function `rate(namespace)` of the propensity of the event, or `None`.
probability:
The function `probability(namespace)` of the probability that the event fires in a step, or `None`.
waiting:
The function `waiting(namespace)` of the time until the next firing, or `None`.
reset:
With a waiting time, a mapping from the name of a state of the target to the function `f(namespace)` of its
new value, one value per component of a vector state, or `None`.
Raises:
ValidationError:
The event has not exactly one clock, changes nothing, or resets states without a waiting time, or a
function is not callable.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import JumpInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> arrivals = JumpInteraction('arrivals', target='queue', changes={'n': 1}, rate=lambda s: 2.0)
>>> system = System(elements=[CompartmentElement('queue')], interactions=[arrivals])
>>> simulation = Simulation(system, dt=1.0, random_state=0).check().prepare()
>>> counts = simulation.run(until=5).trajectory_.values('queue.n')
>>> bool(all(count == round(count) for count in counts))
True
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
changes: Mapping[str, float] | None = None,
rate: Callable[[SimpleNamespace], object] | None = None,
probability: Callable[[SimpleNamespace], object] | None = None,
waiting: Callable[[SimpleNamespace], object] | None = None,
reset: Mapping[str, Callable[[SimpleNamespace], object]] | None = None,
) -> None:
super().__init__(label)
self.target: str | None = target
self.changes = changes
self.rate = rate
self.probability = probability
self.waiting = waiting
self.reset = reset
def _get_clocks(self) -> list[str]:
"""Return the names of the clocks the event is given."""
return [clock for clock in CLOCKS if getattr(self, clock) is not None]
def _check_jump(self) -> None:
"""Raise an error for an event without exactly one clock, without changes, or with a reset of no waiting time.
Raises:
ValidationError:
The event has not exactly one clock, changes nothing, or resets states without a waiting time, or a
function is not callable.
"""
clocks = self._get_clocks()
if len(clocks) != 1:
error_msg = f'The jump {self.label!r} has the clocks {clocks}. Give exactly one of {", ".join(CLOCKS)}.'
raise ValidationError(error_msg)
if not callable(getattr(self, clocks[0])):
error_msg = f'The {clocks[0]} of the jump {self.label!r} must be a function of the namespace.'
raise ValidationError(error_msg)
if self.reset is not None and self.waiting is None:
error_msg = f'The jump {self.label!r} resets states, which only a waiting time does. Give `waiting`.'
raise ValidationError(error_msg)
if not self.changes and not self.reset:
error_msg = f'The jump {self.label!r} changes no state. Give `changes`, or `reset` with a waiting time.'
raise ValidationError(error_msg)
if not all(callable(function) for function in (self.reset or {}).values()):
error_msg = f'The reset of the jump {self.label!r} maps states to {self.reset!r}. Give functions.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for an event without exactly one clock, without changes, or with a reset of no waiting time.
Raises:
ValidationError:
The event has not exactly one clock, changes nothing, or resets states without a waiting time, or a
function is not callable.
"""
super()._check()
self._check_jump()
def _get_changes(self, namespace: SimpleNamespace) -> dict[Expression, float]:
"""Return the amount the event adds to every state it changes."""
return {
getattr(self._get_member(namespace, 'target'), name): amount
for name, amount in (self.changes or {}).items()
}
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the propensity of the event."""
return SR(self.rate(namespace)) if self.rate is not None else SR(0)
def _get_steps(self, namespace: SimpleNamespace, changes: dict[Expression, float]) -> list[BaseLaw]:
"""Return the updates of the states by their amounts in a step, when one draw falls below the probability."""
probability = SR(self.probability(namespace)) if self.probability is not None else SR(0)
name = re.sub(r'\W', '_', f'jump_{self.label}')
fired = where(uniform(name) < probability, 1, 0)
return [UpdateLaw(state, state + amount * fired) for state, amount in changes.items()]
def _get_reset(self, state: object, value: object) -> dict[Expression, Expression]:
"""Return the new value of every component of a state, one per component of a vector state.
Args:
state:
The symbol of a scalar state or the vector of a vector state.
value:
The new value, one expression for every component or a sequence of one per component.
Returns:
reset:
The new value of every component of the state.
Raises:
ValidationError:
The new value of a vector state has not one component per component of the state.
"""
if not isinstance(state, FreeModuleElement):
return {SR(state): SR(value)}
components = list(value) if isinstance(value, FreeModuleElement | list | tuple) else [value] * len(state)
if len(components) != len(state):
error_msg = f'The jump {self.label!r} resets {len(state)} components with {len(components)} values.'
raise ValidationError(error_msg)
return {SR(item): SR(component) for item, component in zip(state, components, strict=True)}
def _get_renewal(self, namespace: SimpleNamespace, changes: dict[Expression, float]) -> JumpLaw:
"""Return the renewal of the event, which adds the changes and sets the reset states after every wait."""
reset: dict[Expression, object] = {state: state + amount for state, amount in changes.items()}
for name, function in (self.reset or {}).items():
reset.update(self._get_reset(getattr(self._get_member(namespace, 'target'), name), function(namespace)))
waiting = SR(self.waiting(namespace)) if self.waiting is not None else SR(0)
return JumpLaw(waiting=waiting, reset=reset, label=self.label)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the jump law of the event at a rate or after a wait, or the updates of a step with a probability.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The laws of the event.
"""
changes = self._get_changes(namespace)
if self.waiting is not None:
return [self._get_renewal(namespace, changes)]
if self.probability is not None:
return self._get_steps(namespace, changes)
return [JumpLaw(self._get_rate(namespace), changes)]
class HoppingInteraction(Interaction):
"""Hops of single organisms between neighbouring patches, each organism moving to each neighbour at the rate `beta`.
The patches are the entries of a variadic role, and every pair `(i, j)` of `links` gives a jump that moves one
organism from the patch of entry `i` to the patch of entry `j` at the rate `beta n_i`.
Args:
label:
The label of the interaction.
patches:
The paths of the patches, such as the members `patches[i]` of a lattice.
beta:
The rate at which an organism hops to one neighbouring patch.
links:
The ordered pairs of entries of `patches` between which organisms hop.
Examples:
>>> from skcomplex.models.biology.ecology import HoppingInteraction
>>> HoppingInteraction('hopping', patches=['a', 'b'], links=[(0, 1), (1, 0)]).beta
0.5
"""
def __init__(
self,
label: str,
*,
patches: Annotated[list[str], RoleSpec(variadic=True)],
beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
links: Sequence[tuple[int, int]] = (),
) -> None:
super().__init__(label)
self.patches = patches
self.beta = beta
self.links = links
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return one hop per ordered pair of neighbouring patches.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
patches = s.patches
return [JumpLaw(s.beta * patches[i].n, {patches[i].n: -1, patches[j].n: 1}) for i, j in self.links]