## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5,
  warning = FALSE,
  message = FALSE
)

library(dendRoAnalyst)

RUN_HEAVY <- FALSE
RUN_INTERACTIVE <- FALSE

## ----data---------------------------------------------------------------------
data(gf_nepa17)
data(nepa17)
data(ktm_clim_hourly)
data(ktm_rain17)

head(gf_nepa17[, 1:3])
head(ktm_clim_hourly)

## ----objects------------------------------------------------------------------
dm_hourly <- gf_nepa17
dm_daily_demo <- nepa17[1:1000, ]
dm_phase_demo <- gf_nepa17[1:800, ]
clim_hourly <- ktm_clim_hourly
clim_rain <- ktm_rain17

## ----read-dendrometer, eval = FALSE-------------------------------------------
# raw_dm <- read.dendrometer(
#   file = "inst/extdata/my_dendrometer_file.csv",
#   detect_resolution = TRUE,
#   return_report = TRUE,
#   quiet = FALSE
# )

## ----read-climate-------------------------------------------------------------
clim_std <- read.climate(
  x = clim_hourly,
  time_col = "TIME",
  vars = c("temp", "prec", "VPD", "RH"),
  verbose = FALSE
)
head(clim_std)

## ----reso---------------------------------------------------------------------
reso_dm(dm_hourly[[1]])

## ----jump-manual--------------------------------------------------------------
jump_free_manual <- jump.locator(
  df = dm_daily_demo[, 1:3],
  detection_method = "manual",
  manual_threshold = 1,
  adjustment_method = "window_median"
)
head(jump_free_manual)

## ----jump-auto, eval = FALSE--------------------------------------------------
# jump_free_auto <- jump.locator(
#   df = dm_daily_demo[, 1:3],
#   detection_method = "auto",
#   auto_method_penalty = 10,
#   adjustment_method = "window_median"
# )

## ----jump-interactive, eval = FALSE-------------------------------------------
# # Interactive graphical workflow
# if (RUN_INTERACTIVE) {
#   i.jump.locator(df = dm_daily_demo[, 1:3])
# }

## ----gap-detect---------------------------------------------------------------
gap_demo <- dm_daily_demo[, 1:3]
gap_demo[150:170, 2] <- NA

gap_info <- dm.na.interpolation(
  df = gap_demo,
  resolution = 60,
  fill = FALSE,
  assess = FALSE
)
head(gap_info$gap_info)

## ----gap-fill-----------------------------------------------------------------
gap_fill <- dm.na.interpolation(
  df = gap_demo,
  resolution = 60,
  fill = TRUE,
  method = "spline",
  assess = TRUE,
  assess_lengths_hours = c(6, 12)
)
head(gap_fill$data)

## ----gap-plots, fig.height=5--------------------------------------------------
plot(gap_fill)
plot(gap_fill)
plot(gap_fill, original = gap_demo)
plot(gap_fill)

## ----resample-----------------------------------------------------------------
resampled_hourly <- dendro.resample(
  df = dm_hourly,
  by = "H",
  value = "mean",
  method = "auto"
)
head(resampled_hourly)

## ----truncate-----------------------------------------------------------------
truncated_demo <- dendro.truncate(dm_hourly, CalYear = 2017, DOY = c(1, 120))
head(truncated_demo)

## ----smooth-------------------------------------------------------------------
res_minutes <- reso_dm(dm_phase_demo[[1]])
smoothed_medmean <- smooth_dm(
  time = dm_phase_demo[[1]],
  dm = dm_phase_demo[[2]],
  resolution_min = res_minutes,
  method = "median_mean",
  window_hours = 3
)
head(smoothed_medmean)

## ----smooth-other, eval = FALSE-----------------------------------------------
# smoothed_pspline <- smooth_dm(
#   time = dm_phase_demo[[1]],
#   dm = dm_phase_demo[[2]],
#   resolution_min = res_minutes,
#   method = "pspline",
#   window_hours = 6
# )

## ----daily-data---------------------------------------------------------------
daily_out <- daily.data(df = dm_daily_demo, TreeNum = 1)
head(daily_out)

## ----daily-plots--------------------------------------------------------------
plot(daily_out)
plot(daily_out, type = "timing")
plot(daily_out, type = "change")

## ----daily-more-plots, eval = FALSE-------------------------------------------
# plot(daily_out, type = "timing_violin", Year = 2017, DOY = c(1, 6))
# plot(daily_out, type = "boxplot", stat = "amplitude", by = "month_of_year")
# plot(daily_out, type = "heatmap")

## ----daily-clim---------------------------------------------------------------
clim_day <- dm_daily_clim(
  clim_std,
  mean_vars = c("temp", "VPD", "RH"),
  max_vars  = c("VPD"),
  sum_vars  = c("prec"),
  lag_vars = c("VPD_max", "temp_mean"),
  lagmean_vars = c("temp_mean", "VPD_mean"),
  lagsum_vars = c("prec_sum"),
  lag_days = c(1, 3, 7)
)
head(clim_day)

## ----daily-join---------------------------------------------------------------
daily_clim_joined <- dm_join_daily_clim(daily_out, clim_day)
head(daily_clim_joined)

## ----add-climate-daily--------------------------------------------------------
daily_clim_added <- dm_add_climate(
  daily_out,
  clim_hourly,
  scale = "daily",
  mean_vars = c("temp", "VPD", "RH"),
  max_vars  = c("temp", "VPD"),
  sum_vars  = c("prec")
)
head(daily_clim_added)

## ----sc-out-------------------------------------------------------------------
sc_out <- phase.sc(df = dm_phase_demo, TreeNum = 1)
sc_out_sm <- phase.sc(df = dm_phase_demo, TreeNum = 1, smoothing = 12)
head(sc_out$SC_phase)
head(sc_out$SC_cycle)

## ----sc-plots-----------------------------------------------------------------
plot(sc_out)
plot(sc_out, type = "ribbon")
plot(sc_out, type = "increment")

## ----sc-more-plots, eval = FALSE----------------------------------------------
# plot(sc_out, type = "transition", x_axis = "doy", Year = 2017, DOY = c(50, 100))
# plot(sc_out, type = "balance", temporal = "daily")
# plot(sc_out, type = "frequency", temporal = "monthly")
# plot(sc_out, type = "heatmap", temporal = "daily")
# plot(sc_out, type = "boxplot", stat = "Magnitude", temporal = "monthly")

## ----subdaily-clim------------------------------------------------------------
clim_sub <- dm_subdaily_clim(
  clim_hourly,
  mean_vars = c("temp", "VPD", "RH"),
  sum_vars  = c("prec"),
  lag_vars  = c("temp", "VPD", "RH"),
  roll_hours = c(1, 3, 6, 24),
  lag_hours  = c(1, 3, 6, 24)
)
head(clim_sub)

## ----sc-joins-----------------------------------------------------------------
sc_phase_clim <- dm_join_phase_clim(
  sc_out,
  clim_hourly,
  mean_vars = c("temp", "VPD", "RH"),
  max_vars  = c("temp", "VPD"),
  sum_vars  = c("prec")
)
head(sc_phase_clim$SC_cycle)

sc_point_clim <- dm_join_subdaily_clim(sc_out_sm, clim_sub)
head(sc_point_clim$SC_phase)

sc_added_clim <- dm_add_climate(
  sc_out_sm,
  clim_hourly,
  scale = "subdaily",
  sub_mean_vars = c("temp", "VPD", "RH"),
  sub_sum_vars  = c("prec"),
  sub_lag_vars  = c("temp", "VPD", "RH"),
  roll_hours = c(1, 3, 6, 24),
  lag_hours  = c(1, 3, 6, 24)
)
head(sc_added_clim$SC_phase)

## ----zg-out-------------------------------------------------------------------
zg_out <- phase.zg(df = dm_phase_demo, TreeNum = 1)
zg_out_sm <- phase.zg(df = dm_phase_demo, TreeNum = 1, smoothing = 6)
head(zg_out$ZG_phase)
head(zg_out$ZG_cycle)

## ----twd-maxima---------------------------------------------------------------
max_twd_tbl <- twd.maxima(dm_phase_demo, TreeNum = 1, smoothing = 5)
head(max_twd_tbl)

## ----zg-plots-----------------------------------------------------------------
plot(zg_out)
plot(zg_out, temporal = "daily")
plot(zg_out, temporal = "monthly", type = "twd")

## ----zg-more-plots, eval = FALSE----------------------------------------------
# plot(zg_out, temporal = "monthly", type = "abr")
# plot(zg_out, type = "transition", x_axis = "doy", Year = 2017, DOY = c(50, 100))
# plot(zg_out, type = "boxplot", stat = "max.twd", temporal = "daily")

## ----zg-clim-joins------------------------------------------------------------
zg_phase_clim <- dm_join_phase_clim(
  zg_out,
  clim_hourly,
  mean_vars = c("temp", "VPD", "RH"),
  max_vars  = c("temp", "VPD"),
  sum_vars  = c("prec")
)
head(zg_phase_clim$ZG_cycle)

zg_point_clim <- dm_join_subdaily_clim(zg_out_sm, clim_sub)
head(zg_point_clim$ZG_phase)

zg_phase_added <- dm_add_climate(
  zg_out,
  clim_hourly,
  scale = "phase",
  mean_vars = c("temp", "VPD", "RH"),
  max_vars  = c("temp", "VPD"),
  sum_vars  = c("prec")
)
head(zg_phase_added$ZG_cycle)

## ----climate-phase-plots------------------------------------------------------
plot(zg_phase_added, climate_var = "VPD_mean_phase", scale = "cycle", type = "boxplot", temporal = 'daily')

## ----climate-phase-compare, eval = FALSE--------------------------------------
# plot(zg_phase_added, climate_var = "temp_mean_phase", scale = "cycle", type = "timeseries")
# 
# plot(
#   zg_phase_added,
#   climate_vars = c("temp_mean_phase", "VPD_mean_phase", "RH_mean_phase"),
#   scale = "cycle",
#   type = "boxplot"
# )
# 
# plot(
#   zg_phase_added,
#   scale = "cycle",
#   type = "cor_heatmap",
#   numeric_vars = c("Duration_h", "max.twd", "ABr.value",
#                    "temp_mean_phase", "VPD_mean_phase", "RH_mean_phase")
# )

## ----epoch--------------------------------------------------------------------
event_tbl <- dm_event_times(zg_out, event = "all")
head(event_tbl)

## ----epoch-extract, eval = FALSE----------------------------------------------
# epoch_extract <- dm_epoch_extract(
#   x = zg_out,
#   clim = clim_hourly,
#   event = "phase_start",
#   phase = "all",
#   response_var = c("temp", "VPD"),
#   before = 24,
#   after = 24,
#   unit = "hours"
# )
# head(epoch_extract)

## ----epoch-test, eval = FALSE-------------------------------------------------
# epoch_test <- dm_epoch_test(
#   x = zg_out,
#   clim = clim_hourly,
#   event = "phase_start",
#   phase = "all",
#   response_var = c("temp", "VPD"),
#   before = 24,
#   after = 24,
#   unit = "hours",
#   n_boot = 99
# )
# 
# epoch_summary <- summary(epoch_test)   # summary.dm_epoch()
# print(epoch_summary)                   # print.summary.dm_epoch()
# plot(epoch_test)

## ----event-climate, eval = FALSE----------------------------------------------
# evt_zg <- dm_event_climate(
#   zg_out,
#   clim_hourly,
#   event = "phase_start",
#   phase = "all",
#   windows = c(0, 3, 6, 12, 24),
#   mean_vars = c("temp", "VPD", "RH"),
#   max_vars  = c("VPD"),
#   sum_vars  = c("prec")
# )
# head(evt_zg)

## ----event-climate-summary, eval = FALSE--------------------------------------
# evt_summary <- dm_event_climate_summary(
#   evt_zg,
#   climate_var = "VPD_mean_prev_6h",
#   group_var = "Phase",
#   by = "month_of_year"
# )
# head(evt_summary)

## ----event-climate-test, eval = FALSE-----------------------------------------
# evt_test <- dm_event_climate_test(
#   evt_zg,
#   climate_var = "VPD_mean_prev_6h",
#   group_var = "Phase",
#   by = "month_of_year",
#   method = "auto"
# )
# head(evt_test$summary)
# head(evt_test$overall_tests)
# head(evt_test$pairwise_tests)

## ----event-climate-plots, eval = FALSE----------------------------------------
# plot(
#   evt_zg,
#   climate_var = "VPD_mean_prev_6h",
#   group_var = "Phase",
#   facet_by = "month_of_year",
#   geom = "both",
#   add_test = TRUE
# )
# 
# evt_maxtwd <- dm_event_climate(
#   zg_out,
#   clim_hourly,
#   event = "max_twd",
#   windows = c(0, 3, 6, 12, 24),
#   mean_vars = c("temp", "VPD", "RH"),
#   max_vars  = c("VPD"),
#   sum_vars  = c("prec")
# )
# 
# plot(
#   evt_maxtwd,
#   climate_var = "VPD_mean_prev_6h",
#   response_var = "max.twd",
#   group_var = "Phase"
# )

## ----legacy-gompertz, eval = FALSE--------------------------------------------
# gomp_legacy <- dm.fit.gompertz(df = dm_hourly, TreeNum = 1, CalYear = 2017, verbose = FALSE)
# head(gomp_legacy)

## ----growth-fit, eval = FALSE-------------------------------------------------
# fit_single <- dm.growth.fit(
#   df = dm_hourly,
#   TreeNum = 1:2,
#   method = "gompertz",
#   year_mode = "yearly",
#   verbose = FALSE
# )
# 
# print(fit_single)                 # print.dm_growth_fit()
# fit_single_summary <- summary(fit_single)   # summary.dm_growth_fit()
# print(fit_single_summary)         # print.summary.dm_growth_fit()

## ----growth-fit-plots, eval = FALSE-------------------------------------------
# plot(fit_single, type = "fit")
# plot(fit_single, type = "season")
# plot(fit_single, type = "residuals")
# plot(fit_single, type = "timing")
# plot(fit_single, type = "overlay")
# plot(fit_single, type = "parameters")

## ----growth-fit-double, eval = FALSE------------------------------------------
# fit_double <- dm.growth.fit.double(
#   df = dm_hourly,
#   TreeNum = 1:2,
#   method = "gompertz",
#   year_mode = "yearly",
#   verbose = FALSE
# )
# print(fit_double)

## ----growth-evaluate, eval = FALSE--------------------------------------------
# growth_eval <- dm.growth.evaluate(
#   df = dm_hourly,
#   TreeNum = 1:2,
#   methods = c("gompertz", "loess", "spline", "double_gompertz"),
#   year_mode = "yearly",
#   verbose = FALSE
# )
# head(growth_eval)

## ----growth-evaluate-plots, eval = FALSE--------------------------------------
# plot(growth_eval, metric = "rmse", type = "boxplot")
# plot(growth_eval, metric = "r2", type = "heatmap")

## ----detrend, eval = FALSE----------------------------------------------------
# detrended_fit <- dm.detrend.fit(fit_single)
# plot(detrended_fit, type = "compare")
# plot(detrended_fit, type = "fit")
# plot(detrended_fit, type = "residual")
# plot(detrended_fit, type = "detrended")
# plot(detrended_fit, type = "boxplot")

## ----mean-detrended, eval = FALSE---------------------------------------------
# mean_det <- mean_detrended.dm(detrended_fit)
# head(mean_det)
# plot(mean_det)

## ----standardize--------------------------------------------------------------
std_nh <- dm_standardize(
  df = dm_hourly,
  season_type = "NH",
  method = "robust_amplitude"
)
head(std_nh$data)
head(std_nh$parameters)

## ----standardize-plot---------------------------------------------------------
plot(std_nh)

## ----mov-cor, eval = FALSE----------------------------------------------------
# mov_out <- mov.cor.dm(
#   df = dm_hourly,
#   Clim = clim_hourly,
#   TreeNum = 1,
#   win_size = 15,
#   cor_method = "pearson",
#   boot = FALSE
# )
# 
# print(mov_out)                # print.mov_cor_dm()
# mov_summary <- summary(mov_out)   # summary.mov_cor_dm()
# print(mov_summary)            # print.summary_mov_cor_dm()

## ----mov-cor-plots, eval = FALSE----------------------------------------------
# plot(mov_out, type = "heatmap")
# plot(mov_out, type = "line", annotate_peak = TRUE)
# plot(mov_out, type = "facet")
# plot(mov_out, type = "peak")
# plot(mov_out, type = "line")
# plot(mov_summary, type = "peak_corr")

## ----mov-cor-boot, eval = FALSE-----------------------------------------------
# mov_boot <- mov.cor.dm(
#   df = dm_hourly,
#   Clim = clim_hourly,
#   TreeNum = 1,
#   win_size = 15,
#   cor_method = "pearson",
#   boot = TRUE,
#   R = 99,
#   set_seed = 1
# )
# plot(mov_boot, type = "line", show_ci = TRUE)

## ----clim-twd, eval = FALSE---------------------------------------------------
# rel_out <- clim.twd(
#   df = dm_hourly,
#   Clim = clim_rain,
#   thresholdClim = "<10",
#   thresholdDays = ">5",
#   showPlot = FALSE
# )

## ----clim-twd-stats, eval = FALSE---------------------------------------------
# clim_stats <- clim.twd.stats(
#   rel_out,
#   response = "full_period_change",
#   group_by = "tree"
# )
# 
# clim_stats_summary <- summary(clim_stats)   # summary.clim_twd_stats()
# print(clim_stats_summary)                   # print.summary.clim_twd_stats()
# plot(clim_stats, type = "trajectory")

## ----clim-twd-test, eval = FALSE----------------------------------------------
# tree_info <- data.frame(
#   tree = c("T2", "T3"),
#   species = c("Pinus", "Pinus"),
#   site = c("Ktm", "Ktm")
# )
# 
# clim_test <- clim.twd.test(
#   rel_out,
#   tree_info = tree_info,
#   parameter = "lag_to_return_zero",
#   compare_by = "species",
#   within = "year"
# )
# 
# clim_test_summary <- summary(clim_test)   # summary.clim_twd_test()
# print(clim_test_summary)                  # print.summary.clim_twd_test()
# plot(clim_test)

## ----network-interp, eval = FALSE---------------------------------------------
# df1 <- dm_hourly
# df1[40:50, "T2"] <- NA
# 
# ref_net <- cbind(dm_hourly, dm_hourly[, 2:3], dm_hourly[, 2:3])
# colnames(ref_net) <- c("Time", "T1", "T2", "T3", "T4", "T5", "T6")
# 
# net_interp <- network.interpolation(
#   df1,
#   ref_net,
#   niMethod = "proportional",
#   n_boot = 100,
#   return_flags = TRUE,
#   return_pi = TRUE,
#   assess = TRUE,
#   assess_lengths_steps = c(1, 2, 4),
#   correct_gap_jumps = TRUE,
#   jump_threshold = 0.05
# )
# head(net_interp)
# # to get the assessment stats we can use attr() function
# attr(out, "network_validation_summary")
# attr(out, "network_validation_points")
# attr(out, "network_validation_seasonal")

## ----network-interp-plots, eval = FALSE---------------------------------------
# plot(net_interp)
# plot(net_interp, type = "uncertainty")
# plot(net_interp, type = "availability")
# plot(net_interp, type = "summary", summary_metric = "n_gap_jumps_removed")
# plot(net_interp, type = "validation", validation_metric = "MAPE")
# plot(net_interp, type = "coverage")
# plot(net_interp, type = "compare")
# plot(net_interp, type = "seasonal_error", seasonal_metric = "mean_abs_error")

## ----wavelet, eval = FALSE----------------------------------------------------
# wv <- dm_wavelet(
#   x = dm_hourly,
#   TreeNum = 1:2,
#   source = "raw",
#   make_pval = TRUE,
#   verbose = FALSE
# )
# 
# wv_summary <- summary(wv)   # summary.dm_wavelet()
# print(wv_summary)           # print.summary.dm_wavelet()

## ----wavelet-plots, eval = FALSE----------------------------------------------
# plot(wv, type = "series")
# plot(wv, type = "average")
# plot(wv, type = "power")
# plot(wv, series = names(wv$results)[1], type = "power")

## ----wavelet-reconstruct, eval = FALSE----------------------------------------
# rec_extract <- dm_wavelet_reconstruct(
#   wv,
#   mode = "extract",
#   lower_hours = 20,
#   upper_hours = 28,
#   only_sig = TRUE
# )
# 
# rec_remove <- dm_wavelet_reconstruct(
#   wv,
#   mode = "remove",
#   lower_hours = 20,
#   upper_hours = 28,
#   only_sig = TRUE
# )
# 
# rec_extract_summary <- summary(rec_extract)   # summary.dm_wavelet_reconstruct()
# print(rec_extract_summary)                    # print.summary.dm_wavelet_reconstruct()

## ----wavelet-reconstruct-plots, eval = FALSE----------------------------------
# plot(rec_extract)
# plot(rec_remove)

## ----wavelet-coherence, eval = FALSE------------------------------------------
# wc <- dm_wavelet_coherence(
#   dm = dm_hourly,
#   clim = clim_hourly,
#   TreeNum = 1,
#   clim_vars = c("temp", "VPD"),
#   source = "raw",
#   dm_fun = "mean",
#   clim_funs = c(temp = "mean", VPD = "mean"),
#   loess_span = 0,
#   make.pval = TRUE,
#   n.sim = 10,
#   verbose = FALSE
# )

## ----wavelet-coherence-plots, eval = FALSE------------------------------------
# plot(wc)
# plot(wc, type = "cross_power")
# plot(wc, type = "average_coherence")
# plot(wc, type = "phase")
# plot(wc, type = "series")