# R packages
library(reticulate)
library(hydroGOF) # for KGE and KGE'
library(dplyr)
library(tidyr)
library(ggplot2)
library(cowplot)

# Python packages
de <- reticulate::import("de") # diaf-eff Python package
hydroErr <- reticulate::import("HydroErr") # HydroEff Python package

# Functions
source("script/utils.R")
source("script/run_eval.R")
source("script/cblind_bw_palette.R")

# Import synthetic time series of flood events
## syn1: one event
## syn2: two consecutive events
source("script/Synthetic time series/syn_flood_event.R")

# Define omega bounds
omega_log_min <- -0.36
omega_log_max <- 0.36

# Calculate omega bounds (min/max)
omega_min <- 10^(omega_log_min)
omega_max <- 10^(omega_log_max)

# Create min and max transformed time series
model_min <- c(syn2$discharge * omega_min)
model_max <- c(syn2$discharge * omega_max)

# Generate plot
syn2 |> 
  mutate(model_min, model_max) |> 
  ggplot(aes(t, discharge)) +
  geom_ribbon(aes(ymin = model_min, 
                  ymax = model_max, 
                  fill = "Set of transformed time series")) +
  geom_line(aes(color = "Reference time series")) +
  scale_color_manual(name = "", values = "black") +
  scale_fill_manual(name = "", values = "lightgrey") +
  xlab("Time [T]") +
  ylab(expression(paste("Discharge [L"^3~T^-1, "]"))) +
  coord_cartesian(xlim = c(0, 40)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  theme_bw() +
  theme(legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        strip.text = element_text(size = 10),
        legend.text = element_text(size = 10))

# Define omega values
omega_1 <- 0.75
omega_2 <- 1.2

# Create two transformed time series
model_bad_bad <- c(syn1$discharge * omega_1, syn1$discharge * omega_2)
model_bad_good <- c(syn1$discharge * omega_1, syn1$discharge)

# Generate plot
syn2 |>
  mutate(model_bad_bad, model_bad_good) |> 
  pivot_longer(c("model_bad_bad", "model_bad_good")) |> 
  ggplot() +
  geom_point(aes(x = t, y = discharge, shape = "Reference time series")) +
  geom_line(aes(x = t, y = value, color = "Transformed time series")) +
  facet_wrap(~name, 
             labeller = as_labeller(c("model_bad_bad" = '"Bad-Bad" model (BB model)', 
                                      "model_bad_good" = '"Bad-Good" model (BG model)'))) +
  coord_cartesian(xlim = c(0, 40)) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  scale_shape_manual(name = "", values = 4) +
  scale_color_manual(name = "", values = "#696969") +
  xlab("Time [T]") +
  ylab(expression(paste("Discharge [L"^3~T^-1, "]"))) +
  theme_bw() +
  theme(legend.position = "bottom",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.title = element_text(size = 12),
        axis.text = element_text(size = 10),
        strip.text = element_text(size = 10),
        legend.text = element_text(size = 10))

source("script/Synthetic time series/compare_score_param.R")

# Define omega values
omega_1 <- 0.75
omega_2 <- 1.2

# Create two transformed time series
model_bad_bad <- c(syn1$discharge * omega_1, syn1$discharge * omega_2)
model_bad_good <- c(syn1$discharge * omega_1, syn1$discharge)

# Generate plot
compare_score_param(syn2$discharge, model_bad_bad, model_bad_good)

source("script/Synthetic time series/plot_model_distribution.R")

# Vector of criteria to study
criteria <- c("kge", "kge_m", "kge_m2", "kgenp", "de", "lme", "nse", "dr", "lce")

# Generate plot
plot_model_distribution(data = syn1,
                        perf_c = criteria,
                        ncol = 3,
                        omega_log_min = omega_log_min,
                        omega_log_max = omega_log_max,
                        step = 0.002)

source("script/Synthetic time series/plot_omega_vs_omega2.R")

# Vector of criteria to study
criteria <- c("kge", "kge_m", "kge_m2", "kgenp", "de", "lme", "nse", "dr", "lce")

# Generate plot
plot_omega_vs_omega2(data = syn1,
                     perf_c = criteria,
                     omega_log_min = omega_log_min,
                     omega_log_max = omega_log_max,
                     step = 0.002)