Group-level and species-level biomass and number of individuals don't add up

[dschlaep]

In the bmassavg.csv output file, rgroup biomass (i.e. p.cool.grass) is not perfectly equal to pssp biomass (despite being the only p.cool.grass species on). The number of individuals that died according to their age is also not identical in a given year between the functional group and only species representing that functional group. This does not make sense.

[kpalmqui]

Yes, this is related to #111, but perhaps was not described well there.

[kpalmqui]

Note also, this problem seems to disappear when _extra_growth is off. At least the discrepancy with relsizes.

[dschlaep]

Above plot was misleading because it was averaged across repeated simulations; here the same plot from one repetition (with seed = 3):

The discrepancies are of one kind: group level biomass is reported as 0 whereas species biomass is repeated at a fixed value.

[dschlaep]

------------------------Repetition/year = 1 / 93
'main' after 'rgroup_Establish': 
	p.cool.grass, relsize = 5.97, est_count = 1
	pssp, relsize = 5.97, est_count = 2
'main' after 'rgroup_PartResources': 
	p.cool.grass, relsize = 5.97, est_count = 1
	pssp, relsize = 5.97, est_count = 2
'main' after 'rgroup_Grow': 
	p.cool.grass, relsize = 6.17, est_count = 1
	pssp, relsize = 6.17, est_count = 2
'main' after 'mort_Main': 
	p.cool.grass, relsize = 4.46, est_count = 0
	pssp, relsize = 4.46, est_count = 0

[dschlaep]

@kpalmqui: This discrepancy arises because the function mort_Main does update est_count of both species and groups, but it does not set relsize to 0.

The output function output_Bmass_Yearly outputs biomass for groups and species (slightly) differently:

• for groups, it calls RGroup_GetBiomass which is a loop ForEachEstSpp summing over species' relsize * max_biomass
• for species, it calls Species_GetBiomass which simply returns that species' relsize * max_biomass

Now, the macro ForEachEstSpp accounts for RGroup[g]->est_count, i.e., it sees that mort_Main set est_count to 0 whereas Species_GetBiomass does not check est_count and reports whatever value remains in relsize.

[dschlaep]

The question now is why does thge function mort_Main set est_count to 0 of both species and groups while not simultaneously setting relsize to 0.

[dschlaep]

The call to _age_independent by mort_Main starts the discrepancy between est_count and relsize:

------------------------Repetition/year = 1 / 93
'mort_Main' after '_no_resources': 
	p.cool.grass, relsize = 6.17, est_count = 1
	pssp, relsize = 6.17, est_count = 2
'mort_Main' after '_age_independent': 
	p.cool.grass, relsize = 4.45, est_count = 0
	pssp, relsize = 4.45, est_count = 0
'mort_Main' after '_slow_growth': 
	p.cool.grass, relsize = 4.45, est_count = 0
	pssp, relsize = 4.45, est_count = 0

[kpalmqui]

Nice find. I suppose I should have reviewed the mortality code more carefuly, given I thought the issue was related to improper killing. Thank you!

I will look more carefully this morning, but my perspective is that why mort_Main does not set relsize to 0 is an oversight that has been missed in the code so far.

[dschlaep]

The species pssp has two individuals in year 93:
The function indiv_Kill_Complete called in a loop over individuals by _age_independent

• reduces relsize of the first individual (relsize drops from 6.17 to 5.41), but does not kill it (est_count stays at 2)
• kills the second individual (est_count drops from 2 to 1 and relsize is reduced from 5.41 to 4.46)

So far so good, but the problem is that upon returning to _age_independent the number of individuals (est_count) miraculously dropped to 0.

------------------------Repetition/year = 1 / 93
'_age_independent' calling 'indiv_Kill_Complete': n=0 <= k=1
'indiv_Kill_Complete' after 'species_Update_Kills': 
	p.cool.grass, relsize = 6.17, est_count = 1
	pssp, relsize = 6.17, est_count = 2
'indiv_Kill_Complete' after 'Species_Update_Newsize': 
	p.cool.grass, relsize = 5.41, est_count = 1
	pssp, relsize = 5.41, est_count = 2
'_age_independent' calling 'indiv_Kill_Complete': n=1 <= k=1
'indiv_Kill_Complete' after 'species_Update_Kills': 
	p.cool.grass, relsize = 5.41, est_count = 1
	pssp, relsize = 5.41, est_count = 1
'indiv_Kill_Complete' after 'Species_Update_Newsize': 
	p.cool.grass, relsize = 4.46, est_count = 1
	pssp, relsize = 4.46, est_count = 1
'_age_independent' after 'indiv_Kill_Complete': 
	p.cool.grass, relsize = 4.46, est_count = 0
	pssp, relsize = 4.46, est_count = 0

[dschlaep]

It seems relsize should not be 0 because est_count should be 1 instead of 0...? Maybe a double count of killed individuals?

[dschlaep]

Oh, no. I may be incorrect in my interpretation in #126 (comment)

Maybe the est_count is correct because the update (i.e., subtraction of 1 = removal of one individual) occurs in the function _delete [rather cryptically hidden in a if statement that checks whether that was the last individual of the species: --s->est_count]. So, indeed two individuals get killed and the est_count is correct, but relsize is incorrect (albeit reduced, but not reduced to a correct 0).

[dschlaep]

The relsize of individuals are subtracted from relsize of a species, but they don't add up:

'_age_independent': individuals of p.cool.grass/pssp with relsize = 6.20
	0, relsize = 0.76
	1, relsize = 0.96
'_age_independent': individuals of p.cool.grass/pssp slated to kill
	0, relsize = 0.76
	1, relsize = 0.96
'_age_independent' calling 'indiv_Kill_Complete': n=0 <= k=1
'indiv_Kill_Complete' after 'species_Update_Kills': 
	p.cool.grass, relsize = 6.20, est_count = 1
	pssp, relsize = 6.20, ndv->relsize = 0.76, est_count = 2
'indiv_Kill_Complete' after 'Species_Update_Newsize': 
	p.cool.grass, relsize = 5.45, est_count = 1
	pssp, relsize = 5.45, est_count = 2
'_age_independent' calling 'indiv_Kill_Complete': n=1 <= k=1
'indiv_Kill_Complete' after 'species_Update_Kills': 
	p.cool.grass, relsize = 5.45, est_count = 1
	pssp, relsize = 5.45, ndv->relsize = 0.96, est_count = 1
'indiv_Kill_Complete' after 'Species_Update_Newsize': 
	p.cool.grass, relsize = 4.49, est_count = 1
	pssp, relsize = 4.49, est_count = 1
'_age_independent' after 'indiv_Kill_Complete': 
	p.cool.grass, relsize = 4.49, est_count = 0
	pssp, relsize = 4.49, est_count = 0

[dschlaep]

Sum of individual relsizes diverges for pssp and its grop p.cool.grass from year 1 step 1:

------------------------Repetition/year = 1 / 1
'main' after 'rgroup_Grow' (1:1): SP: "pssp" size error: SP=0.304738224, ndv=0.230499297
'main' after 'rgroup_Grow' (1:1): RG "p.cool.grass" size error: RG=0.304738224, ndv=0.230499297

At the end of year 100, some species/groups are in sync, others are off; and relsize of groups doesn't get set to 0 if there are no individuals of species of that group:

------------------------Repetition/year = 1 / 100
'main' after '_kill_extra_growth' (1:100): SP: "artr" size error: SP=0.485884070, ndv=0.485884070
'main' after '_kill_extra_growth' (1:100): RG "sagebrush" size error: RG=0.485884070, ndv=0.485884070
'main' after '_kill_extra_growth' (1:100): RG "a.cool.forb" size error: RG=0.000000000, ndv=0.000000000
'main' after '_kill_extra_growth' (1:100): RG "a.warm.forb" size error: RG=0.000000000, ndv=0.000000000
'main' after '_kill_extra_growth' (1:100): RG "p.cool.forb" size error: RG=0.228034645, ndv=0.000000000
'main' after '_kill_extra_growth' (1:100): RG "p.warm.forb" size error: RG=0.001644245, ndv=0.000000000
'main' after '_kill_extra_growth' (1:100): RG "a.cool.grass" size error: RG=0.000000000, ndv=0.000000000
'main' after '_kill_extra_growth' (1:100): RG "p.cool.grass" size error: RG=4.489522934, ndv=0.000000000
'main' after '_kill_extra_growth' (1:100): RG "p.warm.grass" size error: RG=0.423582971, ndv=0.000000000
'main' after '_kill_extra_growth' (1:100): SP: "chvi" size error: SP=0.123962320, ndv=0.019512277
'main' after '_kill_extra_growth' (1:100): RG "shrub" size error: RG=0.123962320, ndv=0.019512277
'main' after '_kill_extra_growth' (1:100): SP: "oppo" size error: SP=0.370149016, ndv=0.370149016
'main' after '_kill_extra_growth' (1:100): RG "succulents" size error: RG=0.370149016, ndv=0.370149016

[kpalmqui]

In regard to the relsize of individuals not adding up to species relsize - that could occur if extra_growth happens in the current year. This extra growth represents increase in biomass that will ultimately be killed at the end of the year. Since it is just killed off at the end of the year anyway, there is no need to update the ndv->relsizes. Instead, the species relsize and rgroup relsizes are updated. This might explain the discrepancy you are seeing (at least before _kill_extra_growth occurs).

Thanks for all of this debugging. So problems so far that required fixing:

1. Making sure relsizes are reduced in _no_resources when individuals are killed (both species and group relsizes). If this is implemented properly, it may take care of the above problem you document in _kill_extra_growth?

[dschlaep]

Hm, ok; I didn't know that about "extra resources".

I don't believe that _no_resources causes the problem because it is not in effect in my case. It is the call to _age_independent that does the killing in year 93 of pssp individuals.

The relsizes are already off before entering _no_resources. The function _age_independent does update the relsizes (but incorrectly so because they don't add up):

'main' after 'rgroup_Grow' (1:93): SP: "pssp" size error: SP=6.203722000, ndv=1.714199066
'main' after 'rgroup_Grow' (1:93): RG "p.cool.grass" size error: RG=6.203722000, ndv=1.714199066
'mort_Main' after '_no_resources': 
	p.cool.grass, relsize = 6.20, est_count = 1
	pssp, relsize = 6.20, est_count = 2
_age_independent': individuals of p.cool.grass/pssp with relsize = 6.20
	0, relsize = 0.76
	1, relsize = 0.96
'_age_independent' after 'indiv_Kill_Complete': 
	p.cool.grass, relsize = 4.49, est_count = 0
	pssp, relsize = 4.49, est_count = 0

[dschlaep]

And if the reasoning about "extra resources" is correct, then this means that pssp has 4.49 units of extra resources while only 1.71 units of "normal" resources in year 93. Does this seem right?

[dschlaep]

And even at the end of a year -- after "extra resources" are removed, there are discrepancies between species/group relsize and the sum of individuals' relsize (both when individuals are present and when all individuals are killed), e.g.,

'main' after '_kill_extra_growth' (1:93): RG "p.cool.grass" size error: RG=4.489522934, ndv=0.000000000
'main' after '_kill_extra_growth' (1:93): SP: "spcr" size error: SP=1.372097373, ndv=1.003715992
'main' after '_kill_extra_growth' (1:93): RG "p.warm.grass" size error: RG=1.372097373, ndv=1.003715992
'main' after '_kill_extra_growth' (1:93): RG "shrub" size error: RG=0.104450047, ndv=0.000000000

[kpalmqui]

I am not sure what you mean by they are off going into _no_resources. They look ok after rgroup_Grow in the above print statements...?

You are right that there is a clearly a problem if ndv = 0.0, but RG or sp relsize > 0.

[dschlaep]

I refer to that pssp relsize is 6.20 as is p.cool.grass relsize before and after_no_resources is called, whereas the sum of the individuals' relsize is only 1.71.

[kpalmqui]

I think that is due to _extra_growth.

[dschlaep]

How do you interpret that these relsizes are still off at the end of a year after extra growth has been removed? For instance,

'main' after '_kill_extra_growth' (1:1): SP: "pssp" size error: SP=0.304738224, ndv=0.230499297
'main' after '_kill_extra_growth' (1:1): RG "p.cool.grass" size error: RG=0.304738224, ndv=0.230499297

[kpalmqui]

In regards to:

And if the reasoning about "extra resources" is correct, then this means that pssp has 4.49 units of extra resources while only 1.71 units of "normal" resources in year 93. Does this seem right?

I have just made extra resources and _extra_growth functional again - it hadn't been functional ever when SOILWAT2 was running. I've been too busy trying to track down these bugs to think carefully about the magnitude of extra resources. 4.49 of extra resources does seem a bit high, but the weather data that is used here is for a wet site. Certainly deserves more thought and evaluation.

[kpalmqui]

In the above case - is there only 1 individual established for each functional group? If not, the species and rgroup relsize will not match the ndv->relsize. Ndv->relsize is 0 to 1, while species->relsize can be >1, as can rgroup->relsize.

[kpalmqui]

If there is only 1 individual, something is clearly wrong.

[dschlaep]

the above output is from

			_kill_extra_growth();
       check_sizes("'main' after '_kill_extra_growth'");
		} /* end model run for this year*/

The existing function check_sizes sums the relsizes of all inviduals of a species (compares that against the relsize of that species), sums up across all species of a resource group (and compares that against relsize of that group).

[dschlaep]

what check_sizes prints out as ndv= is the appropriate sum independent of number of individuals per species and number of species per resource group.

[kpalmqui]

I see - well yes that doesn't make sense.

[kpalmqui]

Can you push up the additional printf and check_sizes statements you have added to you local clone? Thanks!

[dschlaep]

Sure, it's pretty verbose. Here they come

[kpalmqui]

thanks much - no need to me to waste time doing the same thing.

[kpalmqui]

eek getting a segfault with two latest commits

[kpalmqui]

I can just edit locally, not a biggie

[dschlaep]

Do you get a segfault also with commit 6072ad6?

[kpalmqui]

No, 6072ad6 is fine.

[kpalmqui]

So it looks like something might be going wrong in rgroup_Establish for perennial species with individual and species relsizes that is then propagating throughout the year.

------------------------Repetition/year = 1 / 2
WARNING: 'main' after rgroup_Establish (1:2): SP: "pssp" size error: SP=0.299066186, ndv=0.230499297

WARNING: 'main' after rgroup_Establish (1:2): RG "p.cool.grass" size error: RG=0.299066186, ndv=0.230499297

[kpalmqui]

Nevermind, that error is occurring because something went wrong in the previous year and no establishment occurred for pssp in the current year

[kpalmqui]

I think I may have found the problem?

g->res_extra entering _kill_extra_growth is always 0.0, despite it being non-zero in _res_part_extra.

So the difference between the summed ndv->relsize and and species->relsize in the first year appears to be s->extragrowth

[kpalmqui]

Removing:
if (ZRO(g->res_extra))
continue;

At the beginning of _kill_extra_growth solves that problem. However, warnings are still getting thrown after kill_extra_growth:

WARNING: 'main' after kill_extra_growth (1:2): RG "p.warm.grass" size error: RG=0.046365350, ndv=0.000000000

What looks like might be happening here is there are no individuals established anymore (presumably they were killed) and relsize is off because extra_growth is being added back to species-relsize here in the _kill_extra_growth function.

s->extragrowth = 0.046365
, Species = spcr

[kpalmqui]

These small values propagate across years.

Because the _kill_extra_growth loop is set up as ForEachEstSpp(sp, rg, j) { }, when individuals of a species are no longer established, the extragrowth is not removed from the s->relsize because that individual has been killed.

[kpalmqui]

Changing this loop from ForEachEstSpp to ForEachGroupSpp (to loop through all species, skip those where extragrowth=0.0, and do the necessary removal) causes new warnings to the logfile and a logerror:

rgroup_PartResources': Group = sagebrush(0) error with res (relsize = 0.010851, pr = 0.086675):
before correction: res_avail = 0.000000, res_required = 4.340290
after correction: res_avail = 1.000000, res_required = 2.000000
@ generic.c LogError
logfile.log

So that does not seem to be the simple, easy solution to this problem.

[kpalmqui]

This issue has been partially addressed by commits:
c516681
67ff945

In most years, the species biomass and relsize matches the rgroup biomass relsize now. The exception is when individuals belonging to a species are killed, so species->est_count = 0, but species->relsize is not 0 at the end of the year because extra_growth occurred for the species before all the individuals were killed. This does not happen for the rgroup, because there is a check in RGroup_GetBiomass

if (RGroup[rg]->est_count == 0)
return 0.0;

Potential solution is to add similiar logic to Species_GetBiomass.
if (Species[sp]->est_count == 0)
return 0.0;

[kpalmqui]

Resolved by above commit.