get_gamma_default() ignores the external encounter rate that get_f0_default() counts

[gustavdelius]

Split out from the discussion on #577, which deliberately left this case alone.

get_gamma_default() chooses the search volume coefficient gamma that gives a species the target feeding level f0 in a reference state: no fish, and a resource spectrum on the power law $\kappa w^{-\lambda}$. It measures the available energy in that state by giving the species a search volume coefficient of 1, then divides.

mizerEncounter() adds params@ext_encounter (and any component's encounter_fun) to the predation encounter, so that measurement includes them. But it includes them in the unit-gamma normalisation, where they are negligible, and it divides by them where the equation calls for a subtraction. The net effect is that get_gamma_default() ignores the external encounter altogether, while get_f0_default() — which measures with the given gamma, in real units — counts it at full weight. The two stop being inverses of each other, and a model with a substantial external encounter rate does not reach the f0 it asked for.

The arithmetic

In the reference state, with the default q = lambda - 2 + n so that both contributions are the same power of $w$, write the predation encounter of species $i$ at unit gamma as $A_i w^{n_i}$ and the external encounter as $c_i w^{n_i}$. The feeding level is $f_0$ when

$$\gamma_i A_i + c_i = h_i \frac{f_0}{1 - f_0}$$

so the answer wanted is $\gamma_i = \left(h_i f_0/(1-f_0) - c_i\right) / A_i$.

What the code computes is

$$\gamma_i = \frac{h_i}{A_i + c_i}\cdot\frac{f_0}{1 - f_0}$$

The external encounter therefore enters only as a relative correction of order $c_i / A_i$. But $A_i$ is the available energy at unit gamma, so it is about $1/\gamma_i \approx 10^{11}$ times larger than the encounter the species actually experiences, which makes $c_i/A_i \sim \gamma_i \sim 10^{-11}$. The correction is far below the size of the effect it is meant to represent, and gamma comes out as though there were no external encounter at all.

Reproducible example

library(mizer)
packageVersion("mizer")
#> [1] '3.3.0.9000'

# The feeding level actually achieved in the reference state that the `gamma`
# default is derived from: no fish, resource on the power law, read at w_mat
ref_f0 <- function(params) {
    p <- params
    p@initial_n[] <- 0
    p@initial_n_pp[] <- p@resource_params$kappa *
        p@w_full ^ (-p@resource_params$lambda)
    f <- getFeedingLevel(p)
    idx <- vapply(species_params(p)$w_mat,
                  function(w) sum(p@w <= w), integer(1))
    vapply(seq_len(nrow(f)), function(i) f[i, idx[i]], numeric(1))
}

sp <- data.frame(species = c("sp1", "sp2"), w_max = c(150, 1500),
                 w_inf = c(100, 1000), k_vb = c(0.3, 0.2),
                 w_mat = c(10, 100), beta = 100, sigma = 2)
params <- newMultispeciesParams(sp, no_w = 50, info_level = 0)

gamma_0 <- species_params(params)$gamma
ref_f0(params)          # the default delivers the target f0 = 0.6 exactly
#> [1] 0.6 0.6

# Now add an external encounter rate of the same size as the predation
# encounter the species get from the resource
enc <- getEncounter(params)
idx <- vapply(species_params(params)$w_mat,
              function(w) sum(params@w <= w), integer(1))
n <- species_params(params)[["n"]]
E_ext <- vapply(1:2, function(i) enc[i, idx[i]] / params@w[idx[i]]^n[i],
                numeric(1))

with_ext <- params
given_species_params(with_ext)$E_ext <- E_ext
given_species_params(with_ext)$gamma <- NULL   # ask mizer to recalculate gamma

gamma_1 <- species_params(with_ext)$gamma
rbind(without = gamma_0, with = gamma_1, ratio = gamma_1 / gamma_0)
#>                  sp1          sp2
#> without 2.541851e-11 3.746056e-11
#> with    2.541851e-11 3.746056e-11
#> ratio   1.000000e+00 1.000000e+00

ref_f0(with_ext)        # the feeding level actually achieved
#> [1] 0.7500002 0.7417532

# The inverse function, on the very same model
given_species_params(with_ext)$gamma <- species_params(with_ext)$gamma
mizer::get_f0_default(with_ext)
#>       sp1       sp2
#> 0.7500002 0.7417532

gamma is unchanged to every printed digit — the relative difference is about $-2.5\times10^{-11}$, exactly the order the arithmetic above predicts — and the species end up at a feeding level of 0.75 instead of the 0.6 they asked for. get_f0_default() on the same model reports 0.75, the true feeding level with the external encounter included, so f0 → gamma → f0 does not round-trip.

How this differs from #577

The extension case fixed in #577 was worse in kind. An extension scaling search_vol multiplies a quantity that gamma itself determines, so its factor was re-applied on every rebuild, and a factor of zero for some species was a hard error. The external encounter rate is additive and independent of gamma, so it is a single silent offset with a well-defined fixed point. It is also part of mizer's own model rather than something arriving through the extension dispatch chain, which is why #577 left it in place.

Which behaviour is wanted?

There are two self-consistent readings, and mizer currently implements neither:

1. The defaults account for the external encounter. f0 is the feeding level a species actually reaches, so gamma should be solved from $\gamma A + c = h f_0/(1-f_0)$, subtracting the external contribution instead of dividing by it. get_f0_default() already works this way. This needs a guard: where $c_i$ already exceeds $h_i f_0/(1-f_0)$ no positive gamma can reach the target, and that should be an informative error rather than a negative search volume.
2. The defaults ignore it, on the grounds that gamma describes predation on the modelled spectrum and the external encounter is a separate addition on top. Then get_gamma_default() is right as it stands (by accident) and get_f0_default() should do its measuring with ext_encounter zeroed, so that the pair round-trips.

I lean towards 1, since that is what f0 means everywhere else in mizer and what get_f0_default() already reports.

A smaller point that applies under either reading: the measurement is read off at a single size, [, length(params@w)], the largest size on the grid. For the predation part that is exact, because it is a pure power law in $w$ and any size gives the same coefficient. An external encounter rate set by hand with setExtEncounter() need not be a power law at all, and the top of the grid is above the w_max of every species but the largest — a region that mizerEncounter()'s own help page describes as meaningless. Under reading 1 the correction would be taken at a size the species may never occupy, so it would have to be read at a size inside the species' own range, or integrated over it.

All of the above applies equally to the encounter contributed by a component registered with setComponent(..., encounter_fun = ).

[gustavdelius]

I think we want option 2. These defaults should really be calculated under the assumption of power law encounter.