library(stats)
# library(boot)
library(data.table)
# library(DescTools)
library(tidyverse)
# library(mixtools)
library(viridis)

uap = "U3alonePlus"
gup = "G1atU3Plus"

comtheme <- theme_bw() +
  theme(text = element_text(color = "#545454", size = 18),
        legend.position = "none")

#go to the analytical working directory
setwd("C:/USERS/MCS/Desktop/PULLS/R_sim_scripts/data/Analytical_environment")
list.files()




#form list of specific raw data as matricies, perform MSD analysis, plot,
#and merge MSDs for condition/sim
setwd(
  "C:/USERS/MCS/Desktop/PULLS/R_sim_scripts/data/Analytical_environment/raw_data"
)
list.files()
setwd(list.files()[3])

for (i in 1:length(list.files())) {
  if (i == 1) {
    sim_list <- vector(mode = "list", length = length(list.files()))
  }
  sim_list[[i]] <- as.matrix(fread(list.files()[i]))
}

for (i in 1:length(sim_list)) {
  if (i == 1) {
    msd_list <- vector(mode = "list", length = length(sim_list))
  }
  sl_len <- length(sim_list[[i]][, 1])
  msd_list[[i]] <- sim_list[[i]][((sl_len - 199):sl_len), ]
}


tail(msd_list[[1]])

#simMSD(msd_list, "PY_URA3_20_70_FLE")
#kurzplot("PY_URA3_20_70_FLE")

mmsdplot(
  c(
    uap,
    gup,
    "PY_URA3_00_00_FLE",
    "PY_URA3_20_70_FLE",
    "PY_URA3_20_80_FLE",
    "PY_URA3_80_90_FLE",
    "PY_URA3_90_70_FLE",
    "PY_URA3_90_75_FLE",
    "PY_URA3_90_80_FLE",
    "PY_SPB_100_100_FLE"
  )
)

#cut length of sim_list to remove buffer (300 secs or 1500 steps)

debuff_sim_list <- vector("list", length = length(sim_list))
debuff_end <- nrow(sim_list[[1]])
debuff_start <- debuff_end - 3000

for (i in 1:length(sim_list)) {
  debuff_sim_list[[i]] <- sim_list[[i]][debuff_start:debuff_end, ]
  
}

start_x <- c()
start_y <- c()
for (i in 1:length(sim_list)) {
  start_x <- c(start_x, sim_list[[i]][1, 1])
  start_y <- c(start_y, sim_list[[i]][1, 2])
}

start_frame <- data.frame(ID = seq(1,length(start_x),1),x = start_x, y = start_y, rad = sqrt(start_x^2 + start_y^2))

plot(density(start_frame$rad, from = 0, to = 1))

plot(x = start_frame$x, y = start_frame$y)


end_x <- c()
end_y <- c()
for (i in 1:length(sim_list)) {
  end_x <- c(end_x, sim_list[[i]][nrow(sim_list[[i]]), 1])
  end_y <- c(end_y, sim_list[[i]][nrow(sim_list[[i]]), 2])
}

end_frame <- data.frame(ID = seq(1,length(end_x),1),x = end_x, y = end_y, rad = sqrt(end_x^2 + end_y^2))

plot(density(end_frame$rad, from = 0, to = 1))

plot(x = end_frame$x, y = end_frame$y)


#Add time and "distance to" columns (start,previous,origin(aka radius))


TBF <- sim_list
subTBF_nrow <- length(TBF[[1]][, 1])
time_seq <- matrix(seq(0, (subTBF_nrow - 1) * 0.21, 0.21), ncol = 1)
for (i in 1:length(TBF)) {
  TBF[[i]] <- cbind(TBF[[i]], time_seq)
  colnames(TBF[[i]]) <- c("x", "y", "is_bound", "time")
}


t_int <- TBF[[1]][, "time"][2]
d2svec <- c(0)
d2pvec <- c(0)
d2ovec <- c(0)
tot_d2o <- c()
d2_list <- vector("list", length(TBF))

for (i in 1:length(TBF)) {
  float_mat <- TBF[[i]]
  d2ovec[1] <- sqrt(float_mat[1, 1] ^ 2 + float_mat[1, 2] ^ 2)
  for (j in 2:nrow(float_mat)) {
    d2svec[j] <-
      sqrt((float_mat[1, 1] - float_mat[j, 1]) ^ 2 + (float_mat[1, 2] - float_mat[j, 2]) ^
             2)
    d2pvec[j] <-
      sqrt((float_mat[(j - 1), 1] - float_mat[j, 1]) ^ 2 + (float_mat[(j - 1), 2] - float_mat[j, 2]) ^
             2)
    d2ovec[j] <- sqrt(float_mat[j, 1] ^ 2 + float_mat[j, 2] ^ 2)
  }
  float_mat <-
    cbind(float_mat,
          matrix(d2svec, ncol = 1),
          matrix(d2pvec, ncol = 1),
          matrix(d2ovec, ncol = 1))
  colnames(float_mat) <-
    c("x", "y", "is_bound", "time", "D2P", "D2S", "D2O")
  row.names(float_mat) <- c(1:nrow(float_mat))
  d2_list[[i]] <- float_mat
  tot_d2o <- c(tot_d2o, d2ovec)
}

summary(tot_d2o)

plot(density(tot_d2o, from = 0, to = 1))


#Periphery Occupancy

periph_occu_list <- vector("list", length(d2_list))
bound_zone_thickness <-
  (1 - sqrt(2 / 3)) # this is outer third of nucleus
for (i in 1:length(d2_list)) {
  float_po_logic <- d2_list[[i]][, "D2O"] > (1 - bound_zone_thickness)
  float_po_matrix <-
    matrix(seq(0, (length(float_po_logic) - 1) * .21, .21), ncol = 1)
  periph_occu_list[[i]] <-
    cbind(float_po_matrix, matrix(as.numeric(float_po_logic), ncol = 1))
}

preproc_NAs <-
  rep(NA, (length(periph_occu_list[[i]][, 1]) * length(periph_occu_list)))
po_plot_frame <-
  data.frame(Trajectory = preproc_NAs,
             time = preproc_NAs,
             is_periph = preproc_NAs)
for (i in 1:length(periph_occu_list)) {
  po_plot_frame[((i - 1) * nrow(periph_occu_list[[i]]) + 1):(i * nrow(periph_occu_list[[i]])), ] <-
    data.frame(Trajectory = i,
               time = periph_occu_list[[i]][, 1],
               is_periph = periph_occu_list[[i]][, 2])
}

#set how many steps make up one bin
size_of_bin = 50

preproc_NAs <- rep(NA, (length(po_plot_frame[, 1]) / size_of_bin))
bin_po_plot_frame <-
  data.frame(Trajectory = preproc_NAs,
             time = preproc_NAs,
             is_periph = preproc_NAs)
for (i in 1:length(bin_po_plot_frame[, 1])) {
  bin_po_plot_frame[i, 1] <-
    mean(po_plot_frame[((i - 1) * size_of_bin + 1):(i * size_of_bin), 1])
  bin_po_plot_frame[i, 2] <- po_plot_frame[((i - 1) * size_of_bin + 1), 2]
  bin_po_plot_frame[i, 3] <-
    round(mean(po_plot_frame[((i - 1) * size_of_bin + 1):(i * size_of_bin), 3]))
}


bpo_boundfrac_traj <- c()
for (i in 1:length(unique(bin_po_plot_frame[, 1]))) {
  bpo_boundfrac_traj[i] <-
    mean(bin_po_plot_frame[(bin_po_plot_frame[, 1] == unique(bin_po_plot_frame[, 1])[i]), 3])
}

ordered_bpopf <- bin_po_plot_frame
ordered_bpopf <- mutate(ordered_bpopf, rank = NA)
for (i in 1:length(bpo_boundfrac_traj)) {
  float_frame <-
    filter(bin_po_plot_frame, Trajectory == (order(bpo_boundfrac_traj)[i]))
  ordered_bpopf[((i - 1) * nrow(float_frame) + 1):(i * nrow(float_frame)), 1:3] <-
    float_frame
  ordered_bpopf[((i - 1) * nrow(float_frame) + 1):(i * nrow(float_frame)), 4] <-
    i
}

obpo_plot <- ggplot(data = ordered_bpopf) +
  geom_tile(aes(
    x = time,
    y = rank,
    color = is_periph,
    fill = is_periph
  ),
  position = "dodge") +
  theme(legend.position = "none") +
  scale_colour_viridis() +
  scale_fill_viridis() +
  comtheme

print(obpo_plot)

bpo_boundfrac_time <- c()
for (i in 1:length(unique(bin_po_plot_frame[, 2]))) {
  bpo_boundfrac_time[i] <-
    mean(bin_po_plot_frame[(bin_po_plot_frame[, 2] == unique(bin_po_plot_frame[, 2])[i]), 3])
}

plot(density(bpo_boundfrac_time))

plot(y = bpo_boundfrac_time, x = c((1:length(bpo_boundfrac_time)) * (0.21 * size_of_bin)))

bpobf_time_frame <-
  data.frame(bound_frac = bpo_boundfrac_time, time = c((1:length(bpo_boundfrac_time)) * (0.21 * size_of_bin)))

ggplot(data = bpobf_time_frame) +
  stat_smooth(aes(x = time, y = bound_frac * 100)) +
  geom_point(aes(x = time, y = bound_frac * 100), alpha = 0.5) +
  coord_cartesian(ylim = c(0, 100)) +
  comtheme +
  ylab("% of Trajectories at the periphery") +
  xlab("time (s)")


#TODO:create comparative plot for u2b8 and u9b8 FLE sims

#u2b8_bpo_bf_traj <- bpo_boundfrac_traj
#u2b8_bpo_bf_time <- bpo_boundfrac_time

#u9b8_bpo_bf_traj <- bpo_boundfrac_traj
#u9b8_bpo_bf_time <- bpo_boundfrac_time

#u9b9_bpo_bf_traj <- bpo_boundfrac_traj
#u9b9_bpo_bf_time <- bpo_boundfrac_time
#
# u0b0_bpo_bf_traj <- bpo_boundfrac_traj
# u0b0_bpo_bf_time <- bpo_boundfrac_time

comp_traj_frame <-
  data.frame(
    trajperc = c(
      u0b0_bpo_bf_traj,
      u2b8_bpo_bf_traj,
      u9b8_bpo_bf_traj,
      u9b9_bpo_bf_traj
    ),
    ID = c(
      rep("u00/b00", length(u0b0_bpo_bf_traj)),
      rep("u20/b80", length(u2b8_bpo_bf_traj)),
      rep("u90/b80", length(u9b8_bpo_bf_traj)),
      rep("u90/b90", length(u9b9_bpo_bf_traj))
    )
  )


ctjplot <- ggplot(data = comp_traj_frame) +
  geom_violin(
    aes(
      x = ID,
      y = (trajperc * 100),
      fill = ID
    ),
    draw_quantiles = c(0.25, 0.5, 0.75),
    width = 5,
    size = 1,
    position = position_dodge(width = 20)
  ) +
  geom_point(aes(x = ID, y = (trajperc * 100)),
             alpha = 0.2,
             position = position_jitter(height = 0, width = 0.1)) +
  scale_fill_viridis(discrete = T) +
  scale_color_viridis(discrete = T) +
  coord_cartesian(ylim = c(0, 100)) +
  comtheme +
  ylab("% time at the periphery per Trajectory") +
  xlab("Unbound to Bound rate / Bound to Bound rate")

print(ctjplot)

comp_time_frame <-
  data.frame(
    timeperc = c(
      u0b0_bpo_bf_time,
      u2b8_bpo_bf_time,
      u9b8_bpo_bf_time,
      u9b9_bpo_bf_time
    ),
    ID = c(
      rep("u00/b00", length(u0b0_bpo_bf_time)),
      rep("u20/b80", length(u2b8_bpo_bf_time)),
      rep("u90/b80", length(u9b8_bpo_bf_time)),
      rep("u90/b90", length(u9b9_bpo_bf_time))
    )
  )

ctmplot <- ggplot(data = comp_time_frame) +
  geom_violin(
    aes(x = ID, y = timeperc * 100, fill = ID),
    size = 1,
    width = 3,
    draw_quantiles = c(0.25, 0.5, 0.75)
  ) +
  geom_point(aes(x = ID, y = (timeperc * 100)),
             alpha = 0.05,
             position = position_jitter(height = 0, width = 0.1)) +
  scale_fill_viridis(discrete = T) +
  coord_cartesian(ylim = c(0, 100)) +
  comtheme +
  ylab("% of Trajectories at periphery per second")

print(ctmplot)

#cluster trajectory comparison using d2_list from above

#TODO: graphic for anti correlative movement

#TODO:turn d2_list into two cluster pair lists

d2_list

#Repositioning
#REPOSITIONING USES A NON_DEBUFF SIM_LIST passed through the d2_list

Repos_List <- d2_list
sim_lin()

tvv_frame <-
  rbind(tvv_frame, data.frame(velocity <-
                                tvv, ID = rep("u9b8", length(tvv))))
mdt_frame <-
  rbind(mdt_frame,
        data.frame(mean_delta_theta = global_summary_sim[, "mean_delta_theta"], ID = rep("u9b8", nrow(global_summary_sim))))
arrival_frame <-
  rbind(arrival_frame,
        data.frame(arrival_time = global_summary_sim[, "time"], ID = rep("u9b8", nrow(global_summary_sim))))

ggplot(data = tvv_frame) +
  geom_violin(
    aes(x = tvv_frame[, 2], y = tvv_frame[, 1], fill = tvv_frame[, 2]),
    size = 2,
    draw_quantiles = c(0.25, 0.5, 0.75)
  ) +
  scale_fill_viridis(discrete = T) +
  comtheme +
  xlab("Simulation") +
  ylab("velocity (um/s)")

ks.test(tvv_frame[tvv_frame[, 2] == "u2b8", 1], tvv_frame[tvv_frame[, 2] ==
                                                            "u9b8", 1])

ggplot(data = mdt_frame) +
  geom_violin(
    aes(x = mdt_frame[, 2], y = mdt_frame[, 1], fill = mdt_frame[, 2]),
    size = 2,
    draw_quantiles = c(0.25, 0.5, 0.75)
  ) +
  scale_fill_viridis(discrete = T) +
  comtheme +
  xlab("Simulation") +
  ylab("mean delta theta")


ks.test(mdt_frame[mdt_frame[, 2] == "u2b8", 1], mdt_frame[mdt_frame[, 2] ==
                                                            "u9b8", 1])

ggplot(data = arrival_frame) +
  geom_violin(
    aes(x = arrival_frame[, 2], y = arrival_frame[, 1], fill = arrival_frame[, 2]),
    size = 2,
    draw_quantiles = c(0.25, 0.5, 0.75)
  ) +
  scale_fill_viridis(discrete = T) +
  comtheme +
  xlab("Simulation") +
  ylab("Time Before First Contact")

ks.test(arrival_frame[arrival_frame[, 2] == "u2b8", 1], arrival_frame[arrival_frame[, 2] ==
                                                                        "u9b8", 1])