Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update model defaults within the new api. #132

Merged
merged 2 commits into from
Mar 15, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/Internals/model/simulate.jl
Original file line number Diff line number Diff line change
Expand Up @@ -185,7 +185,7 @@ function simulate(
problem,
alg;
callback = callback,
isoutofdomain = (u, p, t) -> any(x -> x < 0, u[species_indices(params)]),
isoutofdomain = (u, p, t) -> any(x -> x < 0, u),
kwargs...,
)
sol
Expand Down
5 changes: 4 additions & 1 deletion src/components/functional_responses.jl
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,10 @@ mutable struct ClassicResponse <: FunctionalResponse
ClassicResponse,
kwargs;
default = (
M = (Z = 1,),
# Don't bring BodyMass by default,
# because it has typically already been added before
# as it is useful to calculate numerous other parameters.
M = nothing,
e = :Miele2019,
h = 2,
w = :homogeneous,
Expand Down
70 changes: 59 additions & 11 deletions src/components/nutrients/nodes.jl
Original file line number Diff line number Diff line change
@@ -1,14 +1,27 @@
# Nutrients nodes compartments, akin to `Nodes`.
# Two possible blueprints because their number can be inferred from producer species.

# ==========================================================================================
# Call it 'Nodes' because the module is already named 'Nutrients'.
mutable struct Nodes <: ModelBlueprint
abstract type Nodes <: ModelBlueprint end
# Don't export, to encourage disambiguated access as `Nutrients.Nodes`.

Nodes(raw) =
if raw isa Symbol
NodesFromFoodweb(raw)
else
RawNodes(raw)
end

# Akin to Species.
mutable struct RawNodes <: Nodes
names::Vector{Symbol}
Nodes(names) = new(Symbol.(names))
Nodes(n::Integer) = new([Symbol(:n, i) for i in 1:n])
Nodes(names::Vector{Symbol}) = new(names)
RawNodes(names) = new(Symbol.(names))
RawNodes(n::Integer) = new([Symbol(:n, i) for i in 1:n])
RawNodes(names::Vector{Symbol}) = new(names)
end

function F.check(_, bp::Nodes)
function F.check(_, bp::RawNodes)
(; names) = bp
# Forbid duplicates (triangular check).
for (i, a) in enumerate(names)
Expand All @@ -19,15 +32,48 @@ function F.check(_, bp::Nodes)
end
end

function F.expand!(model, bp::Nodes)
function add_nutrients!(model, names)
# Store in the scratch, and only alias to model.producer_growth
# if the corresponding component is loaded.
model._scratch[:nutrients_names] = bp.names
model._scratch[:nutrients_index] = OrderedDict(n => i for (i, n) in enumerate(bp.names))
model._scratch[:nutrients_names] = names
model._scratch[:nutrients_index] = OrderedDict(n => i for (i, n) in enumerate(names))
end

F.expand!(model, bp::Nodes) = add_nutrients!(model, bp.names)

@component RawNodes

#-------------------------------------------------------------------------------------------
mutable struct NodesFromFoodweb <: Nodes
method::Symbol
function NodesFromFoodweb(method)
method = Symbol(method)
@check_symbol method (:one_per_producer,)
new(method)
end
end

function F.early_check(_, bp::NodesFromFoodweb)
(; method) = bp
@check_symbol method (:one_per_producer,)
end

@component Nodes
# Don't export to encourage disambiguated access as `Nutrients.Nodes`.
function F.expand!(model, bp::NodesFromFoodweb)
(; method) = bp
(; n_producers) = model
names = @build_from_symbol(
method,
:one_per_producer => [Symbol(:n, i) for i in 1:n_producers],
)
add_nutrients!(model, names)
end

@component NodesFromFoodweb requires(Foodweb)

#-------------------------------------------------------------------------------------------
@conflicts(RawNodes, NodesFromFoodweb)
# Temporary semantic fix before framework refactoring.
F.componentof(::Type{<:Nodes}) = Nodes

# ==========================================================================================
@expose_data graph begin
Expand Down Expand Up @@ -56,7 +102,9 @@ macro nutrients_index()
end

# ==========================================================================================
function F.display(model, ::Type{Nodes})
display_short(bp::Nodes; kwargs...) = display_short(bp, Nodes; kwargs...)
display_long(bp::Nodes; kwargs...) = display_long(bp, Nodes; kwargs...)
function F.display(model, ::Type{<:Nodes})
N = model.n_nutrients
names = model.nutrients_names
"Nutrients: $N ($(join_elided(names, ", ")))"
Expand Down
22 changes: 15 additions & 7 deletions src/components/producer_growth.jl
Original file line number Diff line number Diff line change
Expand Up @@ -50,20 +50,28 @@ mutable struct NutrientIntake <: ProducerGrowth
supply::Option{Nutrients.Supply}
concentration::Option{Nutrients.Concentration}
half_saturation::Option{Nutrients.HalfSaturation}
function NutrientIntake(nodes = nothing; kwargs...)
(!isnothing(nodes) && haskey(kwargs, :nodes)) &&
argerr("Nodes specified once as plain argument ($(repr(nodes))) \
and once as keyword argument (nodes = $(kwargs[:nodes])).")
function NutrientIntake(nodes = missing; kwargs...)
nodes = if haskey(kwargs, :nodes)
iago-lito marked this conversation as resolved.
Show resolved Hide resolved
ismissing(nodes) ||
argerr("Nodes specified once as plain argument ($(repr(nodes))) \
and once as keyword argument (nodes = $(kwargs[:nodes])).")
kwargs[:nodes]
elseif ismissing(nodes)
:one_per_producer
else
nodes
end
fields = fields_from_kwargs(
NutrientIntake,
kwargs;
# Values from Brose2008.
default = (;
r = :Miele2019,
nodes,
turnover = 0.25,
supply = 10,
concentration = 1,
half_saturation = 1,
supply = 4,
concentration = 0.5,
half_saturation = 0.15,
iago-lito marked this conversation as resolved.
Show resolved Hide resolved
),
)
new(fields...)
Expand Down
33 changes: 8 additions & 25 deletions src/default_model.jl
Original file line number Diff line number Diff line change
Expand Up @@ -104,25 +104,6 @@ function default_model(
# println(" collect from caller")
# Pick from caller input.
collect!(take_from_caller!(need))
elseif need <: BodyMass &&
!collected(need) &&
!excluded(need) &&
given(ClassicResponse) &&
any(B <: BodyMass for B in F.brings(input[ClassicResponse]))
# println(" special bodymass case")
# Very special hack-patch to pick body mass
# from a given classic response if given.
# TODO: The proper fix is to reify the tree structure
# of brought sub-blueprints, use dedicated data structures to query it,
# and allow general editions of the blueprints given here
# within the default model.
# This will be better done after the framework has been refactored.
remove_hook!(BodyMass)
bp = deepcopy(input[ClassicResponse])
bm = bp.M
bp.M = nothing
input[ClassicResponse] = bp
collect!(bm)
elseif !collected(need) && !excluded(need) && haskey(hooks, need)
# println(" collect from hooks")
collect!(pop!(hooks, need))
Expand Down Expand Up @@ -210,7 +191,7 @@ function default_model(
#---------------------------------------------------------------------------------------
# Not always required but often missing.

hooks[BodyMass] = BodyMass(; Z = 1)
hooks[BodyMass] = BodyMass(; Z = 10)
hooks[MetabolicClass] = MetabolicClass(:all_invertebrates)

#---------------------------------------------------------------------------------------
Expand All @@ -223,10 +204,12 @@ function default_model(
() -> if nutrients_given
NutrientIntake(;
r = tb!(GrowthRate, temperature_given ? :Binzer2016 : :Miele2019),
# Pick defaults from Brose2008.
nodes = tb!(N.Nodes, :one_per_producer),
turnover = tb!(N.Turnover, 0.25),
supply = tb!(N.Supply, 10),
concentration = tb!(N.Concentration, 1),
half_saturation = tb!(N.HalfSaturation, 1),
supply = tb!(N.Supply, 4),
concentration = tb!(N.Concentration, 0.5),
half_saturation = tb!(N.HalfSaturation, 0.15),
)
elseif temperature_given
LogisticGrowth(;
Expand All @@ -250,7 +233,7 @@ function default_model(
FunctionalResponse,
() -> if temperature_given
ClassicResponse(;
M = tb!(BodyMass, (; Z = 1)),
M = nothing,
e = tb!(Efficiency, :Miele2019),
h = tb!(HillExponent, 2),
w = tb!(ConsumersPreferences, :homogeneous),
Expand All @@ -260,7 +243,7 @@ function default_model(
)
elseif nti_given
ClassicResponse(;
M = tb!(BodyMass, (; Z = 1)),
M = nothing,
e = tb!(Efficiency, :Miele2019),
h = tb!(HillExponent, 2),
w = tb!(ConsumersPreferences, :homogeneous),
Expand Down
30 changes: 27 additions & 3 deletions src/methods/main.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,35 @@
# which is the reason they live after all components specifications.

# Major purpose of the whole model specification: simulate dynamics.
simulate(model, u0; kwargs...) = Internals.simulate(model, u0; kwargs...)
function simulate(model::InnerParms, u0, tmax::Integer; kwargs...)
# Depart from the legacy Internal defaults.
@kwargs_helpers kwargs

# No default simulation time anymore.
given(:tmax) && argerr("Received two values for 'tmax': $tmax and $(take!(:tmax)).")

# Lower threshold.
extinction_threshold = take_or!(:extinction_threshold, 1e-12, Any)
extinction_threshold = @tographdata extinction_threshold {Scalar, Vector}{Float64}

# Shoo.
verbose = take_or!(:verbose, false)
iago-lito marked this conversation as resolved.
Show resolved Hide resolved

# No TerminateSteadyState.
extc = extinction_callback(model, extinction_threshold; verbose)
callback = take_or!(:callbacks, Internals.CallbackSet(extc))

Internals.simulate(model, u0; tmax, extinction_threshold, callback, verbose, kwargs...)
end
@method simulate depends(FunctionalResponse, ProducerGrowth, Metabolism, Mortality)
export simulate

# Re-expose from internals so it works with the new API.
extinction_callback(m, thr; verbose=false) = Internals.ExtinctionCallback(thr, m, verbose)
@method extinction_callback depends(FunctionalResponse, ProducerGrowth, Metabolism, Mortality)
extinction_callback(m, thr; verbose = false) = Internals.ExtinctionCallback(thr, m, verbose)
iago-lito marked this conversation as resolved.
Show resolved Hide resolved
export extinction_callback
@method extinction_callback depends(
FunctionalResponse,
ProducerGrowth,
Metabolism,
Mortality,
)
24 changes: 12 additions & 12 deletions test/user/04-default_model.jl
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@
# Pick another functional response.
m = default_model(fw, ClassicResponse())
@test has_component(m, ClassicResponse)
@test maximum(m.attack_rate) == 50
@test maximum(m.attack_rate) ≈ 334.1719587843073

# Pick yet another functional response.
m = default_model(fw, LinearResponse())
Expand Down Expand Up @@ -71,7 +71,7 @@
@test has_component(m, NutrientIntake)
@test m.n_nutrients == 3
@test m.nutrients_turnover == [1 / 4, 1 / 4, 1 / 4]
@test m.nutrients_concentration == ones(2, 3)
@test m.nutrients_concentration == 0.5 .* ones(2, 3)
@sysfails(
m.K,
Property(K),
Expand All @@ -84,25 +84,25 @@
@test m.K == [0, 0, 1, 1]

# Any nutrient component triggers this default.
m = default_model(fw, Nutrients.Turnover([1, 2, 3, 4]))
m = default_model(fw, Nutrients.Turnover(0.8))
@test !has_component(m, LogisticGrowth)
@test has_component(m, NutrientIntake)
@test m.n_nutrients == 4
@test m.nutrients_turnover == [1, 2, 3, 4]
@test m.nutrients_concentration == ones(2, 4)
@test m.n_nutrients == m.n_producers == 2
@test m.nutrients_turnover == [0.8, 0.8]
@test m.nutrients_concentration == 0.5 .* ones(2, 2)
@sysfails(
m.K,
Property(K),
"A component '$CarryingCapacity' is required to read this property."
)

# Tweak directly from inside the aggregated blueprint.
m = default_model(fw, NutrientIntake(; turnover = [1, 2]))
m = default_model(fw, NutrientIntake(; turnover = 0.8))
@test !has_component(m, LogisticGrowth)
@test has_component(m, NutrientIntake)
@test m.n_nutrients == 2
@test m.nutrients_turnover == [1, 2]
@test m.nutrients_concentration == ones(2, 2)
@test m.n_nutrients == m.n_producers == 2
@test m.nutrients_turnover == [0.8, 0.8]
@test m.nutrients_concentration == 0.5 .* ones(2, 2)
@sysfails(
m.K,
Property(K),
Expand All @@ -111,8 +111,8 @@

# Combine if meaningful.
m = default_model(fw, Temperature(), NutrientIntake(; turnover = [1, 2]))
@test m.nutrients_supply == [10, 10]
@test m.attack_rate[1, 2] == 2.0452306245234897e-6
@test m.nutrients_supply == [4, 4]
@test m.attack_rate[1, 2] == 7.686741690921419e-7

# Add multiplex layers.
m = default_model(fw, CompetitionLayer(; A = (C = 0.2, sym = true), I = 2))
Expand Down
Loading
Loading