BecksLab / EcologicalNetworksDynamics.jl

A simulator for ecological dynamics written in Julia.
GNU Affero General Public License v3.0
18 stars 3 forks source link

Make the package more modular: Definition of a Params abstract type #70

Closed alaindanet closed 1 year ago

alaindanet commented 2 years ago

In order to let the user defines its own version of ModelParameters, we define Params as an abstract type, ModelParameters having now a Params subtype.

All the functions of BEFWM2 that took ModelParameters type as an argument nows take Params type of argument.

I present below a use case where I add an additional mortality term to the model:

using Revise
using BEFWM2

######
# Define your own ModelParameters with an additional mortality term 
#####

# Define the composite type CritDeath
mutable struct CritDeath 
    d::Float64
end

# Add a field crit_death to ModelParameters by defining a new composite type of subtype Params
mutable struct CustomModelParameters <: Params
    network::EcologicalNetwork
    biorates::BioRates
    environment::Environment
    functional_response::FunctionalResponse
    crit_death::CritDeath
end

# Define the function of my custom ModelParameters
function CustomModelParameters(
    network::EcologicalNetwork;
    biorates::BioRates = BioRates(network),
    environment::Environment = Environment(network),
    functional_response::FunctionalResponse = BioenergeticResponse(network),
    crit_death::CritDeath = CritDeath(0.2)
)
    CustomModelParameters(network, biorates, environment, functional_response, crit_death)
end

# My custom dBdt! contains B[i] * params.crit_death.d in metabolic losses
# It takes a object of type Params in input but consumption() should also take a Params input instead of  #ModelParameters type of input
function CustomdBdt!(dB, B, params::Params, t)

    # Set up - Unpack parameters
    S = richness(params.network)
    response_matrix = params.functional_response(B, params.network)
    r = params.biorates.r # vector of intrinsic growth rates
    K = params.environment.K # vector of carrying capacities
    network = params.network

    # Compute ODE terms for each species
    for i in 1:S
        growth = logisticgrowth(i, B, r[i], K[i], network)
        eating, being_eaten = consumption(i, B, params, response_matrix)
        metabolism_loss = metabolic_loss(i, B, params) + B[i] * params.crit_death.d
        net_growth_rate = growth + eating - metabolism_loss
        net_growth_rate = effect_competition(net_growth_rate, i, B, network)
        dB[i] = net_growth_rate - being_eaten
    end
end

# Do the simulation
A = [0 0 0 0; 1 0 0 0; 1 0 0 0; 0 1 1 0] 
foodweb = FoodWeb(A)

params = CustomModelParameters(foodweb, crit_death = CritDeath(0.2))
m = simulate(params, rand(4), diff_code_data=(CustomdBdt!, params))
iago-lito commented 1 year ago

Closing based on https://github.com/BecksLab/BEFWM2/pull/71#issuecomment-1305687226.