ODE Fitting with DiffSL¶

dsl_model = """
in_i {r = 1, k = 1 }
u_i { y = 0.1 }
F_i { (r * y) * (1 - (y / k)) }
"""

print("DiffSL Model Definition:")
print("=" * 50)
print(dsl_model)
print("\nExplanation:")
print("  in_i {r = 1, k = 1 }         → Parameters to fit")
print("  u_i { y = 0.1 }     → Initial condition: y(0) = 0.1")
print("  F_i { ... }         → Right-hand side: dy/dt = ...")

fig, ax = ode_fit(
    t_span,
    y_noisy,
    t_span,
    y_data,
    title="Synthetic Logistic Growth Data",
    xlabel="Time",
    ylabel="Population",
    show=False,
)

# Update legend
ax.legend(["Noisy observations", "True solution"])
plt.show()

# Generate predictions with fitted parameters
r_fit, k_fit = result.x
exp_rt_fit = np.exp(r_fit * t_span)
y_fit = y0 * exp_rt_fit / (1 + y0 * (exp_rt_fit - 1) / k_fit)

# Plot
fig, axes = plt.subplots(2, 1, figsize=(10, 10))

# Top: Data and fits
ax1 = axes[0]
ax1.plot(t_span, y_noisy, "o", label="Noisy data", alpha=0.6, markersize=6)
ax1.plot(t_span, y_fit, "-", label="Fitted model", linewidth=2)
ax1.plot(t_span, y_data, "--", label="True solution", linewidth=2, alpha=0.7)
ax1.set_xlabel("Time")
ax1.set_ylabel("Population")
ax1.set_title("Logistic Growth Model Fit")
ax1.legend()
ax1.grid(True, alpha=0.3)

# Bottom: Residuals
ax2 = axes[1]
residuals = y_noisy - y_fit
ax2.plot(t_span, residuals, "o-", alpha=0.6)
ax2.axhline(0, color="r", linestyle="--", linewidth=2, alpha=0.7)
ax2.set_xlabel("Time")
ax2.set_ylabel("Residuals")
ax2.set_title("Fitting Residuals")
ax2.grid(True, alpha=0.3)

# Add statistics
rmse = np.sqrt(np.mean(residuals**2))
ax2.text(
    0.02,
    0.98,
    f"RMSE = {rmse:.4f}",
    transform=ax2.transAxes,
    verticalalignment="top",
    bbox=dict(boxstyle="round", facecolor="wheat", alpha=0.5),
)

plt.tight_layout()
plt.show()

print("\nFit quality:")
print(f"  RMSE: {rmse:.6f}")
print(f"  Max residual: {np.abs(residuals).max():.6f}")
print(f"  Mean residual: {np.mean(residuals):.6f}")

# Try several initial guesses
initial_guesses = [
    [0.5, 0.5],
    [2.0, 2.0],
    [10.0, 10.0],
    [100.0, 100.0],
]

print("\n" + "=" * 70)
print("INITIAL GUESS SENSITIVITY")
print("=" * 70)
print(f"{'Initial [r, k]':<20} {'Final [r, k]':<30} {'Iterations':<12} {'Success'}")
print("-" * 70)

for guess in initial_guesses:
    result = optimiser.run(problem, guess)
    final_str = f"[{result.x[0]:.4f}, {result.x[1]:.4f}]"
    guess_str = f"{guess}"
    print(f"{guess_str:<20} {final_str:<30} {result.iterations:<12} {result.success}")

print(f"\nTrue values: [r, k] = [{r_true}, {k_true}]")

# Example 1: Multiple state variables (Lotka-Volterra)
predator_prey = """
in_i { alpha = 1, beta = 0.1, gamma = 1.5, delta = 0.075 }
u_i {
    prey = 10.0
    predator = 5.0
}
F_i {
    alpha * prey - beta * prey * predator,
    delta * prey * predator - gamma * predator
}
out_i { prey, predator }
"""

print("Example: Predator-Prey Model (Lotka-Volterra)")
print("=" * 50)
print(predator_prey)

# Example 2: With algebraic equations
with_algebra = """
in_i { k1 = 0.5, k2 = 0.3 }
u_i { A = 1.0 }
F_i {
    -k1 * A,
    k1 * A - k2 * B
}
out_i { A, B }
"""

print("\nExample: Sequential Reactions (A → B → C)")
print("=" * 50)
print(with_algebra)