
args = commandArgs(trailingOnly=TRUE)
if (length(args) != 2) {
  cat("Create a table with afmCC and (normalzed afmCC x total Energy) for each kappa\n")
  cat("Usage:\n")
  cat("       Rscript get_energyTally_3D.R kappa_file frequency_factor\n")
  cat("       kappa_file         file with list of kappa values\n")
  cat("       frequency_factor   interger to print every Nth value from the kappa list\n\n")
  cat("       Output files: energy_tally_3D.txt\n")
  quit()
}

#install.packages("ggplot2")
#library(gdata)
library(dplyr)
library(stringr)
library(bio3d)
#install.packages("gridExtra")
library(gridExtra)


zfile = args[1]
conn = file(zfile,open="r")
zvalues = readLines(conn)
input.folders = paste0("Zscale_", zvalues)
ind = zvalues
every_n = args[2]
#input.folders = c(unlist(strsplit(str_sub(args[1], 2, -2),",")))
#ind <- c(unlist(strsplit(str_sub(args[2], 2, -2),",")))

#print(length(ind))
#print(length(args)-2)
#setwd("/Users/bhandariyr/Documents/projects/afm/cobalamin/miscen/scripts/stats/contact10_5")

## initialize data frames for each energy term
energy.list = list()

## Loop through the given folders
# folders = list("output1_0","output1_1", "output1_2","output1_3","output1_4","output1_5",
#                "output1_6","output1_7","output1_8","output1_9","output1_10")
count = 1 ## counter for list
for (i in 1:(length(input.folders))){

#print(as.integer(i) %% as.integer(every_n))
if ((as.integer(i) %% as.integer(every_n)) == 0){  
ext1 = paste(input.folders[i],"output", "*.ts", sep = "/")
file = Sys.glob(ext1)
data.ts.lines = readLines(file)[grep("#all", readLines(file))]
data.ts =  data.ts.lines %>%
  strsplit("\\s+") %>%
  do.call(rbind, .) %>%
  data.frame()
dim(data.ts)
colnames(data.ts) = unlist(strsplit(readLines(file)[grep("#unit", readLines(file))], "\\s+"))
#head(data.ts)
colnames(data.ts)[12] = "stack"
# View(data.ts)
data.ts$frame = as.numeric(rownames(data.ts)) - 1

data.ts.sort <- data.ts %>% select("frame", everything())
## subset to remove frame zero
data.ts.sort = data.ts.sort[2:nrow(data.ts.sort),]
# View(data.ts.sort)
## print the energy corresponding to top X percent afmcc sturctures
#top_percent = as.numeric(args[2])

## calculate correlation between afmcc and nromCCxEtot
afmcc = as.numeric(data.ts.sort$afmcc) #Afmcc score
toten = as.numeric(data.ts.sort$etot)  #total energy

cc_norm = (afmcc - min(afmcc))/(max(afmcc) - min(afmcc))
ccNormxEtot = cc_norm*toten
col1 = paste0("afmCC_", input.folders[i])
col2 = paste0("normCCxEtot_", input.folders[i])
df_energy = data.frame(cc = afmcc, ce = ccNormxEtot)
colnames(df_energy) = c(col1, col2) 
energy.list[[count]] = df_energy
count = count + 1
#ce.list[[i]] = ccNormxEtot 
} ## end of if loop
} ## end of for loop
file = paste0("energy_tally3D_every", every_n, ".txt")
df_energy = do.call(cbind, energy.list)
write.table(format(data.frame(df_energy)), file = file, sep = "\t", row.names = FALSE, quote = FALSE)

## Reading data table
## data = read.delim(file, header = TRUE, sep = "\t", dec = ".")
