Kramers escape
Particles that escape over the barrier of U = x^4 / 4 - x^2 / 2 by thermal noise of intensity D.
- Name:
kramers_escape - Field:
physics/soft_matter - Formalism: stochastic differential equations
- References: Kramers 1940; Hänggi, Talkner, and Borkovec 1990; Gros 2015
Factory
@register_model(
field='physics/soft_matter',
tags=('noise', 'sde', 'escape', 'double well', 'rates'),
references=('Kramers 1940', 'Hänggi, Talkner, and Borkovec 1990', 'Gros 2015'),
defaults={'until': 100.0, 'dt': 1.0, 'solver': EulerMaruyamaSolver(dt=0.01)},
)
def kramers_escape(
n_particles: int = 200,
*,
D: float = 0.1,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return particles that escape over the barrier of `U = x^4 / 4 - x^2 / 2` by thermal noise of intensity `D`.
The particles start in the left well, and they hop between the wells at the Kramers rate
`sqrt(U''(-1) |U''(0)|) / (2 pi) exp(-1 / (4 D))`, which is `sqrt(2) / (2 pi) exp(-1 / (4 D))`.
Args:
n_particles:
The number of independent particles.
D:
The noise intensity.
random_state:
The seed or generator of the small initial spread in the left well.
Returns:
system:
The model, with the positions `particles.x`.
Examples:
>>> from skcomplex.models.physics.soft_matter import kramers_escape
>>> kramers_escape(10, D=0.2).elements[0].values['D']
0.2
"""
positions = -1.0 + np.random.default_rng(random_state).normal(0.0, 0.1, n_particles)
group = GroupElement(
'particles',
element=DoubleWellElement,
size=n_particles,
values={'x': positions, 'amplitude': 0.0, 'D': D},
)
return System(elements=[group], interactions=_get_double_well_interactions('particles'))
Components
class DoubleWellElement(Element):
"""Particle at the position `x` in the double well `x^4 / 4 - x^2 / 2`, forced periodically and by noise.
Args:
label:
The label of the particle.
visual:
How the particle is drawn, or `None` for none.
x:
The initial position.
amplitude:
The amplitude of the periodic forcing.
frequency:
The frequency of the forcing.
D:
The noise intensity.
Examples:
>>> from skcomplex.models.physics.soft_matter import DoubleWellElement
>>> DoubleWellElement('particle', D=0.2).D
0.2
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec(role='position')] = -1.0,
amplitude: Annotated[float, ParameterSpec()] = 0.1,
frequency: Annotated[float, ParameterSpec()] = 0.01,
D: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.12,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.amplitude = amplitude
self.frequency = frequency
self.D = D
class GradientInteraction(Interaction):
"""Motion of states down the gradient of a potential `U(ns)`, as a gradient flow or as Newtonian motion.
Every state `x` the interaction moves gets the rate `-dU/dx / d` on its highest time derivative, with `d` the
divisor of its member: a state or parameter of the member, an expression of them such as a moment of inertia
`m l^2`, or one divisor per state, such as the masses of two coordinates of one element. On a state of order 1
this is the gradient flow `x' = -dU/dx / d`, with `d` a friction, and on a state of order 2 it is Newton's law
`x'' = -dU/dx / d`, with `d` an inertia. A vector state gets the gradient of every component, and a tensor state
the matrix gradient, the projection of the derivative onto the matrices of its symmetry.
On a field the potential is the density `f` of a functional `F = integral f`, which may read the field, its
derivatives `partial(x, k)` or `gradient(x)`, and its `laplacian(x)`. A variational interaction moves the field by
the functional derivative `dF/dx = df/dx - sum_k partial(df/dx_k, k) + laplacian(df/d laplacian(x))`, as in the
Allen-Cahn equation, and a conserved interaction gives the rate `laplacian(dF/dx) / d` instead, as in the
Cahn-Hilliard equation, so the total of the field stays constant.
The interaction moves the states of every role that binds one element or group, `target` by default. A subclass
may bind several such roles instead, such as both members of a pair potential, and every one of them moves. The
potential function takes the namespace of the interaction, with the namespaces of the roles and of the `sources`,
the time `t`, and the symbols of the parameters a subclass declares.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states move.
states:
The name of the state that moves, or the names of several states.
potential:
A function of the namespace that returns the potential, or the density of the functional of a field.
sources:
The paths of the elements or groups the potential reads but does not move.
divisor:
The name of the state or parameter of every moving member that divides its rate, a function of the
namespace of the member that returns the divisor, a mapping from the name of every moving state to its
own divisor, or `None`.
variational:
Whether the potential is the density of a functional of a field, moved by its functional derivative.
conserved:
Whether the rate is the Laplacian of the derivative of the potential, which conserves the total of a field.
Raises:
ValidationError:
A state or the divisor is not a name, or the potential is not callable.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import GradientInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> well = GradientInteraction('well', target='walker', states='n', potential=lambda s: s.target.n**2 / 2)
>>> system = System(elements=[CompartmentElement('walker', n=1.0)], interactions=[well])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=1.0).trajectory_
>>> round(float(trajectory.values('walker.n')[-1]), 4)
0.3679
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
states: str | Sequence[str] = 'x',
potential: Callable[[SimpleNamespace], object] | None = None,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
divisor: (
str | Callable[[SimpleNamespace], object] | Mapping[str, str | Callable[[SimpleNamespace], object]] | None
) = None,
variational: bool = False,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.states = states
self.potential = potential
self.sources = sources
self.divisor = divisor
self.variational = variational
self.conserved = conserved
def _get_names(self) -> tuple[str, ...]:
"""Return the names of the states the interaction moves."""
return (self.states,) if isinstance(self.states, str) else tuple(self.states)
def _check(self) -> None:
"""Raise an error for a state or a divisor that is not a name, or a potential that is not callable.
Raises:
ValidationError:
The type declares states, a state is not a name, a divisor is neither a name nor a function, or the
potential is not callable.
"""
super()._check()
names = self._get_names()
if not names or not all(isinstance(name, str) and name.isidentifier() for name in names):
error_msg = f'{type(self).__name__} {self.label!r} moves the states {self.states!r}. Give names of states.'
raise ValidationError(error_msg)
divisors = dict(self.divisor) if isinstance(self.divisor, Mapping) else dict.fromkeys(names, self.divisor)
valid = all(
item is None or callable(item) or (isinstance(item, str) and item.isidentifier())
for item in divisors.values()
)
if not valid or not set(divisors) <= set(names):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, a function '
'of the namespace of the member, or a mapping from the moving states to them.'
)
raise ValidationError(error_msg)
if self.potential is not None and not callable(self.potential):
error_msg = f'The potential of {self.label!r} is {self.potential!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_potential(self, namespace: SimpleNamespace) -> object:
"""Return the potential, the hook a subclass overrides to derive it from its parameters."""
if self.potential is None:
error_msg = f'{type(self).__name__} {self.label!r} has no potential. Give a function of the namespace.'
raise ValidationError(error_msg)
return self.potential(namespace)
def _get_divisor(self, member: SimpleNamespace, name: str) -> object:
"""Return the divisor of a moving state of a member, or `None` for none."""
divisor = self.divisor.get(name) if isinstance(self.divisor, Mapping) else self.divisor
if divisor is None:
return None
return getattr(member, divisor) if isinstance(divisor, str) else divisor(member)
def _get_rate_state(self, member: SimpleNamespace, name: str) -> object:
"""Return the state whose highest derivative gets the rate of a moving state, the moving state itself."""
return getattr(member, name)
def _walk(self, expression: Expression) -> Iterator[Expression]:
"""Yield an expression and every subexpression of it."""
yield expression
for operand in expression.operands():
yield from self._walk(operand)
def _get_flux_divergence(
self,
fluxes: dict[int, Expression],
derivatives: dict[int, Expression],
field: Expression,
) -> Expression:
"""Return `sum_k partial(flux_k, k)`, or `c laplacian(x)` when every flux is `c x_k` along every axis read.
The Laplacian keeps the compact stencil of the field, which nested first derivatives would widen.
"""
coefficients = [SR(fluxes[axis] / derivatives[axis]) for axis in sorted(fluxes)]
coefficient = coefficients[0]
uniform = sorted(fluxes) == list(range(len(fluxes))) and all(
item.is_trivially_equal(coefficient) for item in coefficients
)
constant = not coefficient.has(field) and not any(coefficient.has(item) for item in derivatives.values())
if uniform and constant:
return coefficient * laplacian(field)
return sum((partial(flux, SR(axis)) for axis, flux in fluxes.items()), SR(0))
def _get_variation(self, density: Expression, field: Expression) -> Expression:
"""Return the functional derivative of the integral of a density with respect to a field."""
derivatives: dict[int, Expression] = {}
curvature = None
for item in self._walk(density):
kind = item.operator()
if kind == partial and item.operands()[0].is_trivially_equal(field):
derivatives[int(item.operands()[1])] = item
elif kind == laplacian and item.operands()[0].is_trivially_equal(field):
curvature = item
slots = {item: SR.var(f'_slot_{axis}') for axis, item in derivatives.items()}
if curvature is not None:
slots[curvature] = SR.var('_slot_laplacian')
flat = density.subs(slots)
back = {symbol: item for item, symbol in slots.items()}
variation = flat.diff(field).subs(back)
fluxes = {axis: flat.diff(slots[item]).subs(back) for axis, item in derivatives.items()}
if fluxes:
variation -= self._get_flux_divergence(fluxes, derivatives, field)
if curvature is not None:
variation += self._get_laplacian(flat.diff(slots[curvature]).subs(back), field)
return variation
def _get_laplacian(self, expression: Expression, field: Expression) -> Expression:
"""Return the Laplacian of an expression of a field, term by term, with the factors free of the field outside.
The linear terms of the field, such as `laplacian(x)`, stay apart, so a spectral discretization finds them.
"""
expanded = SR(expression).expand()
terms = expanded.operands() if expanded.operator() is add_vararg else [expanded]
total = SR(0)
for term in terms:
factors = term.operands() if term.operator() is mul_vararg else [term]
constant = prod((factor for factor in factors if not factor.has(field)), start=SR(1))
varying = prod((factor for factor in factors if factor.has(field)), start=SR(1))
if varying.has(field):
total += constant * laplacian(varying)
return total
def _get_matrix_derivative(self, potential: Expression, state: Matrix) -> Matrix:
"""Return the matrix gradient of a potential at a tensor state, written in the matrices of its symmetry.
The derivative with respect to the independent components is the product of the matrices of the components, the
metric of the state, and the coefficients of the gradient, so the metric gives the coefficients.
"""
if self.conserved:
error_msg = f'{self.label!r} conserves the tensor state {state}. Give a scalar state.'
raise ValidationError(error_msg)
components = sorted(state.variables(), key=str)
units = [state.derivative(component) for component in components]
metric = matrix(SR, [[(first.transpose() * second).trace() for second in units] for first in units])
density = distribute_operators(potential)
derivatives = [
self._get_variation(density, component) if self.variational else density.diff(component)
for component in components
]
coefficients = metric.inverse() * vector(SR, derivatives)
total = matrix(SR, state.nrows(), state.ncols())
for coefficient, unit in zip(coefficients, units, strict=True):
total += coefficient * unit
return total
def _get_derivative(self, potential: Expression, state: Target | Matrix) -> Target | Matrix:
"""Return the derivative of the potential with respect to a state, every component of a vector state."""
if isinstance(state, Matrix):
return self._get_matrix_derivative(potential, state)
if isinstance(state, FreeModuleElement):
if self.variational or self.conserved:
error_msg = f'{self.label!r} moves the vector state {state} of a field. Give a scalar state.'
raise ValidationError(error_msg)
return vector(SR, [potential.diff(component) for component in state])
if self.variational:
return self._get_variation(potential, state)
return potential.diff(state)
def _get_members(self, namespace: SimpleNamespace) -> list[SimpleNamespace]:
"""Return the namespaces of the members whose states move, those of every role that binds one path."""
return [
getattr(namespace, name)
for name, spec in self._get_specs().items()
if isinstance(spec, RoleSpec) and not spec.variadic and not spec.noise
]
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate of every moving state, minus the derivative of the potential over the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the roles.
Returns:
laws:
The rate laws of the moving states.
"""
potential = SR(self._get_potential(namespace))
laws: list[BasePart] = []
for member in self._get_members(namespace):
for name in self._get_names():
state = getattr(member, name)
derivative = self._get_derivative(potential, state)
rate = self._get_laplacian(derivative, state) if self.conserved else -derivative
divisor = self._get_divisor(member, name)
if divisor is not None:
rate = rate / divisor
laws.append(RateLaw(self._get_rate_state(member, name), rate))
return laws
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
Equations
particles.x[i]: diff(x_particles(t), t) == -x_particles(t)^3 + amplitude_particles*cos(6.283185307179586*frequency_particles*t) + x_particles(t)