library(bio3d)
library(Hmisc)
source("corr_mat_fun.R")
##source("~/Documents/bio3d_new_funs/plot.cij.R")

pdb.A <- read.pdb("1XCK_A.pdb")
sse.A <- dssp(pdb.A)

pdb <- read.pdb("1XCK_chainALM_noWAT_noH.pdb")
sse <- dssp(pdb)

vlines=c(525,1048)
hlines=c(525,1048)


helix.labels = capitalize(letters[1:length(sse.A$helix$start)])
helix.labels=rep(helix.labels, 3)

helix.labels.at=(sse$helix$start+sse$helix$end)/2
axis.at=sse$helix$start[1:54]
axis.labels=rep(sse.A$helix$start,3)

x=525:1572
y=1:524
dim=c(length(x),length(y))
dat <- array(0, dim=c(dim[1], dim[2], 12))
apo <- matrix(0,dim[1],dim[2])
j=1
for ( i in c(1,2,3,4,5,7) ) {
  A <- scan(paste("~/gulrotkake/groel_md/1XCK/116_1XCK_apo/results/correlation_trimer/corr_30-50ns_CA_chain_",i,sep=""))
  n <- sqrt(length(A))
  A <- matrix( A, n, n, byrow = TRUE)
  dat[, ,j]=A[x,y]
  apo=apo+A[x,y]
  j=j+1
}
apo=round(apo/6,2)
##write.table(apo, file="corrmat_apo_avg.dat", row.names=F, col.names=F, quote=F)

holo <- matrix(0,dim[1],dim[2])
##holo <- matrix(0,1585,1585)
for ( i in c(1,2,3,4,5,7) ) {
  A <- scan(paste("~/gulrotkake/groel_md/1XCK/129_1XCK_MGATP/results/correlation_trimer/corr_30-50ns_CA_chain_",i,sep=""))
  n <- sqrt(length(A))
  print(n)
  A <- matrix( A, n, n, byrow = TRUE)
  dat[, ,j]=A[x,y]
  holo=holo+A[x,y]
  j=j+1
}
holo=round(holo/6,2)
##write.table(holo, file="corrmat_holo_avg.dat", row.names=F, col.names=F, quote=F)

diff=apo-holo
diff=filter.corrmat(diff, threshold=-1)

v <- t.list
##v[upper.tri(v, diag=FALSE)] <- NA
t<-which(v<0.025, arr.ind=T)
sig.mat <- matrix(0,dim[1],dim[2])
sig.mat[t]=1
g <- which(sig.mat!=0, arr.ind=T)



png(paste("diff.png",sep=""), width=2000, height=1000,
    pointsize=24)
plot.corrmat(sig.mat, sse,
##plot.corrmat(diff, sse,             
             helix.labels=helix.labels, helix.labels.at=helix.labels.at,
             axis.at=axis.at, axis.labels=axis.labels,
             vlines=vlines, hlines=hlines, main="Cross correlation")
dev.off()

A=filter.corrmat(apo, threshold=.0)
png(paste("apo_avg.png",sep=""), width=2000, height=1000,
    pointsize=24)
plot.corrmat(A, sse,
             helix.labels=helix.labels, helix.labels.at=helix.labels.at,
             axis.at=axis.at, axis.labels=axis.labels,
             vlines=vlines, hlines=hlines, main="Cross correlation")
dev.off()

B=filter.corrmat(holo, threshold=.0)
png(paste("holo_avg.png",sep=""), width=2000, height=1000,
    pointsize=24)
plot.corrmat(B, sse,
             helix.labels=helix.labels, helix.labels.at=helix.labels.at,
             axis.at=axis.at, axis.labels=axis.labels,
             vlines=vlines, hlines=hlines, main="Cross correlation")
dev.off()

   png(paste("corrMat.apo_",i,".png",sep=""), width=1500, height=1500,
      pointsize=24)
  plot.corrmat(A, sse,
               helix.labels=helix.labels, helix.labels.at=helix.labels.at,
               axis.at=axis.at, axis.labels=axis.labels,
               vlines=vlines, hlines=hlines, main="Cross correlation")
  dev.off()
png(paste("corrMat.holo_",i,".png",sep=""), width=1500, height=1500,
      pointsize=24)
  plot.corrmat(A, sse,
               helix.labels=helix.labels, helix.labels.at=helix.labels.at,
               axis.at=axis.at, axis.labels=axis.labels,
               vlines=vlines, hlines=hlines, main="Cross correlation")
  dev.off()




library(multicore)
threads=4
ptm <- proc.time()
rnames=rep(1:threads, each=(nrow(dat)/threads))
print(rnames)

jobs <- list()
for ( k in 1:threads ) {
  r.inds=which(rnames==k)
  mat=dat[r.inds,,]
  q=parallel(t.test.mat(mat))
  jobs[[k]]=q
}

res=collect(jobs, wait=TRUE)
print(proc.time() - ptm)

t.list <- c()
for ( job in res ) {
  t.list=rbind(t.list, job)
}


t.test.mat <- function(mat) {
  n <- dim(mat)[1]
  m <- dim(mat)[2]
  t.list <- matrix(0,n,m)
  
  for ( i in 1:n ) {
    for ( j in 1:m ) {
      a <- mat[i,j,1:6]
      b <- mat[i,j,7:12]

      if (!identical(a,b)) {
        t <- t.test(a,b)
        t.list[i,j]=t$p.value
      }
      else
        t.list[i,j]=1
    }
  }
  return(t.list)
}


filter.corrmat <- function(A, threshold=0.2) {

  h <- which(A>threshold)
  l <- which(A<(threshold*(-1)))
  hmm=A*0
  hmm[h]=A[h]
  hmm[l]=A[l]
  return(hmm)

}