########## 
##########
# Simulate to get equilibrium abundances for the webs
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-06-25
##########
##########

# set up =======================================================================

##include("./src/eqn_component_functions.jl")
#include("./src/species_list_creation.jl")
#include("./src/network_metrics.jl")
using DifferentialEquations
using DelimitedFiles

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

"""

    biomass(dx, x, p, t)

Calling biomass functions, create the actual equation for each species, to be 
passed to the simulation. 
"""
function biomass_basic(du, u, p, t) 
    species_list, network = p[1:2]
    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)
    end 
end 
 

"""
    starting_eq_sims(species_list, network, initial_sim_length)

Function runs through all initial simulations from t=0 to the point of starting 
measurement. 
"""
function starting_eq_sims!(
    species_list,
    network,
    initial_sim_length, 
    network_id
    )

    #initials = species_list, network

    # cut network to size 
    network = network[:,2:end]

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

    # keep track of how many extinctions there were 
    extinct_sp = Float64[]

    # loop through all time values from 0 to 5000
    time_start = 0.0
    while time_start < initial_sim_length

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

        # set time and defie ODE problem 
        tspan = (time_start, time_start+9.9)
        prob = ODEProblem(biomass!, x0, tspan, p)
        sol = solve(prob)

        # keep only last timestep and check for extinction
        x0 = sol.u[end]
        check_extinctions!(x0)
        species_list[:,biomass_pos] = x0 

        # get new network
        network = create_interaction_matrix(species_list)

        # change alive to extinct if species is extinct 
        for i in 1:size(species_list, 1)
            if species_list[i, biomass_pos] == 0
                species_list[i, alive_pos] = 0
                append!(extinct_sp, species_list[i, organism_id_pos])
            end 
        end 

        time_start += 10.0
    end 

    # bind back the network_id and write the species list
    species_list = hcat(species_list, fill(network_id, 50, 1))
    writedlm("./data/intermediate-files/species-lists-pre-perturbation/species_list_$network_id.csv", species_list, ",")

    # create and write network 
    network = create_interaction_matrix(species_list)
    network = hcat(fill(network_id, 50, 1), network) 
    writedlm("./data/intermediate-files/networks-pre-perturbation/network_$network_id.csv", network, ",")
                   
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_pre_perturb(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-06 - 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_pre_perturb!(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