import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

plt.rcParams.update({'font.size': 20})

def simulate_glv(r, A, N0, t_span, t_eval):
    """
    Simulate the generalized Lotka-Volterra (GLV) system:
       dN_i/dt = N_i * (r_i + sum_j A_ij * N_j)
    """
    def glv(t, N):
        return N * (r + A.dot(N))
    sol = solve_ivp(glv, t_span, N0, t_eval=t_eval, rtol=1e-8)
    return sol

#############################
# Case 1: 3 Species, One Species Dies Out
#############################
# Here species 3 is disadvantaged by both a low intrinsic growth rate 
# and strong negative interactions.
r1 = np.array([1.0, 1.0, 0.1])
A1 = np.array([
    [-1.0,  0.2,  0.1],
    [ 0.1, -1.0,  0.2],
    [-1.2, -1.2, -1.0]   # Strong competitive pressure on species 3
])
N0_1 = np.array([0.5, 0.5, 0.5])
t_span1 = (0, 10)
t_eval1 = np.linspace(t_span1[0], t_span1[1], 200)
sol1 = simulate_glv(r1, A1, N0_1, t_span1, t_eval1)

plt.figure(figsize=(10, 6))
for i in range(3):
    plt.plot(sol1.t, sol1.y[i], label=f"Species {i+1}")
plt.xlabel("Time")
plt.ylabel("Population Density")
plt.title("Case 1: 3 Species (One Dies Out)")
plt.legend()
plt.savefig('glv_extinction.pdf')
plt.show()

#############################
# Case 2: 3 Species, Stable Limit Cycle (Cyclic Competitors)
#############################
# This system is set up in a rock--paper--scissors (May--Leonard) style.
# The equations take the form:
#   dN_i/dt = N_i * (1 - N_i - a_ij * N_j - a_ik * N_k)
# with cyclic intransitivity that promotes oscillations.
r2 = np.array([1.0, 1.0, 1.0])
# Note: The off-diagonals are chosen asymmetrically to induce cyclic competition.
A2 = np.array([
    [-10.0, -6.0, -12.0],
    [-14.0, -10.0, -2.0],
    [-8.0, -18.0, -10.0]
])
N0_2 = np.array([0.6, 0.4, 0.2])
t_span2 = (0, 500)
t_eval2 = np.linspace(t_span2[0], t_span2[1], 1500)
sol2 = simulate_glv(r2, A2, N0_2, t_span2, t_eval2)

plt.figure(figsize=(10, 6))
for i in range(3):
    plt.plot(sol2.t, sol2.y[i], label=f"Species {i+1}")
plt.xlabel("Time")
plt.ylabel("Population Density")
plt.title("Case 2: 3 Species with Stable Limit Cycle")
plt.legend()
plt.savefig('glv_limit_cycle.pdf')
plt.show()

