plot_hourly_single_var <- function(data, target, label){
  map_summary_df <- data |>
    group_by(Hour) |>
    summarize(
      mean_value = mean(.data[[target]],na.rm=TRUE),
      ci_lower = t.test(.data[[target]])$conf.int[1],
      ci_upper = t.test(.data[[target]])$conf.int[2],
      n=n(),
      .groups = 'drop'
    )
  overall_mean = mean(map_summary_df[["mean_value"]],na.rm=FALSE)
  
  p <- ggplot(map_summary_df,aes(x=Hour,y=mean_value)) + 
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper),
                fill = "#BFD7ED", alpha = 0.6) +   
    geom_line(color = "#73A9E2", linewidth = 1) +  # main line
    geom_point(color = "#2A5783", size = 2) +      
    geom_hline(yintercept = overall_mean,
               linetype = "dashed", color = "grey40") +
    scale_x_continuous(breaks = seq(0, 23, by = 2)) +
    labs(
      #title = paste("Mean ", label, " by Hour"),
      x = "Hour",
      y = paste("Mean ",  label)) +
    theme_bw(base_family = "Arial", base_size = 10) +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
      axis.title = element_text(size = 12),
      axis.text = element_text(size = 10),
      legend.title = element_text(size = 10),
      legend.text = element_text(size = 10),
      panel.border = element_blank(),              
      axis.line = element_line(),
      axis.line.x.top = element_blank(),          
      axis.line.y.right = element_blank()
    )
  return(p)
}

#function for plotting 2 arousal vars in the same plot
plot_hourly_vars <- function(music_data, targets, labels = NULL) {
  
  if (is.null(labels)) {
    labels <- targets
  }
  names(labels) <- targets
  
  # EXACT COLORS (max 3 variables)
  fixed_colors <- c(
    "#F28E2B",  
    "#4A90E2",  
    "#1b9e77"   
  )
  
  fixed_fills <- c(
    "#F4A261",  
    "#6BAED6",  
    "#a6dba0"   
  )
  
  colors <- fixed_colors[seq_along(targets)]
  fills  <- fixed_fills[seq_along(targets)]
  
  names(colors) <- targets
  names(fills)  <- targets
  
  # Compute summary statistics
  summary_df <- music_data |>
    dplyr::select(Hour, dplyr::all_of(targets)) |>
    tidyr::pivot_longer(
      cols = dplyr::all_of(targets),
      names_to = "variable",
      values_to = "value"
    ) |>
    dplyr::group_by(Hour, variable) |>
    dplyr::summarize(
      mean_value = mean(value, na.rm = TRUE),
      ci_lower = t.test(value)$conf.int[1],
      ci_upper = t.test(value)$conf.int[2],
      .groups = "drop"
    )
  
  overall_mean <- summary_df |>
    dplyr::group_by(variable) |>
    dplyr::summarize(overall_mean = mean(mean_value), .groups = "drop")
  
  p <- ggplot(summary_df, aes(x = Hour, y = mean_value, color = variable)) +
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper, fill = variable),
                alpha = 0.25, color = NA) +
    geom_line(linewidth = 1) +
    geom_point(size = 2) +
    geom_hline(data = overall_mean,
               aes(yintercept = overall_mean, color = variable),
               linetype = "dashed",
               linewidth = 0.8) +
    scale_x_continuous(breaks = seq(0, 23, by = 2)) +
    scale_color_manual(values = colors, labels = labels) +
    scale_fill_manual(values = fills, labels = labels) +
    coord_cartesian(ylim = c(-0.4, 0.4)) +
    labs(
      #title = "Mean Values by Hour",
      x = "Hour",
      y = "Mean value",
      color = "Variable",
      fill = "Variable"
    ) +
    theme_bw(base_family = "Arial", base_size = 10) +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
      axis.title = element_text(size = 12),
      axis.text = element_text(size = 10),
      legend.title = element_text(size = 10),
      legend.text = element_text(size = 10),
      panel.border = element_blank(),              
      axis.line = element_line(),
      axis.line.x.top = element_blank(),          
      axis.line.y.right = element_blank()
    )
  return(p)
}
plot_hourly_freq <- function(data, target, label) {
  freq_df <- data %>%
    dplyr::filter(!is.na(.data[[target]])) %>%
    group_by(Hour, .data[[target]]) %>%
    summarize(freq = n(), .groups = "drop")
  p <- ggplot(freq_df, aes(x = Hour, y = freq, color = as.factor(.data[[target]]))) +
    geom_line(size = 1) +
    scale_color_manual(values = c(
      "#F28E2B",
      "#4A90E2",
      "#4DB77D"
    )) +
    geom_point(size = 2) +
    scale_x_continuous(breaks = seq(0, 23, by = 2)) +
    labs(
      title = paste("Frequency of", label, "by Hour"),
      x = "Hour",
      y = "Frequency",
      color = target
    ) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
  return(p)
}
plot_hourly_percentage <- function(data, target, label, legend_labels = NULL) {
  library(dplyr)
  library(ggplot2)
  
  # Compute proportion per hour
  perc_df <- data %>%
    dplyr::filter(!is.na(.data[[target]])) %>%
    group_by(Hour, .data[[target]]) %>%
    summarize(count = n(), .groups = "drop") %>%
    group_by(Hour) %>%
    mutate(percentage = count / sum(count) * 100) %>%
    ungroup()
  
  # Default colors (up to 3 categories)
  my_colors <- c("#F28E2B", "#4A90E2", "#4DB77D")
  
  p <- ggplot(perc_df, aes(x = Hour, y = percentage, color = as.factor(.data[[target]]))) +
    geom_line(size = 1) +
    geom_point(size = 2) +
    scale_color_manual(
      values = my_colors,
      labels = legend_labels  # custom legend labels
    ) +
    scale_x_continuous(breaks = seq(0, 23, by = 2)) +
    labs(
      #title = paste("Percentage of", label, "by Hour"),
      x = "Hour",
      y = "Percentage (%)",
      color = label  # legend title
    ) +
    theme_bw() +
    theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
  axis.title = element_text(size = 10),
  axis.text = element_text(size = 8),
  legend.title = element_text(size = 8),
  legend.text = element_text(size = 8),
  panel.border = element_blank(),             
  axis.line = element_line(),
  axis.line.x.top = element_blank(),          
  axis.line.y.right = element_blank()
  )
  return(p)
}
plot_hourly_success <- function(data, var1, var2, var3) {
  success_df <- data %>%
    select(Hour, all_of(c(var1, var2, var3))) %>%
    pivot_longer(
      cols = -Hour,
      names_to = "variable",
      values_to = "value"
    ) %>%
    dplyr::filter(!is.na(value)) %>%
    group_by(Hour, variable) %>%
    summarize(
      success_rate = mean(value == 1),
      .groups = "drop"
    )
  p <- ggplot(success_df, aes(x = Hour, y = success_rate, color = variable)) +
    geom_line(linewidth = 1) +
    geom_point(size = 2) +
    scale_color_manual(values = c(
      "#F28E2B",
      "#4A90E2",
      "#4DB77D"
    )) +
    scale_x_continuous(breaks = seq(0, 23, by = 2)) +
    labs(
      title = "Success rate by Hour",
      x = "Hour",
      y = "Proportion of successes (value = 1)",
      color = "Success type"
    ) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
  return(p)
}

fit_cosinor_harmonics <- function(t, y, nHarm, plot = FALSE, label = NULL, color = "#4A90E2") {
  
  # Input checks 
  if (length(t) != length(y)) stop("t and y must have the same length")
  if (!is.numeric(nHarm) || nHarm < 1 || nHarm != round(nHarm)) {
    stop("nHarm must be a positive integer")
  }
  
  # Remove missing/non-finite
  valid_idx <- which(is.finite(t) & is.finite(y))
  t <- as.numeric(t[valid_idx])
  y <- as.numeric(y[valid_idx])
  if (length(t) == 0) stop("No valid (finite) observations after filtering")
  n <- length(t)
  
  # Build design matrix
  omega <- 2 * pi / 24
  X <- matrix(1, nrow = n, ncol = 1)
  for (k in 1:nHarm) {
    X <- cbind(X, cos(k * omega * t), sin(k * omega * t))}
  
  # Fit 
  beta <- solve(t(X) %*% X, t(X) %*% y)
  yhat <- as.vector(X %*% beta)
  resid <- y - yhat
  
  # Parameters
  MESOR <- beta[1]
  amplitude <- phase_rad <- phase_hours <- period_hours <- numeric(nHarm)
  beta_cos <- beta_sin <- numeric(nHarm)
  for (k in 1:nHarm) {
    ak <- beta[2*k]
    bk <- beta[2*k + 1]
    
    beta_cos[k] <- ak
    beta_sin[k] <- bk
    amplitude[k] <- sqrt(ak^2 + bk^2)
    phase_rad[k] <- atan2(bk, ak)
    
    period_hours[k] <- 24 / k
    phase_hours[k] <- (phase_rad[k] / (k * omega)) %% period_hours[k]}
  
  # Fit stats 
  SSE <- sum(resid^2)
  SST <- sum((y - mean(y))^2)
  R2 <- 1 - SSE / SST
  
  p <- ncol(X)
  adjR2 <- if (n > p) 1 - (SSE/(n - p)) / (SST/(n - 1)) else NA
  rmse <- if (n > p) sqrt(SSE / (n - p)) else NA
  
  # Optional plot
  if (plot) {
    df_plot <- data.frame(
      Hour = t,
      Observed = y,
      Fitted = yhat)
    p_plot <- ggplot(df_plot, aes(x = Hour)) +
      geom_line(aes(y = Fitted), color = color, size = 1) +
      labs(
        title = paste("Cosinor Fit (", nHarm, " harmonics)", if (!is.null(label)) paste("-", label) else ""),
        x = "Hour",
        y = "Observed / Fitted"
      ) +
      theme_bw() +
      theme(plot.title = element_text(hjust = 0.5, face = "bold"))
    print(p_plot)}
  
  # Output
  return(list(
    MESOR = MESOR,
    beta_cos = beta_cos,
    beta_sin = beta_sin,
    amplitude = amplitude,
    phase_rad = phase_rad,
    phase_hours = phase_hours,
    period_hours = period_hours,
    beta = as.vector(beta),
    X = X,
    yhat = yhat,
    resid = resid,
    SSE = SSE,
    SST = SST,
    R2 = R2,
    adjR2 = adjR2,
    rmse = rmse,
    n = n,
    p = p,
    nHarm = nHarm,
    t = t,
    y = y))}
fit_cosinor_harmonics_compare_plot <- function(
    t, y, maxHarm = 4, label, plotColors = NULL, criterion = "AIC"
) {
  # Clean data
  valid_idx <- which(is.finite(t) & is.finite(y))
  t <- t[valid_idx]
  y <- y[valid_idx]
  
  # Results storage
  results <- data.frame(
    Harmonics = integer(),
    R2 = numeric(),
    adjR2 = numeric(),
    AIC = numeric(),
    BIC = numeric(),
    stringsAsFactors = FALSE)
  # Colors
  if (is.null(plotColors)) {
    plotColors <- scales::hue_pal()(maxHarm)
  }
  # Store fitted curves
  fit_curves <- data.frame(Hour = t)
  
  # Loop over harmonics
  for (nHarm in 1:maxHarm) {
    fit <- fit_cosinor_harmonics(t, y, nHarm)
    n <- fit$n
    p <- fit$p
    SSE <- fit$SSE
    # AIC & BIC
    AIC_val <- n * log(SSE / n) + 2 * p
    BIC_val <- n * log(SSE / n) + log(n) * p
    
    results <- rbind(results, data.frame(
      Harmonics = nHarm,
      R2 = fit$R2,
      adjR2 = fit$adjR2,
      AIC = AIC_val,
      BIC = BIC_val))
    fit_curves[[paste0("Harm", nHarm)]] <- fit$yhat}
  
  # Reshape 
  fit_long <- fit_curves %>%
    pivot_longer(
      cols = -Hour,
      names_to = "Harmonics",
      values_to = "Fitted")
  
  p <- ggplot(fit_long, aes(x = Hour, y = Fitted, color = Harmonics)) +
    geom_line(size = 1) +
    scale_color_manual(values = plotColors) +
    labs(
      title = paste("Cosinor Fits with Different Harmonics", label),
      x = "Hour",
      y = "Fitted"
    ) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5, face = "bold"))
  
  # Choose best model
  criterion <- toupper(criterion)
  if (!criterion %in% c("AIC", "BIC")) {
    stop("criterion must be 'AIC' or 'BIC'")}
  best_model <- results$Harmonics[which.min(results[[criterion]])]
  message("Best model by ", criterion, ": ", best_model, " harmonics")
  return(list(
    summary = results,
    plot = p,
    best = best_model))}
get_cosinor_24h <- function(t, y, nHarm = 1) {
  
  # Fit cosinor
  fit <- fit_cosinor_harmonics(t, y, nHarm)
  # Create exact hourly grid
  t_grid <- 0:23
  # Rebuild design matrix
  omega <- 2 * pi / 24
  X_grid <- matrix(1, nrow = length(t_grid), ncol = 1)
  for (k in 1:nHarm) {
    X_grid <- cbind(
      X_grid,
      cos(k * omega * t_grid),
      sin(k * omega * t_grid))}
  
  # Predict 
  yhat_24 <- as.vector(X_grid %*% fit$beta)
  
  return(data.frame(
    Hour = t_grid,
    yhat = yhat_24))}
compute_phase_diff <- function(t1, y1, t2, y2, n_points = 1000) {

    # Common time grid
  t_common <- seq(0, 24, length.out = n_points)
  # Interpolate using REAL time 
  y1_interp <- approx(t1, y1, xout = t_common, rule = 2)$y
  y2_interp <- approx(t2, y2, xout = t_common, rule = 2)$y
  
  #  Hilbert phase
  ang1 <- Arg(gsignal::hilbert(y1_interp))
  ang2 <- Arg(gsignal::hilbert(y2_interp))
  
  # Phase difference
  dphi <- (ang2 - ang1 + pi) %% (2*pi) - pi
  
  # Convert to hours
  phase_hours <- 24 * dphi / (2*pi)
  
  return(data.frame(
    Hour = t_common,
    PhaseDiff_hours = phase_hours))}
plot_hourly_single_var_and_cosinor <- function(data, target, label, nHarm = NULL){
  
  # Summary data (hourly means + CI) 
  map_summary_df <- data |>
    group_by(Hour) |>
    summarize(
      mean_value = mean(.data[[target]], na.rm = TRUE),
      ci_lower = t.test(.data[[target]])$conf.int[1],
      ci_upper = t.test(.data[[target]])$conf.int[2],
      n = n(),
      .groups = 'drop')
  overall_mean <- mean(map_summary_df[["mean_value"]], na.rm = FALSE)
  
  #Base plot
  p <- ggplot(map_summary_df, aes(x = Hour, y = mean_value)) + 
    geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper),
                fill = "#BFD7ED", alpha = 0.6) +
    geom_line(color = "#73A9E2", linewidth = 1) +
    geom_point(color = "#2A5783", size = 2) +
    geom_hline(yintercept = overall_mean,
               linetype = "dashed", color = "grey40") +
    scale_x_continuous(breaks = seq(0, 23, by = 2)) +
    labs(
      title = paste("Mean", label, "by Hour"),
      x = "Hour",
      y = paste("Mean", label)
    ) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14))
  # Add cosinor curve 
  if (!is.null(nHarm)) {
    t <- data$Hour
    y <- data[[target]]
    # Fit cosinor
    fit <- fit_cosinor_harmonics(t, y, nHarm)
    # Smooth time grid
    t_smooth <- seq(0, 23, length.out = 1000)
    # Build design matrix
    omega <- 2 * pi / 24
    X_smooth <- matrix(1, nrow = length(t_smooth), ncol = 1)
    for (k in 1:nHarm) {
      X_smooth <- cbind(
        X_smooth,
        cos(k * omega * t_smooth),
        sin(k * omega * t_smooth))}
    
    # Predict smooth curve
    y_smooth <- as.vector(X_smooth %*% fit$beta)
    cos_df <- data.frame(
      Hour = t_smooth,
      Cosinor = y_smooth)
    
    # Add to plot
    p <- p +
      geom_line(
        data = cos_df,
        aes(x = Hour, y = Cosinor),
        color = "#D55E00",
        linewidth = 1.2)}
  return(p)}

plot_multiple_cosinor <- function(spec, data_env) {
  omega <- 2 * pi / 24
  t_smooth <- seq(0, 24, length.out = 1000)
  all_curves <- data.frame()
  for (i in 1:nrow(spec)) {
    df_name <- spec$df[i]
    var     <- spec$var[i]
    label   <- spec$label[i]
    nHarm   <- spec$nHarm[i]
    data <- data_env[[df_name]]
    if (is.null(data)) stop(paste("Dataframe", df_name, "not found"))
    t <- data$Hour
    y <- data[[var]]
    fit <- fit_cosinor_harmonics(t, y, nHarm)
    # Build smooth curve
    X <- matrix(1, nrow = length(t_smooth), ncol = 1)
    for (k in 1:nHarm) {
      X <- cbind(X,
                 cos(k * omega * t_smooth),
                 sin(k * omega * t_smooth))}
    y_smooth <- as.vector(X %*% fit$beta)
    df_out <- data.frame(
      Hour = t_smooth,
      Value = y_smooth,
      Variable = label)
    all_curves <- rbind(all_curves, df_out)}
  ggplot(all_curves, aes(x = Hour, y = Value, color = Variable)) +
    geom_line(linewidth = 1.2) +
    scale_x_continuous(breaks = seq(0, 24, by = 2)) +
    labs(
      #title = "Cosinor Fits",
      x = "Hour",
      y = "Magnitude"
    ) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
          axis.title = element_text(size = 12),
          axis.text = element_text(size = 10),
          legend.title = element_text(size = 10),
          legend.text = element_text(size = 10),
          panel.border = element_blank(),             
          axis.line = element_line(),
          axis.line.x.top = element_blank(),          
          axis.line.y.right = element_blank())}

generate_cosinor_1440 <- function(fit) {
  omega <- 2 * pi / 24
  # Time grid (in hours)
  t_hours <- seq(0, 24, length.out = 1440)
  # Initialize
  yhat <- rep(fit$MESOR, length(t_hours))
  velocity <- rep(0, length(t_hours))
  acceleration <- rep(0, length(t_hours))
  
  # Build signal + derivatives
  for (k in 1:fit$nHarm) {
    ak <- fit$beta_cos[k]
    bk <- fit$beta_sin[k]
    # Signal
    yhat <- yhat +
      ak * cos(k * omega * t_hours) +
      bk * sin(k * omega * t_hours)
    # First derivative (velocity)
    velocity <- velocity +
      (-ak * k * omega * sin(k * omega * t_hours)) +
      ( bk * k * omega * cos(k * omega * t_hours))
    
    # Second derivative (acceleration)
    acceleration <- acceleration +
      (-ak * (k * omega)^2 * cos(k * omega * t_hours)) +
      (-bk * (k * omega)^2 * sin(k * omega * t_hours))}
  
  return(data.frame(
    MinuteOfDay = 0:1439,
    Hour = t_hours,
    Energy = yhat,
    Velocity = velocity,
    Acceleration = acceleration))}

# attach cosinor column to music data
attach_cosinor_energy <- function(music_data, cos_1440) {
  time_to_minute <- function(datetime) {
    lt <- as.POSIXlt(datetime)
    lt$hour * 60 + lt$min}
  music_data$MinuteOfDay <- time_to_minute(music_data$InitStart)
  music_data$CosinorEnergy <- cos_1440$Energy[
    match(music_data$MinuteOfDay, cos_1440$MinuteOfDay)]
  music_data$CosinorVelocity <- cos_1440$Velocity[
    match(music_data$MinuteOfDay, cos_1440$MinuteOfDay)]
  music_data$CosinorAcceleration <- cos_1440$Acceleration[
    match(music_data$MinuteOfDay, cos_1440$MinuteOfDay)]
  return(music_data)}

run_cosinor_first_harmonic <- function(
    ref_df,
    ref_var,
    ref_time = "timeScale",
    ref_harm = 4,
    music_df,
    target_vars,
    harmonics,
    n_points = 1000,   # resolution for cross-corr
    file_name = NULL
) {
  library(dplyr)
  if (length(target_vars) != length(harmonics)) {
    stop("target_vars and harmonics must have same length")
  }
  results_list <- list()
  # REF signal
  ref_agg <- ref_df %>%
    group_by(.data[[ref_time]]) %>%
    summarise(y = mean(.data[[ref_var]], na.rm = TRUE), .groups = "drop")
  
  fit_ref <- fit_cosinor_harmonics(
    ref_agg[[ref_time]],
    ref_agg$y,
    nHarm = ref_harm)
  
  # interpolate reference once
  t_common <- seq(0, 24, length.out = n_points)
  ref_interp <- approx(fit_ref$t, fit_ref$yhat, xout = t_common, rule = 2)$y
  ref_interp <- as.numeric(scale(ref_interp))
  
  ref_row <- data.frame(
    Variable = paste0("Mappiness (", ref_var, ")"),
    Harmonics_Fit = ref_harm,
    MESOR = fit_ref$MESOR,
    Amplitude_24h = fit_ref$amplitude[1],
    Phase_24h = fit_ref$phase_hours[1],
    Period_24h = fit_ref$period_hours[1],
    R2 = fit_ref$R2,
    adjR2 = fit_ref$adjR2,
    RMSE = fit_ref$rmse,
    r = NA)
  results_list[[1]] <- ref_row
  
  # LOOP MUSIC
  for (i in seq_along(target_vars)) {
    var_name <- target_vars[i]
    harm <- harmonics[i]
    var_agg <- music_df %>%
      group_by(.data[[ref_time]]) %>%
      summarise(y = mean(.data[[var_name]], na.rm = TRUE), .groups = "drop")
    
    fit <- fit_cosinor_harmonics(
      var_agg[[ref_time]],
      var_agg$y,
      nHarm = harm)
    
    # interpolate
    mus_interp <- approx(fit$t, fit$yhat, xout = t_common, rule = 2)$y
    mus_interp <- as.numeric(scale(mus_interp))
    
    # cross-correlation
    r_val <- cor(ref_interp, mus_interp, use = "complete.obs")
    df <- data.frame(
      Variable = var_name,
      Harmonics_Fit = harm,
      MESOR = fit$MESOR,
      Amplitude_24h = fit$amplitude[1],
      Phase_24h = fit$phase_hours[1],
      Period_24h = fit$period_hours[1],
      R2 = fit$R2,
      adjR2 = fit$adjR2,
      RMSE = fit$rmse,
      r = r_val)
    results_list[[i + 1]] <- df}
  
  # combine 
  results_table <- do.call(rbind, results_list)
  
  if (!is.null(file_name)) {
    write.csv(results_table, file_name, row.names = FALSE)}
  
  return(results_table)
}