Food security models
A food security system, the rice self-sufficiency of a country under subsidy, research, and extension.
- Name:
food_security - Field:
economics/resources - Formalism: ordinary differential equations
- References: Bala et al. 2017
Factory
@register_model(
field='economics/resources',
tags=('food security', 'agriculture', 'stocks and flows', 'system dynamics', 'ode'),
references=('Bala et al. 2017',),
defaults={'until': 50.0, 'dt': 0.25},
observables=lambda system: [
ObservableSpec(name, expression, element='rice')
for name, expression in (
('input_ratio', _input_ratio),
('productivity', _productivity),
('production', _production),
('requirement', _requirement),
('self_sufficiency', _self_sufficiency),
)
],
)
def food_security(
*,
subsidy: float = 0.67,
bio_target: float = 0.0,
research: float = 0.0,
training_target: float = 0.0,
intensity_target: float = 1.0,
) -> System:
"""Return a food security system, the rice self-sufficiency of a country under subsidy, research, and extension.
The land suitable for rice is developed into rice area over eight years, and the rice area is discarded for
housing, roads, and industry at a small fraction per year. The production is the milled share of the
productivity times the area times the cropping intensity. The productivity is the potential yield, which research
raises towards a limit, times a table of the input ratio, the soil quality, which chemical fertiliser degrades and
fallow regenerates, and the farming practice, which training and extension improve through a learning table. The
input ratio is the share of the recommended fertiliser that the farmers apply, which the subsidy raises, and a
target share of bio-fertiliser replaces the chemical one without degrading the soil. The requirement is the per
capita consumption of a population growing exponentially, and the self-sufficiency level is the production over
the requirement. Land is conserved, so the rice area peaks once the suitable land is developed and then declines,
and without research and extension the self-sufficiency falls as the population grows. The structure and the
initial values follow the food security model of Malaysia of Bala, Arshad and Noh (2017, Chapter 10), and the
table functions and the rates are illustrative, with the land in millions of hectares, the population in
millions, the productivity in tonnes of paddy per hectare, the rice in millions of tonnes, and time in years.
Args:
subsidy:
The share of the recommended fertiliser that the subsidy pays for.
bio_target:
The target share of bio-fertiliser in the inputs.
research:
The funding of research for high-yielding varieties, relative to its normal level.
training_target:
The target coverage of the farmers by training and extension, from 0 to 1.
intensity_target:
The target cropping intensity, the crops per year on the rice area.
Returns:
system:
The model, with the land stocks `rice.available_land` and `rice.rice_area`, the population
`rice.population`, the factors `rice.potential`, `rice.soil_quality`, `rice.practice`, `rice.training`,
`rice.bio_fraction`, and `rice.cropping_intensity`, and the flows `rice.development`, `rice.discard`, and
`rice.population_growth` per year, over years. The observables of the model are `rice.input_ratio`,
`rice.productivity`, `rice.production`, `rice.requirement`, and `rice.self_sufficiency`.
Examples:
>>> from skcomplex.models.economics.resources import food_security
>>> food_security(subsidy=0.0).elements[0].subsidy
0.0
"""
sector = RiceSectorElement(
'rice',
subsidy=subsidy,
bio_target=bio_target,
research=research,
training_target=training_target,
intensity_target=intensity_target,
)
transfers = [
TransferInteraction(
'land_use',
source='rice',
target='rice',
state='available_land',
into='rice_area',
amount=lambda s: s.source.available_land / s.source.conversion_time,
flow='development',
legs=[('rice_area', Sink('urban'), lambda s: s.source.discard_rate * s.source.rice_area, 'discard')],
),
TransferInteraction(
'births',
source='rice',
target='rice',
state=Source('births'),
into='population',
amount=lambda s: s.source.growth_rate * s.source.population,
flow='population_growth',
),
]
relaxations = [
RelaxationInteraction(
'training',
target='rice',
state='training',
goal=lambda s: s.target.training_target,
time=lambda s: s.target.extension_time,
),
RelaxationInteraction(
'bio_fraction',
target='rice',
state='bio_fraction',
goal=lambda s: s.target.bio_target,
time=lambda s: s.target.transition_time,
),
RelaxationInteraction(
'cropping_intensity',
target='rice',
state='cropping_intensity',
goal=lambda s: s.target.intensity_target,
time=lambda s: s.target.implementation_time,
),
]
return System(elements=[sector], interactions=[*transfers, *relaxations])
Components
class RiceSectorElement(Element):
"""Rice sector of a country, its land, population, and the factors of the productivity of its rice land.
The suitable land is developed into rice area and the rice area is discarded, and the population grows
exponentially, through transfers within the sector of the model. Research raises the potential yield towards a
limit, chemical fertiliser degrades the soil and fallow regenerates it, and training improves the farming
practice through a learning table. The training, the share of bio-fertiliser, and the cropping intensity relax
towards their targets through relaxation interactions of the model, which read the targets and times of the
sector.
Args:
label:
The label of the sector.
visual:
How the element is drawn, or `None` for `default_visual_`.
available_land:
The initial land suitable for rice and not yet developed.
rice_area:
The initial rice area.
population:
The initial population.
potential:
The initial potential yield.
soil_quality:
The initial soil quality, from 0 to 1.
practice:
The initial quality of the farming practice, from 0 to 1.
training:
The initial coverage of the farmers by training and extension, from 0 to 1.
bio_fraction:
The initial share of bio-fertiliser in the inputs.
cropping_intensity:
The initial number of crops per year on the rice area.
conversion_time:
The time to develop the suitable land.
discard_rate:
The fraction of the rice area discarded per year.
growth_rate:
The growth rate of the population.
consumption:
The rice consumption per capita, in kilograms per year.
milling_rate:
The share of paddy that milling turns into rice.
subsidy:
The share of the recommended fertiliser that the subsidy pays for.
farmer_inputs:
The share of the unsubsidised fertiliser that the farmers buy.
research:
The funding of research for high-yielding varieties, relative to its normal level.
potential_limit:
The potential yield that research approaches.
research_time:
The time research takes to reach the limit at its normal funding.
regeneration_time:
The time fallow takes to regenerate the soil.
degradation:
The rate at which chemical fertiliser degrades the soil.
learning_time:
The time training takes to improve the practice.
training_target:
The target coverage of training and extension.
extension_time:
The time the training takes to reach its target.
bio_target:
The target share of bio-fertiliser.
transition_time:
The time the share of bio-fertiliser takes to reach its target.
intensity_target:
The target cropping intensity.
implementation_time:
The time the cropping intensity takes to reach its target.
Examples:
>>> from skcomplex.models.economics.resources import RiceSectorElement
>>> RiceSectorElement('rice', subsidy=0.5).subsidy
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
available_land: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo), quantity='land')] = 0.088,
rice_area: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo), quantity='land')] = 0.680647,
population: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo), quantity='people')] = 18.1024,
potential: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 6.0,
soil_quality: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.99,
practice: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.61,
training: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
bio_fraction: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
cropping_intensity: Annotated[object, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
conversion_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 8.0,
discard_rate: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0085,
growth_rate: Annotated[object, ParameterSpec()] = 0.021,
consumption: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 80.0,
milling_rate: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.65,
subsidy: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.67,
farmer_inputs: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
research: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
potential_limit: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
research_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
regeneration_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
degradation: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.002,
learning_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 5.0,
training_target: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
extension_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 5.0,
bio_target: Annotated[object, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
transition_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
intensity_target: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
implementation_time: Annotated[object, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
) -> None:
super().__init__(label, visual=visual)
self.available_land = available_land
self.rice_area = rice_area
self.population = population
self.potential = potential
self.soil_quality = soil_quality
self.practice = practice
self.training = training
self.bio_fraction = bio_fraction
self.cropping_intensity = cropping_intensity
self.conversion_time = conversion_time
self.discard_rate = discard_rate
self.growth_rate = growth_rate
self.consumption = consumption
self.milling_rate = milling_rate
self.subsidy = subsidy
self.farmer_inputs = farmer_inputs
self.research = research
self.potential_limit = potential_limit
self.research_time = research_time
self.regeneration_time = regeneration_time
self.degradation = degradation
self.learning_time = learning_time
self.training_target = training_target
self.extension_time = extension_time
self.bio_target = bio_target
self.transition_time = transition_time
self.intensity_target = intensity_target
self.implementation_time = implementation_time
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the potential, soil, and practice of the sector.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the element.
"""
return _rice_laws(s)
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 RelaxationInteraction(RateInteraction):
"""Relaxation `(goal(ns) - x) / time(ns)` of a state `x` of a target towards a goal, with optional noise.
With a constant goal and time and an additive noise of amplitude `sigma`, the state follows the Ornstein-Uhlenbeck
process `dx = (goal - x) / time dt + sigma dW`, whose stationary variance is `sigma**2 time / 2`. Without noise it
is first-order smoothing, or adaptive expectation, towards the goal. A goal or a time function takes the namespace
of the interaction, as for `RateInteraction`, and a number stands for a constant. A goal that is a number applies
to every component of a vector state.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state relaxes.
state:
The name of the state of the target that relaxes.
goal:
A function of the namespace that returns the goal, or a number.
time:
A function of the namespace that returns the relaxation time, or a positive number.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
Raises:
ValidationError:
The state is not a name, or a function is neither callable nor a number.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RelaxationInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> smoothing = RelaxationInteraction('smoothing', target='belief', state='n', goal=10.0, time=0.5)
>>> system = System(elements=[CompartmentElement('belief', n=0.0)], interactions=[smoothing])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=5.0).trajectory_
>>> round(float(trajectory.values('belief.n')[-1]), 3)
10.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
goal: Callable[[SimpleNamespace], object] | float = 0.0,
time: Callable[[SimpleNamespace], object] | float = 1.0,
noise: Callable[[SimpleNamespace], object] | float | None = None,
) -> None:
super().__init__(label, target=target, state=state, noise=noise)
self.goal = goal
self.time = time
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."""
super()._check_functions()
for name in ('goal', 'time'):
value = getattr(self, name)
if not (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)
def _get_goal(self, namespace: SimpleNamespace) -> object:
"""Return the goal, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.goal, namespace)
def _get_time(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation time, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.time, namespace)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation `(goal - x) / time` of the state.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The rate of the state.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
goal = self._broadcast(self._get_goal(namespace), state)
return (goal - state) / self._get_time(namespace)
Equations
rice.available_land: diff(available_land_rice(t), t) == -available_land_rice(t)/conversion_time_rice
rice.rice_area: diff(rice_area_rice(t), t) == -discard_rate_rice*rice_area_rice(t) + available_land_rice(t)/conversion_time_rice
rice.population: diff(population_rice(t), t) == growth_rate_rice*population_rice(t)
rice.potential: diff(potential_rice(t), t) == (potential_limit_rice - potential_rice(t))*research_rice/research_time_rice
rice.soil_quality: diff(soil_quality_rice(t), t) == ((farmer_inputs_rice*(subsidy_rice - 1) - subsidy_rice)*(bio_fraction_rice(t) - 1) + bio_fraction_rice(t))*degradation_rice*(bio_fraction_rice(t) - 1)*cropping_intensity_rice(t)*soil_quality_rice(t) - (soil_quality_rice(t) - 1)/regeneration_time_rice
rice.practice: diff(practice_rice(t), t) == -(practice_rice(t) - 1)*where(training_rice(t) < 0.0, 0.0, where(training_rice(t) < 0.5, where(training_rice(t) < 0.25, 1.6*training_rice(t), 1.1999999999999997*training_rice(t) + 0.10000000000000009), where(training_rice(t) < 0.75, 0.8000000000000003*training_rice(t) + 0.2999999999999998, where(training_rice(t) < 1.0, 0.3999999999999999*training_rice(t) + 0.6000000000000001, 1.0))))/learning_time_rice
rice.training: diff(training_rice(t), t) == (training_target_rice - training_rice(t))/extension_time_rice
rice.bio_fraction: diff(bio_fraction_rice(t), t) == (bio_target_rice - bio_fraction_rice(t))/transition_time_rice
rice.cropping_intensity: diff(cropping_intensity_rice(t), t) == (intensity_target_rice - cropping_intensity_rice(t))/implementation_time_rice