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 system.	required

Source code in src/skcomplex/core/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/core/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/core/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/core/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.
    """
    system = SystemView(self._stepper)
    dimensions = system.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 system.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/core/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/core/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, or to a uniform point inside a polygon.

Parameters:

Name	Type	Description	Default
location	object	The cell coordinates, network node, or position, or a polygon in a continuous space.	required

Raises:

Type	Description
SpaceError	The location does not fit the space.

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

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

    Args:
        location:
            The cell coordinates, network node, or position, or a polygon in a continuous space.

    Raises:
        SpaceError:
            The location does not fit the space.
    """
    space = self._stepper.problem_.agent_space
    if isinstance(space, ContinuousSpace):
        self._agent.position = self._stepper._position_of(location, space)
        return
    self._agent.cell = self._stepper._cell_at(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/core/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, group=None)

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

Agents of a variable group placed on the nodes of a network make a set of members per node that grows and shrinks during a run, such as the balls of an urn on every node, and an agent on a node reads the sets of the neighbouring nodes as the neighbours of one group.

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
group	str | None	The path of a group whose agents alone are returned, or None for every group.	None

Returns:

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

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

def neighbors(
    self,
    radius: float | None = None,
    layer: str | None = None,
    group: 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.

    Agents of a variable group placed on the nodes of a network make a set of members per node that grows and
    shrinks during a run, such as the balls of an urn on every node, and an agent on a node reads the sets of the
    neighbouring nodes as the neighbours of one group.

    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.

        group:
            The path of a group whose agents alone are returned, or `None` for every group.

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

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/core/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)

regions_here(layer)

Return the indices of the regions of a layer that contain the agent, on their boundaries too.

Parameters:

Name	Type	Description	Default
layer	str	The name of the layer.	required

Returns:

Name	Type	Description
indices	list[int]	The indices of the regions, in increasing order.

Raises:

Type	Description
SimulationError	The space has no such layer.

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

def regions_here(self, layer: str) -> list[int]:
    """Return the indices of the regions of a layer that contain the agent, on their boundaries too.

    Args:
        layer:
            The name of the layer.

    Returns:
        indices:
            The indices of the regions, in increasing order.

    Raises:
        SimulationError:
            The space has no such layer.
    """
    hits = self._stepper._index_of(layer).query(Point(self.location), predicate='intersects')
    return sorted(int(index) for index in hits)

remove()

Remove the agent from the system, 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/core/engine/solvers/_mesa.py

def remove(self) -> None:
    """Remove the agent from the system, 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/core/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/core/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 discrete partition that holds a linear program takes one step per synchronisation interval, whatever its length, so the program solves at every synchronisation point. 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.components.system import System
>>> from skcomplex.systems.biology.common.populations import ExponentialGrowthInteraction, PopulationElement
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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.prepare().run(until=1)
>>> sorted(simulation.stepper_.steppers_)
['a', 'b']

Source code in src/skcomplex/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import RateLaw
>>> from skcomplex.core.specs import ParameterSpec, StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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.prepare().run(until=1.0)
>>> float(round(simulation.trajectory_.values('plant.u')[0], 3))
-0.762

Source code in src/skcomplex/core/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.core.engine.solvers import CrankNicolsonSolver
>>> CrankNicolsonSolver(dt=0.5).dt
0.5

Source code in src/skcomplex/core/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.core.engine.solvers import ETDRK4Solver
>>> ETDRK4Solver(dt=0.25).dt
0.25

Source code in src/skcomplex/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.system import System
>>> from skcomplex.systems.economics.finance import AssetElement, GeometricBrownianMotionInteraction
>>> from skcomplex.core.engine.solvers import EulerMaruyamaSolver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(
...     elements=[AssetElement('stock', price=100)],
...     interactions=[GeometricBrownianMotionInteraction('gbm', target='stock', mu=0.05, sigma=0.2)],
... )
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = EulerMaruyamaSolver(dt=1e-3, random_state=0).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/core/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

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

Bases: Solver

Explicit Euler integrator of ordinary differential equations with a fixed internal step.

Every step is x + h f(t, x), the integration Vensim and Stella use by default, so a run with the time step of a system file reproduces the run of those tools step by step.

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:

>>> from skcomplex.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.system import System
>>> from skcomplex.systems.biology.common.populations import ExponentialGrowthInteraction, PopulationElement
>>> from skcomplex.core.engine.solvers import EulerSolver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(
...     elements=[PopulationElement('a', n=1)],
...     interactions=[ExponentialGrowthInteraction('g', target='a', r=1)],
... )
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = EulerSolver(dt=0.5, random_state=0).prepare(problem)
>>> step = StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)
>>> float(stepper.run(step).state_[0])
2.25

Source code in src/skcomplex/core/engine/solvers/_euler.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/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import RateLaw
>>> from skcomplex.core.specs import StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.compile import PassCompiler
>>> from skcomplex.core.engine.solvers import FractionalABMSolver
>>> from skcomplex.core.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 = PassCompiler().prepare(System(elements=[RelaxingElement('a')])).run().problem_
>>> stepper = FractionalABMSolver(dt=0.001, random_state=0).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/core/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.components.elements import Element, GroupElement
>>> from skcomplex.core.components.laws import FractionalNoiseLaw
>>> from skcomplex.core.specs import StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.compile import PassCompiler
>>> from skcomplex.core.engine.solvers import FractionalEulerMaruyamaSolver
>>> from skcomplex.core.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)]
>>> system = System(elements=[GroupElement('walkers', element=WalkerElement, size=500)])
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = FractionalEulerMaruyamaSolver(dt=0.01, horizon=1.0, random_state=0).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/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.components.elements import CompartmentElement
>>> from skcomplex.systems.biology.epidemiology import RecoveryInteraction
>>> from skcomplex.core.engine.solvers import GillespieSolver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(
...     elements=[CompartmentElement('I', n=10), CompartmentElement('R', n=0)],
...     interactions=[RecoveryInteraction('rec', infected='I', recovered='R', gamma=1.0, stochastic=True)],
... )
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = GillespieSolver(random_state=0).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/core/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.core.engine.solvers import GradientDescentSolver
>>> GradientDescentSolver(precision='float32').precision
'float32'

Source code in src/skcomplex/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import RateLaw
>>> from skcomplex.core.specs import StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.compile import PassCompiler
>>> from skcomplex.core.engine.solvers import GrunwaldLetnikovSolver
>>> from skcomplex.core.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 = PassCompiler().prepare(System(elements=[RelaxingElement('a')])).run().problem_
>>> stepper = GrunwaldLetnikovSolver(dt=0.001, random_state=0).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/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.systems.economics.finance import geometric_brownian_motion
>>> from skcomplex.core.engine.solvers import HeunSolver
>>> from skcomplex.core.specs import StepSpec
>>> problem = PassCompiler().prepare(geometric_brownian_motion()).run().problem_
>>> stepper = HeunSolver(dt=0.01, random_state=0).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/core/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.core.engine.solvers import IMEX2Solver
>>> IMEX2Solver(dt=0.1).dt
0.1

Source code in src/skcomplex/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.systems.physics.heat import heat_diffusion
>>> from skcomplex.core.engine.solvers import IMEXEulerSolver
>>> from skcomplex.core.specs import StepSpec
>>> problem = PassCompiler().prepare(heat_diffusion((16, 16))).run().problem_
>>> stepper = IMEXEulerSolver(dt=1.0, random_state=0).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/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.system import System
>>> from skcomplex.systems.economics.markets import CobwebInteraction, MarketElement
>>> from skcomplex.core.engine.solvers import IterationSolver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(
...     elements=[MarketElement('m', price=1)],
...     interactions=[CobwebInteraction('c', target='m', a=10, b=2, c=1, d=1)],
... )
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = IterationSolver(random_state=0).prepare(problem)
>>> stepper.run(StepSpec(0.0, problem.initial_state, problem.initial_parameters, 1.0)).state_.tolist()
[4.0]

Source code in src/skcomplex/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import DissipatorLaw, GeneratorLaw
>>> from skcomplex.core.components.laws.utils import local_operator
>>> from skcomplex.core.specs import HilbertSpec, ObservableSpec, StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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.core.engine.compile import OperatorCompiler
>>> simulation = Simulation(
...     system, compiler=OperatorCompiler(), solver=solver, dt=1.0, random_state=0, observables=z,
... ).prepare().run(until=20.0)
>>> float(simulation.trajectory_.values('q.z')[-1])
-1.0

Source code in src/skcomplex/core/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.core.engine.solvers import KineticSolver
>>> from skcomplex.core.components.laws import BilinearLaw, EncounterLaw
>>> KineticSolver(n=1000)._get_support().laws == {EncounterLaw, BilinearLaw}
True

Source code in src/skcomplex/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import RateLaw
>>> from skcomplex.core.specs import StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.compile import PassCompiler
>>> from skcomplex.core.engine.solvers import L1Solver
>>> from skcomplex.core.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 = PassCompiler().prepare(System(elements=[RelaxingElement('a')])).run().problem_
>>> stepper = L1Solver(dt=0.001, random_state=0).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/core/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.components.elements import GroupElement, PointElement
>>> from skcomplex.core.components.interactions import PairCouplingInteraction
>>> from skcomplex.core.components.spaces import ContinuousSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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).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/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.systems.computer_science.automata import elementary_automaton
>>> from skcomplex.core.engine.solvers import LatticeSolver
>>> from skcomplex.core.specs import StepSpec
>>> problem = PassCompiler().prepare(elementary_automaton(30, size=11)).run().problem_
>>> stepper = LatticeSolver(random_state=0).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/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import DissipatorLaw, GeneratorLaw
>>> from skcomplex.core.components.laws.utils import local_operator
>>> from skcomplex.core.specs import HilbertSpec, ObservableSpec, StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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.core.engine.compile import OperatorCompiler
>>> compiler = OperatorCompiler()
>>> simulation = Simulation(system, compiler=compiler, solver=LindbladSolver(), dt=1.0, observables=z)
>>> simulation = simulation.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/core/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. It runs a continuous space as a box of two or three sides.

Attributes:

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

Examples:

>>> from typing import Annotated
>>> from skcomplex.core.components.elements import Element, GroupElement
>>> from skcomplex.core.components.laws import RuleLaw
>>> from skcomplex.core.specs import StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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, system):
...         agent.count += 1
>>> system = System(elements=[GroupElement('counters', element=CounterElement, size=3)])
>>> trajectory = Simulation(system, solver=MesaSolver()).prepare().run(until=4).trajectory_
>>> trajectory.agents('counters', 4).values['count'].tolist()
[4.0, 4.0, 4.0]

Source code in src/skcomplex/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.systems.economics.finance import geometric_brownian_motion
>>> from skcomplex.core.engine.solvers import MilsteinSolver
>>> from skcomplex.core.specs import StepSpec
>>> problem = PassCompiler().prepare(geometric_brownian_motion()).run().problem_
>>> stepper = MilsteinSolver(dt=0.01, random_state=0).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/core/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

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.core.engine.compile import PassCompiler
>>> from skcomplex.systems.physics.statistical_mechanics import ising
>>> from skcomplex.core.engine.solvers import MonteCarloSolver
>>> from skcomplex.core.specs import StepSpec
>>> problem = PassCompiler().prepare(ising(size=8, temperature=1.0)).run().problem_
>>> stepper = MonteCarloSolver(random_state=0).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/core/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/core/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/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import RateLaw
>>> from skcomplex.core.specs import ParameterSpec, StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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.prepare().run(until=1.0)
>>> float(round(simulation.trajectory_.values('plant.u')[0], 3))
-0.762

Source code in src/skcomplex/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.system import System
>>> from skcomplex.systems.biology.common.populations import ExponentialGrowthInteraction, PopulationElement
>>> from skcomplex.core.engine.solvers import RK4Solver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(
...     elements=[PopulationElement('a', n=1)],
...     interactions=[ExponentialGrowthInteraction('g', target='a', r=1)],
... )
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = RK4Solver(dt=0.001, random_state=0).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/core/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.components.elements import GroupElement
>>> from skcomplex.core.specs import StateSpec
>>> from skcomplex.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.elements import Element
>>> from skcomplex.core.components.laws import JumpLaw
>>> from skcomplex.core.components.laws.utils import draw
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.solvers import RenewalSolver
>>> from skcomplex.core.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})]
>>> system = System(elements=[GroupElement('counters', element=CounterElement, size=3)])
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = RenewalSolver(random_state=0).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/core/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

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

Bases: BaseSolver[RewritingProblemSpec]

Solver of rewriting systems, which rewrites every structure in parallel once per generation.

The structures live next to the state vector, whose components are their symbol, prototile, or label counts and sizes. Stochastic productions draw from the simulation's random stream.

Examples:

>>> from skcomplex.core.components.spaces import DiscreteTime
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.components.laws import ProductionLaw
>>> from skcomplex.core.engine.solvers import RewritingSolver
>>> from skcomplex.core.components.elements import WordElement
>>> class AlgaeElement(WordElement):
...     def _get_laws(self, namespace):
...         return [ProductionLaw('A', 'AB'), ProductionLaw('B', 'A')]
>>> system = System(elements=[AlgaeElement('algae', axiom='A')], time=DiscreteTime())
>>> from skcomplex.core.engine.compile import RewritingCompiler
>>> compiler = RewritingCompiler()
>>> simulation = Simulation(system, compiler=compiler, solver=RewritingSolver())
>>> simulation = simulation.prepare()
>>> trajectory = simulation.run(until=5).trajectory_
>>> trajectory.values('algae.size').tolist()
[1.0, 2.0, 3.0, 5.0, 8.0, 13.0]

Source code in src/skcomplex/core/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.components.elements import Element
>>> from skcomplex.core.components.laws import EquationLaw
>>> from skcomplex.core.components.spaces import StaticTime
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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).prepare().run(1.0).trajectory_
>>> round(float(trajectory.values('a.x')[0]), 7)
0.7390851

Source code in src/skcomplex/core/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/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.systems.physics.fluids import burgers
>>> from skcomplex.core.engine.solvers import SSPRKSolver
>>> from skcomplex.core.specs import StepSpec
>>> problem = PassCompiler().prepare(burgers(64)).run().problem_
>>> stepper = SSPRKSolver(random_state=0).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/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.system import System
>>> from skcomplex.systems.biology.common.populations import ExponentialGrowthInteraction, PopulationElement
>>> from skcomplex.core.engine.solvers import ScipyODESolver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(
...     elements=[PopulationElement('a', n=1)],
...     interactions=[ExponentialGrowthInteraction('g', target='a', r=1)],
... )
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = ScipyODESolver(rtol=1e-8, random_state=0).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/core/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.components.elements import DensityElement, FieldElement
>>> from skcomplex.core.components.interactions import AdvectionInteraction
>>> from skcomplex.core.components.laws.utils import coordinate
>>> from skcomplex.core.components.spaces import LatticeSpace, PhaseSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.engine.solvers import SemiLagrangianSolver
>>> from skcomplex.core.specs import AxisSpec
>>> space = PhaseSpace(LatticeSpace((16,)), [AxisSpec('v', 8, -1.0, 1.0)]).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).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/core/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/core/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/core/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.components.elements import DensityElement, FieldElement
>>> from skcomplex.core.components.interactions import AdvectionInteraction, RelaxationInteraction
>>> from skcomplex.core.components.laws.utils import coordinate, maximum_entropy
>>> from skcomplex.core.components.spaces import LatticeSpace, PhaseSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.engine.solvers import RK4Solver, SemiLagrangianSolver, StrangSplittingSolver
>>> from skcomplex.core.specs import AxisSpec
>>> space = PhaseSpace(LatticeSpace((8,)), [AxisSpec('v', 16, -6.0, 6.0)]).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).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/core/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

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

System 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 system.	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/core/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 system, 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/core/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/core/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/core/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, a polygon for a uniform point inside it, 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/core/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, a polygon for a uniform point inside it, 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/core/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/core/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/core/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/core/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/core/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/core/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/core/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 system"} 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/core/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

region(layer, index)

Return the polygon of a region of a layer of the continuous space, which a rule can pass to move_to.

Parameters:

Name	Type	Description	Default
layer	str	The name of the layer.	required
index	int	The index of the region in its layer.	required

Returns:

Name	Type	Description
polygon	Polygon	The polygon of the region.

Raises:

Type	Description
SimulationError	The space has no such layer, or the layer has no region of that index.

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

def region(self, layer: str, index: int) -> Polygon:
    """Return the polygon of a region of a layer of the continuous space, which a rule can pass to `move_to`.

    Args:
        layer:
            The name of the layer.

        index:
            The index of the region in its layer.

    Returns:
        polygon:
            The polygon of the region.

    Raises:
        SimulationError:
            The space has no such layer, or the layer has no region of that index.
    """
    polygons = self._stepper._index_of(layer).geometries
    if not 0 <= index < len(polygons):
        error_msg = f'The layer {layer!r} has {len(polygons)} regions, and a rule asks for the region {index}.'
        raise SimulationError(error_msg)
    return polygons[index]

regions(layer)

Return the number of regions of a layer of the continuous space.

Parameters:

Name	Type	Description	Default
layer	str	The name of the layer.	required

Returns:

Name	Type	Description
count	int	The number of regions of the layer.

Raises:

Type	Description
SimulationError	The space has no such layer.

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

def regions(self, layer: str) -> int:
    """Return the number of regions of a layer of the continuous space.

    Args:
        layer:
            The name of the layer.

    Returns:
        count:
            The number of regions of the layer.

    Raises:
        SimulationError:
            The space has no such layer.
    """
    return len(self._stepper._index_of(layer).geometries)

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/core/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 system"} 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/core/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)

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.components.elements import Element
>>> from skcomplex.core.components.laws import DissipatorLaw, GeneratorLaw
>>> from skcomplex.core.components.laws.utils import local_operator
>>> from skcomplex.core.specs import HilbertSpec, ObservableSpec, StateSpec
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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.core.engine.compile import OperatorCompiler
>>> compiler = OperatorCompiler()
>>> simulation = Simulation(system, compiler=compiler, solver=UnitarySolver(), dt=np.pi, observables=z)
>>> simulation = simulation.prepare()
>>> float(round(simulation.run(until=np.pi).trajectory_.values('q.z')[-1], 6))
-1.0

Source code in src/skcomplex/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.spaces import ContinuousSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.components.elements import PointElement
>>> from skcomplex.core.engine.solvers import VelocityVerletSolver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(elements=[PointElement('ball', x_d1=2.0)], space=ContinuousSpace(n_dim=1))
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = VelocityVerletSolver(dt=0.1, random_state=0).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/core/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.core.engine.compile import PassCompiler
>>> from skcomplex.core.components.spaces import ContinuousSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.components.elements import PointElement
>>> from skcomplex.core.engine.solvers import Yoshida4Solver
>>> from skcomplex.core.specs import StepSpec
>>> system = System(elements=[PointElement('ball', x_d1=2.0)], space=ContinuousSpace(n_dim=1))
>>> problem = PassCompiler().prepare(system).run().problem_
>>> stepper = Yoshida4Solver(dt=0.1).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/core/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