# Install and load required packages
if (!require('TreeTools')) install.packages('TreeTools')
library('TreeTools')
if (!require('TreeSearch')) install.packages('TreeSearch')
library('TreeSearch')

# Load data into R
setwd()
Corylophidae_full <- "File S3.tnt"
raw_data <- ReadTntCharacters(Corylophidae_full)
my.phyDat <- MatrixToPhyDat(raw_data)

# Load backbone tree into R
tips <- c('Pharaxonotha', 'Hapalips', 'Holopsis', 'Clypastrea', 'Athrolips', 'Sericoderus', 'Aposericoderus', 'Corylophus', 'Teplinus', 'Orthoperus', 'Gloeosoma', 'Rypobius', 'Foadia', 'Priamima', 'Hyplathrinus', 'Ectinocephalus', 'Aenigmaticum', 'Stanus', 'Cleidostethus', 'Periptyctus', 'Pakalukodes', 'Weirus', 'Sticholotis', 'Coccidophilus', 'Bystus', 'Aphanocephalus', 'Notiophygus', 'Xenostanus')
constraint <- StringToPhyDat('0?1111?1?11111??11?1??1111?? 0?1111?1?11111??11?1??1110?? 0?1111?1?11111??11?1??0010?? 0?0000?0?00000??00?0??1100?? 0?1111?1?11111??11?1??0000?? 0?1111?1?11111??11?0??0000?? 0?0111?1?11111??11?0??0000?? 0?0000?0?00011??00?0??0000?? 0?0111?1?11100??11?0??0000?? 0?0000?0?10000??10?0??0000?? 0?0111?1?01100??01?0??0000?? 0?0001?0?00000??01?0??0000?? 0?0110?1?01100??00?0??0000?? 0?0000?1?01100??00?0??0000?? 0?0000?0?01100??00?0??0000?? 0?0110?0?00000??00?0??0000?? ', tips, FALSE)

# Search tree with inapplicable data
bestTrees <- MaximizeParsimony(my.phyDat, concavity = 12, constraint = constraint, verbosity = -1)[[1]]

bestTrees <- MaximizeParsimony(my.phyDat, tree = bestTrees, concavity = 12, constraint = constraint, ratchIter = 6L, tbrIter = 4L, finalIter = 3L, maxHits = 80L)

outgroup <- c('Pharaxonotha')
rooted.tree <- EnforceOutgroup(ape::consensus(bestTrees), outgroup)

par(mar=rep(0.25, 4), cex=0.75)
plot(rooted.tree)

# Write tree
ape::write.tree(rooted.tree, file = "Full_Constrained_K=12", append = FALSE, digits = 10, tree.names = FALSE)