====== 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 [[http://134.76.18.36/wordpress/courses-and-tutorials/simplefspm/ | "Example08" FSPM model]] from the gallery as an example.
===== Prerequisites =====
Make sure to [[:tutorials:getting-started-with-grolink-and-gror|set up GroR]]. Play a bit around with it and get a feel for how it works and what the [[https://gitlab.com/grogra/groimp-utils/rapilibrary#grolink|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.

==== Downloads ====
  * {{ :tutorials:example08_prepared.zip | Prepared Example08 model}}
  * {{ :tutorials:example08_analysis.zip | R Script of this wiki and model outputs as .rds}}

===== Prepare your model=====
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.

==== Parameter File ====
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:

<code java>
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;}
				}
	...
]

</code>
The parameter values are defined in a seperate file, ''parameters.rgg'':

{{:tutorials:params.png?300|parameters.rgg file}}

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):

<code java>
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;
</code>

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

<code java>
import parameters.*;
</code>

==== POI definition and output communication ====
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:

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

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.





===== Model execution and output gathering in R=====
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:

<code R>
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

</code>

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

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

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.

<code R>
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)
  }
}
</code>







===== Example: Morris Screening using the sensitivity package =====
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:

<code R>
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
</code>

''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:

<code R>
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)
}
</code>

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:

<code R>
plan(multisession, workers = availableCores())

system.time(model_outputs <- future_apply(params, 1, executeModel))

saveRDS(model_outputs, "example08_morris_outputs.rds")
</code>

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:

<code R>
sensitivity::tell(m, model_outputs)

plot(m)
</code>

This should now look something like this:

{{:tutorials:rplot.png?500|}}