# Script to produce data for figures 2(a-c) and 4(d-f)
# Created by Matthew Asker 

import numpy as np
import sys
from next_reaction_method import next_reaction

a = 0.5
s = 0.2
Nth_values = [20, 25, 30]
total_sims = 10
x_values = [0.5]
K_0 = 560
gamma = 11/14
points = 21
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(Nth_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 N_th in range(0, len(Nth_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), Nth_values[N_th], total_sims)
                G.run()
                rgb = np.array([(G.final_N_C, G.final_N_D, G.average_time)])[0]
                coexist_array[x][N_th][delta][nu] = rgb

#index = sys.argv[1]
filename = f'heatmap K_0={K_0}, gamma={round(gamma,2)}, Nth_values={Nth_values}, sims={total_sims}'#' {index}'
results = [[Nth_values, points, switch_rates, bias], coexist_array]                
np.save(filename, results)