Lévy walk
The Lévy walk, flights at constant speed in random directions for times with the tail tau^(-1 - sigma).
- Name:
levy_walk - Field:
physics/soft_matter - Formalism: ordinary differential equations
- References: Argun et al. 2021; Zaburdaev, Denisov, and Klafter 2015; Shlesinger, West, and Klafter 1987
Factory
@register_model(
field='physics/soft_matter',
tags=('anomalous diffusion', 'renewal process', 'superdiffusion', 'random walk'),
references=('Argun et al. 2021', 'Zaburdaev, Denisov, and Klafter 2015', 'Shlesinger, West, and Klafter 1987'),
defaults={'until': 500.0, 'dt': 2.5, 'solver': RK4Solver(dt=0.5)},
)
def levy_walk(
n_walkers: int = 500,
*,
sigma: float = 1.5,
tau0: float = 1.0,
speed: float = 1.0,
n_dim: int = 2,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the Lévy walk, flights at constant speed in random directions for times with the tail `tau^(-1 - sigma)`.
For `1 < sigma < 2` the flight times have a mean but no variance, so long flights dominate and the ensemble mean
squared displacement grows as `t^(3 - sigma)`, faster than Brownian motion and slower than ballistic motion.
Args:
n_walkers:
The number of independent walkers.
sigma:
The tail exponent of the flight times.
tau0:
The shortest flight time.
speed:
The speed of every flight.
n_dim:
The number of dimensions.
random_state:
The seed or generator of the directions of the first flights.
Returns:
system:
The model, with the positions `walkers.x`, the velocities `walkers.v`, and the turn count `turns.renewals`.
Examples:
>>> from skcomplex.models.physics.soft_matter import levy_walk
>>> levy_walk(10, sigma=1.2).elements[0].values['sigma']
1.2
"""
directions = np.random.default_rng(random_state).standard_normal((n_walkers, n_dim))
velocities = speed * directions / np.linalg.norm(directions, axis=1, keepdims=True)
walkers = GroupElement(
'walkers',
element=LevyWalkerElement,
size=n_walkers,
values={'x': np.zeros((n_walkers, n_dim)), 'v': velocities, 'sigma': sigma, 'tau0': tau0, 'speed': speed},
)
flight = RateInteraction('flight', target='walkers', state='x', rate=lambda s: s.target.v)
turns = JumpInteraction(
'turns',
target='walkers',
waiting=lambda s: draw('flight', 'pareto', b=s.target.sigma, scale=s.target.tau0),
reset={'v': lambda s: [s.target.speed * item for item in _direction('direction', len(list(s.target.v)))]},
)
return System(elements=[walkers], interactions=[flight, turns], space=ContinuousSpace(n_dim=n_dim))
Components
class LevyWalkerElement(Element):
"""Walker that flies at constant speed in a random direction for a power-law flight time, then turns.
The flight is a rate interaction of the position, and the turns a jump interaction that resets the velocity after
every flight time.
Args:
label:
The label of the walker or the group.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The initial position.
v:
The initial velocity.
sigma:
The tail exponent of the flight times, between 1 and 2 for superdiffusion.
tau0:
The shortest flight time.
speed:
The speed of every flight.
Examples:
>>> from skcomplex.models.physics.soft_matter import LevyWalkerElement
>>> LevyWalkerElement('walker', speed=2.0).speed
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = POINTS,
*,
x: Annotated[npt.ArrayLike, StateSpec(vector=True, role='position')] = 0.0,
v: Annotated[npt.ArrayLike, StateSpec(vector=True, role='velocity')] = 0.0,
sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.5,
tau0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
speed: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.v = v
self.sigma = sigma
self.tau0 = tau0
self.speed = speed
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
class JumpInteraction(Interaction):
"""Random event that changes states of the target by signed amounts, at a rate, by a probability, or after a wait.
With a `rate`, the event fires at the times of a Poisson process of that propensity in continuous time. With a
`probability`, it fires in a step of discrete time with that probability. With a `waiting` time, usually a draw from
a distribution with nonnegative support, it fires after a time drawn anew at every firing, as in renewal processes,
and every member of a bound group fires on its own clock. When the event fires it adds the `changes` to the states
and, with a waiting time, sets the states of `reset` to their new values. Every function reads `namespace.target`,
the states and parameters of the target, and the time `namespace.t`.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states the event changes.
changes:
A mapping from the name of a state of the target to the signed amount the event adds to it.
rate:
The function `rate(namespace)` of the propensity of the event, or `None`.
probability:
The function `probability(namespace)` of the probability that the event fires in a step, or `None`.
waiting:
The function `waiting(namespace)` of the time until the next firing, or `None`.
reset:
With a waiting time, a mapping from the name of a state of the target to the function `f(namespace)` of its
new value, one value per component of a vector state, or `None`.
Raises:
ValidationError:
The event has not exactly one clock, changes nothing, or resets states without a waiting time, or a
function is not callable.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import JumpInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> arrivals = JumpInteraction('arrivals', target='queue', changes={'n': 1}, rate=lambda s: 2.0)
>>> system = System(elements=[CompartmentElement('queue')], interactions=[arrivals])
>>> simulation = Simulation(system, dt=1.0, random_state=0).check().prepare()
>>> counts = simulation.run(until=5).trajectory_.values('queue.n')
>>> bool(all(count == round(count) for count in counts))
True
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
changes: Mapping[str, float] | None = None,
rate: Callable[[SimpleNamespace], object] | None = None,
probability: Callable[[SimpleNamespace], object] | None = None,
waiting: Callable[[SimpleNamespace], object] | None = None,
reset: Mapping[str, Callable[[SimpleNamespace], object]] | None = None,
) -> None:
super().__init__(label)
self.target: str | None = target
self.changes = changes
self.rate = rate
self.probability = probability
self.waiting = waiting
self.reset = reset
def _get_clocks(self) -> list[str]:
"""Return the names of the clocks the event is given."""
return [clock for clock in CLOCKS if getattr(self, clock) is not None]
def _check_jump(self) -> None:
"""Raise an error for an event without exactly one clock, without changes, or with a reset of no waiting time.
Raises:
ValidationError:
The event has not exactly one clock, changes nothing, or resets states without a waiting time, or a
function is not callable.
"""
clocks = self._get_clocks()
if len(clocks) != 1:
error_msg = f'The jump {self.label!r} has the clocks {clocks}. Give exactly one of {", ".join(CLOCKS)}.'
raise ValidationError(error_msg)
if not callable(getattr(self, clocks[0])):
error_msg = f'The {clocks[0]} of the jump {self.label!r} must be a function of the namespace.'
raise ValidationError(error_msg)
if self.reset is not None and self.waiting is None:
error_msg = f'The jump {self.label!r} resets states, which only a waiting time does. Give `waiting`.'
raise ValidationError(error_msg)
if not self.changes and not self.reset:
error_msg = f'The jump {self.label!r} changes no state. Give `changes`, or `reset` with a waiting time.'
raise ValidationError(error_msg)
if not all(callable(function) for function in (self.reset or {}).values()):
error_msg = f'The reset of the jump {self.label!r} maps states to {self.reset!r}. Give functions.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for an event without exactly one clock, without changes, or with a reset of no waiting time.
Raises:
ValidationError:
The event has not exactly one clock, changes nothing, or resets states without a waiting time, or a
function is not callable.
"""
super()._check()
self._check_jump()
def _get_changes(self, namespace: SimpleNamespace) -> dict[Expression, float]:
"""Return the amount the event adds to every state it changes."""
return {
getattr(self._get_member(namespace, 'target'), name): amount
for name, amount in (self.changes or {}).items()
}
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the propensity of the event."""
return SR(self.rate(namespace)) if self.rate is not None else SR(0)
def _get_steps(self, namespace: SimpleNamespace, changes: dict[Expression, float]) -> list[BaseLaw]:
"""Return the updates of the states by their amounts in a step, when one draw falls below the probability."""
probability = SR(self.probability(namespace)) if self.probability is not None else SR(0)
name = re.sub(r'\W', '_', f'jump_{self.label}')
fired = where(uniform(name) < probability, 1, 0)
return [UpdateLaw(state, state + amount * fired) for state, amount in changes.items()]
def _get_reset(self, state: object, value: object) -> dict[Expression, Expression]:
"""Return the new value of every component of a state, one per component of a vector state.
Args:
state:
The symbol of a scalar state or the vector of a vector state.
value:
The new value, one expression for every component or a sequence of one per component.
Returns:
reset:
The new value of every component of the state.
Raises:
ValidationError:
The new value of a vector state has not one component per component of the state.
"""
if not isinstance(state, FreeModuleElement):
return {SR(state): SR(value)}
components = list(value) if isinstance(value, FreeModuleElement | list | tuple) else [value] * len(state)
if len(components) != len(state):
error_msg = f'The jump {self.label!r} resets {len(state)} components with {len(components)} values.'
raise ValidationError(error_msg)
return {SR(item): SR(component) for item, component in zip(state, components, strict=True)}
def _get_renewal(self, namespace: SimpleNamespace, changes: dict[Expression, float]) -> JumpLaw:
"""Return the renewal of the event, which adds the changes and sets the reset states after every wait."""
reset: dict[Expression, object] = {state: state + amount for state, amount in changes.items()}
for name, function in (self.reset or {}).items():
reset.update(self._get_reset(getattr(self._get_member(namespace, 'target'), name), function(namespace)))
waiting = SR(self.waiting(namespace)) if self.waiting is not None else SR(0)
return JumpLaw(waiting=waiting, reset=reset, label=self.label)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the jump law of the event at a rate or after a wait, or the updates of a step with a probability.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The laws of the event.
"""
changes = self._get_changes(namespace)
if self.waiting is not None:
return [self._get_renewal(namespace, changes)]
if self.probability is not None:
return self._get_steps(namespace, changes)
return [JumpLaw(self._get_rate(namespace), changes)]
Equations
walkers.x[i, 0]: diff(x_walkers_0(t), t) == v_walkers_0(t)
walkers.x[i, 1]: diff(x_walkers_1(t), t) == v_walkers_1(t)
walkers.v[i, 0]: diff(v_walkers_0(t), t) == 0
walkers.v[i, 1]: diff(v_walkers_1(t), t) == 0