Differences

This shows you the differences between two versions of the page.

--- tutorials:sensitivity-analysis-using-grolink-and-gror [2024/06/25 15:29] – [Model execution and output gathering in R] thomas
+++ tutorials:sensitivity-analysis-using-grolink-and-gror [2025/06/06 17:11] (current) – [Example: Morris Screening using the sensitivity package] t
@@ Line 1: / Line 1: @@
 ====== 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.
+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 [[https://gallery.grogra.de/project_gallery/67053a3f256715cda41c867f | "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=====
@@ Line 47: / Line 51: @@
 {{: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 folder and file can be created via **Object** -> **New** in the File explorer.
 The file only contains the parameter definitions (the values are what they were in the original gallery model):
@@ Line 94: / Line 98: @@
 library(sensitivity)
-wb1 <- GroLink.open("http://localhost:58081/api", path="Example08.gsz") # copy gsz to groimp path
+# copy groimp model gsz to groimp path to open it like this:
+wb1 <- GroLink.open("http://localhost:58081/api", path="Example08_prepared.gsz")
+# push model gsz to groimp from any file location like this:
+wb1 <- GroLink.open("http://localhost:58081/api",content = readBin("D:/Example08_prepared.gsz", "raw", 10e6))
 </code>
@@ Line 105: / Line 112: @@
 </code>
-Now here is the function 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.
+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>
@@ Line 119: / Line 126: @@
 </code>
+When you run the while loop, in the end, the ''model_output'' variable will contain the total amount of absorbed light by all leaves at the growth step when flowers emerged (this will be some float) or 0 in case no flowers emerged.
@@ Line 126: / Line 133: @@
 ===== 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,       # Vector containing values for all model parameters
+                         timeout = 200 # max. number of grow function executions
+                         ){
+  # Open workbench with model
+  wb1 <- GroLink.open("http://localhost:58081/api", content = readBin("D:/Example08_prepared.gsz", "raw", 10e6))
+  # Override parameters.rgg with new values from params vector
+  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 = ""))
+  # Compile the workbench
+  WBRef.compile(wb1)
+  model_output <- ""
+  n_grows <- 0
+  # Execute the grow function for as long as there are no flowers
+  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
+  }
+  # Retrun NA in case of timeout
+  if (n_grows == timeout) {
+    model_output <- NA
+  }
+  # Close workbench in the end
+  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|}}