Multiscale tumour angiogenesis
A tumour that outgrows its oxygen supply and recruits new blood vessels from a parent vessel.
- Name:
tumour_angiogenesis - Field:
biology/oncology - Formalism: co-simulation
- References: Carmichael and Hadžikadić 2019; Anderson and Chaplain 1998
Factory
@register_model(
field='biology/oncology',
tags=('agents', 'field', 'multiscale', 'angiogenesis', 'co-simulation'),
references=('Carmichael and Hadžikadić 2019', 'Anderson and Chaplain 1998'),
defaults={'until': 120, 'dt': 1},
)
def tumour_angiogenesis(
width: int = 40,
height: int = 25,
*,
n_sprouts: int = 5,
radius: int = 2,
chemotaxis: float = 30.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return a tumour that outgrows its oxygen supply and recruits new blood vessels from a parent vessel.
The scales are coupled: tumour cells and vessel tips are agents on a grid, and oxygen and VEGF are fields on a
lattice of the same shape. The tumour starts far from the parent vessel along the left edge, so its cells become
hypoxic and secrete VEGF, which draws tips out of the vessel. The new vessels bring oxygen and the tumour grows.
Args:
width:
The number of columns.
height:
The number of rows.
n_sprouts:
The number of tips on the parent vessel.
radius:
The radius of the initial tumour, centred three quarters of the way to the right.
chemotaxis:
How strongly VEGF draws the tips.
random_state:
The seed or generator of the positions of the tips.
Returns:
system:
The model, with the agents `cells/tumour`, `cells/vessels`, and `cells/sprouts`, and the fields
`tissue/tissue.oxygen` and `tissue/tissue.vegf`.
Examples:
>>> from skcomplex.models.biology.oncology import tumour_angiogenesis
>>> [group.label for group in tumour_angiogenesis().systems[0].elements]
['tumour', 'vessels', 'sprouts']
"""
rng = np.random.default_rng(random_state)
center = (3 * width // 4, height // 2)
tumour = [
(x, y)
for x in range(center[0] - radius, center[0] + radius + 1)
for y in range(center[1] - radius, center[1] + radius + 1)
if (x - center[0]) ** 2 + (y - center[1]) ** 2 <= radius**2
]
tips = sorted(int(row) for row in rng.choice(height, size=n_sprouts, replace=False))
groups = [
GroupElement('tumour', element=TumourCellElement, size=len(tumour), variable=True, location=tumour),
GroupElement(
'vessels',
element=VesselElement,
size=height,
variable=True,
location=[(0, row) for row in range(height)],
),
GroupElement(
'sprouts',
element=SproutElement,
size=n_sprouts,
values={'chemotaxis': chemotaxis},
variable=True,
location=[(0, row) for row in tips],
),
]
tissue = FieldElement('tissue', element=TissueElement, space=LatticeSpace((width, height), boundary='reflecting'))
return System(
systems=[
System('cells', elements=groups, space=LatticeSpace((width, height), boundary='fixed')),
System('tissue', elements=[tissue], interactions=_get_tissue()),
],
interactions=_get_cells(),
)
Components
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
class MapInteraction(Interaction):
"""Map of a state of discrete time, whose next value is a function, or the state plus it.
In the mode `'replace'` the next value of the state is `f`, and in the mode `'increment'` it is the state plus `f`.
Under the `'synchronous'` schedule every member of a bound group updates at once, from the states at the start of
the step. Under the `'sequential'` schedule the members update one after the other, in their order, and under the
`'random-sequential'` schedule in a new random order every step, each member reading the values the members before
it set. Under the `'member'` schedule of the agent form, every member updates in its own activation, together with
the rules of its element and the other member interactions, in the order the system lists them.
In the `'symbolic'` form, `function(namespace)` returns a symbolic expression, or a vector of them for a vector
state, where `namespace.target` holds the states and parameters of the target, `namespace.sources` the namespaces
of the sources, and `namespace.t` the time, and may hold draws such as `normal('noise')`. In the `'agent'` form,
the target is an agent group and `function(agent, model)` returns the value as a number, or as an array of the
shape of a vector state, so it can read the agent, its neighbours, and the history and storage of the model view.
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 map sets.
state:
The name of the state of the target.
function:
The function `f(namespace)` in the symbolic form, or `f(agent, model)` in the agent form.
sources:
The paths of the elements, or members of groups such as `cells[0]`, that the function reads.
mode:
`'replace'` or `'increment'`.
schedule:
`'synchronous'`, `'sequential'`, `'random-sequential'`, or `'member'` for the agent form.
form:
`'symbolic'` for update laws of the states of the target, or `'agent'` for a rule over an agent group.
phase:
The phase of the step the agent form runs in, or `None` for the default phase.
Raises:
ValidationError:
The mode, the schedule, or the form is unknown, the symbolic form has the schedule `'member'`, the
function is not callable, or the agent form returns a value that does not match the state, raised by
`run` of the simulation.
Examples:
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.interactions import MapInteraction
>>> from skcomplex.core.spaces import DiscreteTime
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> from typing import Annotated, SupportsFloat
>>> class ValueElement(Element):
... def __init__(self, label, x: Annotated[float, StateSpec()] = 0.5):
... super().__init__(label)
... self.x = x
>>> logistic = MapInteraction(
... 'map', target='value', state='x', function=lambda s: 3 * s.target.x * (1 - s.target.x)
... )
>>> system = System(elements=[ValueElement('value')], interactions=[logistic], time=DiscreteTime())
>>> Simulation(system).check().prepare().run(until=2).trajectory_.values('value.x').tolist()
[0.5, 0.75, 0.5625]
>>> agents = MapInteraction('map', target='walkers', state='x', function=lambda a, m: a.x + 1, form='agent')
>>> type(agents.check().prepare().run(SimpleNamespace()).laws_[0]).__name__
'RuleLaw'
>>> from skcomplex.core.elements import GroupElement
>>> class WeightsElement(Element):
... def __init__(self, label, visual=None, *, w: Annotated[object, StateSpec(vector=True)] = 0.0):
... super().__init__(label, visual=visual)
... self.w = w
>>> group = GroupElement('weights', element=WeightsElement, size=1, values={'w': [1.0, 2.0]})
>>> doubling = MapInteraction(
... 'double', target='weights', state='w', function=lambda a, m: [2 * x for x in a.w], form='agent'
... )
>>> trajectory = Simulation(System(elements=[group], interactions=[doubling])).check().prepare().run(until=2)
>>> trajectory.trajectory_.agents('weights', 2).values['w'].tolist()
[[4.0, 8.0]]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str,
function: Callable[..., object],
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
mode: str = 'replace',
schedule: str = 'synchronous',
form: str = 'symbolic',
phase: str | None = None,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.function = function
self.sources = sources
self.mode = mode
self.schedule = schedule
self.form = form
self.phase = phase
def _check(self) -> None:
"""Raise an error for an unknown mode, schedule, or form, or a function that cannot be called.
Raises:
ValidationError:
The mode, the schedule, or the form is unknown, the symbolic form has the schedule `'member'`, or the
function is not callable.
"""
super()._check()
self._check_variant('mode', self.mode, MODES, (BaseSpec,))
for name, value, variants in (
('schedule', self.schedule, SCHEDULES),
('form', self.form, FORMS),
):
if value not in variants:
error_msg = f'The {name} of map {self.label!r} is {value!r}. Use one of {", ".join(variants)}.'
raise ValidationError(error_msg)
if self.schedule == 'member' and self.form != 'agent':
error_msg = (
f'The map {self.label!r} updates every member in its activation, which only the agent form does.'
)
raise ValidationError(error_msg)
if not callable(self.function):
error_msg = f'The function of map {self.label!r} must be callable, but it is {self.function!r}.'
raise ValidationError(error_msg)
def _apply_mode(self, current: Value, value: Value) -> Value:
"""Return the next value of a state from its current value and the value of the function, by the mode."""
if self.mode == 'increment':
return current + value
if self.mode == 'replace':
return value
raise self._get_variant_error(self.mode, 'mode', '_apply_mode')
def _get_vector(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> Value:
"""Return the next value of a vector state of an agent."""
current = np.asarray(getattr(agent, self.state), dtype=np.float64)
value = np.asarray(self.function(agent, model), dtype=np.float64)
if value.shape != current.shape:
error_msg = (
f'The function of map {self.label!r} returns the shape {value.shape} for the state {self.state!r} of '
f'the shape {current.shape}. Return an array of the shape of the state.'
)
raise ValidationError(error_msg)
return self._apply_mode(current, value)
def _get_value(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> Value:
"""Return the next value of the state of an agent, a number or the array of a vector state."""
if np.ndim(getattr(agent, self.state)) > 0:
return self._get_vector(agent, model)
found = self.function(agent, model)
if not isinstance(found, SupportsFloat):
error_msg = f'The function of map {self.label!r} returns {found!r} for an agent. Return a number.'
raise ValidationError(error_msg)
value = float(found)
return self._apply_mode(float(getattr(agent, self.state)), value)
def _set_value(self, agent: AgentViewProtocol, value: Value) -> None:
"""Set the state of an agent to its next value."""
setattr(agent, self.state, value)
def _update_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Set the state of one member to its next value."""
if agent.alive:
self._set_value(agent, self._get_value(agent, model))
def update(self, model: ModelViewProtocol) -> None:
"""Set the state of every agent of the target to its next value, by the schedule.
Args:
model:
The model view of the interaction.
"""
self._update_agents(model, self.schedule, self._get_value, self._set_value)
def _get_order(self) -> object:
"""Return the order of the members of the symbolic form: none, their own order, or a uniform draw."""
if self.schedule == 'synchronous':
return None
if self.schedule == 'sequential':
return SR(0)
return uniform(re.sub(r'\W', '_', f'order_{self.label}'))
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the update of the state to its next value, or the rule of the agent form.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target` and `namespace.sources`.
Returns:
laws:
The update law of the state, or the rule that runs the map of the agents.
"""
if self.form == 'agent' and self.schedule == 'member':
return [RuleLaw(self._update_member, phase=self.phase, members=self._get_role('target'))]
if self.form == 'agent':
return [RuleLaw(self.update, phase=self.phase)]
value = getattr(self._get_member(namespace, 'target'), self.state)
mapped = cast(Value, self.function(namespace))
next_value = self._apply_mode(value, mapped)
return [UpdateLaw(value, next_value, increment=self.mode == 'increment', order=self._get_order())]
class LifeCycleInteraction(Interaction):
"""Life cycle of the agents of a group, which give birth to offspring and die, every step of discrete time.
Every step the agents act once each, in their order under the schedule `'sequential'`, in a new random order under
`'random-sequential'`, or each in its own activation under `'member'`, together with the rules of its element and
the other member interactions, in the order the system lists them, without the reports of the step and without
arrivals. An agent dies with the probability `death`, and an agent that survives gives birth with the probability
`birth`, or as many times as a Poisson number of the mean `rate`. A probability is a number or a function of the
agent and the model view, and a function that returns a condition gives probability one where it holds. Every birth
creates `offspring` new agents, a number or a function of the parent and the model view, such as a Poisson draw, of
the group `group`, or of the parent's group without one. A new agent stands in the cell of its parent, with the
placement `'here'`, in a neighbouring cell that no agent of the group holds, with the placement `'neighbor'`, or at
the location a function of the parent and the model view returns, where `None` means no place, and a birth without a
place creates fewer agents.
A new agent copies the states `inherit` of its parent and takes the values that `values(parent, model)` returns,
such as a generation one higher or drawn states. The parent splits each state of `split` with its offspring,
keeping the share `keep`, a number or a function of the parent and the model view such as a random cut, or else an
equal share, and the offspring share the rest equally. Then `mutate(values, model)` returns the values of a new
agent from those it would get, and the other states take their defaults. A parent that gave birth takes the values
`reset`, and with `replace` it is removed, so its offspring replace it, also when the birth created none. With
`renew` an agent that dies is renewed in place instead of removed: it keeps its identity, its place, and the states
the new values leave out, and takes the values of an offspring of its own. With the probability `cull` an agent
that does not die removes, after its births, one member of `group`, or of the target without a group, for which
`kin(agent, member)` holds, the first in the order of the group under the pick `'first'`, such as one it owns, or
a random one under `'random'`. Every
step a Poisson number of the mean `arrival`, a number or a function of the model view, of new agents arrive without
a parent, at the location of the placement function or else at a random place, with the values
`values(None, model)`. The life cycle reports the `births`, `deaths`, and `arrivals` of every step. It is an agent
rule, so the system runs in discrete time.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the variable group whose agents give birth and die.
birth:
The probability that an agent gives birth in a step, or `None` for no births by probability.
death:
The probability that an agent dies in a step, or `None` for no deaths.
rate:
The mean number of births of an agent in a step, a number or a function, or `None` for no births by rate.
offspring:
The number of new agents of a birth.
placement:
Where a new agent stands, `'here'`, `'neighbor'`, or a function of the parent, `None` for an arrival, and
the model view that returns its location.
inherit:
The names of the states a new agent copies from its parent.
split:
The names of the states the parent shares with its offspring.
mutate:
A function of the values of a new agent and the model view that returns its values, or `None`.
schedule:
The order of the agents, `'sequential'`, `'random-sequential'`, or `'member'`.
phase:
The phase of the step the life cycle runs in, or `None` for the default phase.
group:
The path of the variable group the new agents join, or `None` for the group of the target.
values:
A function of the parent, `None` for an arrival, and the model view that returns values of a new agent by
state name, or `None`.
keep:
The share of every split state the parent keeps, a number or a function, or `None` for an equal share.
reset:
The values a parent takes after it gives birth, by state name, or `None`.
replace:
Whether a parent that gives birth is removed.
renew:
Whether an agent that dies is renewed in place instead of removed.
cull:
The probability that an agent removes one member of `group`, or of the target without a group, in a
step, or `None` for no removals.
kin:
A function of an agent and a member of `group` that returns whether the agent may remove the member, or
`None` for every member.
arrival:
The mean number of new agents that arrive in a step without a parent, or `None` for no arrivals.
pick:
Which member of `group` a removal takes among those it may remove, `'first'` or `'random'`.
Raises:
ValidationError:
The placement, the schedule, or the pick is unknown, both a birth probability and a rate are given, a
function is not callable, or the schedule `'member'` has arrivals, raised by `check`.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element, GroupElement
>>> from skcomplex.core.interactions import LifeCycleInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, visual=None, *, size: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label, visual=visual)
... self.size = size
>>> cells = GroupElement('cells', element=CellElement, size=1, values={'size': 8.0}, variable=True)
>>> division = LifeCycleInteraction('division', target='cells', birth=1.0, split=('size',))
>>> system = System(elements=[cells], interactions=[division], space=LatticeSpace((3, 3)))
>>> trajectory = Simulation(system, random_state=0).check().prepare().run(until=2).trajectory_
>>> trajectory.agents('cells', 2).values['size'].tolist()
[2.0, 2.0, 2.0, 2.0]
"""
reports: ClassVar[tuple[str, ...]] = ('births', 'deaths', 'arrivals')
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
birth: Chance | None = None,
death: Chance | None = None,
rate: Chance | None = None,
offspring: Count = 1,
placement: Placement = 'here',
inherit: Sequence[str] = (),
split: Sequence[str] = (),
mutate: Mutation | None = None,
schedule: str = 'random-sequential',
phase: str | None = None,
group: Annotated[str | None, RoleSpec(optional=True)] = None,
values: Values | None = None,
keep: Chance | None = None,
reset: Mapping[str, object] | None = None,
replace: bool = False,
renew: bool = False,
cull: Chance | None = None,
kin: Kinship | None = None,
arrival: Arrival | None = None,
pick: str = 'first',
) -> None:
super().__init__(label)
self.target = target
self.birth = birth
self.death = death
self.rate = rate
self.offspring = offspring
self.placement = placement
self.inherit = inherit
self.split = split
self.mutate = mutate
self.schedule = schedule
self.phase = phase
self.group = group
self.values = values
self.keep = keep
self.reset = reset
self.replace = replace
self.renew = renew
self.cull = cull
self.kin = kin
self.arrival = arrival
self.pick = pick
def _check_variants(self) -> None:
"""Raise an error for an unknown placement, schedule, or pick, or arrivals of a member schedule."""
self._check_variant('placement', self.placement, PLACEMENTS, functions=True)
self._check_variant('schedule', self.schedule, SCHEDULES)
self._check_variant('pick', self.pick, PICKS)
if self.schedule == 'member' and self.arrival is not None:
error_msg = f'{self.label!r} runs member by member, so no agent brings arrivals. Give `arrival=None`.'
raise ValidationError(error_msg)
def _check_space(self, space: BaseSpace | None) -> None:
"""Raise an error for a placement in a neighbouring cell or node in a continuous space."""
if self.placement == 'neighbor':
self._check_grid('placement', self.placement, space, "placement='here', or a placement function")
def _check(self) -> None:
"""Raise an error for an unknown placement, a birth by both probability and rate, or a bad function.
Raises:
ValidationError:
The placement, the schedule, or the pick is unknown, both a birth probability and a rate are given,
a value is neither a number nor a function, or the schedule `'member'` has arrivals.
"""
super()._check()
self._check_variants()
if self.birth is not None and self.rate is not None:
error_msg = f'{self.label!r} gives births by a probability and by a rate. Give one of them.'
raise ValidationError(error_msg)
for name in ('birth', 'death', 'rate', 'offspring', 'keep', 'cull', 'arrival'):
value = getattr(self, name)
if value is not None and not callable(value) and not isinstance(value, int | float):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a number or a function.'
raise ValidationError(error_msg)
for name in ('mutate', 'values', 'kin'):
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, or None.'
raise ValidationError(error_msg)
@staticmethod
def _evaluate(value: Chance | Count, agent: AgentViewProtocol, model: ModelViewProtocol) -> float:
"""Return a number, or the value of a function of an agent and the model view, as a float."""
return float(value(agent, model) if callable(value) else value)
def _happens(self, chance: Chance | None, agent: AgentViewProtocol, model: ModelViewProtocol) -> bool:
"""Return whether an event of a probability happens to an agent in this step."""
if chance is None:
return False
probability = self._evaluate(chance, agent, model)
return probability >= 1.0 or (probability > 0.0 and model.rng.random() < probability)
def _count_births(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> int:
"""Return the number of births of a parent in this step."""
if self.rate is not None:
return int(model.rng.poisson(self._evaluate(self.rate, agent, model)))
return int(self._happens(self.birth, agent, model))
def _get_places(
self,
agent: AgentViewProtocol,
count: int,
occupied: set[object],
model: ModelViewProtocol,
) -> list[object]:
"""Return the locations of the new agents of a parent, fewer when places run out."""
if callable(self.placement):
places = [self.placement(agent, model) for _ in range(count)]
return [place for place in places if place is not None]
if self.placement == 'here':
return [agent.location] * count
free = [cell for cell in agent.neighbor_cells() if cell not in occupied]
chosen = [free[int(index)] for index in model.rng.permutation(len(free))[:count]]
occupied.update(chosen)
return chosen
def _get_values(self, agent: AgentViewProtocol | None, model: ModelViewProtocol) -> dict[str, object]:
"""Return the inherited and given values of a new agent of a parent, or of an arrival."""
values: dict[str, object] = {} if agent is None else {name: getattr(agent, name) for name in self.inherit}
if self.values is not None:
values.update(self.values(agent, model))
return values
def _create(
self,
agent: AgentViewProtocol | None,
place: object,
values: dict[str, object],
model: ModelViewProtocol,
) -> None:
"""Create a new agent with mutated values at a place, `None` for a random one, in the group of the cycle."""
if self.mutate is not None:
values = self.mutate(values, model)
if self.group is not None:
model.create(self._get_role('group'), place, **values)
elif agent is not None:
agent.create(None, place, **values)
else:
model.create(self._get_role('target'), place, **values)
def _bear(self, agent: AgentViewProtocol, places: list[object], model: ModelViewProtocol) -> None:
"""Create the new agents of a parent at their places, with inherited, given, split, and mutated values."""
kept = 1 / (len(places) + 1) if self.keep is None else self._evaluate(self.keep, agent, model)
totals = {name: float(getattr(agent, name)) for name in self.split}
for name, amount in totals.items():
setattr(agent, name, amount * kept)
shares = {name: (amount - amount * kept) / max(len(places), 1) for name, amount in totals.items()}
for place in places:
values = self._get_values(agent, model)
values.update(shares)
self._create(agent, place, values, model)
def _renew(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Give an agent that dies the values of an offspring of its own, in place."""
values = self._get_values(agent, model)
if self.mutate is not None:
values = self.mutate(values, model)
for name, value in values.items():
setattr(agent, name, value)
def _remove_kin(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Remove the first or a random member of the group of the cycle that is kin of an agent, if there is one."""
kin = self.kin
group = self._get_role('target' if self.group is None else 'group')
members = (item for item in model.agents(group) if item.alive and (kin is None or kin(agent, item)))
if self.pick == 'random':
candidates = list(members)
if candidates:
candidates[int(model.rng.integers(len(candidates)))].remove()
return
member = next(members, None)
if member is not None:
member.remove()
def _die(self, agent: AgentViewProtocol, occupied: set[object], model: ModelViewProtocol) -> None:
"""Renew an agent in place, or remove it and free its cell."""
if self.renew:
self._renew(agent, model)
return
occupied.discard(agent.location_key)
agent.remove()
def _act(self, agent: AgentViewProtocol, occupied: set[object], model: ModelViewProtocol) -> tuple[int, int]:
"""Let an agent die, or give birth and remove kin, and return its numbers of new agents and deaths."""
if self._happens(self.death, agent, model):
self._die(agent, occupied, model)
return 0, 1
born, died = self._give_birth(agent, occupied, model)
if self._happens(self.cull, agent, model):
self._remove_kin(agent, model)
return born, died
def _give_birth(self, agent: AgentViewProtocol, occupied: set[object], model: ModelViewProtocol) -> tuple[int, int]:
"""Let an agent give birth, and return its number of new agents and whether it was replaced."""
born = 0
events = self._count_births(agent, model)
if events > 0:
births = events * int(self._evaluate(self.offspring, agent, model))
places = self._get_places(agent, births, occupied, model) if births else []
if places or self.split:
self._bear(agent, places, model)
born = len(places)
for name, value in (self.reset or {}).items():
setattr(agent, name, value)
if self.replace:
occupied.discard(agent.location_key)
agent.remove()
return born, 1
return born, 0
def _arrive(self, model: ModelViewProtocol) -> int:
"""Create the agents that arrive without a parent in this step, and return their number."""
if self.arrival is None:
return 0
mean = self.arrival(model) if callable(self.arrival) else self.arrival
count = int(model.rng.poisson(float(mean)))
for _ in range(count):
place = self.placement(None, model) if callable(self.placement) else None
self._create(None, place, self._get_values(None, model), model)
return count
def _live_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Run the life cycle of one agent in its own activation."""
if not agent.alive:
return
neighbor = self.placement == 'neighbor'
occupied = {other.location_key for other in model.agents(self._get_role('target'))} if neighbor else set()
self._act(agent, occupied, model)
def _live(self, model: ModelViewProtocol) -> None:
"""Remove or renew the agents that die and create the offspring of those that give birth, by the schedule."""
agents = model.agents(self._get_role('target'))
occupied = {agent.location_key for agent in agents} if self.placement == 'neighbor' else set()
births = deaths = 0
order = model.rng.permutation(len(agents)) if self.schedule == 'random-sequential' else range(len(agents))
for index in order:
agent = agents[int(index)]
if not agent.alive:
continue
born, died = self._act(agent, occupied, model)
births, deaths = births + born, deaths + died
model.report('births', births)
model.report('deaths', deaths)
model.report('arrivals', self._arrive(model))
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the rule of the life cycle.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rule that removes and creates agents every step.
"""
if self.schedule == 'member':
return [RuleLaw(self._live_member, phase=self.phase, members=self._get_role('target'))]
return [RuleLaw(self._live, phase=self.phase)]
class MovementInteraction(Interaction):
"""Movement of the agents of a group by steps, every step of discrete time.
A step is one of five kinds:
- `'neighbor'` moves an agent to a uniformly drawn cell or node next to its own, on a lattice or a network.
- `'gaussian'` adds to the position of an agent a normal vector whose components have the standard deviation
`scale`, in a continuous space.
- `'levy'` moves an agent in a uniformly drawn direction by a Pareto length, whose least value is `scale` and whose
tail exponent is `exponent`, in a continuous space.
- `'best'` moves an agent to the cell or node next to its own with the highest `score(agent, cell, model)`. With a
`vision`, a number or a function of the agent and the model view, it looks instead along every direction of the
neighbourhood of a lattice up to that many cells, and at the cell of the agent, in the order of the axes, the
positive direction first. Among cells of equal highest score, the ties `'random'` draw one uniformly, `'first'`
take the first in that order, and `'nearest'` take the nearest, the first of them among equal distances. By
default the ties are `'random'` without a vision and `'nearest'` with one.
- A function `step(agent, model)` returns the destination, or `None` to stay.
An agent for which `condition(agent, model)` holds attempts one step every step of time, or a Poisson number of
steps with the mean `rate`, and takes each attempt with the probability `probability` and then with the probability
`accept(agent, destination, model)`. An exclusive movement never moves an agent to a cell another agent of the group
holds. The agents move one after the other, in their order under the schedule `'sequential'`, in a new random order
every step under `'random-sequential'`, or each in its own activation under `'member'`, together with the rules of
its element and the other member interactions, in the order the system lists them. Every step taken charges its
`cost`, a `TransferInteraction` from the agent to itself or a function of the agent and the model view. The
movement is an agent rule, so the system runs in discrete time.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group whose agents move.
step:
The kind of step, `'neighbor'`, `'gaussian'`, `'levy'`, or `'best'`, or a function of the agent and the
model view that returns the destination.
probability:
The probability that an agent takes an attempted step.
rate:
The mean number of attempted steps of an agent in one unit of time, or `None` for one attempt every step.
condition:
A function of the agent and the model view that returns whether the agent moves, or `None` for every agent.
accept:
A function of the agent, its destination, and the model view that returns the probability of taking the
step, or `None` to take every step.
exclusive:
Whether an agent moves only to cells or nodes no other agent of the group holds.
score:
The function of the agent, a neighbouring cell, and the model view that the step `'best'` maximizes.
scale:
The standard deviation of every component of a Gaussian step, or the least length of a Lévy step.
exponent:
The tail exponent of the length of a Lévy step.
cost:
A `TransferInteraction` that the agent pays to itself or to the environment for every step taken, a function
of the agent and the model view run for every step taken, or `None`.
vision:
The number of cells the step `'best'` looks along every direction of a lattice, a number or a function of
the agent and the model view, or `None` for the cells next to the agent.
schedule:
The order of the agents, `'sequential'`, `'random-sequential'`, or `'member'`.
phase:
The phase of the step the movement runs in, or `None` for the default phase.
ties:
The cell the step `'best'` takes among cells of equal highest score, `'random'`, `'first'`, or `'nearest'`,
or `None` for `'random'` without a vision and `'nearest'` with one.
Raises:
ValidationError:
The step, the schedule, or the ties are unknown, the rate is not positive, a function is not callable, or
the step `'best'` has no score, raised by `check`.
Examples:
>>> from skcomplex.core.elements import Element, GroupElement
>>> from skcomplex.core.interactions import MovementInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> walkers = GroupElement('walkers', element=Element, size=2, location=[(2, 2), (2, 2)])
>>> walk = MovementInteraction('walk', target='walkers', step='neighbor')
>>> space = LatticeSpace((5, 5), neighborhood='von_neumann')
>>> system = System(elements=[walkers], interactions=[walk], space=space)
>>> trajectory = Simulation(system, random_state=0).check().prepare().run(until=1).trajectory_
>>> [bool(abs(cell - 2).sum() == 1) for cell in trajectory.agents('walkers', 1).locations]
[True, True]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
step: Step = 'neighbor',
probability: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 1.0,
rate: float | None = None,
condition: Callable[[AgentViewProtocol, ModelViewProtocol], bool] | None = None,
accept: Callable[[AgentViewProtocol, object, ModelViewProtocol], float] | None = None,
exclusive: bool = False,
score: Callable[[AgentViewProtocol, object, ModelViewProtocol], float] | None = None,
scale: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
exponent: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.5,
cost: Cost | None = None,
vision: Vision | None = None,
schedule: str = 'random-sequential',
phase: str | None = None,
ties: str | None = None,
) -> None:
super().__init__(label)
self.target = target
self.step = step
self.probability = probability
self.rate = rate
self.condition = condition
self.accept = accept
self.exclusive = exclusive
self.score = score
self.scale = scale
self.exponent = exponent
self.cost = cost
self.vision = vision
self.schedule = schedule
self.phase = phase
self.ties = ties
def _check_functions(self) -> None:
"""Raise an error for a condition, an acceptance, a score, or a cost that cannot be called."""
for name in ('condition', 'accept', 'score'):
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, or None.'
raise ValidationError(error_msg)
if isinstance(self.cost, TransferInteraction):
self.cost.check()
elif self.cost is not None and not callable(self.cost):
error_msg = f'The cost of {self.label!r} is {self.cost!r}. Give a TransferInteraction, a function, or None.'
raise ValidationError(error_msg)
def _check_space(self, space: BaseSpace | None) -> None:
"""Raise an error for a step to a neighbouring cell or node in a continuous space."""
if isinstance(self.step, str) and self.step in GRID_STEPS:
self._check_grid('step', self.step, space, "step='gaussian' or step='levy', or a step function")
def _check(self) -> None:
"""Raise an error for an unknown step, a rate that is not positive, or a function that cannot be called.
Raises:
ValidationError:
The step, the schedule, or the ties are unknown, the rate is not positive, a function is not callable,
or the step `'best'` has no score.
"""
super()._check()
self._check_variant('step', self.step, STEPS, functions=True)
if self.rate is not None and (isinstance(self.rate, bool) or not self.rate > 0):
error_msg = f'The rate of {self.label!r} is {self.rate!r}. Give a positive number, or None.'
raise ValidationError(error_msg)
if self.schedule not in SCHEDULES:
error_msg = f'Unknown schedule {self.schedule!r} of {self.label!r}. Use one of {", ".join(SCHEDULES)}.'
raise ValidationError(error_msg)
if self.ties is not None:
self._check_variant('ties', self.ties, TIES)
if self.vision is not None and not callable(self.vision) and not isinstance(self.vision, int | float):
error_msg = f'The vision of {self.label!r} is {self.vision!r}. Give a number, a function, or None.'
raise ValidationError(error_msg)
if self.step == 'best' and self.score is None:
error_msg = f'The step `best` of {self.label!r} maximizes a score. Pass `score(agent, cell, model)`.'
raise ValidationError(error_msg)
self._check_functions()
def _count_attempts(self, model: ModelViewProtocol) -> int:
"""Return the number of steps an agent attempts in one step of time."""
return 1 if self.rate is None else int(model.rng.poisson(self.rate))
def _get_cells(self, agent: AgentViewProtocol, occupied: set[object]) -> list[object]:
"""Return the cells or nodes around an agent that it may move to."""
cells = agent.neighbor_cells()
if not self.exclusive:
return cells
return [cell for cell in cells if cell not in occupied]
def _get_best(
self,
agent: AgentViewProtocol,
cells: list[object],
model: ModelViewProtocol,
distances: list[int] | None = None,
) -> object:
"""Return the cell of the highest score, chosen among ties by the tie rule."""
if self.score is None or not cells:
return None
scores = [self.score(agent, cell, model) for cell in cells]
ties = self.ties or ('random' if distances is None else 'nearest')
if ties != 'nearest':
return cells[self._pick_best(scores, ties, model)]
values = np.asarray(scores, dtype=float)
best = np.flatnonzero(values == values.max())
near = np.asarray([1 if distances is None else distances[int(index)] for index in best])
return cells[int(best[near == near.min()][0])]
@staticmethod
def _shift(here: tuple[int, ...], cell: tuple[int, ...], model: ModelViewProtocol) -> tuple[int, ...]:
"""Return the unit step from a cell to a neighbouring cell, across the wrapped edge of a torus."""
shift = [b - a for a, b in zip(here, cell, strict=True)]
if model.torus:
sizes = [int(item) for item in model.dimensions or ()]
shift = [d - round(d / n) * n for d, n in zip(shift, sizes, strict=True)]
return tuple(shift)
def _get_sight(
self,
agent: AgentViewProtocol,
occupied: set[object],
model: ModelViewProtocol,
) -> tuple[list[object], list[int]]:
"""Return the cells in sight of an agent along every direction of a lattice, its own first, with distances."""
vision = self.vision(agent, model) if callable(self.vision) else self.vision
reach = int(vision or 0)
here = tuple(int(item) for item in np.atleast_1d(np.asarray(agent.location_key)))
sizes = [int(item) for item in model.dimensions or ()]
cells: list[object] = [agent.location_key]
distances = [0]
shifts = [
self._shift(here, tuple(int(item) for item in np.atleast_1d(np.asarray(neighbor))), model)
for neighbor in agent.neighbor_cells()
]
for shift in sorted(shifts, key=lambda step: tuple((item == 0, -item) for item in step)):
for distance in range(1, reach + 1):
point = [a + distance * d for a, d in zip(here, shift, strict=True)]
if model.torus:
point = [item % size for item, size in zip(point, sizes, strict=True)]
elif any(not 0 <= item < size for item, size in zip(point, sizes, strict=True)):
break
cell = tuple(point)
if not self.exclusive or cell not in occupied:
cells.append(cell)
distances.append(distance)
return cells, distances
def _get_destination(self, agent: AgentViewProtocol, model: ModelViewProtocol, occupied: set[object]) -> object:
"""Return the location a step of an agent leads to, or `None` to stay."""
rng = model.rng
if callable(self.step):
return self.step(agent, model)
if self.step == 'best' and self.vision is not None:
cells, distances = self._get_sight(agent, occupied, model)
return self._get_best(agent, cells, model, distances)
if self.step in ('neighbor', 'best'):
cells = self._get_cells(agent, occupied)
if self.step == 'best':
return self._get_best(agent, cells, model)
return cells[int(rng.integers(len(cells)))] if cells else None
position = np.asarray(agent.location, dtype=float)
if self.step == 'gaussian':
return position + rng.normal(0.0, self._get_parameter_value('scale', model), position.shape)
direction = rng.normal(size=position.shape)
length = self._get_parameter_value('scale', model) * (
1.0 + rng.pareto(self._get_parameter_value('exponent', model))
)
return position + length * direction / np.linalg.norm(direction)
def _charge(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Charge the cost of one step taken to an agent."""
if isinstance(self.cost, TransferInteraction):
self.cost.transfer(agent, agent, model)
elif self.cost is not None:
self.cost(agent, model)
def _take(self, agent: AgentViewProtocol, destination: object, occupied: set[object]) -> bool:
"""Move an agent to a destination unless an exclusive movement finds it held, and return whether it moved."""
if not self.exclusive:
agent.move_to(destination)
return True
if destination == agent.location_key:
return True
if destination in occupied:
return False
occupied.discard(agent.location_key)
agent.move_to(destination)
occupied.add(agent.location_key)
return True
def _move(
self,
agent: AgentViewProtocol,
model: ModelViewProtocol,
probability: float,
occupied: set[object],
) -> None:
"""Take the attempted steps of one agent, each with its probability and acceptance, charging every step."""
for _ in range(self._count_attempts(model)):
if probability < 1.0 and model.rng.random() >= probability:
continue
destination = self._get_destination(agent, model, occupied)
if destination is None:
continue
if self.accept is not None and model.rng.random() >= self.accept(agent, destination, model):
continue
if self._take(agent, destination, occupied):
self._charge(agent, model)
def _move_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Move one agent whose condition holds, in its own activation."""
if agent.alive and (self.condition is None or self.condition(agent, model)):
occupied = (
{other.location_key for other in model.agents(self._get_role('target'))} if self.exclusive else set()
)
self._move(agent, model, self._get_parameter_value('probability', model), occupied)
def _move_all(self, model: ModelViewProtocol) -> None:
"""Move every agent of the group whose condition holds, in their order or in random order."""
agents = model.agents(self._get_role('target'))
occupied = {agent.location_key for agent in agents} if self.exclusive else set()
probability = self._get_parameter_value('probability', model)
order = model.rng.permutation(len(agents)) if self.schedule == 'random-sequential' else range(len(agents))
for index in order:
agent = agents[int(index)]
if agent.alive and (self.condition is None or self.condition(agent, model)):
self._move(agent, model, probability, occupied)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the rule of the movement.
Args:
namespace:
The namespace of the symbols the laws can use.
Returns:
laws:
The rule that moves the agents every step.
"""
if self.schedule == 'member':
return [RuleLaw(self._move_member, phase=self.phase, members=self._get_role('target'))]
return [RuleLaw(self._move_all, phase=self.phase)]
class TumourCellElement(Element):
"""Tumour cell that consumes oxygen, divides when supplied, secretes VEGF when hypoxic, and dies without oxygen.
Every step the cell takes up `consumption` oxygen at its site, or what is left. Above `proliferation` it divides
with probability `division` into a free neighbouring cell. Below `hypoxia` it secretes `secretion` VEGF, and below
`necrosis` it dies. The model adds the uptake, the hypoxia, the secretion, and the life cycle as member
interactions.
Args:
label:
The label of the agent or group.
visual:
How the element is drawn, or `None` for `default_visual_`.
hypoxic:
Whether the cell is hypoxic, `1`, or not, `0`.
consumption:
The oxygen a cell takes up per step.
proliferation:
The oxygen above which a cell divides.
division:
The probability per step of a division.
hypoxia:
The oxygen below which a cell secretes VEGF.
secretion:
The VEGF a hypoxic cell secretes per step.
necrosis:
The oxygen below which a cell dies.
Examples:
>>> from skcomplex.core.elements import GroupElement
>>> from skcomplex.models.biology.oncology import TumourCellElement
>>> GroupElement('tumour', element=TumourCellElement, size=3, location=(5, 5)).size
3
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = VisualSpec('agents', {'color_by': 'hypoxic', 'color_range': (0.0, 1.0)}),
*,
hypoxic: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 0.0,
consumption: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
proliferation: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
division: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
hypoxia: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
secretion: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
necrosis: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.hypoxic = hypoxic
self.consumption = consumption
self.proliferation = proliferation
self.division = division
self.hypoxia = hypoxia
self.secretion = secretion
self.necrosis = necrosis
@classmethod
def _get_view(cls) -> type[TumourCellView]:
"""Return the view of a tumour cell, `TumourCellView`."""
return TumourCellView
class VesselElement(Element):
"""Segment of a blood vessel, which supplies oxygen to its site every step, a member transfer of the model.
Args:
label:
The label of the agent or group.
visual:
How the element is drawn, or `None` for `default_visual_`.
supply:
The oxygen supplied per step.
Examples:
>>> from skcomplex.core.elements import GroupElement
>>> from skcomplex.models.biology.oncology import VesselElement
>>> GroupElement('vessels', element=VesselElement, size=2, location=[(0, 0), (0, 1)]).size
2
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
supply: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.supply = supply
@classmethod
def _get_view(cls) -> type[VesselView]:
"""Return the view of a vessel, `VesselView`."""
return VesselView
class SproutElement(Element):
"""Endothelial tip cell that climbs the VEGF gradient and leaves a new vessel segment behind it.
Every step the tip moves to a neighbouring cell with a probability that grows with `exp(chemotaxis * VEGF /
VEGF_here)`, so it senses the relative gradient, and it leaves a vessel segment where there was none. With
probability `branching` the tip splits in two, the new tip standing where the tip was. Below the VEGF level
`activation` the tip stays put. The model adds the vessel segment, the branching, and the step as member
interactions.
Args:
label:
The label of the agent or group.
visual:
How the element is drawn, or `None` for `default_visual_`.
chemotaxis:
How strongly VEGF draws the tip.
branching:
The probability per step of a new tip.
activation:
The VEGF above which the tip moves.
Examples:
>>> from skcomplex.core.elements import GroupElement
>>> from skcomplex.models.biology.oncology import SproutElement
>>> GroupElement('sprouts', element=SproutElement, size=2, location=(0, 5)).size
2
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
chemotaxis: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 30.0,
branching: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
activation: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
) -> None:
super().__init__(label, visual=visual)
self.chemotaxis = chemotaxis
self.branching = branching
self.activation = activation
@classmethod
def _get_view(cls) -> type[SproutView]:
"""Return the view of a sprout tip, `SproutView`."""
return SproutView
class TissueElement(Element):
"""Tissue in which oxygen and the vascular endothelial growth factor VEGF diffuse and decay.
Both follow `u' = D Δu - lambda u`, and cells add or remove amounts at their sites: vessels supply oxygen, tumour
cells consume it, and hypoxic tumour cells secrete VEGF. The model adds the diffusion and the decay as
interactions.
Args:
label:
The label of the tissue or the lattice field.
visual:
How the element is drawn, or `None` for `default_visual_`.
oxygen:
The initial oxygen.
vegf:
The initial VEGF.
D_oxygen:
The diffusion of oxygen.
decay_oxygen:
The decay rate of oxygen.
D_vegf:
The diffusion of VEGF.
decay_vegf:
The decay rate of VEGF.
Examples:
>>> from skcomplex.models.biology.oncology import TissueElement
>>> TissueElement('tissue').D_vegf
1.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
oxygen: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
vegf: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
D_oxygen: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
decay_oxygen: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
D_vegf: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
decay_vegf: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.005,
) -> None:
super().__init__(label, visual=visual)
self.oxygen = oxygen
self.vegf = vegf
self.D_oxygen = D_oxygen
self.decay_oxygen = decay_oxygen
self.D_vegf = D_vegf
self.decay_vegf = decay_vegf
class DiffusionInteraction(ConservationInteraction):
"""Diffusion `du/dt = D laplacian(u)` of a state `u` of a field, the conservation law of the flux `-D gradient(u)`.
The coefficient `D` is uniform in space, a number or a function of the namespace that returns an expression of
parameters and time. The interaction lowers to the Laplacian of the field, which its discretization evaluates, so
the total of the state on a periodic field stays constant and a point mass spreads as the heat kernel.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state diffuses.
state:
The name of the state of the field that diffuses.
coefficient:
A function of the namespace that returns the diffusion coefficient, or a number.
Raises:
ValidationError:
The state is not a name, or the coefficient is neither callable nor a number.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import DiffusionInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((40,), spacing=0.1, boundary='periodic')
>>> field = FieldElement('rod', element=CellElement, space=space, values={'u': np.eye(40)[20]})
>>> heat = DiffusionInteraction('heat', target='rod', state='u', coefficient=0.1)
>>> trajectory = Simulation(System(elements=[field], interactions=[heat]), dt=0.01).check().prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('rod.u')[-1].sum()), 6)
1.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
coefficient: Callable[[SimpleNamespace], object] | float = 1.0,
) -> None:
super().__init__(label, target=target, state=state)
self.coefficient = coefficient
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a coefficient that is neither callable nor a number.
Raises:
ValidationError:
The type declares states, the state is not a name, or the coefficient is neither callable nor a number.
"""
super()._check()
if not (callable(self.coefficient) or isinstance(self.coefficient, int | float)):
error_msg = f'The coefficient of {self.label!r} is {self.coefficient!r}. Give a function or a number.'
raise ValidationError(error_msg)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate `D laplacian(u)` of the state, the divergence of the diffusive flux.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The rate law of the state.
"""
coefficient = self.coefficient(namespace) if callable(self.coefficient) else self.coefficient
state = getattr(self._get_member(namespace, 'target'), self.state)
return [RateLaw(state, coefficient * laplacian(state))]
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