Predator-Prey Dynamics

# Rate parameters (arbitrary units).
# PREY_BIRTH: per-capita prey growth rate in the absence of predators.
# PREDATION:  rate at which one predator removes prey (per predator per prey).
# PRED_GROWTH: rate at which predators gain from eating prey (per predator per prey).
# PRED_DEATH: per-capita predator death rate in the absence of prey.
PREY_BIRTH = 1.0
PREDATION = 0.1
PRED_GROWTH = 0.075
PRED_DEATH = 1.5

# Equilibrium populations: the fixed point of the ODE system.
# Setting dx/dt = 0 and dy/dt = 0 gives x* = PRED_DEATH/PRED_GROWTH
# and y* = PREY_BIRTH/PREDATION.
PREY_EQ = PRED_DEATH / PRED_GROWTH  # 20.0
PRED_EQ = PREY_BIRTH / PREDATION  # 10.0

# Initial conditions: start above the prey equilibrium and below the
# predator equilibrium to produce a visible oscillation.
PREY_INIT = 25.0
PRED_INIT = 5.0

# Simulate for ~10 complete cycles.  The small-oscillation period
# (linearised around equilibrium) is 2π / sqrt(PREY_BIRTH * PRED_DEATH).
# The true nonlinear period is ~5.30 for these initial conditions,
# so T_MAX = 53.0 covers roughly 10 full cycles.
PERIOD_APPROX = 2 * np.pi / np.sqrt(PREY_BIRTH * PRED_DEATH)
T_MAX = 53.0

# Dense output resolution: one point per 0.05 time units.
N_EVAL = int(T_MAX / 0.05) + 1

def rhs(t, z):
    """Return [dx/dt, dy/dt] for the Lotka-Volterra equations.

    z[0] = x (prey population)
    z[1] = y (predator population)
    """
    x, y = z
    dxdt = PREY_BIRTH * x - PREDATION * x * y
    dydt = PRED_GROWTH * x * y - PRED_DEATH * y
    return [dxdt, dydt]

def solve():
    """Integrate the Lotka-Volterra equations and return results as a DataFrame."""
    t_eval = np.linspace(0, T_MAX, N_EVAL)
    sol = solve_ivp(
        rhs,
        [0, T_MAX],
        [PREY_INIT, PRED_INIT],
        t_eval=t_eval,
        method="RK45",
        rtol=1e-8,
        atol=1e-10,
    )
    return pl.DataFrame({"t": sol.t, "prey": sol.y[0], "predator": sol.y[1]})

def plot_time_series(df):
    """Return an Altair chart of prey and predator populations over time."""
    long = df.unpivot(
        index="t",
        on=["prey", "predator"],
        variable_name="species",
        value_name="population",
    )
    return (
        alt.Chart(long)
        .mark_line()
        .encode(
            x=alt.X("t", title="Time"),
            y=alt.Y("population", title="Population"),
            color=alt.Color("species:N", title="Species"),
        )
        .properties(width=400, height=250, title="Population over time")
    )


def plot_trajectory(df):
    """Return an Altair chart of the predator-prey population trajectory."""
    return (
        alt.Chart(df)
        .mark_line()
        .encode(
            x=alt.X("prey", title="Prey"),
            y=alt.Y("predator", title="Predator"),
        )
        .properties(width=300, height=300, title="Population trajectory")
    )

import numpy as np
from scipy.signal import find_peaks
from lotka import (
    PRED_EQ,
    PERIOD_APPROX,
    PREY_EQ,
    solve,
)


def test_populations_positive():
    # The Lotka-Volterra equations conserve the positive quadrant: if x > 0
    # and y > 0 initially, both remain positive for all t > 0.  A negative
    # value would mean a population went extinct, which this model cannot
    # produce analytically; it would indicate a numerical error.
    df = solve()
    assert (df["prey"] > 0).all()
    assert (df["predator"] > 0).all()



def test_volterra_principle():
    # Volterra's averaging principle: the time mean of x(t) over any integer
    # number of complete cycles equals x* = PRED_DEATH/PRED_GROWTH, and
    # similarly for y(t).  T_MAX covers ~10 complete cycles, so the fractional-
    # cycle bias is small.  The measured errors are < 0.2%; a tolerance of
    # 1% (relative) gives a 5x safety margin.
    df = solve()
    x_mean = df["prey"].mean()
    y_mean = df["predator"].mean()
    assert abs(x_mean - PREY_EQ) / PREY_EQ < 0.01
    assert abs(y_mean - PRED_EQ) / PRED_EQ < 0.01


def test_period():
    # The small-oscillation period (linearisation around equilibrium) is
    # 2π / sqrt(PREY_BIRTH * PRED_DEATH).  For finite-amplitude oscillations
    # the true period differs by O(amplitude²).  With our initial conditions
    # the measured period is ~5.30 vs the linearised estimate of ~5.13, a
    # 3.3% difference.  A tolerance of 10% gives a 3x safety margin.
    df = solve()
    peaks, _ = find_peaks(df["prey"].to_numpy())
    t_peaks = df["t"].to_numpy()[peaks]
    assert len(t_peaks) >= 3, "too few peaks detected to estimate period"
    measured_period = np.diff(t_peaks).mean()
    assert abs(measured_period - PERIOD_APPROX) / PERIOD_APPROX < 0.10