restart:
read ("DMB.txt"):  # Load the DMB package
with(plots):
with(LinearAlgebra):

# ============================================================
# USER PARAMETERS - CHANGE THESE VALUES TO EXPERIMENT
# ============================================================

# --- SIR MODEL PARAMETERS (Conjecture 1) ---
# Constraints: 0 < a, b, c  and  0 < b + c <= 1
a_sir := 2.0:      # transmission rate
b_sir := 0.3:      # birth/death rate  
c_sir := 0.2:      # recovery rate

# --- SI MODEL STABLE CASE PARAMETERS (Conjecture 2) ---
# Constraints: 0 < a, b  and  a/b > 1 (for endemic equilibrium)
a_si_stable := 2.0:    # transmission rate
b_si_stable := 0.3:    # birth/death rate

# --- SI MODEL UNSTABLE CASE PARAMETERS (Conjecture 2) ---
# Try values where b is large (close to 1) to see instability
a_si_unstable := 3.0:  # transmission rate
b_si_unstable := 0.9:  # birth/death rate

# --- INITIAL CONDITIONS ---
x0_sir := 0.1:     # initial infectives (SIR)
y0_sir := 0.8:     # initial susceptibles (SIR)
x0_si := 0.1:      # initial condition for SI stable
x1_si := 0.15:     # second initial condition for SI (2nd order)
x0_uns := 0.5:     # initial condition for SI unstable

# --- SIMULATION LENGTH ---
n_iter_sir := 150:     # iterations for SIR
n_iter_si := 100:      # iterations for SI stable
n_iter_uns := 200:     # iterations for SI unstable

# --- SPECTRAL RADIUS SWEEP PARAMETERS ---
a_sweep := 2.0:        # fixed 'a' for sweep
b_min := 0.25:         # minimum b to test
b_max := 0.94:         # maximum b to test
b_step := 0.01:        # step size for b

# ============================================================
# END OF USER PARAMETERS - CODE BELOW RUNS AUTOMATICALLY
# ============================================================

print("=== ALLEN (2001) CONJECTURES ANALYSIS ==="):
print(""):

# ============================================================
# CONJECTURE 1: SIR MODEL
# ============================================================

print("=== CONJECTURE 1: SIR MODEL ==="):
print(`Parameters: a =`, a_sir, `, b =`, b_sir, `, c =`, c_sir):

# Check parameter constraints
if b_sir + c_sir > 1 then
    print(`WARNING: b + c > 1 violates model constraints!`):
end if:

R0_sir := a_sir/(b_sir + c_sir):
print(`Basic Reproduction Number R0 =`, R0_sir):

if R0_sir <= 1 then
    print(`WARNING: R0 <= 1, disease will die out (no endemic equilibrium)`):
end if:

# Get transformation using DMB package
T_sir := AllenSIR(a_sir, b_sir, c_sir, x, y):

# Find endemic equilibrium numerically
eq_sir := (b_sir + c_sir)*xs = (1 - xs*(1 + c_sir/b_sir))*(1 - exp(-a_sir*xs)):
xstar_sir := fsolve(eq_sir, xs = 0.01..0.99):
if type(xstar_sir, numeric) then
    ystar_sir := 1 - xstar_sir*(1 + c_sir/b_sir):
    print(`Endemic equilibrium: x* =`, xstar_sir, `, y* =`, ystar_sir):
    
    # Check stability using JAC from DMB
    J_sir := JAC(T_sir, [x, y]):
    J_sir_eval := evalf(subs({x = xstar_sir, y = ystar_sir}, J_sir)):
    print(`Jacobian at equilibrium:`):
    print(Matrix(J_sir_eval)):
    
    is_stable_sir := IsDisStable(J_sir_eval):
    eigs_sir := Eigenvalues(Matrix(J_sir_eval)):
    rho_sir := max(abs(eigs_sir[1]), abs(eigs_sir[2])):
    print(`Spectral radius rho =`, evalf(rho_sir)):
    print(`Is stable?`, is_stable_sir):
else
    print(`Could not find endemic equilibrium`):
end if:

# Simulate using ORB
print(``):
print(`Simulating from initial conditions: x0 =`, x0_sir, `, y0 =`, y0_sir):
full_orbit_sir := ORB(T_sir, [x, y], [x0_sir, y0_sir], 0, n_iter_sir):
orbit_sir := ORB(T_sir, [x, y], [x0_sir, y0_sir], n_iter_sir-10, n_iter_sir):
print(`Final 10 iterates:`):
for i from 1 to nops(orbit_sir) do
    print(orbit_sir[i]):
od:

if type(xstar_sir, numeric) then
    final_pt := orbit_sir[-1]:
    dist := evalf(sqrt((final_pt[1] - xstar_sir)^2 + (final_pt[2] - ystar_sir)^2)):
    print(`Distance to equilibrium:`, dist):
end if:

# Phase portrait
sir_plot := listplot(full_orbit_sir, color = blue, 
    title = cat("SIR Phase Portrait (a=", a_sir, ", b=", b_sir, ", c=", c_sir, ")"),
    labels = ["x (Infectives)", "y (Susceptibles)"],
    style = line, symbol = solidcircle, symbolsize = 5):

# Time series for x (infectives)
sir_x_data := [seq([i, full_orbit_sir[i][1]], i = 1..nops(full_orbit_sir))]:
sir_x_plot := listplot(sir_x_data, color = blue,
    title = cat("SIR Infectives Time Series (a=", a_sir, ", b=", b_sir, ", c=", c_sir, ")"),
    labels = ["n", "x"]):

# Time series for y (susceptibles)
sir_y_data := [seq([i, full_orbit_sir[i][2]], i = 1..nops(full_orbit_sir))]:
sir_y_plot := listplot(sir_y_data, color = red,
    title = cat("SIR Susceptibles Time Series (a=", a_sir, ", b=", b_sir, ", c=", c_sir, ")"),
    labels = ["n", "y"]):

# Combined time series (both x and y)
sir_time_plot := display(
    listplot(sir_x_data, color = blue, legend = "x (Infectives)"),
    listplot(sir_y_data, color = red, legend = "y (Susceptibles)"),
    title = cat("SIR Time Series (a=", a_sir, ", b=", b_sir, ", c=", c_sir, ")"),
    labels = ["n", "population"]):

# ============================================================
# CONJECTURE 2: SI MODEL (STABLE CASE)
# ============================================================

print(``):
print("=== CONJECTURE 2: SI MODEL (STABLE CASE) ==="):
print(`Parameters: a =`, a_si_stable, `, b =`, b_si_stable):

R0_si := a_si_stable/b_si_stable:
print(`R0 =`, R0_si):

if R0_si <= 1 then
    print(`WARNING: R0 <= 1, no endemic equilibrium exists`):
end if:

# Define transformation for 2nd order system
T_si := [z[1]*(1 - b_si_stable) + (1 - z[1])*(1 - exp(-a_si_stable*z[2])), z[1]]:
vars_si := [z[1], z[2]]:

# Find steady state
eq_si := b_si_stable*xs = (1 - xs)*(1 - exp(-a_si_stable*xs)):
xstar_si := fsolve(eq_si, xs = 0.01..0.99):

if type(xstar_si, numeric) then
    print(`Endemic equilibrium x* =`, xstar_si):
    
    # Check stability
    J_si := JAC(T_si, vars_si):
    J_si_eval := evalf(subs({z[1] = xstar_si, z[2] = xstar_si}, J_si)):
    print(`Jacobian at equilibrium:`):
    print(Matrix(J_si_eval)):
    
    is_stable_si := IsDisStable(J_si_eval):
    eigs_si := Eigenvalues(Matrix(J_si_eval)):
    rho_si := max(abs(eigs_si[1]), abs(eigs_si[2])):
    print(`Spectral radius rho =`, evalf(rho_si)):
    print(`Is stable?`, is_stable_si):
else
    print(`Could not find endemic equilibrium`):
end if:

# Simulate
print(``):
print(`Simulating from initial conditions: x0 =`, x0_si, `, x1 =`, x1_si):
full_orbit_si := ORB(T_si, vars_si, [x0_si, x1_si], 0, n_iter_si):
orbit_si := ORB(T_si, vars_si, [x0_si, x1_si], n_iter_si-10, n_iter_si):
print(`Final 10 iterates:`):
for i from 1 to nops(orbit_si) do
    print(orbit_si[i]):
od:

# Time series plot
si_stable_data := [seq([i, full_orbit_si[i][1]], i = 1..nops(full_orbit_si))]:
si_stable_plot := listplot(si_stable_data, color = blue,
    title = cat("SI Stable Time Series (a=", a_si_stable, ", b=", b_si_stable, ")"),
    labels = ["n", "x"]):

# Phase portrait (z[1] = x_n vs z[2] = x_{n-1})
si_stable_phase := listplot(full_orbit_si, color = blue,
    title = cat("SI Stable Phase Portrait (a=", a_si_stable, ", b=", b_si_stable, ")"),
    labels = ["x_n", "x_{n-1}"],
    style = line, symbol = solidcircle, symbolsize = 5):

# ============================================================
# CONJECTURE 2: SI MODEL (UNSTABLE CASE)
# ============================================================

print(``):
print("=== CONJECTURE 2: SI MODEL (UNSTABLE CASE) ==="):
print(`Parameters: a =`, a_si_unstable, `, b =`, b_si_unstable):

R0_uns := a_si_unstable/b_si_unstable:
print(`R0 =`, R0_uns):

T_uns := [z[1]*(1 - b_si_unstable) + (1 - z[1])*(1 - exp(-a_si_unstable*z[2])), z[1]]:

# Find steady state
eq_uns := b_si_unstable*xs = (1 - xs)*(1 - exp(-a_si_unstable*xs)):
xstar_uns := fsolve(eq_uns, xs = 0.01..0.99):

if type(xstar_uns, numeric) then
    print(`Endemic equilibrium x* =`, xstar_uns):
    
    # Check stability
    J_uns := JAC(T_uns, vars_si):
    J_uns_eval := evalf(subs({z[1] = xstar_uns, z[2] = xstar_uns}, J_uns)):
    print(`Jacobian at equilibrium:`):
    print(Matrix(J_uns_eval)):
    
    is_stable_uns := IsDisStable(J_uns_eval):
    eigs_uns := Eigenvalues(Matrix(J_uns_eval)):
    rho_uns := max(abs(eigs_uns[1]), abs(eigs_uns[2])):
    print(`Spectral radius rho =`, evalf(rho_uns)):
    print(`Is stable?`, is_stable_uns):
    if not is_stable_uns then
        print(`UNSTABLE because rho > 1`):
    end if:
else
    print(`Could not find endemic equilibrium`):
end if:

# Simulate
full_orbit_uns := ORB(T_uns, vars_si, [x0_uns, x0_uns], 0, n_iter_uns):
si_unstable_data := [seq([i, full_orbit_uns[i][1]], i = 1..nops(full_orbit_uns))]:
si_unstable_plot := listplot(si_unstable_data, color = blue,
    title = cat("SI Unstable Time Series (a=", a_si_unstable, ", b=", b_si_unstable, ")"),
    labels = ["n", "x"]):

# Phase portrait (z[1] = x_n vs z[2] = x_{n-1})
si_unstable_phase := listplot(full_orbit_uns, color = blue,
    title = cat("SI Unstable Phase Portrait (a=", a_si_unstable, ", b=", b_si_unstable, ")"),
    labels = ["x_n", "x_{n-1}"],
    style = line, symbol = solidcircle, symbolsize = 5):

# ============================================================
# SPECTRAL RADIUS VS B ANALYSIS
# ============================================================

print(``):
print("=== SPECTRAL RADIUS VS B ANALYSIS ==="):
print(`Sweeping b from`, b_min, `to`, b_max, `with a =`, a_sweep):

rho_data := []:
b_critical := NULL:

for bval from b_min by b_step to b_max do
    if a_sweep/bval > 1.01 then
        xsval := fsolve(bval*xx = (1 - xx)*(1 - exp(-a_sweep*xx)), xx = 0.01..0.99):
        
        if type(xsval, numeric) and xsval > 0 and xsval < 1 then
            T_sweep := [z[1]*(1 - bval) + (1 - z[1])*(1 - exp(-a_sweep*z[2])), z[1]]:
            J_sweep := JAC(T_sweep, vars_si):
            J_sweep_eval := evalf(subs({z[1] = xsval, z[2] = xsval}, J_sweep)):
            
            eigsval := Eigenvalues(Matrix(J_sweep_eval)):
            rhoval := evalf(max(abs(eigsval[1]), abs(eigsval[2]))):
            rho_data := [op(rho_data), [bval, rhoval]]:
            
            # Track when stability is lost
            if nops(rho_data) > 1 then
                if rho_data[-2][2] < 1 and rhoval >= 1 and b_critical = NULL then
                    b_critical := (rho_data[-2][1] + bval)/2:
                end if:
            end if:
        end if:
    end if:
end do:

if b_critical <> NULL then
    print(`Critical b value (stability lost) ~ `, b_critical):
else
    print(`No stability transition found in range`):
end if:

# Create plots
rho_plot := listplot(rho_data, color = blue,
    title = cat("Spectral Radius vs b (a=", a_sweep, ")"),
    labels = ["b", "rho"]):

ref_line := listplot([[b_min, 1], [b_max, 1]], color = red, linestyle = dash):
rho_plot_combined := display(rho_plot, ref_line):

# ============================================================
# DISPLAY ALL PLOTS
# ============================================================

print(``):
print("=== PLOTS ==="):

print(`SIR - Phase Portrait:`):
display(sir_plot);

print(`SIR - Time Series (both variables):`):
display(sir_time_plot);

print(`SI Stable - Time Series:`):
display(si_stable_plot);

print(`SI Stable - Phase Portrait:`):
display(si_stable_phase);

print(`SI Unstable - Time Series:`):
display(si_unstable_plot);

print(`SI Unstable - Phase Portrait:`):
display(si_unstable_phase);

print(`Spectral Radius vs b:`):
display(rho_plot_combined);

# ============================================================
# SUMMARY
# ============================================================

print(``):
print("=== SUMMARY ==="):
print(``):
print(`SIR Model (a=`, a_sir, `, b=`, b_sir, `, c=`, c_sir, `):`):
print(`  R0 =`, R0_sir):
if type(rho_sir, numeric) then
    print(`  Spectral radius =`, evalf(rho_sir)):
    print(`  Stable =`, is_stable_sir):
end if:
print(``):
print(`SI Stable (a=`, a_si_stable, `, b=`, b_si_stable, `):`):
print(`  R0 =`, R0_si):
if type(rho_si, numeric) then
    print(`  Spectral radius =`, evalf(rho_si)):
    print(`  Stable =`, is_stable_si):
end if:
print(``):
print(`SI Unstable (a=`, a_si_unstable, `, b=`, b_si_unstable, `):`):
print(`  R0 =`, R0_uns):
if type(rho_uns, numeric) then
    print(`  Spectral radius =`, evalf(rho_uns)):
    print(`  Stable =`, is_stable_uns):
end if: