########## 
##########
# Get all initial populations before perturbation 
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-08-09
##########
##########

# set up =======================================================================
import Pkg
Pkg.instantiate()

using DrWatson
@quickactivate "trait-based-rewiring"

using DifferentialEquations, DelimitedFiles

# logger set up 
#using Logging: global_logger
#using TerminalLoggers: TerminalLogger
#global_logger(TerminalLogger())

include(srcdir("00_species_list_creation.jl"))
include(srcdir("01_eqn_component_functions.jl"))
include(srcdir("03_initial_sims_functions.jl"))
include(srcdir("04_perturbation_sims_functions.jl"))

# set all positional values 
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12 

# read in files 
all_species_lists = readdlm(datadir("intermediate-files/species_lists_for_simulations.csv"), ',', Any)

# little line to set the network id (this happens from the shell line)
network_id = parse(Int64, ARGS[1]) 
 
# set the species list to be worked on 
species_list = all_species_lists[all_species_lists[:,network_id_pos] .== network_id,:]
network = create_interaction_matrix(species_list)

# set up the simulation 
#  allocate u0
u0 = Float64[]
for i in 1:size(species_list, 1)
    append!(u0, species_list[i, biomass_pos])
end

# set the extinction callback
# set vector callback
extinction_callback = 
    VectorContinuousCallback(extinction_threshold_pre_perturb,
                            extinction_affect_pre_perturb!,
                            size(species_list, 1),
                            save_positions = ntuple(x->true,
                                                    size(species_list, 1)))

# put species list and network in same object
p = [species_list, network]

# set  tspan
tspan = (0.0, 1000.0)

# set problem 
prob_setup = ODEProblem(biomass_basic,
                            u0, 
                            tspan, 
                            p) 
@time sol = solve(prob_setup, 
                    Tsit5(),
                    callback = extinction_callback,
                    save_times = [0.0, 1000.0])

# put the values back in the species list
species_list[:, biomass_pos] = sol.u[end]

# change extinction values
for i in 1:size(species_list, 1)
    if species_list[i, biomass_pos] == 0
        species_list[i, alive_pos] = 0
    end 
end

name = lpad(network_id, 4, "0")

# bind back the network_id and write the species list
writedlm(datadir("./intermediate-files/species-lists-pre-perturbation/species_list_$network_id.csv"), species_list, ",")
