Code review request of fecundity function

[ConnectedSystems]

Hi @Rosejoycrocker , please have a look at:

https://github.com/open-AIMS/ADRIA_repo/blob/3f57a89de3cc94de540d769fe7d26c43c4726ae7/ADRIAfunctions/fecundityScope.m#L18

Yes, I can edit :-)

[ConnectedSystems]

Hi @KRNA01 @Rosejoycrocker

Just had a quick look at the declining population issue using all default values (RCP 4.5, etc, etc, etc) with the below changes:

e_mb = coral_params{:, mb_idx}'; %background coral mortality
e_mb(:) = 0.0;  % debug: 0 mortality

% ... snip ...

% competition factor between Small Massives and Acropora
e_comp = sim_params.comp;
e_comp = 0.1;  % debug: very little competition

% ... snip ...

% larval productivity per site
LPs = ADRIA_larvalprod(tstep, assistadapt, natad, past_DHW_stress, ...
    LPdhwcoeff, DHWmaxtot, LPDprm2);
LPs(:) = 1.0;  % debug: maximum larval productivity

% ... snip ...

rec = (Y_pstep * TP_data) .* LPs;
rec(rec > 0) = 1.0;  % debug: unrealistically high recruitment

Only by setting recruitment to 1.0 do we see an increasing population trend for Unenhanced Corymbose Acropora and small/large massives

My gut feeling is that the issue lies in this line:

rec = (Y_pstep * TP_data) .* LPs;

TP_data is the data loaded from MooreTPmean.xlsx, and its values range from: 0.000015135 to 0.14571596

Maximum value in Y_pstep (at least for the first time step) is 0.1214.

Stepping forward in time a bit, it seems recruitment values for small/massives are always very low.

I made a little animation to demonstrate, after removing the debug changes shown above: (first six species are the Enhanced Tabular Acropora, last 12 are the small/large massives)

[Rosejoycrocker]

Hey @ConnectedSystems Great animation. There does seem to be almost no dynamics for these species. As TP_Data is loaded from MooreTPMean.xlsx, do we have any other data sets we could compare with @KRNA01? Or perhaps simulate from expected ranges...

[Rosejoycrocker]

Hi @KRNA01 , Just wanted to clarify because I'm currently piecing together how the new expanded coral species code works. Where does the coral_param.fec parameter come from? (at the following link) https://github.com/open-AIMS/ADRIA_repo/blob/3f57a89de3cc94de540d769fe7d26c43c4726ae7/ADRIAfunctions/fecundityScope.m#L18

Also is coralrunADRIA the new runADRIA function? (Apologies, I'll try and help with debugging once I've got my head around the new code)

[KRNA01]

Hi @ConnectedSystems and @Rosejoycrocker, Thanks for looking into this. @Rosejoycrocker, I am just about to send you an updated comment with pointers to some updated code. I think I have cracked part of the fecundity problem (you be the judge :-). Would then be great if you can review this as well the growthODE4 function and other related functions. More within the hour ..

[Rosejoycrocker]

Thanks @KRNA01 :) I'll review when I receive your pointers

That's great that you might have solved the fecundity problem!

[KRNA01]

@ConnectedSystems and @Rosejoycrocker,

A quick update on where I think we're at plus a humble request for further assistance and advice.

Firstly, I have changed recruitment in 'coralScenario()' and 'growthODE()' such that it now pertains to each of the six coral groups, i.e. it has dimensions ngroups by nsites (6 by 26). Previously recruitment had size 36 by 26, i.e for every size class, which does not make sense - recruitment is by definition the addition of new corals to the smallest size class for each group.

https://github.com/open-AIMS/ADRIA_repo/blob/3a3d4f07875972175cb52e2278cfe81098a4a782/ADRIAfunctions/coralScenario.m#L182

You'll see that there are now more steps added prior to the estimation of recruitment, including the ratio of estimated larval density (in part determined by spawning at the different sites, resulting larval density, and resulting number of settlers actually (potentially) contributing to recruitment. This was missing before, and represents about a factor 1000 in the under-estimate of recruitment - my bad :-(.

@Rosejoycrocker, for now, I have disabled the larval productivity function ('ADRIA_larvalprod()' because it's dimensions are wrong. I would appreciate your help with fixing this. Ideally larval production should have dimensions ngroups by nsites (6 by 26), but is currently 36 by 26. Similar to the 'fecundityScope2()', it should pertain to coral groups, not individual size classes. The larval production function is important because it applies a fecundity penalty to corals that have experienced heat stress in the previous summer.

Reference to relevant paper here: https://www.nature.com/articles/s41586-019-1081-y

Without this penalty function, I suspect recruitment rates following bleaching events are now too high. Here are the results for a single simulation ('single_simulation' from the examples folder in the 40corals branch) with default settings for RCP 4.5. Unenhanced corals now do better than predicted.

Here's the link to the fecundity function:

https://github.com/open-AIMS/ADRIA_repo/blob/3a3d4f07875972175cb52e2278cfe81098a4a782/ADRIAfunctions/ADRIA_larvalprod.m#L1

Lastly, I'd be grateful if you could both review the 'growthODE4()' function. You'll see that I have made a couple of comments.

The most pertinent, I think is this one.
https://github.com/open-AIMS/ADRIA_repo/blob/3a3d4f07875972175cb52e2278cfe81098a4a782/ADRIAfunctions/growthODE4.m#L15

The reason I am flagging this one is that reshaping of arrays (in my limited experience) can upset the row/column order so that you end up with nonsensical data when reshaping. May not be an issue here, and great if it isn't.

@Rosejoycrocker, I apologise in advance for this monstrosity of a function. If you can see ways to simplify, please do. If so, please then rename in a new function file that we can tag when we test and before throwing the old on out.

The most clumsy bit is probably the carrying capacity (P_x), which pertains to coral groups within sites, but is expanded to all 36 groups. I can't see a way around that, but I'm sure you can.

[Rosejoycrocker]

Thanks @KRNA01 , I'll have a look into the dimension issue with the larval productivity function this afternoon. As for the growthODE4 function it does seem that at least the central cases (thinking of it in a finite difference sense) may be able to be done in one operation for each coral group, so I'll have a go at this (just might be messy working it out). In order to test this, is there a script that runs ADRIA using the new functions? I noticed runADRIAscenario still has the 4 groups ode, so is there a script that runs all of ADRIA that you've been using to test the ODEs?

[ConnectedSystems]

@Rosejoycrocker

See the scripts in examples/running_ADRIA

• single_scenario.m
• programmatic_single_scenario.m

The first displays some figures at the end, which I guess could be copied over to the second one

Once we're all happy and this is all fixed up, could we move the constants defined here to simConstants.m?

Unless these parameters are intended to be varied?

[Rosejoycrocker]

Thanks @ConnectedSystems :)

[ConnectedSystems]

Just to clarify:

Ken and I left the runADRIAScenario and runADRIA alone for now - the functions that include representations of the 36 corals are:

runCoralADRIA(), which calls for every scenario coralScenario()

These will replace runADRIA and runADRIAScenario respectively once everything is finalized.

[KRNA01]

Thanks @ConnectedSystems and @Rosejoycrocker,

Happy to move new constants to simConstants.m

I started doing this and got in trouble for two reasons: (1) an error about the dimensions of the table and (2) the dialog box becomes taller than my laptop window :-)

[KRNA01]

@ConnectedSystems, what would be required to update the n-intervention 'run_example.m' on the 40corals branch?

Not urgent, but keen to see how this new 36-species version compares to the old 4-species one across a set of interventions.

Let me know which tasks you need to assign to me.

[Rosejoycrocker]

Hi @KRNA01, I've made a quick fix for the larval prod function. I split the tmp_ad vector into coral groups using reshape (so its a 6 by 6 matrix with groups along the columns), then using the mean function averages along the columns to give a group value for tmp_ad (a 1 by 6 vector). Transposing this in the final calculation gives Y as 6 by 36. I'm not sure how ecologically sound this is using the average for each group so let me know if it needs to be weighted for enhanced vs. unenhanced corals or something (although I think because we start with 36 values the weighting is implicit).

[ConnectedSystems]

@Rosejoycrocker I'm guessing this is the same issue @KRNA01 reported in our catch up yesterday. Do you run into this issue when running single_scenario.m? If so, I just pushed a fix.

The cause was I split/renamed the p variable (an array of 2 values) to gompertz_p1 and gompertz_p2 and forgot to update the example script.

There are three groups of parameters I skip when displaying the dialog box for simplicity (namely those that expect an array of inputs). Because of the above change it skips an incorrect number of parameters.

Tested with a few additional dummy parameters defined in simConstants and it seems to work without issue now.

As for the dialog box extending offscreen, I have the same issue. I've been just pressing "tab" to skip to the "ok" button for now. I've been meaning to ask one of you how to display the dialog options across two columns to save space.

---

@KRNA01

We just need to generate the same intervention and criteria weights and run those through the runADRIA and runCoralADRIA functions. You could cobble together an example looking at run_example and coral_example, the latter being in the sandbox folder (and earmarked to replace run_example when the time comes).

Please note the sandbox directory was my own personal playground for testing out code - I pushed this up by mistake.

[KRNA01]

Thanks @Rosejoycrocker,

Y for the larval productivity function should be a 6 by 26 array, with the second dimension for sites.

Just looked for your suggested edits in the branch, but can't see them. Are you able to point me to the m file please?

[Rosejoycrocker]

Hi @ConnectedSystems, I just pulled and ran single_scenario and ADRIA_larvalprod still outputs a 36 by 26 matrix- are we talking about the same fix? If not I'll push what I've done (which gives the 6 by 26 output)

[Rosejoycrocker]

I just pushed @KRNA01 ... its an easy change if you think it's incorrect, you can just comment out one line in the Larvalprod function... also sorry I meant 6 by 26

[ConnectedSystems]

Hi again @Rosejoycrocker

are we talking about the same fix

I guess not. I thought issue 2 and 1 were related to adding new parameters to simConstants.m in the below but I guess you weren't talking about the parameter table!

I started doing this and got in trouble for two reasons: (1) an error about the dimensions of the table and (2) the dialog box becomes taller than my laptop window :-)

Regarding the larval production/recruitment matrix dimensions, I left a note in coralScenario regarding the discrepancy:

% Note: Matrix format is nspecies * nsites, but the values repeat.
%       Only the first entry for each taxa is intended to be used.
%       Shape of this matrix is simply for convenience.
%rec = (Y_pstep * TP_data) .* LPs;  %old version

Now this may have had other unintended implications but my understanding was that the values produced through these steps are ultimately used in the growthODE* functions, wherein the recruitment values are pulled from specific rows.

In other words, even though a 36*26 matrix is the input, it originally applied only six rows in the calculations (the first class for each taxa).

Note that the indices used for rec increment by 6

%Tabular Acropora Enhanced
Y(1, :) = P_x(1, :) .* rec(1, :) - P_x(2, :) .* X(1, :) .* (r(1)  - X(1,:) .* mb(1));

% ... snip ...

%Tabular Acropora Unenhanced
Y(7, :) = P_x(7, :) .* rec(7, :) - P_x(8, :) .* X(7, :) .* r(7) - X(7,:) .* mb(7);

% ... snip ...

%Corymbose Acropora Enhanced
Y(13, :) = P_x(13, :) .* rec(13, :) - P_x(14, :) .* X(13, :) .* r(13) - X(13,:) .* mb(13);

% ... and so on ...

see change history here

I understand what I did isn't the cleanest approach, it was just the easiest thing I could think of at the time.

If I understand the code correctly, the issue all stems from LPs being in the unexpected shape?

This is caused by assistadapt and natad being defined for each "species" (which I think is correct) so I think the issue you're grappling with could be resolved by defining arrays that hold values for just the smallest size class for each taxa.

In other words, change (circa L 161):

LPs = ADRIA_larvalprod(tstep, assistadapt, natad, past_DHW_stress, ...
            LPdhwcoeff, DHWmaxtot, LPDprm2); % larval productivity ...

to

ADRIA_larvalprod(tstep, assistadapt(1:6:end), natad(1:6:end), ...
    past_DHW_stress, LPdhwcoeff, DHWmaxtot, LPDprm2)

I have a separate branch coral_params in which I'm putting in consideration of coral parameters for each size/taxa so that hardcoded "6" value would be replaced with length(unique(coral_params.class_id)).

I think (emphasis on think) this would solve the issues you're running into, but there are likely other considerations downstream which haven't come out yet.

Testing the above change, I get the following result (default values, run with single_scenario.m):

[Rosejoycrocker]

Thanks @ConnectedSystems . I'm a little confused because I haven't looked at simConstants yet... I'll have a look and try and work it out

[KRNA01]

Hi @ConnectedSystems and @Rosejoycrocker,

I just tested Roses fix and it works, thank you!

Takuya, I get the same results are you and when I'm testing different constellations of inputs, I'm not surprised about what I'm seeing.

Re indices for recruitment (rec) in the growthODE4 function, I have changed these to now take single steps because the input array for recruitment pertains to groups not size classes. @ConnectedSystems, please check that I have not done something that stuffs up the ODE operation here

[ConnectedSystems]

🤦‍♂️ Sorry for the confusion @Rosejoycrocker - I just realised it wasn't you who posted this

I'm going blind or something.

[Rosejoycrocker]

Hi @KRNA01 , Just noticed the index pattern is a bit different on this line of growthode4 for X and r https://github.com/open-AIMS/ADRIA_repo/blob/5e8bd3942d74a270a23e89ecb73cfa370644213e/ADRIAfunctions/growthODE4.m#L78

Is this just part of the ode or a mistake? (apologies if its just part of the ode)

[Rosejoycrocker]

🤦‍♂️ Sorry for the confusion @Rosejoycrocker - I just realised it wasn't you who posted this

I'm going blind or something.

That's all good, I find the structure of Github commenting pretty confusing in terms of tracking issues so I misinterpret things a lot

[KRNA01]

@Rosejoycrocker,

Oops - thank you. I'm going blind too ..

Thank you! Fixing as we speak

[ConnectedSystems]

Hi @ConnectedSystems and @Rosejoycrocker,

I just tested Roses fix and it works, thank you!

Takuya, I get the same results are you and when I'm testing different constellations of inputs, I'm not surprised about what I'm seeing.

Re indices for recruitment (rec) in the growthODE4 function, I have changed these to now take single steps because the input array for recruitment pertains to groups not size classes. @ConnectedSystems, please check that I have not done something that stuffs up the ODE operation here

Just looked at the fix @Rosejoycrocker suggests.

I think my suggestion and Rose's are roughly equivalent (with a caveat I detail below). Obviously I am biased but I prefer my own 😛

As I wrote above, assistadapt and natad specify values for each "species". Not so relevant with assistadapt, but allows us to perturb the natad values across all taxa/size classes for sensitivity analysis or other scenario exploration.

Currently, however, the values are identical for each taxa/size class. In Rose's fix, the mean is taken, which resolves back to the original value.

The caveat here is that the mean may not make sense once we start perturbing natad for individual taxa/sizes. I leave the judgement to @KRNA01.

My own current understanding is that recruitment should only consider the smallest size class, so in my approach I simply pass in the values for the first size class, i.e., with the snippet shown above, and repeated below for convenience.

ADRIA_larvalprod(tstep, assistadapt(1:6:end), natad(1:6:end), ...
    past_DHW_stress, LPdhwcoeff, DHWmaxtot, LPDprm2)

[ConnectedSystems]

the input array for recruitment pertains to groups not size classes

Ah, so it should be the average LP for each taxa rather than just the smallest size class @KRNA01 ?

In which case @Rosejoycrocker suggested fix is more appropriate.

(but yes, in any case the dimensions of the resulting matrix all line up so no issues there)

[KRNA01]

@ConnectedSystems and @Rosejoycrocker,

GitHub now tells me that the growthODE4 function is a conflicted file.

How do we resolve ?

[Rosejoycrocker]

sorry that's probably my fault, I compressed some of the ode onto single lines for ease of reading, can you view conflicts?

[Rosejoycrocker]

I can push the original if that's easier

[KRNA01]

@Rosejoycrocker, Thanks Rose, I've pulled your file :-)

[Rosejoycrocker]

Thanks, sorry @KRNA01 , just did it without thinking :)

[KRNA01]

Not at all, @Rosejoycrocker, nice way to condense the code. And it works 👍

[Rosejoycrocker]

Nice, glad it fits with your code :)

[ConnectedSystems]

The most pertinent, I think is this one.
https://github.com/open-AIMS/ADRIA_repo/blob/3a3d4f07875972175cb52e2278cfe81098a4a782/ADRIAfunctions/growthODE4.m#L15

The reason I am flagging this one is that reshaping of arrays (in my limited experience) can upset the row/column order so that you end up with nonsensical data when reshaping. May not be an issue here, and great if it isn't.

Sorry, @KRNA01 I missed this question initially.

ode45() passes in a single column vector so we have to reshape back to a matrix with the expected dimensions.

This is also why at the end we have to reshape back to single column vector, as otherwise ode45() errors out:

https://github.com/open-AIMS/ADRIA_repo/blob/3a3d4f07875972175cb52e2278cfe81098a4a782/ADRIAfunctions/growthODE4.m#L85

To confirm the transformations don't affect anything I made sure that the reshaped matrix matches the original initial input.

See specification of Yin1 here

global Xtest

Yin1 = ...

global Xtest 
Xtest = Yin1   % initial input in original shape

[~, Y] = ode45(@(t, X) growthODE4(X, e_r, e_P, e_mb, rec, e_comp), tspan, Yin1, non_neg_opt);

%% inside growthODE4()
global Xtest
X = reshape(X, size(rec));

% ensure Xtest == X, otherwise raise an error message "not equal"
assert(isequal(Xtest, X), "not equal")  % this assertion never fails for the initial input
``

[KRNA01]

ode45() passes in a single column vector so we have to reshape back to a matrix with the expected dimensions.
This is also why at the end we have to reshape back to single column vector, as otherwise ode45() errors out:

Thank you, @ConnectedSystems - got it now.

Great that you could test this - and that it works 🥇

[KRNA01]

Hi @ConnectedSystems and @Rosejoycrocker,

Thanks both for your help with this issue. I think we have solved it, so suggest we close. Unless you think there are outstanding bits we need to resolve.

[KRNA01]

@ConnectedSystems, on second thought, please hang off on closing for a day - I'd like to capture a few gems from the discussion first.

[ConnectedSystems]

@KRNA01 yes sure. I've picked up on a few issues too, more related to the metrics though.

[ConnectedSystems]

Hi @KRNA01 @Rosejoycrocker

Completely forgot to raise this in our meeting yesterday but is the fix for larval production acceptable?

The summary of the issue is:

Ken: Previously recruitment had size 36 by 26, i.e for every size class, which does not make sense - recruitment is by definition the addition of new corals to the smallest size class for each group.

My suggested fix is to only pass in the data for smallest class (the hardcoded index could be replaced with the new coralClassIndex() function).

% only the every 6th entry (to be replaced with something more generic)
ADRIA_larvalprod(tstep, assistadapt(1:6:end), natad(1:6:end), ...
    past_DHW_stress, LPdhwcoeff, DHWmaxtot, LPDprm2)

Rose's fix also works as it takes the mean of values (which from memory are identical across size classes), but may cause issues if these are ever varied. Or, it could be that the averaged values incorporates consideration of other processes.

[Rosejoycrocker]

Hey @ConnectedSystems, I wasn't aware of the taxa_id vector when I did this so happy to adopt the more generalised form if you're worried about it.

[KRNA01]

Hi @Rosejoycrocker and @ConnectedSystems ,

Rose, would be great if you could review the 'ADRIA_larvalprod()' function with this in mind (in the 40corals branch). I am in the function now to clean up the documentation, but will be out in a tick.

Both, please note that the larvalprod function should apply to each of the six coral groups individually and to each site, but not to individual size classes. It's okay to collapse (or average) across size classes to get to six groups, but not across groups.

[Rosejoycrocker]

Sure, @KRNA01 , I'll look into it this afternoon :)

[KRNA01]

@Rosejoycrocker and @ConnectedSystems,

Rose, just checked the 'ADRIA_larvalprod()' it works as intended - thank you. We can consider this closed.

[ConnectedSystems]

Hi @Rosejoycrocker

My concern was the averaging more so than the index accessing. But as above, Ken's confirmed it's fine so I will go ahead and close.

Thanks everyone.