R/plot.dm_growth_fit.R

In dendRoAnalyst: A Tool for Processing and Analyzing Dendrometer Data

#' Plot dendrometer growth-fit results
#'
#' @description
#' Creates \pkg{ggplot2}-based plots for objects returned by
#' [dm.growth.fit()] and [dm.growth.fit.double()].
#'
#' The plotting method supports multiple views of fitted dendrometer growth
#' curves, including observed versus fitted trajectories, residuals, timing
#' summaries, overlays of multiple curves, and distributions of fitted model
#' parameters.
#'
#' The x-axis can be displayed as vegetation-season day, calendar day-of-year
#' (DOY), or actual date, depending on the selected \code{x_axis} argument and
#' the type of plot.
#'
#' When \code{x_axis = "date"} and faceting separates individual years or
#' series-year combinations, each facet receives its own x-axis range. This
#' avoids plotting one yearly fitted curve inside the full dendrometer time
#' window when several seasons are present in the input object.
#'
#' @param x An object of class \code{"dm_growth_fit"} returned by
#'   [dm.growth.fit()] or [dm.growth.fit.double()].
#' @param type Character string specifying the plot type. One of
#'   \code{"fit"}, \code{"season"}, \code{"residuals"}, \code{"timing"},
#'   \code{"overlay"}, or \code{"parameters"}.
#' @param series Optional filter selecting one or more dendrometer series to
#'   plot. May be a character vector of series names. Default is \code{NULL},
#'   meaning all available series are used.
#' @param fit_id Optional filter selecting one or more fit identifiers to plot.
#'   May be a character or numeric vector. Numeric values are coerced internally
#'   to character. Default is \code{NULL}, meaning all fits are used.
#' @param facet_by Character string controlling faceting layout. One of
#'   \code{"default"}, \code{"tree"}, \code{"year"}, or \code{"none"}.
#' @param ncol Optional integer giving the number of columns in faceted plots.
#'   Passed to [ggplot2::facet_wrap()].
#' @param normalize Logical. If \code{TRUE}, observed and fitted values are
#'   divided by the maximum observed or fitted value within each curve.
#' @param x_axis Character string controlling the x-axis representation. One of
#'   \code{"default"}, \code{"season_day"}, \code{"doy"}, or \code{"date"}.
#' @param observed_source Character string controlling which observed daily
#'   series is plotted against the fitted curve. One of \code{"processed"} or
#'   \code{"original_daily"}.
#' @param show_observed Logical. If \code{TRUE}, observed values are shown in
#'   plot types where observed data are relevant.
#' @param show_fitted Logical. If \code{TRUE}, fitted values are shown in plot
#'   types where fitted curves are relevant.
#' @param show_timing Logical. If \code{TRUE}, timing markers are added to
#'   \code{"fit"} and \code{"season"} plots whenever timing information is
#'   available in \code{x$fit_statistics}.
#' @param point_alpha Numeric alpha level used for observed points.
#' @param line_width Numeric line width used for fitted curves.
#' @param legend_position Character string specifying legend position, passed to
#'   [ggplot2::theme()].
#' @param ... Further arguments passed to or from other methods.
#'
#' @details
#' The plotting method returns a \code{ggplot} object. The returned plot can be
#' further modified using normal \pkg{ggplot2} syntax.
#'
#' Timing markers are taken from the \code{fit_statistics} table inside the
#' \code{"dm_growth_fit"} object:
#' \itemize{
#'   \item \code{growth_start_*} and \code{growth_end_*} represent growing-season
#'     timing based on cumulative fitted growth.
#'   \item \code{rate_start_*} and \code{rate_end_*} represent active-growth
#'     timing based on the fitted growth-rate curve.
#'   \item For double-growth fits, pulse-specific timing such as
#'     \code{pulse1_start_*}, \code{pulse2_start_*}, and \code{separator_*} is
#'     used in \code{type = "timing"} when available.
#' }
#'
#' For yearly fits in southern hemisphere or cross-year custom seasons,
#' calendar timing variables such as \code{growth_start_day} and
#' \code{rate_start_day} are interpreted as true calendar DOY, while
#' \code{*_season_day} variables represent day counts relative to vegetation
#' season start.
#'
#' In pooled fits, date-based and calendar-DOY timing fields may be unavailable,
#' because pooled fits are not anchored to a single season start date.
#'
#' @return A \code{ggplot} object.
#'
#' @seealso
#' [dm.growth.fit()], [dm.growth.fit.double()], [summary.dm_growth_fit()]
#'
#' @examples
#' \donttest{
#' # fit <- dm.growth.fit(...)
#' # plot(fit, type = "fit")
#' # plot(fit, type = "fit", facet_by = "year", x_axis = "date")
#' # plot(fit, type = "fit", facet_by = "year", x_axis = "season_day")
#' # plot(fit, type = "residuals", facet_by = "year")
#' # plot(fit, type = "timing")
#' # plot(fit, type = "parameters")
#' }
#'
#' @method plot dm_growth_fit
#' @importFrom ggplot2 ggplot aes theme_bw theme labs facet_wrap geom_point
#'   geom_line geom_hline geom_vline geom_text geom_segment geom_boxplot
#'   geom_jitter element_text
#' @importFrom dplyr %>% group_by mutate if_else ungroup filter select all_of
#'   left_join transmute bind_rows
#' @importFrom tidyr pivot_longer
#' @importFrom tibble as_tibble
#' @importFrom stats as.formula na.omit
#' @importFrom rlang .data
#' @export
plot.dm_growth_fit <- function(x,
                               type = c("fit", "season", "residuals", "timing", "overlay", "parameters"),
                               series = NULL,
                               fit_id = NULL,
                               facet_by = c("default", "tree", "year", "none"),
                               ncol = NULL,
                               normalize = FALSE,
                               x_axis = c("default", "season_day", "doy", "date"),
                               observed_source = c("processed", "original_daily"),
                               show_observed = TRUE,
                               show_fitted = TRUE,
                               show_timing = TRUE,
                               point_alpha = 0.7,
                               line_width = 0.8,
                               legend_position = "right",
                               ...) {

  type <- match.arg(type)
  facet_by <- match.arg(facet_by)
  x_axis <- match.arg(x_axis)
  observed_source <- match.arg(observed_source)

  if (!requireNamespace("ggplot2", quietly = TRUE)) {
    stop("Package 'ggplot2' is required for plot.dm_growth_fit().")
  }

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

  if (type %in% c("fit", "season", "residuals", "overlay")) {
    dat <- dmgp_tidy_fit_data(x, observed_source = observed_source)

    if (!is.null(series)) {
      dat <- dat[dat$series %in% series, , drop = FALSE]
    }

    if (!is.null(fit_id)) {
      dat <- dat[dat$fit_id %in% as.character(fit_id), , drop = FALSE]
    }

    if (nrow(dat) == 0) {
      stop("No rows remain after filtering by series and/or fit_id.")
    }

    cfg <- dmgp_prepare_curve_groups(dat, facet_by = facet_by)
    dat <- cfg$data

    x_cfg <- dmgp_resolve_x_axis(type = type, x_axis = x_axis)
    x_var <- x_cfg$var
    x_lab <- x_cfg$label

    facet_scales <- dmgp_resolve_facet_scales(
      x_var = x_var,
      facet_by = facet_by
    )

    if (!x_var %in% names(dat)) {
      stop("Column '", x_var, "' not found in plotting data.")
    }

    if (isTRUE(normalize)) {
      dat <- dat %>%
        dplyr::group_by(.data$group_id) %>%
        dplyr::mutate(
          .norm = suppressWarnings(max(c(.data$observed, .data$fitted), na.rm = TRUE)),
          observed_plot = dplyr::if_else(
            is.finite(.data$.norm) & .data$.norm > 0,
            .data$observed / .data$.norm,
            NA_real_
          ),
          fitted_plot = dplyr::if_else(
            is.finite(.data$.norm) & .data$.norm > 0,
            .data$fitted / .data$.norm,
            NA_real_
          )
        ) %>%
        dplyr::ungroup()

      y_lab <- "Normalized growth"

    } else {
      dat$observed_plot <- dat$observed
      dat$fitted_plot <- dat$fitted
      y_lab <- "Growth"
    }

    plot_title_obs <- if (observed_source == "original_daily") {
      "Original daily and fitted growth"
    } else {
      "Observed and fitted growth"
    }

    if (type == "fit") {
      p <- ggplot2::ggplot() +
        ggplot2::theme_bw() +
        ggplot2::theme(
          legend.position = if (isTRUE(cfg$use_colour)) legend_position else "none"
        )

      if (isTRUE(show_observed)) {
        p <- p +
          ggplot2::geom_point(
            data = dat[is.finite(dat$observed_plot), , drop = FALSE],
            mapping = ggplot2::aes(
              x = .data[[x_var]],
              y = .data$observed_plot,
              colour = .data$colour_group,
              group = .data$group_id
            ),
            alpha = point_alpha,
            size = 1.4
          )
      }

      if (isTRUE(show_fitted)) {
        p <- p +
          ggplot2::geom_line(
            data = dat[is.finite(dat$fitted_plot), , drop = FALSE],
            mapping = ggplot2::aes(
              x = .data[[x_var]],
              y = .data$fitted_plot,
              colour = .data$colour_group,
              group = .data$group_id
            ),
            linewidth = line_width
          )
      }

      p <- p +
        ggplot2::labs(
          x = x_lab,
          y = y_lab,
          colour = cfg$legend_title,
          title = plot_title_obs
        )

      if (isTRUE(cfg$facet_enabled)) {
        p <- p +
          ggplot2::facet_wrap(
            stats::as.formula("~ facet_var"),
            scales = facet_scales,
            ncol = ncol
          )
      }

      if (isTRUE(show_timing) && !is.null(x$fit_statistics)) {
        timing_marks <- dmgp_prepare_boundary_marks(
          fit_statistics = x$fit_statistics,
          facet_by = facet_by,
          x_mode = x_cfg$timing_mode,
          series = series,
          fit_id = fit_id
        )

        if (nrow(timing_marks$cumulative) > 0) {
          p <- p +
            ggplot2::geom_vline(
              data = timing_marks$cumulative,
              mapping = ggplot2::aes(
                xintercept = .data$x,
                linetype = .data$boundary,
                colour = .data$curve_group
              ),
              alpha = 0.5
            )
        }

        if (nrow(timing_marks$rate) > 0) {
          p <- p +
            ggplot2::geom_text(
              data = timing_marks$rate,
              mapping = ggplot2::aes(
                x = .data$x,
                y = .data$y,
                label = .data$triangle,
                colour = .data$curve_group
              ),
              inherit.aes = FALSE,
              show.legend = FALSE,
              size = 4
            )
        }
      }

      return(p)
    }

    if (type == "season") {
      p <- ggplot2::ggplot() +
        ggplot2::theme_bw() +
        ggplot2::theme(
          legend.position = if (isTRUE(cfg$use_colour)) legend_position else "none"
        )

      if (isTRUE(show_observed)) {
        p <- p +
          ggplot2::geom_point(
            data = dat[is.finite(dat$observed_plot), , drop = FALSE],
            mapping = ggplot2::aes(
              x = .data[[x_var]],
              y = .data$observed_plot,
              colour = .data$colour_group,
              group = .data$group_id
            ),
            alpha = point_alpha,
            size = 1.4
          )
      }

      if (isTRUE(show_fitted)) {
        p <- p +
          ggplot2::geom_line(
            data = dat[is.finite(dat$fitted_plot), , drop = FALSE],
            mapping = ggplot2::aes(
              x = .data[[x_var]],
              y = .data$fitted_plot,
              colour = .data$colour_group,
              group = .data$group_id
            ),
            linewidth = line_width
          )
      }

      p <- p +
        ggplot2::labs(
          x = x_lab,
          y = y_lab,
          colour = cfg$legend_title,
          title = plot_title_obs
        )

      if (isTRUE(cfg$facet_enabled)) {
        p <- p +
          ggplot2::facet_wrap(
            stats::as.formula("~ facet_var"),
            scales = facet_scales,
            ncol = ncol
          )
      }

      if (isTRUE(show_timing) && !is.null(x$fit_statistics)) {
        timing_marks <- dmgp_prepare_boundary_marks(
          fit_statistics = x$fit_statistics,
          facet_by = facet_by,
          x_mode = x_cfg$timing_mode,
          series = series,
          fit_id = fit_id
        )

        if (nrow(timing_marks$cumulative) > 0) {
          p <- p +
            ggplot2::geom_vline(
              data = timing_marks$cumulative,
              mapping = ggplot2::aes(
                xintercept = .data$x,
                linetype = .data$boundary,
                colour = .data$curve_group
              ),
              alpha = 0.5
            )
        }

        if (nrow(timing_marks$rate) > 0) {
          p <- p +
            ggplot2::geom_text(
              data = timing_marks$rate,
              mapping = ggplot2::aes(
                x = .data$x,
                y = .data$y,
                label = .data$triangle,
                colour = .data$curve_group
              ),
              inherit.aes = FALSE,
              show.legend = FALSE,
              size = 4
            )
        }
      }

      return(p)
    }

    if (type == "residuals") {
      dat$residual <- dat$observed - dat$fitted
      dat <- dat[is.finite(dat$residual), , drop = FALSE]

      if (nrow(dat) == 0) {
        stop("No finite residuals available to plot.")
      }

      p <- ggplot2::ggplot(
        dat,
        ggplot2::aes(
          x = .data[[x_var]],
          y = .data$residual,
          colour = .data$colour_group,
          group = .data$group_id
        )
      ) +
        ggplot2::geom_hline(yintercept = 0, linetype = 2) +
        ggplot2::geom_point(alpha = point_alpha, size = 1.4) +
        ggplot2::theme_bw() +
        ggplot2::theme(
          legend.position = if (isTRUE(cfg$use_colour)) legend_position else "none"
        ) +
        ggplot2::labs(
          x = x_lab,
          y = "Observed - fitted",
          colour = cfg$legend_title,
          title = if (observed_source == "original_daily") {
            "Residuals of original daily data against fitted curve"
          } else {
            "Residuals"
          }
        )

      if (isTRUE(cfg$facet_enabled)) {
        p <- p +
          ggplot2::facet_wrap(
            stats::as.formula("~ facet_var"),
            scales = facet_scales,
            ncol = ncol
          )
      }

      return(p)
    }

    if (type == "overlay") {
      p <- ggplot2::ggplot() +
        ggplot2::theme_bw() +
        ggplot2::theme(legend.position = legend_position)

      if (isTRUE(show_observed)) {
        p <- p +
          ggplot2::geom_point(
            data = dat[is.finite(dat$observed_plot), , drop = FALSE],
            mapping = ggplot2::aes(
              x = .data[[x_var]],
              y = .data$observed_plot,
              colour = .data$colour_group,
              group = .data$group_id
            ),
            alpha = 0.35,
            size = 1.1
          )
      }

      if (isTRUE(show_fitted)) {
        p <- p +
          ggplot2::geom_line(
            data = dat[is.finite(dat$fitted_plot), , drop = FALSE],
            mapping = ggplot2::aes(
              x = .data[[x_var]],
              y = .data$fitted_plot,
              colour = .data$colour_group,
              group = .data$group_id
            ),
            linewidth = line_width
          )
      }

      p <- p +
        ggplot2::labs(
          x = x_lab,
          y = y_lab,
          colour = cfg$legend_title,
          title = if (observed_source == "original_daily") {
            "Overlay of fitted curves with original daily data"
          } else {
            "Overlay of fitted growth curves"
          }
        )

      if (facet_by == "default") {
        p <- p +
          ggplot2::facet_wrap(
            stats::as.formula("~ series"),
            scales = facet_scales,
            ncol = ncol
          )
      } else if (isTRUE(cfg$facet_enabled)) {
        p <- p +
          ggplot2::facet_wrap(
            stats::as.formula("~ facet_var"),
            scales = facet_scales,
            ncol = ncol
          )
      }

      return(p)
    }
  }

  if (type == "timing") {
    x_cfg <- dmgp_resolve_x_axis(type = "timing", x_axis = x_axis)

    tdat <- dmgp_prepare_timing_plot_data(
      fit_statistics = x$fit_statistics,
      facet_by = facet_by,
      x_mode = x_cfg$timing_mode,
      series = series,
      fit_id = fit_id
    )

    if (nrow(tdat$base_rows) == 0) {
      stop("No timing information available to plot.")
    }

    timing_facet_scales <- dmgp_resolve_timing_facet_scales(
      x_mode = x_cfg$timing_mode,
      facet_by = facet_by
    )

    p <- ggplot2::ggplot() +
      ggplot2::theme_bw() +
      ggplot2::theme(legend.position = legend_position)

    if (nrow(tdat$overall_segments) > 0) {
      p <- p +
        ggplot2::geom_segment(
          data = tdat$overall_segments,
          mapping = ggplot2::aes(
            x = .data$x,
            xend = .data$xend,
            y = .data$row_label,
            yend = .data$row_label,
            colour = .data$window
          ),
          linewidth = 1
        )
    }

    if (nrow(tdat$overall_ticks) > 0) {
      p <- p +
        ggplot2::geom_point(
          data = tdat$overall_ticks,
          mapping = ggplot2::aes(
            x = .data$x,
            y = .data$row_label,
            colour = .data$window
          ),
          shape = 124,
          size = 6,
          stroke = 1.1
        )
    }

    if (nrow(tdat$rate_points) > 0) {
      p <- p +
        ggplot2::geom_text(
          data = tdat$rate_points,
          mapping = ggplot2::aes(
            x = .data$x,
            y = .data$row_label,
            label = .data$triangle,
            colour = .data$window
          ),
          size = 4.5,
          show.legend = FALSE
        )
    }

    if (nrow(tdat$pulse_segments) > 0) {
      p <- p +
        ggplot2::geom_segment(
          data = tdat$pulse_segments,
          mapping = ggplot2::aes(
            x = .data$x,
            xend = .data$xend,
            y = .data$row_label,
            yend = .data$row_label,
            colour = .data$window
          ),
          linewidth = 0.8,
          alpha = 0.9
        )
    }

    if (nrow(tdat$separator_points) > 0) {
      p <- p +
        ggplot2::geom_point(
          data = tdat$separator_points,
          mapping = ggplot2::aes(
            x = .data$x,
            y = .data$row_label,
            colour = .data$window
          ),
          size = 2.2
        )
    }

    p <- p +
      ggplot2::labs(
        x = x_cfg$label,
        y = tdat$y_label,
        colour = "Timing element",
        title = "Timing of overall season and pulses"
      )

    if (isTRUE(tdat$facet_enabled)) {
      p <- p +
        ggplot2::facet_wrap(
          stats::as.formula("~ facet_var"),
          scales = timing_facet_scales,
          ncol = ncol
        )
    }

    return(p)
  }

  if (type == "parameters") {
    if (is.null(x$fit_parameters) || nrow(x$fit_parameters) == 0) {
      stop("No fit_parameters table found in the object.")
    }

    pars <- tibble::as_tibble(x$fit_parameters)

    if (!is.null(series)) {
      pars <- pars[pars$series %in% series, , drop = FALSE]
    }

    if (!is.null(fit_id)) {
      pars <- pars[pars$fit_id %in% as.character(fit_id), , drop = FALSE]
    }

    if (nrow(pars) == 0) {
      stop("No parameter rows remain after filtering by series and/or fit_id.")
    }

    param_candidates <- c(
      "a", "b", "k", "t0", "v",
      "a2", "b2", "k2", "t02", "v2",
      "edf", "span", "spline_df", "spar"
    )

    param_cols <- intersect(param_candidates, names(pars))

    if (length(param_cols) == 0) {
      stop("No parameter columns available to plot.")
    }

    plong <- pars %>%
      dplyr::select(.data$series, .data$fit_id, dplyr::all_of(param_cols)) %>%
      tidyr::pivot_longer(
        cols = dplyr::all_of(param_cols),
        names_to = "parameter",
        values_to = "value"
      ) %>%
      dplyr::filter(is.finite(.data$value))

    if (nrow(plong) == 0) {
      stop("No finite parameter values available to plot.")
    }

    if (facet_by == "tree") {
      plong$group_label <- plong$fit_id
      x_lab <- "Fit ID"
    } else if (facet_by == "year") {
      plong$group_label <- plong$series
      x_lab <- "Series"
    } else {
      if (length(unique(plong$series)) > 1) {
        plong$group_label <- plong$series
        x_lab <- "Series"
      } else {
        plong$group_label <- plong$fit_id
        x_lab <- "Fit ID"
      }
    }

    p <- ggplot2::ggplot(
      plong,
      ggplot2::aes(x = .data$group_label, y = .data$value)
    ) +
      ggplot2::geom_boxplot(outlier.shape = NA) +
      ggplot2::geom_jitter(width = 0.15, alpha = 0.7, size = 1.6) +
      ggplot2::facet_wrap(
        stats::as.formula("~ parameter"),
        scales = "free_y",
        ncol = ncol
      ) +
      ggplot2::theme_bw() +
      ggplot2::theme(
        legend.position = "none",
        axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)
      ) +
      ggplot2::labs(
        x = x_lab,
        y = "Parameter value",
        title = "Distribution of fitted parameters"
      )

    return(p)
  }

  stop("Unknown plot type.")
}


# helpers ----------------------------------------------------------------------

#' Convert growth-fit object to tidy plotting data
#'
#' @description
#' Internal helper used by [plot.dm_growth_fit()] to combine fitted and observed
#' daily dendrometer growth data into a long-format plotting table.
#'
#' @param x An object of class \code{"dm_growth_fit"}.
#' @param observed_source Character string. Either \code{"processed"} or
#'   \code{"original_daily"}.
#'
#' @return A tibble containing fitted values, observed values, series labels,
#'   fit identifiers, and panel identifiers.
#'
#' @keywords internal
dmgp_tidy_fit_data <- function(x,
                               observed_source = c("processed", "original_daily")) {
  observed_source <- match.arg(observed_source)

  fitd <- tibble::as_tibble(x$fitted_data)

  fit_meta_cols <- intersect(
    c(
      "TIME", "doy", "season_label", "season_start", "season_end",
      "season_day", "season_length", "any_observed"
    ),
    names(fitd)
  )

  fit_series_cols <- setdiff(names(fitd), fit_meta_cols)

  if (length(fit_series_cols) == 0) {
    stop("No fitted series columns found in fitted_data.")
  }

  if (observed_source == "processed") {
    obsd <- tibble::as_tibble(x$processed_data)
  } else {
    if (is.null(x$original_daily_data)) {
      stop(
        "This dm_growth_fit object does not contain 'original_daily_data'. ",
        "Refit with the updated fitting function."
      )
    }

    obsd <- tibble::as_tibble(x$original_daily_data)
  }

  obs_meta_cols <- intersect(
    c(
      "TIME", "doy", "season_label", "season_start", "season_end",
      "season_day", "season_length", "any_observed"
    ),
    names(obsd)
  )

  obs_series_cols <- intersect(
    fit_series_cols,
    setdiff(names(obsd), obs_meta_cols)
  )

  if (length(obs_series_cols) == 0) {
    stop("No matching observed series columns found for plotting.")
  }

  fit_long <- fitd %>%
    dplyr::select(
      dplyr::all_of(fit_meta_cols),
      dplyr::all_of(obs_series_cols)
    ) %>%
    tidyr::pivot_longer(
      cols = dplyr::all_of(obs_series_cols),
      names_to = "series",
      values_to = "fitted"
    )

  obs_long <- obsd %>%
    dplyr::select(
      dplyr::all_of(obs_meta_cols),
      dplyr::all_of(obs_series_cols)
    ) %>%
    tidyr::pivot_longer(
      cols = dplyr::all_of(obs_series_cols),
      names_to = "series",
      values_to = "observed"
    )

  join_cols <- intersect(names(fit_long), names(obs_long))
  join_cols <- setdiff(join_cols, c("observed", "fitted"))

  dat <- dplyr::left_join(fit_long, obs_long, by = join_cols)

  if (
    !is.null(x$fit_statistics) &&
    "fit_id" %in% names(x$fit_statistics) &&
    length(unique(stats::na.omit(x$fit_statistics$fit_id))) == 1 &&
    unique(stats::na.omit(x$fit_statistics$fit_id)) == "pooled"
  ) {
    dat$fit_id <- "pooled"
  } else {
    dat$fit_id <- as.character(dat$season_label)
  }

  dat$panel_id <- paste(dat$series, dat$fit_id, sep = " | ")

  dat
}


#' Prepare curve grouping and faceting information
#'
#' @description
#' Internal helper used by [plot.dm_growth_fit()] to define panel identifiers,
#' facet variables, colour groups, and legend titles.
#'
#' @param dat Tidy plotting data returned by [dmgp_tidy_fit_data()].
#' @param facet_by Character string controlling faceting layout.
#'
#' @return A list containing modified plotting data and faceting settings.
#'
#' @keywords internal
dmgp_prepare_curve_groups <- function(dat,
                                      facet_by = c("default", "tree", "year", "none")) {
  facet_by <- match.arg(facet_by)

  n_series <- length(unique(dat$series))
  n_fit <- length(unique(dat$fit_id))

  dat$group_id <- dat$panel_id

  if (facet_by == "default") {
    dat$facet_var <- dat$panel_id
    dat$colour_group <- factor("curve")

    return(list(
      data = dat,
      facet_enabled = TRUE,
      use_colour = FALSE,
      legend_title = NULL
    ))
  }

  if (facet_by == "tree") {
    dat$facet_var <- dat$series
    dat$colour_group <- factor(dat$fit_id)

    return(list(
      data = dat,
      facet_enabled = TRUE,
      use_colour = n_fit > 1,
      legend_title = "Fit ID"
    ))
  }

  if (facet_by == "year") {
    dat$facet_var <- dat$fit_id
    dat$colour_group <- factor(dat$series)

    return(list(
      data = dat,
      facet_enabled = TRUE,
      use_colour = n_series > 1,
      legend_title = "Series"
    ))
  }

  dat$facet_var <- "All"

  if (n_series > 1 && n_fit > 1) {
    dat$colour_group <- factor(dat$panel_id)
    legend_title <- "Series | Fit"
  } else if (n_series > 1) {
    dat$colour_group <- factor(dat$series)
    legend_title <- "Series"
  } else if (n_fit > 1) {
    dat$colour_group <- factor(dat$fit_id)
    legend_title <- "Fit ID"
  } else {
    dat$colour_group <- factor("curve")
    legend_title <- NULL
  }

  list(
    data = dat,
    facet_enabled = FALSE,
    use_colour = length(unique(dat$colour_group)) > 1,
    legend_title = legend_title
  )
}


#' Resolve x-axis mode for dendrometer growth-fit plots
#'
#' @description
#' Internal helper used by [plot.dm_growth_fit()] to determine which column is
#' used on the x-axis and which timing coordinate system should be used.
#'
#' @param type Plot type.
#' @param x_axis Requested x-axis representation.
#'
#' @return A list with elements \code{var}, \code{label}, and
#'   \code{timing_mode}.
#'
#' @keywords internal
dmgp_resolve_x_axis <- function(type,
                                x_axis = c("default", "season_day", "doy", "date")) {
  x_axis <- match.arg(x_axis)

  if (x_axis == "default") {
    if (type %in% c("fit", "timing")) {
      return(list(
        var = "TIME",
        label = "Date",
        timing_mode = "date"
      ))
    }

    return(list(
      var = "season_day",
      label = "Season day",
      timing_mode = "season_day"
    ))
  }

  if (x_axis == "date") {
    return(list(
      var = "TIME",
      label = "Date",
      timing_mode = "date"
    ))
  }

  if (x_axis == "doy") {
    return(list(
      var = "doy",
      label = "DOY",
      timing_mode = "doy"
    ))
  }

  list(
    var = "season_day",
    label = "Season day",
    timing_mode = "season_day"
  )
}


#' Resolve facet scales for curve plots
#'
#' @description
#' Internal helper used by [plot.dm_growth_fit()] to choose facet scale behavior.
#'
#' For date-based plots faceted by year or by individual series-fit panels, the
#' x-axis is freed so that each year or season is shown only over its own
#' temporal window. This prevents yearly fitted curves from being plotted inside
#' the full dendrometer record window when several seasons are present.
#'
#' For season-day and DOY plots, the x-axis remains shared so curves remain
#' directly comparable across facets.
#'
#' @param x_var Name of the x-axis variable.
#' @param facet_by Faceting mode.
#'
#' @return A character string accepted by [ggplot2::facet_wrap()] argument
#'   \code{scales}.
#'
#' @keywords internal
dmgp_resolve_facet_scales <- function(x_var,
                                      facet_by = c("default", "tree", "year", "none")) {
  facet_by <- match.arg(facet_by)

  if (identical(x_var, "TIME") && facet_by %in% c("default", "year")) {
    return("free")
  }

  "free_y"
}


#' Resolve facet scales for timing plots
#'
#' @description
#' Internal helper used by [plot.dm_growth_fit()] to choose facet scales for
#' timing-summary plots.
#'
#' @param x_mode Timing coordinate system. One of \code{"date"}, \code{"doy"},
#'   or \code{"season_day"}.
#' @param facet_by Faceting mode.
#'
#' @return A character string accepted by [ggplot2::facet_wrap()] argument
#'   \code{scales}.
#'
#' @keywords internal
dmgp_resolve_timing_facet_scales <- function(x_mode = c("date", "doy", "season_day"),
                                             facet_by = c("default", "tree", "year", "none")) {
  x_mode <- match.arg(x_mode)
  facet_by <- match.arg(facet_by)

  if (identical(x_mode, "date") && facet_by %in% c("default", "year", "tree")) {
    return("free")
  }

  "free_y"
}


#' Prepare timing boundary markers for fitted-curve plots
#'
#' @description
#' Internal helper used by [plot.dm_growth_fit()] to prepare vertical timing
#' markers for \code{type = "fit"} and \code{type = "season"} plots.
#'
#' @param fit_statistics Fit-statistics table from a \code{"dm_growth_fit"}
#'   object.
#' @param facet_by Faceting mode.
#' @param x_mode Timing coordinate system. One of \code{"date"}, \code{"doy"},
#'   or \code{"season_day"}.
#' @param series Optional series filter.
#' @param fit_id Optional fit identifier filter.
#'
#' @return A list with two data frames: \code{cumulative} and \code{rate}.
#'
#' @keywords internal
dmgp_prepare_boundary_marks <- function(fit_statistics,
                                        facet_by = c("default", "tree", "year", "none"),
                                        x_mode = c("date", "doy", "season_day"),
                                        series = NULL,
                                        fit_id = NULL) {
  facet_by <- match.arg(facet_by)
  x_mode <- match.arg(x_mode)

  fs <- tibble::as_tibble(fit_statistics)

  if (!is.null(series)) {
    fs <- fs[fs$series %in% series, , drop = FALSE]
  }

  if (!is.null(fit_id)) {
    fs <- fs[fs$fit_id %in% as.character(fit_id), , drop = FALSE]
  }

  if (nrow(fs) == 0) {
    return(list(
      cumulative = data.frame(),
      rate = data.frame()
    ))
  }

  fs$panel_id <- paste(fs$series, fs$fit_id, sep = " | ")

  if (facet_by == "tree") {
    fs$facet_var <- fs$series
    fs$curve_group <- as.character(fs$fit_id)
  } else if (facet_by == "year") {
    fs$facet_var <- as.character(fs$fit_id)
    fs$curve_group <- as.character(fs$series)
  } else if (facet_by == "none") {
    fs$facet_var <- "All"
    fs$curve_group <- if (length(unique(fs$panel_id)) > 1) {
      fs$panel_id
    } else {
      "curve"
    }
  } else {
    fs$facet_var <- fs$panel_id
    fs$curve_group <- "curve"
  }

  if (x_mode == "date") {
    gstart <- "growth_start_date"
    gend <- "growth_end_date"
    rstart <- "rate_start_date"
    rend <- "rate_end_date"
  } else if (x_mode == "doy") {
    gstart <- "growth_start_day"
    gend <- "growth_end_day"
    rstart <- "rate_start_day"
    rend <- "rate_end_day"
  } else {
    gstart <- "growth_start_season_day"
    gend <- "growth_end_season_day"
    rstart <- "rate_start_season_day"
    rend <- "rate_end_season_day"
  }

  cumulative <- data.frame()

  if (all(c(gstart, gend) %in% names(fs))) {
    cumulative <- dplyr::bind_rows(
      fs %>%
        dplyr::transmute(
          facet_var = .data$facet_var,
          curve_group = .data$curve_group,
          x = .data[[gstart]],
          boundary = "cumulative_start"
        ),
      fs %>%
        dplyr::transmute(
          facet_var = .data$facet_var,
          curve_group = .data$curve_group,
          x = .data[[gend]],
          boundary = "cumulative_end"
        )
    ) %>%
      dplyr::filter(!is.na(.data$x))
  }

  rate <- data.frame()

  if (all(c(rstart, rend) %in% names(fs))) {
    rate <- dplyr::bind_rows(
      fs %>%
        dplyr::transmute(
          facet_var = .data$facet_var,
          curve_group = .data$curve_group,
          x = .data[[rstart]],
          y = Inf,
          triangle = "\u25C0"
        ),
      fs %>%
        dplyr::transmute(
          facet_var = .data$facet_var,
          curve_group = .data$curve_group,
          x = .data[[rend]],
          y = Inf,
          triangle = "\u25B6"
        )
    ) %>%
      dplyr::filter(!is.na(.data$x))
  }

  list(
    cumulative = cumulative,
    rate = rate
  )
}


#' Prepare timing-summary plot data
#'
#' @description
#' Internal helper used by [plot.dm_growth_fit()] to prepare the line segments,
#' point markers, and labels used in \code{type = "timing"} plots.
#'
#' @param fit_statistics Fit-statistics table from a \code{"dm_growth_fit"}
#'   object.
#' @param facet_by Faceting mode.
#' @param x_mode Timing coordinate system. One of \code{"date"}, \code{"doy"},
#'   or \code{"season_day"}.
#' @param series Optional series filter.
#' @param fit_id Optional fit identifier filter.
#'
#' @return A list of data frames used by the timing plot.
#'
#' @keywords internal
dmgp_prepare_timing_plot_data <- function(fit_statistics,
                                          facet_by = c("default", "tree", "year", "none"),
                                          x_mode = c("date", "doy", "season_day"),
                                          series = NULL,
                                          fit_id = NULL) {
  facet_by <- match.arg(facet_by)
  x_mode <- match.arg(x_mode)

  fs <- tibble::as_tibble(fit_statistics)

  if (!is.null(series)) {
    fs <- fs[fs$series %in% series, , drop = FALSE]
  }

  if (!is.null(fit_id)) {
    fs <- fs[fs$fit_id %in% as.character(fit_id), , drop = FALSE]
  }

  if (nrow(fs) == 0) {
    return(list(
      x_mode = "season_day",
      facet_enabled = FALSE,
      y_label = "Fit",
      base_rows = fs[0, , drop = FALSE],
      overall_segments = data.frame(),
      overall_ticks = data.frame(),
      rate_points = data.frame(),
      pulse_segments = data.frame(),
      separator_points = data.frame()
    ))
  }

  if (facet_by %in% c("default", "tree")) {
    fs$facet_var <- fs$series
    fs$row_label <- as.character(fs$fit_id)
    y_label <- "Fit ID"
    facet_enabled <- length(unique(fs$facet_var)) > 1
  } else if (facet_by == "year") {
    fs$facet_var <- as.character(fs$fit_id)
    fs$row_label <- fs$series
    y_label <- "Series"
    facet_enabled <- length(unique(fs$facet_var)) > 1
  } else {
    fs$facet_var <- "All"
    fs$row_label <- paste(fs$series, fs$fit_id, sep = " | ")
    y_label <- "Series | Fit"
    facet_enabled <- FALSE
  }

  if (x_mode == "date") {
    gstart <- "growth_start_date"
    gend <- "growth_end_date"
    rstart <- "rate_start_date"
    rend <- "rate_end_date"
    p1s <- "pulse1_start_date"
    p1e <- "pulse1_end_date"
    p2s <- "pulse2_start_date"
    p2e <- "pulse2_end_date"
    sepv <- "separator_date"
  } else if (x_mode == "doy") {
    gstart <- "growth_start_day"
    gend <- "growth_end_day"
    rstart <- "rate_start_day"
    rend <- "rate_end_day"
    p1s <- "pulse1_start_day"
    p1e <- "pulse1_end_day"
    p2s <- "pulse2_start_day"
    p2e <- "pulse2_end_day"
    sepv <- "separator_day"
  } else {
    gstart <- "growth_start_season_day"
    gend <- "growth_end_season_day"
    rstart <- "rate_start_season_day"
    rend <- "rate_end_season_day"
    p1s <- "pulse1_start_season_day"
    p1e <- "pulse1_end_season_day"
    p2s <- "pulse2_start_season_day"
    p2e <- "pulse2_end_season_day"
    sepv <- "separator_season_day"
  }

  make_seg <- function(df, x1, x2, label) {
    if (!all(c(x1, x2) %in% names(df))) {
      return(data.frame())
    }

    df %>%
      dplyr::transmute(
        facet_var = .data$facet_var,
        row_label = .data$row_label,
        x = .data[[x1]],
        xend = .data[[x2]],
        window = label
      ) %>%
      dplyr::filter(!is.na(.data$x), !is.na(.data$xend))
  }

  make_ticks <- function(df, xcol, label) {
    if (!xcol %in% names(df)) {
      return(data.frame())
    }

    df %>%
      dplyr::transmute(
        facet_var = .data$facet_var,
        row_label = .data$row_label,
        x = .data[[xcol]],
        window = label
      ) %>%
      dplyr::filter(!is.na(.data$x))
  }

  make_pts <- function(df, xcol, label, tri = NULL) {
    if (!xcol %in% names(df)) {
      return(data.frame())
    }

    out <- df %>%
      dplyr::transmute(
        facet_var = .data$facet_var,
        row_label = .data$row_label,
        x = .data[[xcol]],
        window = label
      ) %>%
      dplyr::filter(!is.na(.data$x))

    if (!is.null(tri)) {
      out$triangle <- tri
    }

    out
  }

  overall_segments <- make_seg(fs, gstart, gend, "overall season")

  overall_ticks <- dplyr::bind_rows(
    make_ticks(fs, gstart, "overall season"),
    make_ticks(fs, gend, "overall season")
  )

  rate_points <- dplyr::bind_rows(
    make_pts(fs, rstart, "rate season", "\u25C0"),
    make_pts(fs, rend, "rate season", "\u25B6")
  )

  pulse_segments <- dplyr::bind_rows(
    make_seg(fs, p1s, p1e, "pulse 1"),
    make_seg(fs, p2s, p2e, "pulse 2")
  )

  separator_points <- make_pts(fs, sepv, "separator")

  list(
    x_mode = x_mode,
    facet_enabled = facet_enabled,
    y_label = y_label,
    base_rows = fs,
    overall_segments = overall_segments,
    overall_ticks = overall_ticks,
    rate_points = rate_points,
    pulse_segments = pulse_segments,
    separator_points = separator_points
  )
}