########## 
##########
# Functions to perform perturbations on the food webs 
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-08-09
##########
##########

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

include("./04_eqn_component_functions.jl")
include("./01_species_list_creation.jl")
#include("./05_network_metrics.jl")
include("./08_perturbation_sims_functions.jl")

#@info{"All prereqs done..."}

# 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 

# pulse perturbations ==========================================================

all_species_lists = readdlm("./data/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]) 

species_list = all_species_lists[all_species_lists[:,network_id_pos] .== network_id,:]
network = create_interaction_matrix(species_list)

# set parameters - pulse
c = 2 # constant amount
x = 4000 # duration
f = 1 # frequency
z = c/(x*(1/f))

# set vector callback
extinction_callback = 
    VectorContinuousCallback(extinction_threshold,
                            extinction_affect!,
                            size(species_list, 1),
                            save_positions = ntuple(x->true,
                                                    size(species_list, 1)))
steady_callback = TerminateSteadyState(1e-12, 
                                        1e-9, 
                                        DiffEqCallbacks.allDerivPass)

# set callback set
callback_set_perturb = extinction_callback
callback_set = CallbackSet(extinction_callback, steady_callback)

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

# set perturbation tspan
tspan = (1000.0, 5000.0)

# set constant stop values
const stop1 = [1000.0]
const stop2 = [5000.0]
const tstop = [1000.;5000.]

# loop through random selections of species to perturb
#@info{"starting loop..."}

for j in 1:10

    # pick set of random species to perturb
    all = Array(1:size(species_list, 1))
    perturb_sp = rand(Array(1:size(species_list, 1)),
                        Int(round(size(species_list, 1)*(1/3))))
    non_perturb = setdiff(Set(all), Set(perturb_sp))

    # set times to save from the simulation 
    save_times = [collect(range(1000.0, 1250.0, step = 0.1));
                collect(range(1250, 4800, step = sample(200:400, 1)[1]));
                collect(range(4800.0, 5000.0, step = 0.1))]

    # set p 
    p = [species_list, network, 0.3, z, perturb_sp, non_perturb]

    prob_press = ODEProblem(biomass_perturb_random,
                                        u0, 
                                        tspan, 
                                        p, 
                                        tstops=tstop) 
    @time sol_prob_press = solve(prob_press, 
                                        Tsit5(),
                                        callback = callback_set,
                                        # save every part for the perturbation and 
                                        # the first 150 time steps after the 
                                        # perturbation, but then just grab a few 
                                        # from there to the end 
                                        saveat = save_times)
    # save output of the solver
    output = hcat(sol_prob_press.t, sol_prob_press.u)
    # write out the biomass outputs
    name = string("press-", network_id, "-", j, "-", "sims.csv") 
    path = "./data/simulation-outputs/press-sims/"
    writedlm(string(path,name), output, ",")
end
