### segment-level code (also applicable to single units):
##' extract segment transitions
##' (1) open the Supplementary_Data_1 (sample data) and download it to your device
##' (2) define the filename in the code below (see #!) by setting the path file for where you 
##'     saved the Supplementary_Data_1 
##' (3) run the function to store it in the environment
##' 

syllable_transition_matrix <- function(names, outfits, lengths, x_positions, y_positions, random_run = FALSE) {
  
  two_letter_codes <- c('CC', 'GR', 'RC', 'KS', 'SH', 'SI', 'RO', 'KH')
  
  if (length(x_positions)!=length(y_positions)) {
    stop("These must be equal length.")
  }
  
  filename <- "PATH_TO/Supplementary_Data_1.csv"
  #! INSERT YOUR PATHFILE FOR Supplementary_Data_1 HERE^
  data <- read.csv(filename)
  data[,1] <- gsub(" ", "", data[,1], fixed = TRUE)
  
  transitions <- data.frame(Transition = character(), Frequency = integer())
  
  for (name in names) {
    for (outfit in outfits) {
      name_match <- data[, 1] == name
      outfit_match <- data[, 2] == outfit
      both_match <- name_match & outfit_match
      if (sum(both_match>0)) {
        relevant_data <- data[both_match, 3]
        
        for (i in 1:length(relevant_data)) {
          entry <- relevant_data[i]
          split_entry <- unlist(strsplit(entry, '_'))
          number_entries <- length(split_entry)
          if (random_run) {
            split_entry <- sample(split_entry,number_entries)
          }
          if (number_entries %in% lengths) { # checks if number_entries is in the lengths vector
            # extracts transitions for specified x and y positions
            for (idx in seq_along(x_positions)) {
              x_pos <- as.numeric(strsplit(x_positions[idx], "_")[[1]])
              y_pos <- as.numeric(strsplit(y_positions[idx], "_")[[1]])
              if (all(x_pos <= number_entries) & all(y_pos <= number_entries)) {
                x <- split_entry[x_pos]
                y <- split_entry[y_pos]
                transition <- paste("x:", paste(x, collapse = "_"), "y:", paste(y, collapse = "_"))
                
                if (all(x %in% two_letter_codes) & all(y %in% two_letter_codes)) {
                  # checks if transition is already in the data frame
                  if (transition %in% transitions$Transition) {
                    transitions$Frequency[transitions$Transition == transition] <- transitions$Frequency[transitions$Transition == transition] + 1
                  } else {
                    transitions <- rbind(transitions, data.frame(Transition = transition, Frequency = 1))
                  }
                } else {
                  print(paste(x,y))
                  stop("Did not think this would happen.")
                }
              } 
            }
          }
        }
      }
    }
  }
  
  # sorts transitions by frequency in descending order
  transitions <- transitions[order(-transitions$Frequency), ]
  
  # prints transitions and their frequencies
  if (!random_run) {
    print(transitions)
  }
  
  if (!random_run) {
    print(paste("The total number of unique transitions is:", nrow(transitions)))
  }
  
  # creates transitions matrix
  # extracts unique x and y values from transitions dataframe
  x_vector <- c()
  y_vector <- c()
  for (transition in transitions$Transition) {
    decoded <- decode_transition(transition)
    x_vector <- c(x_vector,decoded[1])
    y_vector <- c(y_vector,decoded[2])
  }
  unique_x <- unique(x_vector)
  unique_y <- unique(y_vector)
  
  # creates empty transitions matrix
  transition_matrix <- matrix(0, nrow = length(unique_x), ncol = length(unique_y))
  rownames(transition_matrix) <- unique_x
  colnames(transition_matrix) <- unique_y
  
  # populates transitions matrix using frequencies from transitions dataframe
  for (i in 1:nrow(transitions)) {
    transition <- transitions$Transition[i]
    frequency <- transitions$Frequency[i]
    decoded <- decode_transition(transition)
    x <- decoded[1]
    y <- decoded[2]
    transition_matrix[x, y] <- frequency
  }
  
  if (random_run) {
    return(transition_matrix)
  } else {
    return(list(Transitions = transitions, TransitionMatrix = transition_matrix))
  }
}

decode_transition <- function(transition) {
  x_string <- unlist(strsplit(transition,"y:"))[1]
  y_string <- unlist(strsplit(transition,"y:"))[2]
  x_string <- gsub("x:","",x_string)
  x_string <- gsub(" ","",x_string)
  y_string <- gsub(" ","",y_string)
  return(c(x_string,y_string))
}




##' generate random matrix:
##' (1) after completing the previous step, check to see if the syllable_transition_matrix function 
##'     is stored in the environment ->
##' (2) run the get_random_syllable_transition_matrix function below

get_random_syllable_transition_matrix <- function(N = 1000, names, outfits, lengths, x_positions, y_positions) {
  pb = txtProgressBar(min = 0, max = N, initial = 0) 
  all_matrices <- list()
  all_row_names <- c()
  all_col_names <- c()
  for (iteration in 1:N) {
    one_matrix <- syllable_transition_matrix(names, outfits, lengths, x_positions, y_positions, TRUE)
    all_matrices[[iteration]] <- one_matrix
    all_row_names <- c(all_row_names,rownames(one_matrix))
    all_col_names <- c(all_col_names,colnames(one_matrix))
    all_row_names <- unique(all_row_names)
    all_col_names <- unique(all_col_names)
    setTxtProgressBar(pb,iteration)
  }
  new_all_matrices <- list()
  for (i in 1:N) {
    basic_matrix <- all_matrices[[i]]
    new_all_matrices[[i]] <- convert_matrix(basic_matrix,all_row_names,all_col_names)
  }
  one_matrix <- new_all_matrices[[1]]
  number_rows <- nrow(one_matrix)
  number_cols <- ncol(one_matrix)
  result_mean <- one_matrix
  result_mean[,] = 0
  result_sd <- one_matrix
  result_sd[,] = 0
  for (row in 1:number_rows) {
    for (column in 1:number_cols) {
      vector_at_position <- c()
      for (i in 1:N) {
        vector_at_position <- c(vector_at_position,new_all_matrices[[i]][row,column])
      }
      #if ((row==3)&&(column==7)) {
      #hist(vector_at_position)
      #}
      result_mean[row,column] = mean(vector_at_position)
      result_sd[row,column] = sd(vector_at_position)
    }
  }
  overall <- list(result_mean,result_sd)
}

convert_matrix <- function(basic_matrix,rownames,colnames) {
  result <- matrix(0,nrow=length(rownames),ncol=length(colnames))  
  rownames(result) <- rownames
  colnames(result) <- colnames
  if (all(rownames(basic_matrix) %in% rownames) & all(colnames(basic_matrix) %in% colnames)) {
    for (row in rownames(basic_matrix)) {
      for (col in colnames(basic_matrix)) {
        result[row,col] <- basic_matrix[row,col]
      }
    }
    return(result)
  } else {
    stop("Matrices incompatible; increase N of iterations.")
  }
}



##' compare expected vs real transition frequencies
##' (1) check to see if get_random_syllable_transition_matrix is stored in the environment ->
##' (2) run the compare_matrices function below

compare_matrices <- function(syllable_matrix, random_matrix) {
  num_rows <- nrow(syllable_matrix)
  num_cols <- ncol(syllable_matrix)
  comparison_matrix <- array(0, dim = c(num_rows, num_cols, 2))  # Include a third dimension for p-values and expected frequencies
  rownames(comparison_matrix) <- rownames(syllable_matrix)
  colnames(comparison_matrix) <- colnames(syllable_matrix)
  for (i in 1:num_rows) {
    for (j in 1:num_cols) {
      # extracts frequency from syllable_matrix
      frequency <- syllable_matrix[i, j]
      # extracts mean from random_matrix
      mean_value <- random_matrix[[1]][i, j]
      # calculates expected frequency from mean
      expected_frequency <- mean_value
      # performs chi-square test approximation
      chi_square_stat <- (frequency - expected_frequency)^2 / expected_frequency
      # calculates p-value using chi-square distribution
      p_value <- pchisq(chi_square_stat, df = 1, lower.tail = FALSE)
      # stores p-value and expected frequency in comparison_matrix
      comparison_matrix[i, j, 1] <- p_value
      comparison_matrix[i, j, 2] <- expected_frequency
    }
  }
  return(list(p_values = comparison_matrix[,,1], expected_frequencies = comparison_matrix[,,2]))
}


##' combine functions:
##' (1) check all three functions (syllable_transition_matrix, get_random_syllable_transition_matrix, and 
##'     compare_matrices) are stored in the environment ->
##' (2) run the function below to combine all three functions

combined_function <- function(names, outfits, lengths, x_positions, y_positions, N = 10) {
  # runs syllable matrix
  syllable_matrix <- syllable_transition_matrix(names, outfits, lengths, x_positions, y_positions, FALSE)
  # runs randomisation
  random_matrix <- get_random_syllable_transition_matrix(N, names, outfits, lengths, x_positions, y_positions)
  # converts matrix
  converted_syllable_matrix <- convert_matrix(syllable_matrix$TransitionMatrix,
                                              rownames(random_matrix[[1]]),
                                              colnames(random_matrix[[2]]))
  # runs comparison
  comparison_results <- compare_matrices(syllable_matrix$TransitionMatrix, random_matrix)
  # returns syllable matrix and p-values from comparison results
  return(list(syllable_matrix = syllable_matrix, p_values = comparison_results$p_values))
}


##'run the functions to get results:
##'(1) to specify the predator models you'd like to group the data by, adjust the #predator 
##'     model(s) vector, and to the set the individuals you'd like to group by, 
##'     adjust the #individual(s) vector The sample dataset provided only contains Chris' response 
##'     to the Pattern model, so the function input has been pre-set accordingly.
##'(2) to specify the sequence lengths you'd like to include, adjust the #sequence lengths vector,
##'     for example, if you'd only like to look at 3- 4-, 5-, and 6-grams (that is sequences with 3 
##'     to 6 units), you would set it as c(3, 4, 5, 6)
##'(3) decide which positions you'd like to see the call pairs for, for example, if you'd like
##'     to see the call pairs from the first bigram to the  third bigram of an n-gram
##'     (A_B_x_x_x_x to x_x_C_D_x_x where A_B is the first bigram and C_D is the third), 
##'     you would set the #x position(s) vector as c("1_2") and the y position(s) vector as c("4_5")
##'       ! note that the #x position(s) and #y position(s) vectors are pairwise such that if the 
##'         x_positions were c("1_2", "2_3") and the y_positions were c("3_4", "4_5"), the call pairs 
##'         for 1_2 to 3_4 would be extracted, and the call pairs for 2_3 to 4_5 would be extracted
##'(4) for the random_matrix function, you'll also have the option of setting the number of permutations
##'     that are run to create your null matrix, here it has been set to 1000
##'(5) compare the randomised syllable matrix with the actual matrix by running the next function 
##'(6) use the final function to perform all above steps in one go 

# get syllable matrix 
syllable_matrix <- syllable_transition_matrix(
  c("Chris"), #individual(s)
  c("Pattern"), #predator model(s)
  c(3, 4, 5, 6), #sequence length(s)
  c("1_2"), #x position(s)
  c("3_4"), #y position(s)
  FALSE) 
syllable_matrix

# get randomised syllable matrices
random_matrix <- get_random_syllable_transition_matrix(N = 1000, #set permutations
                                                       c("Chris"), #individual(s)
                                                       c("Pattern"), #predator model(s)
                                                       c(3, 4, 5, 6), #sequence length(s)
                                                       c("1_2"), #x position(s)
                                                       c("3_4")) #y position(s)
random_matrix[[1]] # means
random_matrix[[2]] # standard deviation
converted_syllable_matrix <- convert_matrix(syllable_matrix,rownames(random_matrix[[1]]),colnames(random_matrix[[2]]))
converted_syllable_matrix

# compare actual syllable matrices with randomised syllable matrices
comparison_results <- compare_matrices(syllable_matrix$TransitionMatrix, random_matrix)
comparison_results$p_values  # display p-values
comparison_results$expected_frequencies  # display expected frequencies

# run the above three functions to complete all steps in one go using the combined function 
result <- combined_function(
  c("Chris"), #individual(s)         
  c("Pattern"), #predator model(s)
  c(3, 4, 5, 6),   #sequence length(s)
  c("1_2"),  #x position(s)       
  c("3_4"), #y position(s)
  N = 100) #set permutations
result$syllable_matrix
result$p_values

##' see the Supplementary_Data_2 fle for more examples on what vectors to input, e.g.:
##'   an x_position of c("1_2_3", "2_3_4") and a y_position of c("2_3_4", "3_4_5") would show 
##'   the transitions from ABC in A_B_C_x_x_x to DEF in x_D_E_F_x_x, and transitions from ABC in 
##'   x_A_B_C_x_x to DEF in x_x_D_E_F_x