getDiet(proportion = FALSE) overcounts by (1+beta)/2 when second_order_w bin averaging is on

[gustavdelius]

Summary

With second_order_w(params) <- TRUE, summing getDiet(params, proportion = FALSE)
over prey no longer equals getEncounter(params) * (1 - getFeedingLevel(params)).
It is too large by exactly (1 + beta) / 2, where beta = w_full[2] / w_full[1]
is the (constant) grid ratio. For NS_params, beta = 1.1934, so getDiet()
overcounts by 9.67%.

On the default first-order path the identity holds to the usual FFT round-off.
And getDiet(proportion = TRUE) — the default — is unaffected, because the
factor is uniform and divides out in the normalisation. That is presumably why
this has not shown up before.

Reproduce

library(mizer)
p <- NS_params
second_order_w(p) <- TRUE

d <- apply(getDiet(p, proportion = FALSE), c(1, 2), sum)  # sum over prey
e <- getEncounter(p) * (1 - getFeedingLevel(p))
beta <- p@w_full[2] / p@w_full[1]

r <- (d / e)[, 30:80]
r <- r[r > 0]
range(r)          # 1.096681 1.096681
0.5 * (1 + beta)  # 1.096681

The same factor appears on the custom-kernel (non-FFT) branch:

pk <- getPredKernel(p)
comment(pk) <- "set manually"
q <- setPredKernel(p, pred_kernel = pk)

rk <- (apply(getDiet(q, proportion = FALSE), c(1, 2), sum) /
           (getEncounter(q) * (1 - getFeedingLevel(q))))[, 30:80]
range(rk[rk > 0])  # 1.096681 1.096681

Cause

The prey-bin quadrature is applied twice.

getDiet() builds its prey vector with the trapezoidally bin-averaged weight
(R/summary_methods.R:147-148, R/summary_methods.R:170-171):

w_eff      <- bin_average_weight(params@w)
w_full_eff <- bin_average_weight(params@w_full)
...
prey[1:no_sp, idx_sp] <- sweep(n, 2, w_eff * params@dw, "*")
prey[no_sp + 1, ]     <- n_pp * w_full_eff * params@dw_full

and then convolves that with params@ft_pred_kernel_e. But when bin_average
is on, setPredKernel() has already built ft_pred_kernel_e as the kernel
integrated over the prey bin (R/setPredKernel.R:243):

weight_e <- exp(2 * tt) * Delta / Q / (beta_grid - 1)

and the comment there says explicitly that dividing by (beta - 1) cancels the
w * dw factor "that mizerEncounter() and mizerPredRate() already fold into
the prey and predator vectors, so those rate functions need no change".
Accordingly mizerEncounter() keeps the plain point weight
(R/project_methods.R:343):

prey <- sweep(prey, 2, params@w_full * params@dw_full, "*")

So getDiet() applies the second-order prey correction on top of a kernel that
already contains it. On a geometric grid bin_average_weight(w) / w is exactly
(1 + beta) / 2, which is the observed factor.

The custom-kernel branch has the same mismatch for a different reason: there
params@pred_kernel is point-sampled and mizerEncounter() weights it with
plain w, while getDiet() weights it with w_eff.

Suggested fix

In both branches of getDiet(), use params@w / params@w_full rather than
their bin-averaged versions, matching mizerEncounter(). The prey-bin integral
belongs in the kernel, where setPredKernel() already puts it.

Possibly the same issue in getTrophicLevel()

getTrophicLevel() builds its TL-weighted prey mass with
bin_average_summary_weight(params@w, params) and multiplies it by
getPredKernel() (R/summary_methods.R:365-384), while its denominator
accumulator is (1 - feeding_level) * encounter taken from getEncounter()
(R/summary_methods.R:353-357). That is the same numerator/denominator
quadrature mismatch. The net effect is not a uniform factor there — it is a
ratio, and the resource contribution enters through a separately weighted term —
but the trophic levels do move when the flag is flipped:

species	first order	bin averaged
Sprat	4.200	4.153
Sandeel	3.865	3.822
N.pout	5.302	5.358
Herring	3.981	3.938
Dab	5.127	5.166
Cod	5.161	5.204

(mean over getTrophicLevel()[, 60:90].) I have not chased that one down to the
same level of certainty as the getDiet() factor.

Why it matters downstream

In mizerEcopath the diet matrix
is obtained by integrating getDiet(proportion = FALSE) over predator size, and
it is a modelling invariant that summing it over prey reproduces the Ecopath
consumption rate Q. That invariant holds on the default path and breaks by
~10% under second_order_w, which is what led me here.

Version

mizer 3.2.1.9001, commit 9e9b922.