library(ape)
library(phangorn)

chars <- read.nexus.data("char.nex")
chars.ph <- phyDat(chars, type="USER", levels=c("0","1","2","3","4","5","6","7","8","9") )
nachars <- read.nexus.data("nachar.nex")
nachars.ph <- phyDat(nachars, type="USER", levels=c("0","1","2","3","4","5","6","7","8","9") )
#char.nex is a simplified nexus file with absence coded data
#char.nex is a simplified nexus file with each character that has some NA states, coded exclusively with NA or A

ntaxa <- length(attributes(chars)$names)
outgroup <- attributes(chars)$names[1] 
#change the outgroup by changing the number in brackets, to see a list of all taxa use attributes(chars)$names


#Calculates TL for a tree, but then subtracts the number of A<->N transitions inferred from nachar.nex
naparsimony <- function(tree, chars, nachars) {
  parsimony(tree, chars) - parsimony(tree, nachars)
  }

mvec <- c(1,2,3,4,5,6,7,8,9,10)


#attempts is the number of rearrangements you want to try before giving up
searchTrees <- function(ntaxa, chars, chars.phy, nachars.phy, attempts){
  
  outgroup <- attributes(chars)$names[1] #change the outgroup by changing the number in brackets, to see a list of all taxa use attributes(chars)$names
  
  curTree <- rtree(ntaxa,rooted = FALSE, tip.label = attributes(chars)$names)
  curTree <- root(curTree, outgroup, node = NULL, resolve.root = TRUE)

  currentScore<-naparsimony(curTree, chars.ph, nachars.ph)
  targetScore<-currentScore

  for (arrangeAttempt in 1:attempts) {
        candidateTree<-rNNI(curTree, moves = sample(mvec,1), n = 1)
	  candidateScore <- naparsimony(candidateTree,chars.ph, nachars.ph)
	  if (candidateScore < targetScore) {
	  	curTree<-candidateTree
		targetScore<-candidateScore
	  }
  }
resultTrees <- curTree
}



curTree <- rtree(ntaxa,rooted = FALSE, tip.label = attributes(chars)$names)
curTree <- root(curTree, outgroup, node = NULL, resolve.root = TRUE)
resultTrees<- write.tree(curTree, file = "resultTrees", append=FALSE, digits = 3, tree.names = FALSE)
writeScore <- naparsimony(curTree,chars.ph, nachars.ph)
#writeScore <- 21 #Number I feel good about through other evidence/analyses

arrangements<- 200
maxHits<- 7
numHit <- 1

while(numHit < maxHits){
  resultTree <- searchTrees(ntaxa, chars, chars.phy, nachars.phy, arrangements)
  resultScore <- naparsimony(resultTree, chars.ph, nachars.ph)
	if(resultScore < writeScore){
         writeScore<-resultScore
         resultTrees<- write.tree(resultTree, file = "resultTrees", append=FALSE, digits = 3, tree.names = FALSE)
         numHit<-1
         }
      if(resultScore==writeScore){
         resultTrees<- write.tree(resultTree, file = "resultTrees", append=TRUE, digits = 3, tree.names = FALSE)
         numHit <- numHit + 1
         }
}
print(writeScore)
resultTrees<- read.tree("resultTrees")
resultTrees<- unique.multiPhylo(resultTrees)
write.tree(resultTree, file = "resultTrees", append = FALSE, digits = 3, tree.names = FALSE)