Analiz Örnekleri

Genel Bakış

examples/analysis/ klasörü, dynlib'in analiz araçlarının nasıl kullanıldığını gösteren örnekleri içerir. Burada da bu örneklerden bazıları verilmektedir. Kullanılan analiz araçlarının detaylı kullanımı için Analiz rehberine bakınız.

Çekim havzası haritalama ve sınıflandırma

Henon haritası havza keşfi

Yerleşik Henon haritası üzerinde 512×512'lik bir ızgarada, PCR-BM algoritması ile basin_auto analiz aracı ile Henon haritasının çekerleri ve bu çekerlerin çekim havzaları numerik olarak hesaplanmaktadır. basin_auto çeker yapılarını otomatik olarak bulmaya çalışır.

"""
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!")

Bilinen çekerlerin havzası

Bu örnekte ise yine 512×512'lik bir ızgarada yerleşik Henon haritası modeli için basin_known analiz aracı kullanımı gösterilmektedir. Bu aracın kullanılması için çeker yapıları önceden ReferenceRun veya FixedPoint ile tanımlanmalıdır.

"""
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!")

Limit döngü (limit-cycle) havzaları

Bu örnekte basin_auto kullanarak bir limit döngüye ait çekim havzası hesaplanmaktadır. basin_auto limit döngüyü otomatik olarak tespit etmektedir.

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()

Bu örnekte ise yine aynı limit döngü bu sefer ReferenceRun ile tanımlandıktan sonra basin_known aracı ile çekim havzası hesaplanmaktadır.

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()

Lyapunov üstelleri ve spektrumu

Lojistik harita Lyapunov üsteli hesabı

Bu örnekte observer kullanarak maksimum Lyapunov üsteli (MLE) ve Lyapunov spektrumu hesabının nasıl yapılacağı gösterilmektedir. Bir boyutlu lojistik harita için spektrum analizi mantıksız olsa da lyapunov_mle_observer ve lyapunov_spectrum_observer araçlarının numerik olarak ne kadar farklı olduğuna bakılmaktadır.

"""
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.")

Lorenz sistemi Lyapunov üsteli hesabı

Bu örnekte ise Lorenz sistemi için lyapunov_mle_observer ve lyapunov_spectrum_observer kullanarak Lyapunov üsteli hesaplanmaktadır.

"""
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()

Parametre taraması ile Lyapunov üsteli hesaplama

sweep.lyapunov_mle_sweep ile seçilen bir parametre değeri için belli bir aralıkta maksimum Lyapunov üsteli hesabı yapılabilmektedir. Bu örnekte lojistik harita için r ∈ [2.5, 4] aralığında maksimum Lyapunov üsteli hesaplanmakta ve sonuç çizdirilmektedir.

"""
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}")

sweep.lyapunov_spectrum_sweep ise seçilen bir parametre değeri için belli bir aralıkta tüm Lyapunov spektrumunu hesaplamaktadır. Bu örnekte Lorenz sisteminin rho parametresi için [0, 200] aralığında Lyapunov spektrumu hesabı yapılmaktadır.

"""
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()

Manifoldlar ve yörünge bulucular

Manifold izleme

Bu örnekte trace_manifold_1d_map aracı kullanılarak Henon haritasının bir sabit noktası için hem kararlı hem de kararsız 1-boyutlu manifold hesaplanmaktadır. Sonuç plot.manifold ile çizdirilmektedir.

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()

Bu örnekte ise trace_manifold_1d_ode aracı kullanılarak 2-boyutlu bir ODE sisteminin eyer noktasına ait kararlı ve kararsız manifoldlar belirlenerek sonuç çizdirilmektedir.

"""
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()

Homoklinik ve heteroklinik yörünge bulucular

Bu örnekte 2-boyutlu bir ODE sistemi için bir parametre ve denge noktası seçilmekte ve -1.0 değerinden başlayarak [-0.6, -1.2] aralığında seçilen parametre değiştirilerek denge noktası için bir homoklinik yörünge aranmaktadır. Arama için homoclinic_finder kullanılmaktadır. Sonrasında homoklinik yörünge tespit edilen parametre değeri için homoclinic_tracer ve plot.manifold kullanılarak homoklinik yörünge çizdirilmektedir.

"""
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.")

Bu örnekte ise 2-boyutlu başka bir ODE sistemi için iki tane denge noktası seçilmekte ve yine bir parametre için belli bir aralıkta bu iki nokta arasında oluşacak bir heteroklinik yönge aranmaktadır. Arama işlemi heteroclinic_finder aracı ile yapılmaktadır. Sonrasında heteroklinik yörünge tespit edilen parametre değeri için heteroclinic_tracer ve plot.manifold kullanılarak heteroklinik yörünge çizdirilmektedir.

"""
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.")

Parametre taraması ve yörünge analizi

Parametre taramaları

Bu örnekte 1-boyutlu basit bir ODE modelinin a parametresi için farklı değerler seçilerek analysis.sweep.traj_sweep aracı ile her bir parametre için otomatik olarak bir yörünge (trajectory) hesaplanmaktadır. Hesaplanan yörüngeler plot.series.multi ile tek seferde çizdirilmektedir.

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}")

Yörünge analizi

Bu örnekte bir simülasyon yapıldıktan sonra elde edilen sonuçlar üzerinde nasıl analiz yapılabileceği gösterilmektedir. res simülasyon sonucu olmak üzere res.analyze() ile summary(), argmax(), zero_crossings() ve time_above/below gibi post-analiz opsiyonları kullanılabilmektedir.

"""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()