Skip to content

Runtime utilities and diagnostics

Overview

The runtime utilities in this suite exercise dynlib's low-level control hooks: setup() can drive a broad range of steppers, stop conditions, and event observers, while helper APIs reveal what the compiler built or what a trajectory is doing in real time. These scripts demonstrate how to sanity-check numerical accuracy, cancel simulations early, detect interesting transitions, and print the DSL equations so you can audit what is actually running inside Sim.

Example scripts

Stepper accuracy ranking

import numpy as np
import warnings
from dynlib import setup, list_steppers, get_stepper

# Suppress RuntimeWarnings to avoid cluttering output with solver convergence messages
warnings.simplefilter("ignore", RuntimeWarning)

# Define models and their exact solutions
models = [
    {
        "name": "Exponential Decay",
        "model": '''
inline:
[model]
type="ode"

[states]
x=1.0

[params]
a=1.0

[equations]
expr = """
dx = -a * x
"""
''',
        "exact": lambda t: np.exp(-t)
    },
    {
        "name": "Harmonic Oscillator",
        "model": '''
inline:
[model]
type="ode"

[states]
x=1.0
v=0.0

[params]
omega=100.0

[equations]
expr = """
dx = v
dv = -omega**2 * x
"""
''',
        "exact": lambda t: (np.cos(100.0 * t), -100.0 * np.sin(100.0 * t))
    }
]

# Simulation parameters
T = 10  # total time
dt = 1e-4  # time step
N = int(T / dt)  # number of steps

ode_steppers = list_steppers(kind="ode")

for model_info in models:
    model_name = model_info["name"]
    model_str = model_info["model"]
    exact_func = model_info["exact"]

    print(f"\n=== {model_name} ===")

    # Dictionary to store errors
    errors = {}
    failed_steppers = {}

    for name in ode_steppers:
        spec = get_stepper(name)
        try:        
            print(f"Running simulation with stepper: {name}")

            # Setup simulation
            sim = setup(model_str, stepper=name, jit=False)
            sim.config(dt=dt, max_steps=N*10)  # allow more steps if needed

            # Run simulation
            sim.run(T=T)
            res = sim.results()

            # Compute relative error on x
            if model_name == "Exponential Decay":
                x_exact = exact_func(res.t)
                rel_error = np.abs((res["x"] - x_exact) / x_exact)
            else:  # Harmonic Oscillator
                x_exact, v_exact = exact_func(res.t)
                rel_error = np.abs((res["x"] - x_exact) / x_exact)

            rms_rel_error = np.sqrt(np.mean(rel_error**2))

            errors[name] = rms_rel_error
            print(f"RMS relative error for {name}: {rms_rel_error:.2e}")

        except Exception as e:
            print(f"Error with stepper {name}: {e}")
            failed_steppers[name] = str(e)
            continue

    # Sort steppers by accuracy (lowest error first)
    sorted_steppers = sorted(errors.items(), key=lambda x: x[1])

    print(f"\n{model_name} - Stepper accuracy ranking (lowest RMS relative error = most accurate):")
    for rank, (name, error) in enumerate(sorted_steppers, 1):
        print(f"{rank}. {name}: {error:.2e}")

    if failed_steppers:
        print(f"\n{model_name} - Failed steppers:")
        for name, error_msg in failed_steppers.items():
            print(f"- {name}: {error_msg}")
Runs every registered ODE stepper against two analytic solutions (exponential decay and the harmonic oscillator) with a tiny constant dt. Each run computes an RMS relative error so you can sort steppers by accuracy, log failures, and compare how each integrator behaves when your time step is near its stability limit.

Early exit macros

"""
Early Exit Demo
===============

This example demonstrates the early exit feature using the `stop` condition
with built-in DSL macros. The simulation will terminate early when a specified
condition is met.

Available macros for stop conditions:
- cross_up(state, threshold): detects upward threshold crossing
- cross_down(state, threshold): detects downward threshold crossing  
- in_interval(state, lower, upper): checks if state is in range
- enters_interval(state, lower, upper): detects entering a range
- leaves_interval(state, lower, upper): detects leaving a range
- increasing(state): detects if state is increasing
- decreasing(state): detects if state is decreasing
- changed(state): detects any change in state

Note: Event macros that use lagged values (like cross_up) require the
simulation to have at least one prior step, so they won't trigger at t0.
"""

from dynlib import setup


def get_exit_reason(res):
    """Return a user-friendly description of why the simulation ended."""
    if res.exited_early:
        return "Early exit (stop condition met)"
    elif res.ok:
        return "Completed normally (max steps reached)"
    else:
        return "Failed or interrupted"


def demo_simple_threshold():
    """Example 1: Stop using a simple threshold condition"""
    print("\n" + "=" * 70)
    print("Example 1: Simple Threshold (x > 0.8)")
    print("=" * 70)

    model = '''
inline:
[model]
type = "map"
name = "Logistic Map - Simple Threshold"

[states]
x = 0.1

[params]
r = 3.5

[equations.rhs]
x = "r * x * (1 - x)"

[sim]
t0 = 0.0
dt = 1.0
stepper = "map"
record = true
stop = "x > 0.8"
'''

    sim = setup(model, jit=False)
    sim.run(N=100)
    res = sim.results()

    print(f"Exit reason: {get_exit_reason(res)}")
    print(f"Steps executed: {res.step_count_final}")
    print(f"Final state: x = {res['x'][-1]:.6f}")
    print(f"Interpretation: Stopped when x first exceeded 0.8")


def demo_interval_check():
    """Example 2: Stop using in_interval macro"""
    print("\n" + "=" * 70)
    print("Example 2: Interval Check with in_interval()")
    print("=" * 70)

    model = '''
inline:
[model]
type = "map"
name = "Logistic Map - Interval Check"

[states]
x = 0.1

[params]
r = 3.5

[equations.rhs]
x = "r * x * (1 - x)"

[sim]
t0 = 0.0
dt = 1.0
stepper = "map"
record = true
stop = "in_interval(x, 0.82, 0.86)"
'''

    sim = setup(model, jit=False)
    sim.run(N=100)
    res = sim.results()

    print(f"Exit reason: {get_exit_reason(res)}")
    print(f"Steps executed: {res.step_count_final}")
    print(f"Final state: x = {res['x'][-1]:.6f}")
    print(f"Interpretation: Stopped when x entered the interval [0.82, 0.86]")


def demo_cross_up():
    """Example 3: Stop using cross_up macro (requires lag)"""
    print("\n" + "=" * 70)
    print("Example 3: Threshold Crossing with cross_up()")
    print("=" * 70)

    model = '''
inline:
[model]
type = "map"
name = "Logistic Map - Cross Up Detection"

[states]
x = 0.1

[params]
r = 3.5

[equations.rhs]
x = "r * x * (1 - x)"

[sim]
t0 = 0.0
dt = 1.0
stepper = "map"
record = true
stop = "cross_up(x, 0.8)"
'''

    sim = setup(model, jit=False)
    sim.run(N=100)
    res = sim.results()

    print(f"Exit reason: {get_exit_reason(res)}")
    print(f"Steps executed: {res.step_count_final}")
    print(f"Final state: x = {res['x'][-1]:.6f}")
    if res.n >= 2:
        print(f"Previous state: x = {res['x'][-2]:.6f}")
        print(f"Interpretation: Stopped when x crossed 0.8 from below")
        print(f"  (previous: {res['x'][-2]:.6f} <= 0.8, current: {res['x'][-1]:.6f} > 0.8)")


def demo_decreasing():
    """Example 4: Stop when state starts decreasing"""
    print("\n" + "=" * 70)
    print("Example 4: Detect Decreasing Trend with decreasing()")
    print("=" * 70)

    model = '''
inline:
[model]
type = "map"
name = "Logistic Map - Detect Decrease"

[states]
x = 0.1

[params]
r = 3.5

[equations.rhs]
x = "r * x * (1 - x)"

[sim]
t0 = 0.0
dt = 1.0
stepper = "map"
record = true
stop = "decreasing(x)"
'''

    sim = setup(model, jit=False)
    sim.run(N=100)
    res = sim.results()

    print(f"Exit reason: {get_exit_reason(res)}")
    print(f"Steps executed: {res.step_count_final}")
    print(f"Final state: x = {res['x'][-1]:.6f}")
    if res.n >= 2:
        print(f"Previous state: x = {res['x'][-2]:.6f}")
        print(f"Interpretation: Stopped when x started decreasing")
        print(f"  (previous: {res['x'][-2]:.6f}, current: {res['x'][-1]:.6f})")


if __name__ == "__main__":
    print("\n" + "=" * 70)
    print("EARLY EXIT FEATURE DEMONSTRATION")
    print("Using Built-in DSL Macros for Stop Conditions")
    print("=" * 70)

    demo_simple_threshold()
    demo_interval_check()
    demo_cross_up()
    demo_decreasing()

    print("\n" + "=" * 70)
    print("Demo Complete!")
    print("=" * 70)
    print("\nKey Takeaways:")
    print("  • Use res.exited_early to check if stop condition was triggered")
    print("  • res.ok is True for both normal completion and early exit")
    print("  • Simple conditions (x > 0.8) work at any step")
    print("  • in_interval() checks current state without lag")
    print("  • cross_up(), decreasing(), etc. use lag and require prior steps")
    print("  • All conditions are evaluated after each simulation step")
    print("=" * 70 + "\n")
Explores the DSL stop macros such as cross_up, in_interval, and decreasing. For several logistic-map variants the script runs up to 100 steps but halts as soon as the condition becomes true, printing the exit reason, executed steps, and the values just before/after the trigger. It highlights how stopping early keeps the run deterministic while saving time.

Transition detection via events

from dynlib import setup
from dynlib.plot import series, export, theme

"""
This example demonstrate lagging mechanism to detect trasition 
from negative to positive value in one of the state variables 
of a chaotic system.
"""

detect_mod = '''
inline:
[mod]
name = "detect_transition"

[mod.add.events.detect]
cond = "cross_up(x, 0)"
phase = "post"
log = ["t"]
'''

sim = setup("builtin://ode/lorenz", stepper="rk4", jit=False, mods=[detect_mod])
sim.config(dt=0.01)
sim.run(T=50.0, transient=10.0)

res = sim.results()
ev = res.event("detect")

theme.use("paper")

series.plot(x=res.t,
            y=res["x"],
            vlines=ev.t,
            xlabel='Time',
            xlabel_fs=11, # Theme override
            ylabel='x',
            title='Lorenz system: x variable with transition detection')

export.show()
Adds an event that fires when the Lorenz x variable crosses zero from below. The run records the timestamp of each detect event, draws vertical lines on the x(t) series plot, and demonstrates how to inspect res.event("detect") to drive annotations or downstream logic in chaotic regimes.

Printing compiled equations

# examples/print_equations_demo.py
"""
Demonstrate printing model equations from the DSL spec.

Prints the equations for the built-in Henon map and Lorenz system.
"""

from dynlib import setup


def main() -> None:
    print("Henon map equations:")
    sim_henon = setup("builtin://map/henon", stepper="map", jit=False)
    sim_henon.model.print_equations()

    print("\nLorenz system equations:")
    sim_lorenz = setup("builtin://ode/lorenz", stepper="rk4", jit=False)
    sim_lorenz.model.print_equations(tables=("equations", "equations.jacobian"))

if __name__ == "__main__":
    main()
Builds the Henon map and Lorenz system with jit=False and simply prints the RHS and Jacobian tables that the compiler derives from Toml or builtin specs. This is handy when you need to confirm that dynlib sees the model you expect before running expensive simulations.