Refactor emergent parameters code

Project.toml

@@ -3,17 +3,21 @@ 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"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [compat]
 DifferentialEquations = "7.15.0"
 OceanBioME = "0.13.2"
 Oceananigans = "0.93.3"
+StatsPlots = "0.15.7"
 Turing = "0.35.3"
 julia = "1.9"
 

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/NPZD/functions.jl

@@ -1,8 +1,6 @@
 monod_limitation(R, k) = R / (k + R)
 holling_type_2(R, k) = R / (k + R)
 
-Q₁₀(T) = 1.88^(T / 10) # T in °C
-
 smith_light_limitation(PAR, α, μ₀) = α * PAR / sqrt(μ₀^2 + α^2 * PAR^2)
 
 idealized_remineralization(D, r) = r * D
@@ -22,9 +20,3 @@ end
 linear_loss(P, l) = l * P
 
 quadratic_loss(P, l) = l * P^2
-
-function summed_linear_loss(P, l)
-    return sum([
-        linear_loss(P[plankton_index], l[plankton_index]) for plankton_index in eachindex(P)
-    ])
-end

examples/NPZD/tracers.jl

@@ -12,14 +12,13 @@ parameters = (
     kₚ=0.5573,
     β=0.9116,
     lᶻᵈ=0.3395 / day,
-    lⁿ=[0.066, 0.0102] / day,  #[lᵖⁿ, lᶻⁿ]
     rᵈⁿ=0.1213 / day,
     α=0.1953 / day,
 )
 
 tracers = Dict(
     "N" => :(
-        summed_linear_loss([P, Z], lⁿ) + idealized_remineralization(D, rᵈⁿ) -
+        linear_loss(P, lᵖⁿ) + linear_loss(P, lᶻⁿ) + idealized_remineralization(D, rᵈⁿ) -
         idealized_photosynthetic_growth(N, P, PAR, μ₀, kₙ, α)
     ),
     "D" => :(

examples/emergent_4P/emergent_functions.jl

@@ -0,0 +1,78 @@
+
+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
+    # Q: is this the right condition for determining whether to retrieve this value?
+    if prey_data["assimilation_efficiency"] == 0
+        assimilation_efficiency = predator_data["assimilation_efficiency"]
+    end
+    return assimilation_efficiency
+end

examples/emergent_4P/example.jl

@@ -0,0 +1,87 @@
+using Oceananigans.Units
+using Agate
+
+# Q: SHOULD THIS BE A NAMEDTUPLE INSTEAD?
+# TODO: users shouldn't have to specify 0 parameters --> do we just assume 0 and fill it OR can we update tracers to do without them entirely?
+# NOTE: assimilation efficiency has to create an array of values rather than compute anything --> is there a better way to do that?
+# NOTE: don't need to declare units as per second since this is the default time unit i.e., 1 == 1/sec
+
+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),
+        # anything else goes here
+        "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
+        ),
+        # anything else goes here
+        "linear_mortality" => 8e-7 / second,
+        "holling_half_saturation" => 5.0,
+        "quadratic_mortality" => 1e-6 / second,
+        # this should be 0
+        "alpha" => 1e-99 / day,
+    ),
+)
+
+emergent_parameters = compute_darwin_parameters(defined_parameters)
+println(emergent_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_matrix", "palatability_matrix", "volumes"])
+        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 == "volumes")
+        for (i, p) in enumerate(plankton_order)
+            emergent_row = params[p, :]
+            true_row = parameters[Symbol(replace(key, "_matrix" => ""))][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()...