# functions to create parameters

dx_i/dt = r_i*X_i + sum(F_i(X_j)) - sum(F_j(X_i)) - m_i*X_i






dx_i/dt = fun1(X_i) + fun2(X_i) - fun3(X_i) - fun4(X_i)

function lorenz_added!(dx, x, p, t)
    dx[1] = sigma*(x[2]-x[3])
    dx[2] = x[2]*(rho - x[3]) - x[2]
    dx[3] = x[1]*x[2] - beta*x[3]
end

# dx is just a vector
# populate the vector via a loop through the functions where
# you look at each species 
for i in species
    # get abundance from somewhere
    dx[i] = fun1(X_i) + fun2(X_i) - fun3(X_i) - fun4(X_i)
end 

function sigma_fun(x::Array, sigma::Float)

    return (sigma*(x[2] - x[1]))
    
end


x0 = [1.0;0.0;0.0]



# dx1 = r1*x1 + sum(g12*x2, g13*x3) - sum(g21*x1, g31*x1) - m1*x1

x_metrics = [[2, 2], [1, 3], [1, 2.5]]
x_vals = [1, 2, 3]
x0 = [100, 200, 75]

"""
Function defines the density-independent growth of xi. 
"""
function growth(r, i)
    return r*x[i]
end

""" 
Determines for any focal x (i) and any other x (j) if the focal x will gain 
biomass from j based on some arbitrary metric (metric defined in the x_metrics
vector), and if xi does gain from xj, then gives a g value (rate at which the
gain happens)
""" 
function get_g(i, j)
    if x_metrics[i][2] > x_metrics[j][2]
        g = 0.5
    else 
        g = 0
    end 
end 

""" 
Figures out how much biomass is actually taken from any xj and given to a given 
xi.
""" 
function add(i, j)
    return get_g(i, j)*x[j]
end

"""
For a specific xi, gets all other x's that it can consume biomass from, gets the 
biomass from each, and adds it to xi's biomass.  
"""
function add_all(i)
    other_x = filter(x-> x ≠ i, x_vals)
    vals = []
    for j in other_x
        append!(vals, add(i, j))
    end 

    return sum(vals)
end 

"""
For any given xi, finds the other xj's that can remove biomass from xi, gets how
much biomass each xj will take from xi, and removes all of those values from xi. 
"""
function remove_all(i)
    other_x = filter(x -> x ≠ i, x_vals)
    vals = []
    for j in other_x
        append!(vals, add(j, i))
    end 

    return sum(vals)
end 


function death(i)
    m = 1
    m*x[i]
end 
death(i)
function eqns!(dx, x, p, t)
    dx = [[],[],[]]
    for i in x_vals 
        dx[i] = growth(all_x[i][1], i) + add_all(i) - remove_all(i) - death(i)
    end 
end

using DifferentialEquations
x0 = x
tspan = (0.0, 100.0)
prob = ODEProblem(eqns!, x0, tspan)
sol = solve(prob)




# dx1/dt = r1*x1 + x1(a12*x2) - a21*x1 - m1*x1 
# dx2/dt = r2*x2 + x2(a21*x1) - a12*x2 - m2*x2 

#dx/dt = indep_growth + dep_growth - dep_death - indep_death 


function indep_growth(i, x) 
    #define r
    if i == 1
        r = 1.5
        k = 5000
    else 
        r = 1.9
        k = 2000
    end 
    # return product of r and abundance of focal x 
    return (r*x[i]*(1-(x[i]/k)))
end 

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

function dep_death(i, x)
    if i == 1
        a = 0.01
        return x[i]*a*x[2]
    else 
        a = 0.010975
        return x[i]*a*x[1]
    end 
    return death 
end 

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

    return m*x[i]
end 

dx = zeros(2)
function eqns!(dx, x, p, t)
    for i in 1:2
        dx[i] = indep_growth(i, x) + dep_growth(i, x) - dep_death(i, x) - indep_death(i, x)
    end 
end 

r1 = 1.5
r2 = 1.9
a12 = 0.010975
a21 = 0.01
m1 = 1.4
m2 = 1.3
k1 = 5000
k2 = 2000
function eqns!(dx, x, p, t)
    dx[1] = r1*x[1]*(1-(x[1]/k1)) + x[1]*(a12*x[2])*(1-(x[1]/k1)) - 
            x[1]*(a21*x[2]) - m1*x[1]
    dx[2] = r2*x[2]*(1-(x[2]/k2)) + x[2]*(a21*x[1])*(1-(x[2]/k2)) - 
            x[2]*a12*x[1] - m2*x[2]
end

x0 = [100.0; 1000.0]
tspan = (0.0, 10.0)
prob = ODEProblem(eqns!, x0, tspan)
sol2 = solve(prob)
plot(sol2)

sol1[:] == sol2[:]