#### From Lagrangian to Eulerian perspectives
#### Script: Functions to explore the tracking data
#### Last modification: 04 04 2024
#### R version 4.2.3 (2023-03-15 ucrt) -- "Shortstop Beagle"



### 3 - Summary stats of movements ---------------------------------------------
tracksummary_func <- function(trajectories, type = "CPF", n_reloc_per_hour){
  if(type == "CPF"){
    ## distributions of steps and angles
    cli::cli_alert("Computing the step and angle distributions")
    plot_distrib <- ggarrange(
      ggplot(trajectories |> subset(activity != "start")) + 
        geom_density(aes(x = angle, color = activity)) + 
        facet_grid(Colony~Day) + theme_bw() + 
        scale_color_colorblind(name = "Activity", labels = c("Forage", "Travel")) +
        labs(title = "Angle distribution") + 
        theme(axis.title.x = element_blank(),
              panel.grid = element_blank(), 
              strip.background = element_blank(),
              strip.text = element_text(face = "bold", hjust = 0),
              plot.title.position = "plot")
      ,
      ggplot(trajectories |> subset(activity != "start")) + 
        geom_density(aes(x = step, color = activity)) + 
        facet_grid(Colony~Day) + theme_bw() + 
        scale_color_colorblind(name = "Activity", labels = c("Forage", "Travel")) + 
        labs(title = "Step distribution") + 
        theme(axis.title.x = element_blank(),
              panel.grid = element_blank(), 
              strip.background = element_blank(),
              strip.text = element_text(face = "bold", hjust = 0),
              plot.title.position = "plot")
      , nrow = 2, align = "hv", common.legend = T, legend = "bottom") + bgcolor("white") + border("white")
    # ggsave("Simulated_tracks_step_angle_distribution.png", height = 15, width = 15, units = "cm")
    
    
    ## number of fixes per individual
    cli::cli_alert("Summarizing the trajectories per activity")
    smry_track_length <- trajectories |> group_by(BirdId, Day, Colony) |> tally()
    avg_duration <- mean((smry_track_length$n)/n_reloc_per_hour, na.rm = T)
    cli::cli_inform(glue::glue("Average duration of trip is {round(avg_duration, 1)} h"))
      
    
    plot_nfix <- ggplot(smry_track_length) + geom_density(aes(x = n)) + 
      theme_bw() + labs(x = "Number of fixes", y = "Density") + facet_wrap("Day")
    # ggsave("Simulated_tracks_nbr_fixes.png", height = 6, width = 8, units = "cm")
    
    
    ## number of fixes per activity type - relative to track length
    smry_ind <- trajectories |> group_by(Day, Colony, BirdId, activity) |> 
      tally() |>
      ungroup() |>
      left_join(smry_track_length, by = c("BirdId", "Day", "Colony")) |>
      rename(n_tot = n.y, n_act = n.x) |>
      mutate(proportion = n_act / n_tot) |> 
      subset(activity != "start")
    smry_ind$activity <- factor(smry_ind$activity, levels = c("travel", "forage"))
    
    # add individuals which have no foraging steps
    # smry_ind |> subset(Colony == 1) |> subset(activity == "forage") |> group_by(BirdId) |> nrow()
    id.tot <- unique(trajectories$BirdId) # BirdId contient le jour donc c'est bon
    id.sub <- unique(smry_ind$BirdId[which(smry_ind$activity == "forage")])
    id.miss <- setdiff(id.tot, id.sub)
    df.miss <- data.frame(BirdId = id.miss, activity = "forage", n_act = 0, proportion = 0) |>
      left_join(smry_track_length, by = "BirdId") |> rename(n_tot = n) %>%
      mutate(Colony = case_when(
        sapply(strsplit(.$BirdId, "_"), "[[", 1) == "Col1" ~ "1",
        sapply(strsplit(.$BirdId, "_"), "[[", 1) == "Col2" ~ "2"
      ))
    smry_ind <- plyr::rbind.fill(smry_ind, df.miss)
    smry_ind |> group_by(Colony, activity) |> summarize(med = median(proportion, na.rm = T),
                                                        mean = mean(proportion, na.rm = T),
                                                        min = min(proportion, na.rm = T),
                                                        max = max(proportion, na.rm = T),
                                                        sd = sd(proportion, na.rm = T))
    
    
    # plot
    plot_proptrack <- ggplot(smry_ind) + 
      geom_boxplot(aes(y = proportion, 
                       color = activity, 
                       x = as.factor(Colony)),
                       outlier.size = 0.5, lwd = 0.4, fatten = 1.2) + 
      theme_bw() + 
      scale_color_colorblind(name = "Activity", labels = c("Forage", "Travel")) + 
      facet_wrap("Day") +
      labs(x = "Colony", y = "Proportion of time") +
      theme(axis.title = element_text(size = 8),
            axis.text = element_text(size = 6),
            legend.text = element_text(size = 6),
            legend.title = element_text(size = 8))
    # ggsave("Simulated_tracks_nbr_fixes_by_activity.png", height = 7, width = 10, units = "cm")
    
    plot_tripduration <- ggplot(smry_ind) + 
      geom_boxplot(aes(y = n_tot/n_reloc_per_hour, x = as.factor(Colony), color = as.factor(Day))) + scale_color_colorblind()
    
    return(list(plot_distrib = plot_distrib, 
                plot_nfix = plot_nfix, 
                plot_proptrack = plot_proptrack, 
                plot_tripduration = plot_tripduration))
  }
  if(type == "FR"){
    ## distributions of steps and angles
    cli::cli_alert("Computing the step and angle distributions")
    plot_distrib <- ggarrange(
      ggplot(trajectories |> subset(activity != "start")) + 
        geom_density(aes(x = angle, color = activity)) + 
        facet_wrap("Day") + theme_bw() + 
        scale_color_colorblind(name = "Activity", labels = c("Forage", "Travel")) +
        labs(title = "A - Angle distribution") + 
        theme(axis.title.x = element_blank(),
              panel.grid = element_blank(), 
              strip.background = element_blank(),
              strip.text = element_text(face = "bold", hjust = 0))
      ,
      ggplot(trajectories |> subset(activity != "start")) + 
        geom_density(aes(x = step, color = activity)) + 
        facet_wrap("Day") + theme_bw() + 
        scale_color_colorblind(name = "Activity", labels = c("Forage", "Travel")) + 
        labs(title = "B - Step distribution") + 
        theme(axis.title.x = element_blank(),
              panel.grid = element_blank(), 
              strip.background = element_blank(),
              strip.text = element_text(face = "bold", hjust = 0))
      , nrow = 2, align = "hv", common.legend = T, legend = "bottom") + bgcolor("white") + border("white")
    # ggsave("Simulated_tracks_step_angle_distribution.png", height = 15, width = 15, units = "cm")
    
    
    ## number of fixes per individual
    cli::cli_alert("Summarizing the trajectories per activity")
    smry_track_length <- trajectories |> group_by(Day, ind_id) |> tally()
    avg_duration <- mean((smry_track_length$n)/60, na.rm = T)
    cli::cli_inform(glue::glue("Average duration of trip is {round(avg_duration, 1)} h"))
    
    plot_nfix <- ggplot(smry_track_length) + geom_density(aes(x = n)) + 
      theme_bw() + labs(x = "Number of fixes", y = "Density") + facet_wrap("Day")
    # ggsave("Simulated_tracks_nbr_fixes.png", height = 6, width = 8, units = "cm")
    
    
    ## number of fixes per activity type - relative to track length
    smry_ind <- trajectories |> group_by(Day, ind_id, activity) |> 
      tally() |>
      ungroup() |>
      left_join(smry_track_length, by = c("Day", "ind_id")) |>
      rename(n_tot = n.y, n_act = n.x) |>
      mutate(proportion = n_act / n_tot) |> 
      subset(activity != "start") |> 
      mutate(trip_id = paste(ind_id, Day, sep = "_"))
    smry_ind$activity <- factor(smry_ind$activity, levels = c("travel", "forage"))
    
    # add individuals which have no foraging steps
    # smry_ind |> subset(Colony == 1) |> subset(activity == "forage") |> group_by(BirdId) |> nrow()
    id.tot <- unique(trajectories$trip_id) 
    id.sub <- unique(smry_ind$trip_id[which(smry_ind$activity == "forage")])
    id.miss <- setdiff(id.tot, id.sub)
    df.miss <- data.frame(trip_id = id.miss, activity = "forage", n_act = 0, proportion = 0) |>
      left_join(smry_track_length |> mutate(trip_id = paste(ind_id, Day, sep = "_")), by = "trip_id") |> rename(n_tot = n) 
    smry_ind <- plyr::rbind.fill(smry_ind, df.miss) 
    smry_ind |> group_by(activity) |> summarize(med = median(proportion, na.rm = T),
                                                        mean = mean(proportion, na.rm = T),
                                                        min = min(proportion, na.rm = T),
                                                        max = max(proportion, na.rm = T),
                                                        sd = sd(proportion, na.rm = T))
    
    
    # plot
    plot_proptrack <- ggplot(smry_ind) + 
      geom_boxplot(aes(y = proportion, 
                       color = activity),
                   outlier.size = 0.5, lwd = 0.4, fatten = 1.2) + 
      theme_bw() + 
      scale_color_colorblind(name = "Activity", labels = c("Forage", "Travel")) + 
      facet_wrap("Day") +
      labs(y = "Proportion of time") +
      theme(axis.title = element_text(size = 8),
            axis.text = element_text(size = 6),
            legend.text = element_text(size = 6),
            legend.title = element_text(size = 8))
    # ggsave("Simulated_tracks_nbr_fixes_by_activity.png", height = 7, width = 10, units = "cm")
    
    plot_tripduration <- ggplot(smry_ind) + 
      geom_boxplot(aes(y = n_tot/n_reloc_per_hour, color = as.factor(Day))) + scale_fill_colorblind()
    
    return(list(plot_distrib = plot_distrib, 
                plot_nfix = plot_nfix, 
                plot_proptrack = plot_proptrack, 
                plot_tripduration = plot_tripduration))
  }
}


