# Script to produce data for figures 4a-d and 5
# Created by Matthew Asker 

import numpy as np
from next_reaction_method_paper import next_reaction

a = 0.25
s = 0.1
xth_values = [0.1, 0.37, 0.63, 0.9]
total_sims = 1
x_values = [0.5]
K_0 = 550
gamma = 9/11
N_values = np.arange(round(K_0*(1-gamma)), round(K_0*(1+gamma)), 1)
points = 101
switch_rates = np.power(10, np.linspace(np.log10(0.0001), np.log10(100), points))
bias = np.linspace(-1, 1, points)

coexist_array = np.full([len(x_values), len(xth_values), len(bias), len(switch_rates), 3], np.array([(0.0, 0.0, 0.0)])[0])

for x in range(0, len(x_values)):
    for x_th in range(0, len(xth_values)):
        for delta in range(0, len(bias)):
            for nu in range(0, len(switch_rates)):
                G = next_reaction(s, a, switch_rates[nu]*(1-bias[delta]), switch_rates[nu]*(1+bias[delta]), K_0*(1+gamma), K_0*(1-gamma), xth_values[x_th], total_sims)
                G.run()
                rgb = np.array([(G.final_N_R, G.final_N_S, G.average_time)])[0]
                coexist_array[x][x_th][delta][nu] = rgb

filename = f'heatmap K_0={K_0}, gamma={round(gamma,2)}, xth_values={xth_values}, sims={total_sims}'
results = [[xth_values, points, switch_rates, bias], coexist_array]                
np.save(filename, results)