---
title: "Workflow Eisenstadt 2005"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Workflow Eisenstadt 2005}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE, eval=TRUE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
is_ghactions <- tolower(Sys.getenv("GITHUB_ACTIONS")) == "true" ||
  tolower(Sys.getenv("CI")) %in% c("true", "1", "yes")

# base.h5 ships in inst/extdata/models/eisenstadt-2005/. The vignette only
# renders if that file is present; the engine .exe is fetched from the
# kwb.raindrop.binaries Release on demand and only runs on Windows.
path_base <- system.file("extdata/models/eisenstadt-2005/base.h5", package = "kwb.raindrop")
data_available <- nzchar(path_base) && file.exists(path_base)
is_windows <- Sys.info()[["sysname"]] == "Windows"
```

### Input data

The HDF5 model template (`base.h5`) ships with the package under
`inst/extdata/models/eisenstadt-2005/` and is produced with the
Tandler "Regenwasserbewirtschaftung" calculation engine; the engine
itself is downloaded from the `KWB-R/kwb.raindrop.binaries` GitHub
Release on demand via `kwb.raindrop::download_engine()`.

### Define Paths and Scenarios

```{r preparation, eval = data_available}
library(kwb.raindrop)

path_list <- list(
  modelname = "Eisenstadt_2005",
  root_path = file.path(tempdir(), "raindrop_eisenstadt_2005"),
  dir_input  = "<root_path>/models/<modelname>/input",
  dir_output = "<root_path>/models/<modelname>/output",
  dir_target_output = "<dir_output>/<dir_target>",
  file_errors_hdf5 = "Fehlerprotokoll.h5",
  file_results_hdf5_element = "Mulde_Rigole.h5",
  file_results_hdf5_flaeche = "Dach.h5",
  file_results_hdf5_verschaltungen = "<dir_target>_Verschaltungen.h5",
  file_results_txt = "Mulde_Rigole_RAINDROP.txt",
  file_results_txt_multilayer = "Mulde_Rigole_RAINDROP_multi_layer.txt",
  file_target = "<dir_target>.h5",
  path_base = system.file("extdata/models/eisenstadt-2005/base.h5", package = "kwb.raindrop"),
  path_exe  = if (is_windows && !is_ghactions) kwb.raindrop::download_engine() else NA_character_,
  path_errors_hdf5 = "<dir_target_output>/<file_errors_hdf5>",
  path_results_hdf5_element = "<dir_target_output>/<file_results_hdf5_element>",
  path_results_hdf5_flaeche = "<dir_target_output>/<file_results_hdf5_flaeche>",
  path_results_hdf5_verschaltungen = "<dir_target_output>/<file_results_hdf5_verschaltungen>",
  path_results_txt = "<dir_target_output>/<file_results_txt>",
  path_results_txt_multilayer = "<dir_target_output>/<file_results_txt_multilayer>",
  path_target_input = "<dir_input>/<file_target>"
)


parameters <- tibble::tibble(
  para_nama_short = c(
                    "connected_area", 
                    "mulde_area", 
                    "mulde_height",
                    "filter_hydraulicconductivity",
                   # "filter_height",
                    "storage_height"#,
                   # "bottom_hydraulicconductivity"
                    ),
  para_name_long = c(
                   "/Massnahmenelemente/Dach/Allgemein/Flaeche",
                   "/Massnahmenelemente/Mulde_Rigole/Allgemein/Flaeche",
                   "/Massnahmenelemente/Mulde_Rigole/Eigenschaften_Oberflaeche/Ueberlaufhoehe",
                   "Bodenarten/Bodenfilter/Ks_HydraulicConductivity",
                   #"/Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Schichtdicken",
                   "/Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Schichtdicken"#,
                  # "/Massnahmenelemente/Mulde_Rigole/Allgemein/Endversickerungsrate"
                   ),
  index = c(1L,
            1L,
            1L,
            1L,
            2)
)

DT::datatable(parameters, 
              filter = "top",
              options = list(pageLength = 25,
                             autoWidth = TRUE))


connected_area <- 1000
mulde_area <- c(25, 50, 75, 100, 125, 150, 175, 200)
mulde_height <- c(100, 200, 300) 
filter_hydraulicconductivity <-  c(36, 180, 360)
filter_height <- 300
storage_height <- c(100, 500, 1000)
rain_factor <- 1
bottom_hydraulicconductivity <- 12 #c(1,5,10,20,45,90,180,270,360,1860,3600)


# Alle Kombinationen erzeugen
param_grid_all_combinations <- expand.grid(
  connected_area = connected_area,
  mulde_area = mulde_area,
  mulde_height = mulde_height,
  filter_hydraulicconductivity = filter_hydraulicconductivity,
  filter_height = filter_height,
  storage_height = storage_height,
  bottom_hydraulicconductivity = bottom_hydraulicconductivity,
  rain_factor = rain_factor
)

param_grid_all_combinations <- param_grid_all_combinations %>% 
  dplyr::bind_cols(tibble::tibble(scenario_name = sprintf("s%05d", 
                                                          seq_len(nrow(param_grid_all_combinations)))))

ref_scenario <- param_grid_all_combinations %>% 
       dplyr::filter(connected_area == min(unique(param_grid_all_combinations$connected_area)),
                     mulde_area == min(unique(param_grid_all_combinations$mulde_area)), 
                     filter_height == min(filter_height),
                     filter_hydraulicconductivity == min(param_grid_all_combinations$filter_hydraulicconductivity), 
                     bottom_hydraulicconductivity == min(unique(param_grid_all_combinations$bottom_hydraulicconductivity)),
                     mulde_height == min(param_grid_all_combinations$mulde_height),
                     storage_height == min(param_grid_all_combinations$storage_height)) %>% 
       dplyr::pull(scenario_name)

stopifnot(length(ref_scenario)==1)

scenarios_with_single_parameter_variation <- kwb.raindrop::find_single_param_variations(
  data = param_grid_all_combinations,
  ref_scenario = ref_scenario
  ) %>% 
  dplyr::pull(scenario_name) %>% unique()

# param_grid <- param_grid_all_combinations  %>% 
#   dplyr::filter(scenario_name %in% scenarios_with_single_parameter_variation)
param_grid <- param_grid_all_combinations

DT::datatable(param_grid,
              filter = "top",
              options = list(pageLength = 25,
                             autoWidth = TRUE))

htmlwidgets::saveWidget(DT::datatable(parameters,
                                      filter = "top",
                                      options = list(pageLength = 25, 
                                                     autoWidth = TRUE)), "parameters.html")
htmlwidgets::saveWidget(DT::datatable(param_grid, 
                                      filter = "top",
                                      options = list(pageLength = 25, 
                                                     autoWidth = TRUE)), "param_grid.html")


psi_s_mm <- function(kf_mmh) (3.237 * (kf_mmh/25.4)^(-0.328)) * 25.4

paths <- kwb.utils::resolve(path_list)

```

### Run Model


```{r run_model, eval = data_available && is_windows && !is_ghactions}
run_one <- function(i,
                    timestep_hours,
                    debug = FALSE,
                    ...) {

  param_grid_tmp <- param_grid[i, ]

  paths <- kwb.utils::resolve(path_list,
                              dir_target = param_grid_tmp$scenario_name)

      fs::dir_create(paths$dir_input, recurse = TRUE)
      fs::dir_create(paths$dir_output, recurse = TRUE)
      fs::dir_create(paths$dir_target_output, recurse = TRUE)

      fs::file_copy(path = paths$path_base,
                    new_path = paths$path_target_input,
                    overwrite = TRUE)

      h5 <- hdf5r::H5File$new(paths$path_target_input, mode = "a")

      new_path <- stringr::str_c(normalizePath(fs::path_abs(paths$dir_target_output)),
                                 "\\")

      vals <- kwb.raindrop::h5_read_values(h5)

      vals$`//Berechnungsparameter/Ergebnispfad` <- new_path
      vals$`//Berechnungsparameter/Zeitschritt_Infiltration` <- timestep_hours
      vals$`//Berechnungsparameter/Zeitschritt_ET` <- timestep_hours
      vals$`//Berechnungsparameter/Zeitschritt_Verschaltungen` <- timestep_hours
      vals$`//Berechnungsparameter/R-Plots` <- 0
      vals$`//Berechnungsparameter/Ausgabemodus` <- "Optimierung"
      vals$`//Berechnungsparameter/Evapotranspiration_aktiv` <- 1
      
      vals$`//Massnahmenelemente/Dach/Berechnungsparameter/Evapotranspiration_aktiv` <- 1
      vals$`//Massnahmenelemente/Dach/Allgemein/Flaeche` <- param_grid_tmp$connected_area
      
      vals$`//Massnahmenelemente/Mulde_Rigole/Berechnungsparameter/Evapotranspiration_aktiv` <- 1
      vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Regen-Skalierungsfaktor` <- 1
      vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Flaeche` <- param_grid_tmp$mulde_area
      vals$`//Massnahmenelemente/Mulde_Rigole/Eigenschaften_Oberflaeche/Ueberlaufhoehe` <- param_grid_tmp$mulde_height
      vals$`//Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Startwerte_theta_ActualSoilMoisture` <- c(0.3, 0)
      vals$`//Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Schichtdicken` <- c(param_grid_tmp$filter_height,
                                                                                 param_grid_tmp$storage_height)
      vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Endversickerungsrate` <- param_grid_tmp$bottom_hydraulicconductivity
      
      vals$`//Bodenarten/Bodenfilter/Ks_HydraulicConductivity` <- param_grid_tmp$filter_hydraulicconductivity
      vals$`//Bodenarten/Bodenfilter/Psi_Saugspannung_CapillarySuction` <- psi_s_mm(param_grid_tmp$filter_hydraulicconductivity)

# Timeseries (2×N) als tibble?
if (is.data.frame(vals[["//Kurven/Regen"]])) {
  vals[["//Kurven/Regen"]]$value <- vals[["//Kurven/Regen"]]$value * param_grid_tmp$rain_factor
}

      kwb.raindrop::h5_write_values(h5, vals, resize = TRUE,
                                   scalar_strategy = "error",
                                   verbose = FALSE)
      h5$close_all()

      kwb.raindrop::run_model(path_exe = paths$path_exe,
                              path_input = paths$path_target_input,
                              debug = debug)

      invisible(NULL)
}

n_cores <- parallel::detectCores()

system.time(expr = {
kwb.raindrop::run_scenarios(indices = seq_len(nrow(param_grid)),
                            run_one_scenario = run_one,
                            timestep_hours = 0.1,
                            debug = FALSE,
                            parallel = TRUE,
                            workers = n_cores,
                            show_progress = TRUE
                            )
}
)

### Read results for first run
if(FALSE) {
paths <- kwb.utils::resolve(path_list, 
                            dir_target = sprintf("s%05d", i = 1))

#simulation_names <- basename(fs::dir_ls(paths$dir_output))
simulation_names <- param_grid$scenario_name

debug <- TRUE
errors_df <- kwb.raindrop::read_raindrop_errors(simulation_names, path_list)
x <- tidyr::unnest(errors_df, errors) 
x$Fehlerbeschreibung
}
```


### Analyse Results

```{r analyse_results, eval = data_available && is_windows && !is_ghactions}
system.time(
simulation_results <- kwb.raindrop::get_simulation_results_optim_parallel(
  paths = paths, 
  path_list = path_list,
  simulation_names = param_grid$scenario_name,
  debug = FALSE)
)

system.time(
simulation_results_optimisation <- kwb.raindrop::add_overflow_events_and_waterbalance(
  simulation_results = simulation_results,
  event_separation_hours = 4,
  canonical_variables = kwb.raindrop::default_canonical_wb_variables()
  )
)

simulation_results_optimisation <- param_grid %>% 
  dplyr::left_join(simulation_results_optimisation,
                   by = c("scenario_name" = "s_name")) %>% 
  dplyr::relocate(scenario_name, .before = connected_area)

readr::write_csv(simulation_results_optimisation, 
                 file = sprintf("simulation_results_optimisation_%s.csv",
                                paths$modelname)
                 )

htmlwidgets::saveWidget(DT::datatable(simulation_results_optimisation,
                                      filter = "top",
                                      options = list(pageLength = 25, 
                                                     autoWidth = TRUE)), 
                        file = sprintf("simulation_results_optimisation_%s.html",
                                       paths$modelname),
                        title = sprintf("RAINDROP - '%s': Solution Space",
                                        paths$modelname)
)
### Plot results

params <- c(
  #"connected_area",
  "mulde_area",
  "mulde_height",
  "filter_hydraulicconductivity",
  #"filter_height",
  "storage_height"
  #"bottom_hydraulicconductivity",
  #"rain_factor"
)

lang <- "de"
max_n_overflows <- 1

pdff <- sprintf("simulation_results_optimisation_%s_main-effects.pdf",
                paths$modelname)

gg <- kwb.raindrop::plot_main_effects(
  df = simulation_results_optimisation,
  y = "n_overflows",
  params = params,
  lang = lang,
  ylim_lower = 0
)

grDevices::pdf(pdff, width = 9, height = 3, onefile = TRUE)
print(gg)
dev.off()
#kwb.utils::finishAndShowPdf(pdff)


pdff <- sprintf("simulation_results_optimisation_%s_design-space_mulde-area_vs_mulde-height_single.pdf",
                paths$modelname)

# kwb.utils::preparePdf(pdfFile = pdff)
# on.exit(dev.off(), add = TRUE)
grDevices::pdf(pdff, width = 9, height = 4, onefile = TRUE)

  p <- kwb.raindrop::plot_valid_design_space(
    param_grid = param_grid,
    sim_results = simulation_results_optimisation,
    x = "mulde_area",
    y = "mulde_height",
    valid_max = max_n_overflows,
    jitter = TRUE,
    alpha_mode = "duplicates",
    alpha_min = 0.25,
    alpha_max = 1,
    drop_overflow_gt_valid_max = FALSE,
    keep_param_grid_limits = TRUE,
    lang = lang,
    subtitle = ""
  )

  suppressWarnings(print(p))

dev.off()


# --- 2) Interaktiv als HTML: nach dem PDF
plotly_gg <- plotly::ggplotly(gg)

htmlwidgets::saveWidget(
  widget = plotly_gg,
  file = sprintf("simulation_results_optimisation_%s_main-effects.html",
                 paths$modelname),
  selfcontained = TRUE,
  title = sprintf("RAINDROP - '%s': Main Effects",
                  paths$modelname)
)



pdff <- sprintf("simulation_results_optimisation_%s_design-space_mulde-area_vs_parameters.pdf",
                paths$modelname)
kwb.utils::preparePdf(pdfFile = pdff)

for (y in c("mulde_height", "filter_hydraulicconductivity", "storage_height")) {

  p <- kwb.raindrop::plot_valid_design_space(
    param_grid = param_grid,
    sim_results = simulation_results_optimisation,
    x = "mulde_area",
    y = y,
    valid_max = max_n_overflows,
    jitter = TRUE,
    alpha_mode = "duplicates",
    alpha_min = 0.25,
    alpha_max = 1,
    drop_overflow_gt_valid_max = TRUE,
    keep_param_grid_limits = TRUE
  )

  # interaktiv als HTML
  plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
  htmlwidgets::saveWidget(
    widget = plotly_p,
    file = sprintf("simulation_results_optimisation_%s_design-space_mulde-area_vs_%s.html", 
                   paths$modelname,
                   y),
    selfcontained = TRUE,
    title = sprintf("'%s' - Design Space: mulde_area vs. %s", 
                    paths$modelname, 
                    y)
  )

  # statisch ins PDF (WICHTIG!)
  suppressWarnings(print(p))
}
dev.off()
#kwb.utils::finishAndShowPdf(pdff)

pdff <- sprintf("simulation_results_optimisation_%s_water-balance.pdf",
                paths$modelname)
kwb.utils::preparePdf(pdfFile = pdff)

p <- kwb.raindrop::plot_wb_tradeoff_overflows(
  simulation_results_optimisation = simulation_results_optimisation,
  x = max_n_overflows,
  param_grid = param_grid,
  filter_n_gtx = TRUE,
  use_jitter = TRUE
  )

  # interaktiv als HTML
  plotly_p <- suppressWarnings(plotly::ggplotly(p, tooltip = "text"))
  htmlwidgets::saveWidget(
    widget = plotly_p,
    file = sprintf("simulation_results_optimisation_%s_water-balance.html",
                   paths$modelname),
    selfcontained = TRUE,
    title = sprintf("'%s' - Water balance vs overflows", 
                    paths$modelname)
  )

  # statisch ins PDF (WICHTIG!)
  suppressWarnings(print(p))  
dev.off()
#kwb.utils::finishAndShowPdf(pdff)

```