Sensitivity analysis on GroIMP models using GroR

This wiki explains how to do a sensitivity analysis on GroIMP models using the GroR interface using a Morris screening over input parameters of the "Example08" FSPM model from the gallery as an example.

Make sure to set up GroR. Play a bit around with it and get a feel for how it works and what the different functions do and what they return. The approach presented in the wiki here is just one way you could approach a sensitivity analysis. You will likely find your own approach for your specific model, but GroR will be the space you work in.

• Prepared Example08 model
• R Script of this wiki and model outputs as .rds

In any sensitivity analysis, you will need three things:

• Some way of pushing parameter settings to the model
• Some way of knowing that your model has finished or reached a point of interest
• Some way of grabbing the model output at that point of interest

The approach presented here is built on the idea of being able to do all of these things from within R. Therefore, your model will likely have to be adapted so you can feed all of this information easily through GroR.

Parameters will be pushed to the model by modifying a special RGG file which contains only the Parameter definitions, as this is easy to do in GroR. In the Example08 FSPM, I decided to change five hardcoded parameters to variables in the run() and la() methods:

protected void run ()
[ 	Bud(r,p,o),(r<10 && p>0) ==> Bud(r,p-1,o);
	Bud(r,p,o),(r<10 && p==0 && o<4) ==> 
				RV(-0.1) Internode(parameters.NormalInternodeLength,1) NiceNode [ RL(50) Bud(r+1,phyllo+irandom(-5,5),o+1) ] 
				[RL(70) Leaf(parameters.LeafLength, 0.07,0,1,0)] RH(137) RV(-parameters.PlantWideness) NiceInternode Bud(r+1,phyllo+irandom(-5,5), o);
	Bud(r,p,o), (r==10) ==> RV(-0.1) Internode(parameters.FlowerInternodeLength,1) Internode(parameters.FlowerInternodeLength,1) NiceFlower(1);
 
	nf:NiceFlower ::> nf[age]++;
	nf:NiceFlower, (nf[age]>irandom(10,15)) ==> ;
]
 
protected void la ()
[
	lf:Leaf ::> {
		lf[al] = lm.getAbsorbedPower3d(lf).integrate()*2.25;
		lf.(setShader(new AlgorithmSwitchShader(new RGBAShader((float) lf[al]/5.0f, (float) lf[al]*2, (float) lf[al]/100.0f), GREEN)));
		lf[as] = (lf[al]*86400*2)/1000000.0f;
		lf[age]++;
		float lfas = sum((* Leaf *)[as]);
		if (lfas>0) {lf[length] += logistic(2,lf[age],10,0.5);
		lf[width] = lf[length]*parameters.LeafAspectRatio;}
				}
	...
]

The parameter values are defined in a seperate file, parameters.rgg:

The folder and file can be created via <key>Object</key> → <key>New</key> in the File explorer. The file only contains the parameter definitions (the values are what they were in the original gallery model):

static float NormalInternodeLength = 0.1;
static float FlowerInternodeLength = 0.05;
static float LeafLength = 0.1;
static float LeafAspectRatio = 0.7; 
static float PlantWideness = 0.1;

parameters.rgg needs to be imported in your main model file (in this case test.rgg) at the top of the file:

import parameters.*;

The Example08 model in its default way of existing grows a plant and lets it flower. Obvious interesting outputs are the sum of absorbed light by the leaves and the amount of produced assimilates. However, this model has no defined end, it could run forever (even tho nothing really happens in the later stages). For this example, I decided that the interesting output would be the total absorbed light by the leaves at the point when the first flower emerges. Because you can conveniently grab the output on the console in GroR, I added some code in the grow() method to communicate both the existence of flowers and the amount of absorbed light:

if(count((*NiceFlower*)) > 0){
		println(sum((* Leaf *)[al]));
	} else {
		println("no flower");
	}

So, after every grow() iteration, the model prints no flower as long as there are no flowers and a number for the amount of absorbed light once there is a flower.

I will first demonstrate the concept by gathering the output (the amount of absorbed light) from just one model execution (run).

First, load some libraries and set up the workbench with the prepared model:

source("D:\\groimp\\rapilibrary\\R\\GroR.R") # load GroR
 
library(future)
library(future.apply)
library(sensitivity)
 
wb1 <- GroLink.open("http://localhost:58081/api", path="Example08.gsz") # copy gsz to groimp path

You can make sure that everything works by looking up the available model functions or reading the parameter file:

(functions <- WBRef.listRGGFunctions(wb1))
WBRef.getFile(wb1, "param/parameters.rgg")

Now here is the code that gets the model output. It executes the grow() method until it receives a console output that is not no flower, meaning that there is some number for the amount of light available. When creating these kind of functions, it usually makes sense to define some kind of timeout (here 200 grow executions) because you never know if your model will actually work correctly with some weird parameter combinations you might give to it.

model_output <- ""
n_grows <- 0
while (!(is.numeric(model_output)) && (n_grows < 200)) {
  result <- WBRef.runRGGFunction(wb1,"grow")
  model_output <- unlist(result$console)
  if (model_output != "no flower"){
    model_output <- as.numeric(model_output)
  }
}

The Morris Screening will be used to analyze the 5 structural plant growth parameters with regard to their importance (their main and interaction effect) on the total amount of absorbed light by leaves. This is just a random example out of the plethora of available sensitivity analysis methods. Most of the common ones are implemented in the sensitivity R-package. This is not a tutorial on sensitivity, but it is very well documented (e.g. in the R help function). It is worth noting that this example uses the so-called decoupled approach of sensitivity, see ?tell. This basically means that the generation of the parameter set, model execution and output analysis will be done as separate steps.

First, I generate a set of input parameters for the model:

m <-    morris(model = NULL, factors = c("NormalInternodeLength", # baseline is 0.1
                                         "FlowerInternodeLength", # baseline is 0.05
                                         "LeafLength",            # baseline is 0.1
                                         "LeafAspectRatio",       # baseline is 0.7
                                         "PlantWideness"),        # baseline is 0.1
               r = 20,                                                 # 20 repetitions
               binf = c(0.05, 0.01, 0.05, 0.1, 0.1),                   # min value of inputs
               bsup = c(0.8, 0.1, 1, 1, 1),                            # max value of inputs
               design = list(type = "oat", levels = 6, grid.jump = 3)) # 6 levels per parameter (check with e.g. length(unique(x$X[,4])))
                                                                       # grid.jump is recommended to be levels/2, see ?morris
 
 
params <- m$X

params is now a 120-row matrix. In order to be able to loop over the matrix rows, I need a function that takes one parameter set as input, executes the model and returns the output for this set:

executeModel <- function(params, timeout = 200){
  wb1 <- GroLink.open("http://localhost:58081/api", path="Example08.gsz")
  WBRef.updateFile(wb1, "param/parameters.rgg",
                   paste("static float NormalInternodeLength = ", as.character(params[1]),
                         ";\r\nstatic float FlowerInternodeLength = ", as.character(params[2]),
                         ";\r\nstatic float LeafLength = ", as.character(params[3]),
                         ";\r\nstatic float LeafAspectRatio = ", as.character(params[4]),
                         "; \r\nstatic float PlantWideness = ", as.character(params[5]),
                         ";",
                         sep = ""))
 
  WBRef.compile(wb1)
 
  model_output <- ""
  n_grows <- 0
  while (!(is.numeric(model_output)) && (n_grows < timeout)) {
    result <- WBRef.runRGGFunction(wb1,"grow")
    model_output <- unlist(result$console)
    if (model_output != "no flower"){
      model_output <- as.numeric(model_output)
    }
    n_grows <- n_grows + 1
  }
  if (n_grows == timeout) {
    model_output <- NA
  }
  WBRef.close(wb1)
  return(model_output)
}

This function contains the while-loop from above but also modifies (overrides) the parameters.rgg file. The paste() function is used to create the file content. The function also opens and closes a new workbench with every time it is called, this is needed to be able to parallelize model executing using the future framework:

plan(multisession, workers = availableCores())
 
system.time(model_outputs <- future_apply(params, 1, executeModel))
 
saveRDS(model_outputs, "example08_morris_outputs.rds")

Through the use of future_apply and the way the multisession is setup, you will run one instance of GroIMP on every core of your computer. With this Example08 model, this can take some time, on my 20-core machine it took around 2 minutes (and used ~5Gb RAM). This simple parallelization could likely be optimized a lot and will not scale amazingly to e.g. a remote cluster.

The only thing remaining is to analyze the output and plot the results:

sensitivity::tell(m, model_outputs)
 
plot(m)

This should now look something like this: