gadget-framework / gadget3

TMB-based gadget implemtation
GNU General Public License v2.0
8 stars 6 forks source link

Spawning into multiple stocks not splitting evenly #147

Closed lentinj closed 5 months ago

lentinj commented 5 months ago

Initially, __spawnednum contains the stock structure we expect for each newly spawned individuals. We normalise this, then try and scale both by the offspring we expect the parents to make:

fish_imm1__spawnednum <- (fish_imm1__spawnednum/(sum(fish_imm2__spawnednum) + sum(fish_imm1__spawnednum)))
     * sum(fish_mat__offspringnum) * 0.5
fish_imm2__spawnednum <- (fish_imm2__spawnednum/(sum(fish_imm2__spawnednum) + sum(fish_imm1__spawnednum)))
     * sum(fish_mat__offspringnum) * 0.5

...but in the second case fish_imm1__spawnednum has already been scaled, so this sum is nonsense. It's not immediately obvious why we need to scale __spawnednum by the sum of both imm1 & imm2, but according to comments this is what happened in gadget2.

Also, we're not considering area in the above. Offspring should presumably end up in the area their parents were in.

vbartolino commented 5 months ago

This clarifies the current behaviour, thanks. 0.5 is a vector input by the user, so it may be enough to split simply offspringnumby that proportion, which is probably what you're suggesting. Something like?

fish_imm1__spawnednum <- sum(fish_mat__offspringnum) * 0.5
fish_imm2__spawnednum <- sum(fish_mat__offspringnum) * 0.5

I agree that it may be useful to have the area specified, but I'm trying to think what we want it to represent. Thinking general (and aloud), should this cover alternatives like:

bthe commented 5 months ago

Is the formulation in https://github.com/gadget-framework/gadget3/blob/master/R/action_spawn.R#L159 is correct. Should it be just:

output_stock__spawnednum <- stock__offspringnum * output_ratio

We might have multiple spawning events so probably something like:

 output_stock__spawnednum <- stock__offspringnum * output_ratio + output_stock__spawnednum
lentinj commented 5 months ago

0.5 is a vector input by the user

Yeah, this is the output_ratios input.

so it may be enough to split simply offspringnum by that proportion

Not quite. __offspringnum gives us the number of offspring, broken down by parent area/age/length. We need to apply this to a child length distribution using the mean_f, stddev_f inputs. This is what is in fish_imm1__spawnednum. The following seems more like what we want:

fish_imm1__spawnednum <- fish_imm1__spawnednum/sum(fish_imm1__spawnednum)
     * sum(fish_mat__offspringnum) * 0.5
fish_imm2__spawnednum <- fish_imm2__spawnednum/sum(fish_imm2__spawnednum)
     * sum(fish_mat__offspringnum) * 0.5

...but I need to dig around in the gadget2 source to see if I'm missing something.

I agree that it may be useful to have the area specified, but I'm trying to think what we want it to represent.

"spawning in different areas and recruits appearing in each of them" feels like the option of least-surprise to me. But I agree it's not something that we can universally apply. We should leave this as a separate issue for the time being.

lentinj commented 5 months ago

We might have multiple spawning events

i.e. multiple mature stock spawning into a single immature stock at the same timestep?

bthe commented 5 months ago

We might have multiple spawning events

i.e. multiple mature stock spawning into a single immature stock at the same timestep?

Yes, but I could be confusing the order of calculations here. I just wanted to make sure that spawning events would not overwrite one another.

bthe commented 5 months ago

My suggestion would be to normalise the expression directly in https://github.com/gadget-framework/gadget3/blob/master/R/action_spawn.R#L153, i.e.

exp(-(((output_stock__midlen - (mean_f)) * (1 / (stddev_f))) ** 2) * 0.5)

by sqrt(2*pi*stddev_f^2) (that is the normalising constant in the Gaussian density function). Then what will follow becomes a bit simpler, something like

fish_imm1__spawnednum <- fish_imm1__spawnednum* sum(fish_mat__offspringnum) * 0.5
fish_imm2__spawnednum <- fish_imm2__spawnednum* sum(fish_mat__offspringnum) * 0.5
lentinj commented 5 months ago

Yeah, I've misinterpreted gadget2. The relevant lump of code is here: https://github.com/gadget-framework/gadget2/blob/master/src/spawner.cc#L316-L347

calcRecruitNumber() works out the equivalent of __offspringnum in the above, i.e. the total number of recruits for the current area.

Storage is an array of area, minimum child age (spawnAge) & union of child length groups. We fill it with the specified stock structure then normalise it to include calcRecruitNumber recruits.

Finally, we split up the contents of Storage into the output stocks by ratio.

I think I prefer how things are in gadget3; defining the recruit stock structure with a "union stock" would create odd artifacts in the smaller of the 2 (although I'm guessing this basically doesn't happen). It'd also be pretty easy to expand gadget3 to support different mean_f for each output stock, which seems like a fairly obvious thing to want.

Also, since we talked about it, gadget2 considers each area in isolation; parents spawn their offspring into the area they are currently in.

lentinj commented 5 months ago

@vbartolino The above commit should sort this out, and your offspring should be split evenly now. If this all works for you I'll merge it into the main branch.

vbartolino commented 5 months ago

@vbartolino The above commit should sort this out, and your offspring should be split evenly now. If this all works for you I'll merge it into the main branch.

@lentinj just tested also varying the proportions, it works as intended. I think it can be merged