library(tidyverse)
library(leaflet)
library(sf)
library(plotly)
# library(anomalize)
library(slider)
library(anomaly)

# read custom functions
source('../scripts/test_get_tsfeatures.R')

locs <- read_csv2("../data/raw/locatie_mapping.csv")
sensor_lims <- read_csv2("../data/raw/meet_range_sensoren.csv") %>% 
  filter(path_length %in% c(1, NA))
par_ranges <- readxl::read_excel('../data/raw/validation_master_v1.xlsx', sheet = 'realistic_range')
df_dt <-
  arrow::read_parquet("../data/modified/sensor_data_tot.parquet", as_data_frame = FALSE) %>%
  filter(project == "Project RWZI Dinther") %>%
  collect() %>%
  mutate(
    sensor = case_when(
      sensor == 'C4E' ~ 'Ponsel C4E',
      sensor == 'OPTOD' ~ 'Ponsel OPTOD',
      sensor == 'PHEHT' ~ 'Ponsel PHEHT',
      TRUE ~ sensor
    )
  )
# df_extra <- read_csv("../data/intermediate/data_voorbereidingworkshop_extra.csv") %>%
#    mutate(
#     sensor = case_when(
#       sensor == 'C4E' ~ 'Ponsel C4E',
#       sensor == 'PHEHT' ~ 'Ponsel PHEHT',
#       TRUE ~ sensor
#     ),
#     datetime = ymd_hms(datetime, tz = 'CET')
#   )
df <- filter(df_dt, parameter_combi == 'aEFFLDI003 - 2620015709 WA1111-EU-V2 / Dissolved Oxygen (OPTOD)(mg/L)')  

s_lims <- semi_join(sensor_lims, df)
lims <- c(s_lims$ondergrens, s_lims$bovengrens)

# oh <- read_csv('../data/raw/data_oh_voor_workshop.csv') %>% 
#   mutate(datum = ymd_hms(datum, tz = 'CET'))
# 
# oh_saqc <- oh %>% select(datum_start = datum) %>% 
#   mutate(maint = datum_start + hours(3))
# 
# df_lab_cor <- read_csv('../data/raw/lab_cor_data_workshop.csv') %>% 
#   mutate(datetime_compare = ymd_hms(datetime_compare, tz = 'CET'))

sensor_lims %>%
  knitr::kable(caption = glue::glue("Sensor ranges, dit geval: {s_lims$sensor} {s_lims$parameter_code} {s_lims$hoedanigheid_code}"))

par_ranges %>% select(parameter_id, jaar_min, jaar_max) %>% knitr::kable(caption = "Realistische ranges per parameter, dit geval Zuurstof_opgelost_mg/l")

p <- ggplot(df, aes(x = datetime, y = waarde)) +
  geom_point(size = 0.1, shape = 1) +
  geom_line() +
  labs(x = NULL, y = glue::glue("{s_lims$parameter_code} [{s_lims$eenheid_code}]"))
 
#ggplotly(p)
p

lims_real <-  c(0, 12)
df <- df %>% mutate(waarde_realistisch = between(waarde, left = lims_real[1], right = lims_real[2]))

ggplot(df, aes(x = datetime, y = waarde)) +
  geom_line() +
      geom_point(aes(col = waarde_realistisch, shape = waarde_realistisch)) +
      labs(x = NULL, y = "NO3-N [mg/l]") +
      scale_color_manual(values = c("red", "black")) +
      scale_shape_manual(values = c(1, NA))

df <- get_flatline(dataset = df, value_col = "waarde", threshold = 0, n_measurements = 10)


ggplot(filter(df, !is.na(flatline)), aes(x = datetime, y = waarde)) +
  geom_line() +
      geom_point(aes(col = flatline, shape = flatline)) +
      labs(x = NULL, y = "NO3-N [mg/l]") +
      scale_color_manual(values = c("black", "red")) +
      scale_shape_manual(values = c(NA, 1))

df_trans <- filter(df, flatline == FALSE, waarde_realistisch == TRUE) %>% 
  distinct(datetime, waarde) %>% 
    group_by(datetime) %>%
    summarise(waarde = mean(waarde, na.rm = TRUE)) %>%
    ungroup() %>%
    transform_data_uv(time_col = 'datetime')

  df_trans$features <- slide(df_trans, ~get_tsfeatures(.x$waarde), .before = 9, .complete = FALSE)
  df_trans <- df_trans %>% unnest_wider(features)
  
col_keep <- is.na(df_trans) %>% colMeans()  
col_keep <- names(col_keep[col_keep < 1])
col_keep <- col_keep[!col_keep %in%  c('datetime', 'waarde', 'time')] 
df_testout <- df_trans[, col_keep]

out_forest <- isotree::isolation.forest(data = df_testout, ntrees = 500)

df_testout$datetime <- df_trans$datetime
df_testout$waarde <- df_trans$waarde
df_testout$iforest_score <- predict(out_forest, df_testout)
iforest_th <- quantile(df_testout$iforest_score, probs = 0.99, na.rm = TRUE, names = FALSE)
df_testout$iforest_type <- df_testout$iforest_score >= iforest_th

#plot
ggplot(df_testout, aes(x = datetime, y = waarde)) +
  geom_line() +
      geom_point(aes(col = iforest_type, shape = iforest_type)) +
      # labs(title = par) +
      scale_color_manual(values = c("black", "red")) +
      scale_shape_manual(values = c(NA, 1))