
args=commandArgs(TRUE)

# Packages to load #######
library(gtools)
library(data.table)

indir <- args[1]#"brain-hprd/dyn-net-processed"#args[1]
outdir <- args[2]#"brain-hprd/dyn-net-norm-diff"#args[2]
epsilon <- as.numeric(as.character(args[3]))#0.01#as.numeric(as.character(args[3]))
jump <- as.numeric(as.character(args[4]))#2#as.numeric(as.character(args[4]))
start <- as.numeric(as.character(args[5]))+1

# Min. amd max of edge weights among all network snapshots
adirs <- mixedsort(list.files(indir, full.names=TRUE))

temp <- fread(adirs[1], sep='\t')
MIN <- min(temp[[3]])
MAX <- max(temp[[3]])

for(k in 1:length(adirs)){

  temp <- fread(adirs[k], sep='\t')
  MIN <- min(MIN, min(temp[[3]]))
  MAX <- max(MAX, max(temp[[3]]))

}

ll <- length(adirs)-1

if( jump > ll ) stop('higher-order relationship should be at most 1 less than the maximum number of snapshots')

# normalize weights in every network snapshot 
loop <- length(adirs)-jump

for(k in 1:loop){

  j <- k+jump

  temp1 <- fread(adirs[k], sep='\t')
  # scale the edge weights between 0 and 1
  temp1$V3 <- ((temp1$V3-MIN)/(MAX-MIN))*(1-epsilon)+epsilon


  temp2 <- fread(adirs[j], sep='\t')
  # scale the edge weights between epsilon and 1
  temp2$V3 <- ((temp2$V3-MIN)/(MAX-MIN))*(1-epsilon)+epsilon


  temp11 <- data.table(temp1, key = c("V1", "V2"))
  temp22 <- data.table(temp2, key = c("V1", "V2"))

  temp33 <- data.frame(temp11$V3, temp22$V3)
  temp11$V4 <- apply(temp33, 1, FUN=sum)
  # save percentage change relative to the higher weight
  temp22$V3 <- ((temp22$V3-temp11$V3)/temp11$V4)*100

  #write.table(temp22, paste0(outdir, '/', strsplit(basename(adirs[k]),'[.]')[[1]][1],'.txt'),  col.names=FALSE, row.names=FALSE, quote=FALSE, sep='\t')
  write.table(temp22, paste0(outdir, '/',start,'.txt'),  col.names=FALSE, row.names=FALSE, quote=FALSE, sep='\t')

  start <- start + 1

}

