Refactor emergent parameters code

Project.toml

@@ -3,9 +3,10 @@ uuid = "41889b7c-c35c-4f16-abd7-94b299d9fd29"
 version = "0.0.1"
 
 [deps]
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+NamedArrays = "86f7a689-2022-50b4-a561-43c23ac3c673"
 OceanBioME = "a49af516-9db8-4be4-be45-1dad61c5a376"
 Oceananigans = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

examples/N2P2ZD/tracers.jl

@@ -9,7 +9,7 @@ parameters = (
     linear_mortality=[8e-7, 8e-7, 8e-7, 8e-7] / second,
     quadratic_mortality=[0, 0, 1e-6, 1e-6] / second,
     maximum_predation_rate=[0, 0, 8.86e-5, 4.88e-5] / second,
-    alpha=[0.5, 0.1953, 1e-99, 1e-99] / day,
+    alpha=[0.1953, 0.1953, 1e-99, 1e-99] / day,
     feeding_export_poc_doc_fraction=0.5,
     mortality_export_fraction=0.5,
     palatability=[

examples/differential_equations.jl

@@ -1,7 +1,7 @@
 # This example shows how to integrate a Biogeochemistry (BGC) box model with DifferentialEquations.
 using Agate.Library.Light
 
-using DifferentialEquations
+using OrdinaryDiffEq
 using Plots
 
 using Oceananigans.Units

examples/emergent_4P/emergent_functions.jl

@@ -0,0 +1,91 @@
+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_palat(
+    prey_data,
+    predator_data;
+    volume_key::String="volume",
+    optimum_ratio_key::String="optimum_predator_prey_ratio",
+    protection_key::String="protection",
+)
+    prey_volume = prey_data[volume_key]
+    predator_volume = predator_data[volume_key]
+    optimum_predator_prey_ratio = predator_data[optimum_ratio_key]
+    protection = prey_data[protection_key]
+
+    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_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
+
+"""
+Default fall-back function if no emergent function is defined
+"""
+function default_emergent(my_parameter::Real)
+    return my_parameter
+end
+
+function dummy_emergent_assimilation_efficiency(
+    prey_data, predator_data; assimilation_efficiency_key::String="assimilation_efficiency"
+)
+    assimilation_efficiency = 0
+    # predators don't eat other predators
+    if prey_data[assimilation_efficiency_key] == 0
+        assimilation_efficiency = predator_data[assimilation_efficiency_key]
+    end
+    return assimilation_efficiency
+end

examples/emergent_4P/example.jl

@@ -0,0 +1,78 @@
+using Oceananigans.Units
+using Agate
+
+defined_parameters = Dict(
+    "P" => Dict(
+        "n" => 2,
+        "min_volume" => 1,
+        "max_volume" => 10,
+        "splitting" => "log_splitting",
+        "growth_a" => 1,
+        "growth_b" => 1,
+        "protection" => 0,
+        "optimum_predator_prey_ratio" => 0,
+        "nitrogen_half_saturation_a" => 1,
+        "nitrogen_half_saturation_b" => 1,
+        "predation_rate_a" => 0,
+        "predation_rate_b" => 0,
+        "linear_mortality" => 8e-7 / second,
+        "holling_half_saturation" => 0,
+        "quadratic_mortality" => 0,
+        "alpha" => 0.1953 / day,
+        "assimilation_efficiency" => 0,
+    ),
+    "Z" => Dict(
+        "n" => 2,
+        "min_volume" => 10,
+        "max_volume" => 100,
+        "splitting" => "linear_splitting",
+        "growth_a" => 0,
+        "growth_b" => 0,
+        "protection" => 1,
+        "optimum_predator_prey_ratio" => 10,
+        "nitrogen_half_saturation_a" => 0,
+        "nitrogen_half_saturation_b" => 0,
+        "predation_rate_a" => 1,
+        "predation_rate_b" => 1,
+        "linear_mortality" => 8e-7 / second,
+        "holling_half_saturation" => 5.0,
+        "quadratic_mortality" => 1e-6 / second,
+        "alpha" => 1e-99 / day,
+        "assimilation_efficiency" => 0.32,
+    ),
+)
+
+emergent_parameters = compute_darwin_parameters(defined_parameters)
+
+# sanity check --> compare generated values to parameters in N2P2ZD example
+# NOTE: eventually the emergent functions will be updated to true ones
+# that's why we don't include this in tests
+
+include(joinpath("..", "N2P2ZD", "tracers.jl"))
+
+plankton_order = ["P1", "P2", "Z1", "Z2"]
+
+for (key, params) in emergent_parameters
+    # start with arrays of values
+    if !(key ∈ ["assimilation_efficiency", "palatability", "volume"])
+        comparison = all(parameters[Symbol(key)] .== [params[p] for p in plankton_order])
+        println(key, " values are the same: ", comparison)
+        # matrices of values -> compare row at a time
+    elseif !(key == "volume")
+        for (i, p) in enumerate(plankton_order)
+            emergent_row = params[p, :]
+            true_row = parameters[Symbol(key)][i, :]
+            comparison = all(true_row .== [emergent_row[p] for p in plankton_order])
+            println(key, " ", p, " values are the same: ", comparison)
+        end
+    end
+end
+
+# for simplicity define the biogeochemistry dict seperately
+biogeochemistry_parameters = Dict(
+    "detritus_remineralization" => 0.1213 / day,
+    "feeding_export_poc_doc_fraction" => 0.5,
+    "mortality_export_fraction" => 0.5,
+)
+created_parameters = merge(biogeochemistry_parameters, emergent_parameters)
+#note that this dictionary would need to be converted to a named tuple to work with create_bgc_struc()...

examples/inference_hmc.jl

@@ -5,7 +5,7 @@
 using Agate
 using Agate.Library.Light
 
-using DifferentialEquations
+using OrdinaryDiffEq
 using LinearAlgebra
 using Plots
 using StatsPlots

src/Agate.jl

@@ -5,6 +5,7 @@ include("Models/Models.jl")
 
 using .Models
 
+export compute_darwin_parameters
 export define_tracer_functions
 
 end # module

src/Models/Dynamic.jl

@@ -18,7 +18,7 @@ export define_tracer_functions
         parameters, tracers; auxiliary_fields=[:PAR], helper_functions=nothing
     ) -> DataType
 
-Creates an Oceananigans.Biogeochemistry model type.
+Create an Oceananigans.Biogeochemistry model type.
 
 # Arguments
 - `parameters`: named sequence of values of the form ((<field name> = <default value>, ...)
@@ -61,8 +61,8 @@ end
 Create a subtype of Oceananigans.Biogeochemistry with field names defined in `parameters`.
 
 # Arguments
-- `struct_name`: name for the new struct passed as a Symbol. The struct will be accessible
-   as: `Agate.Models.Dynamic.<struct_name>`
+- `struct_name`: name for the struct to create passed as a Symbol (the new struct will be
+   accessible as `Agate.Models.Dynamic.<struct_name>`)
 - `parameters`: named sequence of values of the form (<field name> = <default value>, ...)
 
 Note that the field names defined in `parameters` can't be any of [:x, :y, :z, :t] as these
@@ -183,7 +183,7 @@ end
 """
     parse_expression(f_expr) -> Vector
 
-Returns all symbols (argument names and method names) called in expression.
+Return all symbols (argument names and method names) called in expression.
 
 # Example
 ```julia
@@ -211,7 +211,7 @@ end
 """
     expression_check(args, f_expr) -> nothing
 
-Checks that all methods and arguments are defined. Specifically:
+Check that all methods and arguments are defined. Specifically:
     - vector `args` contains all arguments of expression `f_expr`
     - all methods called in `f_expr` are defined in module (e.g., Base, Main, Agate)
 If not, throws an UnderVarError.

src/Models/Models.jl

@@ -1,8 +1,12 @@
 module Models
 
 include("Dynamic.jl")
+include("Parameters.jl")
+
 using .Dynamic
+using .Parameters
 
+export compute_darwin_parameters
 export define_tracer_functions
 
 end # module