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 [2024/07/01 12:39] (current) – [Example: Morris Screening using the sensitivity package] thomas
@@ Line 3: / Line 3: @@
 ===== 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 105: / Line 109: @@
 </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 126: / Line 130: @@
 ===== 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|}}