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

using DifferentialEquations
using DelimitedFiles

# functions ====================================================================

"""
    biomass_perturb_all(du, u, p, t)

Calling biomass functions, create the actual equation for each species, to be 
passed to the simulation. Do this by perturbing all species. 
"""
function biomass_perturb_all(du, u, p, t) 
    # pick set of random 
    # extract all values to be used w/in function from the p value
    species_list, network, base_mortality, z = p
    for i in 1:size(species_list, 1)
        du[i] = growth_rate(i, species_list, u) + 
                biomass_gain(species_list, network, u, 2, i) -
                biomass_loss(species_list, network, u, 2, i) -
                mortality_change(i, u, base_mortality, z)
    end 
end 

"""
    biomass_perturb_random(du, u, p, t)

Calling biomass functions, create the actual equation for each species, to be 
passed to the simulation. Do this by randomly picking ~1/3 of the species and 
changing the mortality rate. The species to be perturbed are defined before this
function is used, as are the ones NOT to be perturbed. 
"""
function biomass_perturb_random!(du, u, p, t) 
    # extract all values to be used w/in function from the p value
    # extract all values to be used w/in function from the p value
    species_list, network, base_mortality, amount, perturb_sp, non_perturb = p
    if amount === nothing
        for i in 1:size(species_list, 1)
            du[i] = growth_rate(i, species_list, u) + 
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, amount)
        end 
    else 
        for i in perturb_sp
            du[i] = growth_rate(i, species_list, u) + 
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, amount) 
        end 
        for i in non_perturb
            du[i] = growth_rate(i, species_list, u) + 
                    biomass_gain(species_list, network, u, 2, i) -
                    biomass_loss(species_list, network, u, 2, i) -
                    mortality_change(i, u, base_mortality, amount)
        end
    end 
end

"""
    extinction_cond(out, x, t, integrator)

Condition asking if any of the values in the integrator go below the extinction
threshold. 
"""
function extinction_threshold(out,u,t,integrator)
    # loop through all species to make the condition check all of them
    for i in 1:size(species_list, 1)
        out[i] = 1e-05 - u[i]
    end
end 

"""
    extinction_affect!(integrator, event_idx)

Affect function to use in callback vector, this will allow us to check 
if any of the values of the threshold condition are true, and if so, 
set the value to zero in the integrator.
"""
function extinction_affect!(integrator, event_idx)
    # loop again through all species
    for i in 1:size(species_list, 1)
        if event_idx == i
            integrator.u[i] = 0 
        end
    end 
end

"""
    inc_mort_cond(u,t,integrator)

Condition to find the stop value and when it's hit, increase the parameter value
specified. Do this whenever a value of stop1 is reached. 
"""
function inc_mort_cond(u,t,integrator)
  t in stop1
end

"""
    dec_mort_cond(u,t,integrator)

Condition to find the stop value and when it's hit, decrease the parameter value
specified. Do this whenever a value of stop2 is reached. 
"""
function dec_mort_cond(u,t,integrator)
  t in stop2
end

"""
    affect_inc!(integrator)

Affect for the condition to work on. Increase the mortality by making the amount
(the fourth element in p) equal to z.  
"""
function affect_inc!(integrator, z)
  integrator.p[4] = z
end

"""
    affect_dec!(integrator)

Affect for the condition to work on. Decrease the fourth element in p, 
returning it to 1.0.  
"""
function affect_dec!(integrator)
  integrator.p[4] = nothing 
end
