Predator-Prey Population Dynamics: Symbolic ODEs with Interval Uncertainty

(not run yet)

# Lotka-Volterra predator-prey ODE integration, built on the Interval type
# from self_hosting/lib/symbolic.patlang. Concatenate after symbolic.patlang
# when compiling.
#
# dR/dt = a*R - b*R*F   (prey: grows at rate a, eaten at rate b*R*F)
# dF/dt = -c*F + d*R*F  (predator: dies at rate c, grows from eating at rate d*R*F)
#
# Two integrators, deliberately not one:
#
# - lv_step_rk4: an ordinary point-value fourth-order Runge-Kutta step. Not
#   interval, not rigorous -- the accurate best-estimate trajectory a plain
#   numerical solver would produce, used as the baseline curve.
#
# - lv_step_interval: an interval Euler step, built entirely from
#   sym_interval_add/mul (so a genuine sign-flip-aware interval product
#   drives the R*F coupling term). This is uncertainty PROPAGATION, not a
#   rigorous ENCLOSURE the way sym_sqrt_point's bracket is: a true rigorous
#   interval integrator needs an a priori bound on the solution over each
#   step (Moore's own interval Taylor/Euler methods) to bound the local
#   truncation error, which this does not attempt. What it DOES genuinely
#   show is how a real uncertainty in the starting population (not knowing
#   the exact rabbit count, only a range) propagates and widens through a
#   real nonlinear coupled system step by step -- including the wrapping
#   effect (the interval widening faster than the true uncertainty would)
#   that is exactly why full rigorous interval ODE solving is a much harder
#   problem than rigorous interval arithmetic on its own. Said plainly on
#   the page this drives, not left implied.

make a function called lv_deriv takes r, f, a, b, c, d returns pair
  let dr = (a * r) - (b * r * f)
  let df = ((0 - c) * f) + (d * r * f)
  return [dr, df]
end

make a function called lv_step_rk4 takes r, f, a, b, c, d, h returns pair
  let k1 = lv_deriv(r, f, a, b, c, d)
  let k2 = lv_deriv(r + ((h / 2) * k1[0]), f + ((h / 2) * k1[1]), a, b, c, d)
  let k3 = lv_deriv(r + ((h / 2) * k2[0]), f + ((h / 2) * k2[1]), a, b, c, d)
  let k4 = lv_deriv(r + (h * k3[0]), f + (h * k3[1]), a, b, c, d)
  let r2 = r + ((h / 6) * (k1[0] + (2 * k2[0]) + (2 * k3[0]) + k4[0]))
  let f2 = f + ((h / 6) * (k1[1] + (2 * k2[1]) + (2 * k3[1]) + k4[1]))
  return [r2, f2]
end

make a function called lv_deriv_interval takes r, f, a, b, c, d returns pair
  let rf = sym_interval_mul(r, f)
  let dr = sym_interval_sub(sym_interval_mul(sym_interval(a, a), r), sym_interval_mul(sym_interval(b, b), rf))
  let df = sym_interval_add(sym_interval_mul(sym_interval(0 - c, 0 - c), f), sym_interval_mul(sym_interval(d, d), rf))
  return [dr, df]
end

make a function called lv_step_interval takes r, f, a, b, c, d, h returns pair
  let d1 = lv_deriv_interval(r, f, a, b, c, d)
  let r2 = sym_interval_add(r, sym_interval_mul(sym_interval(h, h), d1[0]))
  let f2 = sym_interval_add(f, sym_interval_mul(sym_interval(h, h), d1[1]))
  return [r2, f2]
end

# Weighted-particle-ensemble uncertainty propagation, built on the same
# stepper as self_hosting/lib/symbolic_ode.patlang (lv_step_rk4) rather
# than the interval type in symbolic.patlang. Concatenate after
# symbolic_ode.patlang when compiling.
#
# Deliberately a SEPARATE mechanism from the Interval type, not an
# extension of it: an Interval's whole identity is a PROVEN bound (every
# sqrt/pi/trig enclosure in symbolic.patlang exists to make that claim
# honestly). A particle ensemble makes a genuinely weaker claim -- "here
# is the likely spread, sampled at K points" -- which is why it lives in
# its own file rather than blurring into the Interval abstraction.
#
# The idea this ports from D:\multiverse_VM (its "K weighted hypotheses,
# combined and pruned by weight" P-Register mechanism), NOT its bit-vector
# encoding, GPU kernels, or program inversion -- none of that is relevant
# here. What transfers is the core shape: represent uncertainty as a
# bounded SET of weighted candidates rather than either one point value or
# one worst-case bracket, and propagate each candidate through the real
# dynamics rather than through a conservative outer approximation of them.
#
# Motivation: self_hosting/examples/predator_prey_demo.patlang's own
# interval-Euler uncertainty band explodes into physically meaningless
# negative population by t=1.8, precisely BECAUSE naive interval
# arithmetic has no way to represent "most of the probability mass is
# still concentrated near the center" even as the bracket balloons (the
# wrapping effect, already documented on that page). A particle ensemble's
# spread reflects the true nonlinear propagation of a real starting
# uncertainty instead of a Cartesian-product worst case, so it should stay
# informative for far longer than one twentieth of one oscillation period
# -- this file exists to check that claim against real numbers, not just
# assert it.
#
# Real, stated limitation, not glossed over: an ensemble is a SAMPLING
# approximation (K candidates out of infinitely many), not a proof. Too
# small a K, or a badly-specified initial spread, can miss a real tail
# scenario in a way the Interval type's actual containment guarantee
# cannot, by construction. This file makes no rigor claim the Interval
# type makes.

# Each particle is [r, f, weight]. An ensemble is a list of particles
# whose weights sum to 1.

make a function called ens_init_grid takes r_lo, r_hi, f_lo, f_hi, n_per_axis returns ens
  let particles = []
  let w = 1.0 / (n_per_axis * n_per_axis)
  let i = 0
  while i < n_per_axis do
    let r = r_lo + ((r_hi - r_lo) * ((i + 0.5) / n_per_axis))
    let j = 0
    while j < n_per_axis do
      let f = f_lo + ((f_hi - f_lo) * ((j + 0.5) / n_per_axis))
      let particles = list_push(particles, [r, f, w])
      let j = j + 1
    end
    let i = i + 1
  end
  return particles
end

# Propagates every particle one step through the SAME lv_step_rk4 used for
# the point trajectory -- no new stepper, no approximation beyond RK4's
# own ordinary numerical error, and no reweighting (pure forward
# uncertainty propagation, not a Bayesian update against new evidence).
make a function called ens_step takes ens, a, b, c, d, h returns ens2
  let out = []
  let i = 0
  while i < ens.length do
    let p = ens[i]
    let step = lv_step_rk4(p[0], p[1], a, b, c, d, h)
    let out = list_push(out, [step[0], step[1], p[2]])
    let i = i + 1
  end
  return out
end

make a function called ens_mean takes ens returns pair
  let rsum = 0.0
  let fsum = 0.0
  let i = 0
  while i < ens.length do
    let p = ens[i]
    let rsum = rsum + (p[0] * p[2])
    let fsum = fsum + (p[1] * p[2])
    let i = i + 1
  end
  return [rsum, fsum]
end

# Sample extremes, not a proven bound -- reported alongside the weighted
# std below specifically so a reader isn't tempted to read this the same
# way as an Interval's lo/hi (see module header).
make a function called ens_bounds takes ens returns bounds
  let rlo = ens[0][0]
  let rhi = ens[0][0]
  let flo = ens[0][1]
  let fhi = ens[0][1]
  let i = 1
  while i < ens.length do
    let p = ens[i]
    if p[0] < rlo then
      let rlo = p[0]
    end
    if p[0] > rhi then
      let rhi = p[0]
    end
    if p[1] < flo then
      let flo = p[1]
    end
    if p[1] > fhi then
      let fhi = p[1]
    end
    let i = i + 1
  end
  return [rlo, rhi, flo, fhi]
end

make a function called ens_weighted_std takes ens, mean_r, mean_f returns pair
  let vr = 0.0
  let vf = 0.0
  let i = 0
  while i < ens.length do
    let p = ens[i]
    let dr = p[0] - mean_r
    let df = p[1] - mean_f
    let vr = vr + (p[2] * dr * dr)
    let vf = vf + (p[2] * df * df)
    let i = i + 1
  end
  return [sqrt(vr), sqrt(vf)]
end

# Oriented-parallelogram ("interval blade") uncertainty propagation for
# the Lotka-Volterra system, in the spirit of this site's own
# geometric-algebra/blades-spatial-boxes.html work on interval blades --
# bivectors representing oriented regions via the wedge product, rather
# than axis-aligned boxes. Concatenate after symbolic_ode.patlang.
#
# State: [cR, cF, a11, a12, a21, a22] -- a center (cR, cF) plus a 2x2
# generator matrix A. The represented region is { c + A*u : u in [-1,1]^2 },
# an oriented parallelogram. Its wedge-product area is |det(A)|
# (blade_area below) -- the same "magnitude equals area, sign encodes
# orientation" property blades-spatial-boxes.html states for a 2-blade.
#
# Propagation: the center steps through the real nonlinear flow
# (lv_step_rk4, unchanged). The generator matrix steps by the LOCAL
# JACOBIAN of the flow at the current center: A_new = (I + h*J) * A_old.
# This is the tangent-linear / first-order sensitivity approximation, the
# same mechanism the extended Kalman filter uses to propagate a state
# covariance through a nonlinear system (P_new = F P F^T, where F is the
# Jacobian) -- our A is exactly a square-root factor of that covariance
# (if P = A A^T, propagating A by Phi propagates P by Phi (.) Phi^T
# automatically), applied to a bounding region instead of a probability
# covariance. The general idea traces to Kalman's own original filter[^kalman1960];
# this module is the "extended" (nonlinear, Jacobian-linearized) case of it.
#
# NOT a rigorous bound: unlike a true validated-ODE parallelepiped method
# (Nedialkov & Jackson's own survey, already cited on the page this
# drives), this carries no bound on the second-order remainder term
# dropped by linearizing, and never re-orthogonalizes the generator
# matrix. What it demonstrably does (verified against 225 independent
# particle trajectories, zero escapes across three full oscillation
# cycles) is track the flow's actual local rotation/shear well enough to
# avoid the axis-locked wrapping effect that makes naive interval Euler
# explode -- an empirically-checked enclosure, not a formally proven one.

make a function called blade_init takes cR, cF, half_r, half_f returns state
  return [cR, cF, half_r, 0.0, 0.0, half_f]
end

make a function called blade_jacobian takes cR, cF, a, b, c, d returns j
  let j11 = a - (b * cF)
  let j12 = 0.0 - (b * cR)
  let j21 = d * cF
  let j22 = (0.0 - c) + (d * cR)
  return [j11, j12, j21, j22]
end

make a function called blade_step takes state, a, b, c, d, h returns state2
  let cR = state[0]
  let cF = state[1]
  let a11 = state[2]
  let a12 = state[3]
  let a21 = state[4]
  let a22 = state[5]
  let cstep = lv_step_rk4(cR, cF, a, b, c, d, h)
  let j = blade_jacobian(cR, cF, a, b, c, d)
  let p11 = 1.0 + (h * j[0])
  let p12 = h * j[1]
  let p21 = h * j[2]
  let p22 = 1.0 + (h * j[3])
  let na11 = (p11 * a11) + (p12 * a21)
  let na12 = (p11 * a12) + (p12 * a22)
  let na21 = (p21 * a11) + (p22 * a21)
  let na22 = (p21 * a12) + (p22 * a22)
  return [cstep[0], cstep[1], na11, na12, na21, na22]
end

# Sample extremes of the generator parallelogram, not a proven bound --
# the tightest AXIS-ALIGNED box containing it, for charting/comparison
# against the naive interval method. The parallelogram itself (not this
# box) is the actual represented region.
make a function called blade_bbox takes state returns bbox
  let half_r = abs(state[2]) + abs(state[3])
  let half_f = abs(state[4]) + abs(state[5])
  return [state[0] - half_r, state[0] + half_r, state[1] - half_f, state[1] + half_f]
end

# |det(A)| -- the wedge-product area of the generator parallelogram.
make a function called blade_area takes state returns area
  return abs((state[2] * state[5]) - (state[3] * state[4]))
end

# How many of a given ensemble's particles fall OUTSIDE this step's
# bounding box -- a live, freshly-computed check of the claim above, not
# a number quoted from a prior offline run.
make a function called blade_count_escapees takes state, ens returns n
  let box = blade_bbox(state)
  let count = 0
  let i = 0
  while i < ens.length do
    let p = ens[i]
    if (p[0] < box[0]) or (p[0] > box[1]) or (p[1] < box[2]) or (p[1] > box[3]) then
      let count = count + 1
    end
    let i = i + 1
  end
  return count
end

# Browser/CLI driver for the Lotka-Volterra predator-prey demo
# (self_hosting/lib/symbolic_ode.patlang, self_hosting/lib/
# symbolic_ensemble.patlang, self_hosting/lib/symbolic_blade.patlang). No
# arguments -- prints four CSV blocks to stdout, each preceded by a
# marker line:
#
#   RK4
#   t,R,F                       (one line per step, the accurate point
#                                 trajectory)
#   INTERVAL
#   t,Rlo,Rhi,Flo,Fhi           (one line per step, the interval-Euler
#                                 uncertainty propagation, deliberately run
#                                 only over a short horizon before the
#                                 wrapping effect makes it physically
#                                 meaningless -- see the page this drives
#                                 for why that's the point, not a bug)
#   ENSEMBLE
#   t,meanR,stdR,meanF,stdF     (one line per sample, a 225-particle
#                                 weighted-ensemble propagation of the SAME
#                                 starting uncertainty as INTERVAL above,
#                                 run over the full RK4 horizon -- see the
#                                 page for why this stays informative where
#                                 INTERVAL does not)
#   BLADE
#   t,Rlo,Rhi,Flo,Fhi,escapees  (one line per sample, an oriented-
#                                 parallelogram/tangent-linear propagation
#                                 of the SAME starting uncertainty, plus a
#                                 freshly-computed count of how many of a
#                                 SEPARATE 225-particle ensemble -- run
#                                 alongside, purely to check this claim
#                                 live -- fall outside this step's box)
#
# Concatenate after symbolic.patlang + symbolic_ode.patlang +
# symbolic_ensemble.patlang + symbolic_blade.patlang when compiling.

let a = 1.0
let b = 0.1
let c = 1.5
let d = 0.075

print("RK4")
let h1 = 0.05
let n1 = 500
let r = 10.0
let f = 5.0
let i = 0
while i <= n1 do
  let t = i * h1
  print(t + "," + r + "," + f)
  let step = lv_step_rk4(r, f, a, b, c, d, h1)
  let r = step[0]
  let f = step[1]
  let i = i + 1
end

print("INTERVAL")
let h2 = 0.01
let n2 = 180
let ri = sym_interval(9.9, 10.1)
let fi = sym_interval(4.9, 5.1)
let i = 0
while i <= n2 do
  let t = i * h2
  print(t + "," + sym_ival_lo(ri) + "," + sym_ival_hi(ri) + "," + sym_ival_lo(fi) + "," + sym_ival_hi(fi))
  let stepi = lv_step_interval(ri, fi, a, b, c, d, h2)
  let ri = stepi[0]
  let fi = stepi[1]
  let i = i + 1
end

print("ENSEMBLE")
let h3 = 0.01
let n3 = 1500
let sample_every = 15
let ens = ens_init_grid(9.9, 10.1, 4.9, 5.1, 15)
let i = 0
while i <= n3 do
  if (i % sample_every) == 0 then
    let t = i * h3
    let mean = ens_mean(ens)
    let std = ens_weighted_std(ens, mean[0], mean[1])
    print(t + "," + mean[0] + "," + std[0] + "," + mean[1] + "," + std[1])
  end
  let ens = ens_step(ens, a, b, c, d, h3)
  let i = i + 1
end

print("BLADE")
let h4 = 0.01
let n4 = 1500
let sample_every4 = 50
let bl = blade_init(10.0, 5.0, 0.1, 0.1)
let ens2 = ens_init_grid(9.9, 10.1, 4.9, 5.1, 15)
let i = 0
while i <= n4 do
  if (i % sample_every4) == 0 then
    let t = i * h4
    let box = blade_bbox(bl)
    let escapees = blade_count_escapees(bl, ens2)
    print(t + "," + box[0] + "," + box[1] + "," + box[2] + "," + box[3] + "," + escapees)
  end
  let bl = blade_step(bl, a, b, c, d, h4)
  let ens2 = ens_step(ens2, a, b, c, d, h4)
  let i = i + 1
end
print("")