import Pkg
Pkg.add("DifferentialEquations")
using DifferentialEquations, Plots
function indep_growth(u, i) 
    #define r
    if i == 1
        r = 1.5
        k = 5000
    else 
        r = 1.8
        k = 2000
    end 
    # return product of r and abundance of focal u 
    return (r*u[i]*(1-(u[i]/k)))
end 

function dep_growth(u, i, a_pass)
    # determine if a12 or a21 should be used 
    if i == 1 # use a12 
        #a = 0.01
        a = a_pass
        k = 5000
        return u[i]*(a*u[2])*(1-(u[i]/k))
    else 
        a = 0.01
        k = 2000
        return u[i]*(a*u[1])*(1-(u[i]/k))
    end 
end 

function dep_death(u, i, a_pass)
    if i == 1
        death = dep_growth(u, 2, a_pass)
    else 
        death = dep_growth(u, 1, a_pass)
    end 
    return death 
end 

function indep_death(u, i)
    if i == 1
        m = 1.4
    else 
        m = 1.3
    end 

    return m*u[i]
end 


function eqns(du, u, p, t)
    a_pass = p
    for i in 1:1
        du[i] = indep_growth(u, i) + 
        dep_growth(u, i, a_pass) - 
        dep_death(u, i, a_pass) - 
        indep_death(u, i)
    end
    for i in 2:2
        du[i] = indep_growth(u, i) + 
        dep_growth(u, i, a_pass) - 
        dep_death(u, i, a_pass) - 
        indep_death(u, i)
    end
end 

function extinction_threshold(out,u,t,integrator)
    for i in 1:2
        out[i] = 10 - u[i]
    end
end 

function extinction_affect!(integrator, event_idx)
    # loop again through all species
    for i in 1:2
        if event_idx == i
            integrator.u[i] = 0 
        end
    end 
end 

# function eqns_test!(du, u, p, t)
#     for i in 1:2
#         du[i] = indep_growth(u, i) + 
#         dep_growth(u, i) - 
#         dep_death(u, i) - 
#         dep_death(u, i)
#     end 
# end 

# set constant stop values
const stop1 = [5.0, 60.0]
const stop2 = [8.0, 75.0]

function inc_mort_cond(u,t,integrator)
  t in stop1
end

# decrease mortality (stop pulse)
function dec_mort_cond(u,t,integrator)
  t in stop2
end

# apply an affect when condition is reached
function affect_inc!(integrator)
  integrator.p += 0.001
end

function affect_dec!(integrator)
  integrator.p -= 0.001
end
# set vector callback
extinction_callback = 
    VectorContinuousCallback(extinction_threshold,
                             extinction_affect!,
                             2,
                             save_positions = ntuple(x->true,2))

save_positions = (true, true)
cb_inc = DiscreteCallback(inc_mort_cond, 
                            affect_inc!, 
                            save_positions=save_positions)
save_positions = (true, true)
cb_dec = DiscreteCallback(dec_mort_cond, 
                            affect_inc!, 
                            save_positions=save_positions)
steady_callback = TerminateSteadyState(1e-8, 
                            1e-6, 
                            DiffEqCallbacks.allDerivPass)
                             
callback_set = CallbackSet(cb_inc, cb_dec, extinction_callback, steady_callback)

const tstop = [5.;8.;60.;75.;]
save_values = [range(5.0, 20.0, step = 0.1); 
                range(20.0, 59.9, step = 5); 
                range(60, 100, step = 1)]
save_values = [range(0,100, step = 10)]

network_id = parse(Int64, ARGS[1]) 
p = 0.01
u0 = [100.0; 100.0]
tspan = (0.0, 1000.0)
prob = ODEProblem(eqns, 
                    u0, 
                    tspan,
                    p, 
                    tstops=tstop)
@time sol = solve(prob, 
                Tsit5(),
                reltol = 1e-08,
                callback = callback_set)
object = hcat(sol.t, sol.u)
writedlm("./out/out$network_id.csv", species_list, ",")
plot(sol)
