The RStudio developers recently released reticulate, an R package for interfacing directly with Python from R. It’s a significant improvement over PythonInR and rPython, and deprecated at least four of my packages basically overnight (including my own attempt at a Python-R interface. I’ve been playing around with reticulate and it works great—the only thing that I feel is missing right now is conversion between Pandas dataframes and R dataframes, and it looks like the reticulate developers plan to add support.

There aren’t many examples for reticulate yet, so that’s what this post is for. UC Davis Professor John Herman developed the Python module SALib (also available on pypi) which implements commonly used sensitivity analysis methods. I’ve been curious about trying the module out for a while, but I didn’t have any existing models I wrote in Python that would be a good fit for playing with it. By using reticulate, I can use SALib with my R code directly. To be fair, many of the same methods are probably included in the sensitivity package for R, but it’s a good excuse make a reticulate example.

Consider a hypothetical trapezoidal channel that is designed to pass a flood wave. The local engineering district is considering design options and wants to know what will have the strongest impact on the wave crest height and travel time to a specific location, subject to certain constraints. The current options available to the district are (1) channel bedslope design, (2) channel bottom width, (3) channel sidewall design, and (4) channel bed material and vegetation removal strategies (reflected in the value of the roughness coefficient).

We can use the rivr package to model the progression of a floodwave in R. We’ll investigate the sensitivty of the flood wave crest height and travel time by varying the four parameters within the bounds shown in the table below.

I’ll use SALib to do a sensitivity of the analysis of the four parameters. The first step is to connect to Python from R. Load the library and specify the Python environment you want to use.

library(reticulate)
# I will use my base Python 3 installation
use_python("C:/python/python36")

I already installed SALib (via pip) so I can immediately import the module.

# import numpy module
numpy = import('numpy')

# import some functions from SALib
saltelli = import('SALib.sample.saltelli')
sobol = import('SALib.analyze.sobol')

I’m importing two submodules from SALib: sample.saltelli and analyze.sobol. I’ll use saltelli to generate a collection of parameter values to use with rivr. I’ll then use sobol to analyze the results and generate the sensitivity analysis.

To use sample.saltelli, you have to provide the problem definition in a specific format. We have to be explicit with our types (e.g. integers vs. floats, vectors vs. lists) here, hence the explicit integer specification of num_vars and representing bounds as a list of lists.

problem = list(
  num_vars = 4L,
  names = list('bottom.width', 'bed.slope', 'side.slope', 'roughness'),
  bounds = list(
    list(40.0, 60.0),
    list(0.0005, 0.002),
    list(0.0, 2.0),
    list(0.003, 0.008)
  )
)

param.values = as.data.frame(saltelli$sample(problem, 1000L))
names(param.values) = problem$names

now that we’ve generated our parameter sets with saltelli.sample, it’s time to run the simulations. I first define a function to run a simulation and extract the time it takes for the flood crest to pass a point 30,000 ft downstream, and the crest height.

library(rivr)
#  values are in imperial units (ft, lbs, etc.)
Cm = 1.486
gravity = 32.2

# model resolution
numnodes = 301
xresolution = 250
tresolution = 2.5; 
times = seq(0, 30000, by = tresolution)

# baseflow and floodwave 
baseflow = 2500
floodwave = ifelse(times >= 9000, baseflow,
  baseflow + (7500/pi)*(1 - cos(pi*times/(60*75))))

ds.distance = 3.0e4
obs.location = ds.distance/xresolution + 1
obs.times = seq_along(floodwave)
run_model = function(...) {
  route_wave(..., Cm = Cm, g = gravity, 
    initial.condition = baseflow, boundary.condition = floodwave, 
    downstream.condition = rep(-1, length(floodwave)), 
    timestep = tresolution, spacestep = xresolution, numnodes = numnodes, 
    monitor.nodes = obs.location, monitor.times = obs.times,
    engine = "Dynamic", scheme = "MacCormack", boundary.type = "QQ") 
}

Now we’re ready to run the simulations:

outputs = vector("list", nrow(param.values))
results = vector("list", nrow(param.values))
# this will take a while...
for(i in seq_along(results)) {
  if (i %% 50 < 1)
	print(sprintf("loop %s", i))
  outputs[[i]] = as.data.frame(subset(
    with(param.values[i,], 
      run_model(So = bed.slope, n = roughness, 
        B = bottom.width, SS = side.slope)
    ),
	node == obs.location))
  peak.index = with(outputs[[i]], which.max(depth))
  results[[i]] = outputs[[i]][peak.index, c("time", "depth")]
}
crest.times = sapply(results, function(x) x$time)
crest.heights = sapply(results, function(x) x$depth)

Once that’s done, we use sobol.analyze to analyze the results.

sensitivity.height = sobol$analyze(problem,
  numpy$asarray(crest.heights, dtype = numpy$float))

sensitivity.time = sobol$analyze(problem,
  numpy$asarray(crest.times, dtype = numpy$float))

sobol.analyze returns a named list structure. Each list element holds a different set of sensitivity indices; these indices tell you something about how the model output variance relates to the model inputs. The first-order indices (list element S1) measure the contribution to the output variance by a single model input alone. Second-order indices (list element S2) measur the contribution to the output variance caused by the interaction of two model inputs. Finally, the total-order indices (list element ST) measures the contribution to the output variance caused by a model input, including both its first-order effects (the input varying alone) and all higher-order interactions. Confidences intervals for these indices are also returned (list elements S1_conf, S2_conf, and ST_conf).

With a little bit of help from dplyr and ggplot2, we can visualize these indices:

library(dplyr)
sens.ind = bind_rows(
  "Crest Height" = bind_rows(
    first = data.frame(
      parameter = unlist(problem$names),
      value = sensitivity.height$S1,
      confint = sensitivity.height$S1_conf
    ),
    total = data.frame(
      parameter = unlist(problem$names),
      value = sensitivity.height$ST,
      confint = sensitivity.height$ST_conf
    ),
	.id = "index"
  ),
  "Travel Time" = bind_rows(
    first = data.frame(
      parameter = unlist(problem$names),
      value = sensitivity.time$S1,
      confint = sensitivity.time$S1_conf
    ),
    total = data.frame(
      parameter = unlist(problem$names),
      value = sensitivity.time$ST,
      confint = sensitivity.time$ST_conf
    ),
	.id = "index"
  ),
  .id = "Output"
)

library(ggplot2)
sens.ind %>% mutate(
  conf.upper = value + confint,
  conf.lower = if_else(confint > value, 0, value - confint),
  parameter = gsub(".", " ", parameter, fixed = TRUE)
) %>%
  ggplot() + aes(x = parameter, y = value, fill = index,
  ymin = conf.lower, ymax = conf.upper) + 
  geom_col(position = "dodge", width = 0.5) + 
  geom_errorbar(position = position_dodge(0.5), width = 0.25) +
  facet_wrap(~Output, scales = "free_y", ncol = 1) +
  theme(legend.position = "top") + 
  scale_fill_brewer(NULL) + theme_bw() +
  ylab("Sensitivity") + xlab(NULL)

Note that the sensitivity indices cannot be negative, so a confidence interval that results in a negative value for the lower bound of a sensitivity index can be interpreted as equal to zero. The results are mostly as expected: bed slope and roughness have the most significant effect both on crest height and travel time. Side slope appears to have a low first-order effect on crest height, but becomes more significant when higher-order effects are taken into account (we’ll get to that in a minute). Bottom width has a significant effect on crest height, which is to be expected because a wider channel will by definition result in lower depths for the same flow. More surprising is that bottom width does not have a significant first-order effect on travel time, and potentially no significant higher-order effects. This is surprising, considering that the velocity of the flood wave is directly proportional to the area of flow under one-dimensional assumptions. Furthermore, there are significant higher-order effects from side slope, which similarly impacts flow velocities (since a larger side-slope results in larger flow areas). The uncertainty in the effect of bottom width is probably due to the fact that we explored a very large parameter range (50–200 feet). It might turn out that more iterations would clarify the effect of bottom width.

The “total” effects returned by sobol.analyze combine both the first-order effect of the parameter and all higher-order effects. We can also explore second-order effects of the parameters directly. Visualizing this is a little harder, but I took a shot at it by mapping point sizes to the indices and confidence intervals.

library(tidyr)
second.ind = list(
  "Crest Height" = setNames(
    as_data_frame(sensitivity.height$S2), 
	unlist(problem$names)
  ) %>% 
    mutate(var1 = unlist(problem$names)) %>% 
    gather(var2, value, -var1),
  "Travel Time" = setNames(
    as_data_frame(sensitivity.time$S2), 
	unlist(problem$names)
  ) %>% 
    mutate(var2 = unlist(problem$names)) %>% 
    gather(var1, value, -var2)
) %>%
  bind_rows(.id = "Output") %>%
  mutate(
    var1 = factor(var1, levels = unlist(problem$names), 
	  labels = gsub(".", " ", unlist(problem$names), fixed = TRUE)),
    var2 = factor(var2, levels = unlist(problem$names), 
	  labels = gsub(".", " ", unlist(problem$names), fixed = TRUE))
  ) 
second.conf = list(
  "Crest Height" = setNames(
    as_data_frame(sensitivity.height$S2_conf), 
	unlist(problem$names)
  ) %>% 
    mutate(var1 = unlist(problem$names)) %>% 
    gather(var2, confint, -var1),
  "Travel Time" = setNames(
    as_data_frame(sensitivity.time$S2_conf), 
	unlist(problem$names)
  ) %>% 
    mutate(var2 = unlist(problem$names)) %>% 
    gather(var1, confint, -var2)
) %>%
  bind_rows(.id = "Output") %>%
  mutate(
    var1 = factor(var1, levels = unlist(problem$names), 
	  labels = gsub(".", " ", unlist(problem$names), fixed = TRUE)),
    var2 = factor(var2, levels = unlist(problem$names), 
	  labels = gsub(".", " ", unlist(problem$names), fixed = TRUE))
  ) 

second.ind %>% left_join(second.conf) %>%
  mutate(
    lower = if_else(confint > value, 0, value - confint),
	upper = value + confint
  ) %>%
  ggplot() + aes(x = var1, y = var2) + 
  geom_point(aes(size = upper, fill = Output), shape = 22) + 
  geom_point(aes(size = value), shape = 22) + 
  geom_point(aes(size = lower), shape = 22) + 
  theme_bw() + xlab(NULL) + ylab(NULL) + 
  scale_size_area("Effect", max_size = 20) +
  scale_fill_brewer(NULL, palette = "Accent") +
  coord_fixed()

The above plot shows the second-order effects, with the inner/outer squares indicating the confidence intervals. Truthfully, I’m not that excited about this visualization, but it does illustrate relative effect sizes and uncertainty. Except for bed slope:roughness effect on travel time, the second-order effects of all interactions are potentially negligible (since the confidence interval includes zero). Also, while this doesn’t show up in the plot (because of the way I defined the size scale), some of my second-order indices are negative. This usually indicates a problem with sample size, i.e. the parameter space wasn’t sufficiently explored to get a good estimate of the effect.

Okay, that’s all I have for SALib and reticulate today. Now to rewrite arcpyrextra to use reticulate…