simulator

Hybrid system simulation utilities.

ensure_state(state)

Validate and coerce a state vector into a 1D array.

Source code in src/flowcean/ode/simulator.py

def ensure_state(state: Iterable[float]) -> State:
    """Validate and coerce a state vector into a 1D array."""
    array = np.asarray(state, dtype=float)
    if array.ndim != 1:
        message = "State must be a 1D array."
        raise ValueError(message)
    return array

simulate(system, t_span, x0=None, location0=None, *, input_stream=None, capture_inputs=None, capture_derivatives=False, max_jumps=256, rtol=1e-07, atol=1e-09, max_step=None, dense_output=False, sample_times=None, sample_dt=None)

Simulate a hybrid system and return a trace.

Parameters:

Name	Type	Description	Default
system	HybridSystem	Hybrid system to simulate.	required
t_span	tuple[float, float]	Start and end time for integration.	required
x0	Iterable[float] | None	Optional initial state override.	None
location0	Location | None	Optional initial location override.	None
input_stream	InputStream | None	Optional input stream accessor for callbacks.	None
capture_inputs	bool | None	Input capture mode. If None, capture iff an input stream is provided.	None
capture_derivatives	bool	Whether to re-evaluate Location.dynamics.flow on the returned trace grid and store the sampled derivatives in Trace.dx. This assumes pure flow callbacks under repeated evaluation. Scalar derivative returns are accepted only for single-state systems.	False
max_jumps	int	Maximum number of transitions allowed.	256
rtol	float	Relative tolerance for the solver.	1e-07
atol	float	Absolute tolerance for the solver.	1e-09
max_step	float | None	Optional maximum step size.	None
dense_output	bool	Whether to build a continuous solution per segment.	False
sample_times	Iterable[float] | None	Monotone time grid to sample from the dense solution.	None
sample_dt	float | None	Fixed sampling interval to generate a time grid.	None

Returns:

Name	Type	Description
Trace	Trace	The simulation trace with location labels and events.

Source code in src/flowcean/ode/simulator.py

def simulate(  # noqa: C901, PLR0912, PLR0915
    system: HybridSystem,
    t_span: tuple[float, float],
    x0: Iterable[float] | None = None,
    location0: Location | None = None,
    *,
    input_stream: InputStream | None = None,
    capture_inputs: bool | None = None,
    capture_derivatives: bool = False,
    max_jumps: int = 256,
    rtol: float = 1e-7,
    atol: float = 1e-9,
    max_step: float | None = None,
    dense_output: bool = False,
    sample_times: Iterable[float] | None = None,
    sample_dt: float | None = None,
) -> Trace:
    """Simulate a hybrid system and return a trace.

    Args:
        system: Hybrid system to simulate.
        t_span: Start and end time for integration.
        x0: Optional initial state override.
        location0: Optional initial location override.
        input_stream: Optional input stream accessor for callbacks.
        capture_inputs: Input capture mode. If ``None``, capture iff an
            input stream is provided.
        capture_derivatives: Whether to re-evaluate ``Location.dynamics.flow``
            on the returned trace grid and store the sampled derivatives in
            ``Trace.dx``. This assumes pure flow callbacks under repeated
            evaluation. Scalar derivative returns are accepted only for
            single-state systems.
        max_jumps: Maximum number of transitions allowed.
        rtol: Relative tolerance for the solver.
        atol: Absolute tolerance for the solver.
        max_step: Optional maximum step size.
        dense_output: Whether to build a continuous solution per segment.
        sample_times: Monotone time grid to sample from the dense solution.
        sample_dt: Fixed sampling interval to generate a time grid.

    Returns:
        Trace: The simulation trace with location labels and events.
    """
    start_location = (
        system.initial_location if location0 is None else location0
    )
    if not isinstance(start_location, Location):
        message = "location0 must be a Location."
        raise TypeError(message)
    location_ids = {id(location) for location in system.locations}
    if id(start_location) not in location_ids:
        message = "location0 must be included in system.locations."
        raise ValueError(message)

    state = ensure_state(x0 if x0 is not None else system.initial_state)
    location = start_location
    effective_input_stream = input_stream or _missing_input_stream
    should_capture = _resolve_capture_inputs(
        capture_inputs=capture_inputs,
        input_stream=input_stream,
    )

    t_segments: list[np.ndarray] = []
    x_segments: list[np.ndarray] = []
    location_segments: list[np.ndarray] = []
    sol_segments: list[Callable[[np.ndarray], np.ndarray] | None] = []
    events: list[Event] = []

    t_current = float(t_span[0])
    t_final = float(t_span[1])
    jumps = 0

    time_epsilon = 1e-12
    sample_grid = _prepare_sample_times(t_span, sample_times, sample_dt)
    needs_dense = dense_output or sample_grid is not None

    while t_current < t_final:
        transitions = system.transitions_from(location)
        event_fns = _build_event_functions(
            transitions,
            system.parameters,
            location.parameters,
            effective_input_stream,
        )

        solve_kwargs = {
            "fun": _wrap_flow(
                location,
                system.parameters,
                effective_input_stream,
            ),
            "t_span": (t_current, t_final),
            "y0": state,
            "events": event_fns or None,
            "rtol": rtol,
            "atol": atol,
            "dense_output": needs_dense,
        }
        if max_step is not None:
            solve_kwargs["max_step"] = max_step

        result = solve_ivp(**solve_kwargs)
        if not result.success:
            message = f"ODE integration failed: {result.message}"
            raise RuntimeError(message)

        t_segments.append(result.t)
        x_segments.append(result.y.T)
        location_segments.append(
            np.full(result.t.shape, location, dtype=object),
        )
        sol_segments.append(result.sol)

        if not result.t_events or all(
            len(event_list) == 0 for event_list in result.t_events
        ):
            break

        triggered_index, event_time, event_state = _first_event(
            result.t_events,
            result.y_events,
        )
        is_zero_time_event = event_time - t_current <= time_epsilon
        transition = transitions[triggered_index]
        jumps += 1
        if jumps > max_jumps:
            message = "Maximum number of transitions exceeded."
            raise RuntimeError(message)

        state, event = _apply_transition(
            transition,
            event_time,
            event_state,
            system.parameters,
            effective_input_stream,
        )
        events.append(event)

        location = transition.target
        while True:
            transitions = system.transitions_from(location)
            immediate_index = _immediate_transition_index(
                transitions,
                location,
                event_time,
                state,
                system.parameters,
                effective_input_stream,
                time_epsilon,
            )
            if immediate_index is None:
                break

            transition = transitions[immediate_index]
            jumps += 1
            if jumps > max_jumps:
                message = "Maximum number of transitions exceeded."
                raise RuntimeError(message)

            state, event = _apply_transition(
                transition,
                event_time,
                state,
                system.parameters,
                effective_input_stream,
            )
            events.append(event)

            location = transition.target

        if is_zero_time_event:
            t_current = min(t_final, t_current + time_epsilon)
        else:
            t_current = min(t_final, event_time + time_epsilon)

    if sample_grid is None:
        t_all = _concat_segments(t_segments)
        x_all = _concat_segments(x_segments)
        location_objects = _concat_segments(location_segments)
        location_all = _location_labels(location_objects)
        u_all = None
        dx_all = None
        if should_capture:
            if input_stream is None:
                message = "Internal error: expected input_stream for capture."
                raise RuntimeError(message)
            u_all = _capture_inputs(t_all, input_stream)
        if capture_derivatives:
            dx_all = _capture_derivatives(
                system=system,
                times=t_all,
                states=x_all,
                locations=location_objects,
                input_stream=effective_input_stream,
            )
        return Trace(
            t=t_all,
            x=x_all,
            location=location_all,
            events=tuple(events),
            u=u_all,
            dx=dx_all,
        )

    rolled = _rollout_segments(
        sample_grid,
        t_segments,
        x_segments,
        location_segments,
        sol_segments,
    )
    u_all = None
    dx_all = None
    if should_capture:
        if input_stream is None:
            message = "Internal error: expected input_stream for capture."
            raise RuntimeError(message)
        u_all = _capture_inputs(rolled.t, input_stream)
    if capture_derivatives:
        dx_all = _capture_derivatives(
            system=system,
            times=rolled.eval_t,
            states=rolled.x,
            locations=rolled.location,
            input_stream=effective_input_stream,
        )
    location_all = _location_labels(rolled.location)
    return Trace(
        t=rolled.t,
        x=rolled.x,
        location=location_all,
        events=tuple(events),
        u=u_all,
        dx=dx_all,
    )

generate_traces(system, t_span, initial_states, *, input_stream=None, capture_inputs=None, capture_derivatives=False, max_jumps=256, rtol=1e-07, atol=1e-09, max_step=None, dense_output=False, sample_times=None, sample_dt=None)

Simulate a batch of traces for a set of initial states.

The input stream and capture semantics match :func:simulate, including the requirement that capture_derivatives=True assumes pure flow callbacks under repeated evaluation on the returned trace grid. Scalar derivative returns are accepted only for single-state systems.

Source code in src/flowcean/ode/simulator.py

def generate_traces(
    system: HybridSystem,
    t_span: tuple[float, float],
    initial_states: Iterable[Iterable[float]],
    *,
    input_stream: InputStream | None = None,
    capture_inputs: bool | None = None,
    capture_derivatives: bool = False,
    max_jumps: int = 256,
    rtol: float = 1e-7,
    atol: float = 1e-9,
    max_step: float | None = None,
    dense_output: bool = False,
    sample_times: Iterable[float] | None = None,
    sample_dt: float | None = None,
) -> list[Trace]:
    """Simulate a batch of traces for a set of initial states.

    The input stream and capture semantics match :func:`simulate`, including
    the requirement that ``capture_derivatives=True`` assumes pure flow
    callbacks under repeated evaluation on the returned trace grid. Scalar
    derivative returns are accepted only for single-state systems.
    """
    return [
        simulate(
            system,
            t_span,
            x0=state,
            input_stream=input_stream,
            capture_inputs=capture_inputs,
            capture_derivatives=capture_derivatives,
            max_jumps=max_jumps,
            rtol=rtol,
            atol=atol,
            max_step=max_step,
            dense_output=dense_output,
            sample_times=sample_times,
            sample_dt=sample_dt,
        )
        for state in initial_states
    ]