Problems with numberfleet

[lentinj]

Whilst trying to construct a model using tagging data, we tripped over numberfleet causing abundance to become NaN.

@willbutler42 can you add anything else about it? Hopefully turning igfs to a numberfleet will be enough to recreate our woes.

[willbutler42]

Just to add to this, I was trying numberfleet yesterday & in pretty much all circumstances (different models and setups) NaN's were produced

[lentinj]

Test case from Will:

library(tidyverse)
library(gadget3)
library(mfdb)

year_range <- 1952:1953

areas <- c('1' = 1)
timestep <- mfdb_timestep_quarterly

# Timekeeping for the model, i.e. how long we run for
time_actions <- list(
  g3a_time(start_year = min(year_range),
           end_year = max(year_range),
           timestep),
  list())

stock <-
  g3_stock('stock', seq(10, 20, 5)) %>%
  g3s_livesonareas(areas[c('1')]) %>%
  g3s_age(minage = 0, maxage = 1)

stock_actions <-
  list(
    g3a_initialconditions_normalparam(stock),
    g3a_age(stock),
    NULL
  )

survey <-
  g3_fleet('survey') %>%
  g3s_livesonareas(areas[c('1')])

survey_landings <-
  structure(expand.grid(year = year_range, step = 1:4) %>%
              mutate(area = 1, total_weight = 1),
            area_group = mfdb_group(`1` = 1))

numberfleet_actions <-
  list(
    g3a_predate_fleet(survey,
                      list(stock),
                      suitabilities =  g3_suitability_exponentiall50(g3_parameterized('alpha'),
                                                                     g3_parameterized('l50')),
                      catchability_f = g3a_predate_catchability_numberfleet(10)
                      )
  )

totalfleet_actions <-
  list(
    g3a_predate_fleet(survey,
                      list(stock),
                      suitabilities =  g3_suitability_exponentiall50(g3_parameterized('alpha'),
                                                                     g3_parameterized('l50')),
                      catchability_f = g3a_predate_catchability_totalfleet(10)
    )
  )

nactions <- c(time_actions, stock_actions, numberfleet_actions)
tactions <- c(time_actions, stock_actions, totalfleet_actions)

nactions <- c(nactions, list(g3a_report_history(nactions, '^survey__catch$')))
tactions <- c(tactions, list(g3a_report_history(tactions, '^survey__catch$')))

nummod <- g3_to_r(c(nactions))
totalmod <- g3_to_r(c(tactions))

par <- attr(nummod, 'parameter_template')

par[] <- 1
par$retro_years <- par$project_years <- 0

nres <- nummod(par)
tres <- totalmod(par)

attr(nres, 'stock__num') %>% as.data.frame.table() %>% pull(Freq) %>% sum()
attr(tres, 'stock__num') %>% as.data.frame.table() %>% pull(Freq) %>% sum()

attr(nres, 'hist_survey__catch')
attr(tres, 'hist_survey__catch')

You can see the differences in the generated modes with:

unattr <- function (x) { attributes(x) <- NULL ; x }
unittest::ut_cmp_identical(unattr(totalmod), unattr(nummod))

... which gives some very likely candidates of where things are going wrong.

[lentinj]

Predictably, the problem is the division by weight:

https://github.com/gadget-framework/gadget3/blob/a8fda621b26d2fe593ecc106fe25a8f6cda177c6/R/action_predate.R#L7-L11

Both stock__predby_fleet_stock and stock__wgt are zero, resulting in NaN. Quick answer is avoid_zero_vec, probably the right one but I wonder if there's scope for a general stock__wgt <- avoid_zero_vec(stock__wgt), and if that'd be quicker.

[lentinj]

The vague comment above is probably best interpreted as doing avoid_zero_vec(stock__wgt) when setting weight, i.e. in growth, transitions, initialconditions, renewal, .. Whilst it would simplify some code (a lot of the uses of avoid_zero_vec are to duck exactly this), I'm not convinced it'd be a performance gain, given the number of times we'd then have to do avoid_zero_vec. Worse, forgetting to do so will result in much harder to track bugs.

[willbutler42]

The zeros for stock__wgt seem to be entering during ageing at stock__minage, so would re-populating the stock__wgt with non-zeros here solve the problem? https://github.com/gadget-framework/gadget3/blob/master/R/action_age.R#L97

[lentinj]

Ouch, that's a particularly poor choice for default value. Thanks for spotting that! I wouldn't be surprised if that makes a few avoid_zero_vec instances redundant.

There's no real reason to wipe __wgt here anyway. Sure, it's not a useful value but wiping it to 1 isn't any better.

[lentinj]

If you're feeling brave, the avoid_zero_vec in g3l_distribution might well be redundant too. Would be interesting to give this a try with some larger models:

diff --git a/R/likelihood_distribution.R b/R/likelihood_distribution.R
index 39c066e..52528cc 100644
--- a/R/likelihood_distribution.R
+++ b/R/likelihood_distribution.R
@@ -303,7 +303,7 @@ g3l_distribution <- function (
                 if (compare_num) {
                     debug_trace("Take prey_stock__predby_fleet_stock weight, convert to individuals, add to our count")
                     stock_ss(modelstock__num) <- stock_ss(modelstock__num) +
-                        stock_reshape(modelstock, stock_ss(prey_stock__predby_fleet_stock) / avoid_zero_vec(stock_ss(prey_stock__wgt)))
+                        stock_reshape(modelstock, stock_ss(prey_stock__predby_fleet_stock) / stock_ss(prey_stock__wgt))
                 }
                 if (compare_wgt) {
                     debug_trace("Take prey_stock__predby_fleet_stock weight, add to our count")

[willbutler42]

Have not managed to test g3l_distribution yet, but I have been having some problems with larger numberfleet models despite the above commit (bcaa9b1), so feeling that avoid_zero_vec(ss_stock(stock__wgt)) is the way forward here

[lentinj]

Fair do's, done. This issue is drifting off-topic anyway. Added #89 if we want to pick this up at a later date.