R/plot_SC_output.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' @title Plot method for stem-cycle output
#'
#' @description
#' Unified S3 plotting method for objects returned by \code{phase.sc()}.
#' Supports raw phase timelines, phase ribbons, transition diagnostics,
#' daily/monthly phase balance, cumulative increment, event frequency,
#' phase heatmaps, and daily/monthly boxplots of phase statistics.
#'
#' For \code{type = "balance"}, two plotting modes are available. With
#' \code{balance_mode = "time_of_day"}, phases are drawn as continuous
#' within-day intervals, and the y-axis represents hour of day from 0 to 24.
#' This makes the timing of shrinkage, expansion, and increment visible. For
#' example, if shrinkage occurs from 00:00 to 05:00, expansion from 05:00 to
#' 15:00, and shrinkage again from 15:00 to 24:00, these intervals are shown
#' as continuous coloured blocks along the y-axis.
#'
#' With \code{balance_mode = "duration"}, the older stacked-bar behaviour is
#' used. In that case, the y-axis represents the amount of time spent in each
#' phase. With \code{temporal = "daily"}, the y-axis is hours per day. With
#' \code{temporal = "monthly"}, the y-axis is hours per month.
#'
#' @param x Object of class \code{"SC_output"} returned by \code{phase.sc()}.
#' @param y Unused.
#' @param DOY Optional numeric vector of length 2 giving start and end
#'   day-of-year.
#' @param Year Optional numeric year used together with \code{DOY}.
#' @param type Plot type. One of
#'   \code{"points"}, \code{"ribbon"}, \code{"transition"},
#'   \code{"balance"}, \code{"increment"}, \code{"frequency"},
#'   \code{"heatmap"}, or \code{"boxplot"}.
#' @param temporal Temporal scale: \code{"raw"}, \code{"daily"}, or
#'   \code{"monthly"}. For \code{type = "boxplot"}, use \code{"daily"} or
#'   \code{"monthly"}. For \code{type = "balance"}, \code{"frequency"}, and
#'   \code{"heatmap"}, \code{"raw"} is internally treated as daily summary and
#'   a warning is emitted.
#' @param x_axis X-axis style for time-based plots. One of \code{"time"} or
#'   \code{"doy"}. When \code{"doy"} is used, timestamps are converted to
#'   day-of-year on the x-axis. For \code{temporal = "monthly"}, \code{"doy"}
#'   is not meaningful and time is used instead.
#' @param balance_mode For \code{type = "balance"}, controls how phase balance
#'   is drawn. One of \code{"time_of_day"} or \code{"duration"}.
#'   \code{"time_of_day"} draws continuous phase intervals along the y-axis
#'   from 0 to 24 h, preserving when phases occurred within each day.
#'   \code{"duration"} draws stacked bars showing total time spent in each
#'   phase per day or month.
#' @param stat For \code{type = "boxplot"}, one of
#'   \code{"Duration_h"}, \code{"Duration_m"}, \code{"Magnitude"}, or
#'   \code{"rate"}.
#' @param phase For \code{type = "boxplot"}, one of
#'   \code{"all"}, \code{"Shrinkage"}, \code{"Expansion"}, or
#'   \code{"Increment"}.
#' @param cols Vector of three colours for Shrinkage, Expansion, and Increment.
#' @param phNames Vector of three labels for the three phase names.
#' @param transition_linetype Line type used for transition markers in
#'   \code{type = "transition"}.
#' @param transition_alpha Transparency of transition markers.
#' @param singleton_as_points Logical. If \code{TRUE}, the boxplot panel falls
#'   back to a point plot when every group contains only one value.
#' @param ... Unused.
#'
#' @return A \code{ggplot2} object, returned invisibly.
#'
#' @details
#' Plot types:
#' \itemize{
#'   \item \code{"points"}: dendrometer series with points coloured by phase.
#'   \item \code{"ribbon"}: dendrometer series with background ribbons for
#'         contiguous phases.
#'   \item \code{"transition"}: dendrometer series with vertical lines at phase
#'         changes.
#'   \item \code{"balance"}: phase balance through time. With
#'         \code{balance_mode = "time_of_day"}, phase intervals are drawn
#'         continuously along the y-axis from 0 to 24 h, so the timing of
#'         shrinkage, expansion, and increment within each day is visible.
#'         With \code{balance_mode = "duration"}, stacked bars show total time
#'         spent in each phase per day or month.
#'   \item \code{"increment"}: cumulative increment, phase 3, over time.
#'   \item \code{"frequency"}: number of phase events starting per day or month.
#'   \item \code{"heatmap"}: dominant phase by hour-of-day across dates or months.
#'   \item \code{"boxplot"}: distribution of selected cycle statistics after
#'         assembling cycle summaries to the chosen daily or monthly time scale.
#' }
#'
#' The \code{"frequency"} plot uses \code{SC_cycle}.
#' The \code{"heatmap"} plot uses \code{SC_phase}.
#' The \code{"transition"} plot is especially useful for checking whether
#' smoothing reduced spurious phase switching.
#'
#' @examples
#' \donttest{
#' # data(gf_nepa17)
#' # sc <- phase.sc(df = gf_nepa17, TreeNum = 1, smoothing = 12)
#'
#' # plot(sc)
#' # plot(sc, type = "ribbon")
#' # plot(sc, type = "transition")
#'
#' # Daily within-day phase timing
#' # plot(sc, type = "balance", temporal = "daily")
#'
#' # Daily within-day phase timing for selected DOY range
#' # plot(sc, type = "balance", temporal = "daily",
#' #      DOY = c(150, 160), Year = 2023)
#'
#' # Old stacked-duration balance plot
#' # plot(sc, type = "balance", temporal = "daily",
#' #      balance_mode = "duration")
#'
#' # Monthly stacked duration
#' # plot(sc, type = "balance", temporal = "monthly",
#' #      balance_mode = "duration")
#'
#' # Boxplots
#' # plot(sc, type = "boxplot", stat = "Magnitude", temporal = "monthly")
#' # plot(sc, type = "boxplot", stat = "Duration_h",
#' #      temporal = "daily", phase = "Shrinkage")
#' }
#'
#' @method plot SC_output
#' @importFrom dplyr mutate rename select group_by summarise arrange filter
#'   slice_max ungroup n bind_rows %>%
#' @importFrom tibble tibble
#' @importFrom lubridate floor_date hour minute second yday month year
#' @importFrom ggplot2 ggplot aes geom_line geom_point geom_rect geom_col
#'   geom_boxplot geom_tile geom_vline facet_wrap theme_minimal theme_bw theme
#'   element_text labs scale_color_manual scale_fill_manual scale_x_date
#'   scale_x_discrete scale_y_continuous position_jitter
#' @export
plot.SC_output <- function(x,
                           y = NULL,
                           DOY = NULL,
                           Year = NULL,
                           type = c("points", "ribbon", "transition",
                                    "balance", "increment", "frequency",
                                    "heatmap", "boxplot"),
                           temporal = c("raw", "daily", "monthly"),
                           x_axis = c("time", "doy"),
                           balance_mode = c("time_of_day", "duration"),
                           stat = c("Duration_h", "Duration_m", "Magnitude", "rate"),
                           phase = c("all", "Shrinkage", "Expansion", "Increment"),
                           cols = c("#fee8c8", "#fdbb84", "#e34a33"),
                           phNames = c("Shrinkage", "Expansion", "Increment"),
                           transition_linetype = "dashed",
                           transition_alpha = 0.55,
                           singleton_as_points = TRUE,
                           ...) {

  TIME <- dm <- Phases <- Phase <- period <- hours <- value <- group_label <- n <- NULL
  Start <- End <- Duration_h <- Duration_m <- Magnitude <- rate <- hour <- xval <- NULL
  xmin <- xmax <- cum_increment <- X <- group_date <- ymin <- ymax <- NULL

  if (!inherits(x, "SC_output")) {
    stop("Input must be an object of class 'SC_output'.")
  }

  if (!is.character(phNames) || length(phNames) != 3) {
    stop("'phNames' must be a character vector of length 3.")
  }

  if (length(cols) != 3) {
    stop("'cols' must be a vector of length 3.")
  }

  type <- match.arg(type)
  temporal <- match.arg(temporal)
  x_axis <- match.arg(x_axis)
  balance_mode <- match.arg(balance_mode)
  stat <- match.arg(stat)
  phase <- match.arg(phase)

  if (type == "boxplot" && temporal == "raw") {
    stop("For type = 'boxplot', temporal must be 'daily' or 'monthly'.")
  }

  .time_unit <- function(x) {
    switch(
      x,
      raw = NA_character_,
      daily = "day",
      monthly = "month"
    )
  }

  .ensure_posixct <- function(tt) {
    if (inherits(tt, "POSIXct")) {
      return(tt)
    }

    if (inherits(tt, "Date")) {
      return(as.POSIXct(tt))
    }

    out <- suppressWarnings(lubridate::ymd_hms(tt, quiet = TRUE))

    if (all(is.na(out))) {
      out <- suppressWarnings(
        lubridate::parse_date_time(
          tt,
          orders = c(
            "ymd HMS", "ymd HM", "ymd",
            "dmy HMS", "dmy HM", "dmy",
            "mdy HMS", "mdy HM", "mdy"
          ),
          quiet = TRUE
        )
      )
    }

    out
  }

  .resolution_hours <- function(tt) {
    tt <- tt[!is.na(tt)]
    tt <- sort(unique(tt))

    if (length(tt) < 2) {
      stop("Could not determine time resolution because fewer than two timestamps are available.")
    }

    res_h <- stats::median(
      as.numeric(diff(tt), units = "hours"),
      na.rm = TRUE
    )

    if (!is.finite(res_h) || res_h <= 0) {
      stop("Could not determine a valid positive time resolution from the stem-cycle data.")
    }

    res_h
  }

  .hour_decimal <- function(tt) {
    lubridate::hour(tt) +
      lubridate::minute(tt) / 60 +
      lubridate::second(tt) / 3600
  }

  .split_phase_interval_by_day <- function(t0, t1, ph) {
    if (is.na(t0) || is.na(t1) || is.na(ph) || t1 <= t0) {
      return(NULL)
    }

    tz_use <- attr(t0, "tzone")

    if (is.null(tz_use) || length(tz_use) == 0 || is.na(tz_use)) {
      tz_use <- "UTC"
    }

    pieces <- list()
    ii <- 1L
    cur <- t0

    while (cur < t1) {
      cur_date <- as.Date(cur, tz = tz_use)
      next_midnight <- as.POSIXct(cur_date + 1, tz = tz_use)

      seg_end <- if (t1 < next_midnight) {
        t1
      } else {
        next_midnight
      }

      ymin <- .hour_decimal(cur)

      ymax <- if (as.Date(seg_end, tz = tz_use) > cur_date &&
                  abs(.hour_decimal(seg_end)) < 1e-8) {
        24
      } else {
        .hour_decimal(seg_end)
      }

      if (is.finite(ymin) && is.finite(ymax) && ymax > ymin) {
        pieces[[ii]] <- data.frame(
          period = cur_date,
          ymin = ymin,
          ymax = ymax,
          Phases = ph
        )
        ii <- ii + 1L
      }

      cur <- seg_end
    }

    if (length(pieces) == 0) {
      return(NULL)
    }

    dplyr::bind_rows(pieces)
  }

  .build_phase_segments <- function(td_raw) {
    td2 <- td_raw %>%
      dplyr::filter(!is.na(TIME), !is.na(Phases)) %>%
      dplyr::arrange(TIME)

    if (nrow(td2) < 1) {
      stop("No non-missing phase records available for the balance plot.")
    }

    res_h <- .resolution_hours(td2$TIME)
    res_s <- res_h * 3600

    next_time <- c(
      td2$TIME[-1],
      td2$TIME[nrow(td2)] + res_s
    )

    gap_s <- as.numeric(next_time - td2$TIME, units = "secs")

    bad_gap <- !is.finite(gap_s) |
      gap_s <= 0 |
      gap_s > 1.5 * res_s

    if (any(bad_gap)) {
      next_time[bad_gap] <- td2$TIME[bad_gap] + res_s
    }

    seg_list <- vector("list", nrow(td2))

    for (ii in seq_len(nrow(td2))) {
      seg_list[[ii]] <- .split_phase_interval_by_day(
        t0 = td2$TIME[ii],
        t1 = next_time[ii],
        ph = td2$Phases[ii]
      )
    }

    seg <- dplyr::bind_rows(seg_list)

    if (nrow(seg) == 0) {
      stop("Could not construct phase intervals for the balance plot.")
    }

    seg$Phase <- .phase_factor(seg$Phases)

    seg
  }

  .x_values <- function(tt, temporal, x_axis) {
    if (x_axis == "time") {
      return(list(x = tt, xlab = "Time"))
    }

    if (temporal == "monthly") {
      warning("x_axis = 'doy' is not meaningful for temporal = 'monthly'; using time instead.")
      return(list(x = tt, xlab = "Time"))
    }

    list(
      x = lubridate::yday(tt),
      xlab = "Day of year"
    )
  }

  .phase_factor <- function(x) {
    factor(x, levels = c(1, 2, 3), labels = phNames)
  }

  .safe_sum <- function(z) {
    if (all(is.na(z))) {
      NA_real_
    } else {
      sum(z, na.rm = TRUE)
    }
  }

  .aggregate_cycles <- function(cyc, temporal) {
    if (temporal == "raw") {
      return(cyc)
    }

    unit <- .time_unit(temporal)

    cyc %>%
      dplyr::filter(!is.na(Start), !is.na(Phases)) %>%
      dplyr::mutate(period = lubridate::floor_date(Start, unit = unit)) %>%
      dplyr::group_by(period, Phases) %>%
      dplyr::summarise(
        Start = min(Start, na.rm = TRUE),
        End = max(End, na.rm = TRUE),
        Duration_h = .safe_sum(Duration_h),
        Duration_m = .safe_sum(Duration_m),
        Magnitude = .safe_sum(Magnitude),
        rate = {
          mag <- .safe_sum(Magnitude)
          dur <- .safe_sum(Duration_h)

          if (is.na(mag) || is.na(dur) || dur == 0) {
            NA_real_
          } else {
            mag * 1000 / dur
          }
        },
        .groups = "drop"
      )
  }

  .truncate_sc <- function(obj, DOY = NULL, Year = NULL) {
    if (is.null(DOY) && is.null(Year)) {
      td <- obj$SC_phase %>%
        dplyr::rename(TIME = 1, dm = 2, Phases = 3)
    } else {
      if (is.null(DOY) || is.null(Year) || length(DOY) != 2) {
        stop("Provide both DOY = c(start, end) and Year, or neither.")
      }

      td <- dendro.truncate(obj$SC_phase, CalYear = Year, DOY = DOY) %>%
        dplyr::rename(TIME = 1, dm = 2, Phases = 3)
    }

    td$TIME <- .ensure_posixct(td$TIME)

    if (any(is.na(td$TIME))) {
      stop("Some timestamps in SC_phase could not be parsed.")
    }

    if (nrow(td) == 0) {
      stop("No data available for the requested plotting period.")
    }

    cyc <- obj$SC_cycle
    cyc$Start <- .ensure_posixct(cyc$Start)
    cyc$End <- .ensure_posixct(cyc$End)

    cyc <- cyc %>%
      dplyr::filter(
        End >= min(td$TIME, na.rm = TRUE),
        Start <= max(td$TIME, na.rm = TRUE)
      )

    list(td = td, cyc = cyc)
  }

  .aggregate_points <- function(td, temporal) {
    if (temporal == "raw") {
      td$Phase <- .phase_factor(td$Phases)
      return(td)
    }

    unit <- .time_unit(temporal)

    td %>%
      dplyr::filter(!is.na(Phases)) %>%
      dplyr::mutate(period = lubridate::floor_date(TIME, unit = unit)) %>%
      dplyr::group_by(period) %>%
      dplyr::summarise(
        dm = mean(dm, na.rm = TRUE),
        Phases = {
          ptab <- table(Phases)
          as.integer(names(ptab)[which.max(ptab)])
        },
        .groups = "drop"
      ) %>%
      dplyr::rename(TIME = period) %>%
      dplyr::mutate(Phase = .phase_factor(Phases))
  }

  out <- .truncate_sc(x, DOY = DOY, Year = Year)

  td_raw <- out$td
  cyc_raw <- out$cyc

  cyc <- .aggregate_cycles(cyc_raw, temporal = temporal)
  td <- .aggregate_points(td_raw, temporal = temporal)

  if (type %in% c("points", "ribbon", "transition")) {
    if (temporal == "raw") {
      td_plot <- td_raw %>%
        dplyr::mutate(Phase = .phase_factor(Phases))
    } else {
      td_plot <- td
    }
  }

  if (type == "points") {
    xv <- .x_values(td_plot$TIME, temporal, x_axis)
    td_plot$X <- xv$x
    td_plot <- dplyr::filter(td_plot, !is.na(X))

    p <- suppressWarnings(
      ggplot2::ggplot(td_plot, ggplot2::aes(x = X, y = dm)) +
        ggplot2::geom_line(color = "grey50") +
        ggplot2::geom_point(ggplot2::aes(color = Phase), size = 2.4) +
        ggplot2::scale_color_manual(
          values = cols,
          na.value = "grey70",
          name = "Phases",
          limits = phNames
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme_bw() +
        ggplot2::labs(
          title = "Stem-cycle phases over time",
          x = xv$xlab,
          y = if (temporal == "raw") {
            "Stem size variation"
          } else {
            paste0("Stem size variation (", temporal, ")")
          }
        ) +
        ggplot2::theme(
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 12)
        )
    )

    print(p)
    return(invisible(p))
  }

  if (type == "ribbon") {
    td_r <- td_plot %>%
      dplyr::filter(!is.na(Phases))

    xv <- .x_values(td_r$TIME, temporal, x_axis)

    td_r$X <- xv$x
    td_r <- dplyr::filter(td_r, !is.na(X))

    rr <- rle(td_r$Phases)
    idx_end <- cumsum(rr$lengths)
    idx_start <- c(1, head(idx_end + 1, -1))

    ribbons <- tibble::tibble(
      xmin = td_r$X[idx_start],
      xmax = td_r$X[idx_end],
      Phase = factor(rr$values, levels = c(1, 2, 3), labels = phNames)
    )

    p <- ggplot2::ggplot(td_r, ggplot2::aes(x = X, y = dm)) +
      ggplot2::geom_rect(
        data = ribbons,
        ggplot2::aes(
          xmin = xmin,
          xmax = xmax,
          ymin = -Inf,
          ymax = Inf,
          fill = Phase
        ),
        inherit.aes = FALSE,
        alpha = 0.18
      ) +
      ggplot2::geom_line(color = "black", linewidth = 0.45) +
      ggplot2::scale_fill_manual(
        values = cols,
        name = "Phases",
        limits = phNames
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme_bw() +
      ggplot2::labs(
        title = "Stem-cycle phases as ribbons",
        x = xv$xlab,
        y = if (temporal == "raw") {
          "Stem size variation"
        } else {
          paste0("Stem size variation (", temporal, ")")
        }
      ) +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "transition") {
    td_tr <- td_plot %>%
      dplyr::filter(!is.na(Phases)) %>%
      dplyr::arrange(TIME)

    xv <- .x_values(td_tr$TIME, temporal, x_axis)

    td_tr$X <- xv$x
    td_tr <- dplyr::filter(td_tr, !is.na(X))

    trans_idx <- which(diff(td_tr$Phases) != 0) + 1L

    transitions <- if (length(trans_idx) > 0) {
      td_tr[trans_idx, , drop = FALSE]
    } else {
      td_tr[0, , drop = FALSE]
    }

    p <- ggplot2::ggplot(td_tr, ggplot2::aes(x = X, y = dm)) +
      ggplot2::geom_line(color = "grey35", linewidth = 0.5) +
      ggplot2::geom_point(ggplot2::aes(color = Phase), size = 1.8) +
      ggplot2::geom_vline(
        data = transitions,
        ggplot2::aes(xintercept = X, color = Phase),
        linetype = transition_linetype,
        alpha = transition_alpha,
        show.legend = FALSE
      ) +
      ggplot2::geom_point(
        data = transitions,
        ggplot2::aes(x = X, y = dm, color = Phase),
        size = 3,
        inherit.aes = FALSE
      ) +
      ggplot2::scale_color_manual(
        values = cols,
        na.value = "grey70",
        name = "New phase",
        limits = phNames
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme_bw() +
      ggplot2::labs(
        title = "Stem-cycle transitions",
        x = xv$xlab,
        y = if (temporal == "raw") {
          "Stem size variation"
        } else {
          paste0("Stem size variation (", temporal, ")")
        }
      ) +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "balance") {
    unit <- if (temporal == "monthly") {
      "month"
    } else {
      "day"
    }

    if (temporal == "raw") {
      warning("For type = 'balance', raw data are summarised to daily balance.")
      unit <- "day"
    }

    if (balance_mode == "duration") {
      res_h <- .resolution_hours(td_raw$TIME)

      bal <- td_raw %>%
        dplyr::filter(!is.na(Phases)) %>%
        dplyr::mutate(
          period = lubridate::floor_date(TIME, unit = unit),
          Phase = .phase_factor(Phases)
        ) %>%
        dplyr::group_by(period, Phase) %>%
        dplyr::summarise(
          hours = dplyr::n() * res_h,
          .groups = "drop"
        )

      if (x_axis == "doy" && unit == "day") {
        bal$X <- lubridate::yday(bal$period)
        x_lab <- "Day of year"
      } else {
        if (x_axis == "doy" && unit == "month") {
          warning("x_axis = 'doy' is not meaningful for monthly balance; using time instead.")
        }

        bal$X <- bal$period
        x_lab <- if (unit == "month") {
          "Month"
        } else {
          "Date"
        }
      }

      bal <- dplyr::filter(bal, !is.na(X))

      p <- ggplot2::ggplot(
        bal,
        ggplot2::aes(x = X, y = hours, fill = Phase)
      ) +
        ggplot2::geom_col() +
        ggplot2::scale_fill_manual(
          values = cols,
          name = "Phases",
          limits = phNames
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme_bw() +
        ggplot2::labs(
          title = "Phase balance",
          x = x_lab,
          y = if (unit == "month") {
            "Hours per month"
          } else {
            "Hours per day"
          }
        ) +
        ggplot2::theme(
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 12)
        )

      print(p)
      return(invisible(p))
    }

    if (balance_mode == "time_of_day") {
      if (unit == "month") {
        warning(
          "balance_mode = 'time_of_day' is designed for daily phase timing. ",
          "Using daily phase timing even though temporal = 'monthly'."
        )
        unit <- "day"
      }

      seg <- .build_phase_segments(td_raw)

      if (x_axis == "doy") {
        seg$X <- lubridate::yday(seg$period)
        seg$xmin <- seg$X - 0.45
        seg$xmax <- seg$X + 0.45
        x_lab <- "Day of year"

        p <- ggplot2::ggplot(
          seg,
          ggplot2::aes(
            xmin = xmin,
            xmax = xmax,
            ymin = ymin,
            ymax = ymax,
            fill = Phase
          )
        ) +
          ggplot2::geom_rect(color = NA)
      } else {
        seg$X <- as.Date(seg$period)
        seg$xmin <- seg$X - 0.45
        seg$xmax <- seg$X + 0.45
        x_lab <- "Date"

        p <- ggplot2::ggplot(
          seg,
          ggplot2::aes(
            xmin = xmin,
            xmax = xmax,
            ymin = ymin,
            ymax = ymax,
            fill = Phase
          )
        ) +
          ggplot2::geom_rect(color = NA) +
          ggplot2::scale_x_date()
      }

      p <- p +
        ggplot2::scale_fill_manual(
          values = cols,
          name = "Phases",
          limits = phNames
        ) +
        ggplot2::scale_y_continuous(
          limits = c(0, 24),
          breaks = seq(0, 24, by = 4),
          expand = c(0, 0)
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme_bw() +
        ggplot2::labs(
          title = "Daily phase timing",
          x = x_lab,
          y = "Hour of day"
        ) +
        ggplot2::theme(
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 12),
          axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5)
        )

      print(p)
      return(invisible(p))
    }
  }

  if (type == "increment") {
    inc <- cyc %>%
      dplyr::filter(Phases == 3L) %>%
      dplyr::arrange(Start)

    if (nrow(inc) == 0) {
      stop("No increment phases available in the selected period.")
    }

    if (temporal == "daily") {
      inc <- inc %>%
        dplyr::mutate(period = lubridate::floor_date(Start, "day")) %>%
        dplyr::group_by(period) %>%
        dplyr::summarise(
          Magnitude = sum(Magnitude, na.rm = TRUE),
          .groups = "drop"
        ) %>%
        dplyr::mutate(
          xval = period,
          cum_increment = cumsum(Magnitude)
        )
    } else if (temporal == "monthly") {
      inc <- inc %>%
        dplyr::mutate(period = lubridate::floor_date(Start, "month")) %>%
        dplyr::group_by(period) %>%
        dplyr::summarise(
          Magnitude = sum(Magnitude, na.rm = TRUE),
          .groups = "drop"
        ) %>%
        dplyr::mutate(
          xval = period,
          cum_increment = cumsum(Magnitude)
        )
    } else {
      inc <- inc %>%
        dplyr::mutate(
          xval = Start,
          cum_increment = cumsum(Magnitude)
        )
    }

    xv <- .x_values(inc$xval, temporal, x_axis)

    inc$X <- xv$x
    inc <- dplyr::filter(inc, !is.na(X))

    p <- ggplot2::ggplot(inc, ggplot2::aes(x = X, y = cum_increment)) +
      ggplot2::geom_line(color = cols[3], linewidth = 0.7) +
      ggplot2::geom_point(color = cols[3], size = 2) +
      ggplot2::theme_minimal() +
      ggplot2::theme_bw() +
      ggplot2::labs(
        title = "Cumulative increment",
        x = xv$xlab,
        y = "Cumulative increment"
      ) +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "frequency") {
    unit <- if (temporal == "monthly") {
      "month"
    } else {
      "day"
    }

    if (temporal == "raw") {
      warning("For type = 'frequency', raw data are summarised to daily event counts.")
    }

    freq <- cyc %>%
      dplyr::filter(!is.na(Phases)) %>%
      dplyr::mutate(
        period = lubridate::floor_date(Start, unit = unit),
        Phase = .phase_factor(Phases)
      ) %>%
      dplyr::group_by(period, Phase) %>%
      dplyr::summarise(n = dplyr::n(), .groups = "drop")

    if (x_axis == "doy" && unit == "day") {
      freq$X <- lubridate::yday(freq$period)
      x_lab <- "Day of year"
    } else {
      if (x_axis == "doy" && unit == "month") {
        warning("x_axis = 'doy' is not meaningful for monthly frequency; using time instead.")
      }

      freq$X <- freq$period
      x_lab <- if (unit == "month") {
        "Month"
      } else {
        "Date"
      }
    }

    freq <- dplyr::filter(freq, !is.na(X))

    p <- ggplot2::ggplot(freq, ggplot2::aes(x = X, y = n, fill = Phase)) +
      ggplot2::geom_col() +
      ggplot2::scale_fill_manual(
        values = cols,
        name = "Phases",
        limits = phNames
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme_bw() +
      ggplot2::labs(
        title = "Frequency of stem-cycle events",
        x = x_lab,
        y = if (unit == "month") {
          "Number of events per month"
        } else {
          "Number of events per day"
        }
      ) +
      ggplot2::theme(
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 12)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "heatmap") {
    unit <- if (temporal == "monthly") {
      "month"
    } else {
      "day"
    }

    if (temporal == "raw") {
      warning("For type = 'heatmap', raw data are summarised to daily dominant phase by hour.")
    }

    hm <- td_raw %>%
      dplyr::filter(!is.na(Phases)) %>%
      dplyr::mutate(
        period = if (unit == "month") {
          format(lubridate::floor_date(TIME, unit = "month"), "%Y-%m")
        } else {
          as.Date(TIME)
        },
        hour = factor(lubridate::hour(TIME), levels = 23:0),
        Phase = .phase_factor(Phases)
      ) %>%
      dplyr::group_by(period, hour, Phase) %>%
      dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
      dplyr::group_by(period, hour) %>%
      dplyr::slice_max(order_by = n, n = 1, with_ties = FALSE) %>%
      dplyr::ungroup()

    if (unit == "day" && x_axis == "doy") {
      hm$X <- lubridate::yday(hm$period)
      x_lab <- "Day of year"
    } else {
      if (unit == "month" && x_axis == "doy") {
        warning("x_axis = 'doy' is not meaningful for monthly heatmap; using time instead.")
      }

      hm$X <- if (unit == "month") {
        format(hm$period, "%Y-%m")
      } else {
        format(hm$period, "%Y-%m-%d")
      }

      x_lab <- if (unit == "month") {
        "Month"
      } else {
        "Date"
      }
    }

    p <- ggplot2::ggplot(hm, ggplot2::aes(x = X, y = hour, fill = Phase)) +
      ggplot2::geom_tile(color = "white", linewidth = 0.15) +
      ggplot2::scale_fill_manual(
        values = cols,
        name = "Phases",
        limits = phNames
      ) +
      ggplot2::theme_minimal() +
      ggplot2::theme_bw() +
      ggplot2::labs(
        title = "Dominant phase heatmap",
        x = x_lab,
        y = "Hour of day"
      ) +
      ggplot2::theme(
        axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5),
        axis.title = ggplot2::element_text(size = 14),
        axis.text = ggplot2::element_text(size = 11)
      )

    print(p)
    return(invisible(p))
  }

  if (type == "boxplot") {
    cyc_box <- cyc %>%
      dplyr::mutate(Phase = .phase_factor(Phases))

    if (phase != "all") {
      cyc_box <- cyc_box %>%
        dplyr::filter(Phase == phase)
    }

    use_continuous_time <- FALSE

    if (temporal == "daily") {
      if (x_axis == "doy") {
        full_levels <- as.character(1:366)

        cyc_box$group_label <- factor(
          as.character(lubridate::yday(cyc_box$Start)),
          levels = full_levels
        )

        x_lab <- "Day of year"
      } else {
        use_continuous_time <- TRUE

        full_dates <- seq(
          as.Date(min(td_raw$TIME, na.rm = TRUE)),
          as.Date(max(td_raw$TIME, na.rm = TRUE)),
          by = "day"
        )

        cyc_box$group_date <- as.Date(cyc_box$Start)

        x_lab <- "Date"
        x_limits <- range(full_dates)
        box_width <- 0.8
        point_jitter_width <- 0.2
        date_breaks <- if (length(full_dates) > 120) "3 months" else "1 month"
        date_labels <- if (length(unique(format(full_dates, "%Y"))) > 1) "%Y-%m" else "%m-%d"
      }
    } else {
      if (x_axis == "doy") {
        warning("x_axis = 'doy' is not meaningful for monthly boxplots; using month instead.")
      }

      use_continuous_time <- TRUE

      full_dates <- seq(
        as.Date(lubridate::floor_date(min(td_raw$TIME, na.rm = TRUE), "month")),
        as.Date(lubridate::floor_date(max(td_raw$TIME, na.rm = TRUE), "month")),
        by = "month"
      )

      cyc_box$group_date <- as.Date(lubridate::floor_date(cyc_box$Start, "month"))

      x_lab <- "Date"
      x_limits <- range(full_dates)
      box_width <- 20
      point_jitter_width <- 5
      date_breaks <- if (length(full_dates) > 24) "6 months" else "1 month"
      date_labels <- if (length(unique(format(full_dates, "%Y"))) > 1) "%Y-%m" else "%b"
    }

    cyc_box$value <- cyc_box[[stat]]

    if (use_continuous_time) {
      cyc_box <- cyc_box %>%
        dplyr::filter(!is.na(value), !is.na(group_date))
    } else {
      cyc_box <- cyc_box %>%
        dplyr::filter(!is.na(value), !is.na(group_label))
    }

    if (nrow(cyc_box) == 0) {
      stop("No non-missing values available for the selected statistic/phase.")
    }

    if (use_continuous_time) {
      count_key <- format(cyc_box$group_date, "%Y-%m-%d")

      if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
        count_key <- paste(count_key, cyc_box$Phase, sep = "_")
      }
    } else {
      count_key <- as.character(cyc_box$group_label)

      if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
        count_key <- paste(count_key, cyc_box$Phase, sep = "_")
      }
    }

    counts <- table(count_key)
    singleton_groups <- all(counts <= 1)

    if (singleton_groups) {
      warning("Each group has only one value; boxplots collapse to lines.")
    }

    if (singleton_groups && isTRUE(singleton_as_points)) {
      if (use_continuous_time) {
        p <- ggplot2::ggplot(
          cyc_box,
          ggplot2::aes(x = group_date, y = value, color = Phase)
        ) +
          ggplot2::geom_point(
            position = ggplot2::position_jitter(width = point_jitter_width, height = 0),
            size = 2,
            alpha = 0.8
          ) +
          ggplot2::scale_color_manual(
            values = cols,
            name = "Phases",
            limits = phNames
          ) +
          ggplot2::scale_x_date(
            limits = x_limits,
            date_breaks = date_breaks,
            date_labels = date_labels
          ) +
          ggplot2::theme_minimal() +
          ggplot2::theme_bw() +
          ggplot2::labs(
            title = paste("Values of", stat),
            x = x_lab,
            y = stat
          ) +
          ggplot2::theme(
            axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5),
            axis.title = ggplot2::element_text(size = 14),
            axis.text = ggplot2::element_text(size = 11)
          )
      } else {
        p <- ggplot2::ggplot(
          cyc_box,
          ggplot2::aes(x = group_label, y = value, color = Phase)
        ) +
          ggplot2::geom_point(
            position = ggplot2::position_jitter(width = 0.12, height = 0),
            size = 2,
            alpha = 0.8
          ) +
          ggplot2::scale_color_manual(
            values = cols,
            name = "Phases",
            limits = phNames
          ) +
          ggplot2::theme_minimal() +
          ggplot2::theme_bw() +
          ggplot2::labs(
            title = paste("Values of", stat),
            x = x_lab,
            y = stat
          ) +
          ggplot2::theme(
            axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5),
            axis.title = ggplot2::element_text(size = 14),
            axis.text = ggplot2::element_text(size = 11)
          ) +
          ggplot2::scale_x_discrete(drop = FALSE)
      }

      if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
        p <- p + ggplot2::facet_wrap(~ Phase, scales = "free_y")
      }

      print(p)
      return(invisible(p))
    }

    if (use_continuous_time) {
      p <- ggplot2::ggplot(
        cyc_box,
        ggplot2::aes(
          x = group_date,
          y = value,
          fill = Phase,
          group = interaction(group_date, Phase, drop = TRUE)
        )
      ) +
        ggplot2::geom_boxplot(width = box_width) +
        ggplot2::scale_fill_manual(
          values = cols,
          name = "Phases",
          limits = phNames
        ) +
        ggplot2::scale_x_date(
          limits = x_limits,
          date_breaks = date_breaks,
          date_labels = date_labels
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme_bw() +
        ggplot2::labs(
          title = paste("Boxplot of", stat),
          x = x_lab,
          y = stat
        ) +
        ggplot2::theme(
          axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5),
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 11)
        )
    } else {
      p <- ggplot2::ggplot(
        cyc_box,
        ggplot2::aes(x = group_label, y = value, fill = Phase)
      ) +
        ggplot2::geom_boxplot() +
        ggplot2::scale_fill_manual(
          values = cols,
          name = "Phases",
          limits = phNames
        ) +
        ggplot2::theme_minimal() +
        ggplot2::theme_bw() +
        ggplot2::labs(
          title = paste("Boxplot of", stat),
          x = x_lab,
          y = stat
        ) +
        ggplot2::theme(
          axis.text.x = ggplot2::element_text(angle = 90, hjust = 1, vjust = 0.5),
          axis.title = ggplot2::element_text(size = 14),
          axis.text = ggplot2::element_text(size = 11)
        ) +
        ggplot2::scale_x_discrete(drop = FALSE)
    }

    if (phase == "all" && length(unique(cyc_box$Phase)) > 1) {
      p <- p + ggplot2::facet_wrap(~ Phase, scales = "free_y")
    }

    print(p)
    return(invisible(p))
  }

  stop("Unknown plot type.")
}