diff --git a/Project.toml b/Project.toml index 5b4fd875..8c86756b 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,9 @@ uuid = "41889b7c-c35c-4f16-abd7-94b299d9fd29" version = "0.0.1" [deps] +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" +NamedArrays = "86f7a689-2022-50b4-a561-43c23ac3c673" OceanBioME = "a49af516-9db8-4be4-be45-1dad61c5a376" Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" diff --git a/examples/emergent_4P/emergent_functions.jl b/examples/emergent_4P/emergent_functions.jl new file mode 100644 index 00000000..cbd85825 --- /dev/null +++ b/examples/emergent_4P/emergent_functions.jl @@ -0,0 +1,77 @@ + +function dummy_emergent_growth(growth_a::Real, growth_b::Real, volume::Real) + rate = 0.0 + if growth_a == 0 + return rate + else + if volume == 1 + rate = 7.190e-6 + elseif volume == 10 + rate = 2.216e-5 + end + end + + return rate +end + +function dummy_emergent_predation_rate( + predation_rate_a::Real, predation_rate_b::Real, volume::Real +) + rate = 0 + if predation_rate_a == 0 + return rate = 0 # Early return if volume_a is zero + else + # Set rate based on the value of volume + if volume == 10 + rate = 8.86e-5 + elseif volume == 100 + rate = 4.88e-5 + end + end + + return rate +end + +function dummy_emergent_nitrogen_half_saturation( + nitrogen_half_saturation_a::Real, nitrogen_half_saturation_b::Real, volume::Real +) + if nitrogen_half_saturation_a == 0 + return rate = 0 # Early return if volume_a is zero + else + # Set rate based on the value of volume + if volume == 1 + rate = 6.73e-3 + elseif volume == 10 + rate = 0.12 + end + end + + return rate +end + +function dummy_emergent_palat(prey_data, predator_data) + prey_volume = prey_data["volumes"] + predator_volume = predator_data["volumes"] + optimum_predator_prey_ratio = predator_data["optimum_predator_prey_ratio"] + protection = prey_data["protection"] + + ratio = predator_volume / prey_volume + + if optimum_predator_prey_ratio == 0 + palat = 0.0 + elseif ratio == optimum_predator_prey_ratio + palat = 1 * (1 - protection) + else + palat = 0.3 * (1 - protection) + end + return palat +end + +function dummy_emergent_assimilation_efficiency(prey_data, predator_data) + assimilation_efficiency = 0 + # predators don't eat other predators + if prey_data["assimilation_efficiency"] == 0 + assimilation_efficiency = predator_data["assimilation_efficiency"] + end + return assimilation_efficiency +end diff --git a/src/Agate.jl b/src/Agate.jl index c29bac31..f22ea006 100644 --- a/src/Agate.jl +++ b/src/Agate.jl @@ -6,6 +6,7 @@ include("Models/Models.jl") using .Library using .Models +export compute_allometric_parameters export define_tracer_functions end # module diff --git a/src/Models/Models.jl b/src/Models/Models.jl index 80d79ff7..db3342ac 100644 --- a/src/Models/Models.jl +++ b/src/Models/Models.jl @@ -1,8 +1,12 @@ module Models include("Biogeochemistry.jl") +include("Parameters.jl") + using .Biogeochemistry +using .Parameters +export compute_allometric_parameters export define_tracer_functions end # module diff --git a/src/Models/Parameters.jl b/src/Models/Parameters.jl new file mode 100644 index 00000000..a9c1752e --- /dev/null +++ b/src/Models/Parameters.jl @@ -0,0 +1,219 @@ +module Parameters + +using DataStructures: DefaultDict +using NamedArrays + +export compute_allometric_parameters + +# TODO: the real palatability and assimilation functions should eventually be defined here +include(joinpath("..", "..", "examples", "emergent_4P", "emergent_functions.jl")) +emergent_palatability_f = dummy_emergent_palat +emergent_assimilation_efficiency_f = dummy_emergent_assimilation_efficiency + +# TODO: update this placeholder function (should only take in `a`, `b` and `volume`) +function allometry_f(param, a, b, volume) + if param == "maximum_growth_rate" + return dummy_emergent_growth(a, b, volume) + elseif param == "nitrogen_half_saturation" + return dummy_emergent_nitrogen_half_saturation(a, b, volume) + elseif param == "maximum_predation_rate" + return dummy_emergent_predation_rate(a, b, volume) + end +end + +""" + compute_allometric_parameters(plankton::Dict) -> Dict + +This function: + - generates `n` names for each `plankton` group (e.g., "P", "Z" or "cocco") of the form: + `["P1", ..., "P", "Z1", ...., "Z", ...]` + - generates `n` volumes for each `plankton` group using either a linear or a log + splitting scale + - optionally computes emergent parameters (allometric functions, assimilation matrix, + palatability matrix) for each group-volume combination + - reshapes any other group specific parameters (e.g., `linear_mortality`) to length `n` +All parameters are returned as: + `Dict( => , ....)` +using names generated in the first step. + +# Arguments +- `plankton`: a Dictionary of plankton groups' specific parameters of the form: + `Dict( => Dict( => , ....), ...)` + + The Dictionary for each group (e.g., "P", "Z" or "cocco") has to contain at least the + keys "n" and "volumes" and have the following form: + ``` + Dict( + "P" => Dict( + "n" => 2, + "volumes" => + Dict("min_volume" => 1, "max_volume" => 10, "splitting" => "log_splitting"), + ... + ), + "Z" => ... + ) + ``` + + The Dictionary for each group can optionally include the following keys and values: + - key: "allometry" (compute volume dependent parameters) + - value: Dictionary of the form + `Dict( => Dict("a" => , "b" => ), ...)` + - key: "palatability" (generate a palatability matrix) + - value: Dictionary with "optimum_predator_prey_ratio" and "protection" keys + - key: "assimilation_efficiency" (generate assimilation efficiency matrix) + - value: Dictionary with "assimilation_efficiency" key + For example: + ``` + Dict( + "P" => Dict( + ..., + "allometry" => Dict( + "maximum_growth_rate" => Dict("a" => 1, "b" => 1), + "nitrogen_half_saturation" => Dict("a" => 1, "b" => 1), + ), + ... + ), + "Z" => Dict( + ..., + "palatability" => Dict("optimum_predator_prey_ratio" => 10, "protection" => 1), + "assimilation_efficiency" => Dict("assimilation_efficiency" => 0.32), + ... + ) + ) + ``` +""" +function compute_allometric_parameters(plankton::Dict) + + # intermediate representations for parameters that are output as matrices + intermediate_palatability = DefaultDict{AbstractString,NamedArray}(NamedArray([], [])) + intermediate_assimilation = DefaultDict{AbstractString,NamedArray}(NamedArray([], [])) + # final outputs here + results = DefaultDict{AbstractString,NamedArray}(NamedArray([], [])) + + for (plankton_name, params) in plankton + n = params["n"] + + plankton_names = ["$plankton_name$i" for i in 1:n] + + # 1. compute volumes + splitting_function = getfield(Parameters, Symbol(params["volumes"]["splitting"])) + volumes = splitting_function( + params["volumes"]["min_volume"], params["volumes"]["max_volume"], n + ) + results["volumes"] = vcat(results["volumes"], NamedArray(volumes, plankton_names)) + + # 2. compute allometric functions (if any specified by user) + if "allometry" ∈ keys(params) + for (param, args) in params["allometry"] + # TODO: once `allometry_f` is updated, it should not have `param` as argument + values = [allometry_f(param, args["a"], args["b"], v) for v in volumes] + results[param] = vcat(results[param], NamedArray(values, plankton_names)) + end + end + + # 3. reshape remaining parameters + for (param, value) in params + if !(param ∈ ["n", "volumes", "allometry"]) + # NOTE: for palatability and assimilation_efficiency, `value` is a Dictionary + # store the values in intermediate form to use later when creating matrices + if param ∈ ["palatability", "assimilation_efficiency"] + for (arg, val) in value + values = NamedArray(repeat([val], n), plankton_names) + if param == "palatability" + intermediate_palatability[arg] = vcat( + intermediate_palatability[arg], values + ) + elseif param == "assimilation_efficiency" + intermediate_assimilation[arg] = vcat( + intermediate_assimilation[arg], values + ) + end + end + else + # NOTE: expect here that in all other cases `value` is a single number + results[param] = vcat( + results[param], NamedArray(repeat([value], n), plankton_names) + ) + end + end + end + end + + # 4. palatability & assimilation efficiency matrices (if specified by user) + if !(isempty(intermediate_palatability)) + intermediate_palatability["volumes"] = results["volumes"] + results["palatability_matrix"] = emergent_2D_array( + intermediate_palatability, emergent_palatability_f + ) + end + + if !(isempty(intermediate_assimilation)) + intermediate_assimilation["volumes"] = results["volumes"] + results["assimilation_efficiency_matrix"] = emergent_2D_array( + intermediate_assimilation, emergent_assimilation_efficiency_f + ) + end + + return results +end + +""" +Log splitting function to generate a set of volumes. +""" +function log_splitting(min_volume::Real, max_volume::Real, n::Int) + log_min = log10(min_volume) + log_max = log10(max_volume) + log_step = (log_max - log_min) / (n - 1) + return [10^(log_min + i * log_step) for i in 0:(n - 1)] +end + +""" +Linear splitting function to generate a set of volumes. +""" +function linear_splitting(min_volume::Real, max_volume::Real, n::Int) + linear_step = (max_volume - min_volume) / (n - 1) + return [min_volume + i * linear_step for i in 0:(n - 1)] +end + +""" +Apply function `func` to every every pair of `plankton` groups, returning a 2D NamedArray. + +The `plankton` Dictionary keys have to contain argument names of the function to calculate +and be of the form `Dict( => NamedArray(, ), ...)` + +# Example +```julia +function add_ab(prey, pred) + return prey["a"] + pred["b] +end + +plankton = Dict( + "a" => NamedArray([1, 2, 3], ["A1", "B2", "C3"]), + "b" => NamedArray([9, 8, 7], ["A1", "B2", "C3"]), +) + +values_matrix = emergent_2D_array(plankton, add_ab) +``` +""" +function emergent_2D_array(plankton, func) + plankton_names = names(plankton["volumes"])[1] + arg_names = keys(plankton) + values = zeros(Float64, length(plankton_names), length(plankton_names)) + plankton_matrix = NamedArray(values, (predator=plankton_names, prey=plankton_names)) + + # Q: is there a better way to populate this ?! + # Populate the NamedArray with calculated values + for pred_name in plankton_names + predator_data = Dict(arg => plankton[arg][pred_name] for arg in arg_names) + for prey_name in plankton_names + prey_data = Dict(arg => plankton[arg][prey_name] for arg in arg_names) + # Pass prey and predator data dictionaries to the function + plankton_matrix[pred_name, prey_name] = func(prey_data, predator_data) + end + end + + # Create and return a NamedArray with the values and names + return plankton_matrix +end + +end #module diff --git a/test/runtests.jl b/test/runtests.jl index 158a77b7..c8647ef4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,10 @@ using Agate using Test +# modules include("test_biogeochemistry.jl") +include("test_parameters.jl") + +# examples include("test_box_model.jl") include("test_n2p2zd.jl") diff --git a/test/test_parameters.jl b/test/test_parameters.jl new file mode 100644 index 00000000..c5d422b0 --- /dev/null +++ b/test/test_parameters.jl @@ -0,0 +1,86 @@ +using Agate + +using Oceananigans.Units + +# TODO: this will eventually be replaced with the actual emergent functions +include(joinpath("..", "examples", "emergent_4P", "emergent_functions.jl")) + +@testset "Models.Parameters" begin + @testset "create_allometric_parameters" begin + + # include the 0 parameters to make it comparable to the N2P2ZD example + defined_parameters = Dict( + "P" => Dict( + "n" => 2, + "volumes" => Dict( + "min_volume" => 1, + "max_volume" => 10, + "splitting" => "log_splitting", + ), + "allometry" => Dict( + "maximum_growth_rate" => Dict("a" => 1, "b" => 1), + "nitrogen_half_saturation" => Dict("a" => 1, "b" => 1), + "maximum_predation_rate" => Dict("a" => 0, "b" => 0), + ), + "palatability" => + Dict("optimum_predator_prey_ratio" => 0, "protection" => 0), + "assimilation_efficiency" => Dict( + "can_be_eaten" => 1, + "can_eat" => 0, + "assimilation_efficiency" => 0, + ), + "linear_mortality" => 8e-7 / second, + "holling_half_saturation" => 0, + "quadratic_mortality" => 0, + "alpha" => 0.1953 / day, + ), + "Z" => Dict( + "n" => 2, + "volumes" => Dict( + "min_volume" => 10, + "max_volume" => 100, + "splitting" => "linear_splitting", + ), + "allometry" => Dict( + "maximum_growth_rate" => Dict("a" => 0, "b" => 0), + "nitrogen_half_saturation" => Dict("a" => 0, "b" => 0), + "maximum_predation_rate" => Dict("a" => 1, "b" => 1), + ), + "palatability" => + Dict("optimum_predator_prey_ratio" => 10, "protection" => 1), + "assimilation_efficiency" => Dict( + "can_be_eaten" => 0, + "can_eat" => 1, + "assimilation_efficiency" => 0.32, + ), + "linear_mortality" => 8e-7 / second, + "holling_half_saturation" => 5.0, + "quadratic_mortality" => 1e-6 / second, + "alpha" => 1e-99 / day, + ), + ) + + emergent_parameters = compute_allometric_parameters(defined_parameters) + + # compare against hand computer `parameters` in examples + include(joinpath("..", "examples", "N2P2ZD", "tracers.jl")) + + plankton_order = ["P1", "P2", "Z1", "Z2"] + + for (key, emerge_params) in emergent_parameters + # start with arrays of values + if !(key ∈ ["assimilation_efficiency_matrix", "palatability_matrix", "volumes"]) + @test all( + parameters[Symbol(key)] .== [emerge_params[p] for p in plankton_order] + ) + # matrices of values -> compare row at a time + elseif !(key == "volumes") + for (i, p) in enumerate(plankton_order) + emergent_row = emerge_params[p, :] + true_row = parameters[Symbol(replace(key, "_matrix" => ""))][i, :] + @test all(true_row .== [emergent_row[p] for p in plankton_order]) + end + end + end + end +end