Forest dieback with replicator dynamics of mitigators
The forest dieback model coupled with the replicator dynamics of mitigators who perceive its warming.
- Name:
social_forest_dieback - Field:
physics/climate - Formalism: ordinary differential equations
- References: Liu and Bauch 2023; Bury, Bauch, and Anand 2019
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('social_forest_dieback')
entry = describe('social_forest_dieback')
trajectory = create_simulation('social_forest_dieback').prepare().run().trajectory_
Factory
@register_system(
field='physics/climate',
tags=('tipping points', 'vegetation', 'social dynamics', 'replicator', 'delay', 'dde'),
references=('Liu and Bauch 2023', 'Bury, Bauch, and Anand 2019'),
defaults={'until': 200.0, 'dt': 0.5},
)
def social_forest_dieback(
v: float = 0.8,
x: float = 0.5,
*,
Tv: float = 30.0,
Tc: float = 1.0,
alpha: float = 1.0,
delta: float = 1.0,
fmax: float = 5.0,
w: float = 3.0,
tp: float = 10.0,
tf: float = 15.0,
) -> System:
"""Return the forest dieback model coupled with the replicator dynamics of mitigators who perceive its warming.
The forest grows as in `forest_dieback`, with its growth rate scaled by `0.2 + 0.4 x`, so that mitigation by the
share `x` of the population helps it. The share follows the replicator equation
`x' = x (1 - x) (-alpha + f(Tf) + delta (2 x - 1))`, where the cost of warming `f = fmax / (1 + exp(-w (Tf - Tc)))`
reads the warming `Tf = tf / tp (T(t) - T(t - tp))` projected from the last `tp` years, a delay differential
equation. Before the start the cover keeps its initial value. With no warming the shares 0 and 1 are stable and the
share `x* = 1 - (fmax / 2) / (1 + exp(w Tc))` between them is unstable, for the defaults of `alpha`, `delta`, and
`w`. The population converges to all or no mitigators, and the forest to a positive cover or to bare soil, the
green, brown, and barren outcomes, with long swings for some temperatures.
Args:
v:
The initial vegetation cover, and its value before the start.
x:
The initial share of mitigators.
Tv:
The ambient temperature under full forest cover, in degrees Celsius.
Tc:
The critical warming.
alpha:
The cost of mitigation.
delta:
The strength of the social norm.
fmax:
The largest cost of warming.
w:
The steepness of the cost of warming.
tp:
The years of past warming that the projection reads.
tf:
The years ahead of the projection.
Returns:
system:
The system, with the vegetation cover `forest.v` and the share of mitigators `people.x`.
Examples:
>>> from skcomplex.systems.physics.climate import social_forest_dieback
>>> social_forest_dieback(Tc=2.5).elements[1].Tc
2.5
"""
forest = ForestCoverElement('forest', v=v, Tv=Tv)
people = MitigatorsElement('people', x=x, alpha=alpha, delta=delta, fmax=fmax, w=w, Tc=Tc, tp=tp, tf=tf)
growth = RateInteraction(
'growth',
target='forest',
state='v',
sources=['people'],
rate=lambda s: _logistic(s.target, 0.2 + 0.4 * s.sources[0].x),
)
mitigation = RateInteraction('mitigation', target='people', state='x', sources=['forest'], rate=_replicator)
return System(elements=[forest, people], interactions=[growth, mitigation])
Components
Elements:
ForestCoverElement, fromElementMitigatorsElement, fromElement
Interactions:
RateInteraction
class ForestCoverElement(Element):
"""Vegetation cover of a tropical forest, whose growth rate falls on both sides of an optimal temperature.
The cover `v` runs from 0 for bare soil to 1 for full forest. Bare soil is warmer than forest by `a`, so the local
temperature is `T = Tv + a (1 - v)`, and the growth rate `g = g0 (1 - ((T - Topt) / beta)²)` peaks at `Topt`. A
rate interaction gives the logistic growth `g v (1 - v)` and the disturbance `gamma v`, so the element holds the
state and the parameters it reads.
Args:
label:
The label of the forest.
visual:
How the forest is drawn, or `None` for `default_visual_`.
v:
The initial vegetation cover.
Tv:
The ambient temperature, the temperature under full forest cover, in degrees Celsius.
g0:
The largest growth rate.
beta:
The half width of the growth rate against the temperature.
Topt:
The optimal temperature of growth.
a:
The warming of bare soil over forest.
gamma:
The disturbance rate.
Attributes:
default_visual_:
The visual of the forest when it sets none.
Examples:
>>> from skcomplex.systems.physics.climate import ForestCoverElement
>>> ForestCoverElement('forest').Topt
28.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
v: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.8,
Tv: Annotated[float, ParameterSpec()] = 30.0,
g0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
Topt: Annotated[float, ParameterSpec()] = 28.0,
a: Annotated[float, ParameterSpec()] = 5.0,
gamma: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
) -> None:
super().__init__(label, visual=visual)
self.v = v
self.Tv = Tv
self.g0 = g0
self.beta = beta
self.Topt = Topt
self.a = a
self.gamma = gamma
class MitigatorsElement(Element):
"""Share of mitigators in a population that weighs the cost of mitigation against the perceived cost of warming.
The share `x` follows the replicator equation `x' = x (1 - x) (-alpha + f + delta (2 x - 1))`, with the cost of
mitigation `alpha`, the strength `delta` of the social norm, and the cost of warming
`f = fmax / (1 + exp(-w (Tf - Tc)))`. The perceived warming `Tf = tf / tp (T(t) - T(t - tp))` projects the warming
of the last `tp` years over `tf` years. A rate interaction reads the temperature it perceives, so the element holds
the state and the parameters it reads.
Args:
label:
The label of the population.
visual:
How the population is drawn, or `None` for `default_visual_`.
x:
The initial share of mitigators.
alpha:
The cost of mitigation.
delta:
The strength of the social norm.
fmax:
The largest cost of warming.
w:
The steepness of the cost of warming.
Tc:
The critical warming.
tp:
The years of past warming that the projection reads.
tf:
The years ahead of the projection.
Attributes:
default_visual_:
The visual of the population when it sets none.
Examples:
>>> from skcomplex.systems.physics.climate import MitigatorsElement
>>> MitigatorsElement('people').tp
10.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec(domain=RealSet.closed(0, 1))] = 0.5,
alpha: Annotated[float, ParameterSpec()] = 1.0,
delta: Annotated[float, ParameterSpec()] = 1.0,
fmax: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0,
w: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
Tc: Annotated[float, ParameterSpec()] = 1.0,
tp: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
tf: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 15.0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.alpha = alpha
self.delta = delta
self.fmax = fmax
self.w = w
self.Tc = Tc
self.tp = tp
self.tf = tf
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.components.elements import CompartmentElement
>>> from skcomplex.core.components.interactions import RateInteraction
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.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).prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (RateLaw, NoiseLaw, FractionalNoiseLaw)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state', 'divisor'))],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True, elements=(ElementKind.ANY,))] = (),
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
Equations
forest.v: diff(v_forest(t), t) == g0_forest*((a_forest*(v_forest(t) - 1) + Topt_forest - Tv_forest)^2/beta_forest^2 - 1)*(v_forest(t) - 1)*(0.4*x_people(t) + 0.2)*v_forest(t) - gamma_forest*v_forest(t)
people.x: diff(x_people(t), t) == -(delta_people*(2*x_people(t) - 1) - alpha_people + fmax_people/(e^((Tc_people - (a_forest*(v_forest(t - tp_people) - 1) - a_forest*(v_forest(t) - 1))*tf_people/tp_people)*w_people) + 1))*(x_people(t) - 1)*x_people(t)