"trim" <- function(s) {                                                          
        s <- sub("^ +", "", s)                                                     
        s <- sub(" +$", "", s)                                                     
        s[(s == "")] <- NA                                                         
        s                                                                          
      }    





"plot.zoom" <- function(apo, holo, inds, n, seq, ylim=NULL, ylab="Energy (kcal/mol)", mtext.line=0, ...) {
  s.inds <- ((inds-1)%%524)+1
  hmm <- matrix(c(apo$total[inds], holo$total[inds]), nrow=2, byrow=T)
  print(head(hmm))
  mp <- barplot(hmm, ylim=ylim, ylab=ylab, beside=T, ...)

  stderr=(apo$totals.stds[inds]/sqrt(n))
  errbar(mp[1,], apo$total[inds], apo$total[inds]+stderr, apo$total[inds]-stderr, add=T, cex=0.5, cap=.0075)

  stderr=(holo$totals.stds[inds]/sqrt(n))
  errbar(mp[2,], holo$total[inds], holo$total[inds]+stderr, holo$total[inds]-stderr, add=T, cex=0.65, cap=.0075)
  mtext(1, at=colMeans(mp), text=paste(seq[s.inds], s.inds+1), las=2, line=mtext.line, ...)

  return(mp)
}


"plot.zoom.three" <- function(apo, holo, xtra, inds, n, seq, ylim=NULL, ylab="Energy (kcal/mol)", mtext.line=0, ...) {
  s.inds <- ((inds-1)%%524)+1
  hmm <- matrix(c(apo$total[inds], holo$total[inds], xtra$total[inds]), nrow=3, byrow=T)
  print(head(hmm))
  mp <- barplot(hmm, ylim=ylim, ylab=ylab, beside=T, ...)

  stderr=(apo$totals.stds[inds]/sqrt(n))
  errbar(mp[1,], apo$total[inds], apo$total[inds]+(stderr*1.96), apo$total[inds]-(stderr*1.96), add=T, cex=0.5, cap=.0075)

  stderr=(holo$totals.stds[inds]/sqrt(n))
  errbar(mp[2,], holo$total[inds], holo$total[inds]+(stderr*1.96), holo$total[inds]-(stderr*1.96), add=T, cex=0.65, cap=.0075)

  stderr=(xtra$totals.stds[inds]/sqrt(n))
  errbar(mp[3,], xtra$total[inds], xtra$total[inds]+(stderr*1.96), xtra$total[inds]-(stderr*1.96), add=T, cex=0.65, cap=.0075)
  
  mtext(1, at=colMeans(mp), text=paste(seq[s.inds], s.inds+1), las=2, line=mtext.line, ...)

  return(mp)
}


"plot.zoom.four" <- function(apo, holo, xtra.a, xtra.b, inds, n, seq,
                             ylim=NULL, ylab="Energy (kcal/mol)", mtext.line=0, ...) {
  s.inds <- ((inds-1)%%524)+1
  hmm <- matrix(c(apo$total[inds], holo$total[inds], xtra.a$total[inds], xtra.b$total[inds]), nrow=4, byrow=T)
  print(head(hmm))
  mp <- barplot(hmm, ylim=ylim, ylab=ylab, beside=T, ...)

  stderr=(apo$totals.stds[inds]/sqrt(n))
  errbar(mp[1,], apo$total[inds], apo$total[inds]+stderr, apo$total[inds]-stderr, add=T, cex=0.5, cap=.0075)

  stderr=(holo$totals.stds[inds]/sqrt(n))
  errbar(mp[2,], holo$total[inds], holo$total[inds]+stderr, holo$total[inds]-stderr, add=T, cex=0.65, cap=.0075)

  stderr=(xtra.a$totals.stds[inds]/sqrt(n))
  errbar(mp[3,], xtra.a$total[inds], xtra.a$total[inds]+stderr, xtra.a$total[inds]-stderr, add=T, cex=0.65, cap=.0075)

  stderr=(xtra.b$totals.stds[inds]/sqrt(n))
  errbar(mp[4,], xtra.b$total[inds], xtra.b$total[inds]+stderr, xtra.b$total[inds]-stderr, add=T, cex=0.65, cap=.0075)
  
  mtext(1, at=colMeans(mp), text=paste(seq[s.inds], s.inds+1), las=2, line=mtext.line, ...)

  return(mp)
}





"read.mmpbsa.decomp" <- function(filename) {
  
  zz <- readLines(filename, n=-1)

  head <- head(zz, n=3)
  i <- 4

  ## check what information is given
  backbone.inds <- grep("Backbone", zz)
  sidechain.inds <- grep("Sidechain", zz)
  total.inds <- grep("Total", zz)

  com.inds <- grep("Complex", zz)
  rec.inds <- grep("Receptor", zz)
  lig.inds <- grep("Ligand", zz)
  delta.inds <- grep("DELTAS", zz)

  seps <- grep("-------------", zz)
  seps=c(seps, length(zz))
  empts <- grep("         ", zz)


  m <- 1 ## lets keep track of how far in the seps-list we have gotten
  complex <- NULL
  receptor <- NULL
  ligand <- NULL
  deltas <- NULL
  if ( length(com.inds) > 0 ) {

    # if complex is written, then Total is written
    # thus, read until 2nd seperator
    inds <- c( (com.inds[1]+4) : (seps[m+1]-2) )
    complex$total = zz[inds]
    m=m+1
    
    if ( length(backbone.inds) > 0 ) {
      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      complex$sidechain = zz[inds]
      m=m+1

      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      complex$backbone = zz[inds]
      m=m+1
    }
  
    # in fact, lets assume that receptor is given
    inds <- c( (rec.inds[1]+4) : (seps[m+1]-2) )
   
    receptor$total = zz[inds]
    m=m+1

    if ( length(backbone.inds) > 0 ) {
      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      receptor$sidechain = zz[inds]
      m=m+1

      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      receptor$backbone = zz[inds]
      m=m+1
    }

    ## and then ligand should be given
    inds <- c( (lig.inds[1]+4) : (seps[m+1]-2) )
    ligand$total = zz[inds]
    m=m+1

    if ( length(backbone.inds) > 0 ) {
      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      ligand$sidechain = zz[inds]
      m=m+1

      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      ligand$backbone = zz[inds]
      m=m+1
    }

  }

  if ( length(delta.inds) > 0 ) {
    if ( m==1 )
      m=m+1

    inds <- c( (delta.inds[1]+4) : (seps[m+1]-2) )
    deltas$total = zz[inds]
    m=m+1
    
    if ( length(backbone.inds) > 0 ) {
      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      deltas$sidechain = zz[inds]
      m=m+1

      inds <- c( (seps[m]+1) : (seps[m+1]-2) )
      deltas$backbone = zz[inds]
      m=m+1
    }
  
  }

  complex$total = parse.all(complex$total)
  complex$sidechain = parse.all(complex$sidechain)
  complex$backbone = parse.all(complex$backbone)
  
  receptor$total = parse.all(receptor$total)
  receptor$sidechain = parse.all(receptor$sidechain)
  receptor$backbone = parse.all(receptor$backbone)

  ligand$total = parse.all(ligand$total)
  ligand$sidechain = parse.all(ligand$sidechain)
  ligand$backbone = parse.all(ligand$backbone)

  deltas$total = parse.all(deltas$total)
  deltas$sidechain = parse.all(deltas$sidechain)
  deltas$backbone = parse.all(deltas$backbone)

  deltas2=calc.deltas(complex, receptor, ligand)
  #deltas2 <- NULL
  
  out <- list( complex=complex, receptor=receptor, ligand=ligand, deltas=deltas, deltas2=deltas2)
  return(out)

}


"calc.deltas" <- function(complex, receptor, ligand) {

  deltas2 <- NULL
  c = matrix(as.numeric( complex$total[,c(3,5,7,9,11,13)] ), ncol=6)
  c.std = matrix(as.numeric( complex$total[,c(4,6,8,10,12,14)] ), ncol=6)
  r = matrix(as.numeric( receptor$total[,c(3,5,7,9,11,13)] ), ncol=6)
  r.std = matrix(as.numeric( receptor$total[,c(4,6,8,10,12,14)] ), ncol=6)
  l = matrix(as.numeric( ligand$total[,c(3,5,7,9,11,13)] ), ncol=6)
  l.std = matrix(as.numeric( ligand$total[,c(4,6,8,10,12,14)] ), ncol=6)


  print(paste(dim(c), dim(r), dim(l)))
  
  
  means=round(c-rbind(r,l),3)
  stds=round(sqrt( c.std**2 + rbind(r.std, l.std)**2 ), 3)

  deltas2$total = cbind(means[,1], stds[,1],
    means[,2], stds[,2],
    means[,3], stds[,3],
    means[,4], stds[,4],
    means[,5], stds[,5],
    means[,6], stds[,6]
    )
  
  colnames(deltas2$total) <- c("int",  "std", "vdw", "std", "ele", "std", "polsol", "std", "nonpolsol", "std", "total", "std")

  if ( !is.null(complex$sidechain) ) {
    c = matrix(as.numeric( complex$sidechain[,c(3,5,7,9,11,13)] ), ncol=6)
    c.std = matrix(as.numeric( complex$sidechain[,c(4,6,8,10,12,14)] ), ncol=6)
    r = matrix(as.numeric( receptor$sidechain[,c(3,5,7,9,11,13)] ), ncol=6)
    r.std = matrix(as.numeric( receptor$sidechain[,c(4,6,8,10,12,14)] ), ncol=6)
    l = matrix(as.numeric( ligand$sidechain[,c(3,5,7,9,11,13)] ), ncol=6)
    l.std = matrix(as.numeric( ligand$sidechain[,c(4,6,8,10,12,14)] ), ncol=6)

    print(paste(dim(c), dim(r), dim(l)))
    
    means=round(c-rbind(r,l),3)
    stds=round(sqrt( c.std**2 + rbind(r.std, l.std)**2 ), 3)
      
    deltas2$sidechain = cbind(means[,1], stds[,1],
      means[,2], stds[,2],
      means[,3], stds[,3],
      means[,4], stds[,4],
      means[,5], stds[,5],
      means[,6], stds[,6]
      )
    
    c = matrix(as.numeric( complex$backbone[,c(3,5,7,9,11,13)] ), ncol=6)
    c.std = matrix(as.numeric( complex$backbone[,c(4,6,8,10,12,14)] ), ncol=6)
    r = matrix(as.numeric( receptor$backbone[,c(3,5,7,9,11,13)] ), ncol=6)
    r.std = matrix(as.numeric( receptor$backbone[,c(4,6,8,10,12,14)] ), ncol=6)
    l = matrix(as.numeric( ligand$backbone[,c(3,5,7,9,11,13)] ), ncol=6)
    l.std = matrix(as.numeric( ligand$backbone[,c(4,6,8,10,12,14)] ), ncol=6)

    print(paste(dim(c), dim(r), dim(l)))
    
    means=round(c-rbind(r,l),3)
    stds=round(sqrt( c.std**2 + rbind(r.std, l.std)**2 ), 3)
    
    deltas2$backbone = cbind(means[,1], stds[,1],
      means[,2], stds[,2],
      means[,3], stds[,3],
      means[,4], stds[,4],
      means[,5], stds[,5],
      means[,6], stds[,6]
      )

    colnames(deltas2$sidechain) <- c("int",  "std", "vdw", "std", "ele", "std", "polsol", "std", "nonpolsol", "std", "total", "std")
    colnames(deltas2$backbone) <- c("int",  "std", "vdw", "std", "ele", "std", "polsol", "std", "nonpolsol", "std", "total", "std")
  }

  return(deltas2)

}




"parse.all" <- function(obj) {

  if ( length(obj) == 0 )
    return (NULL)

  hmm <- NULL
  for ( l in 1:length(obj) ) {
    m <- parse.energies(obj[l])

    if ( !is.null(m) )
      hmm=rbind(hmm, m)
  }

  if ( !is.null(hmm) ) {
    if ( ncol(hmm)==14 )
      colnames(hmm) <- c("aa", "resno", "int", "std", "vdw", "std", "ele", "std", "polsol", "std", "nonpolsol", "std", "total", "std")
    else if ( ncol(hmm)==15 )
      colnames(hmm) <- c("aa", "resno", "location", "int",  "std", "vdw", "std", "ele", "std", "polsol", "std", "nonpolsol", "std", "total", "std")
  }
    return(hmm)
}



"parse.energies" <- function(line) {

  split=strsplit(line, " | ", fixed=T)
  d <- unlist(split)

  ## if 7 elements, this is not deltas
  j <- 0
  if ( length( d ) > 6 ) {
    res=trim(d[j+1])
    resnum=trim( substr(res, 4, nchar(res)) )
    res=trim( substr(res, 1, 3) )

    loc=NULL
    if ( length( d ) == 8 ) {
      j=j+1
      loc=trim(d[j+1])
    }
            
    int=strsplit( trim(d[j+2]), "+/-", fixed=T )
    int=trim(unlist(int))
    vdw=strsplit( trim(d[j+3]), "+/-", fixed=T )
    vdw=trim(unlist(vdw))
    ele=strsplit( trim(d[j+4]), "+/-", fixed=T )
    ele=trim(unlist(ele))
    pol=strsplit( trim(d[j+5]), "+/-", fixed=T )
    pol=trim(unlist(pol))
    nonpol=strsplit( trim(d[j+6]), "+/-", fixed=T )
    nonpol=trim(unlist(nonpol))
    tot=strsplit( trim(d[j+7]), "+/-", fixed=T )
    tot=trim(unlist(tot))

    if ( !is.null(loc) )
      m <- matrix( c(res, resnum, loc, int, vdw, ele, pol, nonpol, tot), nrow=1 )
    else
      m <- matrix( c(res, resnum, int, vdw, ele, pol, nonpol, tot), nrow=1 )

    return(m)
  }
  else
    return(NULL)


}
  
"read.mmpbsa.decomp.old" <- function(filename) {


  ## read in residue contributions
  i=0
  k=0
  res.cont <- NULL
  while ( i < length(zz) ) {
    i=i+1
    line = zz[i]

    if ( length( grep("------", line) ) > 0 ) {
      k=k+1
      j=0
    }

    split=strsplit(line, " | ", fixed=T)
    d <- unlist(split)
    
    if ( length( d ) == 7 && d[1] != "Residue" ){

      res=trim(d[1])
      tot=strsplit(d[length(d)], "+/-", fixed=T)
      tot=unlist(tot)
      tot=as.numeric(tot[1])

      if ( j == 0 ) {
        res.cont[[k]] <- matrix(c(res, tot), ncol=2)
        #res.cont[[k]] = as.data.frame(res.cont[[k]])
        j=1
      }
      else
        res.cont[[k]] = rbind(res.cont[[k]], c(res, tot))
    }
   
  }
  #print(res.cont[[k]]$total)

  #print(total.inds)

  complex=list(
    total=res.cont[[1]],
    sidechain=res.cont[[2]],
    backbone=res.cont[[3]] )

  receptor=list(
    total=res.cont[[4]],
    sidechain=res.cont[[5]],
    backbone=res.cont[[6]] )

  ligand=list(
    total=res.cont[[7]],
    sidechain=res.cont[[8]],
    backbone=res.cont[[9]] )


  rec.lig.tot <- c(receptor$total[,2], ligand$total[,2])
  deltas <- NULL
  deltas$total = as.numeric(complex$total[,2]) - as.numeric(rec.lig.tot)

  #if ( length(delta.inds)>0 ) {
  #  deltas=list(
  #    total=res.cont[[10]],
  #    sidechain=res.cont[[11]],
  #    backbone=res.cont[[12]] )
  #}
  #else {
  #total.rec = as.numeric(complex$total[1:524,2]) - as.numeric(receptor$total[,2])
  #total.lig = as.numeric(complex$total[525:1048,2]) - as.numeric(ligand$total[,2])
  #  sidechain=as.numeric(complex$sidechain[,2])-as.numeric(receptor$sidechainl[,2])
  #  backbone=as.numeric(complex$backbone[,2])-as.numeric(receptor$backbone[,2])
  #deltas=list(total=total)
    #, sidechain=sidechain, backbone=backbone)
  #}

  out <- list(complex=complex, receptor=receptor, ligand=ligand, deltas=deltas)
  return(out)
}



"calc.means" <- function(z, n) {
  tots <- NULL
  vdw <- NULL
  ele <- NULL
  polsol <- NULL
  sas <- NULL
  stds <- NULL
  for ( j in 1:n ) {
    k <- z[[j]]
    tots=c(tots, k$deltas2$total[,"total"])
    vdw=c(vdw, k$deltas2$total[,"vdw"])
    ele=c(ele, k$deltas2$total[,"ele"])
    polsol=c(polsol, k$deltas2$total[,"polsol"])
    sas=c(sas, k$deltas2$total[,"nonpolsol"])
    stds=c(stds, k$deltas2$total[,12])
  }

  tots=matrix(as.numeric(tots), ncol=n)
  tot.means=rowMeans(tots)

  vdw=matrix(as.numeric(vdw), ncol=n)
  vdw=rowMeans(vdw)

  ele=matrix(as.numeric(ele), ncol=n)
  ele=rowMeans(ele)

  polsol=matrix(as.numeric(polsol), ncol=n)
  polsol=rowMeans(polsol)

  ele2 <- cbind(ele, polsol)
  ele2 = rowSums(ele2)
 
  sas=matrix(as.numeric(sas), ncol=n)
  sas=rowMeans(sas)  

  tots=matrix(as.numeric(tots), ncol=n)
  tot.means=rowMeans(tots)
    
  stds=matrix(as.numeric(stds), ncol=n)
  stds=stds**2
  stds=rowSums(stds)
  stds=sqrt(stds)

  out <- list(total=tot.means, totals.stds=stds, vdw=vdw, ele=ele2, sas=sas )
  return(out)
}