Remove discrepancies between FFT method and old method for calculating integrals about mizer HOT 16 CLOSED

Hi Gustav,

thanks for looking into this issue. Have you had any insights about the causes of discrepancies? Should we assume the discrepancies are just normal discretisation errors? We would like to proceed with developing parameter optimisation in mizer, but want to make sure the two integration methods are comparable. The FFT method is much faster, so it naturally would be preferable for the optimisation.

from mizer.

Unfortunately there is more than just rounding errors going on. This can be seen clearly when one looks at the predation rate. That rate should be zero at sizes greater than the predator, but the fft code allows feeding also in one size class above the predator's size. It seems that the truncation of the feeding kernel is not implemented correctly. I will try to fix that.

Another discrepancy arises when a species has a feeding kernel that extends considerably below the smallest phytoplankton. Allowing such a scenario is a bit weird, because it seems to artificially restrict the growth rate of small larvae of that species compared to slightly larger individuals. Unfortunately the default value for min_w_pp is such that it happens for herring and sandeel in the North Sea model, where both species have been given enormous values for beta. I guess this shows the difficulty of modelling filter feeders as having a constant predator/prey mass ratio.

The fft code was written under the assumption that that does not happen. If the feeding kernel extends to below the smallest phytoplankton then this is wrapped around to the largest fish. Because the abundance of large fish is very small compared to the abundance of the smallest phytoplankton, this error only has a very small effect on the energy available to the fish, but the error can accumulate during longer simulations. Again this is something I can fix.

from mizer.

I guess the reason for very large beta values for some species is that they should remain feeding on small prey even if they grow large. I have the same challenge with the Tasmanian model, and was wondering how to make sure that larvae still have food but adults stay feeding on small prey (e.g. planktivores). This will be partly remedied once we implement variable PPMR option (#41 ), although that will compromise the speed quite substantially, because FFT method won't be available.

The way I currently tried to solve this challenge in the Tasmanian model, is to increase the minimum body size of my functional groups. Currently the w_min for all groups is 10x the egg size, and perhaps it should be increased even further. This means that I don't track dynamics and mortality of eggs and tiny larvae, but I am not sure many multi-species models in mizer are equipped to do it either.

It would be interesting to hear what others think about it

from mizer.

from mizer.

@juliablanchard , I like your suggestion to set the minimum size of the plankton spectrum to coincide with the smallest size at which consumers feed. There is of course the slight problem that a Gaussian kernel has no lower cutoff. The fft code already sets a lower cutoff at 3 sigma. If there is no objection, I would like to impose the same cutoff also for the default kernel when not using fft. We could allow users to change this cutoff to another value if they add a min_ppmr slot to the species_params data frame.

While we are at it, I would also like to propose to set the minimum size of the consumer spectrum automatically to the smallest egg size, if provided in species_params, and to set the maximum size of the consumer spectrum to the largest asymptotic size. Any other choice just would slow down the simulations unnecessarily, and it makes it easier for the user if they do not have to think about how to choose these parameters.

from mizer.

from mizer.

I think I have identified one source of discrepancy between the old and new discretisation schemes. The issue has to do with the restriction that predators should not feed on prey larger than themselves. This is a restriction that both the old and new code tried to impose. However in order to know how to correctly implement that restriction, one has to be clear about how to approximate the integrals.

I believe that the idea in mizer was that it stores the abundance densities at grid points and then approximates the integrals for calculating the growth and death rates at these gridpoints by left Riemann sums. So for example, the total number of individuals between w[i] and w[i+1] is approximated by n[i] * dw[i] (where in mizer dw[i] = w[i+1]-w[i]). So it is the values at the left end of a size bracket that are used to estimate the integral over that size bracket.

Unfortunately, for the calculation of the predation rate the old code did not follow this convention. The integral over the contribution from all predators, if approximated by a left Riemann sum should get its first contribution from the size bracket just above the prey size. The contribution should be approximated in terms of the values at the left of this size bracket, i.e., at the prey size. The old code dropped that contribution. One way to say this is that the old code used a left Riemann sum for the growth rate but a right Riemann sum for the death rate. This is a quirk in the old code that I had not spotted when I wrote the fft code. I will now change the fft code to replicate that behaviour.

from mizer.

What I wrote at the end of the last post was wrong: the missing term in the old code can not be explained by saying that the code is using a right Riemann sum. It is unfortunately simply a left Riemann sum with a missing term. Unfortunately the old code is designed in a way that makes it difficult to fix the error. Of course as the size of the brackets gets smaller, the error from the missing term gets smaller. And we should always keep in mind that these discretisation errors are neglegible compared to all the other effects that the model does not capture. Mizer is not supposed to be a quantitatively precise model. Anyway, after I found my way through this I will write down carefully the details of how the discretisation is currently working in mizer.

from mizer.

I have now changed the fft code to agree with the old code. Also the code is now setting min_w_pp, min_w and max_w automatically. min_w_pp is set so that even the smallest fish has their entire feeding kernel inside the plankton spectrum. All that is left to do is to document the fft method.

from mizer.

Hi Gustav,

I am working with Asta and Julia on a version of Mizer that includes evolution in baldrech/mizerEvoMerging. Is it possible to merge your fix of the feeding kernel in my branch without having all the other commits you have in sizespectrum/mizer/dev ?

from mizer.

from mizer.

Hi @gustavdelius and thanks for the help with the merge. I have not yet had a chance to look, but was going to try and merge new updates next week. Perhaps it is easiest (given the time difference) if we go through the changes with Romain first, and contact you if things are not clear.

On a separate, but related issue - if the new fft code changes results, is there a way to keep earlier mizer versions available, so that people can reproduce results they published with the earlier version? Presumably with version control like git people can install whatever version they want or used for their old analyses. As long as they know which version to install...

from mizer.

Hi Asta, thanks for bringing up the issue of keeping track of which mizer version to use. I have put that into a separate issue #45

from mizer.

Thanks for your help Gustav.
Can we do the merge next week? I can do it with Asta's repository as I am collaborator and Asta could join if she has time (as with the time difference it's going to be pretty late in Tasmania).

from mizer.

Hi Romain, yes we can, but did you see Asta's suggestion in her post above? You can always reach me via Google hangouts on [email protected].

from mizer.

Yes, we met last Friday and decided to go with your plan of merging Asta's repository first while screen-sharing.

from mizer.