add test, nicer output

R/pml_generics.R

@@ -144,32 +144,30 @@ print.pml <- function(x, ...) {
 #' unlink(c("woodmouse_pml.txt", "woodmouse_tree.nwk", "woodmouse.rds"))
 #' @importFrom utils citation
 #' @export
-write.pml <- function(x, file="pml", save_rds=TRUE, ...){
+write.pml <- function(x, file="pml", save_rds=FALSE, ...){
   digits <- -1
   if (hasArg("digits")) digits <- list(...)$digits
   write.tree(x$tree, file=paste0(file, "_tree.nwk"))
   if(save_rds) saveRDS(x, file=paste0(file, ".rds"))
-  else write.phyDat(x$data, file=paste0(file, "_align.fasta"), format="fasta")
+  write.phyDat(x$data, file=paste0(file, "_align.fasta"), format="fasta")
   if(!is.null(x$bs)) write.nexus(x$bs, file=paste0(file, "_bs.nex"),
                                  digits=digits)
   sink(paste0(file, ".txt"))
   cat("phangorn", packageDescription("phangorn", fields = "Version"), "\n\n")
   print(x)
-  cat("\n\n")
-  cat("You can (re-)create the pml object using:\n\n")
+  cat("\n\nThe following lines (re-)creates the pml object up to numerical inaccuracies:\n\n")
+  call <- x$call
+  call$data <- quote(align)
+  call$tree <- quote(tree)
+  cat("tree <- read.tree(\"", file, "_tree.nwk\")\n", sep="")
+  cat("align <- read.phyDat(\"", file, "_align.fasta\", format=\"fasta\")",
+      sep="")
+  cat( "\nfit <- ")
+  print(call)
   if(save_rds){
+    cat("\nAnd the following reproduces the exact pml object:\n\n")
     cat("fit <- readRDS(\"", file,".rds\")", sep="")
   }
-  else {
-    call <- x$call
-    call$data <- quote(align)
-    call$tree <- quote(tree)
-    cat("tree <- read.tree(\"", file, "_tree.nwk\")\n", sep="")
-    cat("align <- read.phyDat(\"", file, "_align.fasta\", format=\"fasta\")",
-        sep="")
-    cat( "\nfit <- ")
-    print(call)
-  }
   cat("\n\nREFERENCES\n\n")
   cat("To cite phangorn please use:\n\n")
   print(citation("phangorn") [[1]], style="text")

inst/tinytest/test_pml_generics.R

@@ -1,5 +1,5 @@
 data(woodmouse)
-
+set.seed(42)
 fit <- pml_bb(woodmouse, "HKY+I", rearrangement = "NNI",
               control=pml.control(trace=0))
 
@@ -14,18 +14,16 @@ write.pml(fit, save_rds = TRUE)
 fit2 <- readRDS("pml.rds")
 expect_equal(fit, fit2)
 
+# Add more checks (best expect_equal(fit, fit3))
+# order of tip label changes (and of data object)
+# attribute method might missing (when optim.pml was not called)
+tree <- read.tree("pml_tree.nwk")
+tree <- relabel(tree, fit$tree$tip.label)
+align <- read.phyDat("pml_align.fasta", format="fasta")
+fit3 <- pml(tree = tree, data = align, k = 1L, site.rate = "gamma", ASC = FALSE,
+           bf = c(0.306541405706333, 0.261308281141267, 0.126026443980515,
+                  0.306123869171886), Q = c(1, 23.0401029286372, 1, 1, 23.0401029286372,
+                                            1), inv = 0.850681184138605, model = "HKY")
+expect_equal(logLik(fit), logLik(fit3))
 
-
-#write.pml(fit, save_rds = FALSE)
-
-#tree <- read.tree("pml_tree.nwk")
-#align <- read.phyDat("pml_align.fasta", format="fasta")
-#fit3 <- pml(tree = tree, data = align, k = 1L, site.rate = "gamma", ASC = FALSE,
-#           bf = c(0.306541405706333, 0.261308281141267, 0.126026443980515,
-#                  0.306123869171886),
-#           Q = c(1, 23.0401288315219, 1, 1, 23.0401288315219, 1),
-#           inv = 0.850681761772363, model = "HKY")
-
-#fit$method <- NULL
-#expect_equal(fit, fit3)
 unlink(c("pml_align.fasta", "pml_tree.nwk", "pml.rds", "pml.txt"))