gadget-framework / gadget3

TMB-based gadget implemtation
GNU General Public License v2.0
8 stars 6 forks source link

list_to_stock_switch() means we can't use sub-formulas in suitability definitions #160

Open lentinj opened 2 weeks ago

lentinj commented 2 weeks ago
--- a/R/action_predate.R
+++ b/R/action_predate.R
@@ -210,7 +210,9 @@ g3a_predate <- function (

         # Main predation step, iterate over prey and pull out everything this fleet needs
         catchability <- f_substitute(catchability_f$suit, list(suit_f = quote(suitability)))
-        environment(catchability)$suitability <- list_to_stock_switch(suitabilities)
+        environment(catchability)$suitability <- suitabilities[[stock$name]]
+        print(catchability)
+        print(environment(environment(catchability)$suitability))
         out[[step_id(run_at, 1, predstock, stock, action_name)]] <- g3_step(f_substitute(~{
             debug_label("g3a_predate_fleet for ", stock)
             debug_trace("Zero ", predstock, "-", stock, " biomass-consuming counter")
@@ -301,7 +303,7 @@ g3a_predate <- function (
         do.call(c, lapply(prey_stocks, function (stock) g3a_suitability_report(
             predstock,
             stock,
-            list_to_stock_switch(suitabilities) )))))
+            suitabilities[[stock$name]] )))))
     return(as.list(out))
 }

diff --git a/R/suitability.R b/R/suitability.R
index 055fff7..d4b836b 100644
--- a/R/suitability.R
+++ b/R/suitability.R
@@ -13,18 +13,18 @@ g3_suitability_exponentiall50 <- function (
 }

 g3_suitability_andersen <- function (p0, p1, p2, p3 = p4, p4, p5 = quote(predator_length)) {
-  f_substitute(~p0 +
-                 avoid_zero(p2) * exp(-(log(avoid_zero_vec(p5/stock__midlen)) - p1)**2/avoid_zero(p3)) *
-                 bounded_vec(1000*(p1 - log(avoid_zero_vec(p5/stock__midlen))),0,1) +
-                 avoid_zero(p2) * exp(-(log(avoid_zero_vec(p5/stock__midlen)) - p1)**2/avoid_zero(p4)) *
-                 bounded_vec(1000*(log(avoid_zero_vec(p5/stock__midlen)) - p1),0,1),
-               list(
-                 p0 = p0,
-                 p1 = p1,
-                 p2 = p2,
-                 p3 = p3,
-                 p4 = p4,
-                 p5 = p5))
+  g3_formula(
+    quote(p0 +
+      p2 * exp(-(log(p5 - p1)**2/p3)) *
+      bounded_vec(1000*(p1 - log(p5)),0,1) +
+      p2 * exp(-(log(p5 - p1)**2/p4)) *
+      bounded_vec(1000*(log(p5) - p1),0,1) ),
+    p0 = f_substitute(quote( p0 ), list(p0 = p0)),
+    p1 = f_substitute(quote( p1 ), list(p1 = p1)),
+    p2 = f_substitute(quote( avoid_zero(p2) ), list(p2 = p2)),
+    p3 = f_substitute(quote( avoid_zero(p3) ), list(p3 = p3)),
+    p4 = f_substitute(quote( avoid_zero(p4) ), list(p4 = p4)),
+    p5 = f_substitute(quote( avoid_zero_vec(p5/stock__midlen) ), list(p5 = p5)) )
 }

 g3_suitability_andersenfleet <- function (
@@ -42,6 +42,19 @@ g3_suitability_andersenfleet <- function (
         by_stock = TRUE,
         by_predator = TRUE,
         exponentiate = TRUE) {
+  g3_formula(
+    quote(p0 +
+      p2 * exp(-(log(p5 - p1)**2/p3)) *
+      bounded_vec(1000*(p1 - log(p5)),0,1) +
+      p2 * exp(-(log(p5 - p1)**2/p4)) *
+      bounded_vec(1000*(log(p5) - p1),0,1) ),
+    p0 = f_substitute(quote( p0 ), list(p0 = p0)),
+    p1 = f_substitute(quote( p1 ), list(p1 = p1)),
+    p2 = f_substitute(quote( avoid_zero(p2) ), list(p2 = p2)),
+    p3 = f_substitute(quote( avoid_zero(p3) ), list(p3 = p3)),
+    p4 = f_substitute(quote( avoid_zero(p4) ), list(p4 = p4)),
+    p5 = f_substitute(quote( p5/stock__midlen ), list(p5 = p5)) )
+
     f_substitute(~p0 +
         p2 * exp(-(log(p5/stock__midlen) - p1)**2/p3) *
         bounded_vec(1000*(p1 - log(p5/stock__midlen)),0,1) +
@@ -100,7 +113,7 @@ g3a_suitability_report <- function (

   suit_f <- g3_step(f_substitute(~stock_with(stock, suit_f), list(suit_f = suit_f)), recursing = TRUE)  # Resolve stock_switch

-  suit_dims <- all.vars(suit_f)
+  suit_dims <- all_undefined_vars(suit_f, recursive = TRUE)
   if ("cur_step" %in% suit_dims) stop("Can't generate a time-varying suitability report")
   # Special case, swap use of stock__midlen with general iterator name
   suit_dims[suit_dims == paste0(stock$name, "__midlen")] <- "length"
lentinj commented 2 weeks ago

g3_eval() doesn't support these either, which is another problem to doing the above.

lentinj commented 2 weeks ago

Alternatively, we're already generating suit_...__report, why doesn't our suitability become suit_...__report[vars]?