Automating X number of ensemble runs with unique parameter values

[Scoobys-keeper]

Hi all!

I have scoured through the documentation to find an example, but my apologies if I have missed the answer to this.

For my masters project, I am modeling global temperature change in 2100 relative to a baseline from 1960-1990 in response to 4 key perturbed parameters (Climate sensitivity, ocean heat diffusivity, Volcanic and aerosol forcing). I used a latin hypercube in r to generate values for these parameters within these ranges, and my issue now is how to automate the running of these different values and saving their outputs. Here is a photo of my current code trying to create a loop, but I feel I'm missing a major component here.

Ideally I'll save and plot all the ensembles (in this trial case 35 but hoping to increase this to around 50) and then conduct a sensitivity analysis with historical observations, ruling out those that don't meet my implausibility threshold. My ideal output would be a file with all the ensembles in csv format and a plot showing how they compare to HADCrut historical observational data.

Thank you so much for your time and help!!

[kdorheim]

Hi @Scoobys-keeper exciting to see Hector being used this way! I think you're really close, but the for-loop might need some tweaks. Do you mind sharing your R script?

[Scoobys-keeper]

@kdorheim Wow- Very honored to have your assistance! Thank you very much for your time and bandwidth here.

Could I also highlight an issue I have when trying to use RCPs, being the error around 'could not parse var: rho - unknown variable name'. Ill paste the whole script here and the error I receive when trying to use any of the RCPs (doesn't happen when I use the SSPs in my hector core, however)

The portion of interest is down towards the bottom, under the final set of hashmarks titled #trial attempt to change...

ini_file <- system.file("input/hector_rcp45.ini", package = "hector")
core <- newcore(ini_file)
core

ecs <- fetchvars(core, NA, ECS())
OceanDiffusivity <- fetchvars(core, NA, DIFFUSIVITY())
AerosolForcing <- fetchvars(core, NA, AERO_SCALE())
VolcanicForcing <- fetchvars(core, NA, VOLCANIC_SCALE())
beta <- fetchvars(core, NA, BETA())
PreindustrialC02 <- fetchvars(core, NA, PREINDUSTRIAL_CO2())
Hetero_respiration_tempSensi <- fetchvars(core, NA, Q10_RH())

#################################
#calibrated Values from doreheim et al
#showing difference from hector default to calibrated values 
#trying to troubleshoot as of 07.04.24- why this only returning one line on graph 
setvar(core, NA, DIFFUSIVITY(), 1.16, "cm2/s")
fetchvars(core, NA, DIFFUSIVITY())
core
reset(core)
core
run(core)
results_OD_1.16 <- fetchvars(core, 1850:2024)
results[["DIFFUSIVITY"]] <- 2.38
results_OD_1.16[["DIFFUSIVITY"]] <- 1.16
compare_results <- rbind(results, results_OD_1.16)

ggplot(compare_results) +
  aes(x = year, y = value, color = factor((DIFFUSIVITY()))) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y") +
  guides(color = guide_legend(title = expression(DIFFUSIVITY())))

###########################################
#lhs set up

install.packages("lhs")
library(lhs)
set.seed(2455)
#35 ensembles and 4 parameters as trial
lhs_design = maximinLHS(35, 4)

colnames(lhs_design) <- c("ECS", "OceanDiff","AerosolF","VolcanicF")
ranges <- 
  list(lower = c(1, 1, 0, 0),
       upper = c(4.5, 5, 2.44, 3.66))

cube2box <- function(X, lower, upper, check = TRUE, eps = 1e-6) {
  is.axial <- inherits(X, "axial")
  stopifnot(is.matrix(X))
  p <- ncol(X)
  stopifnot(length(lower) == p, length(upper) == p)
  X <- t(X)
  if (check) stopifnot(0 - eps < X, X < 1 + eps)
  # This bit:
  robj <- t(X * (upper - lower) + lower)
  if (is.axial)
    class(robj) <- c("matrix", "axial")
  robj
}

lhs_des <- maximinLHS(35, 4)
ens_params <- cube2box(lhs_des, ranges$lower, ranges$upper)
ens_params
write.csv(ens_params, "ens_params.csv")

##############################################
#trial attempt to change Ocean diffusivity and run.
#Keep in mind calibrated parameters from Dorheim et al., 2023 to match observation data.

parameter_values <- read.csv("ens_params.csv")  

for (i in 1:nrow(parameter_values)) {
  # Extract parameter values for the current ensemble run
  Climate.Sensitivity <- parameter_values[i, "Climate.Sensitivity"]
  Ocean.Heat.diffusivity <- parameter_values[i, "Ocean.Heat.diffusivity"]
  Aerosol.forcing <- parameter_values[i, "Aerosol.forcing"]
  Volcanic.forcing <- parameter_values[i, "Volcanic.forcing"]
}

setvar(core, NA, DIFFUSIVITY(), Ocean.Heat.diffusivity, "cm2/s")
run(core)
OD_attempt <- fetchvars(core, 1850:2024)

OD_attempt

ggplot(OD_attempt) +
  aes(x = year, y = value) +
  geom_line() +
  facet_wrap(~variable, scales = "free_y")

The error I mentioned from the initial try with rcp 4.5

Error: While parsing hector input file /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/hector/input/hector_rcp45.ini: msg:     Could not parse var: rho - Unknown variable name while parsing CH3Br_halocarbon: rho
func:   setData
file:   halocarbon_component.cpp
ffile:  halocarbon_component.cpp

line:   155

Cant thank you enough for your time!!

Willy Phillips

[kdorheim]

Hi @Scoobys-keeper, awesome thanks, that is really helpful!

Re the RCP issue - As of V3 the RCPs are no longer shipped with the Hector package, the error message you are seeing there is related to the fact that the old RCP ini files are out of date (missing parameters). Would you be able to use the SSP scenarios? Or you need to use the RCPs?

Re the sensitivity runs - as it was originally written the for loop was running over your data frame of the different parameters and passing that information to the Hector core, but the lines of code actually running the Hector instance and then fetching/formatting the results also need to be in the for loop as well. See below your code with some modifications/comments that should get you the kind of behavior you are interested in.

# because of the hector version we will be running the ssps instead of the RCPS
ini_file <- system.file("input/hector_ssp245.ini", package = "hector")
core <- newcore(ini_file)
core

ecs <- fetchvars(core, NA, ECS())
OceanDiffusivity <- fetchvars(core, NA, DIFFUSIVITY())
AerosolForcing <- fetchvars(core, NA, AERO_SCALE())
VolcanicForcing <- fetchvars(core, NA, VOLCANIC_SCALE())
beta <- fetchvars(core, NA, BETA())
PreindustrialC02 <- fetchvars(core, NA, PREINDUSTRIAL_CO2())
Hetero_respiration_tempSensi <- fetchvars(core, NA, Q10_RH())

#################################
#calibrated Values from doreheim et al
#showing difference from hector default to calibrated values
#trying to troubleshoot as of 07.04.24- why this only returning one line on graph
setvar(core, NA, DIFFUSIVITY(), 1.16, "cm2/s")
fetchvars(core, NA, DIFFUSIVITY())
core
reset(core)
core
run(core)
results_OD_1.16 <- fetchvars(core, 1850:2024)
results[["DIFFUSIVITY"]] <- 2.38
results_OD_1.16[["DIFFUSIVITY"]] <- 1.16
compare_results <- rbind(results, results_OD_1.16)

ggplot(compare_results) +
    aes(x = year, y = value, color = factor((DIFFUSIVITY()))) +
    geom_line() +
    facet_wrap(~variable, scales = "free_y") +
    guides(color = guide_legend(title = expression(DIFFUSIVITY())))

###########################################
#lhs set up

install.packages("lhs")
library(lhs)
set.seed(2455)
#35 ensembles and 4 parameters as trial
lhs_design = maximinLHS(35, 4)

colnames(lhs_design) <- c("ECS", "OceanDiff","AerosolF","VolcanicF")
ranges <-
    list(lower = c(1, 1, 0, 0),
         upper = c(4.5, 5, 2.44, 3.66))

cube2box <- function(X, lower, upper, check = TRUE, eps = 1e-6) {
    is.axial <- inherits(X, "axial")
    stopifnot(is.matrix(X))
    p <- ncol(X)
    stopifnot(length(lower) == p, length(upper) == p)
    X <- t(X)
    if (check) stopifnot(0 - eps < X, X < 1 + eps)
    # This bit:
    robj <- t(X * (upper - lower) + lower)
    if (is.axial)
        class(robj) <- c("matrix", "axial")
    robj
}

lhs_des <- maximinLHS(35, 4)
ens_params <- cube2box(lhs_des, ranges$lower, ranges$upper)
ens_params
write.csv(ens_params, "ens_params.csv")

##############################################
#trial attempt to change Ocean diffusivity and run.
#Keep in mind calibrated parameters from Dorheim et al., 2023 to match observation data.

parameter_values <- read.csv("ens_params.csv")

# Adding the parameter names names to the data frame so the subsetting in the for
# loop works.
names(parameter_values) <- c("Climate.Sensitivity", "Ocean.Heat.diffusivity",
                             "Aerosol.forcing", "Volcanic.forcing")

# Make an empty data frame to store the results in!
output <- data.frame()

# This for loop is running over each of the rows in the parameter_values data frame
for (i in 1:nrow(parameter_values)) {
    # Adding some print statements can be really helpful when debugging a for loop!
    print(i)

    # Extract parameter values for the current ensemble run
    Climate.Sensitivity <- parameter_values[i, "Climate.Sensitivity"]
    Ocean.Heat.diffusivity <- parameter_values[i, "Ocean.Heat.diffusivity"]
    Aerosol.forcing <- parameter_values[i, "Aerosol.forcing"]
    Volcanic.forcing <- parameter_values[i, "Volcanic.forcing"]

    # Setting the new diffusivity parameter value, running the model, and
    # fetching the results is going to need to happen in the for loop.
    print(cat("Ocean.Heat.diffusivity ", Ocean.Heat.diffusivity))
    setvar(core = core, dates = NA, var = DIFFUSIVITY(), values = Ocean.Heat.diffusivity, unit = "cm2/s")
    # TODO when you are ready add more setvar calls for the other parameters you'd like to vary.

    # reset the core to make sure the core is ready to run and uses
    # the new parameter value, this is an important step when using a for loop
    # to rerun the same hector core over and over again.
    reset(core)

    run(core)

    # Get the variables!
    OD_attempt <- fetchvars(core, 1850:2024)

    # Add information about which parameters were used to for the Hector run
    OD_attempt[, "scenario"] <- i

    # Since every time the for loop is running the OD_attempt data frame
    # is going to be overwritten so we will want to concatenate the new
    # results to previous run results.
    output <- rbind(output, OD_attempt)

}

library(ggplot2)

ggplot(output) +
    aes(x = year, y = value, color = factor(scenario)) +
    geom_line() +
    facet_wrap(~variable, scales = "free_y")

[Scoobys-keeper]

@kdorheim Thank you so very much for your assistance. I’ll give this a bash and report back with my results! Thank you again. And it is just fine to use the SSPs- just wanted to ask! Thanks!!

Willy

[bpbond]

FYI: you might be interested in @jk-brown 's "matilda" R package that provides tools for conducting sensitivity and probabilistic analyses using Hector.

[Scoobys-keeper]

@bpbond Very well received. will be working through that today! Thank you for your help.

[kdorheim]

@Scoobys-keeper, just following up. How are things? Is it okay to close this issue?

[Scoobys-keeper]

Hi @kdorheim !

Thanks for following up- yes, all is working as expected and I Appreciate All your help!