From 338247149e395593ba5b1ef670f1b0f9e00ce79b Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Wed, 17 Jul 2024 15:53:25 +0200
Subject: [PATCH 01/17] small improvements to plot2D (plot.networx)

---
 R/plot_networx.R | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/R/plot_networx.R b/R/plot_networx.R
index 70de35f9..a48d80a1 100644
--- a/R/plot_networx.R
+++ b/R/plot_networx.R
@@ -325,7 +325,8 @@ plot.networx <- function(x, type = "equal angle", use.edge.length = TRUE,
     plot2D(coord, x, show.tip.label = show.tip.label,
            show.edge.label = show.edge.label, edge.label = edge.label,
            show.node.label = show.node.label, node.label = node.label,
-           show.nodes = show.nodes, tip.color = tip.color, edge.color = edge.color,
+           #show.nodes = show.nodes,
+           tip.color = tip.color, edge.color = edge.color,
            edge.width = edge.width, edge.lty = edge.lty, font = font, cex = cex,
            cex.node.label = cex.node.label, cex.edge.label = cex.edge.label,
            col.node.label = col.node.label, col.edge.label = col.edge.label,
@@ -408,24 +409,27 @@ plot2D <- function(coords, net, show.tip.label = TRUE, show.edge.label = FALSE,
   xx <- coords[, 1]
   yy <- coords[, 2]
   nTips <- length(label)
+  offset <- max(nchar(label)) * 0.018 * cex * diff(range(xx))
   if(is.null(xlim)){
     xlim <- range(xx)
-    offset <- max(nchar(label)) * 0.018 * cex * diff(xlim)
+#    offset <- max(nchar(label)) * 0.018 * cex * diff(xlim)
     if (show.tip.label) xlim <- c(xlim[1] - offset, xlim[2] + offset)
   }
   if(is.null(ylim)){
     ylim <- range(yy)
     if (show.tip.label){
       if(direction=="axial"){
-        offset <- max(nchar(label)) * 0.018 * cex * diff(ylim)
+#        offset <- max(nchar(label)) * 0.018 * cex * diff(ylim)
         ylim <- c(ylim[1] - offset, ylim[2] + offset)
       } else ylim <- c(ylim[1] - 0.03 * cex * diff(ylim),
                      ylim[2] + 0.03 * cex * diff(ylim))
     }
   }
   if (!add) {
-    plot.new()
-    plot.window(xlim, ylim, asp = 1)
+#    plot.new()
+#    plot.window(xlim, ylim, asp = 1, ...)
+    plot.default(0, type = "n", xlim = xlim, ylim = ylim, xlab = "",
+                 ylab = "", axes = FALSE, asp = 1, ...)
   }
   cladogram.plot(edge, xx, yy, edge.color, edge.width, edge.lty)
   if (show.tip.label) {

From c3ce9117508ca5de203795d09635f34ad9b96ebc Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Fri, 26 Jul 2024 11:16:08 +0200
Subject: [PATCH 02/17] add alternative to acctran. Small improvements to
 vignette.

---
 R/joint_ASR.R           | 67 +++++++++++++++++++++++++++++++++++++++++
 vignettes/Ancestral.Rmd | 11 ++++---
 2 files changed, 73 insertions(+), 5 deletions(-)

diff --git a/R/joint_ASR.R b/R/joint_ASR.R
index e279b46d..df332b1d 100644
--- a/R/joint_ASR.R
+++ b/R/joint_ASR.R
@@ -70,3 +70,70 @@ joint_pml <- function(x){
   res
 }
 
+
+# Joint reconstruction for parsimony
+joint_sankoff <- function(tree, data, cost=NULL){
+  stopifnot(inherits(data, "phyDat"))
+  stopifnot(inherits(tree, "phylo"))
+  tree <- reorder(tree, "postorder")
+  edge <- tree$edge
+  ntip <- Ntip(tree)
+  if (is.null(cost)) {
+    levels <- attr(data, "levels")
+    l <- length(levels)
+    cost <- matrix(1, l, l)
+    cost <- cost - diag(l)
+  }
+  P <- cost
+
+  desc <- Descendants(tree, type="children")
+  nr <- attr(data, "nr")
+  nc <- attr(data, "nc")
+  levels <- attr(data, "levels")
+  allLevels <- attr(data, "allLevels")
+
+  l <- nrow(tree$edge)
+  C <- array(NA, c(nr, nc, max(tree$edge)))
+  tmp <- sankoff(tree, data, cost, "data")
+  pscore <- tmp$pscore
+  L <- tmp$dat
+  nn <- unique(edge[,1])
+  pa <- 0
+  root <- edge[l, 1]
+  lnr <- seq_len(nr)
+
+  for(i in seq_len(length(nn))){
+    pa <- nn[i]
+    if(pa==root) break()
+    for(j in 1:nc){
+      pp <- L[[pa]] + rep(P[j,], each=nr)
+      pos <- max.col(-pp)
+      C[,j,pa] <- pos
+    }
+  }
+  pp <- L[[pa]]
+  pos <- max.col(-pp)
+  C[,,pa] <- pos
+  tree <- reorder(tree)
+  if(is.null(tree$node.label)) tree <- makeNodeLabel(tree)
+  res <- vector("list", length(tree$node.label))
+  names(res) <- tree$node.label
+  res[[1]] <- pos
+  att <- attributes(data)
+  att$names <- tree$node.label
+  labels <- c(tree$tip.label, tree$node.label)
+  edge <- tree$edge
+  nrw <- seq_len(nr)
+  for(i in seq_along(edge[,1])){
+    ch_i <- edge[i,2]
+    pa_i <- edge[i,1]
+    if(ch_i > ntip){
+      pos <-res[[labels[pa_i]]]
+      res[[labels[ch_i]]] <- C[cbind(nrw,pos,ch_i)]
+    }
+  }
+  ind <- match(levels, allLevels)
+  for(i in length(res))  res[[i]] <- ind[res[[i]]]
+  attributes(res) <- att
+  res
+}
diff --git a/vignettes/Ancestral.Rmd b/vignettes/Ancestral.Rmd
index caafacb8..6240d35b 100644
--- a/vignettes/Ancestral.Rmd
+++ b/vignettes/Ancestral.Rmd
@@ -38,16 +38,17 @@ parsimony(tree, primates)
 ```
 
 For parsimony analysis of the edge length represent the observed number of changes. Reconstructing ancestral states therefore defines also the edge lengths of a tree. However there can exist several equally parsimonious reconstructions or states can be ambiguous and therefore edge length can differ.
+In _phangorn_ all the ancestral reconstructions for parsimony used to be based on the fitch algorithm (below version 3.0) and needed bifurcating trees. However trees can get pruned afterwards using the function `multi2di` from _ape_. Recently we replaced the acctran routine with a method based on the sankoff algorithm adopting the algorithm for joint reconstruction [@Pupko2000] and breaking ties at random. This has the additional benefit that to infer phylogenies with multifurcations. 
+
 "MPR" reconstructs the ancestral states for each (internal) node as if the tree would be rooted in that node. However the nodes are not independent of each other. If one chooses one state for a specific node, this can restrict the choice of neighboring nodes (figures 2 and 3).
-The function acctran (accelerated transformation) assigns edge length and internal nodes to the tree [@Swofford1987].
+There is also an option "POSTORDER" which is only a one pass algotrithm, which we use for teaching purposes. 
+The function acctran (accelerated transformation) assigns edge length and internal nodes to the tree [@Swofford1987]. 
 ```{r ancestral_reconstruction}
 anc.acctran <- ancestral.pars(tree, primates, "ACCTRAN")
 anc.mpr <- ancestral.pars(tree, primates, "MPR")
 ```
 
-All the ancestral reconstructions for parsimony are based on the fitch algorithm and so far only bifurcating trees are allowed. However trees can get pruned afterwards using the function `multi2di` from _ape_.
-
-The `plotSeqLogo` function is a wrapper around the from the _ggseqlogo_ function the the _ggseqlogo_ package [@ggseqlogo]  and provides a simple way to show proportions of a nucleotides of ancestral states (see figure 1). 
+The `plotSeqLogo` function is a wrapper around the from the _ggseqlogo_ function in the _ggseqlogo_ package [@ggseqlogo]  and provides a simple way to show proportions of a nucleotides of ancestral states (see figure 1). 
 
 ```{r seqLogo, fig.cap="Fig 1. Ancestral reconstruction for a node.", eval=TRUE}
 #library(seqLogo)
@@ -68,7 +69,7 @@ title("ACCTRAN")
 
 _phangorn_ also offers the possibility to estimate ancestral states using ML.
 The advantages of ML over parsimony is that the reconstruction accounts for different edge lengths.
-So far only a marginal construction is implemented (see [@Yang2006][@Koshi1996]) and no joint reconstruction [@Pupko2000].
+Currently a marginal construction is implemented (see [@Yang2006][@Koshi1996]), but the joint reconstruction [@Pupko2000] only for models without rate variation (e.g. gamma models) or invariant sites.
 ```{r fit_ML}
 fit <- pml(tree, primates)
 fit <- optim.pml(fit, model="F81", control = pml.control(trace=0))

From 5d73e732fe29135f3a8c0e28b8ded20863c0bd78 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Mon, 26 Aug 2024 18:05:11 +0200
Subject: [PATCH 03/17] add support for rhub

---
 .github/workflows/rhub.yaml | 95 +++++++++++++++++++++++++++++++++++++
 1 file changed, 95 insertions(+)
 create mode 100644 .github/workflows/rhub.yaml

diff --git a/.github/workflows/rhub.yaml b/.github/workflows/rhub.yaml
new file mode 100644
index 00000000..74ec7b05
--- /dev/null
+++ b/.github/workflows/rhub.yaml
@@ -0,0 +1,95 @@
+# R-hub's generic GitHub Actions workflow file. It's canonical location is at
+# https://github.com/r-hub/actions/blob/v1/workflows/rhub.yaml
+# You can update this file to a newer version using the rhub2 package:
+#
+# rhub::rhub_setup()
+#
+# It is unlikely that you need to modify this file manually.
+
+name: R-hub
+run-name: "${{ github.event.inputs.id }}: ${{ github.event.inputs.name || format('Manually run by {0}', github.triggering_actor) }}"
+
+on:
+  workflow_dispatch:
+    inputs:
+      config:
+        description: 'A comma separated list of R-hub platforms to use.'
+        type: string
+        default: 'linux,windows,macos'
+      name:
+        description: 'Run name. You can leave this empty now.'
+        type: string
+      id:
+        description: 'Unique ID. You can leave this empty now.'
+        type: string
+
+jobs:
+
+  setup:
+    runs-on: ubuntu-latest
+    outputs:
+      containers: ${{ steps.rhub-setup.outputs.containers }}
+      platforms: ${{ steps.rhub-setup.outputs.platforms }}
+
+    steps:
+    # NO NEED TO CHECKOUT HERE
+    - uses: r-hub/actions/setup@v1
+      with:
+        config: ${{ github.event.inputs.config }}
+      id: rhub-setup
+
+  linux-containers:
+    needs: setup
+    if: ${{ needs.setup.outputs.containers != '[]' }}
+    runs-on: ubuntu-latest
+    name: ${{ matrix.config.label }}
+    strategy:
+      fail-fast: false
+      matrix:
+        config: ${{ fromJson(needs.setup.outputs.containers) }}
+    container:
+      image: ${{ matrix.config.container }}
+
+    steps:
+      - uses: r-hub/actions/checkout@v1
+      - uses: r-hub/actions/platform-info@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+      - uses: r-hub/actions/setup-deps@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+      - uses: r-hub/actions/run-check@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+
+  other-platforms:
+    needs: setup
+    if: ${{ needs.setup.outputs.platforms != '[]' }}
+    runs-on: ${{ matrix.config.os }}
+    name: ${{ matrix.config.label }}
+    strategy:
+      fail-fast: false
+      matrix:
+        config: ${{ fromJson(needs.setup.outputs.platforms) }}
+
+    steps:
+      - uses: r-hub/actions/checkout@v1
+      - uses: r-hub/actions/setup-r@v1
+        with:
+          job-config: ${{ matrix.config.job-config }}
+          token: ${{ secrets.RHUB_TOKEN }}
+      - uses: r-hub/actions/platform-info@v1
+        with:
+          token: ${{ secrets.RHUB_TOKEN }}
+          job-config: ${{ matrix.config.job-config }}
+      - uses: r-hub/actions/setup-deps@v1
+        with:
+          job-config: ${{ matrix.config.job-config }}
+          token: ${{ secrets.RHUB_TOKEN }}
+      - uses: r-hub/actions/run-check@v1
+        with:
+          job-config: ${{ matrix.config.job-config }}
+          token: ${{ secrets.RHUB_TOKEN }}

From 7fa487020860a02641792d0cde6a3eadb2e0a7b9 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Mon, 26 Aug 2024 20:08:40 +0200
Subject: [PATCH 04/17] fix error

---
 inst/tinytest/test_ancestral.R | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/inst/tinytest/test_ancestral.R b/inst/tinytest/test_ancestral.R
index bd01eb41..14fa79d3 100644
--- a/inst/tinytest/test_ancestral.R
+++ b/inst/tinytest/test_ancestral.R
@@ -29,7 +29,8 @@ test_acctran <- ancestral.pars(tree, dna, "ACCTRAN")
 
 
 data(Laurasiatherian)
-fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none")
+fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none",
+              control=pml.control(trace=0))
 anc_ml <- ancestral.pml(fit)
 write.ancestral(anc_ml)
 align <- read.phyDat("ancestral_align.fasta")
@@ -37,7 +38,7 @@ tree <- read.tree("ancestral_tree.nwk")
 df <- read.table("ancestral.state", header=TRUE)
 anc_ml_disc <- ancestral(tree, align, df)
 expect_equal(anc_ml[[1]], anc_ml_disc[[1]])
-expect_equal(anc_ml[[2]], anc_ml_disc[[2]])
+expect_equal(anc_ml[[2]][tree$tip.label], anc_ml_disc[[2]])
 expect_equal(anc_ml[[3]], anc_ml_disc[[3]])
 expect_equal(anc_ml[[4]], anc_ml_disc[[4]])
 unlink(c("ancestral_align.fasta", "ancestral_tree.nwk", "ancestral.state"))

From 902800363a7570ab128567442d39d4a6af2aaa99 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 27 Aug 2024 10:46:58 +0200
Subject: [PATCH 05/17] fix warning

---
 man/add.tips.Rd | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/man/add.tips.Rd b/man/add.tips.Rd
index 7689364f..bdd80ab3 100644
--- a/man/add.tips.Rd
+++ b/man/add.tips.Rd
@@ -31,7 +31,7 @@ tree1 <- add.tips(tree, c("A", "B", "C"), c(1,2,15))
 plot(tree1)
 }
 \seealso{
-\code{\link{bind.tree}}
+\code{\link[ape]{bind.tree}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}

From 6ab3adc1a74035c79e1d3f5698746e7aecef3ddf Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 27 Aug 2024 10:59:37 +0200
Subject: [PATCH 06/17] avoid warnings

---
 R/addConfidences.R     | 2 +-
 man/addConfidences.Rd  | 2 +-
 man/add_ci.Rd          | 2 +-
 man/add_edge_length.Rd | 5 +++--
 man/as.splits.Rd       | 2 +-
 5 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/R/addConfidences.R b/R/addConfidences.R
index 37b33788..6f6d1c08 100644
--- a/R/addConfidences.R
+++ b/R/addConfidences.R
@@ -15,7 +15,7 @@
 #' @return The object \code{x} with added bootstrap / MCMC support values.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{as.splits}}, \code{\link{as.networx}},
-#' \code{\link{RF.dist}}, \code{\link{plot.phylo}}
+#' \code{\link{RF.dist}}, \code{\link[ape]{plot.phylo}}
 #' @references Schliep, K., Potts, A. J., Morrison, D. A. and Grimm, G. W.
 #' (2017), Intertwining phylogenetic trees and networks.
 #' \emph{Methods Ecol Evol}.\bold{8}, 1212--1220. doi:10.1111/2041-210X.12760
diff --git a/man/addConfidences.Rd b/man/addConfidences.Rd
index 06ba9911..29c31797 100644
--- a/man/addConfidences.Rd
+++ b/man/addConfidences.Rd
@@ -65,7 +65,7 @@ Schliep, K., Potts, A. J., Morrison, D. A. and Grimm, G. W.
 }
 \seealso{
 \code{\link{as.splits}}, \code{\link{as.networx}},
-\code{\link{RF.dist}}, \code{\link{plot.phylo}}
+\code{\link{RF.dist}}, \code{\link[ape]{plot.phylo}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/add_ci.Rd b/man/add_ci.Rd
index 9a3c8535..83cf11a5 100644
--- a/man/add_ci.Rd
+++ b/man/add_ci.Rd
@@ -52,7 +52,7 @@ plot(tree, direction="downwards")
 add_boxplot(tree, trees, boxwex=.7)
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{plotBS}},
+\code{\link[ape]{plot.phylo}}, \code{\link{plotBS}},
 \code{\link{add_edge_length}}, \code{\link{maxCladeCred}}
 }
 \author{
diff --git a/man/add_edge_length.Rd b/man/add_edge_length.Rd
index 424cba2a..4b25d978 100644
--- a/man/add_edge_length.Rd
+++ b/man/add_edge_length.Rd
@@ -38,8 +38,9 @@ plot(tree_compat)
 add_boxplot(tree_compat, bs, boxwex=.7)
 }
 \seealso{
-\code{\link{node.depth.edgelength}}, \code{\link{consensus}},
-\code{\link{maxCladeCred}}, \code{\link{add_boxplot}}
+\code{\link[ape::node.depth]{node.depth.edgelength}},
+\code{\link[ape]{consensus}}, \code{\link{maxCladeCred}},
+\code{\link{add_boxplot}}
 }
 \author{
 Klaus Schliep
diff --git a/man/as.splits.Rd b/man/as.splits.Rd
index 4f82b6b5..bccb480e 100644
--- a/man/as.splits.Rd
+++ b/man/as.splits.Rd
@@ -105,7 +105,7 @@ plot(as.networx(spl))
 
 }
 \seealso{
-\code{\link{prop.part}}, \code{\link{lento}},
+\code{\link[ape]{prop.part}}, \code{\link{lento}},
 \code{\link{as.networx}}, \code{\link{distanceHadamard}},
 \code{\link{read.nexus.splits}}
 }

From 3b440d0c3963d88c7470d5f7bd4b39eb787d08cd Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 27 Aug 2024 11:04:57 +0200
Subject: [PATCH 07/17] fix notes

---
 R/ancestral.R         | 10 ++++++++--
 R/draw_CI.R           |  7 ++++---
 R/phyDat_conversion.R | 27 ++++++++++++++++++++++++++-
 R/splits.R            |  2 +-
 R/treeManipulation.R  |  2 +-
 man/ancestral.pml.Rd  |  2 +-
 man/as.phyDat.Rd      |  5 ++++-
 7 files changed, 45 insertions(+), 10 deletions(-)

diff --git a/R/ancestral.R b/R/ancestral.R
index ce2a45d9..203d75d2 100644
--- a/R/ancestral.R
+++ b/R/ancestral.R
@@ -42,7 +42,7 @@
 ## returned.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{pml}}, \code{\link{parsimony}}, \code{\link[ape]{ace}},
-#' \code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
+#' \code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link[ape]{latag2n}},
 #' \code{\link{gap_as_state}}, \code{\link[ape]{root}},
 #' \code{\link[ape]{makeNodeLabel}}
 #' @references Felsenstein, J. (2004). \emph{Inferring Phylogenies}. Sinauer
@@ -107,6 +107,9 @@ ancestral.pml <- function(object, type = "marginal", ...) {
   tree$node.label <- node_label
   tmp <- length(data)
 
+  joint <- TRUE
+  if(length(w) > 1 || object$inv > 0) joint <- FALSE
+
   eig <- object$eig
 
   bf <- object$bf
@@ -165,12 +168,15 @@ ancestral.pml <- function(object, type = "marginal", ...) {
       result2[[k]][ind] <- data[[1]][ind]
     }
   }
+  result3 <- NULL
+  if(joint) result3 <- joint_pml(object)
   attrib$names <- node_label
   attributes(result2) <- attrib
   attributes(result) <- attrib
   result <- list2df_ancestral(result, result2)
   result2 <- compress.phyDat(result2)
-  erg <- list(tree=tree, data=data, prob=result, state=result2)
+  erg <- list(tree=tree, data=data, prob=result, state=result2,
+              joint_state = result3)
   class(erg) <- "ancestral"
   erg
 }
diff --git a/R/draw_CI.R b/R/draw_CI.R
index 06b89042..3d6d34ce 100644
--- a/R/draw_CI.R
+++ b/R/draw_CI.R
@@ -58,8 +58,9 @@ edge_length_matrix <- function(tree, trees, rooted=TRUE){
 ##'               add_edge_length(bs)
 ##' plot(tree_compat)
 ##' add_boxplot(tree_compat, bs, boxwex=.7)
-##' @seealso \code{\link{node.depth.edgelength}}, \code{\link{consensus}},
-##' \code{\link{maxCladeCred}}, \code{\link{add_boxplot}}
+##' @seealso \code{\link[ape::node.depth]{node.depth.edgelength}},
+##' \code{\link[ape]{consensus}}, \code{\link{maxCladeCred}},
+##' \code{\link{add_boxplot}}
 ##' @keywords aplot
 ##' @export
 add_edge_length <- function(tree, trees, fun=\(x)median(na.omit(x)),
@@ -109,7 +110,7 @@ add_edge_length <- function(tree, trees, fun=\(x)median(na.omit(x)),
 ##' add_ci(tree, trees)
 ##' plot(tree, direction="downwards")
 ##' add_boxplot(tree, trees, boxwex=.7)
-##' @seealso \code{\link{plot.phylo}}, \code{\link{plotBS}},
+##' @seealso \code{\link[ape]{plot.phylo}}, \code{\link{plotBS}},
 ##' \code{\link{add_edge_length}}, \code{\link{maxCladeCred}}
 ##' @keywords aplot
 ##' @rdname add_ci
diff --git a/R/phyDat_conversion.R b/R/phyDat_conversion.R
index de664302..1484b84f 100644
--- a/R/phyDat_conversion.R
+++ b/R/phyDat_conversion.R
@@ -30,7 +30,7 @@
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso [DNAbin()], [as.DNAbin()],
 #' \code{\link{baseFreq}}, \code{\link{glance.phyDat}},
-#' \code{\link{read.dna}}, \code{\link{read.nexus.data}}
+#' \code{\link[ape]{read.dna}}, \code{\link[ape]{read.nexus.data}}
 #' and the chapter 1 in the \code{vignette("phangorn-specials",
 #' package="phangorn")} and the example of \code{\link{pmlMix}} for the use of
 #' \code{allSitePattern}
@@ -158,6 +158,31 @@ as.phyDat.DNAStringSet <- function(x, ...){
 }
 
 
+# @rdname phyDat
+#' @export
+as.StringSet <- function (x, ...){
+  if (inherits(x, "StringSet")) return(x)
+  UseMethod("as.StringSet")
+}
+
+
+
+#' @rdname as.phyDat
+#' @export
+as.StringSet.phyDat <- function(x, ...){
+  if (requireNamespace("Biostrings")){
+    z <- as.character(x)
+    type <- attr(x, "type")
+    seq <- switch(type,
+                  DNA = tolower(apply(z, 1, paste, collapse="")),
+                  AA = toupper(apply(z, 1, paste, collapse="")))
+    if(type=="DNA") return(Biostrings::DNAStringSet(seq))
+    if(type=="AA") return(Biostrings::AAStringSet(seq))
+  }
+  return(NULL)
+}
+
+
 # @rdname phyDat
 #' @export
 as.MultipleAlignment <- function (x, ...){
diff --git a/R/splits.R b/R/splits.R
index 88d7abc9..c1fefc63 100644
--- a/R/splits.R
+++ b/R/splits.R
@@ -26,7 +26,7 @@
 #' two splits are incompatible.
 #' @note The internal representation is likely to change.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{prop.part}}, \code{\link{lento}},
+#' @seealso \code{\link[ape]{prop.part}}, \code{\link{lento}},
 #' \code{\link{as.networx}}, \code{\link{distanceHadamard}},
 #' \code{\link{read.nexus.splits}}
 #' @keywords cluster
diff --git a/R/treeManipulation.R b/R/treeManipulation.R
index fd0e4837..a98541f2 100644
--- a/R/treeManipulation.R
+++ b/R/treeManipulation.R
@@ -430,7 +430,7 @@ addOneTree <- function(tree, subtree, i, node) {
 #' @param edge.length optional numeric vector with edge length
 #' @return an object of class phylo
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{bind.tree}}
+#' @seealso \code{\link[ape]{bind.tree}}
 #' @keywords cluster
 #' @examples
 #' tree <- rcoal(10)
diff --git a/man/ancestral.pml.Rd b/man/ancestral.pml.Rd
index b84f3a60..b41a1189 100644
--- a/man/ancestral.pml.Rd
+++ b/man/ancestral.pml.Rd
@@ -92,7 +92,7 @@ Press, Oxford.
 }
 \seealso{
 \code{\link{pml}}, \code{\link{parsimony}}, \code{\link[ape]{ace}},
-\code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
+\code{\link{plotAnc}}, \code{\link{ltg2amb}}, \code{\link[ape]{latag2n}},
 \code{\link{gap_as_state}}, \code{\link[ape]{root}},
 \code{\link[ape]{makeNodeLabel}}
 }
diff --git a/man/as.phyDat.Rd b/man/as.phyDat.Rd
index afba3199..81c902d7 100644
--- a/man/as.phyDat.Rd
+++ b/man/as.phyDat.Rd
@@ -18,6 +18,7 @@
 \alias{as.phyDat.MultipleAlignment}
 \alias{as.phyDat.AAStringSet}
 \alias{as.phyDat.DNAStringSet}
+\alias{as.StringSet.phyDat}
 \alias{as.character.phyDat}
 \alias{as.data.frame.phyDat}
 \alias{as.DNAbin.phyDat}
@@ -45,6 +46,8 @@ phyDat2alignment(x)
 
 \method{as.phyDat}{DNAStringSet}(x, ...)
 
+\method{as.StringSet}{phyDat}(x, ...)
+
 \method{as.MultipleAlignment}{phyDat}(x, ...)
 
 \method{as.character}{phyDat}(x, allLevels = TRUE, ...)
@@ -113,7 +116,7 @@ all.equal(Laurasiatherian, as.phyDat(LauraDNAbin))
 \seealso{
 [DNAbin()], [as.DNAbin()],
 \code{\link{baseFreq}}, \code{\link{glance.phyDat}},
-\code{\link{read.dna}}, \code{\link{read.nexus.data}}
+\code{\link[ape]{read.dna}}, \code{\link[ape]{read.nexus.data}}
 and the chapter 1 in the \code{vignette("phangorn-specials",
 package="phangorn")} and the example of \code{\link{pmlMix}} for the use of
 \code{allSitePattern}

From 469b46976fb4f2a54f76e4a844e680122c931015 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 27 Aug 2024 11:37:42 +0200
Subject: [PATCH 08/17] fix notes

---
 R/Coalescent.R               | 2 +-
 R/Densi.R                    | 4 ++--
 R/baseFreq.R                 | 2 +-
 R/bootstrap.R                | 4 ++--
 R/cladePar.R                 | 2 +-
 R/codon.R                    | 2 +-
 R/gap_as_state.R             | 2 +-
 R/image_phyDat.R             | 2 +-
 R/ltg2amb.R                  | 2 +-
 R/maxCladeCred.R             | 8 ++++----
 R/phyDat_conversion.R        | 2 +-
 R/phylo.R                    | 3 ++-
 R/plotBS.R                   | 8 ++++----
 R/plot_pml.R                 | 4 ++--
 R/pml_bb.R                   | 4 ++--
 R/read.nexus.partitions.R    | 2 +-
 R/read.nexus.splits.R        | 2 +-
 R/transferBootstrap.R        | 2 +-
 R/upgma.R                    | 4 ++--
 man/as.phyDat.Rd             | 2 +-
 man/baseFreq.Rd              | 2 +-
 man/bootstrap.pml.Rd         | 4 ++--
 man/cladePar.Rd              | 2 +-
 man/coalSpeciesTree.Rd       | 2 +-
 man/densiTree.Rd             | 4 ++--
 man/dna2codon.Rd             | 2 +-
 man/gap_as_state.Rd          | 2 +-
 man/image.phyDat.Rd          | 2 +-
 man/ltg2amb.Rd               | 2 +-
 man/maxCladeCred.Rd          | 8 ++++----
 man/plot.pml.Rd              | 4 ++--
 man/plotBS.Rd                | 8 ++++----
 man/pml.Rd                   | 3 ++-
 man/pml_bb.Rd                | 4 ++--
 man/read.nexus.partitions.Rd | 2 +-
 man/read.nexus.splits.Rd     | 2 +-
 man/transferBootstrap.Rd     | 2 +-
 man/upgma.Rd                 | 4 ++--
 38 files changed, 62 insertions(+), 60 deletions(-)

diff --git a/R/Coalescent.R b/R/Coalescent.R
index 21925ad1..d572c43c 100644
--- a/R/Coalescent.R
+++ b/R/Coalescent.R
@@ -64,7 +64,7 @@ comp2 <- function(x, y) {
 #' @param sTree A species tree which fixes the topology.
 #' @return The function returns an object of class \code{phylo}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com} Emmanuel Paradies
-#' @seealso \code{\link{speciesTree}}
+#' @seealso \code{\link[ape]{speciesTree}}
 #' @references Liu, L., Yu, L. and Pearl, D. K. (2010) Maximum tree: a
 #' consistent estimator of the species tree. \emph{Journal of Mathematical
 #' Biology}, \bold{60}, 95--106.
diff --git a/R/Densi.R b/R/Densi.R
index 8708fe78..ec48e18a 100644
--- a/R/Densi.R
+++ b/R/Densi.R
@@ -97,8 +97,8 @@ add_tiplabels <- function(xy, tip.label, direction, adj, font, srt = 0, cex = 1,
 #' @param \dots further arguments to be passed to plot.
 #' @returns \code{densiTree} returns silently x.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{plot.phylo}}, \code{\link{plot.networx}},
-#' \code{\link{jitter}}, \code{\link{rtt}}
+#' @seealso \code{\link[ape]{plot.phylo}}, \code{\link{plot.networx}},
+#' \code{\link{jitter}}, \code{\link[ape]{rtt}}
 #' @references densiTree is inspired from the great
 #' \href{https://www.cs.auckland.ac.nz/~remco/DensiTree/}{DensiTree} program of
 #' Remco Bouckaert.
diff --git a/R/baseFreq.R b/R/baseFreq.R
index 8e8629d7..fd32dac5 100644
--- a/R/baseFreq.R
+++ b/R/baseFreq.R
@@ -17,7 +17,7 @@
 #'
 #' @return  \code{baseFreq} returns a named vector and \code{glance} a one row
 #' \code{data.frame}.
-#' @seealso \code{\link{phyDat}, \link{base.freq}, \link{glance}}
+#' @seealso \code{\link{phyDat}, \link[ape]{base.freq}, \link{glance}}
 #' @author Klaus Schliep
 #' @examples
 #'
diff --git a/R/bootstrap.R b/R/bootstrap.R
index b0d6ac67..6c2105d1 100644
--- a/R/bootstrap.R
+++ b/R/bootstrap.R
@@ -29,8 +29,8 @@
 #' if not supplied the tree with labels supplied in the \code{node.label} slot.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{optim.pml}}, \code{\link{pml}},
-#' \code{\link{plot.phylo}}, \code{\link{maxCladeCred}}
-#' \code{\link{nodelabels}},\code{\link{consensusNet}} and
+#' \code{\link[ape]{plot.phylo}}, \code{\link{maxCladeCred}}
+#' \code{\link[ape]{nodelabels}},\code{\link{consensusNet}} and
 #' \code{\link{SOWH.test}} for parametric bootstrap
 #' @references Felsenstein J. (1985) Confidence limits on phylogenies. An
 #' approach using the bootstrap. \emph{Evolution} \bold{39}, 783--791
diff --git a/R/cladePar.R b/R/cladePar.R
index fc28b477..103652f8 100644
--- a/R/cladePar.R
+++ b/R/cladePar.R
@@ -14,7 +14,7 @@
 #' @param \dots Further arguments passed to or from other methods.
 #' @return A list containing the information about the edges and tips.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{plot.phylo}}
+#' @seealso \code{\link[ape]{plot.phylo}}
 #' @keywords plot
 #' @examples
 #'
diff --git a/R/codon.R b/R/codon.R
index 64769a3e..e2fa28e4 100644
--- a/R/codon.R
+++ b/R/codon.R
@@ -45,7 +45,7 @@
 #' @return The functions return an object of class \code{phyDat}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @references \url{https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes}
-#' @seealso \code{\link{trans}}, \code{\link{phyDat}} and the chapter 4 in the
+#' @seealso \code{\link[ape]{trans}}, \code{\link{phyDat}} and the chapter 4 in the
 #' \code{vignette("phangorn-specials", package="phangorn")}
 #' @keywords cluster
 #' @examples
diff --git a/R/gap_as_state.R b/R/gap_as_state.R
index c6917a11..5d0897f1 100644
--- a/R/gap_as_state.R
+++ b/R/gap_as_state.R
@@ -10,7 +10,7 @@
 #' @param ambiguous a character which codes for the ambiguous state
 #' @return The functions return an object of class \code{phyDat}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{phyDat}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
+#' @seealso \code{\link{phyDat}}, \code{\link[ape]{ltg2amb}}, \code{\link{latag2n}},
 #' \code{\link{ancestral.pml}}, \code{\link{gap_as_state}}
 #' @keywords cluster
 #' @examples
diff --git a/R/image_phyDat.R b/R/image_phyDat.R
index 8a686477..33ea8156 100644
--- a/R/image_phyDat.R
+++ b/R/image_phyDat.R
@@ -7,7 +7,7 @@
 #' @param x	 an object containing sequences, an object of class \code{phyDat}.
 #' @param ... further arguments passed to or from other methods.
 #' @returns Nothing. The function is called for plotting.
-#' @seealso \code{\link{image.DNAbin}}, \code{\link{image.AAbin}}
+#' @seealso \code{\link[ape]{image.DNAbin}}, \code{\link[ape]{image.AAbin}}
 #' @examples
 #' data("chloroplast")
 #' image(chloroplast[, 1:50], scheme="Clustal", show.aa = TRUE)
diff --git a/R/ltg2amb.R b/R/ltg2amb.R
index 09c5aa1a..7b984f72 100644
--- a/R/ltg2amb.R
+++ b/R/ltg2amb.R
@@ -7,7 +7,7 @@
 #' @param amb character of the ambiguous state t replace the gaps.
 #' @param gap gap parameter to replace.
 #' @returns returns an object of class \code{phyDat}.
-#' @seealso \code{\link{latag2n}}, \code{\link{ancestral.pml}},
+#' @seealso \code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
 #' \code{\link{gap_as_state}}
 #' @examples
 #' x <- phyDat(matrix(c("-", "A", "G", "-", "T", "C"), 2, 3))
diff --git a/R/maxCladeCred.R b/R/maxCladeCred.R
index 42a3a0c4..ab17d088 100644
--- a/R/maxCladeCred.R
+++ b/R/maxCladeCred.R
@@ -24,10 +24,10 @@
 #' @return a tree (an object of class \code{phylo}) with the highest clade
 #' credibility or a numeric vector of clade credibilities for each tree.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{consensus}}, \code{\link{consensusNet}},
-#' \code{\link{prop.part}}, \code{\link{bootstrap.pml}}, \code{\link{plotBS}},
-#' \code{\link{transferBootstrap}}, \code{\link{add_edge_length}},
-#' \code{\link{add_boxplot}}
+#' @seealso \code{\link[ape]{consensus}}, \code{\link{consensusNet}},
+#' \code{\link[ape]{prop.part}}, \code{\link{bootstrap.pml}},
+#' \code{\link{plotBS}}, \code{\link{transferBootstrap}},
+#' \code{\link{add_edge_length}}, \code{\link{add_boxplot}}
 #' @keywords cluster
 #' @importFrom fastmatch fmatch
 #' @examples
diff --git a/R/phyDat_conversion.R b/R/phyDat_conversion.R
index 1484b84f..9b9a3665 100644
--- a/R/phyDat_conversion.R
+++ b/R/phyDat_conversion.R
@@ -28,7 +28,7 @@
 #' @param ... further arguments passed to or from other methods.
 #' @return The functions return an object of class \code{phyDat}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso [DNAbin()], [as.DNAbin()],
+#' @seealso \code{\link[ape]{DNAbin}}, \code{\link[ape]{as.DNAbin}},
 #' \code{\link{baseFreq}}, \code{\link{glance.phyDat}},
 #' \code{\link[ape]{read.dna}}, \code{\link[ape]{read.nexus.data}}
 #' and the chapter 1 in the \code{vignette("phangorn-specials",
diff --git a/R/phylo.R b/R/phylo.R
index 9f3fa581..9c4e13cc 100644
--- a/R/phylo.R
+++ b/R/phylo.R
@@ -1247,7 +1247,8 @@ pml.fit <- function(tree, data, bf = rep(1 / length(levels), length(levels)),
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{pml_bb}}, \code{\link{bootstrap.pml}},
 #' \code{\link{modelTest}}, \code{\link{pmlPart}}, \code{\link{pmlMix}},
-#' \code{\link{plot.phylo}}, \code{\link{SH.test}}, \code{\link{ancestral.pml}}
+#' \code{\link[ape]{plot.phylo}}, \code{\link{SH.test}},
+#' \code{\link{ancestral.pml}}
 #' @references Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a
 #' maximum likelihood approach. \emph{Journal of Molecular Evolution},
 #' \bold{17}, 368--376.
diff --git a/R/plotBS.R b/R/plotBS.R
index 5588f297..1d2ce7d7 100644
--- a/R/plotBS.R
+++ b/R/plotBS.R
@@ -75,10 +75,10 @@ support <- function(tree, trees, method="FBP", tol=2e-8, scale=TRUE){
 #' \code{trees} is optional and if not supplied the labels supplied
 #' in the \code{node.label} slot will be used.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso  \code{\link{plot.phylo}}, \code{\link{add_ci}},
-#' \code{\link{nodelabels}},
-#' \code{\link{prop.clades}}, \code{\link{maxCladeCred}},
-#' \code{\link{transferBootstrap}}, \code{\link{consensus}},
+#' @seealso  \code{\link[ape]{plot.phylo}}, \code{\link{add_ci}},
+#' \code{\link[ape]{nodelabels}},
+#' \code{\link[ape]{prop.clades}}, \code{\link{maxCladeCred}},
+#' \code{\link{transferBootstrap}}, \code{\link[ape]{consensus}},
 #' \code{\link{consensusNet}}
 #' @references Felsenstein J. (1985) Confidence limits on phylogenies. An
 #' approach using the bootstrap. \emph{Evolution} \bold{39}, 783--791
diff --git a/R/plot_pml.R b/R/plot_pml.R
index 553eb864..66507283 100644
--- a/R/plot_pml.R
+++ b/R/plot_pml.R
@@ -14,8 +14,8 @@
 #' @return \code{plot.pml} returns invisibly a list with arguments dexcribing the plot.
 #' For further details see the \code{plot.phylo}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{plot.phylo}}, \code{\link{axisPhylo}},
-#' \code{\link{add.scale.bar}}
+#' @seealso \code{\link[ape]{plot.phylo}}, \code{\link[ape]{axisPhylo}},
+#' \code{\link[ape]{add.scale.bar}}
 #' @keywords plot
 #' @examples
 #' fdir <- system.file("extdata/trees", package = "phangorn")
diff --git a/R/pml_bb.R b/R/pml_bb.R
index 9e2d1a7a..aeb28734 100644
--- a/R/pml_bb.R
+++ b/R/pml_bb.R
@@ -36,8 +36,8 @@
 #' @param \dots Further arguments passed to or from other methods.
 #' @return \code{pml_bb} returns an object of class pml.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{optim.pml}}, \code{\link{modelTest}}, \code{\link{rtt}},
-#' \code{\link{pml.control}}
+#' @seealso \code{\link{optim.pml}}, \code{\link{modelTest}},
+#' \code{\link[ape]{rtt}}, \code{\link{pml.control}}
 #' @keywords cluster
 #' @examples
 #'
diff --git a/R/read.nexus.partitions.R b/R/read.nexus.partitions.R
index b23ada8d..f68e1442 100644
--- a/R/read.nexus.partitions.R
+++ b/R/read.nexus.partitions.R
@@ -48,7 +48,7 @@ read.nexus.charset <- function(file){
 #' @return a list where each element is a 'phyDat' object or an object of class
 #' 'multiphyDat'.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{read.nexus.data}}, \code{\link{read.phyDat}}
+#' @seealso \code{\link[ape]{read.nexus.data}}, \code{\link{read.phyDat}}
 #' @keywords cluster
 #' @examples
 #' tree <- rtree(10)
diff --git a/R/read.nexus.splits.R b/R/read.nexus.splits.R
index 9f0e9f57..528abb4b 100644
--- a/R/read.nexus.splits.R
+++ b/R/read.nexus.splits.R
@@ -30,7 +30,7 @@
 #' assumes that different co-variables are tab delimited and the bipartition
 #' are separated with white-space. Comments in square brackets are ignored.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{prop.part}}, \code{\link{lento}},
+#' @seealso \code{\link[ape]{prop.part}}, \code{\link{lento}},
 #' \code{\link{as.splits}}, \code{\link{as.networx}}
 #' @keywords cluster
 #' @examples
diff --git a/R/transferBootstrap.R b/R/transferBootstrap.R
index e01533bb..9675ddbf 100644
--- a/R/transferBootstrap.R
+++ b/R/transferBootstrap.R
@@ -13,7 +13,7 @@
 #' the node labels.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso  \code{\link{plotBS}}, \code{\link{maxCladeCred}},
-#' \code{\link{drawSupportOnEdges}}
+#' \code{\link[ape]{drawSupportOnEdges}}
 #' @references Lemoine, F., Entfellner, J. B. D., Wilkinson, E., Correia, D.,
 #' Felipe, M. D., De Oliveira, T., & Gascuel, O. (2018). Renewing Felsenstein’s
 #' phylogenetic bootstrap in the era of big data. \emph{Nature},
diff --git a/R/upgma.R b/R/upgma.R
index 7ac5275d..5607696c 100644
--- a/R/upgma.R
+++ b/R/upgma.R
@@ -20,8 +20,8 @@
 #' @return A phylogenetic tree of class \code{phylo}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
 #' @seealso \code{\link{hclust}}, \code{\link{dist.hamming}}, \code{\link{NJ}},
-#' \code{\link{as.phylo}}, \code{\link{fastme}}, \code{\link{nnls.tree}},
-#' \code{\link{rtt}}
+#' \code{\link[ape]{as.phylo}}, \code{\link[ape]{fastme}},
+#' \code{\link{nnls.tree}}, \code{\link[ape]{rtt}}
 #' @references Sneath, P. H., & Sokal, R. R. (1973). \emph{Numerical taxonomy.
 #' The principles and practice of numerical classification.}
 #'
diff --git a/man/as.phyDat.Rd b/man/as.phyDat.Rd
index 81c902d7..4b69ddf6 100644
--- a/man/as.phyDat.Rd
+++ b/man/as.phyDat.Rd
@@ -114,7 +114,7 @@ all.equal(Laurasiatherian, as.phyDat(LauraDNAbin))
 
 }
 \seealso{
-[DNAbin()], [as.DNAbin()],
+\code{\link[ape]{DNAbin}}, \code{\link[ape]{as.DNAbin}},
 \code{\link{baseFreq}}, \code{\link{glance.phyDat}},
 \code{\link[ape]{read.dna}}, \code{\link[ape]{read.nexus.data}}
 and the chapter 1 in the \code{vignette("phangorn-specials",
diff --git a/man/baseFreq.Rd b/man/baseFreq.Rd
index 5ba3339e..822796a1 100644
--- a/man/baseFreq.Rd
+++ b/man/baseFreq.Rd
@@ -51,7 +51,7 @@ glance(chloroplast)
 composition_test(Laurasiatherian)[1:10,]
 }
 \seealso{
-\code{\link{phyDat}, \link{base.freq}, \link{glance}}
+\code{\link{phyDat}, \link[ape]{base.freq}, \link{glance}}
 }
 \author{
 Klaus Schliep
diff --git a/man/bootstrap.pml.Rd b/man/bootstrap.pml.Rd
index 1616bc4e..71613afc 100644
--- a/man/bootstrap.pml.Rd
+++ b/man/bootstrap.pml.Rd
@@ -100,8 +100,8 @@ trees. \emph{Molecular Biology and Evolution} \bold{3}, 403--417
 }
 \seealso{
 \code{\link{optim.pml}}, \code{\link{pml}},
-\code{\link{plot.phylo}}, \code{\link{maxCladeCred}}
-\code{\link{nodelabels}},\code{\link{consensusNet}} and
+\code{\link[ape]{plot.phylo}}, \code{\link{maxCladeCred}}
+\code{\link[ape]{nodelabels}},\code{\link{consensusNet}} and
 \code{\link{SOWH.test}} for parametric bootstrap
 }
 \author{
diff --git a/man/cladePar.Rd b/man/cladePar.Rd
index 0fc2d612..fbfacca8 100644
--- a/man/cladePar.Rd
+++ b/man/cladePar.Rd
@@ -42,7 +42,7 @@ cladePar(tree, 18, "blue", "blue", x=x, plot=TRUE)
 
 }
 \seealso{
-\code{\link{plot.phylo}}
+\code{\link[ape]{plot.phylo}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/coalSpeciesTree.Rd b/man/coalSpeciesTree.Rd
index 878be9ed..0af0610b 100644
--- a/man/coalSpeciesTree.Rd
+++ b/man/coalSpeciesTree.Rd
@@ -40,7 +40,7 @@ consistent estimator of the species tree. \emph{Journal of Mathematical
 Biology}, \bold{60}, 95--106.
 }
 \seealso{
-\code{\link{speciesTree}}
+\code{\link[ape]{speciesTree}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com} Emmanuel Paradies
diff --git a/man/densiTree.Rd b/man/densiTree.Rd
index a6d88b1a..7b4e63ad 100644
--- a/man/densiTree.Rd
+++ b/man/densiTree.Rd
@@ -117,8 +117,8 @@ Remco R. Bouckaert (2010) DensiTree: making sense of sets of phylogenetic
 trees \emph{Bioinformatics}, \bold{26 (10)}, 1372-1373.
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{plot.networx}},
-\code{\link{jitter}}, \code{\link{rtt}}
+\code{\link[ape]{plot.phylo}}, \code{\link{plot.networx}},
+\code{\link{jitter}}, \code{\link[ape]{rtt}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/dna2codon.Rd b/man/dna2codon.Rd
index f9077f4c..88e6c958 100644
--- a/man/dna2codon.Rd
+++ b/man/dna2codon.Rd
@@ -75,7 +75,7 @@ dna2codon(Laurasiatherian)
 \url{https://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=cgencodes}
 }
 \seealso{
-\code{\link{trans}}, \code{\link{phyDat}} and the chapter 4 in the
+\code{\link[ape]{trans}}, \code{\link{phyDat}} and the chapter 4 in the
 \code{vignette("phangorn-specials", package="phangorn")}
 }
 \author{
diff --git a/man/gap_as_state.Rd b/man/gap_as_state.Rd
index 0cb61e83..6f0fc687 100644
--- a/man/gap_as_state.Rd
+++ b/man/gap_as_state.Rd
@@ -36,7 +36,7 @@ rownames(contr) <- attr(tmp, "allLevels")
 contr
 }
 \seealso{
-\code{\link{phyDat}}, \code{\link{ltg2amb}}, \code{\link{latag2n}},
+\code{\link{phyDat}}, \code{\link[ape]{ltg2amb}}, \code{\link{latag2n}},
 \code{\link{ancestral.pml}}, \code{\link{gap_as_state}}
 }
 \author{
diff --git a/man/image.phyDat.Rd b/man/image.phyDat.Rd
index 637c9e1b..3e7b1cd6 100644
--- a/man/image.phyDat.Rd
+++ b/man/image.phyDat.Rd
@@ -26,5 +26,5 @@ data("chloroplast")
 image(chloroplast[, 1:50], scheme="Clustal", show.aa = TRUE)
 }
 \seealso{
-\code{\link{image.DNAbin}}, \code{\link{image.AAbin}}
+\code{\link[ape]{image.DNAbin}}, \code{\link[ape]{image.AAbin}}
 }
diff --git a/man/ltg2amb.Rd b/man/ltg2amb.Rd
index 4e7b9655..63e9dc1e 100644
--- a/man/ltg2amb.Rd
+++ b/man/ltg2amb.Rd
@@ -28,7 +28,7 @@ image(x)
 image(y)
 }
 \seealso{
-\code{\link{latag2n}}, \code{\link{ancestral.pml}},
+\code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
 \code{\link{gap_as_state}}
 }
 \keyword{cluster}
diff --git a/man/maxCladeCred.Rd b/man/maxCladeCred.Rd
index 1de7ead8..82b07ce9 100644
--- a/man/maxCladeCred.Rd
+++ b/man/maxCladeCred.Rd
@@ -81,10 +81,10 @@ maxCladeCred(c(tree, bs[[1]]), tree=FALSE, part = pp)
 
 }
 \seealso{
-\code{\link{consensus}}, \code{\link{consensusNet}},
-\code{\link{prop.part}}, \code{\link{bootstrap.pml}}, \code{\link{plotBS}},
-\code{\link{transferBootstrap}}, \code{\link{add_edge_length}},
-\code{\link{add_boxplot}}
+\code{\link[ape]{consensus}}, \code{\link{consensusNet}},
+\code{\link[ape]{prop.part}}, \code{\link{bootstrap.pml}},
+\code{\link{plotBS}}, \code{\link{transferBootstrap}},
+\code{\link{add_edge_length}}, \code{\link{add_boxplot}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/plot.pml.Rd b/man/plot.pml.Rd
index 35a212e5..ed0a8f57 100644
--- a/man/plot.pml.Rd
+++ b/man/plot.pml.Rd
@@ -48,8 +48,8 @@ plot(fit_unrooted, cex=.5)
 
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{axisPhylo}},
-\code{\link{add.scale.bar}}
+\code{\link[ape]{plot.phylo}}, \code{\link[ape]{axisPhylo}},
+\code{\link[ape]{add.scale.bar}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/plotBS.Rd b/man/plotBS.Rd
index 04d50ea2..4f80e832 100644
--- a/man/plotBS.Rd
+++ b/man/plotBS.Rd
@@ -91,10 +91,10 @@ Penny D. and Hendy M.D. (1986) Estimating the reliability of evolutionary
 trees. \emph{Molecular Biology and Evolution} \bold{3}, 403--417
 }
 \seealso{
-\code{\link{plot.phylo}}, \code{\link{add_ci}},
-\code{\link{nodelabels}},
-\code{\link{prop.clades}}, \code{\link{maxCladeCred}},
-\code{\link{transferBootstrap}}, \code{\link{consensus}},
+\code{\link[ape]{plot.phylo}}, \code{\link{add_ci}},
+\code{\link[ape]{nodelabels}},
+\code{\link[ape]{prop.clades}}, \code{\link{maxCladeCred}},
+\code{\link{transferBootstrap}}, \code{\link[ape]{consensus}},
 \code{\link{consensusNet}}
 }
 \author{
diff --git a/man/pml.Rd b/man/pml.Rd
index 9d79876b..451b3319 100644
--- a/man/pml.Rd
+++ b/man/pml.Rd
@@ -232,7 +232,8 @@ discrete morphological character data. \emph{Systematic Biology} \bold{50},
 \seealso{
 \code{\link{pml_bb}}, \code{\link{bootstrap.pml}},
 \code{\link{modelTest}}, \code{\link{pmlPart}}, \code{\link{pmlMix}},
-\code{\link{plot.phylo}}, \code{\link{SH.test}}, \code{\link{ancestral.pml}}
+\code{\link[ape]{plot.phylo}}, \code{\link{SH.test}},
+\code{\link{ancestral.pml}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/pml_bb.Rd b/man/pml_bb.Rd
index 333386f2..fb79d66e 100644
--- a/man/pml_bb.Rd
+++ b/man/pml_bb.Rd
@@ -69,8 +69,8 @@ fit_HKY_R2 <- pml_bb(woodmouse, model="HKY+R(2)")
 }
 }
 \seealso{
-\code{\link{optim.pml}}, \code{\link{modelTest}}, \code{\link{rtt}},
-\code{\link{pml.control}}
+\code{\link{optim.pml}}, \code{\link{modelTest}},
+\code{\link[ape]{rtt}}, \code{\link{pml.control}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/read.nexus.partitions.Rd b/man/read.nexus.partitions.Rd
index 468746c5..2000a44c 100644
--- a/man/read.nexus.partitions.Rd
+++ b/man/read.nexus.partitions.Rd
@@ -42,7 +42,7 @@ tmp
 unlink(zz)
 }
 \seealso{
-\code{\link{read.nexus.data}}, \code{\link{read.phyDat}}
+\code{\link[ape]{read.nexus.data}}, \code{\link{read.phyDat}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/read.nexus.splits.Rd b/man/read.nexus.splits.Rd
index 5b9d9660..87c75b78 100644
--- a/man/read.nexus.splits.Rd
+++ b/man/read.nexus.splits.Rd
@@ -75,7 +75,7 @@ write.splits(spl, print.labels = FALSE)
 
 }
 \seealso{
-\code{\link{prop.part}}, \code{\link{lento}},
+\code{\link[ape]{prop.part}}, \code{\link{lento}},
 \code{\link{as.splits}}, \code{\link{as.networx}}
 }
 \author{
diff --git a/man/transferBootstrap.Rd b/man/transferBootstrap.Rd
index fb714273..71989c32 100644
--- a/man/transferBootstrap.Rd
+++ b/man/transferBootstrap.Rd
@@ -45,7 +45,7 @@ phylogenetic bootstrap in the era of big data. \emph{Nature},
 }
 \seealso{
 \code{\link{plotBS}}, \code{\link{maxCladeCred}},
-\code{\link{drawSupportOnEdges}}
+\code{\link[ape]{drawSupportOnEdges}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}
diff --git a/man/upgma.Rd b/man/upgma.Rd
index ddf6be05..94d94f7c 100644
--- a/man/upgma.Rd
+++ b/man/upgma.Rd
@@ -56,8 +56,8 @@ samples under the assumption of a molecular clock using serial-sample UPGMA.
 }
 \seealso{
 \code{\link{hclust}}, \code{\link{dist.hamming}}, \code{\link{NJ}},
-\code{\link{as.phylo}}, \code{\link{fastme}}, \code{\link{nnls.tree}},
-\code{\link{rtt}}
+\code{\link[ape]{as.phylo}}, \code{\link[ape]{fastme}},
+\code{\link{nnls.tree}}, \code{\link[ape]{rtt}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}

From a7971d9b403938c2033f2e50735610b889d00239 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 27 Aug 2024 11:56:59 +0200
Subject: [PATCH 09/17] bug fixes

---
 NEWS                   | 13 ++++++-------
 R/draw_CI.R            |  2 +-
 R/gap_as_state.R       |  5 +++--
 man/add_edge_length.Rd |  2 +-
 man/gap_as_state.Rd    |  5 +++--
 5 files changed, 14 insertions(+), 13 deletions(-)

diff --git a/NEWS b/NEWS
index 9416a4ee..d555445b 100644
--- a/NEWS
+++ b/NEWS
@@ -1,4 +1,4 @@
-    CHANGES in PHANGORN VERSION 3.0.0
+    CHANGES in PHANGORN VERSION 2.12.0
 
 NEW FEATURES
 
@@ -18,15 +18,15 @@ NEW FEATURES
 
       character vector for the node argument and not only integers.
 
-    o Ancestral state reconstructionwas in parts rewritten. Output is now an
+    o Ancestral state reconstruction was rewritten. Output is now an object of
 
-      object of class ancestral. E.g. ancestral states of constant sites are now
+      class ancestral. E.g. ancestral states of constant sites are now always
 
-      always this state.
+      this state.
 
     o Improvements to several plot functions to get nicer plots out of the box.
 
-    o Nicer defaults for plot.pml, mainly for rooted trees.
+      Nicer defaults for plot.pml, mainly for rooted trees.
 
     o plotAnc got an argument scheme allowing to use different color schemes,
 
@@ -38,8 +38,7 @@ NEW FEATURES
 
     o new function composition_test comparing indicating possible problems
 
-      with base / state composition.
-
+      with base / state composition as in iqtree.
 
 OTHER CHANGES
 
diff --git a/R/draw_CI.R b/R/draw_CI.R
index 3d6d34ce..b33f9437 100644
--- a/R/draw_CI.R
+++ b/R/draw_CI.R
@@ -58,7 +58,7 @@ edge_length_matrix <- function(tree, trees, rooted=TRUE){
 ##'               add_edge_length(bs)
 ##' plot(tree_compat)
 ##' add_boxplot(tree_compat, bs, boxwex=.7)
-##' @seealso \code{\link[ape::node.depth]{node.depth.edgelength}},
+##' @seealso \code{\link[ape]{node.depth}},
 ##' \code{\link[ape]{consensus}}, \code{\link{maxCladeCred}},
 ##' \code{\link{add_boxplot}}
 ##' @keywords aplot
diff --git a/R/gap_as_state.R b/R/gap_as_state.R
index 5d0897f1..b7d1555e 100644
--- a/R/gap_as_state.R
+++ b/R/gap_as_state.R
@@ -10,8 +10,9 @@
 #' @param ambiguous a character which codes for the ambiguous state
 #' @return The functions return an object of class \code{phyDat}.
 #' @author Klaus Schliep \email{klaus.schliep@@gmail.com}
-#' @seealso \code{\link{phyDat}}, \code{\link[ape]{ltg2amb}}, \code{\link{latag2n}},
-#' \code{\link{ancestral.pml}}, \code{\link{gap_as_state}}
+#' @seealso \code{\link{phyDat}}, \code{\link{ltg2amb}},
+#' \code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
+#' \code{\link{gap_as_state}}
 #' @keywords cluster
 #' @examples
 #' data(Laurasiatherian)
diff --git a/man/add_edge_length.Rd b/man/add_edge_length.Rd
index 4b25d978..cfc75e67 100644
--- a/man/add_edge_length.Rd
+++ b/man/add_edge_length.Rd
@@ -38,7 +38,7 @@ plot(tree_compat)
 add_boxplot(tree_compat, bs, boxwex=.7)
 }
 \seealso{
-\code{\link[ape::node.depth]{node.depth.edgelength}},
+\code{\link[ape]{node.depth}},
 \code{\link[ape]{consensus}}, \code{\link{maxCladeCred}},
 \code{\link{add_boxplot}}
 }
diff --git a/man/gap_as_state.Rd b/man/gap_as_state.Rd
index 6f0fc687..76148104 100644
--- a/man/gap_as_state.Rd
+++ b/man/gap_as_state.Rd
@@ -36,8 +36,9 @@ rownames(contr) <- attr(tmp, "allLevels")
 contr
 }
 \seealso{
-\code{\link{phyDat}}, \code{\link[ape]{ltg2amb}}, \code{\link{latag2n}},
-\code{\link{ancestral.pml}}, \code{\link{gap_as_state}}
+\code{\link{phyDat}}, \code{\link{ltg2amb}},
+\code{\link[ape]{latag2n}}, \code{\link{ancestral.pml}},
+\code{\link{gap_as_state}}
 }
 \author{
 Klaus Schliep \email{klaus.schliep@gmail.com}

From d22268024703edb1a5323d928d206391d0cdfbb8 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 27 Aug 2024 15:33:37 +0200
Subject: [PATCH 10/17] allow w, g as input for pml

---
 R/phylo.R | 45 +++++++++++++++++++++++++--------------------
 1 file changed, 25 insertions(+), 20 deletions(-)

diff --git a/R/phylo.R b/R/phylo.R
index 9c4e13cc..ae41ff78 100644
--- a/R/phylo.R
+++ b/R/phylo.R
@@ -1341,7 +1341,7 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
                 rate = 1, model = NULL, site.rate = "gamma", ASC = FALSE, ...) {
   call <- match.call()
   extras <- match.call(expand.dots = FALSE)$...
-  pmla <- c("wMix", "llMix", "dnds", "tstv")
+  pmla <- c("wMix", "llMix", "dnds", "tstv", "w", "g")
   existing <- match(pmla, names(extras))
   wMix <- ifelse(is.na(existing[1]), 0,
     eval(extras[[existing[1]]], parent.frame()))
@@ -1349,14 +1349,18 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
 #    eval(extras[[existing[2]]], parent.frame()))
   llMix <- 0
   if(!is.na(existing[2])) llMix <- eval(extras[[existing[2]]], parent.frame())
-
-  # allow
   dnds <- tstv <- 1
   dnds <- ifelse(is.na(existing[3]), 1,
     eval(extras[[existing[3]]], parent.frame()))
   tstv <- ifelse(is.na(existing[4]), 1,
     eval(extras[[existing[4]]], parent.frame()))
-
+  w <- g <- NULL
+  if(!is.na(existing[5]) & !is.na(existing[6])) {
+    w <- eval(extras[[existing[5]]], parent.frame())
+    g <- eval(extras[[existing[6]]], parent.frame())
+    stopifnot(length(g) == length(w))
+    k <- length(w)
+  }
   if (!inherits(tree, "phylo")) stop("tree must be of class phylo")
   if (is.null(tree$edge.length)) stop("tree must have edge weights")
   if (any(duplicated(tree$tip.label))) stop("tree must have unique labels!")
@@ -1420,22 +1424,23 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
     Q <- rep(1, length(levels) * (length(levels) - 1) / 2)
   m <- 1
   eig <- edQt(bf = bf, Q = Q)
-  if(site.rate=="free_rate"){
-    w <- rep(1/k, k)
-    g <- rep(1, k)
-  }
-  else{
-    rw <- rates_n_weights(shape, k, site.rate)
-    g <- rw[, 1]
-    w <- rw[, 2]
-  }
-  if (inv > 0){
-    w <- (1 - inv) * w
-    g <- g / (1 - inv)
-  }
+  if(is.null(g) | is.null(w)){
+    if(site.rate=="free_rate"){
+      w <- rep(1/k, k)
+      g <- rep(1, k)
+    }
+    else{
+      rw <- rates_n_weights(shape, k, site.rate)
+      g <- rw[, 1]
+      w <- rw[, 2]
+    }
+    if (inv > 0){
+      w <- (1 - inv) * w
+      g <- g / (1 - inv)
+    }
 #  if (wMix > 0) w <- (1 - wMix) * w
-  g <- g * rate
-
+    g <- g * rate
+  }
   inv0 <- inv
   kmax <- k
   if (any(g < .gEps)) {
@@ -1461,7 +1466,7 @@ pml <- function(tree, data, bf = NULL, Q = NULL, inv = 0, k = 1, shape = 1,
   tmp <- pml.fit(tree, data, bf, shape = shape, k = k, Q = Q,
     levels = attr(data, "levels"), inv = inv, rate = rate, g = g, w = w,
     eig = eig, INV = INV, ll.0 = ll.0, llMix = llMix, wMix = wMix, site = TRUE,
-    ASC=ASC)
+    ASC=ASC, site.rate = site.rate)
   df <- ifelse(is.ultrametric(tree), tree$Nnode, length(tree$edge.length))
 
   df <- switch(type,

From f6da31c7046ab50c2402142aa7f9d25f9af35fb9 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 27 Aug 2024 18:03:25 +0200
Subject: [PATCH 11/17] small updates

---
 src/dist.c    | 2 +-
 src/ml.c      | 2 +-
 src/sankoff.c | 6 +++---
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/dist.c b/src/dist.c
index e6d7c770..08a061ba 100644
--- a/src/dist.c
+++ b/src/dist.c
@@ -1,7 +1,7 @@
 /*
  * dist.c
  *
- * (c) 2008-2021 Klaus Schliep (klaus.schliep@gmail.com)
+ * (c) 2008-2024 Klaus Schliep (klaus.schliep@gmail.com)
  *
  *
  * This code may be distributed under the GNU GPL
diff --git a/src/ml.c b/src/ml.c
index e7678e04..7ed105d7 100644
--- a/src/ml.c
+++ b/src/ml.c
@@ -1,7 +1,7 @@
 /*
  * ml.c
  *
- * (c) 2008-2021  Klaus Schliep (klaus.schliep@gmail.com)
+ * (c) 2008-2024  Klaus Schliep (klaus.schliep@gmail.com)
  *
  *
  * This code may be distributed under the GNU GPL
diff --git a/src/sankoff.c b/src/sankoff.c
index 0c1edbe4..214e8e37 100644
--- a/src/sankoff.c
+++ b/src/sankoff.c
@@ -1,7 +1,7 @@
 /*
  * dist.c
  *
- * (c) 2008-2023  Klaus Schliep (klaus.schliep@gmail.com)
+ * (c) 2008-2024  Klaus Schliep (klaus.schliep@gmail.com)
  *
  *
  * This code may be distributed under the GNU GPL
@@ -31,7 +31,7 @@ SEXP C_rowMin(SEXP sdat, SEXP sn, SEXP sk){
     return(result);
 }
 
-
+/*
 void rowMin2(double *dat, int n,  int k, double *res){
     int i, h;
     double x;
@@ -41,7 +41,7 @@ void rowMin2(double *dat, int n,  int k, double *res){
         res[i] = x;
     }
 }
-
+*/
 
 double get_ps(double *dat, int n,  int k, double *weight){
   int i, h;

From f3d45751f08e697be922122bb53459613962abd1 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Wed, 4 Sep 2024 12:09:25 +0200
Subject: [PATCH 12/17] test and try

---
 .github/workflows/release-file-system-image.yml | 17 +++++++++++++++++
 1 file changed, 17 insertions(+)
 create mode 100644 .github/workflows/release-file-system-image.yml

diff --git a/.github/workflows/release-file-system-image.yml b/.github/workflows/release-file-system-image.yml
new file mode 100644
index 00000000..87bf570d
--- /dev/null
+++ b/.github/workflows/release-file-system-image.yml
@@ -0,0 +1,17 @@
+# Workflow derived from https://github.com/r-wasm/actions/tree/v1/examples
+# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
+on:
+  release:
+    # Must republish release to update assets
+    types: [ published ]
+
+name: Build and deploy wasm R package image
+
+jobs:
+  release-file-system-image:
+    uses: r-wasm/actions/.github/workflows/release-file-system-image.yml@v1
+    permissions:
+      # For publishing artifact files to the release
+      contents: write
+      # To download GitHub Packages within action
+      repository-projects: read

From 8f25070251f68f054e166e2643daca9b47ee88e8 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Wed, 4 Sep 2024 16:46:32 +0200
Subject: [PATCH 13/17] testing

---
 src/fitch64.cpp | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/fitch64.cpp b/src/fitch64.cpp
index 2b20bf2b..cf71a805 100644
--- a/src/fitch64.cpp
+++ b/src/fitch64.cpp
@@ -846,6 +846,7 @@ double pscore(Fitch* obj, const IntegerMatrix & orig){
 
 NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
   int i,j;
+  size_t ij;
   int states = obj->nStates;
   int nBits = obj->nBits;
   std::vector< std::vector<uint64_t> > vector = obj->X;
@@ -867,6 +868,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
   uint64_t ones = ~0ull;
   uint64_t tmp = 0ull;
   for(int k=0; k<nl; k+=2){
+    ij = (size_t) anc[k] - 1L;
     child1 = vector[desc[k] - 1L].data();
     child2 = vector[desc[k+1] - 1L].data();
     parent = vector[anc[k] - 1L].data();
@@ -884,7 +886,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       parent += states;
       tmp = ~orvand & ones;
       for(int l=0; l<BIT_SIZE; ++l){
-        if((tmp >> l) & 1ull) pars[anc[k] - 1L]+= obj->weight[i*BIT_SIZE + l];
+        if((tmp >> l) & 1ull) pars[ij]+= obj->weight[i*BIT_SIZE + l];
       }
     }
     for (int i = obj->wBits; i < nBits; ++i){
@@ -897,12 +899,13 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       child2 += states;
       parent += states;
       tmp = ~orvand & ones;
-      pars[anc[k] - 1L] += popcnt64(tmp);
+      pars[ij] += popcnt64(tmp);
     }
   }
   if(unrooted){
     child1 = vector[desc[nl] - 1].data();
     parent = vector[anc[nl] - 1].data();
+    ij = (size_t) anc[nl] - 1L;
     for (i = 0; i < obj->wBits; ++i){
       // OR the ANDs of all states
       uint64_t orvand = 0ull;
@@ -915,7 +918,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       parent += states;
       uint64_t tmp = ~orvand & ones;
       for(int l=0; l<BIT_SIZE; ++l){
-        if((tmp >> l) & 1ull) pars[anc[nl] - 1L]+= obj->weight[i*BIT_SIZE + l];
+        if((tmp >> l) & 1ull) pars[ij]+= obj->weight[i*BIT_SIZE + l];
       }
     }
     for (int i = obj->wBits; i < nBits; ++i){
@@ -929,7 +932,7 @@ NumericVector pscore_node(Fitch* obj, const IntegerMatrix & orig){
       child1 += states;
       parent += states;
       uint64_t tmp = ~orvand & ones;
-      pars [anc[nl] - 1L]+= popcnt64(tmp);
+      pars [ij]+= popcnt64(tmp);
     }
   }
   return(pars);

From 0ddb6304211a1b711d7f9baef1946ee862f33c55 Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Wed, 4 Sep 2024 23:02:04 +0200
Subject: [PATCH 14/17] fix wasm errors

---
 src/phangorn_utils.cpp | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/phangorn_utils.cpp b/src/phangorn_utils.cpp
index de7e80d0..d8f3dde0 100644
--- a/src/phangorn_utils.cpp
+++ b/src/phangorn_utils.cpp
@@ -88,14 +88,14 @@ int countCycle_cpp(IntegerMatrix M){
 // [[Rcpp::export]]
 IntegerVector countCycle2_cpp(IntegerMatrix M){
   int tmp;
-  int l = M.nrow();
-  int m = M.ncol();
+  size_t l = M.nrow();
+  size_t m = M.ncol();
   IntegerVector res(l);
-  for (int i=0; i<l; i++) {
+  for (size_t i=0; i<l; i++) {
     tmp = 0L;
     if(M[i] != M[i + (m -1) * l])tmp=1L;
-    for (int j=1; j<m; j++) {
-      if(M[i + (j-1L)* l] != M[i + j * l])tmp++;
+    for (size_t j=1; j<m; j++) {
+      if(M[i + (j-1L)* l] != M[i + j * l])tmp++;    ///
     }
     res[i]=tmp;
   }
@@ -395,12 +395,12 @@ Fast cophenetic distance
 */
 
 
-void copheneticHelpCpp(std::vector<int> left, std::vector<int> right, int h, NumericVector nh, int nTips, NumericVector dm){
+void copheneticHelpCpp(std::vector<int> left, std::vector<int> right, size_t h, NumericVector nh, int nTips, NumericVector dm){
     int ind;
     for(std::size_t i=0; i<left.size(); i++){
         for(std::size_t j=0; j<right.size(); j++){
             ind = give_index3(left[i], right[j], nTips);
-            dm[ind] = 2.0*nh[h] - nh[left[i]-1L] - nh[right[j]-1L];
+            dm[ind] = 2.0*nh[h] - nh[left[i]-1] - nh[right[j]-1];
         }
     }
 }
@@ -416,7 +416,7 @@ NumericVector cophenetic_cpp(IntegerMatrix edge, NumericVector edge_length,
     std::vector< std::vector<int> > bip = bipCPP(edge, nTips);
     NumericVector dm( nTips * (nTips-1) /2 );
     int l=0, left, right;
-    for(int h=nNode; h<(nNode + nTips); h++){
+    for(size_t h=nNode; h<(nNode + nTips); h++){
         // changed from NumericVector to IntegerVector
         IntegerVector tmp_ch = ch[h];
         l = tmp_ch.size();

From 087b4e442b8df313043e247c55cbe7b0a34900ea Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Fri, 6 Sep 2024 18:17:59 +0200
Subject: [PATCH 15/17] Update NEWS

---
 NEWS | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/NEWS b/NEWS
index d555445b..a4aa1058 100644
--- a/NEWS
+++ b/NEWS
@@ -18,11 +18,17 @@ NEW FEATURES
 
       character vector for the node argument and not only integers.
 
-    o Ancestral state reconstruction was rewritten. Output is now an object of
+    o Ancestral state reconstruction was rewritten. anc_pml and anc_pars are now
 
-      class ancestral. E.g. ancestral states of constant sites are now always
+      prefered over ancestral.pml and ańcestral.pars.
 
-      this state.
+      Output is now an object of class ancestral, which contains tree and
+
+      reconstructions.
+
+      Joint reconstruction for ML is now available.
+
+      Ancestral states of constant sites are now always this state.
 
     o Improvements to several plot functions to get nicer plots out of the box.
 

From 50be4d8aaaf9509ca8a5c2035467b68193495dec Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 10 Sep 2024 00:37:14 +0200
Subject: [PATCH 16/17] improve rbind.phyDat

---
 R/phyDat.R | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/R/phyDat.R b/R/phyDat.R
index 64431b65..935cca99 100644
--- a/R/phyDat.R
+++ b/R/phyDat.R
@@ -177,6 +177,12 @@ rbind.phyDat <- function(...){
   for(i in seq_along(x)){
     res[(mcs[i]+1):mcs[i+1], ] <- as.character(x[[i]])
   }
+  if(types[1]=="USER"){
+    contrast <- attr(x[[1]], "contrast")
+    dimnames(contrast) <- list(attr(x[[1]], "allLevels"),
+                               attr(x[[1]], "levels"))
+    return(phyDat(res, type="USER", contrast=contrast))
+  }
   phyDat(res, type=types[1])
 }
 

From b49b151f90a65737918bbc99cdf144ac049f219b Mon Sep 17 00:00:00 2001
From: Klaus Schliep <klaus.schliep@gmail.com>
Date: Tue, 10 Sep 2024 00:45:28 +0200
Subject: [PATCH 17/17] update ancestral code

---
 R/ancestral.R                        | 259 ++++++++++++++++++---------
 R/phyDat_conversion.R                |   1 +
 R/plotAnc.R                          |  28 +--
 inst/tinytest/test_ancestral.R       |  20 +--
 man/ancestral.pml.Rd                 |  23 ++-
 man/as.phyDat.Rd                     |   3 +-
 man/plot.ancestral.Rd                |   4 +-
 man/write.ancestral.Rd               |  10 +-
 tests/testthat/test_plot_ancestral.R |   2 +-
 vignettes/Ancestral.Rmd              |  26 +--
 10 files changed, 239 insertions(+), 137 deletions(-)

diff --git a/R/ancestral.R b/R/ancestral.R
index 203d75d2..96ec2b48 100644
--- a/R/ancestral.R
+++ b/R/ancestral.R
@@ -29,7 +29,7 @@
 #' @param data an object of class phyDat
 #' @param type method used to assign characters to internal nodes, see details.
 #' @param cost A cost matrix for the transitions between two states.
-## @param return return a \code{phyDat} object or matrix of probabilities.
+#' @param return return a \code{phyDat} object or matrix of probabilities.
 ##  @param x an object of class ancestral.
 #' @param \dots Further arguments passed to or from other methods.
 #' @return An object of class ancestral. This is a list containing the tree with
@@ -61,13 +61,8 @@
 #' # generate node labels to ensure plotting will work
 #' tree <- makeNodeLabel(tree)
 #' fit <- pml(tree, Laurasiatherian)
-#' anc.ml <- ancestral.pml(fit, type = "ml")
-#' anc.p <- ancestral.pars(tree, Laurasiatherian)
-## \dontrun{
-## require(seqLogo)
-## seqLogo( t(subset(anc.ml, 48, 1:20)[[1]]), ic.scale=FALSE)
-## seqLogo( t(subset(anc.p, 48, 1:20)[[1]]), ic.scale=FALSE)
-## }
+#' anc.ml <- anc_pml(fit)
+#' anc.p <- anc_pars(tree, Laurasiatherian)
 #' # plot ancestral sequences at the root
 #' plotSeqLogo( anc.ml, 48, 1, 20)
 #' plotSeqLogo( anc.p, 48, 1, 20)
@@ -78,9 +73,10 @@
 #'
 #' @rdname ancestral.pml
 #' @export
-ancestral.pml <- function(object, type = "marginal", ...) {
+ancestral.pml <- function(object, type = "marginal", return = "prob", ...) {
   call <- match.call()
   pt <- match.arg(type, c("marginal", "ml", "bayes")) # "joint",
+  rt <- match.arg(return, c("prob", "phyDat", "ancestral"))
   tree <- object$tree
   INV <- object$INV
   inv <- object$inv
@@ -168,20 +164,34 @@ ancestral.pml <- function(object, type = "marginal", ...) {
       result2[[k]][ind] <- data[[1]][ind]
     }
   }
-  result3 <- NULL
-  if(joint) result3 <- joint_pml(object)
+#  browser()
   attrib$names <- node_label
   attributes(result2) <- attrib
+  if(rt == "phyDat") return(rbind(data, result2))
+  if(rt == "prob") {
+    tmp <- c(unclass(new2old.phyDat(data)), result)
+    attrib$names <- c(tree$tip.label, node_label)
+    attributes(tmp) <- attrib
+    return(tmp)
+  }
   attributes(result) <- attrib
-  result <- list2df_ancestral(result, result2)
-  result2 <- compress.phyDat(result2)
-  erg <- list(tree=tree, data=data, prob=result, state=result2,
-              joint_state = result3)
+  if(joint) result2 <- joint_pml(object)
+#  result <- list2df_ancestral(result, result2)
+#  result2 <- compress.phyDat(result2)
+  erg <- list(tree=tree, data=data, prob=result, state=result2)
   class(erg) <- "ancestral"
   erg
 }
 
 
+#' @rdname ancestral.pml
+#' @export
+anc_pml <- function(object, type = "marginal", ...) {
+  ancestral.pml(object, type=type, return="ancestral")
+}
+
+
+
 #' Export and convenience functions for  ancestral reconstructions
 #'
 #' \code{write.ancestral} allows to export ancestral reconstructions. It writes
@@ -198,21 +208,22 @@ ancestral.pml <- function(object, type = "marginal", ...) {
 #' @examples
 #' data(Laurasiatherian)
 #' fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none")
-#' anc_ml <- ancestral.pml(fit)
+#' anc_ml <- anc_pml(fit)
 #' write.ancestral(anc_ml)
 #' # Can be also results from iqtree
-#' align <- read.phyDat("ancestral_align.fasta")
+#' align <- read.phyDat("ancestral_align.fasta", format="fasta")
 #' tree <- read.tree("ancestral_tree.nwk")
 #' df <- read.table("ancestral.state", header=TRUE)
-#' anc_ml_disc <- ancestral(tree, align, df)
+#' anc_ml_disc <- as.ancestral(tree, align, df)
 #' plotAnc(anc_ml_disc, 20)
 #' unlink(c("ancestral_align.fasta", "ancestral_tree.nwk", "ancestral.state"))
 #' @rdname write.ancestral
 #' @export
 write.ancestral <- function(x, file="ancestral"){
   stopifnot(inherits(x, "ancestral"))
-  write.phyDat(x$data, file=paste0(file, "_align.fasta"))
-  write.table(x$prob, file=paste0(file, ".state"), quote=FALSE, row.names=FALSE,
+  write.phyDat(x$data, file=paste0(file, "_align.fasta"), format="fasta")
+  tab <- as.data.frame(x)
+  write.table(tab, file=paste0(file, ".state"), quote=FALSE, row.names=FALSE,
               sep="\t")
   write.tree(x$tree, file=paste0(file, "_tree.nwk"))
   invisible(x)
@@ -226,13 +237,14 @@ write.ancestral <- function(x, file="ancestral"){
 #' @importFrom utils head
 #' @rdname write.ancestral
 #' @export
-ancestral <- function(tree, align, prob){
+as.ancestral <- function(tree, align, prob){
   stopifnot(inherits(tree, "phylo"))
   stopifnot(inherits(align, "phyDat"))
   stopifnot(inherits(prob, "data.frame"))
   if(is.null(tree$node.label))stop("tree needs node.label")
   state <- extract_states(prob, attr(align, "type"),
                           levels=attr(align, "levels"))
+  prob <- extract_prob(prob, align, tree$node.label)
   erg <- list(tree=tree, data=align[tree$tip.label], prob=prob,
               state=state[tree$node.label])
   class(erg) <- "ancestral"
@@ -248,6 +260,21 @@ extract_states <- function(x, type, levels=NULL){
   phyDat(y, type=type, levels=levels)
 }
 
+extract_prob <- function(x, align, node_label){
+  index <- attr(align, "index")
+  idx <- which(!duplicated(index))
+  att <- attributes(align)
+  res <- vector("list", length(node_label))
+  for(i in seq_along(node_label)){
+     tmp <- x[x$"Node"==node_label[i], -c(1:3)] |> as.matrix()
+     dimnames(tmp) <- NULL
+     res[[i]] <- as.matrix(tmp)[idx,]
+  }
+  att$names <- node_label
+  attributes(res) <- att
+  res
+}
+
 #' @rdname write.ancestral
 #' @export
 print.ancestral <- function(x, ...){
@@ -255,8 +282,8 @@ print.ancestral <- function(x, ...){
   print(x$tree)
   cat("\n")
   print(x$data)
-  cat("\n")
-  print(head(x$prob))
+#  cat("\n")
+#  print(head(x$prob))
 }
 
 #' @rdname ancestral.pml
@@ -317,6 +344,46 @@ list2df_ancestral <- function(x, y=NULL, ...) {
   res
 }
 
+#' @rdname ancestral.pml
+#' @export
+as.data.frame.ancestral <- function(x, ...) {
+  stopifnot(inherits(x, "ancestral"))
+  y <- rbind(x$data, x$state)
+  contrast <- attr(x$data, "contrast")
+  Y <- unlist(as.data.frame(y))
+  n_node <- Nnode(x$tree)
+  n_tip <- Ntip(x$tree)
+  l <- length(y)
+  nr <- attr(x$data, "nr")
+  nc <- attr(x$data, "nc")
+  index <- attr(x$data, "index")
+  nr <- length(index)
+  nam <- names(y)
+  X <- matrix(0, l*length(index), nc)
+  j <- 0
+  for(i in seq_len(n_tip)){
+    X[(j+1):(j+nr), ] <- contrast[x$data[[i]][index], ]
+    j <- j + nr
+  }
+  for(i in seq_len(n_node)){
+    X[(j+1):(j+nr), ] <- x$prob[[i]][index, ]
+    j <- j + nr
+  }
+#  if(!is.null(y) & inherits(y, "phyDat")){
+#    y <- y[names(x)]
+#    Y <- unlist(as.data.frame(y))
+  res <- data.frame(Node=rep(nam, each=nr), Site=rep(seq_len(nr), l), Y, X)
+  colnames(res) <- c("Node", "Site", "State", paste0("p_", attr(x$data, "levels")))
+#  } else{
+#    res <- data.frame(Node=rep(nam, each=nr), Site=rep(seq_len(nr), l), X)
+#    colnames(res) <- c("Node", "Site", paste0("p_", attr(x, "levels")))
+#  }
+  rownames(res) <- NULL
+  res
+}
+
+
+
 
 fitchCoding2ambiguous <- function(x, type = "DNA") {
   y <- c(1L, 2L, 4L, 8L, 8L, 3L, 5L, 9L, 6L, 10L, 12L, 7L, 11L, 13L,
@@ -325,46 +392,59 @@ fitchCoding2ambiguous <- function(x, type = "DNA") {
 }
 
 
+
 #' @rdname ancestral.pml
 #' @export
 ancestral.pars <- function(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"),
-                           cost = NULL, ...) {
-  #, return = "prob"
+                           cost = NULL, return = "prob") {
   call <- match.call()
   type <- match.arg(type)
-  tree$node.label <- makeAncNodeLabel(tree, ...)
-  contrast <- attr(data, "contrast")
-  data <- data[tree$tip.label,]
+  tree$nodel.label <- makeAncNodeLabel(tree)
   if (type == "ACCTRAN" || type=="POSTORDER") {
-    #, return = return
-    res <- ptree(tree, data, acctran=(type == "ACCTRAN"))
-#    attr(res, "call") <- call
+    res <- ptree(tree, data, return = return, acctran=(type == "ACCTRAN"))#, tips=TRUE)
+    attr(res, "call") <- call
   }
   if (type == "MPR") {
-    res <- mpr(tree, data, cost = cost)
-    #, return = return
-#    attr(res, "call") <- call
+    res <- mpr(tree, data, cost = cost, return = return) #, tips=TRUE)
+    attr(res, "call") <- call
   }
-  result <- res[[1]]
-  result2 <- res[[2]]
+  res
+}
+
+
+
+#' @rdname ancestral.pml
+#' @export
+anc_pars <- function(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"),
+                           cost = NULL, ...) {
+  #
+  call <- match.call()
+  type <- match.arg(type)
+  tree$node.label <- makeAncNodeLabel(tree, ...)
+  contrast <- attr(data, "contrast")
+  data <- data[tree$tip.label,]
+
+  prob <- ancestral.pars(tree, data, type, cost)
+  joint <- joint_sankoff(tree, data, cost)
   ind <- identical_sites(data)
   if(length(ind)>0){
     for(k in seq_len(Nnode(tree))){
-      result[[k]][ind,] <- contrast[data[[1]][ind],]
-      result2[[k]][ind] <- data[[1]][ind]
+      prob[[k]][ind,] <- contrast[data[[1]][ind],]
+      joint[[k]][ind] <- data[[1]][ind]
     }
   }
 #  attrib$names <- node_label
 #  attributes(result2) <- attrib
 #  attributes(result) <- attrib
-  result <- list2df_ancestral(result, result2)
-  result2 <- compress.phyDat(result2) # needed?
-  erg <- list(tree=tree, data=data, prob=result, state=result2)
+#  result <- list2df_ancestral(result, result2)
+#  result2 <- compress.phyDat(result2) # needed?
+  erg <- list(tree=tree, data=data, prob=prob, state=joint, call=call)
   class(erg) <- "ancestral"
   erg
 }
 
 
+
 #' @rdname ancestral.pml
 #' @export
 pace <- ancestral.pars
@@ -392,11 +472,12 @@ mpr.help <- function(tree, data, cost = NULL) {
   res
 }
 
-# , return = "prob"
-mpr <- function(tree, data, cost = NULL, ...) {
+
+mpr <- function(tree, data, cost = NULL, return="prob", tips=FALSE, ...) {
   data <- subset(data, tree$tip.label)
   att <- attributes(data)
-  att$names <- tree$node.label
+  if(tips) att$names <- c(tree$tip.label, tree$node.label)
+  else att$names <- tree$node.label
   type <- att$type
   nr <- att$nr
   nc <- att$nc
@@ -416,31 +497,31 @@ mpr <- function(tree, data, cost = NULL, ...) {
   }
 #  for (i in 1:ntips) res[[i]] <- contrast[data[[i]], , drop = FALSE]
   for (i in (ntips + 1):m) res[[i]][] <- as.numeric(res[[i]] < RM)
-  res <- res[(ntips + 1):m]
-#  if (return == "prob") {
-    #        for(i in 1:ntips) res[[i]] <- contrast[data[[i]],,drop=FALSE]
-    res_prob <- lapply(res, fun)
-    attributes(res_prob) <- att
-    class(res_prob) <- c("ancestral", "phyDat")
-#  }
+  if(tips)for(i in seq_len(ntips)) res[[i]] <- contrast[data[[i]],,drop=FALSE]
+  else  res <- res[(ntips + 1):m]
+  res_prob <- lapply(res, fun)
+  attributes(res_prob) <- att
+  if(return=="prob") return(res_prob)
+#    class(res_prob) <- c("ancestral", "phyDat")
   #    else res[1:ntips] <- data[1:ntips]
-  fun2 <- function(x) {
-    x <- p2dna(x)
-    fitchCoding2ambiguous(x)
-  }
+  #fun2 <- function(x) {
+  #  x <- p2dna(x)
+  #  fitchCoding2ambiguous(x)
+  #}
 #  if (return != "prob") {
-  if (type == "DNA") {
-    res_state <- lapply(res, fun2)
-    attributes(res_state) <- att
-  }
-  else {
-    attributes(res) <- att
-    res_state <- highest_state(res)
-    attributes(res_state) <- att
-  }
+#  if (type == "DNA") {
+#    res_state <- lapply(res, fun2)
+#    attributes(res_state) <- att
+#  }
+#  else {
+  attributes(res) <- att
+  res_state <- highest_state(res)
+  attributes(res_state) <- att
+  if(tips)res_state[seq_len(ntips)] <- data
+#  }
 #    res[1:ntips] <- data
 #  }
-  list(res_prob, res_state)
+  res_state # list(res_prob, res_state)
 }
 
 
@@ -489,12 +570,13 @@ acctran <- function(tree, data) {
 }
 
 #, return = "prob"
-ptree <- function(tree, data, acctran=TRUE, ...) {
+ptree <- function(tree, data, acctran=TRUE, return = "prob", tips=FALSE, ...) {
   tree <- reorder(tree, "postorder")
   data <- subset(data, tree$tip.label)
   edge <- tree$edge
   att <- attributes(data)
-  att$names <- tree$node.label
+  att$names <- c(tree$tip.label, tree$node.label)
+  #else att$names <- tree$node.label
   nr <- att$nr
   type <- att$type
   m <- max(edge)
@@ -506,33 +588,44 @@ ptree <- function(tree, data, acctran=TRUE, ...) {
   tmp <- tmp[tmp[,2]>Ntip(tree),]
   if(length(tmp)>0 && acctran==TRUE)f$acctran_traverse(tmp)
   res <- res_state <- vector("list", nNode)
-#  res <- vector("list", m)
-  att$names <- tree$node.label #makeAncNodeLabel(tree, ...)
+  res <- vector("list", m)
+  att$names <- c(tree$tip.label, tree$node.label) #makeAncNodeLabel(tree, ...)
 #  else {
     fun <- function(X) {
       rs <- rowSums(X)
       X / rs
     }
     contrast <- att$contrast
-#    for(i in seq_len(nTip)) res[[i]] <- contrast[data[[i]], , drop=FALSE]
+    for(i in seq_len(nTip)) res[[i]] <- contrast[data[[i]], , drop=FALSE]
 #    for(i in (nTip+1):m) res[[i]] <- f$getAnc(i)[1:nr, , drop=FALSE]
-    for(i in seq_len(nNode)) res[[i]] <- f$getAnc(i+nTip)[seq_len(nr), , drop=FALSE]
+    for(i in seq_len(nNode)) res[[i+nTip]] <- f$getAnc(i+nTip)[seq_len(nr), , drop=FALSE]
     res <- lapply(res, fun)
     attributes(res) <- att
-    class(res) <- c("ancestral", "phyDat")
+#    class(res) <- c("ancestral", "phyDat")
 #  }
-  if(type=="DNA"){
-    indx <- c(1, 2, 6, 3, 7, 9, 12, 4, 8, 10, 13, 11, 14, 15, 16)
-    #    res[1:nTip] <- data[1:nTip]
-    for(i in seq_len(nNode)) #  (nTip+1):m)
-      res_state[[i]] <- indx[f$getAncAmb(i+nTip)[1:nr]]
-    attributes(res_state) <- att
-    #    return(res)
-  } else{
-    res_state <- highest_state(res)
-    class(res_state) <- "phyDat"
+  if(!tips) res <- res[tree$node.label]
+
+  if(return=="prob"){
+#    if(tips)
+    return(res)
+#    else return(res[tree$node.label])
   }
-  list(res, res_state)
+
+#  if(type=="DNA"){
+#    indx <- c(1, 2, 6, 3, 7, 9, 12, 4, 8, 10, 13, 11, 14, 15, 16)
+#    res_state[seq_len(nTip)] <- data
+#    for(i in (nTip+1):m)
+#      res_state[[i]] <- indx[f$getAncAmb(i+nTip)[1:nr]]
+#    attributes(res_state) <- att
+#    #    return(res)
+ # } else{
+    res_state <- highest_state(res)
+    attributes(res_state) <-  attributes(res)
+#    class(res_state) <- "phyDat"
+#  }
+#  if(tips) return(res_state)
+#  res_state[tree$node.label]
+    res_state
 }
 
 
diff --git a/R/phyDat_conversion.R b/R/phyDat_conversion.R
index 9b9a3665..b46792c7 100644
--- a/R/phyDat_conversion.R
+++ b/R/phyDat_conversion.R
@@ -15,6 +15,7 @@
 #' @aliases
 #' as.phyDat.character as.phyDat.data.frame as.phyDat.matrix
 #' as.MultipleAlignment as.MultipleAlignment.phyDat
+#' as.StringSet as.StringSet.phyDat
 #' acgt2ry phyDat2MultipleAlignment
 #' @param data An object containing sequences.
 #' @param x An object containing sequences.
diff --git a/R/plotAnc.R b/R/plotAnc.R
index 9ea7fca4..c0edf625 100644
--- a/R/plotAnc.R
+++ b/R/plotAnc.R
@@ -65,7 +65,7 @@ getAncDF <- function(x){
 #' example(NJ)
 #' # generate node labels to ensure plotting will work
 #' tree <- makeNodeLabel(tree)
-#' anc.p <- ancestral.pars(tree, Laurasiatherian)
+#' anc.p <- anc_pars(tree, Laurasiatherian)
 #' # plot the third character
 #' plotAnc(anc.p, 3, pos="topright")
 #' plotSeqLogo(anc.p, node="Node10", 1, 25)
@@ -73,7 +73,7 @@ getAncDF <- function(x){
 #' data(chloroplast)
 #' tree <- pratchet(chloroplast,  maxit=10, trace=0)
 #' tree <- makeNodeLabel(tree)
-#' anc.ch <- ancestral.pars(tree, chloroplast)
+#' anc.ch <- anc_pars(tree, chloroplast)
 #' image(as.phyDat(anc.ch)[, 1:25])
 #' plotAnc(anc.ch, 21, scheme="Ape_AA", pos="topleft")
 #' plotAnc(anc.ch, 21, scheme="Clustal", pos="topleft")
@@ -88,7 +88,7 @@ plotAnc <- function(x, i = 1, type="phylogram", ..., col = NULL,
   type <- match.arg(type, c("phylogram", "cladogram", "fan", "unrooted",
                             "radial", "tidy"))
   phylo_clado <- type %in% c("phylogram", "cladogram")
-  df <- getAncDF(x)
+  df <- as.data.frame(x)
   data <- x$data
   tree <- x$tree
   subset <- df[,"Site"] == i
@@ -158,32 +158,36 @@ plotAnc <- function(x, i = 1, type="phylogram", ..., col = NULL,
 my_ggseqlogo <-function (data, facet = "wrap", scales = "free_x", ncol = NULL,
           nrow = NULL, start=NULL, end=NULL, ...)
 {
-  x <- geom_logo(data = data, ...)
-  x[[2]] <- scale_x_continuous(limits = c(start-0.5, end+.5) ,
+  x <- ggseqlogo::geom_logo(data = data, ...)
+  x[[2]] <- ggplot2::scale_x_continuous(limits = c(start-0.5, end+.5) ,
                                breaks=pretty(seq(start, end)))
-  p <- ggplot() + x + theme_logo()
+  p <- ggplot2::ggplot() + x + ggseqlogo::theme_logo()
   if (!"list" %in% class(data)) return(p)
   facet <- match.arg(facet, c("grid", "wrap"))
   if (facet == "grid") {
-    p <- p + facet_grid(~seq_group, scales = scales)
+    p <- p + ggplot2::facet_grid(~seq_group, scales = scales)
   }
   else if (facet == "wrap") {
-    p <- p + facet_wrap(~seq_group, scales = scales, nrow = nrow, ncol = ncol)
+    p <- p + ggplot2::facet_wrap(~seq_group, scales = scales, nrow = nrow,
+                                 ncol = ncol)
   }
   return(p)
 }
 
 
 #' @rdname plot.ancestral
-#' @importFrom ggplot2 scale_x_continuous ggplot facet_grid facet_wrap
-#' @importFrom ggseqlogo geom_logo theme_logo
+## @importFrom ggplot2 scale_x_continuous ggplot facet_grid facet_wrap
+## @importFrom ggseqlogo geom_logo theme_logo
 #' @returns \code{plotSeqLogo} returns a ggplot object.
 #' @export
-plotSeqLogo <- function(x, node=getRoot(x$tree), start=1, end=10, scheme="Ape_NT", ...){
+plotSeqLogo <- function(x, node=getRoot(x$tree), start=1, end=10,
+                        scheme="Ape_NT", ...){
   stopifnot(inherits(x, "ancestral"))
+  chk <- requireNamespace("ggseqlogo", quietly = TRUE)
+  if (!chk) stop("package ggseqlogo needs to be not installed!\n")
   type <- attr(x$data, "type")
   tree <- x$tree
-  df <- getAncDF(x)
+  df <- as.data.frame(x)
   nodes <- c(tree$tip.label, tree$node.label)
   if(is.numeric(node)) node <- nodes[node]
   subset <- df[,"Node"] == node
diff --git a/inst/tinytest/test_ancestral.R b/inst/tinytest/test_ancestral.R
index 14fa79d3..b0bc2fb1 100644
--- a/inst/tinytest/test_ancestral.R
+++ b/inst/tinytest/test_ancestral.R
@@ -15,9 +15,9 @@ fit <- pml(tree, dna)
 # dna tests differs from other data types as it may returns ambiguous data
 # test ancestral generics
 #    test.ml1 <- ancestral.pml(fit, type = "ml")
-test_ml <- ancestral.pml(fit, type = "ml")
-test_mpr <- ancestral.pars(tree, dna, "MPR")
-test_acctran <- ancestral.pars(tree, dna, "ACCTRAN")
+#test_ml <- ancestral.pml(fit, type = "ml")
+#test_mpr <- ancestral.pars(tree, dna, "MPR")
+#test_acctran <- ancestral.pars(tree, dna, "ACCTRAN")
 
 #expect_equal(as.character(test_ml), as.character(test_acctran))
 #expect_equal(as.character(test_ml), as.character(test_mpr))
@@ -31,16 +31,16 @@ test_acctran <- ancestral.pars(tree, dna, "ACCTRAN")
 data(Laurasiatherian)
 fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none",
               control=pml.control(trace=0))
-anc_ml <- ancestral.pml(fit)
+anc_ml <- anc_pml(fit)
 write.ancestral(anc_ml)
-align <- read.phyDat("ancestral_align.fasta")
+align <- read.phyDat("ancestral_align.fasta", format = "fasta")
 tree <- read.tree("ancestral_tree.nwk")
 df <- read.table("ancestral.state", header=TRUE)
-anc_ml_disc <- ancestral(tree, align, df)
-expect_equal(anc_ml[[1]], anc_ml_disc[[1]])
-expect_equal(anc_ml[[2]][tree$tip.label], anc_ml_disc[[2]])
-expect_equal(anc_ml[[3]], anc_ml_disc[[3]])
-expect_equal(anc_ml[[4]], anc_ml_disc[[4]])
+anc_ml_disc <- as.ancestral(tree, align, df)
+expect_equal(anc_ml$tree, anc_ml_disc$tree)
+expect_equal(anc_ml$data[tree$tip.label], anc_ml_disc$data)
+#expect_equal(anc_ml[[3]], anc_ml_disc[[3]])
+expect_equal(anc_ml$state, anc_ml_disc$state)
 unlink(c("ancestral_align.fasta", "ancestral_tree.nwk", "ancestral.state"))
 
 
diff --git a/man/ancestral.pml.Rd b/man/ancestral.pml.Rd
index b41a1189..033b2afc 100644
--- a/man/ancestral.pml.Rd
+++ b/man/ancestral.pml.Rd
@@ -2,25 +2,38 @@
 % Please edit documentation in R/ancestral.R
 \name{ancestral.pml}
 \alias{ancestral.pml}
+\alias{anc_pml}
 \alias{as.phyDat.ancestral}
+\alias{as.data.frame.ancestral}
 \alias{ancestral.pars}
+\alias{anc_pars}
 \alias{pace}
 \title{Ancestral character reconstruction.}
 \usage{
-ancestral.pml(object, type = "marginal", ...)
+ancestral.pml(object, type = "marginal", return = "prob", ...)
+
+anc_pml(object, type = "marginal", ...)
 
 \method{as.phyDat}{ancestral}(x, ...)
 
+\method{as.data.frame}{ancestral}(x, ...)
+
 ancestral.pars(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"),
-  cost = NULL, ...)
+  cost = NULL, return = "prob")
 
-pace(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"), cost = NULL, ...)
+anc_pars(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"), cost = NULL,
+  ...)
+
+pace(tree, data, type = c("MPR", "ACCTRAN", "POSTORDER"), cost = NULL,
+  return = "prob")
 }
 \arguments{
 \item{object}{an object of class pml}
 
 \item{type}{method used to assign characters to internal nodes, see details.}
 
+\item{return}{return a \code{phyDat} object or matrix of probabilities.}
+
 \item{\dots}{Further arguments passed to or from other methods.}
 
 \item{x}{an object of class ancestral}
@@ -69,8 +82,8 @@ example(NJ)
 # generate node labels to ensure plotting will work
 tree <- makeNodeLabel(tree)
 fit <- pml(tree, Laurasiatherian)
-anc.ml <- ancestral.pml(fit, type = "ml")
-anc.p <- ancestral.pars(tree, Laurasiatherian)
+anc.ml <- anc_pml(fit)
+anc.p <- anc_pars(tree, Laurasiatherian)
 # plot ancestral sequences at the root
 plotSeqLogo( anc.ml, 48, 1, 20)
 plotSeqLogo( anc.p, 48, 1, 20)
diff --git a/man/as.phyDat.Rd b/man/as.phyDat.Rd
index 4b69ddf6..6bc1412d 100644
--- a/man/as.phyDat.Rd
+++ b/man/as.phyDat.Rd
@@ -7,6 +7,8 @@
 \alias{as.phyDat.matrix}
 \alias{as.MultipleAlignment}
 \alias{as.MultipleAlignment.phyDat}
+\alias{as.StringSet}
+\alias{as.StringSet.phyDat}
 \alias{acgt2ry}
 \alias{phyDat2MultipleAlignment}
 \alias{as.phyDat}
@@ -18,7 +20,6 @@
 \alias{as.phyDat.MultipleAlignment}
 \alias{as.phyDat.AAStringSet}
 \alias{as.phyDat.DNAStringSet}
-\alias{as.StringSet.phyDat}
 \alias{as.character.phyDat}
 \alias{as.data.frame.phyDat}
 \alias{as.DNAbin.phyDat}
diff --git a/man/plot.ancestral.Rd b/man/plot.ancestral.Rd
index e6a5bcf2..9b699bf5 100644
--- a/man/plot.ancestral.Rd
+++ b/man/plot.ancestral.Rd
@@ -57,7 +57,7 @@ For further details see vignette("Ancestral").
 example(NJ)
 # generate node labels to ensure plotting will work
 tree <- makeNodeLabel(tree)
-anc.p <- ancestral.pars(tree, Laurasiatherian)
+anc.p <- anc_pars(tree, Laurasiatherian)
 # plot the third character
 plotAnc(anc.p, 3, pos="topright")
 plotSeqLogo(anc.p, node="Node10", 1, 25)
@@ -65,7 +65,7 @@ plotSeqLogo(anc.p, node="Node10", 1, 25)
 data(chloroplast)
 tree <- pratchet(chloroplast,  maxit=10, trace=0)
 tree <- makeNodeLabel(tree)
-anc.ch <- ancestral.pars(tree, chloroplast)
+anc.ch <- anc_pars(tree, chloroplast)
 image(as.phyDat(anc.ch)[, 1:25])
 plotAnc(anc.ch, 21, scheme="Ape_AA", pos="topleft")
 plotAnc(anc.ch, 21, scheme="Clustal", pos="topleft")
diff --git a/man/write.ancestral.Rd b/man/write.ancestral.Rd
index 923ee026..2a0daa7c 100644
--- a/man/write.ancestral.Rd
+++ b/man/write.ancestral.Rd
@@ -2,13 +2,13 @@
 % Please edit documentation in R/ancestral.R
 \name{write.ancestral}
 \alias{write.ancestral}
-\alias{ancestral}
+\alias{as.ancestral}
 \alias{print.ancestral}
 \title{Export and convenience functions for  ancestral reconstructions}
 \usage{
 write.ancestral(x, file = "ancestral")
 
-ancestral(tree, align, prob)
+as.ancestral(tree, align, prob)
 
 \method{print}{ancestral}(x, ...)
 }
@@ -41,13 +41,13 @@ plotting capabilities of R.
 \examples{
 data(Laurasiatherian)
 fit <- pml_bb(Laurasiatherian[,1:100], "JC", rearrangement = "none")
-anc_ml <- ancestral.pml(fit)
+anc_ml <- anc_pml(fit)
 write.ancestral(anc_ml)
 # Can be also results from iqtree
-align <- read.phyDat("ancestral_align.fasta")
+align <- read.phyDat("ancestral_align.fasta", format="fasta")
 tree <- read.tree("ancestral_tree.nwk")
 df <- read.table("ancestral.state", header=TRUE)
-anc_ml_disc <- ancestral(tree, align, df)
+anc_ml_disc <- as.ancestral(tree, align, df)
 plotAnc(anc_ml_disc, 20)
 unlink(c("ancestral_align.fasta", "ancestral_tree.nwk", "ancestral.state"))
 }
diff --git a/tests/testthat/test_plot_ancestral.R b/tests/testthat/test_plot_ancestral.R
index 3b79e5b2..d384ca89 100644
--- a/tests/testthat/test_plot_ancestral.R
+++ b/tests/testthat/test_plot_ancestral.R
@@ -7,7 +7,7 @@ dat <- matrix(c("a", "a",
               dimnames = list(c("t1", "t2", "t3", "t4"), NULL))
 dna <- phyDat(dat)
 fit <- pml(tree, dna)
-test_ml <- ancestral.pml(fit, type = "ml")
+test_ml <- anc_pml(fit)
 
 
 test_that("Pie_plots works", {
diff --git a/vignettes/Ancestral.Rmd b/vignettes/Ancestral.Rmd
index 6240d35b..ccaa0f72 100644
--- a/vignettes/Ancestral.Rmd
+++ b/vignettes/Ancestral.Rmd
@@ -44,26 +44,18 @@ In _phangorn_ all the ancestral reconstructions for parsimony used to be based o
 There is also an option "POSTORDER" which is only a one pass algotrithm, which we use for teaching purposes. 
 The function acctran (accelerated transformation) assigns edge length and internal nodes to the tree [@Swofford1987]. 
 ```{r ancestral_reconstruction}
-anc.acctran <- ancestral.pars(tree, primates, "ACCTRAN")
-anc.mpr <- ancestral.pars(tree, primates, "MPR")
+anc.pars <- anc_pars(tree, primates)
 ```
 
 The `plotSeqLogo` function is a wrapper around the from the _ggseqlogo_ function in the _ggseqlogo_ package [@ggseqlogo]  and provides a simple way to show proportions of a nucleotides of ancestral states (see figure 1). 
 
 ```{r seqLogo, fig.cap="Fig 1. Ancestral reconstruction for a node.", eval=TRUE}
-#library(seqLogo)
-#seqLogo( t(subset(anc.mpr, getRoot(tree), 1:20)[[1]]), ic.scale=FALSE)
-plotSeqLogo(anc.mpr, node=getRoot(tree), 1, 20)
+plotSeqLogo(anc.pars, node=getRoot(tree), 1, 20)
 ```
 ```{r MPR, fig.cap="Fig 2. Ancestral reconstruction using MPR."}
-plotAnc(anc.mpr, 17)
+plotAnc(anc.pars, 17)
 title("MPR")
 ```
-```{r ACCTRAN, fig.cap="Fig 3. Ancestral reconstruction using ACCTRAN."}
-plotAnc(anc.acctran, 17)
-title("ACCTRAN")
-```
-
 
 # Likelihood reconstructions
 
@@ -82,11 +74,9 @@ and the highest posterior probability ("bayes") criterion:
 \[
 P(x_r=A) =  \frac{\pi_A L(x_r=A)}{\sum_{k \in \{A,C,G,T\}}\pi_k L(x_r=k)},
 \]
-where $L(x_r)$ is the joint probability of states at the tips and the state at the root $x_r$ and $\pi_i$ are the estimated base frequencies of state $i$.
-Both methods agree if all states (base frequencies) have equal probabilities.
+where $L(x_r)$ is the joint probability of states at the tips and the state at the root $x_r$ and $\pi_i$ are the estimated base frequencies of state $i$. Both methods agree if all states (base frequencies) have equal probabilities.
 ```{r ML_reconstruction}
-anc.ml <- ancestral.pml(fit, "ml")
-anc.bayes <- ancestral.pml(fit, "bayes")
+anc.ml <- anc_pml(fit)
 ```
 The differences of the two approaches for a specific site (17) are represented in the following figures.
 ```{r plotML, fig.cap="Fig 4. Ancestral reconstruction the using the maximum likelihood."}
@@ -94,8 +84,8 @@ plotAnc(anc.ml, 17)
 title("ML")
 ```
 ```{r plotB, fig.cap="Fig 5. Ancestral reconstruction using (empirical) Bayes."}
-plotAnc(anc.bayes, 17)
-title("Bayes")
+#plotAnc(anc.bayes, 17)
+#title("Bayes")
 ```
 
 # Fitting for discrete comparative data
@@ -143,7 +133,7 @@ all.equal(fit_SYM$logLik, fit_ace$loglik+log(1/3))
 
 
 ```{r SYM_reconstruction}
-anc_SYM <- ancestral.pml(fit_SYM, "ml")
+anc_SYM <- anc_pml(fit_SYM)
 plotAnc(anc_SYM)
 ```