R/all.equal.phylo.R

   1 ## all.equal.phylo.R (2009-07-05)
   2 ##
   3 ##     Global Comparison of two Phylogenies
   4
   5 ## Copyright 2006 Benoît Durand
   6 ##    modified by EP for the new coding of "phylo" (2006-10-04)
   7
   8 ## This file is part of the R-package `ape'.
   9 ## See the file ../COPYING for licensing issues.
  10
  11 ## Recherche de la correspondance entre deux arbres
  12 ## Parcours en profondeur et en parallèle des deux arbres (current et target)
  13 ## current, target: les deux arbres à comparer
  14 ## use.edge.length: faut-il comparer les longueurs de branches ?
  15 ## use.tip.label: faut-il comparer les étiquettes de feuilles ou seulement la
  16 ##      topologie des deux arbres ?
  17 ## index.return: si TRUE, retourner la matrice de correspondance entre noeuds
  18 ##      et feuilles, une matrice à deux colonnes (current et target) avec pour
  19 ##      chaque ligne des paires d'identifiants de noeuds/feuilles, tels qu'ils
  20 ##      apparaissent dans l'attribut 'edge' des objets phylo
  21 ## tolerance, scale: paramètres de comparaison des longueurs de branches
  22 ##      (voir 'all.equal')
  23 all.equal.phylo <- function(target, current,
  24                         use.edge.length = TRUE,
  25                         use.tip.label = TRUE,
  26                         index.return = FALSE,
  27                         tolerance = .Machine$double.eps ^ 0.5,
  28                         scale = NULL, ...) {
  29
  30         same.node <- function(i, j) {
  31                 # Comparaison de un noeud et une feuille
  32                 if (xor(i > Ntip1, j > Ntip2)) return(NULL)
  33                 # Comparaison de deux feuilles
  34                 if (i <= Ntip1) {
  35                         if (!use.tip.label) return(c(i, j))
  36                         if (current$tip.label[i] == target$tip.label[j])
  37                                 return(c(i, j))
  38                         return(NULL)
  39                 }
  40                 # Comparaison de deux noeuds
  41                 i.children <- which(current$edge[, 1] == i)
  42                 j.children <- which(target$edge[, 1] == j)
  43                 if (length(i.children) != length(j.children)) return(NULL)
  44                 correspondance <- NULL
  45                 for (i.child in i.children) {
  46                         corresp <- NULL
  47                         for (j.child in j.children) {
  48                                 if (!use.edge.length ||
  49                                     isTRUE(all.equal(current$edge.length[i.child],
  50                                                      target$edge.length[j.child],
  51                                                      tolerance = tolerance,
  52                                                      scale = scale)))
  53                                     corresp <- same.node(current$edge[i.child, 2],
  54                                                          target$edge[j.child, 2])
  55                                 if (!is.null(corresp)) break
  56                         }
  57                         if (is.null(corresp)) return(NULL)
  58                         correspondance <- c(correspondance, i, j, corresp)
  59                         j.children <- j.children[j.children != j.child]
  60                 }
  61                 return(correspondance)
  62         }
  63
  64         Ntip1 <- length(target$tip.label)
  65         Ntip2 <- length(current$tip.label)
  66         root1 <- Ntip1 + 1
  67         root2 <- Ntip2 + 1
  68         if (root1 != root2) return(FALSE)
  69         ## Fix by EP so that unrooted trees are correctly compared:
  70         if (!is.rooted(target) && !is.rooted(current)) {
  71             outg <- target$tip.label[1]
  72             if (! outg %in% current$tip.label) return(FALSE)
  73             target <- root(target, outg)
  74             current <- root(current, outg)
  75         }
  76         ## End
  77         result <- same.node(root1, root2)
  78         if (!isTRUE(index.return)) return(!is.null(result))
  79         if (is.null(result)) return(result)
  80         result <- t(matrix(result, nrow = 2))
  81       colnames(result) = c('current', 'target')
  82       return(result)
  83 }