#!/usr/bin/env catnip
# Intégration d'un système proie-prédateur (Lotka-Volterra) par solve_ivp.
#
# dx/dt = alpha*x - beta*x*y   (proies)
# dy/dt = delta*x*y - gamma*y  (prédateurs)
#
# Catnip porte le modèle, les paramètres et les oracles ; SciPy intègre et
# localise les maxima. Deux invariants indépendants contrôlent la trajectoire :
# la quantité conservée V = delta*x - gamma*ln(x) + beta*y - alpha*ln(y), et la
# régularité des périodes entre maxima de proies. Ces maxima sont trouvés par
# le mécanisme d'événements de solve_ivp — une seconde lambda Catnip, dont le
# solveur cherche les racines — puis triés côté Catnip, faute de pouvoir régler
# leur direction (voir plus bas).
#
# Référence : Lotka (1925), Volterra (1926) — voir Wikipedia « Lotka-Volterra
# equations » pour la dérivation de l'invariant.
#
# DEPS: numpy scipy matplotlib
# OFFLINE: données synthétiques déterministes

numpy = import('numpy')
import('scipy.integrate', 'solve_ivp')
mpl = import('matplotlib')
mpl.use('Agg')
plt = import('matplotlib.pyplot')
import('pathlib', 'Path')

script_dir = Path(META.file).parent
output_dir = script_dir / 'output'
output_dir.mkdir(exist_ok=True)

ALPHA = 1.0    # croissance des proies sans prédateur
BETA = 0.1     # prédation par rencontre
DELTA = 0.075  # conversion proies → prédateurs
GAMMA = 1.5    # mortalité des prédateurs sans proie

# Équilibre non trivial : (gamma/delta, alpha/beta) = (20, 10).

EQUIL_X = GAMMA / DELTA
EQUIL_Y = ALPHA / BETA

rhs = (_t, y) => {
    list(
        ALPHA * y[0] - BETA * y[0] * y[1],
        DELTA * y[0] * y[1] - GAMMA * y[1],
    )
}

# Un maximum de proies est un zéro de dx/dt = x*(alpha - beta*y). Passé en
# événement, il est localisé par recherche de racine à la précision de
# l'intégration, et non quantifié au pas de t_eval.

prey_peak = (_t, y) => { ALPHA * y[0] - BETA * y[0] * y[1] }

sol = solve_ivp(rhs, tuple(0.0, 100.0), list(40.0, 9.0), t_eval=numpy.linspace(0.0, 100.0, 4001), rtol=1e-9, atol=1e-12,
    events=prey_peak)

prey = sol.y[0]
predators = sol.y[1]
times = sol.t

print(f"⇒ Intégration : {sol.nfev} évaluations du champ, statut {sol.status}")
print(f"  Équilibre théorique : ({EQUIL_X}, {EQUIL_Y})")

# Oracle 1 : invariant de Lotka-Volterra conservé le long de la trajectoire.

invariant = DELTA * prey - GAMMA * numpy.log(prey) + BETA * predators - ALPHA * numpy.log(predators)
drift = float(numpy.max(numpy.abs(invariant - invariant[0])))
print(f"⇒ Dérive de l'invariant V : {drift:.2e} (tolérance 1e-6)")

# Oracle 2 : maxima de proies localisés par le solveur, périodes régulières.
#
# Ce que le solveur ne peut pas recevoir de Catnip, c'est le *réglage* de
# l'événement : scipy lit `direction` et `terminal` comme attributs de la
# fonction, or une lambda Catnip est un objet Rust sans attributs assignables
# (`AttributeError: no __dict__`). Sans direction, t_events porte les 37 zéros
# de dx/dt — maxima et minima confondus. Le tri revient donc à Catnip, et il
# tient en une comparaison : un maximum de proies est le seul zéro où la
# population dépasse l'équilibre.

peak_times = sol.t_events[0]
peak_states = sol.y_events[0]

peaks = list()
i = 0
while i < len(peak_times) {
    if peak_states[i][0] > EQUIL_X { peaks.append(float(peak_times[i])) }
    i = i + 1
}

periods = list()
i = 1
while i < len(peaks) {
    periods.append(peaks[i] - peaks[i - 1])
    i = i + 1
}

period_min = periods[0]
period_max = periods[0]
for p in periods {
    if p < period_min { period_min = p }
    if p > period_max { period_max = p }
}
spread = period_max - period_min

print(f"⇒ Maxima de proies : {len(peaks)} pics sur {len(peak_times)} zéros de dx/dt")
print(f"  Périodes : min {round(period_min, 6)} s, max {round(period_max, 6)} s, écart {spread:.2e} s")

# Moyennes sur un nombre entier de périodes : l'équilibre est la moyenne
# temporelle des populations (propriété classique du système).

last_peak = peaks[len(peaks) - 1]
first_peak = peaks[0]
mask = numpy.logical_and(times >= first_peak, times <= last_peak)
mean_prey = float(numpy.mean(prey[mask]))
mean_pred = float(numpy.mean(predators[mask]))
print(f"⇒ Moyennes sur {len(peaks) - 1} périodes : proies {round(mean_prey, 3)} (th. {EQUIL_X}), prédateurs {round(mean_pred, 3)} (th. {EQUIL_Y})")

oracle_ok = drift < 1e-6
oracle_ok = oracle_ok and len(peaks) >= 5 and spread < 1e-6
oracle_ok = oracle_ok and abs(mean_prey - EQUIL_X) < 0.1 and abs(mean_pred - EQUIL_Y) < 0.1

print()
print(f"⇒ Oracle (invariant + périodes + moyennes d'équilibre) : {oracle_ok}")

# Planche : séries temporelles en haut, portrait de phase en bas.

fig, axes = plt.subplots(2, 1, figsize=tuple(10, 8))

axes[0].plot(times, prey, color='#4c78a8', label='proies')
axes[0].plot(times, predators, color='#f58518', label='prédateurs')
axes[0].axhline(EQUIL_X, color='#4c78a8', linestyle=':', linewidth=0.8)
axes[0].axhline(EQUIL_Y, color='#f58518', linestyle=':', linewidth=0.8)
axes[0].set_xlabel('temps')
axes[0].set_ylabel('population')
axes[0].legend(loc='upper right')
axes[0].set_title('Lotka-Volterra : proies et prédateurs')

axes[1].plot(prey, predators, color='#54a24b', linewidth=0.8)
axes[1].scatter(list(EQUIL_X), list(EQUIL_Y), color='#e45756', marker='x', s=60, label='équilibre')
axes[1].set_xlabel('proies')
axes[1].set_ylabel('prédateurs')
axes[1].legend(loc='upper right')
axes[1].set_title('Portrait de phase')

fig.tight_layout()
out_path = output_dir / 'lotka_volterra.png'
fig.savefig(str(out_path), dpi=130, bbox_inches='tight')
print(f"⇒ Planche écrite : {out_path}")

# --no-browser garde un chemin headless (le PNG reste écrit ci-dessus).
if '--no-browser' not in import('sys').argv {
    http = import('http')
    b64 = import('base64').b64encode(out_path.read_bytes()).decode('ascii')
    http.serve(f'<!doctype html><meta charset="utf-8"><title>Lotka-Volterra (solve_ivp)</title><body style="margin:0;background:#0d1117"><img style="max-width:100%;display:block;margin:0 auto" src="data:image/png;base64,{b64}">',
        0, 'text/html; charset=utf-8', True)
}