Fix double-counted prey-bin quadrature in getDiet() and getTrophicLevel() - #476
Open
gustavdelius wants to merge 1 commit into
Open
Fix double-counted prey-bin quadrature in getDiet() and getTrophicLevel()#476gustavdelius wants to merge 1 commit into
gustavdelius wants to merge 1 commit into
Conversation
…el() Under `second_order_w(params) <- TRUE`, `setPredKernel()` builds `ft_pred_kernel_e` as the predation kernel integrated over the prey bin, dividing by `beta - 1` so that the plain `w * dw` supplied by the prey vector cancels. `mizerEncounter()` accordingly keeps the point weight. `getDiet()` did not: it weighted its prey vectors with the bin-averaged `w`, applying the prey-bin quadrature a second time. The result was uniformly too large by `(1 + beta) / 2` — 9.7% for `NS_params`. Both branches (FFT and custom kernel) were affected. The default `proportion = TRUE` hid it, because a uniform factor divides out in the normalisation. `getTrophicLevel()` had the same mismatch plus a second one: its trophic-level-weighted numerator was built from the point-sampled `getPredKernel()` while its denominator came from `getEncounter()`, which uses the bin-integrated kernel once the flag is on. The two are then no longer the same integral and the trophic levels moved by up to 0.06. Both now use the plain point weight `w * dw`. New internal `encounter_kernel()` returns the kernel `mizerEncounter()` actually uses — recovered from `ft_pred_kernel_e` by inverse FFT, so it cannot drift from whatever quadrature `setPredKernel()` chose — and `getTrophicLevel()` builds its numerator from that. `getPredKernel()` is refactored onto the shared `expand_kernel_offsets()` helper; its behaviour is unchanged. Tests cover the `encounter_kernel()`/`getEncounter()` identity, the diet-sums-to-consumption identity, and a predator whose prey all have trophic level 1 coming out at exactly 2 — each across both settings of the flag and both kernel branches. All fail against the previous code. Fixes #474 Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #474.
The bug
Under
second_order_w(params) <- TRUE,setPredKernel()buildsft_pred_kernel_eas the predation kernel integrated over the prey bin, dividing bybeta - 1so that the plainw * dwsupplied by the prey vector cancels.mizerEncounter()accordingly keeps the point weightparams@w_full * params@dw_full.getDiet()did not. It weighted its prey vectors with the trapezoidally bin-averagedw, applying the prey-bin quadrature a second time on top of a kernel that already contained it. On a geometric gridbin_average_weight(w) / wis exactly(1 + beta) / 2, so the diet came out uniformly too large by that factor — 9.7% forNS_params. Both branches were affected: the FFT path for the reason above, and the custom-kernel path because thereparams@pred_kernelis point-sampled whilegetDiet()weighted it with the bin average.The default
proportion = TRUEhid this, since a uniform factor divides out in the normalisation.getTrophicLevel()had the same mismatch plus a second one. Its trophic-level-weighted numerator was built from the point-sampledgetPredKernel(), while its denominator came fromgetEncounter(), which uses the bin-integrated kernel once the flag is on. The two are then no longer the same integral, and the trophic levels moved by up to 0.06 when the flag was flipped.The fix
Both functions now use the plain point weight
w * dw, matchingmizerEncounter(). The prey-bin integral belongs in the kernel, wheresetPredKernel()already puts it.New internal
encounter_kernel()returns the kernel thatmizerEncounter()actually uses: the point-sampled kernel on the default path, the bin-integrated one undersecond_order_w, and the stored array when a custom kernel is set. It recovers the bin-integrated weights fromft_pred_kernel_eby inverse FFT rather than duplicating the quadrature code, so it cannot drift from whateversetPredKernel()chose.getTrophicLevel()builds its numerator from that.getPredKernel()is refactored onto a sharedexpand_kernel_offsets()helper; its behaviour is unchanged. It remains the point-sampled kernel — the right object for plotting and for supplying a custom kernel — which is why the internal helper is separate.Effect
Summing
getDiet(proportion = FALSE)over prey now reproducesgetEncounter() * (1 - getFeedingLevel())in all four combinations of flag and kernel branch. Trophic levels under the flag now sit ~0.014 below the first-order values, consistently signed, instead of moving by up to ±0.06 with mixed signs.Nothing changes on the default path.
Tests
Four new tests in
test-second_order_summary.R, each run under both settings of the flag and both kernel branches:encounter_kernel()contracted with the prey vector reproducesgetEncounter()getDiet(proportion = FALSE)summed over prey equals the consumption rateEach fails against the previous code (the trophic-level numerator/denominator ratio was 0.97 rather than 1). Full suite passes: 3545 tests, 0 failures.
🤖 Generated with Claude Code