---
title: "dendRoAnalyst: end-to-end workflow and function tour"
author: "Package vignette"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{dendRoAnalyst: end-to-end workflow and function tour}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r 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
```

## Overview

This vignette presents a complete workflow for **dendRoAnalyst** and also acts
as a compact function tour for the user-facing functions covered by the uploaded
source files. The sequence follows a typical dendrometer analysis pipeline:

1. read dendrometer and climate data,
2. check time resolution, jumps, and gaps,
3. resample and truncate,
4. smooth the signal,
5. compute daily summaries,
6. classify stem-cycle phases,
7. classify zero-growth phases,
8. perform superposed epoch analysis,
9. fit growth curves and compare methods,
10. detrend and standardize,
11. compute moving climate correlations,
12. analyze harsh-climate episodes,
13. interpolate missing values from a network, and
14. inspect periodicity with wavelet tools.

The vignette is written in a CRAN-friendly style. Lightweight sections are
executed; heavier, interactive, or computationally expensive sections are shown
as code examples with `eval = FALSE`. All the summary and plot functions can also
be executed by simply applying `summary()` and `plot()` functions. 

## Example datasets

```{r data}
data(gf_nepa17)
data(nepa17)
data(ktm_clim_hourly)
data(ktm_rain17)

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

Throughout this vignette we use:

- `gf_nepa17` for subdaily dendrometer workflows,
- `nepa17[1:1000, ]` for compact daily examples,
- `ktm_clim_hourly` for climate-feature examples, and
- `ktm_rain17` for harsh-climate examples.

```{r 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
```

## 1. Reading dendrometer and climate data

### `read.dendrometer()`

When users start from raw files, `read.dendrometer()` is the recommended entry
point.

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

This call returns an imported dendrometer table and, optionally, an import
report describing parsing and time-resolution diagnostics.

### `read.climate()`

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

The standardized climate object can later be aggregated to daily scale,
converted to subdaily climate features, or joined directly to dendrometer
outputs.

## 2. Looking for jumps or gaps in dendrometer data

### `reso_dm()`

```{r reso}
reso_dm(dm_hourly[[1]])
```

This returns the dominant sampling interval in minutes and helps diagnose
irregular logging intervals before subsequent processing.

### `jump.locator()` and `i.jump.locator()`

```{r 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)
```

```{r 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"
)
```

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

### `dm.na.interpolation()` and its plot methods

```{r 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)
```

```{r 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)
```

```{r gap-plots, fig.height=5}
plot(gap_fill)
plot(gap_fill)
plot(gap_fill, original = gap_demo)
plot(gap_fill)
```

These functions provide both the filled series and visual diagnostics of where
and how the gaps were handled.

## 3. Resample and truncate the dataset

### `dendro.resample()`

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

### `dendro.truncate()`

```{r truncate}
truncated_demo <- dendro.truncate(dm_hourly, CalYear = 2017, DOY = c(1, 120))
head(truncated_demo)
```

Resampling harmonizes time steps; truncation helps focus on one year or one
season.

## 4. Smoothing

### `smooth_dm()`

```{r 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)
```

```{r 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
)
```

Smoothing is often applied before phase classification to suppress short spikes.

## 5. Daily approach and daily climate functions

### `daily.data()` and `plot.daily_output()`

```{r daily-data}
daily_out <- daily.data(df = dm_daily_demo, TreeNum = 1)
head(daily_out)
```

```{r daily-plots}
plot(daily_out)
plot(daily_out, type = "timing")
plot(daily_out, type = "change")
```

```{r 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")
```

### `dm_daily_clim()` and `dm_join_daily_clim()`

```{r 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)
```

```{r daily-join}
daily_clim_joined <- dm_join_daily_clim(daily_out, clim_day)
head(daily_clim_joined)
```

### `dm_add_climate()` for daily outputs

```{r 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)
```

## 6. Stem-cycle approach and associated functions

### `phase.sc()` and `plot.SC_output()`

```{r 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)
```

```{r sc-plots}
plot(sc_out)
plot(sc_out, type = "ribbon")
plot(sc_out, type = "increment")
```

```{r 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")
```

### `dm_subdaily_clim()`, `dm_join_subdaily_clim()`, `dm_join_phase_clim()`, and `dm_add_climate()`

```{r 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)
```

```{r 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)
```

## 7. Zero-growth approach and associated functions

### `phase.zg()`, `twd.maxima()`, and `plot.ZG_output()`

```{r 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)
```

```{r twd-maxima}
max_twd_tbl <- twd.maxima(dm_phase_demo, TreeNum = 1, smoothing = 5)
head(max_twd_tbl)
```

```{r zg-plots}
plot(zg_out)
plot(zg_out, temporal = "daily")
plot(zg_out, temporal = "monthly", type = "twd")
```

```{r 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")
```

### Climate joins and climate-phase plots

```{r 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)
```

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

```{r 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")
)
```

## 8. Superposed epoch and associated functions

### `dm_event_times()`, `dm_epoch_extract()`, `dm_epoch_test()`, and method dispatch

```{r epoch}
event_tbl <- dm_event_times(zg_out, event = "all")
head(event_tbl)
```

```{r 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)
```

```{r 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)
```

## 9. Event-based climate functions

```{r 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)
```

```{r 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)
```

```{r 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)
```

```{r 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"
)
```

## 10. Growth fit and evaluation

### `dm.fit.gompertz()`

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

### `dm.growth.fit()`, `plot.dm_growth_fit()`, and print/summary methods

```{r 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()
```

```{r 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")
```

### `dm.growth.fit.double()` and `dm.growth.evaluate()`

```{r 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)
```

```{r 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)
```

```{r growth-evaluate-plots, eval = FALSE}
plot(growth_eval, metric = "rmse", type = "boxplot")
plot(growth_eval, metric = "r2", type = "heatmap")
```

## 11. Detrending and associated functions

```{r 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")
```

```{r mean-detrended, eval = FALSE}
mean_det <- mean_detrended.dm(detrended_fit)
head(mean_det)
plot(mean_det)
```

### `dm_standardize()` and `plot.dm_standardized()`

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

```{r standardize-plot}
plot(std_nh)
```

## 12. Moving correlations and associated functions

```{r 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()
```

```{r 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")
```

```{r 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)
```

## 13. Harsh climate using `clim.twd()` and associated functions

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

```{r 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")
```

```{r 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)
```

## 14. Network interpolation

### `network.interpolation()` and `plot.network_interpolation()`

```{r 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")
```

```{r 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")
```

## 15. Wavelet analysis and associated functions

### `dm_wavelet()`, `summary.dm_wavelet()`, and `plot.dm_wavelet()`

```{r 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()
```

```{r 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")
```

### `dm_wavelet_reconstruct()` and its methods

```{r 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()
```

```{r wavelet-reconstruct-plots, eval = FALSE}
plot(rec_extract)
plot(rec_remove)
```

### `dm_wavelet_coherence()` and `plot.dm_wavelet_coherence()`

```{r 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
)
```

```{r 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")
```

## Closing note

This vignette is intended to be both a workflow guide and a compact reference
for the main user-facing functions in the current **dendRoAnalyst** 
analysis pipeline. For classroom teaching, the companion demonstration script
can be run section by section; for package documentation, this vignette keeps
only the lighter sections evaluated and leaves heavy or interactive examples in
place for users to run locally.

## For suggestions, comments and questions please contact : <sugam.aryal@fau.de>