#dv <- diff.vector(xyz, c(29,99),  gaps.pos$f.inds)
#o <- overlap(sim116$pc, dv)

"rmsip" <-
function(pca1, pca2, subset = 10, name.prefix=NULL) {

  if (ncol(pca1$U)!=ncol(pca2$U))
    stop("unequal dimensions")

  o <- matrix(0,subset,subset)
  for ( i in 1:subset ) {
    for ( j in 1:subset ) {
      s <- pca1$U[,i] %*% pca2$U[,j]
      o[i,j] = (s)^2
    }
  }

  rmsip <- sqrt(sum(o)/subset)

  if (!is.null(name.prefix)) {
    rownames(o) <- paste(name.prefix[1], c(1:10))
    colnames(o) <- paste(name.prefix[2], c(1:10))
  }
  out <- list(matrix=round(o,3), rmsip=rmsip)
  
  return( out )
}



"diff.vector" <- function(xyz, pdb.inds, xyz.inds) {
  
  a <- xyz[pdb.inds[1],xyz.inds]
  b <- xyz[pdb.inds[2],xyz.inds]

  if (length(a)!=length(b))
    stop("unequal lengths")
  
  diff <- normalizedVector(b-a)
  return( diff )

}

"overlap" <-
  function(pca, diff, num.modes=20) {

    if (ncol(pca$U)!=length(diff))
      stop("unequal vector lengths")
      
    overlap.values <- c()
    for ( i in 1:num.modes ) {
      p <- pca$U[,i]
      
      o <- dotProduct( diff, normalizedVector(p) )
      o <- c(o*o)
      overlap.values <- c( overlap.values, c(o) )
    }
    
    cum <- cumsum(overlap.values)
    #print(colSums(as.matrix(overlap.values)))

    out <- list(v=overlap.values, c=cum) 
    
    return(out)
    
  }



"normalizedVector" <-
  function(v) {
    return( v/sqrt(dotProduct(v,v)) )
  }

"dotProduct" <-
  function(a,b) {
    o <- colSums( matrix( a*b, length(a) ) )
    return(o)
}

#print(rmsip(pc.traj4, pc.traj4))