Skip to content

Analysis examples

Overview

The examples/analysis/ folder showcases how to apply dynlib's diagnostic and geometric toolkits. You will find scripts that (1) map basins of attraction with persistence or targeted classifiers, (2) compute Lyapunov exponents and spectra, (3) trace stable/unstable manifolds and homoclinic/heteroclinic orbits, and (4) sweep parameters while atomically analyzing every recorded trajectory.

Basin mapping and classification

Henon map basin discovery

"""
Auto basin of attraction calculation for the Henon Map.

Demonstrates the basin_auto function (Persistent Cell-Recurrence Basin Mapping algorithm)
to identify basins of attraction for the 2D Henon map.

basin_auto tries to automatically discover attractors and their basins. Its success highly
depends on the choice of parameters and the nature of the dynamical system. 
"""

from __future__ import annotations
import numpy as np
from dynlib import setup
from dynlib.analysis import basin_auto, print_basin_summary
from dynlib.plot import export, theme, fig, basin_plot
from dynlib.utils import Timer

# Setup Henon map model
print("Setting up Henon map model...")
sim = setup("builtin://map/henon", stepper="map", jit=True)

# Set standard Henon parameters
a_param = 1.4
b_param = 0.3
sim.assign(a=a_param, b=b_param)

print(f"  Parameters: a={a_param}, b={b_param}")
print(f"  State variables: x, y")

# Domain for basin calculation
# Focusing on a region near the Henon attractor
x_min, x_max = -2.5, 2.5
y_min, y_max = -2.5, 2.5

# Grid resolution (increase for higher detail)
grid_nx = 512
grid_ny = 512

print(f"  Grid size: {grid_nx} × {grid_ny} = {grid_nx * grid_ny} points")
print(f"  Domain: x ∈ [{x_min}, {x_max}], y ∈ [{y_min}, {y_max}]")

# Compute basin of attraction
print("\nComputing basin of attraction...")
print("  Algorithm: PCR-BM (Persistent Cell-Recurrence Basin Mapping)")
print(f"  Grid size: {grid_nx} × {grid_ny} = {grid_nx * grid_ny} points")
print(f"  Domain: x ∈ [{x_min}, {x_max}], y ∈ [{y_min}, {y_max}]")

with Timer("Basin computation time"):
    result = basin_auto(
        sim,
        ic_grid=[grid_nx, grid_ny],
        ic_bounds=[(x_min, x_max), (y_min, y_max)],
        observe_vars=["x", "y"],  # Use both state variables
        grid_res=[128, 128],  # Cell resolution for recurrence detection
        merge_downsample=4,
        post_detect_samples=512,
        max_samples=2000,  # Number of iterations per initial condition
        transient_samples=500,  # Skip initial transient
        window=128,  # Recurrence detection window size
        u_th=0.95,  # Uniqueness threshold
        recur_windows=3,  # Required recurrence windows
        s_merge=0.3,  # Similarity threshold for merging attractors
        p_in=30,  # Persistence threshold
        b_max=100.0,  # Blowup threshold (lower to catch diverging trajectories)
        outside_limit=50,  # Maximum consecutive outside steps
        parallel_mode="auto",
    )

print("Done!")

# Analyze results
print_basin_summary(result)

# Visualization
print("\nCreating visualization...")
theme.use("notebook")
theme.update(grid=False)

# Create figure
ax = fig.single(size=(10, 8))

basin_plot(
    result,
    ax=ax,
    xlabel="x",
    ylabel="y",
    ylabel_rot=0,
    ypad=15,
    title=f"Basin of Attraction - Henon Map (a={a_param}, b={b_param})",
    titlepad=15,
)

export.show()
print("\nVisualization complete!")
Runs basin_auto on the built-in Henon map with a 512×512 grid, PCR-BM detection across both x/y, and timer instrumentation. The script prints attractor metadata, summarizes the result, and then hands the BasinResult to basin_plot so you can see how the automatic fingerprinting divides the plane.

Basin of known attractors

"""
Basin of Attraction for Henon Map using Known Attractors.

Demonstrates the simplified basin_known function to identify basins for the
2D Henon map. basin_known requires reference(s) of the attractor(s) to be 
provided beforehand.
"""

from __future__ import annotations

import numpy as np

from dynlib import setup
from dynlib.analysis import basin_known, print_basin_summary, ReferenceRun
from dynlib.plot import export, theme, fig, basin_plot
from dynlib.utils import Timer

# Setup Henon map model
print("Setting up Henon map model...")
sim = setup("builtin://map/henon", stepper="map", jit=True)

# Set standard Henon parameters
a_param = 1.4
b_param = 0.3
sim.assign(a=a_param, b=b_param)

print(f"  Parameters: a={a_param}, b={b_param}")
print(f"  State variables: x, y")

# Define grid of initial conditions
print("\nGenerating initial condition grid...")

# Domain for basin calculation (region near Henon attractor)
x_min, x_max = -2.5, 2.5
y_min, y_max = -2.5, 2.5

# Grid resolution
grid_nx = 512
grid_ny = 512

print(f"  Grid size: {grid_nx} × {grid_ny} = {grid_nx * grid_ny} points")
print(f"  Domain: x ∈ [{x_min}, {x_max}], y ∈ [{y_min}, {y_max}]")

# Compute basin of known attractors
print("\nComputing basin of known attractors...")

with Timer("basin_known computation"):
    result = basin_known(
        sim,
        attractors=[
            ReferenceRun(name="Henon attractor", ic=[0.1, 0.1]),
        ],
        ic_grid=[grid_nx, grid_ny],
        ic_bounds=[(x_min, x_max), (y_min, y_max)],
        max_samples=500,
        transient_samples=200,
        signature_samples=500,
        tolerance=0.05,      # 5% of attractor range
        min_match_ratio=0.8,  # 80% of points must match
        escape_bounds=[(-5.0, 5.0), (-5.0, 5.0)],  # Wide bounds for escape detection
        b_max=1e6, # Blowup threshold / None means literal NaN/Inf
        parallel_mode="auto",
        refine=True, # refine is faster for high resolution grids
    )

print("Done!")

# Analyze results
print_basin_summary(result)

# Visualization
print("\nCreating visualization...")
theme.use("notebook")
theme.update(grid=False)

ax = fig.single(size=(10, 8))

basin_plot(
    result,
    ax=ax,
    xlabel="x",
    ylabel="y",
    ylabel_rot=0,
    ypad=15,
    title=f"Basin of Attraction - Henon Map (a={a_param}, b={b_param})",
    titlepad=15,
)

export.show()

print("\nVisualization complete!")
Demonstrates basin_known + ReferenceRun to classify the Henon attractor during refinement. It feeds pre-computed attractor fingerprints into the known-attractor library, runs on a 512×512 grid, and shows how signature_samples, tolerance, and min_match_ratio influence the assignment.

Limit-cycle basins

from dynlib import setup
from dynlib.analysis import basin_auto
from dynlib.plot import export, basin_plot, fig, phase

# Energy Template Oscillator (ETO) with Circular L Curve
sim = setup("builtin://ode/eto-circular", jit=True, disk_cache=True, stepper="rk4")
# Globally stable limit cycle parameters
sim.assign(mu=0.8, a=2.0)

results = basin_auto(
    sim,
    ic_grid=[128, 128],
    ic_bounds=[(-3, 3), (-3, 3)],
    dt_obs=0.01,
    # run long enough to converge + see the cycle
    max_samples=4000,          # 40 time units
    transient_samples=2000,    # 20 time units of transient
    # recurrence / evidence
    window=512,                # ~5.12 time units > 1 period
    recur_windows=3,
    u_th=0.8,                  # less strict; avoids missing clean periodic motion
    post_detect_samples=1024,  # add extra evidence (aim: >= 1–2 periods total)
    # merging
    merge_downsample=2,        # coarser fingerprint grid => more overlap
    s_merge=0.4,               # easier merge for same attractor under phase shifts
    # optional: make persistence assignment more stable
    p_in=12,
    online_max_cells=8192,     # avoid truncating the stored cell set for a large cycle
)

# Create an attractor
sim.run(T=300, transient=100)
attr = sim.results()

# Plotting
ax = fig.single()
phase.xy(x=attr["x"], y=attr["y"], ax=ax)
basin_plot(results, ax=ax)

print("Globally stable limit cycle.")
print("Edges can go outside for high mu values but they should converge back.")

export.show()


from dynlib import setup
from dynlib.analysis import basin_known, ReferenceRun
from dynlib.plot import export, basin_plot, fig, phase

# Energy Template Oscillator (ETO) with Circular L Curve
sim = setup("builtin://ode/eto-circular", jit=True, disk_cache=True, stepper="rk4")
# Globally stable limit cycle parameters
sim.assign(mu=0.8, a=2.0)

results = basin_known(
    sim,
    attractors=[
        ReferenceRun(name="Limit Cycle Attractor", ic=[1.0, 0.0]),
    ],
    ic_grid=[128, 128],
    ic_bounds=[(-3, 3), (-3, 3)],
    dt_obs=0.01,
    # run long enough to converge + see the cycle
    signature_samples=8000, # samples are in steps not time
    max_samples=4000,
    transient_samples=2000, 
    tolerance=0.05,      # 5% of attractor range
    min_match_ratio=0.8,  # 80% of points must match
    escape_bounds=[(-5.0, 5.0), (-5.0, 5.0)],  # Wide bounds for escape detection
    b_max=1e6, # Blowup threshold / None means literal NaN/Inf
)

# Create an attractor
sim.run(T=300, transient=100)
attr = sim.results()

# Plotting
ax = fig.single()
phase.xy(x=attr["x"], y=attr["y"], ax=ax)
basin_plot(results, ax=ax)

export.show()
The first script uses basin_auto on the Energy Template Oscillator to capture a globally stable limit cycle via recurrence detection and post-detection persistence. The second repeats the experiment with basin_known/ReferenceRun, showcasing how to compare phase-space points with a reference trajectory while still refining the grid.

Refinement benchmark

"""
Benchmark for basin_known with refine=False and refine=True.

Tests performance on different grid sizes: 64x64, 128x128, 256x256, 512x512.
Uses the same parameters and model as basin_henon_known.py.
"""

from __future__ import annotations

import time
import numpy as np

from dynlib import setup
from dynlib.analysis import basin_known, ReferenceRun

# Setup Henon map model
print("Setting up Henon map model...")
sim = setup("builtin://map/henon", stepper="map", jit=True)

# Set standard Henon parameters
a_param = 1.4
b_param = 0.3
sim.assign(a=a_param, b=b_param)

print(f"  Parameters: a={a_param}, b={b_param}")
print(f"  State variables: x, y")

# Domain for basin calculation (region near Henon attractor)
x_min, x_max = -2.5, 2.5
y_min, y_max = -2.5, 2.5

# Common parameters for basin_known
attractors = [
    ReferenceRun(name="Henon attractor", ic=[0.1, 0.1]),
]
ic_bounds = [(x_min, x_max), (y_min, y_max)]
max_samples = 500
transient_samples = 200
signature_samples = 500
tolerance = 0.05      # 5% of attractor range
min_match_ratio = 0.8  # 80% of points must match
escape_bounds = [(-5.0, 5.0), (-5.0, 5.0)]  # Wide bounds for escape detection
b_max = 1e6  # Blowup threshold / None means literal NaN/Inf
parallel_mode = "auto"

# Grid sizes to test
grid_sizes = [64, 128, 256, 512, 1024]

print("\nBenchmarking basin_known with different grid sizes and refine options...")
grid_width = 12
time_width = 12
match_width = 10
print(f"{'Grid Size':<{grid_width}} | {'refine = False':>{time_width}} | {'refine = True':>{time_width}} | {'Match %':>{match_width}}")
print("-" * (grid_width + time_width * 2 + match_width + 9))

# Warm-up
basin_known(
    sim,
    attractors=attractors,
    ic_grid=[grid_sizes[0], grid_sizes[0]],
    ic_bounds=ic_bounds,
    max_samples=max_samples,
    transient_samples=transient_samples,
    signature_samples=signature_samples,
    tolerance=tolerance,
    min_match_ratio=min_match_ratio,
    escape_bounds=escape_bounds,
    b_max=b_max,
    parallel_mode=parallel_mode,
    refine=False,
)

for size in grid_sizes:
    ic_grid = [size, size]

    # Measure refine=False
    start_time = time.perf_counter()
    result_false = basin_known(
        sim,
        attractors=attractors,
        ic_grid=ic_grid,
        ic_bounds=ic_bounds,
        max_samples=max_samples,
        transient_samples=transient_samples,
        signature_samples=signature_samples,
        tolerance=tolerance,
        min_match_ratio=min_match_ratio,
        escape_bounds=escape_bounds,
        b_max=b_max,
        parallel_mode=parallel_mode,
        refine=False,
    )
    false_time = time.perf_counter() - start_time

    # Measure refine=True
    start_time = time.perf_counter()
    result_true = basin_known(
        sim,
        attractors=attractors,
        ic_grid=ic_grid,
        ic_bounds=ic_bounds,
        max_samples=max_samples,
        transient_samples=transient_samples,
        signature_samples=signature_samples,
        tolerance=tolerance,
        min_match_ratio=min_match_ratio,
        escape_bounds=escape_bounds,
        b_max=b_max,
        parallel_mode=parallel_mode,
        refine=True,
    )
    true_time = time.perf_counter() - start_time

    # Compute match percentage
    match_percentage = np.sum(result_false.labels == result_true.labels) / result_false.labels.size * 100

    print(f"{f'{size}x{size}':<{grid_width}} | {false_time:>{time_width}.2f} | {true_time:>{time_width}.2f} | {match_percentage:>{match_width}.1f}")

print("\nBenchmark complete!")
Bi-directional benchmark toggles refine on basin_known across grid sizes from 64×64 up to 1024×1024, printing wall-clock time and label agreement so you can see the performance payoff before enabling coarse-to-fine passes.

Lyapunov exponents and spectral sweeps

Logistic-map Lyapunov demo

"""
Lyapunov exponent calculation for the logistic map.

Demonstrates the runtime analysis system for computing maximum Lyapunov exponents (MLE)
and Lyapunov spectrum in discrete dynamical systems using the high-level Sim.run() API.

The Lyapunov exponent λ characterizes divergence of nearby trajectories:
    - λ > 0: Chaotic behavior 
    - λ = 0: Periodic orbits
    - λ < 0: Stable fixed points

For the logistic map at r=4: λ = ln(2) ≈ 0.6931
"""

from __future__ import annotations
import numpy as np
from dynlib import setup
from dynlib.runtime.observers import lyapunov_mle_observer, lyapunov_spectrum_observer
from dynlib.plot import series, export, theme, fig


# Single run with multiple observers
sim = setup("builtin://map/logistic", jit=True, disk_cache=True)
model = sim.model

record_every = 1

# Run simulation using the high-level Sim API with multiple observers
print("\nComputing Lyapunov exponents with Sim.run()...")
print(f"  Parameter: r = 4.0")
print(f"  Iterations: 5000")
print(f"  Initial condition: x = 0.4")
print(f"  Observers: MLE and Spectrum")

sim.assign(x=0.4, r=4.0)
sim.run(
    N=5000,
    dt=1.0,
    record_interval=record_every,
    observers=[
        lyapunov_mle_observer(model=sim.model, record_interval=record_every),
        lyapunov_spectrum_observer(model=sim.model, k=1, record_interval=record_every),
    ],
)
result = sim.results()

# Extract runtime observer results
print("\n" + "="*60)
print("RESULTS")
print("="*60)

# Get observer result with ergonomic named access
lyap = result.observers["lyapunov_mle"]
spectrum = result.observers["lyapunov_spectrum"]

# Direct access to final MLE (auto-computed from trace)
mle = lyap.mle  # Final converged value from trace
log_growth = lyap.log_growth
n_steps = int(lyap.steps)

# Get spectrum results (lyapunov_spectrum stores the final value in the trace)
spectrum_steps = int(spectrum.steps)

theoretical_mle = np.log(2.0)
x_trajectory = result["x"]
lyap_trace = lyap["mle"]  # Full trace array (use bracket for arrays)
spectrum_trace = spectrum["lyap0"]  # Spectrum trace for first exponent

# Note: trace may have one less point than trajectory if recording starts at t0
n_points = min(len(x_trajectory), len(lyap_trace), len(spectrum_trace))
iterations = np.arange(n_points)

# Get final spectrum value from trace (after data is loaded)
# spectrum trace has alternating zeros, so find the last non-zero value
nonzero_indices = np.nonzero(spectrum_trace)[0]
spectrum_mle = spectrum_trace[nonzero_indices[-1]] if len(nonzero_indices) > 0 else 0

print("\nMaximum Lyapunov Exponent (MLE):")
print(f"  Computed MLE:      {mle:.10f}")
print(f"  Theoretical MLE:   {theoretical_mle:.10f} (ln(2))")
print(f"  Relative error:    {abs(mle - theoretical_mle)/theoretical_mle * 100:.4f}%")
print(f"  Total iterations:  {n_steps}")

print("\nLyapunov Spectrum:")
print(f"  Spectrum λ₀:       {spectrum_mle:.10f}")
print(f"  Theoretical λ₀:    {theoretical_mle:.10f} (ln(2))")
print(f"  Relative error:    {abs(spectrum_mle - theoretical_mle)/theoretical_mle * 100:.4f}%")
print(f"  Total iterations:  {spectrum_steps}")

# Visualize
theme.use("notebook")
theme.update(grid=True)

fig_obj = fig.grid(rows=3, cols=1, size=(10, 12))

# Plot trajectory (first 500 iterations)
series.plot(
    x=iterations[:500],
    y=x_trajectory[:500],
    style="line",
    ax=fig_obj[0, 0],
    xlabel="Iteration (n)",
    ylabel="$x_n$",
    title=f"Logistic Map Trajectory (r=4.0)",
    lw=1.0,
    color="C0",
)

# Plot Lyapunov exponent convergence
series.plot(
    x=iterations,
    y=lyap_trace[:n_points],
    style="line",
    ax=fig_obj[1, 0],
    xlabel="Iteration (n)",
    ylabel="$\\lambda$ (MLE)",
    title="Lyapunov Exponent Convergence",
    lw=1.5,
    color="C1",
    hlines=[(theoretical_mle, f'Theoretical: λ = ln(2) ≈ {theoretical_mle:.4f}')],
    hlines_kwargs={"color": "gray", "linestyle": "--", "linewidth": 2, "label_pad": 0.10},
)

# Annotate final value
fig_obj[1, 0].text(
    0.98, 0.05,
    f'Computed: λ = {mle:.6f}',
    transform=fig_obj[1, 0].transAxes,
    fontsize=11,
    verticalalignment='bottom',
    horizontalalignment='right',
    bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)
)

# Plot Lyapunov spectrum convergence
series.plot(
    x=iterations,
    y=spectrum_trace[:n_points],
    style="line",
    ax=fig_obj[2, 0],
    xlabel="Iteration (n)",
    ylabel="$\\lambda$ (Spectrum)",
    title="Lyapunov Spectrum Convergence",
    lw=1.5,
    color="C2",
    hlines=[(theoretical_mle, f'Theoretical: λ = ln(2) ≈ {theoretical_mle:.4f}')],
    hlines_kwargs={"color": "red", "linestyle": "--", "linewidth": 2},
)

# Annotate final value from spectrum
spectrum_final = spectrum_trace[n_points-1] if n_points > 0 else 0
fig_obj[2, 0].text(
    0.98, 0.05,
    f'Computed: λ = {spectrum_final:.6f}',
    transform=fig_obj[2, 0].transAxes,
    fontsize=11,
    verticalalignment='bottom',
    horizontalalignment='right',
    bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5)
)

export.show()

print("\n" + "="*60)
print("Parameter Scan")
print("="*60)
print(f"{'r':>6s} {'λ (MLE)':>12s} {'Regime':>15s}")
print("-"*60)

test_r_values = [2.5, 3.2, 3.5, 3.83, 4.0]
scan_sim = setup("builtin://map/logistic", stepper="map", jit=False)

for r_test in test_r_values:
    params_test = np.array([r_test], dtype=model.dtype)
    scan_sim.run(
        N=2000,
        dt=1.0,
        record=False,
        record_interval=1,
        max_steps=2000,
        cap_rec=1,
        cap_evt=1,
        ic=np.array([0.4], dtype=model.dtype),
        params=params_test,
        observers=lyapunov_mle_observer,  # Factory mode: Sim injects model
    )

    lyap_result = scan_sim.results().observers["lyapunov_mle"]
    mle_test = lyap_result.mle  # Direct final value access

    if mle_test < -0.01:
        regime = "Stable"
    elif abs(mle_test) < 0.01:
        regime = "Periodic"
    else:
        regime = "Chaotic"

    print(f"{r_test:6.2f} {mle_test:12.6f} {regime:>15s}")
    scan_sim.reset()

print("="*60)
print("\nThe logistic map exhibits the period-doubling route to chaos.")
print("Chaos emerges around r ≈ 3.57, fully developed at r = 4.")
Runs the logistic map at r=4, attaches both the MLE and spectrum runtime observers, and plots the trajectory, Lyapunov convergence, and final spectrum trace while printing the error against ln(2). The script also prints a quick scan over several r values to classify stability vs. chaos.

Lorenz-system Lyapunov demo

"""
Lyapunov exponent calculation for the Lorenz system.

Demonstrates the runtime analysis system for computing maximum Lyapunov exponents (MLE)
and Lyapunov spectrum in continuous dynamical systems using the high-level Sim.run() API.
"""

from __future__ import annotations
from dynlib import setup
from dynlib.runtime.observers import lyapunov_mle_observer, lyapunov_spectrum_observer
from dynlib.plot import fig, export, phase, series

# 1. Setup simulation
# -------------------
sim = setup("builtin://ode/lorenz", stepper="rk4", jit=True)

# Parameters for standard chaotic regime
sigma, rho, beta = 10.0, 28.0, 8.0/3.0
initial_state = {"x": 1.0, "y": 1.0, "z": 1.0}

# Simulation control
dt = 0.001
total_time = 400
transient = 20
record_interval = 4

print(f"Running Lorenz simulation (T={total_time}, dt={dt})...")

# 2. Run MLE Analysis
# -------------------
print("Computing Maximum Lyapunov Exponent (MLE)...")
sim.assign(**initial_state, sigma=sigma, rho=rho, beta=beta)
sim.run(
    transient=transient,
    T=total_time,
    dt=dt,
    record_interval=record_interval,
    observers=lyapunov_mle_observer(record_interval=record_interval)
)
res_mle = sim.results()
mle_analysis = res_mle.observers["lyapunov_mle"]

# 3. Run Spectrum Analysis
# ------------------------
print("Computing Lyapunov Spectrum (3 exponents)...")
sim.reset()
sim.assign(**initial_state, sigma=sigma, rho=rho, beta=beta)
sim.run(
    transient=transient,
    T=total_time,
    dt=dt,
    record_interval=record_interval,
    observers=lyapunov_spectrum_observer(k=3, record_interval=record_interval)
)
res_spec = sim.results()
spec_analysis = res_spec.observers["lyapunov_spectrum"]

# 4. Results & Validation
# -----------------------
print("\n" + "="*60)
print("RESULTS & VALIDATION")
print("="*60)

# Theoretical values for Lorenz (standard parameters)
ref_lambda1 = 0.9056
ref_lambda2 = 0.0000
ref_lambda3 = -14.5723

# MLE Results
mle_calc = mle_analysis.mle
mle_err = abs(mle_calc - ref_lambda1) / abs(ref_lambda1) * 100
print(f"\nMaximum Lyapunov Exponent (MLE):")
print(f"  Calculated:  {mle_calc:.4f}")
print(f"  Theoretical: {ref_lambda1:.4f}")
print(f"  Error:       {mle_err:.2f}%")

# Spectrum Results
# Get final values from traces (using attribute access for final value)
l1_calc = spec_analysis.lyap0
l2_calc = spec_analysis.lyap1
l3_calc = spec_analysis.lyap2

# Calculate errors (handle division by zero for lambda2)
l1_err = abs(l1_calc - ref_lambda1) / abs(ref_lambda1) * 100
l3_err = abs(l3_calc - ref_lambda3) / abs(ref_lambda3) * 100

print(f"\nLyapunov Spectrum:")
print(f"  λ₁: {l1_calc:.4f} (Ref: {ref_lambda1:.4f}, Err: {l1_err:.2f}%)")
print(f"  λ₂: {l2_calc:.4f} (Ref: {ref_lambda2:.4f}, Err: N/A)")
print(f"  λ₃: {l3_calc:.4f} (Ref: {ref_lambda3:.4f}, Err: {l3_err:.2f}%)")
print(f"  Sum: {l1_calc + l2_calc + l3_calc:.4f} (Ref: {ref_lambda1 + ref_lambda2 + ref_lambda3:.4f})")

# 5. Visualization
# ----------------
print("\nPlotting results...")

# Create a 3-row, 1-column grid
grid = fig.grid(rows=3, cols=1, size=(8, 12), title="Lorenz System Lyapunov Analysis")
ax_phase, ax_mle, ax_spec = grid[0][0], grid[1][0], grid[2][0]

# Plot 1: Phase Portrait (x vs z)
# Use the trajectory from the spectrum run (it has the same parameters)
x = res_spec["x"]
z = res_spec["z"]
phase.xy(
    x=x,
    y=z,
    ax=ax_phase,
    lw=0.5,
    alpha=0.8,
    title="Phase Portrait (Lorenz Attractor)",
    xlabel="x",
    ylabel="z",
)

# Plot 2: MLE Convergence
# Use explicit trace_time from analysis result
t_mle = mle_analysis.trace_time
mle_trace = mle_analysis["mle"]
series.plot(
    x=t_mle,
    y=mle_trace,
    ax=ax_mle,
    label=f"MLE (λ₁) ≈ {mle_analysis.mle:.3f}",
    title="Maximum Lyapunov Exponent Convergence",
    ylabel="MLE",
    legend=True,
    hlines=[(ref_lambda1, f"Theoretical λ₁ ≈ {ref_lambda1:.3f}")],
    hlines_kwargs={"alpha": 0.5},
)
ax_mle.grid(True, alpha=0.3)

# Plot 3: Lyapunov Spectrum Convergence
# Use explicit trace_time from analysis result
t_spec = spec_analysis.trace_time
# Spectrum trace columns: lyap0, lyap1, lyap2
l1 = spec_analysis["lyap0"]
l2 = spec_analysis["lyap1"]
l3 = spec_analysis["lyap2"]

series.multi(
    x=t_spec,
    y=[l1, l2, l3],
    names=[f"λ₁ ≈ {l1[-1]:.3f}", f"λ₂ ≈ {l2[-1]:.3f}", f"λ₃ ≈ {l3[-1]:.3f}"],
    ax=ax_spec,
    title="Lyapunov Spectrum Convergence",
    xlabel="Time",
    ylabel="Exponents",
    legend=True,
    hlines=[ref_lambda1, ref_lambda2, ref_lambda3],
    hlines_kwargs={"alpha": 0.3},
)

ax_spec.grid(True, alpha=0.3)

# Save figure
# export.savefig(grid, "lyapunov_lorenz_demo.png")
export.show()
Repeats the approach for the continuous Lorenz attractor. It computes an rk4 run with lyapunov_mle_observer and another with lyapunov_spectrum_observer, prints the error to the reference (0.9056, 0.0, -14.57) spectrum, and draws the attractor plus convergence traces.

Parameter sweeps with Lyapunov statistics

"""
Lyapunov exponent parameter sweep for the logistic map.

Demonstrates sweep.lyapunov_mle_sweep for computing maximum Lyapunov exponents (MLE)
across a range of parameter values. This reveals the transition from order to chaos
as a continuous function of the control parameter.

For the logistic map x_{n+1} = r*x_n*(1-x_n):
    - r < 3.0: Stable fixed point (λ < 0)
    - 3.0 < r < ~3.57: Periodic orbits (λ = 0)
    - r > ~3.57: Chaotic regime (λ > 0)
    - r = 4.0: Fully chaotic (λ = ln(2) ≈ 0.6931)
"""

from __future__ import annotations
import numpy as np
from dynlib import setup
from dynlib.analysis import sweep
from dynlib.plot import series, export, theme, fig, bifurcation_diagram


# Setup simulation
sim = setup("builtin://map/logistic", jit=True, disk_cache=False)

# Parameter sweep configuration (1000 points won't trigger parallelisation in auto mode)
r_values = np.linspace(2.5, 4.0, 1000)  
record_every = 1

print(f"Computing bifurcation diagram and Lyapunov exponents for {len(r_values)} parameter values...")
print(f"  Parameter range: r ∈ [{r_values[0]:.2f}, {r_values[-1]:.2f}]")
print(f"  Iterations per run: 3000")
print(f"  Transient: 500")
print(f"  Trace recording interval: {record_every}")

# First, compute bifurcation diagram
print("\nComputing bifurcation diagram...")
sweep_bif = sweep.traj_sweep(
    sim,
    param="r",
    values=r_values,
    record_vars=["x"],
    N=200,
    transient=300,
)
result_bif = sweep_bif.bifurcation("x")

# Then, run parameter sweep with MLE analysis
print("\nComputing Lyapunov exponents...")
sim.assign(x=0.4)  # Initial condition
res = sweep.lyapunov_mle_sweep(
    sim,
    param="r",
    values=r_values,
    N=3000,
    transient=500,
    dt=1.0,
    record_interval=record_every,
    parallel_mode="auto",  # Enable parallelization
)

print(f"\nSweep completed!")
print(f"  Sweep kind: {res.kind}")
print(f"  MLE range: [{res.mle.min():.4f}, {res.mle.max():.4f}]")
print(f"  Chaos onset (λ=0): r ≈ {r_values[np.argmin(np.abs(res.mle))]:.4f}")
print(f"  MLE at r=4.0: {res.mle[-1]:.4f} (theoretical: {np.log(2):.4f})")

# Optional: inspect recorded trace per-parameter (list of arrays)
trace_runs = res.traces.get("mle")
if trace_runs:
    print(f"  Recorded {len(trace_runs)} convergence traces (first length={trace_runs[0].shape[0]})")

# ===== Visualization =====
theme.use("notebook")
theme.update(grid=True, vline_label_placement_pad=0.16)

# Create figure with 2 rows: bifurcation diagram (top) and MLE (bottom)
ax = fig.grid(rows=2, cols=1, size=(12, 10))

# Top panel: Bifurcation diagram
bifurcation_diagram(
    result_bif,
    alpha=0.8,
    color="black",
    ax=ax[0, 0],
    xlim=(2.5, 4.0),
    ylim=(0, 1),
    ylabel="x*",
    title="Bifurcation Diagram and Maximum Lyapunov Exponent",
    title_fs=14,
    ylabel_fs=12,
    vlines=[(3.0, 'r=3\n(period-2)'), (3.57, 'r≈3.57\n(chaos)')],
    vlines_kwargs={'color': 'darkred', 
                   'linestyle': '--', 
                   'alpha': 0.3, 
                   'linewidth': 1,
                   'label_position': 'bottom',},
)

# Bottom panel: MLE vs parameter (sharing same x-axis)
series.plot(
    x=res.values,
    y=res.mle,
    style="continuous",
    ax=ax[1, 0],
    xlabel="r",
    ylabel="λ (MLE)",
    lw=1.5,
    color="darkred",
    xlabel_fs=12,
    ylabel_fs=12,
    vlines=[3.0, 3.57],
    vlines_kwargs={'color': 'orange', 'linestyle': '--', 'alpha': 0.3, 'linewidth': 1},
)

# Add horizontal line at λ=0 (chaos boundary)
ax[1, 0].axhline(0, color="gray", ls="--", lw=1, alpha=0.7)
ax[1, 0].text(2.6, 0.05, "λ = 0 (chaos boundary)", fontsize=10, color="gray")

# Set x-limits to match
ax[1, 0].set_xlim(2.5, 4.0)

export.show()

# ===== Analysis Summary =====
print("\n" + "="*70)
print("ANALYSIS SUMMARY")
print("="*70)

# Find chaos onset (where λ crosses zero from negative to positive)
zero_crossings = np.where(np.diff(np.sign(res.mle)) > 0)[0]
if len(zero_crossings) > 0:
    chaos_onset_idx = zero_crossings[0]
    print(f"Chaos onset (first λ=0 crossing): r ≈ {res.values[chaos_onset_idx]:.4f}")

# Count negative, zero, and positive MLE regions
n_negative = np.sum(res.mle < -0.01)
n_zero = np.sum(np.abs(res.mle) <= 0.01)
n_positive = np.sum(res.mle > 0.01)

print(f"\nRegion distribution:")
print(f"  Stable (λ < 0):   {n_negative:3d} points ({n_negative/len(res.mle)*100:.1f}%)")
print(f"  Neutral (λ ≈ 0):  {n_zero:3d} points ({n_zero/len(res.mle)*100:.1f}%)")
print(f"  Chaotic (λ > 0):  {n_positive:3d} points ({n_positive/len(res.mle)*100:.1f}%)")

print(f"\nExtreme values:")
print(f"  Min MLE: {res.mle.min():.6f} at r={res.values[np.argmin(res.mle)]:.4f}")
print(f"  Max MLE: {res.mle.max():.6f} at r={res.values[np.argmax(res.mle)]:.4f}")
uses sweep.lyapunov_mle_sweep for the logistic map, pairing the bifurcation extraction with an MLE sweep over 1000 r values so you can see chaos onset from λ crossing zero.

"""
Lyapunov spectrum parameter sweep for the Lorenz system.

Demonstrates sweep.lyapunov_spectrum_sweep for computing the full Lyapunov spectrum
across a range of parameters, and overlays the results with a Lorenz
bifurcation diagram.

For the Lorenz system (sigma=10, beta=8/3):
    - rho < ~24.74: stable fixed points
    - rho > ~24.74: chaotic regime
"""

from __future__ import annotations

import numpy as np
from dynlib import setup
from dynlib.analysis import sweep
from dynlib.plot import series, export, theme, fig, bifurcation_diagram


# Setup simulation
sim = setup("builtin://ode/lorenz", stepper="rk4", jit=True)

sigma, beta = 10.0, 8.0/3.0
initial_state = {"x": 1.0, "y": 1.0, "z": 1.0}

# Sweep configuration
rho_values = np.linspace(0.0, 200.0, 4000)
dt = 0.01
total_time = 50.0
transient = 50.0
record_interval = 1

print(f"Computing Lorenz bifurcation + spectrum for {len(rho_values)} values...")
print(f"  rho range: [{rho_values[0]:.2f}, {rho_values[-1]:.2f}]")
print(f"  T={total_time}, dt={dt}, transient={transient}")

# Bifurcation diagram (z vs rho)
print("\nComputing bifurcation diagram...")
sim.assign(**initial_state, sigma=sigma, rho=rho_values[0], beta=beta)
bif_sweep = sweep.traj_sweep(
    sim,
    param="rho",
    values=rho_values,
    record_vars=["z"],
    T=total_time,
    dt=dt,
    transient=transient,
    record_interval=record_interval,
)
# Local extrema give clean bifurcation diagrams for continuous systems
bif_result = bif_sweep.bifurcation("z").extrema(
    max_points=50,  # Limit points in chaotic regime
    min_peak_distance=8
)

# Lyapunov spectrum sweep
print("\nComputing Lyapunov spectrum...")
sim.reset()
sim.assign(**initial_state, sigma=sigma, rho=rho_values[0], beta=beta)
res = sweep.lyapunov_spectrum_sweep(
    sim,
    param="rho",
    values=rho_values,
    k=3,
    T=total_time,
    dt=dt,
    transient=transient,
    record_interval=record_interval,
    parallel_mode="auto",
)

print("\nSweep completed!")
print(
    "  Spectrum min/max (lambda1): "
    f"[{res.lyap0.min():.4f}, {res.lyap0.max():.4f}]"
)

# ===== Visualization =====
theme.use("notebook")
theme.update(grid=True)

ax = fig.grid(rows=2, cols=1, size=(12, 10), sharex=True)

# Upper panel: Bifurcation diagram
bifurcation_diagram(
    bif_result,
    alpha=0.6,
    color="black",
    ax=ax[0, 0],
    xlabel="rho",
    ylabel="z",
    title="Lorenz Bifurcation Diagram and Lyapunov Spectrum",
    title_fs=14,
    ylabel_fs=12,
)

# Lower panel: Lyapunov spectrum vs parameter
series.multi(
    x=res.values,
    y=[res.lyap0, res.lyap1, res.lyap2],
    names=["$\\lambda_1$", "$\\lambda_2$", "$\\lambda_3$"],
    ax=ax[1, 0],
    xlabel="rho",
    ylabel="Lyapunov exponents",
    legend=True,
    hlines=[0.0],
    hlines_kwargs={"color": "gray", "ls": "--", "lw": 1, "alpha": 0.7},
    xlim=(rho_values[0], rho_values[-1]),
)

export.show()
sweeps the Lorenz rho parameter with sweep.lyapunov_spectrum_sweep, overlays the z-plane bifurcation, and plots all three Lyapunov exponents across the range.

Manifolds and orbit finders

Manifold tracing

from dynlib import setup
from dynlib.analysis import trace_manifold_1d_map
from dynlib.plot import manifold, theme, export, fig


a, b = 1.4, 0.3
bounds = ((-2.5, 2.5), (-2.5, 2.5))

sim = setup("builtin://map/henon2", jit=True)
sim.assign(a=1.4, b=0.3)

fp = sim.model.fixed_points(seeds=[[0.8, 0.8]])
stable_results = trace_manifold_1d_map( sim, 
                                        fp=fp.points[0],
                                        bounds=bounds,
                                        kind="stable",
                                        steps=80,
                                        hmax=2e-3,
                                        clip_margin=0.35,
                                        max_points_per_segment=40000,
                                        max_segments=300,
                                    )

unstable_results = trace_manifold_1d_map( sim, 
                                          fp=fp.points[0],
                                          bounds=bounds,
                                          kind="unstable",
                                          steps=80,
                                          hmax=2e-3,
                                          clip_margin=0.35,
                                          max_points_per_segment=120000,
                                          max_segments=400,
                                        )

theme.use("paper")
ax = fig.single(size=(5, 5))
manifold(
    result=stable_results,
    ax=ax,
    color="blue",
    label="stable",
    xlim=bounds[0],
    ylim=bounds[1],
    xlabel="$x$",
    ylabel="$y$",
    title="Henon map manifolds",
    aspect="equal",
    lw=0.7,
)
manifold(result=unstable_results, ax=ax, color="red", label="unstable", lw=0.7)
ax.plot(fp.points[0][0], fp.points[0][1], "k+", ms=12)
export.show()
traces the 1D stable and unstable manifolds of the Henon map fixed point using trace_manifold_1d_map, shows both branches on a single plot, and overlays the fixed point.

"""
Example: Tracing stable and unstable manifolds of a 2D saddle point.

Model: x' = x, y' = -y + x^2
Equilibrium at origin (0, 0) with:
  - Unstable eigenvalue λ=+1 (eigenvector along x-axis)
  - Stable eigenvalue λ=-1 (eigenvector along y-axis)

Known analytic solutions:
  - Stable manifold: x = 0 (the y-axis)
  - Unstable manifold: y = x^2/3
"""
from dynlib import setup
from dynlib.analysis import trace_manifold_1d_ode
from dynlib.plot import manifold, theme, export, fig
import numpy as np

# Define the saddle model inline
model_toml = """
inline:
[model]
type = "ode"
name = "saddle2d"

[sim]
t0 = 0.0
T = 10.0
dt = 0.01

[states]
x = 0.1
y = 0.1

[equations.rhs]
x = "x"
y = "-y + x*x"
"""

# Setup simulation
sim = setup(model_toml, jit=True)
bounds = ((-2.5, 2.5), (-2.5, 2.5))

# Trace stable manifold (backward in time)
stable_result = trace_manifold_1d_ode(
    sim,
    fp={"x": 0.0, "y": 0.0},
    kind="stable",
    bounds=bounds,
    dt=0.01,
    max_time=50.0,
    resample_h=0.01,
)

# Trace unstable manifold (forward in time)
unstable_result = trace_manifold_1d_ode(
    sim,
    fp={"x": 0.0, "y": 0.0},
    kind="unstable",
    bounds=bounds,
    dt=0.01,
    max_time=50.0,
    resample_h=0.01,
)

# Validate against analytic solutions
print(f"Stable eigenvalue: {stable_result.eigenvalue:.4f} (expected: -1)")
print(f"Unstable eigenvalue: {unstable_result.eigenvalue:.4f} (expected: +1)")

# Check stable manifold (should be x=0)
stable_x_err = 0.0
for branch in stable_result.branches[0] + stable_result.branches[1]:
    if branch.size > 0:
        stable_x_err = max(stable_x_err, np.max(np.abs(branch[:, 0])))
print(f"Stable manifold max |x| error: {stable_x_err:.3e}")

# Check unstable manifold (should be y=x^2/3)
unstable_err = 0.0
for branch in unstable_result.branches[0] + unstable_result.branches[1]:
    if branch.shape[0] >= 2:
        x_vals = branch[:, 0]
        y_vals = branch[:, 1]
        y_analytic = x_vals**2 / 3.0
        unstable_err = max(unstable_err, np.max(np.abs(y_vals - y_analytic)))
print(f"Unstable manifold max |y - x²/3| error: {unstable_err:.3e}")

# Plot
theme.use("paper")
ax = fig.single(size=(6, 6))

# Stable manifold (blue)
manifold(
    result=stable_result,
    ax=ax,
    color="C0",
    label="W^s (stable)",
    lw=1.5,
)

# Unstable manifold (red)
manifold(
    result=unstable_result,
    ax=ax,
    color="C3",
    label="W^u (unstable)",
    lw=1.5,
)

# Analytic unstable manifold for comparison (dashed black)
x_analytic = np.linspace(bounds[0][0], bounds[0][1], 500)
y_analytic = x_analytic**2 / 3.0
ax.plot(x_analytic, y_analytic, "k--", lw=1.0, label="y = x²/3 (analytic)")

# Mark the equilibrium
ax.plot(0, 0, "ko", ms=8)

ax.set_xlim(bounds[0])
ax.set_ylim(bounds[1])
ax.set_aspect("equal")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_title("Stable/Unstable Manifolds: x'=x, y'=-y+x²")
ax.legend(loc="upper left")

export.show()
builds a simple saddle ODE, traces the manifolds via trace_manifold_1d_ode, compares against analytic formulas (x=0, y=x²/3), and reports eigenvalue errors.

Homoclinic & heteroclinic finders

"""
Example: Homoclinic finder/tracer (manifold finder/tracer) on an inline model.

Workflow:
1) Use homoclinic_finder to locate a parameter value with a homoclinic orbit.
2) Use homoclinic_tracer to trace and plot the orbit at the found parameter.

This example demonstrates the simplified API with:
- List inputs for equilibrium guess (no need for np.array)
- window parameter as list of tuples
- preset parameter for common configurations ("fast", "default", "precise")

Important: Matching Finder and Tracer Settings
-----------------------------------------------
The tracer should use the same section/radius settings as the finder:
- r_sec (section radius) and t_min_event (minimum event time) should match
  the finder's configuration to ensure consistent return-section logic.
- If tracing fails with status 'no_cross', try increasing t_max or loosening
  the finder's gap_tol/x_tol tolerances.

Advanced Configuration
----------------------
For fine-grained control, you can use the full configuration classes:

    from dynlib.analysis import (
        HomoclinicRK45Config,      # RK45 integrator settings (dt0, atol, rtol, etc.)
        HomoclinicBranchConfig,    # Branch tracing settings (eps, t_max, r_sec, etc.)
        HomoclinicFinderConfig,    # Finder config
    )

    # Example: Custom RK45 and branch configuration
    rk_cfg = HomoclinicRK45Config(dt0=1e-4, atol=1e-12, rtol=1e-9, max_steps=5_000_000)
    branch_cfg = HomoclinicBranchConfig(
        t_max=2000.0,
        r_blow=200.0,
        r_sec=1e-2,
        window_min=[-3.0, -3.0],
        window_max=[3.0, 3.0],
        rk=rk_cfg,
    )
    finder_cfg = HomoclinicFinderConfig(
        trace=branch_cfg,
        scan_n=81,
        max_bisect=60,
        gap_tol=1e-6,
        x_tol=1e-4,
    )

    result = homoclinic_finder(sim, ..., cfg=finder_cfg)

You can also use `trace_cfg` for unified branch configuration:

    result = homoclinic_finder(
        sim, ...,
        trace_cfg=branch_cfg,  # applies to the unstable branch
        gap_tol=1e-6,
    )
"""
from dynlib import setup
from dynlib.analysis import (
    homoclinic_finder,
    homoclinic_tracer,
)
from dynlib.plot import theme, fig, export, manifold


MODEL = """
inline:
[model]
type = "ode"
name = "homoclinic-demo"

[sim]
t0 = 0.0
T = 10.0
dt = 0.01

[states]
x = 0.0
y = 0.0

[params]
mu = -0.86

[equations.rhs]
x = "y"
y = "mu*y + x - x*x + x*y"

[equations.jacobian]
expr = [
    ["0", "1"],
    ["1 - 2*x + y", "mu + x"]
]
"""


sim = setup(MODEL, jit=True, disk_cache=False)

# ============================================================================
# SIMPLIFIED API: Using preset and flattened kwargs
# ============================================================================

# Simple case: just specify essential parameters, uses "default" preset
result = homoclinic_finder(
    sim,
    param="mu",
    param_min=-1.2,
    param_max=-0.6,
    param_init=-1.0,
    eq_guess=[0.0, 0.0],  # Lists work! No need for np.array
    # Simplified kwargs:
    preset="default",  # or "fast", "precise"
    window=[(-3, 3), (-3, 3)],  # state-space bounds as (min, max) tuples
    scan_n=61,
    max_bisect=60,
    gap_tol=1e-4,
    x_tol=1e-4,
    t_max=2000.0,
    r_blow=200.0,
)

print("success:", result.success)
print("mu_found:", result.param_found)
print("fail:", result.info["fail"])
print("best_by_abs_gap:", result.info["best_by_abs_gap"])
print("best_by_q:", result.info["best_by_q"])
if result.success and result.miss is not None:
    print(
        "gap_found:",
        result.info["gap_found"],
        "q_found:",
        result.info["q_found"],
        "t_cross:",
        result.info["t_cross"],
    )

if (not result.success) or result.miss is None:
    raise SystemExit("Homoclinic finder did not converge; see info above.")

# Trace the orbit using the simplified API
trace = homoclinic_tracer(
    sim,
    param="mu",
    param_value=result.param_found,
    eq=result.miss.eq,          # Can also use lists: [0.0, 0.0]
    sign_u=result.miss.sign_u,
    preset="default",  # Use same preset as finder for consistency
    window=[(-3, 3), (-3, 3)],
    t_max=2000.0,
    r_blow=200.0,
    r_sec=result.miss.r_sec,
    t_min_event=result.miss.t_min,
)

print("trace success: ", trace.success)
print("status: ", trace.meta.status)
print("event: ", None if trace.meta.event is None else trace.meta.event.kind)

if not trace.success:
    print(f"\nTracing failed with status '{trace.meta.status}'.")
    print("Consider increasing t_max or loosening gap_tol/x_tol.")

# Only plot if we have valid trajectory data
traj = trace.branch_pos[0] if trace.branch_pos and len(trace.branch_pos) > 0 else None
if traj is not None and len(traj) > 1:
    theme.use("paper")
    ax = fig.single()

    title_suffix = "" if trace.success else " (partial, no return-section hit)"

    manifold(
        result=trace,
        ax=ax,
        color="C0",
        lw=1.2,
        label="Homoclinic excursion" + (" (closed)" if trace.success else " (incomplete)"),
        xlim=(-1.6, 1.6),
        ylim=(-1.0, 1.0),
        xlabel="x",
        ylabel="y",
        title=f"Homoclinic orbit at mu={result.param_found:.10f}{title_suffix}",
        aspect="equal",
    )
    ax.scatter([result.miss.eq[0]], [result.miss.eq[1]], s=60, label="Equilibrium")
    ax.scatter([trace.meta.x_cross[0]], [trace.meta.x_cross[1]], s=30, label="Section cross")
    ax.legend(loc="best")

    """
    Section cross: This is the point where the unstable manifold trajectory
    re-enters the return section defined by r_sec around the equilibrium.
    When a homoclinic orbit exists, the crossing should align with the unstable
    direction, making the signed return value g close to zero.
    """

    export.show()
else:
    print("\nNo trajectory data to plot.")


"""
Example: Heteroclinic finder/tracer (manifold finder/tracer) on an inline model.

Workflow:
1) Use heteroclinic_finder to locate a parameter value with a heteroclinic orbit.
2) Use heteroclinic_tracer to trace and plot the orbit at the found parameter.

This example demonstrates the simplified API with:
- List inputs for fixed points (no need for np.array)
- window parameter as list of tuples
- preset parameter for common configurations ("fast", "default", "precise")

Important: Matching Finder and Tracer Settings
-----------------------------------------------
The tracer's hit_radius should be chosen carefully based on the finder's accuracy:
- If finder reports gap_found ~ 1e-4, use hit_radius ~ 0.01 to 0.05
- If finder reports gap_found ~ 1e-6, use hit_radius ~ 0.001 to 0.01
- Rule of thumb: hit_radius ≈ 5-50x the finder's gap_found value
- Use the same preset for both finder and tracer for consistency

The finder and tracer may have slightly different numerical behavior, so the
tracer might not achieve the exact same precision as the finder's gap measurement.
If the tracer fails with status 'window' or gets close but doesn't hit the target,
either increase hit_radius or tighten the finder's gap_tol/x_tol tolerances.

Advanced Configuration
----------------------
For fine-grained control, you can use the full configuration classes:

    from dynlib.analysis import (
        HeteroclinicRK45Config,      # RK45 integrator settings (dt0, atol, rtol, etc.)
        HeteroclinicBranchConfig,    # Branch tracing settings (eps, t_max, window, etc.)
        HeteroclinicFinderConfig2D,  # Finder config for 2D systems
        HeteroclinicFinderConfigND,  # Finder config for N-dimensional systems
    )

    # Example: Custom RK45 and branch configuration
    rk_cfg = HeteroclinicRK45Config(dt0=1e-4, atol=1e-12, rtol=1e-9, max_steps=5_000_000)
    branch_cfg = HeteroclinicBranchConfig(
        t_max=500.0,
        r_blow=500.0,
        window_min=[-20.0, -20.0],
        window_max=[20.0, 20.0],
        rk=rk_cfg,
    )
    finder_cfg = HeteroclinicFinderConfigND(
        trace_u=branch_cfg,
        trace_s=branch_cfg,
        scan_n=121,
        max_bisect=80,
        gap_tol=1e-6,
    )

    result = heteroclinic_finder(sim, ..., cfg=finder_cfg)

You can also use `trace_cfg` for unified branch configuration:

    result = heteroclinic_finder(
        sim, ...,
        trace_cfg=branch_cfg,  # applies to both unstable and stable branches
        gap_tol=1e-6,
    )
"""
import numpy as np

from dynlib import setup
from dynlib.analysis import (
    heteroclinic_finder,
    heteroclinic_tracer,
)
from dynlib.plot import theme, fig, export, manifold


MODEL = """
inline:
[model]
type = "ode"
name = "heteroclinic-demo"

[sim]
t0 = 0.0
T = 10.0
dt = 0.01

[states]
u = 0.0
v = 0.0

[params]
c = 0.0
a = 0.3

[equations.rhs]
u = "v"
v = "-c*v - u*(1-u)*(u-a)"

[equations.jacobian]
expr = [
    ["0", "1"],
    ["3*u*u - 2*(1+a)*u + a", "-c"]
]
"""


sim = setup(MODEL, jit=True, disk_cache=False)

# NOTE: Stricter hit_radius requires stricter config values!
#       With 0.02 you will see a gap near fixed points.
hit_radius = 0.02

# ============================================================================
# SIMPLIFIED API: Using preset and flattened kwargs
# ============================================================================

# Simple case: just specify essential parameters, uses "default" preset
result = heteroclinic_finder(
    sim,
    param="c",
    param_min=-0.5,
    param_max=1.0,
    param_init=0.0,
    source_eq_guess=[0.0, 0.0],  # Lists work! No need for np.array
    target_eq_guess=[1.0, 0.0],
    # Simplified kwargs:
    preset="default",  # or "fast", "precise"
    window=[(-10, 10), (-10, 10)],  # state-space bounds as (min, max) tuples
    scan_n=61,
    max_bisect=60,
    gap_tol=1e-4,
    x_tol=1e-4,
    t_max=200.0,
    r_blow=200.0,
)

print("success:", result.success)
print("c_found:", result.param_found)
print("fail:", result.info["fail"])
print("sign_u, sign_s:", result.info["sign_u"], result.info["sign_s"])
print("best_by_abs_gap:", result.info["best_by_abs_gap"])
print("best_by_q:", result.info["best_by_q"])
if result.success and result.miss is not None:
    print(
        "gap_found:",
        result.info["gap_found"],
        "gap_tol_eff_found:",
        result.info["gap_tol_eff_found"],
        "q_found:",
        result.info["q_found"],
    )

if (not result.success) or result.miss is None:
    raise SystemExit("Heteroclinic finder did not converge; see info above.")

# Trace the orbit using the simplified API
# Note: hit_radius should be adjusted based on the accuracy of the finder's gap (gap_found).
# If the finder's gap is ~1e-4, a hit_radius of 1e-3 may be too strict for the tracer.
# Reasonable values: 0.01 to 0.05 for typical cases, or 5-10x the finder's gap_found.
trace = heteroclinic_tracer(
    sim,
    param="c",
    param_value=result.param_found,
    source_eq=result.miss.source_eq,  # Can also use lists: [0.0, 0.0]
    target_eq=result.miss.target_eq,
    sign_u=result.miss.sign_u,
    preset="default",  # Use same preset as finder for consistency
    window=[(-10, 10), (-10, 10)],
    t_max=200.0,
    hit_radius=hit_radius,  # Adjusted to be more realistic given gap_found ~ 2e-4
)

print("trace success: ", trace.success)
print("status: ", trace.meta.status)
print("event: ", None if trace.meta.event is None else trace.meta.event.kind)

if not trace.success:
    # Diagnose why tracing failed
    if trace.branch_pos and len(trace.branch_pos) > 0:
        traj = trace.branch_pos[0]
        dist_to_B = np.linalg.norm(traj[-1] - result.miss.target_eq)
        print(f"\nTracing failed with status '{trace.meta.status}'.")
        print(f"Trajectory reached distance {dist_to_B:.6f} from target_eq.")
        print(f"This exceeds hit_radius={hit_radius}.")
        print("Consider increasing hit_radius or tightening finder tolerances.")

if not trace.success:
    # Diagnose why tracing failed
    if trace.branch_pos and len(trace.branch_pos) > 0:
        traj = trace.branch_pos[0]
        dist_to_B = np.linalg.norm(traj[-1] - result.miss.target_eq)
        print(f"\nTracing failed with status '{trace.meta.status}'.")
        print(f"Trajectory reached distance {dist_to_B:.6f} from target_eq.")
        print(f"This exceeds hit_radius={hit_radius}.")
        print("Consider increasing hit_radius or tightening finder tolerances.")

# Only plot if we have valid trajectory data
if trace.branch_pos and len(trace.branch_pos) > 0:
    theme.use("paper")
    ax = fig.single()

    # Add title indicating whether full heteroclinic connection was achieved
    title_suffix = "" if trace.success else " (partial, did not reach B within hit_radius)"

    manifold(
        result=trace,
        ax=ax,
        color="C0",
        lw=1.2,
        label="Wu(A)" + (" → B" if trace.success else " (incomplete)"),
        xlim=(-0.5, 1.5),
        ylim=(-0.2, 0.4),
        xlabel="u",
        ylabel="v",
        title=f"Heteroclinic orbit at c={result.param_found:.10f}{title_suffix}",
        aspect="equal",
    )
    ax.scatter([result.miss.source_eq[0]], [result.miss.source_eq[1]], s=60, label="A")
    ax.scatter([result.miss.target_eq[0]], [result.miss.target_eq[1]], s=60, label="B")
    ax.scatter([result.miss.x_u_cross[0]], [result.miss.x_u_cross[1]], s=30, label="cross (Wu)")
    ax.scatter([result.miss.x_s_cross[0]], [result.miss.x_s_cross[1]], s=30, label="cross (Ws)")
    ax.legend(loc="best")

    """
    cross (Wu) (or x_u_cross): This is the point where the unstable manifold (Wu) of
    equilibrium point A crosses a section plane centered at equilibrium point B.

    cross (Ws) (or x_s_cross): This is the point where the stable manifold (Ws) of
    equilibrium point B crosses the same section plane.

    When a heteroclinic orbit exists, these two crossing points should be very close
    (ideally the same point), meaning the unstable manifold of A connects to the stable
    manifold of B. The algorithm searches for parameter values where this "miss distance"
    is minimized.
    """

    export.show()
else:
    print("\nNo trajectory data to plot.")
Each script runs the corresponding finder/tracer pair on an inline 2D model. They demonstrate the simplified API (preset, tuple-style window, optional x_tol/gap_tol), print finder diagnostics (gap size, status), and plot the traced orbit in the return section so you can visualize whether the unstable manifold connects back.

Sweeps and trajectory helpers

Parameter sweeps

import numpy as np
from dynlib import setup
from dynlib.plot import series, export
from dynlib.analysis import sweep


DSL = """
inline:
[model]
type="ode"

[states]
x=10.0

[params]
a=-1.0

[equations.rhs]
x = "a * x"
"""

sim = setup(DSL, 
            jit=True,
            disk_cache=False)

values = np.arange(-5.0, -1.0, 1.0)
res=sweep.traj_sweep(sim, param="a", record_vars=["x"], values=values, T=5)

# Plot all trajectories at once (stacked access)
series.multi(x=res.t, y=res["x"], legend=False)
export.show()

# Access individual runs for custom processing
print(f"\nSweep has {len(res.runs)} runs:")
for run in res.runs:
    final_x = run["x"][-1]
    print(f"  a={run.param_value:.1f}: x(T)={final_x:.6f}")

# Access a specific run
run = res.runs[2]  # Third run (a=-3.0)
print(f"\nRun details for a={run.param_value}:")
print(f"  Time points: {len(run.t)}")
print(f"  Initial x: {run['x'][0]:.6f}")
print(f"  Final x: {run['x'][-1]:.6f}")
Uses analysis.sweep.traj_sweep to run the toy ODE x' = a x for several a values, stacks the trajectories with series.multi, and accesses each SweepRun to print the final state and step details.

Trajectory analyzer

"""Demonstration of TrajectoryAnalyzer helpers on a damped oscillator."""

from __future__ import annotations
import numpy as np
from dynlib import setup


MODEL = """
inline:
[model]
type="ode"
name="Damped oscillator with aux energy"

[states]
x=1.0
v=0.0

[params]
omega=2.0
zeta=0.15

[equations.rhs]
x = "v"
v = "-2*zeta*omega*v - omega**2 * x"

[aux]
energy = "0.5 * (v**2 + (omega * x)**2)"
"""


def _round_dict(values: dict[str, float], digits: int = 3) -> dict[str, float]:
    return {k: round(float(v), digits) for k, v in values.items()}


def main() -> None:
    sim = setup(MODEL, stepper="rk4", jit=False, disk_cache=False)
    sim.config(dt=0.01, max_steps=150_000)
    sim.run(T=15.0, record_vars=["x", "v", "aux.energy"])
    res = sim.results()

    print("\n--- Single-variable analysis (x) ---")
    x = res.analyze("x")
    print("summary:", _round_dict(x.summary()))
    t_peak, peak = x.argmax()
    print(f"peak x at t={t_peak:.3f}: {peak:.3f}")
    crossings = x.zero_crossings(direction="up")
    print("first three zero up-crossings (s):", np.round(crossings[:3], 3))
    print(f"time with x >= 0.5: {x.time_above(0.5):.3f}s")
    print(f"time with x < -0.5: {x.time_below(-0.5):.3f}s")

    print("\n--- Multi-variable analysis (x, v) ---")
    mv = res.analyze(["x", "v"])
    print("argmax by var:", {k: (round(t, 3), round(v, 3)) for k, (t, v) in mv.argmax().items()})
    print("range by var:", _round_dict(mv.range()))
    print("mean by var:", _round_dict(mv.mean()))

    print("\n--- Aux variable analysis (energy) ---")
    energy = res.analyze("energy")  # auto-detected aux without explicit prefix
    print("energy min/max:", (round(energy.min(), 3), round(energy.max(), 3)))
    print("energy range:", round(energy.range(), 3))
    below_times = energy.crossing_times(threshold=0.1, direction="down")
    print("first drop below 0.1:", round(below_times[0], 3) if len(below_times) else "never")

    print("\n--- Automatic variable selection ---")
    auto = res.analyze()
    print("recorded vars picked by res.analyze():", auto.vars)


if __name__ == "__main__":
    main()
Runs a damped oscillator with an auxiliary energy variable, then calls res.analyze() to show summary(), argmax(), zero_crossings(), and time_above/Below for single- and multi-variable selections. It also inspects the automatically selected auxiliary variables to demonstrate how TrajectoryAnalyzer finds recorded names.