Re-evaluate calculation of Rgroup relsize via getRGroupRelsize

[kpalmqui]

We track individual relsizes, species relsizes and rgroup relsizes.

individual relsize: between 0 and 1, where 1 represents a full-sized individual of the species.

species relsize: summed individual relsizes for all individuals belonging to the species, can exceed 1.

rgroup relsize: via getRGroupRelsize currently calculated as,

ForEachEstSpp( sp, rg, n){
		sum += getSpeciesRelsize(sp);
    }

    if(RGroup[rg]->est_count > 0){
        return (RealF) sum / (RealF) RGroup[rg]->est_count;

where getSpeciesRelsize

RealF getSpeciesRelsize(SppIndex sp)
{
	IndivType *p = Species[sp]->IndvHead;
    double sum = 0;

	if(p)
	{
    	while(p)
    	{
        	sum += p->relsize;
			p = p->Next;
    	}
	}

	return (RealF) (sum + Species[sp]->extragrowth);
}

The only major place getRGroupRelsize is called within the code is within RGroup_Update_GrpResProp to calculate grp_res_prop, which is then utilized in _extra_growth.

I am opening this issue because I am concerned this is not the desired calculation of Rgroup relsize. When there is only 1 species established per rgroup, species relsizes are equivalent to rgroup relsizes and we are good to go. However, when more than one species per rgroup is established, we are understimating rgroup relsize and therefore overestimating the relative contribution of each individual to Rgroup relsize (i.e. grp_res_prop). I do not understand why Rgroup relsize is calculated this way.

@dschlaep what do you think? Should we potentially change the calculation of Rgroup relsize to:

ForEachEstSpp( sp, rg, n){
		sum += getSpeciesRelsize(sp);
    }

[dschlaep]

_extra_growth() calculates "increment in size due to res_extra for each individual" in units [(g / plant) / (g / plant)] as:

  extra_ndv = ndv->res_extra * g->xgrow * indivpergram;

where

• indivpergram = 1.0 / s->mature_biomass with units [plant / g],
• g->xgrow is a user input with unitless values in [0, 1] representing ability to capture extra resources,
• ndv->res_extra = ndv->grp_res_prop * g->res_extra.

The "contribution of each indiv to the group's size" is calculated by function getRGroupRelsize() as

  ndv->grp_res_prop = ndv->relsize / sumsize

where

• sumsize = getRGroupRelsize(rg) in questionable units
• and where individual relsize is equivalent to ndv->relsize = ndv->biomass / s->mature_biomass with units [(g / plant) / (g / plant)]

The group's extra resources g->res_extra in units [sum(g / plant)] were previously calculated by function _res_part_extra() that was called by rgroup_ResPartIndiv() that was called by rgroup_PartResources() that was called in the main loop immediately after Env_Generate() but before rgroup_Grow():

  g->res_extra = req_prop * extra

where

• req_prop represents the "proportional biomass of each group out of the total biomass" in units [sum(g / plant) / sum(g / plant)],
• extra represents the extra resources of the group in units [sum(g / plant)]

Thus,

  extra_ndv = ndv->biomass / s->mature_biomass / getRGroupRelsize(rg) * g->res_extra * g->xgrow / s->mature_biomass

and in units

  [(g / plant) / (g / plant)] = [g / plant] / [g / plant] / [getRGroupRelsize(rg)] * [sum(g / plant)] * [-] / [g / plant]

We solve for getRGroupRelsize(rg):

  [getRGroupRelsize(rg)] = ndv->biomass / s->mature_biomass / extra_ndv * g->res_extra * g->xgrow / s->mature_biomass

and in units

  [getRGroupRelsize(rg)] = [g / plant] / [g / plant] / [(g / plant) / (g / plant)] * [sum(g / plant)] * [-] / [g / plant]

which simplifies to

  [getRGroupRelsize(rg)] = [sum(g / plant) / (g / plant)]

The current function getRGroupRelsize() calculates, however,

  [getRGroupRelsize(rg)] = [sum{(g / plant) / (g / plant)} / plant]

The division by number of plants was added to the code on Feb 22, 2019 with commit 068923d.
Unfortunately, neither the commit message nor the issue #232 describe why that division was added.

Therefore, I agree with @kpalmqui's suggested solution which has getRGroupRelsize() calculate as before

  [getRGroupRelsize(rg)] = [sum{(g / plant) / (g / plant)}]

Funnily, 1st time, the division by number of plants was added to the code on Feb 22, 2019 with commit 068923d.
Unfortunately, neither the commit message nor the issue #232 describe why that division was added; however, it appears that this
division was added in response to question 202c74b#commitcomment-32402683
which pointed out that previously the code divided, i.e.,

  //For calculating rgroup relSize, sumsize should be divide by no of current established species in rgroup rather than total no of species in rgroup.
  RGroup[rg]->relsize = sumsize / (RealF) RGroup[rg]->est_count;

This was from a previous version from Oct 22, 2015 commit deed157 which did

  RGroup[rg]->relsize = sumsize / (RealF) RGroup[rg]->max_spp;

which again was changed from a previous version from Oct 21, 2015 commit 8fd0a76 which had

  //For calculating rgroup relSize, sumsize should be divide by no of current established species in rgroup rather than total no of species in rgroup.
  RGroup[rg]->relsize = sumsize / (RealF) RGroup[rg]->est_count;

This in turn was changed on Oct 20, 2015 commit 5a043cc which had

  RGroup[rg]->relsize = sumsize / (RealF) RGroup[rg]->max_spp;

which was the state when the code was first added to github on Dec 5, 2013 commit 4a2ec31.

Funnily, 2nd time, there is a function check_sizes(), which is never used, but it tests whether getRGroupRelsize(rg) is calculated as the sum of getSpeciesRelsize(sp) -- which is currently not the case, but would do so with @kpalmqui's solution!

[kpalmqui]

@dschlaep thanks for the detailed summary of how this has been handled since 2013!

A bit more backstory:

You are right that originally the code divided sumsize by RGroup[rg]->max_spp. I deemed this as incorrect and changed it to RGroup[rg]->est_count because RGroup[rg]->max_spp is fixed during runs and represents the maximum number of species that can be established in a RGroup, not the actual number of species established in the RGroup, which is represented by RGroup[rg]->est_count.

Chandler replaced all occurrences of RGroup[rg]->relsize = sumsize / (RealF) RGroup[rg]->est_count in the Data_Precision_Improvments branch/milestone via getRGroupRelsize. But he did not change the code logic.

So despite changing RGroup[rg]->max_spp to RGroup[rg]->est_count, Rgroup relsize has always been calculated by dividing by the number of species, rather than just summing the Species relsizes.

It sounds like you are on board with updating to:

ForEachEstSpp( sp, rg, n){
		sum += getSpeciesRelsize(sp);
    }

Perhaps I will open a new branch for this, unless you think I should just fold this into feature_space? This isn't impacting my runs at all because I am only using 1 species to represent each RGroup.

[dschlaep]

@kpalmqui - I am fine both with a new branch or with just a commit directly to feature_space.

[kpalmqui]

@dschlaep as expected, this change resulted in a slight reduction in STEPPE biomass for RGroups where I turned on two species, presumably due to a reduction in extragrowth. This was the case with the default sagebrush parameters. It will be interesting to see if this has any impact with the Swiss grassland inputs. I do believe you looked at this before, but it didn't really matter?

[kpalmqui]

closed by 43f180f