Solvers

Advance compiled problems through time, with a default solver per form.

`AgentView(agent, stepper)`

Agent as rules see it: its states as attributes, its group's parameters, its neighbours, and its location.

The stepper gives every agent one view of the class of its group, whose states and parameters are properties, so two views of an agent are the same object. Reading a state reads the agent, and setting it sets the agent, at the end of the phase in a simultaneous phase. Parameters are read only, and an unknown name raises Python's own AttributeError.

Parameters:

Name	Type	Description	Default
`agent`	`Agent`	The Mesa agent.	required
`stepper`	`MesaSolver`	The stepper that runs the model.	required

Source code in src/skcomplex/engine/solvers/_mesa.py

def __init__(self, agent: mesa.Agent, stepper: 'MesaSolver') -> None:
    self._agent = agent
    self._stepper = stepper

`alive` `property`

Whether the agent has not been removed.

`group` `property`

Path of the agent's group.

`group_label` `property`

Label of the agent's group, the last part of its path.

`id` `property`

Unique id of the agent.

`location` `property`

Grid coordinates, network node, or position of the agent.

`location_key` `property`

Location of the agent as a key: a tuple of integers for a cell, a tuple of floats for a position, or a node.

The key equals the coordinates neighbor_cells gives for the agent's cell or node, so rules compare the cells around an agent with the places of other agents.

`agents_here(group=None)`

Return the other agents in the agent's cell or node, of one group or of every group.

Parameters:

Name	Type	Description	Default
`group`	`str \| None`	The path of a group, or `None` for every group.	`None`

Returns:

Name	Type	Description
`agents`	`list[AgentView]`	The other agents at the same place.

Source code in src/skcomplex/engine/solvers/_mesa.py

def agents_here(self, group: str | None = None) -> list['AgentView']:
    """Return the other agents in the agent's cell or node, of one group or of every group.

    Args:
        group:
            The path of a group, or `None` for every group.

    Returns:
        agents:
            The other agents at the same place.
    """
    stepper = self._stepper
    path = stepper._resolve_group(group, self.group) if group is not None else None
    return [
        MesaSolver._view_of(agent)
        for agent in self._agent.cell.agents
        if agent is not self._agent and (path is None or stepper._spec_of(agent).path == path)
    ]

`create(group=None, location=None, **values)`

Create an agent of a group, by default the agent's own group at the agent's location.

Parameters:

Name	Type	Description	Default
`group`	`str \| None`	The path of the group, or `None` for the agent's group.	`None`
`location`	`object`	The location of the new agent, or `None` for the agent's location.	`None`
`**values`	`object`	The initial values of states, with the element type's defaults for the others.	`{}`

Returns:

Name	Type	Description
`agent`	`AgentView`	The new agent.

Raises:

Type	Description
`SimulationError`	The group is not variable, so its size cannot change.

Source code in src/skcomplex/engine/solvers/_mesa.py

def create(self, group: str | None = None, location: object = None, **values: object) -> 'AgentView':
    """Create an agent of a group, by default the agent's own group at the agent's location.

    Args:
        group:
            The path of the group, or `None` for the agent's group.

        location:
            The location of the new agent, or `None` for the agent's location.

        **values:
            The initial values of states, with the element type's defaults for the others.

    Returns:
        agent:
            The new agent.

    Raises:
        SimulationError:
            The group is not variable, so its size cannot change.
    """
    stepper = self._stepper
    target = stepper._resolve_group(group, self.group) if group is not None else self.group
    place = location if location is not None else self.location
    return MesaSolver._view_of(stepper._create_agent(target, place, values))

`displacement_to(other)`

Return the vector from the agent to another agent or a location, the shortest one on a torus.

Parameters:

Name	Type	Description	Default
`other`	`AgentView \| object`	An agent, or grid coordinates or a position.	required

Returns:

Name	Type	Description
`displacement`	`NDArray[float64]`	One component per axis of the space.

Raises:

Type	Description
`SpaceError`	The agents live on a network or in no space.

Source code in src/skcomplex/engine/solvers/_mesa.py

def displacement_to(self, other: 'AgentView | object') -> npt.NDArray[np.float64]:
    """Return the vector from the agent to another agent or a location, the shortest one on a torus.

    Args:
        other:
            An agent, or grid coordinates or a position.

    Returns:
        displacement:
            One component per axis of the space.

    Raises:
        SpaceError:
            The agents live on a network or in no space.
    """
    model = ModelView(self._stepper)
    dimensions = model.dimensions
    if dimensions is None:
        error_msg = (
            'Displacements need a grid or a continuous space, and the agents live on a network or in no space. '
            'Use `distance_to` for the number of hops on a network.'
        )
        raise SpaceError(error_msg)
    target = other.location if isinstance(other, AgentView) else other
    delta = np.asarray(target, dtype=float) - np.asarray(self.location, dtype=float)
    if not model.torus:
        return delta
    size = np.asarray(dimensions, dtype=float)
    return delta - size * np.round(delta / size)

`distance_to(other)`

Return the distance to another agent or a location, the number of hops on a network.

Parameters:

Name	Type	Description	Default
`other`	`AgentView \| object`	An agent, or a network node, grid coordinates, or a position.	required

Returns:

Name	Type	Description
`distance`	`float`	The Euclidean length of the shortest displacement, or the length of the shortest path on a network.

Raises:

Type	Description
`SpaceError`	The agents live in no space.

Source code in src/skcomplex/engine/solvers/_mesa.py

def distance_to(self, other: 'AgentView | object') -> float:
    """Return the distance to another agent or a location, the number of hops on a network.

    Args:
        other:
            An agent, or a network node, grid coordinates, or a position.

    Returns:
        distance:
            The Euclidean length of the shortest displacement, or the length of the shortest path on a network.

    Raises:
        SpaceError:
            The agents live in no space.
    """
    space = self._stepper.problem_.agent_space
    if isinstance(space, NetworkSpace):
        target = other.location if isinstance(other, AgentView) else other
        return float(nx.shortest_path_length(space.graph, self.location, target))
    return float(np.linalg.norm(self.displacement_to(other)))

`lags(state, window=None)`

Return the past values of a state of an agent, the values at the ends of the last steps, the last one last.

An agent remembers a state declared with lags, StateSpec(lags=m): the agent stepper keeps its last m values in the companion states <state>_lag1 to <state>_lag<m>, which rules may also read one by one. Before the first steps the lags hold their initial values, the initial of the spec or else the initial value of the state.

Parameters:

Name	Type	Description	Default
`state`	`str`	The name of a state with lags.	required
`window`	`int \| None`	The number of last values, or `None` for every kept value.	`None`

Returns:

Name	Type	Description
`values`	`NDArray[float64]`	The values, the oldest first.

Raises:

Type	Description
`KeyError`	The state of the agent keeps no lags.

Source code in src/skcomplex/engine/solvers/_mesa.py

def lags(self, state: str, window: int | None = None) -> npt.NDArray[np.float64]:
    """Return the past values of a state of an agent, the values at the ends of the last steps, the last one last.

    An agent remembers a state declared with lags, `StateSpec(lags=m)`: the agent stepper keeps its last `m` values
    in the companion states `<state>_lag1` to `<state>_lag<m>`, which rules may also read one by one. Before the
    first steps the lags hold their initial values, the `initial` of the spec or else the initial value of the
    state.

    Args:
        state:
            The name of a state with lags.

        window:
            The number of last values, or `None` for every kept value.

    Returns:
        values:
            The values, the oldest first.

    Raises:
        KeyError:
            The state of the agent keeps no lags.
    """
    spec = self._stepper._spec_of(self._agent)
    declaration = declarations_of(MesaSolver._agent_type(spec), VariableKind.STATE).get(state)
    if declaration is None or not declaration.lags:
        error_msg = (
            f'The agents of {spec.path} do not remember {state!r}. Declare the state with lags, such as '
            '`StateSpec(lags=3)`, to keep its last values.'
        )
        raise KeyError(error_msg)
    lags = declaration.lags if window is None else min(window, declaration.lags)
    return np.asarray([getattr(self, f'{state}_lag{index}') for index in range(lags, 0, -1)], dtype=float)

`move_to(location)`

Move the agent to a cell, a node, or a position.

Parameters:

Name	Type	Description	Default
`location`	`object`	The cell coordinates, network node, or position.	required

Source code in src/skcomplex/engine/solvers/_mesa.py

def move_to(self, location: object) -> None:
    """Move the agent to a cell, a node, or a position.

    Args:
        location:
            The cell coordinates, network node, or position.
    """
    if isinstance(self._stepper.problem_.agent_space, ContinuousSpace):
        self._agent.position = np.asarray(location, dtype=float)
        return
    self._agent.cell = self._stepper._cells()[MesaSolver._cell_key(location)]

`neighbor_cells()`

Return the coordinates of the cells or nodes around the agent's cell or node.

Returns:

Name	Type	Description
`locations`	`list[object]`	The neighbouring cell coordinates or network nodes.

Source code in src/skcomplex/engine/solvers/_mesa.py

def neighbor_cells(self) -> list[object]:
    """Return the coordinates of the cells or nodes around the agent's cell or node.

    Returns:
        locations:
            The neighbouring cell coordinates or network nodes.
    """
    return [cell.coordinate for cell in self._agent.cell.neighborhood]

`neighbors(radius=None, layer=None)`

Return the agents in the cells or nodes around the agent, within a radius, or around it in a layer.

Parameters:

Name	Type	Description	Default
`radius`	`float \| None`	The distance within which agents are neighbours in a continuous space.	`None`
`layer`	`str \| None`	The layer of a network space whose edges join the agent to its neighbours, or `None` for the graph of the network.	`None`

Returns:

Name	Type	Description
`agents`	`list[AgentView]`	The neighbouring agents.

Source code in src/skcomplex/engine/solvers/_mesa.py

def neighbors(self, radius: float | None = None, layer: str | None = None) -> list['AgentView']:
    """Return the agents in the cells or nodes around the agent, within a radius, or around it in a layer.

    Args:
        radius:
            The distance within which agents are neighbours in a continuous space.

        layer:
            The layer of a network space whose edges join the agent to its neighbours, or `None` for the graph of
            the network.

    Returns:
        agents:
            The neighbouring agents.
    """
    if layer is not None:
        graph = getattr(self._stepper.problem_.agent_space, 'graphs_', {})[layer]
        cells = self._stepper._cells()
        return [
            MesaSolver._view_of(agent)
            for node in graph.neighbors(self.location_key)
            for agent in cells[MesaSolver._cell_key(node)].agents
        ]
    if radius is not None:
        agents, _ = self._agent.get_neighbors_in_radius(radius)
        return [MesaSolver._view_of(agent) for agent in agents if agent is not self._agent]
    return [MesaSolver._view_of(agent) for agent in self._agent.cell.neighborhood.agents]

`referenced(state)`

Return the member of another group that a reference state of the agent names by its id.

A state declared with StateSpec(reference=group) holds the id of a member of the group, or a negative value for none, such as the bank of a customer or the employer of a worker.

Parameters:

Name	Type	Description	Default
`state`	`str`	The name of a reference state.	required

Returns:

Name	Type	Description
`agent`	`AgentView \| None`	The member the state names, or `None` when it names none or the member was removed.

Raises:

Type	Description
`KeyError`	The state of the agent names no member of a group.

Source code in src/skcomplex/engine/solvers/_mesa.py

def referenced(self, state: str) -> 'AgentView | None':
    """Return the member of another group that a reference state of the agent names by its id.

    A state declared with `StateSpec(reference=group)` holds the id of a member of the group, or a negative value
    for none, such as the bank of a customer or the employer of a worker.

    Args:
        state:
            The name of a reference state.

    Returns:
        agent:
            The member the state names, or `None` when it names none or the member was removed.

    Raises:
        KeyError:
            The state of the agent names no member of a group.
    """
    stepper = self._stepper
    path = stepper._reference_of(stepper._spec_of(self._agent), state)
    value = float(getattr(self._agent, state))
    if value < 0 or path is None:
        return None
    agent = stepper.members_.get(int(value))
    if agent is None or getattr(agent, 'removed_', False) or stepper._spec_of(agent).path != path:
        return None
    return MesaSolver._view_of(agent)

`remove()`

Remove the agent from the model, after checking that it holds none of a conserved quantity.

Raises:

Type	Description
`ConservationError`	A state tagged with a quantity is not zero, so removing the agent would destroy it. Settle the holdings with transfers first.
`SimulationError`	The agent's group is not variable, so its size cannot change.

Source code in src/skcomplex/engine/solvers/_mesa.py

def remove(self) -> None:
    """Remove the agent from the model, after checking that it holds none of a conserved quantity.

    Raises:
        ConservationError:
            A state tagged with a quantity is not zero, so removing the agent would destroy it. Settle the
            holdings with transfers first.

        SimulationError:
            The agent's group is not variable, so its size cannot change.
    """
    spec = self._stepper._spec_of(self._agent)
    MesaSolver._check_variable(spec, 'remove')
    held = {
        name: value
        for name in self._stepper._quantities[spec.path]
        if (value := float(getattr(self._agent, name))) != 0.0
    }
    if self.alive and held:
        error_msg = (
            f'Agent {self.id} of {spec.path} is removed while holding {held}. Settle its holdings with transfers '
            'to other agents before removing it, so conserved quantities are not destroyed.'
        )
        raise ConservationError(error_msg)
    if self.alive:
        self._stepper._recount(self._agent, -1)
        self._agent.removed_ = True
        self._agent.remove()

`BaseSolver(*, random_state=None, until=math.inf)`

Bases: BaseConfig[P, StepSpec]

Base of the solvers, which advance a problem of one form through time.

A solver stores its settings in its constructor, among them the random state of stochastic solvers and the end time of the run it prepares for. prepare(problem) prepares it for a problem and run(step) advances the state of a StepSpec from its time to its end time and keeps the state vector there in state_. The last call run(step, final=True) releases what preparing acquired, such as an external process.

Subclasses implement the hooks _get_form, the form they solve, _get_support, what they solve, and _advance, and solvers that keep events, agents, or other state outside the state vector also implement _collect. Solvers that also solve other forms list them in _get_forms. Solvers that treat diffusion implicitly answer _is_implicit, so they do not warn about the stability limit of explicit steps. Solvers that fire renewal laws answer _is_renewing, solvers that integrate fractional noise _is_fractional, and solvers that read external inputs _reads_inputs. A solver is generic over the problem it advances, so one written for the compiler's problem reads its fields.

Parameters:

Name	Type	Description	Default
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers. A seed starts a stream spawned from `np.random.SeedSequence(random_state)`, so the dynamics never repeat the numbers that `np.random.default_rng(random_state)` gives a factory or a user for initial states; a generator is used as it is.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped, which solvers that plan or draw a whole run read.	`inf`

Attributes:

Name	Type	Description
`problem_`	`P`	The problem the solver is prepared to.
`rng_`	`Generator`	The random number generator of stochastic solvers.
`initial_state_`	`Array`	The state vector of the first frame, which solvers that derive states fill in.
`events_`	`list[BaseRecordSpec]`	The events of the last run, in the order they happened.
`stopped_at_`	`float \| None`	The time at which a terminal event ended the run, or `None` while it goes on.
`snapshot_`	`object \| None`	The state the solver keeps outside the state vector, such as agents, after the last run, or `None`.
`agents_`	`Mapping[str, list[object]]`	The views of the agents alive after the last run, as rules receive them, by path of their agent group.
`state_`	`Array`	The state vector at the end time of the last step.

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`CaputoSolver(*, random_state=None, until=math.inf)`

Bases: Solver

Base class of steppers of Caputo fractional differential equations on a fixed grid from the start of a run.

A state of order alpha follows D^alpha x = f(t, x), and a state of order 1 follows dx/dt = f(t, x). The stepper keeps the history of the run on the grid t0 + k dt. A change of a fractional state between frames starts a new history from the changed state, and a change of a state of order 1 carries on.

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`CoSimulationSolver(solvers=None, sync_dt=None, scheme='gauss-seidel', split=False, *, random_state=None, until=math.inf)`

Bases: Solver

Scheduler that advances every partition of a co-simulation to each synchronisation point with its own solver.

Between synchronisation points a partition holds the states of the other partitions at their last values, so a split adds an error of the order of sync_dt. A law of one partition that reads a reduction of the agents of another, such as count(s.farmers), reads the value the scheduler computes at every synchronisation point and holds until the next. Gauss-Seidel advances the partitions in system order, each reading the values the partitions before it already advanced. Jacobi advances every partition from the same values.

Parameters:

Name	Type	Description	Default
`solvers`	`Mapping[str, BaseSolver] \| None`	The solver of a partition by its name, with the default solver of its form for the others.	`None`
`sync_dt`	`float \| None`	The time between synchronisation points, which must divide the frame duration, or `None` for one synchronisation per frame.	`None`
`scheme`	`str`	`'gauss-seidel'` or `'jacobi'`.	`'gauss-seidel'`
`split`	`bool`	Whether to split a problem of one form into one partition per subsystem.	`False`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`steppers_`		The stepper of every partition by partition name.

Examples:

>>> from skcomplex.core.system import System
>>> from skcomplex.models.biology.common.populations import ExponentialGrowthInteraction, PopulationElement
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import CoSimulationSolver
>>> system = System(
...     systems=[
...         System(
...             'a',
...             elements=[PopulationElement('x', n=1)],
...             interactions=[ExponentialGrowthInteraction('g', target='x')],
...         ),
...         System(
...             'b',
...             elements=[PopulationElement('y', n=1)],
...             interactions=[ExponentialGrowthInteraction('g', target='y')],
...         ),
...     ],
... )
>>> simulation = Simulation(system, solver=CoSimulationSolver(split=True), dt=0.5)
>>> simulation = simulation.check().prepare().run(until=1)
>>> sorted(simulation.stepper_.steppers_)
['a', 'b']

Source code in src/skcomplex/engine/solvers/_cosimulation.py

def __init__(
    self,
    solvers: Mapping[str, BaseSolver] | None = None,
    sync_dt: float | None = None,
    scheme: str = 'gauss-seidel',
    split: bool = False,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.solvers = solvers
    self.sync_dt = sync_dt
    self.scheme = scheme
    self.split = split
    self.random_state = random_state
    self.until = until

`CollocationSolver(objective=None, *, terminal=0.0, sense='min', nodes=101, scheme='hermite-simpson', method='SLSQP', tolerance=1e-10, max_iterations=1000, random_state=None, until=math.inf)`

Bases: OptimalControlSolver

Optimal control by direct collocation, a transcription of the problem solved by nonlinear programming.

The states and the controls at the points of a grid are the variables of a nonlinear programme. The trapezoidal scheme requires the states of neighbouring nodes to differ by the trapezoidal integral of their rates, and the Hermite-Simpson scheme adds a midpoint to every interval, with a cubic interpolation of the states and Simpson's integral. The cost is the quadrature of the running cost over the points plus the terminal cost. The controls lie within their bounds and the states of nonnegative domains stay nonnegative. SciPy's minimize solves the programme with the gradient of the cost and the Jacobian of the constraints derived symbolically. An objective to maximise is minimised with its sign changed.

Parameters:

Name	Type	Description	Default
`objective`	`Cost \| None`	The running cost, a function of the namespace `s` of the states, parameters, and time of the system by dotted name, such as `lambda s: s.plant.x**2`, or a number.	`None`
`terminal`	`Cost`	The terminal cost at the horizon, a function of the same namespace or a number.	`0.0`
`sense`	`str`	`'min'` to minimise the objective or `'max'` to maximise it.	`'min'`
`nodes`	`int`	The number of times of the grid, from the start to the horizon.	`101`
`scheme`	`str`	`'trapezoidal'` or `'hermite-simpson'`.	`'hermite-simpson'`
`method`	`str`	The method of `scipy.optimize.minimize`, `'SLSQP'` or `'trust-constr'`.	`'SLSQP'`
`tolerance`	`float`	The tolerance of the optimiser.	`1e-10`
`max_iterations`	`int`	The largest number of iterations of the optimiser.	`1000`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`solution_`		The optimal solution of the run.

Raises:

Type	Description
`ValidationError`	The sense, the scheme, or the method is unknown, or the grid has fewer than 2 nodes.
`SolverError`	The solver has no objective, or the optimiser does not converge.

Examples:

>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.specs import ParameterSpec, StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import CollocationSolver
>>> class PlantElement(Element):
...     def __init__(
...         self,
...         label,
...         x: Annotated[float, StateSpec()] = 1.0,
...         u: Annotated[float, ParameterSpec(adaptive=True, bounds=(-10, 10))] = 0.0,
...     ):
...         super().__init__(label)
...         self.x = x
...         self.u = u
...
...     def _get_laws(self, s):
...         return [RateLaw(s.x, s.u)]
>>> cost = lambda s: (s.plant.x**2 + s.plant.u**2) / 2
>>> solver = CollocationSolver(objective=cost, nodes=51, scheme='hermite-simpson')
>>> simulation = Simulation(System(elements=[PlantElement('plant')]), solver=solver)
>>> simulation = simulation.check().prepare().run(until=1.0)
>>> float(round(simulation.trajectory_.values('plant.u')[0], 3))
-0.762

Source code in src/skcomplex/engine/solvers/_collocation.py

def __init__(
    self,
    objective: Cost | None = None,
    *,
    terminal: Cost = 0.0,
    sense: str = 'min',
    nodes: int = 101,
    scheme: str = 'hermite-simpson',
    method: str = 'SLSQP',
    tolerance: float = 1e-10,
    max_iterations: int = 1000,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.objective = objective
    self.terminal = terminal
    self.sense = sense
    self.nodes = nodes
    self.scheme = scheme
    self.method = method
    self.tolerance = tolerance
    self.max_iterations = max_iterations
    self.random_state = random_state
    self.until = until

`CrankNicolsonSolver(dt=0.01, *, random_state=None, until=math.inf)`

Bases: SplitSolver

Crank-Nicolson for the linear part with a Heun predictor-corrector for the remainder, of second order.

The predictor takes a Crank-Nicolson step with the remainder at the start, and the corrector repeats it with the remainder averaged over the start and the predicted end. For linear fields it is the trapezoidal rule.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed internal step.	`0.01`

Examples:

>>> from skcomplex.engine.solvers import CrankNicolsonSolver
>>> CrankNicolsonSolver(dt=0.5).dt
0.5

Source code in src/skcomplex/engine/solvers/_split.py

def __init__(
    self,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`ETDRK4Solver(dt=0.01, *, random_state=None, until=math.inf)`

Bases: SplitSolver

Fourth-order exponential time differencing Runge-Kutta of Cox and Matthews, for spectral fields.

The linear part is integrated exactly per Fourier mode and the remainder by a fourth-order scheme, with the coefficients computed by contour integrals as Kassam and Trefethen propose. Coupled states of a field use a small matrix exponential per mode. States outside spectral fields follow classical RK4.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed internal step.	`0.01`

Examples:

>>> from skcomplex.engine.solvers import ETDRK4Solver
>>> ETDRK4Solver(dt=0.25).dt
0.25

Source code in src/skcomplex/engine/solvers/_split.py

def __init__(
    self,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`EulerMaruyamaSolver(dt=0.001, *, random_state=None, until=math.inf)`

Bases: SDESolver

Euler-Maruyama integrator of stochastic differential equations with a fixed internal step.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest internal time step.	`0.001`

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.system import System
>>> from skcomplex.models.economics.finance import AssetElement, GeometricBrownianMotionInteraction
>>> from skcomplex.engine.solvers import EulerMaruyamaSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(
...     System(
...         elements=[AssetElement('stock', price=100)],
...         interactions=[GeometricBrownianMotionInteraction('gbm', target='stock', mu=0.05, sigma=0.2)],
...     ),
... )
>>> stepper = EulerMaruyamaSolver(dt=1e-3, random_state=0).check().prepare(problem)
>>> bool(stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_[0] > 0)
True

Source code in src/skcomplex/engine/solvers/_sde.py

def __init__(
    self,
    dt: float = 1e-3,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`ExternalSolver(*, random_state=None, until=math.inf)`

Bases: Solver

Base class of solvers that translate a compiled problem into a third-party library.

Subclasses implement the hooks _get_name, _get_form, _unsupported, _build, and _advance, and list the other forms the library runs in _get_forms. Before the first frame, prepare checks every element, interaction, and law, and raises one error that lists everything the library cannot run.

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`FractionalABMSolver(dt=0.01, *, random_state=None, until=math.inf)`

Bases: CaputoSolver

Fractional Adams-Bashforth-Moulton predictor-corrector of Diethelm, Ford, and Freed, of order 1 + alpha.

The scheme solves the equivalent Volterra equation x(t) = x_0 + (1 / Gamma(alpha)) int_0^t (t - s)^(alpha - 1) f(s, x(s)) ds. A rectangle rule over the drifts of the whole run predicts the next state, and a trapezoidal rule with the drift at the prediction corrects it. States of order 1 take Heun steps.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed step. Frames must fall on multiples of it.	`0.01`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`orders_`		The order of every state component.

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.specs import StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.engine.solvers import FractionalABMSolver
>>> from skcomplex.specs import StepSpec
>>> class RelaxingElement(Element):
...     def __init__(self, label, x: Annotated[float, StateSpec(order=0.5)] = 1.0):
...         super().__init__(label)
...         self.x = x
...
...     def _get_laws(self, s):
...         return [RateLaw(s.x, -s.x)]
>>> problem = compile_system(System(elements=[RelaxingElement('a')]))
>>> stepper = FractionalABMSolver(dt=0.001, random_state=0).check().prepare(problem)
>>> step = StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)
>>> round(float(stepper.run(step).state_[0]), 2)
0.43

Source code in src/skcomplex/engine/solvers/_fractional_abm.py

def __init__(
    self,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`FractionalEulerMaruyamaSolver(dt=0.001, horizon=None, max_increments=_MAX_INCREMENTS, *, random_state=None, until=math.inf)`

Bases: Solver

Euler-Maruyama integrator of stochastic differential equations driven by fractional and ordinary noise.

The stepper draws the correlated increments of every fractional noise row for the whole run when it first advances, on the grid t0 + k dt up to the horizon, and adds the stored increment of each step times the amplitude. Ordinary noise rows take independent Wiener increments. A simulation passes the end time of its run as the horizon when the stepper has none.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed internal step. Frames must fall on multiples of it.	`0.001`
`horizon`	`float \| None`	The last time the run reaches, or `None` for the end time of the simulation's run.	`None`
`max_increments`	`int`	The largest number of stored increments, one per step, fractional noise row, and group member.	`_MAX_INCREMENTS`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`horizon_`		The horizon of the drawn increments.
`hurst_`		The Hurst index of every noise row the increments were drawn with.

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, GroupElement
>>> from skcomplex.core.laws import FractionalNoiseLaw
>>> from skcomplex.specs import StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.engine.solvers import FractionalEulerMaruyamaSolver
>>> from skcomplex.specs import StepSpec
>>> class WalkerElement(Element):
...     def __init__(self, label, x: Annotated[float, StateSpec()] = 0.0):
...         super().__init__(label)
...         self.x = x
...
...     def _get_laws(self, s):
...         return [FractionalNoiseLaw(s.x, 1, hurst=0.75)]
>>> problem = compile_system(System(elements=[GroupElement('walkers', element=WalkerElement, size=500)]))
>>> stepper = FractionalEulerMaruyamaSolver(dt=0.01, horizon=1.0, random_state=0).check().prepare(problem)
>>> state = stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_
>>> round(float(np.var(state)), 1)
1.0

Source code in src/skcomplex/engine/solvers/_fractional_euler_maruyama.py

def __init__(
    self,
    dt: float = 1e-3,
    horizon: float | None = None,
    max_increments: int = _MAX_INCREMENTS,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.horizon = horizon
    self.max_increments = max_increments
    self.random_state = random_state
    self.until = until

`GillespieSolver(*, random_state=None, until=math.inf)`

Bases: Solver

Exact stochastic simulation of jump processes by the Gillespie direct method.

An event that would fire after until is discarded, which is exact because waiting times are memoryless.

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.system import System
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.models.biology.epidemiology import RecoveryInteraction
>>> from skcomplex.engine.solvers import GillespieSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(
...     System(
...         elements=[CompartmentElement('I', n=10), CompartmentElement('R', n=0)],
...         interactions=[RecoveryInteraction('rec', infected='I', recovered='R', gamma=1.0, stochastic=True)],
...     ),
... )
>>> stepper = GillespieSolver(random_state=0).check().prepare(problem)
>>> stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 100.0)).state_.tolist()
[0.0, 10.0]

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`GradientDescentSolver(precision='float64', *, random_state=None, until=math.inf)`

Bases: BaseSolver[LearningProblemSpec]

Solver of the learning form, which trains the parameters of a nested system by gradient steps on JAX and Optax.

Every unit of time is one epoch. The solver shuffles the sequences of the dataset with its random number generator, splits them into batches, and for every batch unrolls the nested system from its initial states over each sequence, writing the values of the dataset into its parameters before every inner step. A discrete map unrolls with jax.lax.scan, and ordinary differential equations with one diffrax.diffeqsolve per inner step, with Tsit5, a PID step size controller, and recursive checkpointing. The loss is the weighted mean of the loss observable over the steps and sequences, its gradient comes from jax.value_and_grad, and an Optax optimizer updates the trained parameters. Truncated backpropagation through time stops the gradient of the carried inner states every truncation steps. The optimizer reads its hyperparameters from the parameter vector before every frame, so a schedule takes effect at the next frame, and its state persists across frames until the solver is prepared again. The solver records the mean batch loss of every epoch and every metric over the whole dataset at the end of the epoch as reporters <interaction path>/loss and <interaction path>/<metric>.

Parameters:

Name	Type	Description	Default
`precision`	`str`	`'float64'` to compute in double precision, inside `jax.enable_x64(True)`, so the unrolled runs agree with NumPy within 1e-10, or `'float32'` for speed.	`'float64'`
`random_state`	`int \| Generator \| None`	The seed or generator of the shuffling of the sequences.	`None`
`until`	`float`	The end time of the run the solver prepares for.	`inf`

Attributes:

Name	Type	Description
`theta_`	`TracedArray`	The trained parameters after the last epoch, as a JAX array.
`optimizer_state_`	`OptState`	The state of the Optax optimizer, such as the moments of Adam.

Raises:

Type	Description
`ValidationError`	The precision is unknown.

Examples:

>>> from skcomplex.engine.solvers import GradientDescentSolver
>>> GradientDescentSolver(precision='float32').check().precision
'float32'

Source code in src/skcomplex/engine/solvers/_gradient_descent.py

def __init__(
    self,
    precision: str = 'float64',
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    super().__init__(random_state=random_state, until=until)
    self.precision = precision

`GrunwaldLetnikovSolver(dt=0.01, memory=None, *, random_state=None, until=math.inf)`

Bases: ConvolutionSolver

Explicit Grünwald-Letnikov scheme of Caputo fractional differential equations, of order 1.

A step sets dt^alpha f(t_(k-1), x_(k-1)) = sum_(j=0)^k w_j (x_(k-j) - x_0) with the weights w_j = (-1)^j binomial(alpha, j). States of order 1 take explicit Euler steps.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed step. Frames must fall on multiples of it.	`0.01`
`memory`	`float \| None`	The length of time over which the history enters with its exact weights, or `None` for the whole run. The history older than the window enters through exponential modes fitted to the weights, so a step costs the same however long the run.	`None`

Attributes:

Name	Type	Description
`orders_`		The order of every state component.

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.specs import StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.engine.solvers import GrunwaldLetnikovSolver
>>> from skcomplex.specs import StepSpec
>>> class RelaxingElement(Element):
...     def __init__(self, label, x: Annotated[float, StateSpec(order=0.5)] = 1.0):
...         super().__init__(label)
...         self.x = x
...
...     def _get_laws(self, s):
...         return [RateLaw(s.x, -s.x)]
>>> problem = compile_system(System(elements=[RelaxingElement('a')]))
>>> stepper = GrunwaldLetnikovSolver(dt=0.001, random_state=0).check().prepare(problem)
>>> step = StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)
>>> round(float(stepper.run(step).state_[0]), 2)
0.43

Source code in src/skcomplex/engine/solvers/_convolution.py

def __init__(
    self,
    dt: float = 0.01,
    memory: float | None = None,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.memory = memory
    self.random_state = random_state
    self.until = until

`HeunSolver(dt=0.001, *, random_state=None, until=math.inf)`

Bases: SDESolver

Predictor-corrector scheme for SDEs that averages the drift over the step and keeps the Itô noise.

The predictor X~ = X + f(X) h + g(X) dW gives the corrector X + (f(X) + f(X~)) h / 2 + g(X) dW, which has strong order 1 for additive noise and integrates the Itô form for any noise.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed internal step.	`0.001`

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.models.economics.finance import geometric_brownian_motion
>>> from skcomplex.engine.solvers import HeunSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(geometric_brownian_motion())
>>> stepper = HeunSolver(dt=0.01, random_state=0).check().prepare(problem)
>>> bool(stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_[0] > 0)
True

Source code in src/skcomplex/engine/solvers/_sde.py

def __init__(
    self,
    dt: float = 1e-3,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`HorizonError`

Bases: SolverError

A stepper with a fixed horizon is asked to advance beyond it.

`IMEX2Solver(dt=0.01, *, random_state=None, until=math.inf)`

Bases: SplitSolver

Crank-Nicolson Adams-Bashforth scheme of second order: trapezoidal linear part and extrapolated remainder.

Every step solves (I - h L / 2) x_{n+1} = (I + h L / 2) x_n + h (3 N_n - N_{n-1}) / 2, starting with one Crank-Nicolson step with an Euler remainder. A call that does not continue the time, state, and step of the previous call starts the scheme again.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed internal step.	`0.01`

Examples:

>>> from skcomplex.engine.solvers import IMEX2Solver
>>> IMEX2Solver(dt=0.1).dt
0.1

Source code in src/skcomplex/engine/solvers/_split.py

def __init__(
    self,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`IMEXEulerSolver(dt=0.01, *, random_state=None, until=math.inf)`

Bases: SplitSolver

Implicit-explicit Euler: the linear part implicit and the nonlinear remainder explicit, of first order.

Every step solves (I - h L) x_{n+1} = x_n + h N(x_n), with sparse solves for finite differences and finite elements and division per Fourier mode for spectral fields.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed internal step.	`0.01`

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.models.physics.heat import heat_diffusion
>>> from skcomplex.engine.solvers import IMEXEulerSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(heat_diffusion((16, 16)))
>>> stepper = IMEXEulerSolver(dt=1.0, random_state=0).check().prepare(problem)
>>> state = stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 10.0)).state_
>>> bool(np.isclose(state.sum(), problem.initial_state.sum()))
True

Source code in src/skcomplex/engine/solvers/_split.py

def __init__(
    self,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`IterationSolver(*, random_state=None, until=math.inf)`

Bases: Solver

Iteration of a discrete map, one application per unit of time.

A stochastic map takes its random draws from the stepper's random number generator at every step.

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.system import System
>>> from skcomplex.models.economics.markets import CobwebInteraction, MarketElement
>>> from skcomplex.engine.solvers import IterationSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(
...     System(
...         elements=[MarketElement('m', price=1)],
...         interactions=[CobwebInteraction('c', target='m', a=10, b=2, c=1, d=1)],
...     ),
... )
>>> stepper = IterationSolver(random_state=0).check().prepare(problem)
>>> stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_.tolist()
[4.0]

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`JumpTrajectoriesSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)`

Bases: OperatorSolver

Solver that samples the master equation by jumps of a pure state with QuTiP's MCSolver.

Every run is one trajectory, whose jumps are drawn from the simulation's random number generator, so a seeded run is reproducible. A problem stored as a density matrix records the projector on the pure state of every frame.

Parameters:

Name	Type	Description	Default
`atol`	`float`	The absolute tolerance of the integrator.	`1e-12`
`rtol`	`float`	The relative tolerance of the integrator.	`1e-10`
`nsteps`	`int`	The largest number of internal steps between two frames.	`100000`

Raises:

Type	Description
`SolverError`	The initial state is a mixed density matrix.

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from sage.all__sagemath_symbolics import CC, SR
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import DissipatorLaw, GeneratorLaw
>>> from skcomplex.core.laws.utils import local_operator
>>> from skcomplex.specs import HilbertSpec, ObservableSpec, StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import JumpTrajectoriesSolver
>>> operators = {'x': (np.array([[0, 1], [1, 0]]), True), 'z': (np.diag([1, -1]), True)}
>>> operators['minus'] = (np.array([[0, 0], [1, 0]]), False)
>>> QUBIT = HilbertSpec(2, operators, basis={'up': 0, 'down': 1})
>>> class TwoLevelElement(Element):
...     def __init__(self, label, psi: Annotated[object, StateSpec(domain=CC**2)] = 'up', hilbert=QUBIT):
...         super().__init__(label)
...         self.psi = psi
...         self.hilbert = hilbert
...
...     def _get_laws(self, s):
...         return [DissipatorLaw(local_operator(s.psi, SR.symbol('minus')), 1.0)]
>>> z = [ObservableSpec('z', element='q', variable='z')]
>>> system = System(elements=[TwoLevelElement('q')])
>>> solver = JumpTrajectoriesSolver()
>>> from skcomplex.engine.compile import OperatorCompiler
>>> simulation = Simulation(
...     system, compiler=OperatorCompiler(), solver=solver, dt=1.0, random_state=0, observables=z,
... ).check().prepare().run(until=20.0)
>>> float(simulation.trajectory_.values('q.z')[-1])
-1.0

Source code in src/skcomplex/engine/solvers/_operator.py

def __init__(
    self,
    atol: float = 1e-12,
    rtol: float = 1e-10,
    nsteps: int = 100_000,
    max_dimension: int = DEFAULT_MAX_DIMENSION,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.atol = atol
    self.rtol = rtol
    self.nsteps = nsteps
    self.max_dimension = max_dimension
    self.random_state = random_state
    self.until = until

`KineticSolver(n=10000, dt=0.01, *, random_state=None, until=math.inf)`

Bases: Solver

Direct simulation Monte Carlo of encounter laws with finite particles on activity grids.

Every density is held as counts of particles of one common mass, chosen so that the largest initial density has n particles. Every step of length at most dt applies the encounter laws in turn. A candidate particle at activity h encounters with probability 1 - exp(-lambda_h dt), where lambda_h = sum_k eta[h, k] g_k over the field density g, and picks its field activity and outcome in proportion to eta[h, k] A[h, k, i] g_k. A proliferation adds a Poisson number of particles and a destruction removes a binomial number. The state holds the empirical densities, so the trajectories match those of the gain and loss equations, which the particles approach as n grows.

Parameters:

Name	Type	Description	Default
`n`	`int`	The number of particles of the largest initial density.	`10000`
`dt`	`float`	The largest time step.	`0.01`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Examples:

>>> from skcomplex.engine.solvers import KineticSolver
>>> from skcomplex.core.laws import BilinearLaw, EncounterLaw
>>> KineticSolver(n=1000)._get_support().laws == {EncounterLaw, BilinearLaw}
True

Source code in src/skcomplex/engine/solvers/_kinetics.py

def __init__(
    self,
    n: int = 10_000,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.n = n
    self.dt = dt
    self.random_state = random_state
    self.until = until

`L1Solver(dt=0.01, memory=None, *, random_state=None, until=math.inf)`

Bases: ConvolutionSolver

Explicit L1 scheme of Caputo fractional differential equations, of order 1.

The L1 scheme interpolates the state linearly between grid points inside the Caputo integral. A step sets Gamma(2 - alpha) dt^alpha f(t_(k-1), x_(k-1)) = sum_(j=0)^(k-1) b_j (x_(k-j) - x_(k-j-1)) with the weights b_j = (j + 1)^(1 - alpha) - j^(1 - alpha). States of order 1 take explicit Euler steps.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed step. Frames must fall on multiples of it.	`0.01`
`memory`	`float \| None`	The length of time over which the history enters with its exact weights, or `None` for the whole run. The history older than the window enters through exponential modes fitted to the weights, so a step costs the same however long the run.	`None`

Attributes:

Name	Type	Description
`orders_`		The order of every state component.

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.specs import StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.engine.solvers import L1Solver
>>> from skcomplex.specs import StepSpec
>>> class RelaxingElement(Element):
...     def __init__(self, label, x: Annotated[float, StateSpec(order=0.5)] = 1.0):
...         super().__init__(label)
...         self.x = x
...
...     def _get_laws(self, s):
...         return [RateLaw(s.x, -s.x)]
>>> problem = compile_system(System(elements=[RelaxingElement('a')]))
>>> stepper = L1Solver(dt=0.001, random_state=0).check().prepare(problem)
>>> step = StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)
>>> round(float(stepper.run(step).state_[0]), 2)
0.43

Source code in src/skcomplex/engine/solvers/_convolution.py

def __init__(
    self,
    dt: float = 0.01,
    memory: float | None = None,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.memory = memory
    self.random_state = random_state
    self.until = until

`LammpsFixSpec(name, group, value, mass, seed)` `dataclass`

Bases: BaseSpec

Values the native fix of an interaction reads when the LAMMPS backend writes its command.

Parameters:

Name	Type	Description	Default
`name`	`str`	The unique name of the fix.	required
`group`	`str`	The LAMMPS group of the atoms the fix acts on.	required
`value`	`Callable[[str], Array]`	The function that returns the current value of a parameter of the interaction by name.	required
`mass`	`Callable[[], float]`	The function that returns the mass of the atoms of the interaction's target.	required
`seed`	`Callable[[], int]`	The function that draws a random seed.	required

`LammpsSolver(dt=0.005, pairs=None, fixes=None, *, random_state=None, until=math.inf)`

Bases: ExternalSolver

Solver that runs systems of points in LAMMPS, with LAMMPS velocity Verlet integration.

An element or group whose members declare a vector position state of order 2, such as PointElement, becomes an atom type, whose atoms move with the first derivative of the position as their velocity. The mass of an atom type is the parameter that the divisor of the interactions on its members names, and 1 without a divisor. Boxes become LAMMPS boundaries and fix wall/reflect. A PairCouplingInteraction with a potential becomes the native pair style pairs gives for its class, or a Lepton expression potential otherwise. An interaction whose class fixes gives a native fix becomes that fix. Any other interaction that only changes the derivatives of the positions runs through fix external, evaluated by the framework.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest LAMMPS time step.	`0.005`
`pairs`	`Mapping[type[BaseInteraction], tuple[str, PairArguments]] \| None`	The native LAMMPS pair style of pair potential classes, with the function of the parameter values that returns the coefficients of `pair_coeff` without the cutoff, which the solver appends.	`None`
`fixes`	`Mapping[type[BaseInteraction], FixCommand] \| None`	The function that returns the native LAMMPS fix command of interaction classes from the values of the fix.	`None`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`lammps_`		The LAMMPS instance, which the final run `run(..., final=True)` closes.
`atom_elements_`		The elements and groups of points, each a LAMMPS atom type.
`atom_numbers_`		The LAMMPS atom type of every element and group of points, by path.
`position_indices_`		The flat state indices of the position of every atom, of shape `(atoms, dimensions)`.
`velocity_indices_`		The flat state indices of the velocity of every atom, of shape `(atoms, dimensions)`.
`atom_types_`		The LAMMPS atom type of every atom.
`n_dim_`		The number of dimensions of the space of the atoms.
`released_`		Whether the LAMMPS instance is closed.

Examples:

>>> from skcomplex.core.elements import GroupElement, PointElement
>>> from skcomplex.core.interactions import PairCouplingInteraction
>>> from skcomplex.core.spaces import ContinuousSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import LammpsSolver
>>> potential = lambda namespace, r: 4 * (r**-12 - r**-6)
>>> system = System(
...     elements=[
...         GroupElement('gas', element=PointElement, size=2, values={'x': [[0.0, 0.0], [1.5, 0.0]]}),
...     ],
...     interactions=[PairCouplingInteraction('lj', 'gas', 'gas', cutoff=2.5, potential=potential)],
...     space=ContinuousSpace(n_dim=2),
... )
>>> solver = LammpsSolver(dt=0.005)
>>> trajectory = Simulation(system, solver=solver, dt=0.1).check().prepare().run(until=0.5).trajectory_
>>> positions = trajectory.values('gas.x')
>>> bool(positions[-1, 1, 0] - positions[-1, 0, 0] < 1.5)
True

Source code in src/skcomplex/engine/solvers/_lammps.py

def __init__(
    self,
    dt: float = 0.005,
    pairs: Mapping[type[BaseInteraction], tuple[str, PairArguments]] | None = None,
    fixes: Mapping[type[BaseInteraction], FixCommand] | None = None,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.pairs = pairs
    self.fixes = fixes
    self.random_state = random_state
    self.until = until

`LatticeSolver(*, random_state=None, until=math.inf)`

Bases: Solver

Solver of lattice dynamics that applies every local rule, drive, cascade, and extremal law once per step.

Steps apply in the order the laws are declared, and the boundary laws apply after every step. Cascades record their size, duration, and area as events, and extremal laws record the site and the value of every application.

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.models.computer_science.automata import elementary_automaton
>>> from skcomplex.engine.solvers import LatticeSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(elementary_automaton(30, size=11))
>>> stepper = LatticeSolver(random_state=0).check().prepare(problem)
>>> step = StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)
>>> stepper.run(step).state_.astype(int).tolist()
[0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0]

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`LindbladSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)`

Bases: OperatorSolver

Solver of the Lindblad master equation with QuTiP's MESolver, and of pure closed systems with SESolver.

It is the default solver of the operator form. A system without dissipators that starts in a pure state follows the linear evolution dx/dt = -i H x, which is the master equation without dissipators.

Parameters:

Name	Type	Description	Default
`atol`	`float`	The absolute tolerance of the integrator.	`1e-12`
`rtol`	`float`	The relative tolerance of the integrator.	`1e-10`
`nsteps`	`int`	The largest number of internal steps between two frames.	`100000`

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from sage.all__sagemath_symbolics import CC, SR
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import DissipatorLaw, GeneratorLaw
>>> from skcomplex.core.laws.utils import local_operator
>>> from skcomplex.specs import HilbertSpec, ObservableSpec, StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import LindbladSolver
>>> operators = {'x': (np.array([[0, 1], [1, 0]]), True), 'z': (np.diag([1, -1]), True)}
>>> operators['minus'] = (np.array([[0, 0], [1, 0]]), False)
>>> QUBIT = HilbertSpec(2, operators, basis={'up': 0, 'down': 1})
>>> class TwoLevelElement(Element):
...     def __init__(self, label, psi: Annotated[object, StateSpec(domain=CC**2)] = 'up', hilbert=QUBIT):
...         super().__init__(label)
...         self.psi = psi
...         self.hilbert = hilbert
...
...     def _get_laws(self, s):
...         return [DissipatorLaw(local_operator(s.psi, SR.symbol('minus')), 0.5)]
>>> z = [ObservableSpec('z', element='q', variable='z')]
>>> system = System(elements=[TwoLevelElement('q')])
>>> from skcomplex.engine.compile import OperatorCompiler
>>> compiler = OperatorCompiler()
>>> simulation = Simulation(system, compiler=compiler, solver=LindbladSolver(), dt=1.0, observables=z)
>>> simulation = simulation.check().prepare()
>>> values = simulation.run(until=2.0).trajectory_.values('q.z')
>>> bool(abs((1 + values[-1]) / 2 - np.exp(-1.0)) < 1e-6)
True

Source code in src/skcomplex/engine/solvers/_operator.py

def __init__(
    self,
    atol: float = 1e-12,
    rtol: float = 1e-10,
    nsteps: int = 100_000,
    max_dimension: int = DEFAULT_MAX_DIMENSION,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.atol = atol
    self.rtol = rtol
    self.nsteps = nsteps
    self.max_dimension = max_dimension
    self.random_state = random_state
    self.until = until

`MesaSolver(*, random_state=None, until=math.inf)`

Bases: ExternalSolver

Mesa backend of agent-based problems, one Mesa step per unit of time.

The agents are the members of the elements and groups whose laws include RuleLaws, and of the variable groups. Every step advances Mesa's clock by one, activates all agents in random order, or by the phases of the system's discrete time, and each agent runs the rules of its element type in declaration order. Interaction rules run once after the agents. Rules add and remove members of variable groups only. The class variable supports declares what the backend runs, which validation checks systems against: discrete time, no space or a lattice, network, or continuous space, rule, boundary, transfer, and matching laws, and elements and fixed and variable groups with any domain but Hilbert spaces.

Attributes:

Name	Type	Description
`model_`		The Mesa model.
`members_`		Every agent created in the run by its id, which reference states name.

Examples:

>>> from typing import Annotated
>>> from skcomplex.core.elements import Element, GroupElement
>>> from skcomplex.core.laws import RuleLaw
>>> from skcomplex.specs import StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import MesaSolver
>>> class CounterElement(Element):
...     def __init__(self, label, count: Annotated[float, StateSpec()] = 0.0):
...         super().__init__(label)
...         self.count = count
...
...     def _get_laws(self, namespace):
...         return [RuleLaw(self.tick)]
...
...     def tick(self, agent, model):
...         agent.count += 1
>>> system = System(elements=[GroupElement('counters', element=CounterElement, size=3)])
>>> trajectory = Simulation(system, solver=MesaSolver()).check().prepare().run(until=4).trajectory_
>>> trajectory.agents('counters', 4).values['count'].tolist()
[4.0, 4.0, 4.0]

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`MilsteinSolver(dt=0.001, *, random_state=None, until=math.inf)`

Bases: SDESolver

Derivative-free Milstein scheme of Platen for SDEs with diagonal noise, of strong order 1.

It takes a supporting value Y = X + f h + g sqrt(h) and corrects the noise with (g(Y) - g(X)) (dW^2 - h) / (2 sqrt(h)), which needs no derivative of the noise.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The fixed internal step.	`0.001`

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.models.economics.finance import geometric_brownian_motion
>>> from skcomplex.engine.solvers import MilsteinSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(geometric_brownian_motion())
>>> stepper = MilsteinSolver(dt=0.01, random_state=0).check().prepare(problem)
>>> bool(stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_[0] > 0)
True

Source code in src/skcomplex/engine/solvers/_sde.py

def __init__(
    self,
    dt: float = 1e-3,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`ModelView(stepper, owner='', roles=None)`

Model as rules see it: parameters, random numbers, time, agents, and the space.

Every rule receives the view scoped to the interaction or agent group that owns it, so the rule reads the owner's parameters by name and the groups bound to the owner's roles by role, also when its system is nested in another.

Parameters:

Name	Type	Description	Default
`stepper`	`MesaSolver`	The stepper that runs the model.	required
`owner`	`str`	The path of the interaction or agent group whose rules receive the view, or `''` for no owner.	`''`
`roles`	`Mapping[str, str \| tuple[str, ...]] \| None`	The paths of the elements and groups the owner's roles bind to, by role.	`None`

Source code in src/skcomplex/engine/solvers/_mesa.py

def __init__(
    self,
    stepper: 'MesaSolver',
    owner: str = '',
    roles: Mapping[str, str | tuple[str, ...]] | None = None,
) -> None:
    self._stepper = stepper
    self._owner = owner
    self._roles = dict(roles or {})

`dimensions` `property`

Extent of the agent space per axis, cells of a grid or lengths of a continuous space, None otherwise.

`inputs` `property`

Values of the states of other partitions at the last synchronisation point, by dotted name.

`parameters` `property`

Current values of every parameter by dotted name.

`rng` `property`

Random number generator of the model, seeded by the simulation's random state.

`round` `property`

Index of the current round of a repeated phase, 0 outside repeated phases.

`storage` `property`

Mutable mapping private to the owner that persists across the steps of one run, such as an order book.

`time` `property`

Time of Mesa's clock, k during the step that produces frame k, plus i / n at visit i of n.

`torus` `property`

Whether the agent space wraps around at its edges.

`agents(name, view=None)`

agents(name: str) -> list[AgentView]

agents(name: str, view: type[V]) -> list[V]

Return the views of the agents of a group bound to a role of the owner, or of a group by path.

Parameters:

Name	Type	Description	Default
`name`	`str`	A role of the owner, the path of a group relative to the owner's system, or the path of a group.	required
`view`	`type[V] \| None`	The view class of the group, which its element names with `_get_view`, or `None` for `AgentView`.	`None`

Returns:

Name	Type	Description
`agents`	`list[V] \| list[AgentView]`	The views of the agents alive in the group, or in every group of a role that binds several.

Raises:

Type	Description
`SimulationError`	No role or group has that name, or the views of a group are not of the view class.

Source code in src/skcomplex/engine/solvers/_mesa.py

def agents[V: AgentView](self, name: str, view: type[V] | None = None) -> list[V] | list['AgentView']:
    """Return the views of the agents of a group bound to a role of the owner, or of a group by path.

    Args:
        name:
            A role of the owner, the path of a group relative to the owner's system, or the path of a group.

        view:
            The view class of the group, which its element names with `_get_view`, or `None` for `AgentView`.

    Returns:
        agents:
            The views of the agents alive in the group, or in every group of a role that binds several.

    Raises:
        SimulationError:
            No role or group has that name, or the views of a group are not of the view class.
    """
    stepper = self._stepper
    views = [MesaSolver._view_of(agent) for path in self._paths(name) for agent in stepper._group_agents(path)]
    wanted = AgentView if view is None else view
    wrong = [item for item in views if not isinstance(item, wanted)]
    if wrong:
        viewed = type(wrong[0]).__mro__[1].__name__
        error_msg = f'The agents of {name} are viewed as {viewed}, not {wanted.__name__}.'
        raise SimulationError(error_msg)
    if view is None:
        return views
    return [item for item in views if isinstance(item, view)]

`composition(name, state)`

Return the number of agents of a group that hold every value of a state, for the values some agent holds.

The stepper keeps the counts as agents are created, removed, and change the state, so reading them costs no pass over the agents.

Parameters:

Name	Type	Description	Default
`name`	`str`	A role of the owner, the path of a group relative to the owner's system, or the path of a group.	required
`state`	`str`	The name of a scalar state of the agents of the group.	required

Returns:

Name	Type	Description
`counts`	`dict[Hashable, int]`	The number of agents alive in the group, or in every group of a role that binds several, by value.

Raises:

Type	Description
`SimulationError`	No role or group has that name, or its agents have no such state.

Source code in src/skcomplex/engine/solvers/_mesa.py

def composition(self, name: str, state: str) -> dict[Hashable, int]:
    """Return the number of agents of a group that hold every value of a state, for the values some agent holds.

    The stepper keeps the counts as agents are created, removed, and change the state, so reading them costs no
    pass over the agents.

    Args:
        name:
            A role of the owner, the path of a group relative to the owner's system, or the path of a group.

        state:
            The name of a scalar state of the agents of the group.

    Returns:
        counts:
            The number of agents alive in the group, or in every group of a role that binds several, by value.

    Raises:
        SimulationError:
            No role or group has that name, or its agents have no such state.
    """
    counts: dict[Hashable, int] = {}
    for path in self._paths(name):
        for value, count in self._stepper._composition(path, state).items():
            counts[value] = counts.get(value, 0) + count
    return {value: count for value, count in counts.items() if count}

`count(name)`

Return the number of agents of a group bound to a role of the owner, or of a group by path.

Parameters:

Name	Type	Description	Default
`name`	`str`	A role of the owner, the path of a group relative to the owner's system, or the path of a group.	required

Returns:

Name	Type	Description
`count`	`int`	The number of agents alive in the group, or in every group of a role that binds several.

Raises:

Type	Description
`SimulationError`	No role or group has that name.

Source code in src/skcomplex/engine/solvers/_mesa.py

def count(self, name: str) -> int:
    """Return the number of agents of a group bound to a role of the owner, or of a group by path.

    Args:
        name:
            A role of the owner, the path of a group relative to the owner's system, or the path of a group.

    Returns:
        count:
            The number of agents alive in the group, or in every group of a role that binds several.

    Raises:
        SimulationError:
            No role or group has that name.
    """
    return sum(len(self._stepper._group_agents(path)) for path in self._paths(name))

`create(group, location=None, **values)`

Create an agent of a group, at a location or at a random place.

Parameters:

Name	Type	Description	Default
`group`	`str`	A role of the owner that binds the group, or the path of the group.	required
`location`	`object`	The cell coordinates, network node, or position, or `None` for a random place.	`None`
`**values`	`object`	The initial values of states, with the element type's defaults for the others.	`{}`

Returns:

Name	Type	Description
`agent`	`AgentView`	The new agent.

Raises:

Type	Description
`SimulationError`	The group is not variable, so its size cannot change.

Source code in src/skcomplex/engine/solvers/_mesa.py

def create(self, group: str, location: object = None, **values: object) -> 'AgentView':
    """Create an agent of a group, at a location or at a random place.

    Args:
        group:
            A role of the owner that binds the group, or the path of the group.

        location:
            The cell coordinates, network node, or position, or `None` for a random place.

        **values:
            The initial values of states, with the element type's defaults for the others.

    Returns:
        agent:
            The new agent.

    Raises:
        SimulationError:
            The group is not variable, so its size cannot change.
    """
    return MesaSolver._view_of(self._stepper._create_agent(self._paths(group)[0], location, values))

`deposit(name, cell, amount)`

Add an amount to a field at a cell at the next synchronisation point.

Parameters:

Name	Type	Description	Default
`name`	`str`	The dotted path of the field state, such as `ground/pheromone.level`, or `<role>.<state>` for a field bound to a role of the owner.	required
`cell`	`object`	The lattice coordinates of the cell or a network node.	required
`amount`	`float`	The amount to add.	required

Raises:

Type	Description
`KeyError`	The name is not a state of a field.

Source code in src/skcomplex/engine/solvers/_mesa.py

def deposit(self, name: str, cell: object, amount: float) -> None:
    """Add an amount to a field at a cell at the next synchronisation point.

    Args:
        name:
            The dotted path of the field state, such as `ground/pheromone.level`, or `<role>.<state>` for a field
            bound to a role of the owner.

        cell:
            The lattice coordinates of the cell or a network node.

        amount:
            The amount to add.

    Raises:
        KeyError:
            The name is not a state of a field.
    """
    name = self._field_name(name)
    site = self._site(name, cell)
    entry = self._stepper.problem_.states.entry(name)
    self._stepper.deposits_[entry.start + site * entry.columns] += amount

`field(name, cell)`

Return the value of a field at a cell, as it was at the last synchronisation point.

Parameters:

Name	Type	Description	Default
`name`	`str`	The dotted path of the field state, such as `ground/pheromone.level`, or `<role>.<state>` for a field bound to a role of the owner.	required
`cell`	`object`	The lattice coordinates of the cell, which match the coordinates of a grid of the same shape, or a network node.	required

Returns:

Name	Type	Description
`value`	`float`	The value of the field there.

Raises:

Type	Description
`KeyError`	The name is not a state of a field.

Source code in src/skcomplex/engine/solvers/_mesa.py

def field(self, name: str, cell: object) -> float:
    """Return the value of a field at a cell, as it was at the last synchronisation point.

    Args:
        name:
            The dotted path of the field state, such as `ground/pheromone.level`, or `<role>.<state>` for a field
            bound to a role of the owner.

        cell:
            The lattice coordinates of the cell, which match the coordinates of a grid of the same shape, or a
            network node.

    Returns:
        value:
            The value of the field there.

    Raises:
        KeyError:
            The name is not a state of a field.
    """
    name = self._field_name(name)
    site = self._site(name, cell)
    values = self._stepper.inputs_.get(name)
    if values is None:
        error_msg = f'The agents read the field {name}, which no other partition provides. Run it in co-simulation.'
        raise KeyError(error_msg)
    return float(np.ravel(np.asarray(values))[site])

`groups(name)`

Return the paths of the groups a role of the owner binds, or of a group by relative or absolute path.

Parameters:

Name	Type	Description	Default
`name`	`str`	A role of the owner, the path of a group relative to the owner's system, or the path of a group.	required

Returns:

Name	Type	Description
`paths`	`list[str]`	The paths of the groups.

Raises:

Type	Description
`SimulationError`	No role or group has that name.

Source code in src/skcomplex/engine/solvers/_mesa.py

def groups(self, name: str) -> list[str]:
    """Return the paths of the groups a role of the owner binds, or of a group by relative or absolute path.

    Args:
        name:
            A role of the owner, the path of a group relative to the owner's system, or the path of a group.

    Returns:
        paths:
            The paths of the groups.

    Raises:
        SimulationError:
            No role or group has that name.
    """
    return self._paths(name)

`history(name, window=None)`

Return the past values of a reporter, or of an observable at the end of each earlier step, the last one last.

Parameters:

Name	Type	Description	Default
`name`	`str`	The name of a value the owner reports, the full name of any reporter, such as `bar.attendance`, or the path of an observable of the run, such as `attendance`.	required
`window`	`int \| None`	The number of last values, or `None` for every recorded value.	`None`

Returns:

Name	Type	Description
`values`	`NDArray[float64]`	The values, at most `window` of them.

Source code in src/skcomplex/engine/solvers/_mesa.py

def history(self, name: str, window: int | None = None) -> npt.NDArray[np.float64]:
    """Return the past values of a reporter, or of an observable at the end of each earlier step, the last one last.

    Args:
        name:
            The name of a value the owner reports, the full name of any reporter, such as `bar.attendance`, or the
            path of an observable of the run, such as `attendance`.

        window:
            The number of last values, or `None` for every recorded value.

    Returns:
        values:
            The values, at most `window` of them.
    """
    stepper = self._stepper
    history = stepper.reported_
    values = history.get(self._report_name(name), history.get(name, stepper.observed_.get(name, ())))
    array = (
        np.asarray(values, dtype=float)
        if all(isinstance(item, numbers.Real) for item in values)
        else np.fromiter(values, dtype=object, count=len(values))
    )
    return array if window is None else array[len(array) - min(window, len(array)) :]

`observable(name)`

Return the current value of an observable of the system, such as an aggregate of a group.

Parameters:

Name	Type	Description	Default
`name`	`str`	The path of the observable.	required

Returns:

Name	Type	Description
`value`	`object`	Its value for the current agents and parameters.

Raises:

Type	Description
`SimulationError`	The system has no such observable.

Source code in src/skcomplex/engine/solvers/_mesa.py

def observable(self, name: str) -> object:
    """Return the current value of an observable of the system, such as an aggregate of a group.

    Args:
        name:
            The path of the observable.

    Returns:
        value:
            Its value for the current agents and parameters.

    Raises:
        SimulationError:
            The system has no such observable.
    """
    return self._stepper._observable(name)

`occupancy(name)`

Return the number of agents of a group in every cell or node that holds one, by the key of the location.

Parameters:

Name	Type	Description	Default
`name`	`str`	A role of the owner, the path of a group relative to the owner's system, or the path of a group.	required

Returns:

Name	Type	Description
`counts`	`dict[Hashable, int]`	The number of agents alive in the group, or in every group of a role that binds several, by the `location_key` of their cell, node, or position.

Raises:

Type	Description
`SimulationError`	No role or group has that name.

Source code in src/skcomplex/engine/solvers/_mesa.py

def occupancy(self, name: str) -> dict[Hashable, int]:
    """Return the number of agents of a group in every cell or node that holds one, by the key of the location.

    Args:
        name:
            A role of the owner, the path of a group relative to the owner's system, or the path of a group.

    Returns:
        counts:
            The number of agents alive in the group, or in every group of a role that binds several, by the
            `location_key` of their cell, node, or position.

    Raises:
        SimulationError:
            No role or group has that name.
    """
    counts: dict[Hashable, int] = {}
    for path in self._paths(name):
        for agent in self._stepper._group_agents(path):
            key = self._stepper._view_of(agent).location_key
            counts[key] = counts.get(key, 0) + 1
    return counts

`parameter(name)`

Return the current value of a parameter of the owner, or of one component as name[i].

Parameters:

Name	Type	Description	Default
`name`	`str`	The name of a scalar parameter, such as `threshold`, or of one component, such as `weights[1]`.	required

Returns:

Name	Type	Description
`value`	`float`	The value.

Raises:

Type	Description
`SimulationError`	The owner has no such parameter, or the parameter has several components and no index.

Source code in src/skcomplex/engine/solvers/_mesa.py

def parameter(self, name: str) -> float:
    """Return the current value of a parameter of the owner, or of one component as `name[i]`.

    Args:
        name:
            The name of a scalar parameter, such as `threshold`, or of one component, such as `weights[1]`.

    Returns:
        value:
            The value.

    Raises:
        SimulationError:
            The owner has no such parameter, or the parameter has several components and no index.
    """
    stepper = self._stepper
    layout = stepper.problem_.parameters
    try:
        entry, indices = layout.locate(f'{self._owner}.{name}')
    except KeyError as error:
        prefix = f'{self._owner}.'
        names = [item.name[len(prefix) :] for item in layout.entries if item.name.startswith(prefix)]
        error_msg = (
            f'{self._owner or "The model"} has no parameter {name!r}. Its parameters are: '
            f'{", ".join(names) or "none"}.'
        )
        raise SimulationError(error_msg) from error
    if entry.size != 1:
        error_msg = f'{entry.name} has {entry.size} components. Read one of them as {name}[i].'
        raise SimulationError(error_msg)
    return float(stepper.parameters_[indices][0])

`random_empty_cell()`

Return the coordinates of a random empty cell of the grid or network.

Returns:

Name	Type	Description
`location`	`object`	The cell coordinates or the network node.

Source code in src/skcomplex/engine/solvers/_mesa.py

def random_empty_cell(self) -> object:
    """Return the coordinates of a random empty cell of the grid or network.

    Returns:
        location:
            The cell coordinates or the network node.
    """
    return self._stepper._cells().select_random_empty_cell().coordinate

`report(name, value)`

Publish a value of the owner as the reporter <owner>.<name> and add it to the reporter's history.

The owner declares the names it reports in its class variable reports, so every frame records them, NaN before the first report.

Parameters:

Name	Type	Description	Default
`name`	`str`	The name of the value, such as `price`.	required
`value`	`float`	The value.	required

Raises:

Type	Description
`SimulationError`	The owner does not declare the name in `reports`.

Source code in src/skcomplex/engine/solvers/_mesa.py

def report(self, name: str, value: float) -> None:
    """Publish a value of the owner as the reporter `<owner>.<name>` and add it to the reporter's history.

    The owner declares the names it reports in its class variable `reports`, so every frame records them, `NaN`
    before the first report.

    Args:
        name:
            The name of the value, such as `price`.

        value:
            The value.

    Raises:
        SimulationError:
            The owner does not declare the name in `reports`.
    """
    full = self._report_name(name)
    if full not in self._stepper.reports_:
        error_msg = (
            f'{self._owner or "The model"} reports {name!r}, which it does not declare. List it in the class '
            'variable `reports` of the interaction or element type.'
        )
        raise SimulationError(error_msg)
    kept = float(value) if isinstance(value, numbers.Real) else value
    self._stepper.reports_[full] = kept
    self._stepper.reported_.setdefault(full, deque(maxlen=_HISTORY_LENGTH)).append(kept)

`set_parameter(name, value)`

Set a parameter of any group or interaction during the run, recorded in the trajectory like any parameter.

Parameters:

Name	Type	Description	Default
`name`	`str`	The dotted name of the parameter, such as `bank.rate`, or the name of a parameter of the owner.	required
`value`	`float`	The new value.	required

Raises:

Type	Description
`SimulationError`	No parameter has that name.

Source code in src/skcomplex/engine/solvers/_mesa.py

def set_parameter(self, name: str, value: float) -> None:
    """Set a parameter of any group or interaction during the run, recorded in the trajectory like any parameter.

    Args:
        name:
            The dotted name of the parameter, such as `bank.rate`, or the name of a parameter of the owner.

        value:
            The new value.

    Raises:
        SimulationError:
            No parameter has that name.
    """
    layout = self._stepper.problem_.parameters
    full = name if layout.has(name) else f'{self._owner}.{name}'
    if not layout.has(full):
        error_msg = f'No parameter named {name!r} to set. Available names are: {", ".join(layout.names)}.'
        raise SimulationError(error_msg)
    self._stepper._set_parameter(full, value)

`MonteCarloSolver(*, random_state=None, until=math.inf)`

Bases: LatticeSolver

Solver of Monte Carlo sampling that sweeps the lattice once per unit of time.

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.models.physics.statistical_mechanics import ising
>>> from skcomplex.engine.solvers import MonteCarloSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(ising(size=8, temperature=1.0))
>>> stepper = MonteCarloSolver(random_state=0).check().prepare(problem)
>>> state = stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 50.0)).state_
>>> bool(abs(state.mean()) > 0.9)
True

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`OperatorSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)`

Bases: BaseSolver[OperatorProblemSpec]

Base class of solvers that advance operator problems with a stepwise QuTiP solver.

The solver is created at the first frame and again whenever the parameters change, from the current state, so a parameter update or a schedule applies from the next frame on. Coefficients that depend on time stay functions of time, and the others are folded into constant matrices.

Parameters:

Name	Type	Description	Default
`atol`	`float`	The absolute tolerance of the integrator.	`1e-12`
`rtol`	`float`	The relative tolerance of the integrator.	`1e-10`
`nsteps`	`int`	The largest number of internal steps between two frames.	`100000`
`max_dimension`	`int`	The largest dimension of the Hilbert space the solver accepts, which compiling checks.	`DEFAULT_MAX_DIMENSION`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Source code in src/skcomplex/engine/solvers/_operator.py

def __init__(
    self,
    atol: float = 1e-12,
    rtol: float = 1e-10,
    nsteps: int = 100_000,
    max_dimension: int = DEFAULT_MAX_DIMENSION,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.atol = atol
    self.rtol = rtol
    self.nsteps = nsteps
    self.max_dimension = max_dimension
    self.random_state = random_state
    self.until = until

`OptimalControlSolver(*, random_state=None, until=math.inf)`

Bases: Solver

Base class of solvers that compute the controls minimising or maximising an objective over a run of a system.

The solver reads the bounded adaptive parameters of the system as its controls. Its objective is the running cost objective integrated over the run plus the terminal cost terminal at the horizon, the time until of the run. Each cost is a function of the namespace s of the states, parameters, and time of the system by dotted name, such as lambda s: s.plant.x**2 + s.plant.u**2, or a number. The solver minimises the objective for the sense 'min' and maximises it for 'max'. It solves the problem over the whole run from the initial state and the initial parameters at the start of the run, and then follows the solution frame by frame. Subclasses take objective, terminal, and sense as arguments, implement _solve, and may implement _prepare.

Attributes:

Name	Type	Description
`control_problem_`		The optimal control problem of the compiled system.
`horizon_`		The end of the run, the horizon of the objective.
`solution_`		The optimal solution of the run.

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`PontryaginSolver(objective=None, *, terminal=0.0, sense='min', nodes=201, relaxation=0.5, tolerance=1e-06, max_iterations=1000, random_state=None, until=math.inf)`

Bases: OptimalControlSolver

Optimal control by Pontryagin's principle, with costates and the stationarity condition derived symbolically.

The solver builds the Hamiltonian H = L + λ · f from the running cost L and the rates f of the states, and derives the costate equations λ' = -∂H/∂x with λ(T) = ∂φ/∂x at the horizon T for the terminal cost φ. The controls solve ∂H/∂u = 0, projected onto their bounds, or take the bound that minimises H where H is linear in them. A forward-backward sweep integrates the states forward and the costates backward with fourth-order Runge-Kutta steps on a grid, and relaxes the controls towards their update until they converge. The weight of the update halves whenever the distance of the controls to their update stops falling, which breaks the cycles of the sweep. An objective to maximise is minimised with its sign changed. The solution holds the costates.

Parameters:

Name	Type	Description	Default
`objective`	`Cost \| None`	The running cost, a function of the namespace `s` of the states, parameters, and time of the system by dotted name, such as `lambda s: s.plant.x**2`, or a number.	`None`
`terminal`	`Cost`	The terminal cost at the horizon, a function of the same namespace or a number.	`0.0`
`sense`	`str`	`'min'` to minimise the objective or `'max'` to maximise it.	`'min'`
`nodes`	`int`	The number of times of the grid, from the start to the horizon.	`201`
`relaxation`	`float`	The initial weight of the update of the controls at every sweep, between 0 and 1.	`0.5`
`tolerance`	`float`	The largest relative change of the states, costates, and controls between sweeps at convergence.	`1e-06`
`max_iterations`	`int`	The largest number of sweeps.	`1000`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`solution_`		The optimal solution of the run.

Raises:

Type	Description
`SolverError`	The solver has no objective, the sweep diverges or does not converge, or the stationarity condition has no unique explicit solution.
`ValidationError`	The sense is unknown, the grid has fewer than 2 nodes, or the relaxation is not in `(0, 1]`.

Examples:

>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import RateLaw
>>> from skcomplex.specs import ParameterSpec, StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import PontryaginSolver
>>> class PlantElement(Element):
...     def __init__(
...         self,
...         label,
...         x: Annotated[float, StateSpec()] = 1.0,
...         u: Annotated[float, ParameterSpec(adaptive=True, bounds=(-10, 10))] = 0.0,
...     ):
...         super().__init__(label)
...         self.x = x
...         self.u = u
...
...     def _get_laws(self, s):
...         return [RateLaw(s.x, s.u)]
>>> solver = PontryaginSolver(objective=lambda s: (s.plant.x**2 + s.plant.u**2) / 2)
>>> simulation = Simulation(System(elements=[PlantElement('plant')]), solver=solver)
>>> simulation = simulation.check().prepare().run(until=1.0)
>>> float(round(simulation.trajectory_.values('plant.u')[0], 3))
-0.762

Source code in src/skcomplex/engine/solvers/_pontryagin.py

def __init__(
    self,
    objective: Cost | None = None,
    *,
    terminal: Cost = 0.0,
    sense: str = 'min',
    nodes: int = 201,
    relaxation: float = 0.5,
    tolerance: float = 1e-6,
    max_iterations: int = 1000,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.objective = objective
    self.terminal = terminal
    self.sense = sense
    self.nodes = nodes
    self.relaxation = relaxation
    self.tolerance = tolerance
    self.max_iterations = max_iterations
    self.random_state = random_state
    self.until = until

`RK4Solver(dt=0.001, *, random_state=None, until=math.inf)`

Bases: Solver

Classical fourth-order Runge-Kutta integrator of ordinary differential equations with a fixed internal step.

With the Numba backend every internal step of a frame runs in compiled code.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest internal time step.	`0.001`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.system import System
>>> from skcomplex.models.biology.common.populations import ExponentialGrowthInteraction, PopulationElement
>>> from skcomplex.engine.solvers import RK4Solver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(
...     System(
...         elements=[PopulationElement('a', n=1)],
...         interactions=[ExponentialGrowthInteraction('g', target='a', r=1)],
...     ),
... )
>>> stepper = RK4Solver(dt=0.001, random_state=0).check().prepare(problem)
>>> step = StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)
>>> round(float(stepper.run(step).state_[0]), 10)
2.7182818285

Source code in src/skcomplex/engine/solvers/_rk4.py

def __init__(
    self,
    dt: float = 1e-3,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`RenewalSolver(*, random_state=None, until=math.inf)`

Bases: Solver

Exact stepper of renewal processes, which jumps from renewal to renewal without a time grid.

Every member renews on its own clock, and a frame holds the states the last renewal at or before its time set.

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import GroupElement
>>> from skcomplex.specs import StateSpec
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import JumpLaw
>>> from skcomplex.core.laws.utils import draw
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.solvers import RenewalSolver
>>> from skcomplex.specs import StepSpec
>>> class CounterElement(Element):
...     def __init__(self, label, n: Annotated[float, StateSpec()] = 0.0):
...         super().__init__(label)
...         self.n = n
...
...     def _get_laws(self, s):
...         return [JumpLaw(waiting=draw('wait', 'expon'), reset={s.n: s.n + 1})]
>>> problem = compile_system(System(elements=[GroupElement('counters', element=CounterElement, size=3)]))
>>> stepper = RenewalSolver(random_state=0).check().prepare(problem)
>>> state = stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 100.0)).state_
>>> bool(np.all(state[:3] > 50))
True

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`RootSolver(*, method='hybr', tolerance=1e-09, random_state=None, until=np.inf)`

Bases: Solver

Solver of the equations of a static system, which it solves once for the one frame of the run.

The compiled solution gives the states when Sage finds a unique explicit one. Otherwise scipy.optimize.root solves the residuals of the equations with their exact Jacobian, from the initial state as the guess. A system without equations solves to its initial state.

Parameters:

Name	Type	Description	Default
`method`	`str`	The method of `scipy.optimize.root`.	`'hybr'`
`tolerance`	`float`	The largest residual of a solution.	`1e-09`
`random_state`	`int \| Generator \| None`	The seed or generator of the solver, which draws nothing.	`None`
`until`	`float`	The end time of the run the solver prepares for, which a static system does not read.	`inf`

Attributes:

Name	Type	Description
`initial_state_`	`Array`	The solved state vector, the one frame of the run.

Examples:

>>> from sage.all__sagemath_symbolics import cos
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import EquationLaw
>>> from skcomplex.core.spaces import StaticTime
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class FixedPointElement(Element):
...     def __init__(self, label, x: Annotated[float, StateSpec()] = 1.0):
...         super().__init__(label)
...         self.x = x
...
...     def _get_laws(self, s):
...         return [EquationLaw(s.x == cos(s.x), target=s.x)]
>>> system = System(elements=[FixedPointElement('a')], time=StaticTime())
>>> trajectory = Simulation(system).check().prepare().run(1.0).trajectory_
>>> round(float(trajectory.values('a.x')[0]), 7)
0.7390851

Source code in src/skcomplex/engine/solvers/_root.py

def __init__(
    self,
    *,
    method: str = 'hybr',
    tolerance: float = 1e-9,
    random_state: int | np.random.Generator | None = None,
    until: float = np.inf,
) -> None:
    super().__init__(random_state=random_state, until=until)
    self.method = method
    self.tolerance = tolerance

`SDESolver(dt=0.001, *, random_state=None, until=math.inf)`

Bases: Solver

Base class of fixed-step schemes of stochastic differential equations, which take one increment per step.

Subclasses set supports and implement _increment. The steps stop at breakpoints, read delays by linear interpolation of the steps, and locate events in them.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest internal time step.	`0.001`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Source code in src/skcomplex/engine/solvers/_sde.py

def __init__(
    self,
    dt: float = 1e-3,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`SSPRKSolver(order=3, cfl=0.4, dt=0.01, *, random_state=None, until=math.inf)`

Bases: Solver

Strong-stability-preserving Runge-Kutta of order 2 or 3, with steps from the CFL condition of conservation laws.

The schemes of Shu and Osher are convex combinations of Euler steps, so they keep the total variation bound of a limited finite-volume scheme. With finite-volume fields the step is the smaller of cfl * h / s, for the cell size h and the largest wave speed s, and the stability limit of the Laplacian rates at the current parameters, and otherwise dt.

Parameters:

Name	Type	Description	Default
`order`	`int`	2 or 3.	`3`
`cfl`	`float`	The Courant number.	`0.4`
`dt`	`float`	The largest step, and the step without conservation laws.	`0.01`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.models.physics.fluids import burgers
>>> from skcomplex.engine.solvers import SSPRKSolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(burgers(64))
>>> stepper = SSPRKSolver(random_state=0).check().prepare(problem)
>>> state = stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_
>>> bool(np.isclose(state.sum(), problem.initial_state.sum()))
True

Source code in src/skcomplex/engine/solvers/_ssprk.py

def __init__(
    self,
    order: int = 3,
    cfl: float = 0.4,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.order = order
    self.cfl = cfl
    self.dt = dt
    self.random_state = random_state
    self.until = until

`ScipyODESolver(method='auto', rtol=1e-06, atol=1e-09, jacobian='auto', *, random_state=None, until=math.inf)`

Bases: Solver

Integrator of ordinary differential equations based on the step solvers of scipy.integrate.

One integrator runs across the frames of a run, and every frame reads its state from the dense output of the step that covers the frame's time, and it restarts at every frame when a condition on time jumps at a time that depends on parameters. With delays, steps are at most the shortest delay, or a hundredth of the frame when the delay is shorter, and a delay shorter than the step reads the extension of the previous step.

The implicit methods BDF, Radau, and LSODA receive the structure of the Jacobian that the compiler derives from the laws, at every start and restart of the integrator. BDF and Radau estimate the Jacobian by finite differences over groups of independent columns and factor it with sparse LU, which makes them the stiff path for finite-difference and finite element fields. An explicit method that spans more than 1000 stable steps of the diffusion of a field in one call warns once.

Parameters:

Name	Type	Description	Default
`method`	`str`	The integration method: `RK45`, `RK23`, or `DOP853` for non-stiff systems, `BDF` or `Radau` for stiff ones, `LSODA`, which starts with Adams methods and switches to BDF methods while the system is stiff, or the default `'auto'`. It chooses `LSODA` when its dense work array fits, for at most 4096 states or a Jacobian with a narrow band, since `LSODA` stores a dense matrix of the size of the state squared. Larger systems get `BDF` with the sparsity pattern of the Jacobian when fields limit the stable step by diffusion, and `RK45` otherwise.	`'auto'`
`rtol`	`float`	The relative tolerance.	`1e-06`
`atol`	`float`	The absolute tolerance.	`1e-09`
`jacobian`	`str`	How implicit methods obtain the Jacobian. `'auto'` passes the sparsity pattern to `BDF` and `Radau` and the half bandwidths of a pattern narrower than a tenth of the size to `LSODA`, and passes nothing when the pattern is denser than one half. `'exact'` passes the exact sparse Jacobian of the laws, which needs laws with a symbolic form. `'dense'` passes nothing, so the Jacobian is estimated densely. Explicit methods ignore it.	`'auto'`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`method_`		The method the stepper integrates with, which `'auto'` resolves for the problem.
`jacobian_pattern_`		The sparsity pattern of the Jacobian that an implicit method uses, or `None`.
`jacobian_options_`		The keyword arguments of the pattern or the half bandwidths that the integrator receives. The exact Jacobian is passed as `jac` at every start of the integrator, with the parameters of the call.

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.system import System
>>> from skcomplex.models.biology.common.populations import ExponentialGrowthInteraction, PopulationElement
>>> from skcomplex.engine.solvers import ScipyODESolver
>>> from skcomplex.specs import StepSpec
>>> problem = compile_system(
...     System(
...         elements=[PopulationElement('a', n=1)],
...         interactions=[ExponentialGrowthInteraction('g', target='a', r=1)],
...     ),
... )
>>> stepper = ScipyODESolver(rtol=1e-8, random_state=0).check().prepare(problem)
>>> x = stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_
>>> round(float(x[0]), 4)
2.7183

Source code in src/skcomplex/engine/solvers/_scipy_ode.py

def __init__(
    self,
    method: str = 'auto',
    rtol: float = 1e-6,
    atol: float = 1e-9,
    jacobian: str = 'auto',
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.method = method
    self.rtol = rtol
    self.atol = atol
    self.jacobian = jacobian
    self.random_state = random_state
    self.until = until

`SemiLagrangianSolver(dt=0.01, *, random_state=None, until=math.inf)`

Bases: Solver

Directional splitting of Cheng and Knorr, which shifts every slice of a field along each axis in turn.

Every flux of the field must be a linear drift a u, with a free of the state and constant along its own axis. A shift along a periodic axis multiplies the Fourier modes of the slices through scipy.fft, which conserves the total, and a shift along an interval evaluates a cubic spline through scipy.ndimage.map_coordinates. The solver is stable at any Courant number. Rate laws advance by one midpoint step between two half transport steps.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest internal time step.	`0.01`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Examples:

>>> import numpy as np
>>> from skcomplex.core.elements import DensityElement, FieldElement
>>> from skcomplex.core.interactions import AdvectionInteraction
>>> from skcomplex.core.laws.utils import coordinate
>>> from skcomplex.core.spaces import LatticeSpace, PhaseSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import SemiLagrangianSolver
>>> from skcomplex.specs import AxisSpec
>>> space = PhaseSpace(LatticeSpace((16,)), [AxisSpec('v', 8, -1.0, 1.0)]).check().prepare()
>>> start = np.outer(2 + np.sin(2 * np.pi * np.arange(16) / 16), np.ones(8))
>>> gas = FieldElement(
...     'gas', element=DensityElement, space=space, values={'f': start}, discretization='finite-volume'
... )
>>> stream = AdvectionInteraction('s', target='gas', state='f', drift=lambda s: [coordinate('v'), 0])
>>> system = System(elements=[gas], interactions=[stream])
>>> solver = SemiLagrangianSolver(dt=1.0)
>>> final = Simulation(system, solver=solver).check().prepare().run(until=2.0).trajectory_.values('gas.f')[-1]
>>> abs(float(final.sum() / start.sum()) - 1.0) < 1e-12
True

Source code in src/skcomplex/engine/solvers/_semi_lagrangian.py

def __init__(
    self,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`Solver(*, random_state=None, until=math.inf)`

Bases: BaseSolver[ProblemSpec]

Base of the solvers of the compiler's problem, which run its events, renewals, delays, and boundary laws.

Every call to run receives the current parameters, so parameter changes apply at the next call, and applies the boundary laws of the problem after every internal step.

Source code in src/skcomplex/engine/solvers/_base.py

def __init__(self, *, random_state: int | np.random.Generator | None = None, until: float = math.inf) -> None:
    self.random_state = random_state
    self.until = until

`SplitSolver(dt=0.01, *, random_state=None, until=math.inf)`

Bases: Solver

Base of steppers that treat the linear part of the rates of fields apart from the nonlinear remainder.

The linear part is held in Fourier blocks for spectral fields, a small matrix per Fourier mode, and in a sparse matrix over the whole state vector for the other fields. Every block is the shape of the field, the state indices of its states, and its matrices.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest internal time step.	`0.01`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`sparse_`		The sparse matrix of the linear terms outside spectral fields, or `None` without them.
`blocks_`		The spectral blocks.

Source code in src/skcomplex/engine/solvers/_split.py

def __init__(
    self,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`StrangSplittingSolver(transport=None, collision=None, dt=0.01, *, random_state=None, until=math.inf)`

Bases: Solver

Strang splitting of a half transport step, a full collision step, and a half transport step.

The transport solver advances the flux differences of the finite-volume fields and the collision solver the source rates, which the problem keeps apart.

Parameters:

Name	Type	Description	Default
`transport`	`Solver \| None`	The solver of the flux differences, or `None` for `RK4Solver`.	`None`
`collision`	`Solver \| None`	The solver of the source rates, or `None` for `RK4Solver`.	`None`
`dt`	`float`	The largest internal time step, the length of one split step.	`0.01`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`transport_`		The prepared solver of the flux differences.
`collision_`		The prepared solver of the source rates.

Examples:

>>> import numpy as np
>>> from skcomplex.core.elements import DensityElement, FieldElement
>>> from skcomplex.core.interactions import AdvectionInteraction, RelaxationInteraction
>>> from skcomplex.core.laws.utils import coordinate, maximum_entropy
>>> from skcomplex.core.spaces import LatticeSpace, PhaseSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import RK4Solver, SemiLagrangianSolver, StrangSplittingSolver
>>> from skcomplex.specs import AxisSpec
>>> space = PhaseSpace(LatticeSpace((8,)), [AxisSpec('v', 16, -6.0, 6.0)]).check().prepare()
>>> start = np.outer(2 + np.sin(2 * np.pi * np.arange(8) / 8), np.exp(-space.coordinates_['v'] ** 2))
>>> gas = FieldElement(
...     'gas', element=DensityElement, space=space, values={'f': start}, discretization='finite-volume'
... )
>>> stream = AdvectionInteraction('s', target='gas', state='f', drift=lambda s: [coordinate('v'), 0])
>>> relax = RelaxationInteraction(
...     'r', target='gas', state='f', goal=lambda s: maximum_entropy(s.target.f, axes=('v',)), time=0.5
... )
>>> system = System(elements=[gas], interactions=[stream, relax])
>>> solver = StrangSplittingSolver(SemiLagrangianSolver(dt=0.1), RK4Solver(dt=0.1), dt=0.1)
>>> final = Simulation(system, solver=solver).check().prepare().run(until=0.2).trajectory_.values('gas.f')[-1]
>>> abs(float(final.sum() / start.sum()) - 1.0) < 1e-12
True

Source code in src/skcomplex/engine/solvers/_strang_splitting.py

def __init__(
    self,
    transport: Solver | None = None,
    collision: Solver | None = None,
    dt: float = 0.01,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.transport = transport
    self.collision = collision
    self.dt = dt
    self.random_state = random_state
    self.until = until

`UnitarySolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)`

Bases: OperatorSolver

Solver of the linear evolution dx/dt = -i H x of a pure state with QuTiP's SESolver.

Parameters:

Name	Type	Description	Default
`atol`	`float`	The absolute tolerance of the integrator.	`1e-12`
`rtol`	`float`	The relative tolerance of the integrator.	`1e-10`
`nsteps`	`int`	The largest number of internal steps between two frames.	`100000`

Examples:

>>> from typing import Annotated
>>> import numpy as np
>>> from sage.all__sagemath_symbolics import CC, SR
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.laws import DissipatorLaw, GeneratorLaw
>>> from skcomplex.core.laws.utils import local_operator
>>> from skcomplex.specs import HilbertSpec, ObservableSpec, StateSpec
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.engine.solvers import UnitarySolver
>>> operators = {'x': (np.array([[0, 1], [1, 0]]), True), 'z': (np.diag([1, -1]), True)}
>>> operators['minus'] = (np.array([[0, 0], [1, 0]]), False)
>>> QUBIT = HilbertSpec(2, operators, basis={'up': 0, 'down': 1})
>>> class TwoLevelElement(Element):
...     def __init__(self, label, psi: Annotated[object, StateSpec(domain=CC**2)] = 'up', hilbert=QUBIT):
...         super().__init__(label)
...         self.psi = psi
...         self.hilbert = hilbert
...
...     def _get_laws(self, s):
...         return [GeneratorLaw(local_operator(s.psi, SR.symbol('x')) / 2)]
>>> z = [ObservableSpec('z', element='q', variable='z')]
>>> system = System(elements=[TwoLevelElement('q')])
>>> from skcomplex.engine.compile import OperatorCompiler
>>> compiler = OperatorCompiler()
>>> simulation = Simulation(system, compiler=compiler, solver=UnitarySolver(), dt=np.pi, observables=z)
>>> simulation = simulation.check().prepare()
>>> float(round(simulation.run(until=np.pi).trajectory_.values('q.z')[-1], 6))
-1.0

Source code in src/skcomplex/engine/solvers/_operator.py

def __init__(
    self,
    atol: float = 1e-12,
    rtol: float = 1e-10,
    nsteps: int = 100_000,
    max_dimension: int = DEFAULT_MAX_DIMENSION,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.atol = atol
    self.rtol = rtol
    self.nsteps = nsteps
    self.max_dimension = max_dimension
    self.random_state = random_state
    self.until = until

`VelocityVerletSolver(dt=0.001, *, random_state=None, until=math.inf)`

Bases: Solver

Velocity Verlet integrator of particle systems, whose positions move with their velocities.

Every step kicks the velocities by half a step, moves the positions, applies the boundary laws, and kicks the velocities again with the accelerations at the new positions, so a step costs one drift evaluation. It conserves energy over long runs of conservative forces. With the Numba backend every internal step of a frame runs in compiled code.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest internal time step.	`0.001`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Examples:

>>> import numpy as np
>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.spaces import ContinuousSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.core.elements import PointElement
>>> from skcomplex.engine.solvers import VelocityVerletSolver
>>> from skcomplex.specs import StepSpec
>>> system = System(elements=[PointElement('ball', x_d1=2.0)], space=ContinuousSpace(n_dim=1))
>>> problem = compile_system(system)
>>> stepper = VelocityVerletSolver(dt=0.1, random_state=0).check().prepare(problem)
>>> stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_.round(10).tolist()
[2.0, 2.0]

Source code in src/skcomplex/engine/solvers/_velocity_verlet.py

def __init__(
    self,
    dt: float = 1e-3,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

`Yoshida4Solver(dt=0.001, *, random_state=None, until=math.inf)`

Bases: Solver

Yoshida integrator of fourth order, three velocity Verlet steps over two blocks of states.

The solver splits the states into two blocks whose rates each read only the other block: a state of order 2 into its position and velocity, and states of order 1 by the sparsity of the Jacobian of the drift, as the coordinates and the momenta of a separable Hamiltonian split. Every step composes three velocity Verlet steps of the lengths w1 h, w0 h, and w1 h, with w1 = 1 / (2 - 2^(1/3)) and w0 = 1 - 2 w1, which kick the second block by half a step, drift the first, and kick the second again. The scheme is symplectic, so the energy of a separable Hamiltonian stays bounded over long runs. With the Numba backend every internal step of a frame runs in compiled code.

Parameters:

Name	Type	Description	Default
`dt`	`float`	The largest internal time step.	`0.001`
`random_state`	`int \| Generator \| None`	The seed or generator of stochastic solvers.	`None`
`until`	`float`	The end time of the run the solver prepares for, infinite for a run that goes on until it is stopped.	`inf`

Attributes:

Name	Type	Description
`first_`		The flat indices of the states of the block that drifts.
`second_`		The flat indices of the states of the block that kicks.

Raises:

Type	Description
`SolverError`	The problem is not an ODE, or its states admit no split into two blocks whose rates read only the other.

Examples:

>>> from skcomplex.engine.compile import compile_system
>>> from skcomplex.core.spaces import ContinuousSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.core.elements import PointElement
>>> from skcomplex.engine.solvers import Yoshida4Solver
>>> from skcomplex.specs import StepSpec
>>> system = System(elements=[PointElement('ball', x_d1=2.0)], space=ContinuousSpace(n_dim=1))
>>> problem = compile_system(system)
>>> stepper = Yoshida4Solver(dt=0.1).check().prepare(problem)
>>> stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_.round(10).tolist()
[2.0, 2.0]

Source code in src/skcomplex/engine/solvers/_yoshida4.py

def __init__(
    self,
    dt: float = 1e-3,
    *,
    random_state: int | np.random.Generator | None = None,
    until: float = math.inf,
) -> None:
    self.dt = dt
    self.random_state = random_state
    self.until = until

Solvers

AgentView(agent, stepper)

alive property

group property

group_label property

id property

location property

location_key property

agents_here(group=None)

create(group=None, location=None, **values)

displacement_to(other)

distance_to(other)

lags(state, window=None)

move_to(location)

neighbor_cells()

neighbors(radius=None, layer=None)

referenced(state)

remove()

BaseSolver(*, random_state=None, until=math.inf)

CaputoSolver(*, random_state=None, until=math.inf)

CoSimulationSolver(solvers=None, sync_dt=None, scheme='gauss-seidel', split=False, *, random_state=None, until=math.inf)

CollocationSolver(objective=None, *, terminal=0.0, sense='min', nodes=101, scheme='hermite-simpson', method='SLSQP', tolerance=1e-10, max_iterations=1000, random_state=None, until=math.inf)

CrankNicolsonSolver(dt=0.01, *, random_state=None, until=math.inf)

ETDRK4Solver(dt=0.01, *, random_state=None, until=math.inf)

EulerMaruyamaSolver(dt=0.001, *, random_state=None, until=math.inf)

ExternalSolver(*, random_state=None, until=math.inf)

FractionalABMSolver(dt=0.01, *, random_state=None, until=math.inf)

FractionalEulerMaruyamaSolver(dt=0.001, horizon=None, max_increments=_MAX_INCREMENTS, *, random_state=None, until=math.inf)

GillespieSolver(*, random_state=None, until=math.inf)

GradientDescentSolver(precision='float64', *, random_state=None, until=math.inf)

GrunwaldLetnikovSolver(dt=0.01, memory=None, *, random_state=None, until=math.inf)

HeunSolver(dt=0.001, *, random_state=None, until=math.inf)

HorizonError

IMEX2Solver(dt=0.01, *, random_state=None, until=math.inf)

IMEXEulerSolver(dt=0.01, *, random_state=None, until=math.inf)

IterationSolver(*, random_state=None, until=math.inf)

JumpTrajectoriesSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)

KineticSolver(n=10000, dt=0.01, *, random_state=None, until=math.inf)

L1Solver(dt=0.01, memory=None, *, random_state=None, until=math.inf)

LammpsFixSpec(name, group, value, mass, seed) dataclass

LammpsSolver(dt=0.005, pairs=None, fixes=None, *, random_state=None, until=math.inf)

LatticeSolver(*, random_state=None, until=math.inf)

LindbladSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)

MesaSolver(*, random_state=None, until=math.inf)

MilsteinSolver(dt=0.001, *, random_state=None, until=math.inf)

ModelView(stepper, owner='', roles=None)

dimensions property

inputs property

parameters property

rng property

round property

storage property

time property

torus property

agents(name, view=None)

composition(name, state)

count(name)

create(group, location=None, **values)

deposit(name, cell, amount)

field(name, cell)

groups(name)

history(name, window=None)

observable(name)

occupancy(name)

parameter(name)

random_empty_cell()

report(name, value)

set_parameter(name, value)

MonteCarloSolver(*, random_state=None, until=math.inf)

OperatorSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)

OptimalControlSolver(*, random_state=None, until=math.inf)

PontryaginSolver(objective=None, *, terminal=0.0, sense='min', nodes=201, relaxation=0.5, tolerance=1e-06, max_iterations=1000, random_state=None, until=math.inf)

RK4Solver(dt=0.001, *, random_state=None, until=math.inf)

RenewalSolver(*, random_state=None, until=math.inf)

RootSolver(*, method='hybr', tolerance=1e-09, random_state=None, until=np.inf)

SDESolver(dt=0.001, *, random_state=None, until=math.inf)

SSPRKSolver(order=3, cfl=0.4, dt=0.01, *, random_state=None, until=math.inf)

ScipyODESolver(method='auto', rtol=1e-06, atol=1e-09, jacobian='auto', *, random_state=None, until=math.inf)

SemiLagrangianSolver(dt=0.01, *, random_state=None, until=math.inf)

Solver(*, random_state=None, until=math.inf)

`AgentView(agent, stepper)`

`alive` `property`

`group` `property`

`group_label` `property`

`id` `property`

`location` `property`

`location_key` `property`

`agents_here(group=None)`

`create(group=None, location=None, **values)`

`displacement_to(other)`

`distance_to(other)`

`lags(state, window=None)`

`move_to(location)`

`neighbor_cells()`

`neighbors(radius=None, layer=None)`

`referenced(state)`

`remove()`

`BaseSolver(*, random_state=None, until=math.inf)`

`CaputoSolver(*, random_state=None, until=math.inf)`

`CoSimulationSolver(solvers=None, sync_dt=None, scheme='gauss-seidel', split=False, *, random_state=None, until=math.inf)`

`CollocationSolver(objective=None, *, terminal=0.0, sense='min', nodes=101, scheme='hermite-simpson', method='SLSQP', tolerance=1e-10, max_iterations=1000, random_state=None, until=math.inf)`

`CrankNicolsonSolver(dt=0.01, *, random_state=None, until=math.inf)`

`ETDRK4Solver(dt=0.01, *, random_state=None, until=math.inf)`

`EulerMaruyamaSolver(dt=0.001, *, random_state=None, until=math.inf)`

`ExternalSolver(*, random_state=None, until=math.inf)`

`FractionalABMSolver(dt=0.01, *, random_state=None, until=math.inf)`

`FractionalEulerMaruyamaSolver(dt=0.001, horizon=None, max_increments=_MAX_INCREMENTS, *, random_state=None, until=math.inf)`

`GillespieSolver(*, random_state=None, until=math.inf)`

`GradientDescentSolver(precision='float64', *, random_state=None, until=math.inf)`

`GrunwaldLetnikovSolver(dt=0.01, memory=None, *, random_state=None, until=math.inf)`

`HeunSolver(dt=0.001, *, random_state=None, until=math.inf)`

`HorizonError`

`IMEX2Solver(dt=0.01, *, random_state=None, until=math.inf)`

`IMEXEulerSolver(dt=0.01, *, random_state=None, until=math.inf)`

`IterationSolver(*, random_state=None, until=math.inf)`

`JumpTrajectoriesSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)`

`KineticSolver(n=10000, dt=0.01, *, random_state=None, until=math.inf)`

`L1Solver(dt=0.01, memory=None, *, random_state=None, until=math.inf)`

`LammpsFixSpec(name, group, value, mass, seed)` `dataclass`

`LammpsSolver(dt=0.005, pairs=None, fixes=None, *, random_state=None, until=math.inf)`

`LatticeSolver(*, random_state=None, until=math.inf)`

`LindbladSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)`

`MesaSolver(*, random_state=None, until=math.inf)`

`MilsteinSolver(dt=0.001, *, random_state=None, until=math.inf)`

`ModelView(stepper, owner='', roles=None)`

`dimensions` `property`

`inputs` `property`

`parameters` `property`

`rng` `property`

`round` `property`

`storage` `property`

`time` `property`

`torus` `property`

`agents(name, view=None)`

`composition(name, state)`

`count(name)`

`create(group, location=None, **values)`

`deposit(name, cell, amount)`

`field(name, cell)`

`groups(name)`

`history(name, window=None)`

`observable(name)`

`occupancy(name)`

`parameter(name)`

`random_empty_cell()`

`report(name, value)`

`set_parameter(name, value)`

`MonteCarloSolver(*, random_state=None, until=math.inf)`

`OperatorSolver(atol=1e-12, rtol=1e-10, nsteps=100000, max_dimension=DEFAULT_MAX_DIMENSION, *, random_state=None, until=math.inf)`

`OptimalControlSolver(*, random_state=None, until=math.inf)`

`PontryaginSolver(objective=None, *, terminal=0.0, sense='min', nodes=201, relaxation=0.5, tolerance=1e-06, max_iterations=1000, random_state=None, until=math.inf)`

`RK4Solver(dt=0.001, *, random_state=None, until=math.inf)`

`RenewalSolver(*, random_state=None, until=math.inf)`

`RootSolver(*, method='hybr', tolerance=1e-09, random_state=None, until=np.inf)`

`SDESolver(dt=0.001, *, random_state=None, until=math.inf)`

`SSPRKSolver(order=3, cfl=0.4, dt=0.01, *, random_state=None, until=math.inf)`

`ScipyODESolver(method='auto', rtol=1e-06, atol=1e-09, jacobian='auto', *, random_state=None, until=math.inf)`

`SemiLagrangianSolver(dt=0.01, *, random_state=None, until=math.inf)`

`Solver(*, random_state=None, until=math.inf)`

`SplitSolver(dt=0.01, *, random_state=None, until=math.inf)`