# python script for generating data for Figure 2 in 
# Parametric resonance, chaos and spatial structure in the Lotka--Volterra model
# Mohamed Swailem and Alastair M. Rucklidge
# Copyright 2025 University of Leeds
# Unless otherwise stated, this dataset is licensed under a Creative Commons 
# Attribution 4.0 International Licence: 
# https://creativecommons.org/licenses/by/4.0/.

# mean field ODE stability calculation, varying alpha
# run the ODEs for a transient and then see what we have left
# 0 = eqm point (alpha=1 only), 1, 2, 3... = periodic orbit, -1 = chaos
# vary alpha and output a series of values of (alpha, a, b) for the attractor
# based on LV_mean_field_stability.py verion 1.03

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

version=1.00 # 25 August 2025
# edited for publication

# parameters -- equilibrium point scaled to (astar,1)
# carrying capacity is 1/(kk0 + kk1*cos(omega*time)), where omega=omega0/n
# and omega0 is the imaginary part of the eigenvalue of the (astar,1) equilibrium
# alpha is the homotopy parameter: alpha=0 has time-dependent predation rate

kappa0 = 0.25
kappa1 = 0.24    # instability with alpha=1 is around kappa1=0.225
astar = 1
n = 0.70
alpha = 0.0  # stability calculations need alpha = 1

tol = 1e-6

# controlling the calculation

# first tranche
n_forcing_periods_transient = 1000
n_forcing_periods_keep = 100
alpha_start = 0.15
alpha_end   = 0.12
n_alpha = 151

# second tranche
# n_forcing_periods_transient = 10000
# n_forcing_periods_keep = 100
# alpha_start = 0.12
# alpha_end   = 0.09
# n_alpha = 151

# third tranche
# n_forcing_periods_transient = 10000
# n_forcing_periods_keep = 100
# alpha_start = 0.09
# alpha_end   = 0.06
# n_alpha = 151

# fourth tranche
# n_forcing_periods_transient = 10000
# n_forcing_periods_keep = 100
# alpha_start = 0.06
# alpha_end   = 0.03
# n_alpha = 151

# fifth tranche
# n_forcing_periods_transient = 10000
# n_forcing_periods_keep = 100
# alpha_start = 0.03
# alpha_end   = 0.00
# n_alpha = 151

alpha = alpha_start

# derived quantities
lamda0 = (1 - (1 + astar) * kappa0) / astar
# lamda1 = (1 - (1 + astar) * (kappa0 + kappa1*cos(omega*time)) / astar
omega0squared = (1-kappa0)*(1-kappa0)/astar - kappa0*(1-kappa0) - kappa0*kappa0/4
omega0 = math.sqrt(omega0squared)
omega = omega0 / n

a_ic = astar * 1.001
b_ic =     1 * 1.001

save_figs_as_pdfs = 1

# construct a filename
file_name_1 = "LV_ODE_alpha"
file_name_2 = "_kappa0_" + "{:06.4f}".format(kappa0) + \
              "_kappa1_" + "{:06.4f}".format(kappa1) + \
              "_n" + "{:06.4f}".format(n)
file_name_2 = file_name_2.replace(".", "p")

verbose = 0

LV_mean_field_alpha_output_filename = file_name_1 + file_name_2 + ".txt"
LV_mean_field_alpha_output_file = open(LV_mean_field_alpha_output_filename, mode = "wt", buffering=1)  # buffering=1 write to disk at the end of each line

print("astar, kappa0, kappa1, n:", "{:06.4f}".format(astar),  \
                                   "{:06.4f}".format(kappa0), \
                                   "{:06.4f}".format(kappa1), \
                                   "{:06.4f}".format(n))
print("alpha_start, alpha_end, n_alpha", "{:06.4f}".format(alpha_start), \
                                         "{:06.4f}".format(alpha_end), 
                                         "{:5d}".format(n_alpha))
print("lamda0, omega0, omega:", "{:06.4f}".format(lamda0), \
                                "{:06.4f}".format(omega0), \
                                "{:06.4f}".format(omega))
print("n_forcing_periods_transient, n_forcing_periods_keep, tol:", "{:5d}".format(n_forcing_periods_transient), \
                                                                   "{:5d}".format(n_forcing_periods_keep), \
                                                                   "{:10.3g}".format(tol))
print("Version: ", "{:6.2f}".format(version))
print(" ")
print("    alpha              a                    b       count  period n_transient")

print("# astar, kappa0, kappa1, n:", "{:06.4f}".format(astar),  \
                                     "{:06.4f}".format(kappa0), \
                                     "{:06.4f}".format(kappa1), \
                                     "{:06.4f}".format(n), file = LV_mean_field_alpha_output_file)
print("# alpha_start, alpha_end, n_alpha", "{:06.4f}".format(alpha_start), \
                                           "{:06.4f}".format(alpha_end), 
                                           "{:5d}".format(n_alpha), file = LV_mean_field_alpha_output_file)
print("# lamda0, omega0, omega:", "{:06.4f}".format(lamda0), \
                                  "{:06.4f}".format(omega0), \
                                  "{:06.4f}".format(omega), file = LV_mean_field_alpha_output_file)
print("# n_forcing_periods_transient, n_forcing_periods_keep, tol:", "{:5d}".format(n_forcing_periods_transient), \
                                                                     "{:5d}".format(n_forcing_periods_keep), \
                                                                     "{:10.3g}".format(tol), file = LV_mean_field_alpha_output_file)
print("# Version: ", "{:6.2f}".format(version), file = LV_mean_field_alpha_output_file)
print("# ", file = LV_mean_field_alpha_output_file)
print("#   alpha              a                    b       count  period n_transient", file = LV_mean_field_alpha_output_file)

for alpha_count in range(n_alpha):
    alpha = alpha_start + alpha_count * (alpha_end - alpha_start) / (n_alpha-1)
    if (verbose == 1):
        print("alpha_count, alpha:", alpha_count, alpha)
    file_name_3 = "_alpha_" + "{:06.4f}".format(alpha)
    file_name_3 = file_name_3.replace(".", "p")

    def LV(t, ab):
        a, b = ab
        cosomt = math.cos(omega*t)
        kappa = kappa0 + kappa1*cosomt
        lamda1 = (1 - (1 + astar) * kappa) / astar
        lamda = (1-alpha) * lamda0 + alpha * lamda1
        if (lamda<0) or (kappa<0) or (omega0squared<0):
            print("Problem with parameters: lamda, kappa, omega0squared are:", lamda, kappa, omega0squared)
            error
    
        return [lamda * a*(b - 1), b*(1 - (a + b) * kappa) - lamda * a * b]
    
    def LV_log(t, ab):
        loga, logb = ab
        a = math.exp(loga)
        b = math.exp(logb)
        cosomt = math.cos(omega*t)
        kappa = kappa0 + kappa1*cosomt
        lamda1 = (1 - (1 + astar) * kappa) / astar
        lamda = (1-alpha) * lamda0 + alpha * lamda1
        if (lamda<0) or (kappa<0) or (omega0squared<0):
            print("Problem with parameters: lamda, kappa, omega0squared are:", lamda, kappa, omega0squared)
            error
    
        return [lamda * (b - 1), (1 - (a + b) * kappa) - lamda * a]

    # keep track of the number of transient periods we've run
    n_transients_run = 0
    
    # run the mean field for many periods to get over any transients
    if (verbose == 1):
        print("Running transient: ", n_forcing_periods_transient, " periods of the forcing")
    t_eval = [2 * math.pi / omega * nn for nn in [0, n_forcing_periods_transient-5, n_forcing_periods_transient-4, n_forcing_periods_transient-3, n_forcing_periods_transient-2, n_forcing_periods_transient-1]]
    sol = solve_ivp(LV, [0,n_forcing_periods_transient * 2 * math.pi / omega], [a_ic, b_ic], t_eval=t_eval, method='RK45', rtol=1e-10, atol=1e-10)
    n_transients_run = n_transients_run + n_forcing_periods_transient
    if sol.status != 0:
        print("Status: ", sol.status)
    mf_stroboscopic_a = sol.y[0]
    mf_stroboscopic_b = sol.y[1]
    a_ic = mf_stroboscopic_a[-1]
    b_ic = mf_stroboscopic_b[-1]

    # check to see if we might be close to a period one or two or four orbit, if so, run longer
    check_near_periodic = min([abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-2]), \
                               abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-3]), \
                               abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-5])])
    if (verbose == 1):
        print(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-2], \
              mf_stroboscopic_a[-1] - mf_stroboscopic_a[-3], \
              mf_stroboscopic_a[-1] - mf_stroboscopic_a[-5], check_near_periodic)
    next_tol = tol
    if ((check_near_periodic < 1000*tol) and (check_near_periodic > tol)):
        next_tol = check_near_periodic
        if (verbose == 1):
            print("Running extra transient: ", mf_stroboscopic_a[-2]-mf_stroboscopic_a[-1], \
                                               mf_stroboscopic_a[-3]-mf_stroboscopic_a[-1], \
                                               mf_stroboscopic_a[-5]-mf_stroboscopic_a[-1], \
                                               check_near_periodic);
        n_extra_transients = 2 * n_forcing_periods_transient
        t_eval = [2 * math.pi / omega * nn for nn in [0, n_extra_transients-5, n_extra_transients-4, n_extra_transients-3, n_extra_transients-2, n_extra_transients-1]]
        sol = solve_ivp(LV, [0, n_extra_transients * 2 * math.pi / omega], [a_ic, b_ic], t_eval=t_eval, method='RK45', rtol=1e-10, atol=1e-10)
        n_transients_run = n_transients_run + n_extra_transients
        if sol.status != 0:
            print("Status: ", sol.status)
        mf_stroboscopic_a = sol.y[0]
        mf_stroboscopic_b = sol.y[1]
        a_ic = mf_stroboscopic_a[-1]
        b_ic = mf_stroboscopic_b[-1]
       
    # if we are getting closer, run for a bit longer
    check_near_periodic = min([abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-2]), \
                               abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-3]), \
                               abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-5])])
    if (verbose == 1):
        print(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-2], \
              mf_stroboscopic_a[-1] - mf_stroboscopic_a[-3], \
              mf_stroboscopic_a[-1] - mf_stroboscopic_a[-5], check_near_periodic)
    if ((check_near_periodic < next_tol) and (check_near_periodic > tol)):
        next_tol = check_near_periodic
        if (verbose == 1):
            print("Running extra transient: ", mf_stroboscopic_a[-2]-mf_stroboscopic_a[-1], \
                                               mf_stroboscopic_a[-3]-mf_stroboscopic_a[-1], \
                                               mf_stroboscopic_a[-5]-mf_stroboscopic_a[-1], \
                                               check_near_periodic);
        n_extra_transients = 4 * n_forcing_periods_transient
        t_eval = [2 * math.pi / omega * nn for nn in [0, n_extra_transients-5, n_extra_transients-4, n_extra_transients-3, n_extra_transients-2, n_extra_transients-1]]
        sol = solve_ivp(LV, [0, n_extra_transients * 2 * math.pi / omega], [a_ic, b_ic], t_eval=t_eval, method='RK45', rtol=1e-10, atol=1e-10)
        n_transients_run = n_transients_run + n_extra_transients
        if sol.status != 0:
            print("Status: ", sol.status)
        mf_stroboscopic_a = sol.y[0]
        mf_stroboscopic_b = sol.y[1]
        a_ic = mf_stroboscopic_a[-1]
        b_ic = mf_stroboscopic_b[-1]

    # if we are still getting closer, run for a bit longer
    check_near_periodic = min([abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-2]), \
                               abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-3]), \
                               abs(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-5])])
    if (verbose == 1):
        print(mf_stroboscopic_a[-1] - mf_stroboscopic_a[-2], \
              mf_stroboscopic_a[-1] - mf_stroboscopic_a[-3], \
              mf_stroboscopic_a[-1] - mf_stroboscopic_a[-5], check_near_periodic)
    if ((check_near_periodic < next_tol) and (check_near_periodic > tol)):
        next_tol = check_near_periodic
        if (verbose == 1):
            print("Running extra transient: ", mf_stroboscopic_a[-2]-mf_stroboscopic_a[-1], \
                                               mf_stroboscopic_a[-3]-mf_stroboscopic_a[-1], \
                                               mf_stroboscopic_a[-5]-mf_stroboscopic_a[-1], \
                                               check_near_periodic);
        n_extra_transients = 8 * n_forcing_periods_transient
        t_eval = [2 * math.pi / omega * nn for nn in [0, n_extra_transients-5, n_extra_transients-4, n_extra_transients-3, n_extra_transients-2, n_extra_transients-1]]
        sol = solve_ivp(LV, [0, n_extra_transients * 2 * math.pi / omega], [a_ic, b_ic], t_eval=t_eval, method='RK45', rtol=1e-10, atol=1e-10)
        n_transients_run = n_transients_run + n_extra_transients
        if sol.status != 0:
            print("Status: ", sol.status)
        mf_stroboscopic_a = sol.y[0]
        mf_stroboscopic_b = sol.y[1]
        a_ic = mf_stroboscopic_a[-1]
        b_ic = mf_stroboscopic_b[-1]

    a_ic = mf_stroboscopic_a[-1] # keep these for the next alpha
    b_ic = mf_stroboscopic_b[-1]
    
    # run the mean field for n_forcing_periods_keep periods
    # we will look for repeats in the stroboscopic map, going forward a further n_forcing_periods_keep steps
    if (verbose == 1):
        print("Running data to keep: ", n_forcing_periods_keep, " periods of the forcing")
    t_eval = 2 * math.pi / omega * np.arange(n_forcing_periods_keep)
    sol = solve_ivp(LV, [0,n_forcing_periods_keep * 2 * math.pi / omega], [a_ic, b_ic], t_eval=t_eval, method='RK45', rtol=1e-10, atol=1e-10, max_step=0.1)
    if sol.status != 0:
        print("Status: ", sol.status)
    mf_stroboscopic_time = sol.t
    mf_stroboscopic_a = sol.y[0]
    mf_stroboscopic_b = sol.y[1]
    
    period = 0
    
    # hunt back to see if there is periodicity
    if (abs(a_ic - astar) > 100 * tol) and (abs(b_ic - 1) > 100 * tol):
        if (verbose == 1):
            print("Not an eqm point, hunt for the period")
        period = -1
    
        for i_hunt_back in range(n_forcing_periods_keep-1):
            if (verbose == 1):
                print("Checking 0 and ",i_hunt_back + 1," : ", mf_stroboscopic_a[0], mf_stroboscopic_a[i_hunt_back + 1], abs(mf_stroboscopic_a[0] - mf_stroboscopic_a[i_hunt_back + 1]))
            if (abs(mf_stroboscopic_a[0] - mf_stroboscopic_a[i_hunt_back + 1]) < tol):
                period = i_hunt_back + 1
                break
    
    # after the break, we have period = 1, 2, ... if it is periodic
    # period = -1: chaotic 
    # period =  0: we shouldn't be in this if block

    # print alpha, the values of a and b and the period
    if (period == 0) or (period == 1):
        print("  ", "{:06.4f}".format(alpha), "{:20.12g}".format(mf_stroboscopic_a[0]), "{:20.12g}".format(mf_stroboscopic_b[0]), "{:5d}".format(0), "  {:5d}".format(period), " {:10d}".format(n_transients_run)) 
        print("  ", "{:06.4f}".format(alpha), "{:20.12g}".format(mf_stroboscopic_a[0]), "{:20.12g}".format(mf_stroboscopic_b[0]), "{:5d}".format(0), "  {:5d}".format(period), " {:10d}".format(n_transients_run), file = LV_mean_field_alpha_output_file) 
    if (period > 1) or (period == -1):
        n_to_print = period
        if (n_to_print == -1):
            n_to_print = n_forcing_periods_keep
        for i_period in range(n_to_print):
            print("  ", "{:06.4f}".format(alpha), "{:20.12g}".format(mf_stroboscopic_a[i_period]), "{:20.12g}".format(mf_stroboscopic_b[i_period]), "{:5d}".format(i_period), "  {:5d}".format(period), " {:10d}".format(n_transients_run)) 
            print("  ", "{:06.4f}".format(alpha), "{:20.12g}".format(mf_stroboscopic_a[i_period]), "{:20.12g}".format(mf_stroboscopic_b[i_period]), "{:5d}".format(i_period), "  {:5d}".format(period), " {:10d}".format(n_transients_run), file = LV_mean_field_alpha_output_file) 

LV_mean_field_alpha_output_file.close()
    
#    # save a phase portrait and the data, but only if not an equilibrium point
#    # and save the phase portrait only if the period > 2
#    
#    # generate and save trajectory data for n_forcing_periods_keep steps, or period steps if it is periodic
#        t_max = t_eval[-1]
#        if period > 0:
#            t_max = period * 2 * math.pi / omega
#        # it shouldn't matter which initial conditions we choose, if we are in the attractor
#        sol = solve_ivp(LV, [0,t_max], [a_ic, a_ic], method='RK45', rtol=1e-10, atol=1e-10, max_step=0.1)
#        time = sol.t
#        a = sol.y[0]
#        b = sol.y[1]
#    
#        n_points = len(time)
#        timeab = np.zeros((n_points, 3))
#        timeab[:,0] = time
#        timeab[:,1] = a
#        timeab[:,2] = b
#        np.savetxt(file_name_1 + file_name_2 + "_timeab.txt.gz", timeab)
#    
#    # extract and save stroboscopic data for n_forcing_periods_keep steps, or period steps if it is periodic
#        n_points_stroboscopic = len(mf_stroboscopic_time)
#        if period > 0:
#            n_points_stroboscopic = period
#        timeab_stroboscopic = np.zeros((n_points_stroboscopic, 3))
#        timeab_stroboscopic[:,0] = mf_stroboscopic_time[0 : n_points_stroboscopic]
#        timeab_stroboscopic[:,1] = mf_stroboscopic_a   [0 : n_points_stroboscopic]
#        timeab_stroboscopic[:,2] = mf_stroboscopic_b   [0 : n_points_stroboscopic]
#        np.savetxt(file_name_1 + file_name_2 + "_timeab_stroboscopic.txt.gz", timeab_stroboscopic)
#    
#        print(mf_stroboscopic_a[0:period])
#    
#    # plot a phase portrait, but only if the period > 2, or n = 0.7 or we are doing only one run
#        if (period > 2) or (period < 0) or (abs(n - 0.7) < 1.0e-6) or (one_run_only == 1):
#            fig_phaseportrait, fig_phaseportrait_axes = plt.subplots(nrows = 1, ncols = 1, figsize=(10, 10), tight_layout=True)
#            fig_phaseportrait_axes.set_xlabel("$a$")
#            fig_phaseportrait_axes.set_ylabel("$b$")
#            fig_phaseportrait_axes.plot(a, b, 'k', linewidth=2)
#            if period > 0:
#                fig_phaseportrait_axes.plot(mf_stroboscopic_a[0:period], mf_stroboscopic_b[0:period], 'm+')
#            if period == -1:
#                fig_phaseportrait_axes.plot(mf_stroboscopic_a, mf_stroboscopic_b, 'm+', linewidth=2)
#            fig_phaseportrait.tight_layout()
#            if save_figs_as_pdfs == 1:
#                fig_phaseportrait.savefig(file_name_1 + file_name_2 + "_phase_portrait" + pdf_png_extension, bbox_inches='tight')
#            else:
#                plt.show(block = True)
#            plt.close(fig_phaseportrait)
#    
#    # exit the "if not at a fixed point" block
#    print("Summary: ", file_name_2, " ", \
#          "{:06.4f}".format(astar), "{:06.4f}".format(kappa0), "{:06.4f}".format(kappa1), "{:06.4f}".format(n), "{:06.4f}".format(alpha), \
#          "{:06.4f}".format(lamda0), "{:06.4f}".format(omega0), "{:06.4f}".format(omega), \
#          "{:06.4f}".format(a_ic), "{:06.4f}".format(b_ic), \
#          "{:5.2f}".format(version), \
#          "{:6d}".format(period))
#    
