## function to pack a vcov back into a tree
phylo.vcv <- function(v,reorder=TRUE) {
  ## no polytomies allowed
  va <- v
  tipnames <- rownames(v)
  ntip <- nrow(v)
  dimnames(v) <- list(as.character(1:ntip),
                      as.character(1:ntip))
  diag(va) <- 0
  edgemat <- matrix(ncol=2,nrow=0)
  ## termlens <- diag(v)-colSums(va)
  edgelens <- numeric(0)
  maxnode <- ntip
  while (nrow(v)>1) {
    mva <- max(va)  ## find pair with max shared evolution
    nextpr <- if (nrow(v)==2) c(1,2) else which(va==mva,arr.ind=TRUE)[1,]
    maxnode <- maxnode+1  ## new node
    ## points to both of current identified nodes
    ##   (indexed by names)
    edgemat <- rbind(edgemat,
                     c(maxnode,as.numeric(rownames(v)[nextpr[1]])),
                     c(maxnode,as.numeric(rownames(v)[nextpr[2]])))
    ## descending edges are amount of *unshared* evolution
    edgelens <- c(edgelens,
                  
                  diag(v)[nextpr]-mva)
    ## this clade has total evolution = shared evolution
    diag(v)[nextpr] <- mva
    ## assign new node name
    rownames(v)[nextpr[1]] <- colnames(v)[nextpr[1]] <- maxnode
    ## drop rows/cols from matrix
    v <- v[-nextpr[2],-nextpr[2],drop=FALSE]
    va <- va[-nextpr[2],-nextpr[2],drop=FALSE]
  }
  ## switch order of node numbers to put root in the right place:
  ##  much plotting code seems to assume root = node # (ntips+1)
  nn <- nrow(edgemat)
  nnode <- nn-ntip+1
  newedge <- edgemat
  for (i in 1:nnode) {
    newedge[edgemat==(ntip+i)] <- nn-i+2
  }
  names(edgelens) <- NULL
  r <- list(edge=newedge,
        tip.label=tipnames,
        edge.length=edgelens,
        Nnode =  nnode)
   class(r) <- "phylo"
  if (!reorder) r else reorder(r) ## for plotting safety
 }