Cache Oblivious Algorithms for Computing the Triplet Distance Between Trees\footnote{Research supported by the Danish National Research Foundation, grant DNRF84, Center for Massive Data Algorithmics (MADALGO).
We study the problem of computing the triplet distance between two rooted unordered trees with $n$ labeled leafs. Introduced by Dobson 1975, the triplet distance is the number of leaf triples that induce different topologies in the two trees. The current theoretically best algorithm is an $\mathrm{O(n \log n)$ time algorithm by Brodal et al. (SODA 2013). Recently Jansson et al. proposed a new algorithm that, while slower in theory, requiring $\mathrm{O(n \log^3 n)$ time, in practice it outperforms the theoretically faster $\mathrm{O(n \log n)$ algorithm. Both algorithms do not scale to external memory. We present two cache oblivious algorithms that combine the best of both worlds. The first algorithm is for the case when the two input trees are binary trees and the second a generalized algorithm for two input trees of arbitrary degree. Analyzed in the RAM model, both algorithms require $\mathrm{O(n \log n)$ time, and in the cache oblivious model $\mathrm{O(\frac{n{B \log_{2 \frac{n{M)$ I/Os. Their relative simplicity and the fact that they scale to external memory makes them achieve the best practical performance. We note that these are the first algorithms that scale to external memory, both in theory and practice, for this problem.
lipics-v2016-sample-article.tex
\documentclass[a4paper,UKenglish]{lipics-v2016}
%This is a template for producing LIPIcs articles. 
%See lipics-manual.pdf for further information.
%for A4 paper format use option "a4paper", for US-letter use option "letterpaper"
%for british hyphenation rules use option "UKenglish", for american hyphenation rules use option "USenglish"
% for section-numbered lemmas etc., use "numberwithinsect"
 
\usepackage{microtype}%if unwanted, comment out or use option "draft"

%%%%%%%%%%%ADDITIONAL%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{booktabs}
\usepackage{enumitem}

%\graphicspath{{./graphics/}}%helpful if your graphic files are in another directory
\graphicspath{{./figures/}}

%%%%%%%%%%%%%COMMANDS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\ietal}{\textit{et al.}~}
\newcommand{\floor}[1]{\lfloor #1 \rfloor}
\bibliographystyle{plainurl}% the recommended bibstyle

% Author macros::begin %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\title{Cache Oblivious Algorithms for Computing the Triplet Distance Between Trees\footnote{Research supported by the Danish National Research Foundation, grant DNRF84, Center for Massive Data Algorithmics (MADALGO).}}
\titlerunning{Cache Oblivious Algorithms for Computing the Triplet Distance Between Trees} %optional, in case that the title is too long; the running title should fit into the top page column

%% Please provide for each author the \author and \affil macro, even when authors have the same affiliation, i.e. for each author there needs to be the  \author and \affil macros
\author[1]{Gerth Stølting Brodal}
\author[2]{Konstantinos Mampentzidis}
\affil[1]{Department of Computer Science, Aarhus University, Aarhus, Denmark\\
  \texttt{gerth@cs.au.dk}}
\affil[2]{Department of Computer Science, Aarhus University, Aarhus, Denmark\\
  \texttt{kmampent@cs.au.dk}}
\authorrunning{G.\,S. Brodal and K. Mampentzidis} %mandatory. First: Use abbreviated first/middle names. Second (only in severe cases): Use first author plus 'et. al.'

\Copyright{Gerth Stølting Brodal and Konstantinos Mampentzidis}%mandatory, please use full first names. LIPIcs license is "CC-BY";  http://creativecommons.org/licenses/by/3.0/

\subjclass{G.2.2 Trees, G.2.1 Combinatorial Algorithms}% mandatory: Please choose ACM 1998 classifications from http://www.acm.org/about/class/ccs98-html . E.g., cite as "F.1.1 Models of Computation". 
\keywords{Phylogenetic tree, tree comparison, triplet distance, cache oblivious algorithm}% mandatory: Please provide 1-5 keywords
% Author macros::end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Editor-only macros:: begin (do not touch as author)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\EventEditors{Kirk Pruhs and Christian Sohler}
\EventNoEds{2}
\EventLongTitle{25th Annual European Symposium on Algorithms (ESA 2017)}
\EventShortTitle{ESA 2017}
\EventAcronym{ESA}
\EventYear{2017}
\EventDate{September 4--6, 2017}
\EventLocation{Vienna, Austria}
\EventLogo{}
\SeriesVolume{87}
\ArticleNo{1} % Set article-no=1
% Editor-only macros::end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\maketitle

\begin{abstract}
We study the problem of computing the triplet distance between two rooted unordered trees with $n$ labeled leafs. Introduced by Dobson 1975, the triplet distance is the number of leaf triples that induce different topologies in the two trees. The current theoretically best algorithm is an $\mathrm{O}(n \log n)$ time algorithm by Brodal \textit{et al.} (SODA 2013). Recently Jansson \textit{et al.} proposed a new algorithm that, while slower in theory, requiring $\mathrm{O}(n \log^3 n)$ time, in practice it outperforms the theoretically faster $\mathrm{O}(n \log n)$ algorithm. Both algorithms do not scale to external memory.

We present two cache oblivious algorithms that combine the best of both worlds. The first algorithm is for the case when the two input trees are binary trees and the second a generalized algorithm for two input trees of arbitrary degree. Analyzed in the RAM model, both algorithms require $\mathrm{O}(n \log n)$ time, and in the cache oblivious model $\mathrm{O}(\frac{n}{B} \log_{2} \frac{n}{M})$ I/Os. Their relative simplicity and the fact that they scale to external memory makes them achieve the best practical performance. We note that these are the first algorithms that scale to external memory, both in theory and practice, for this problem. \end{abstract}

\section{Introduction}
\label{sec:introduction}

\subparagraph*{Background.}
Trees are data structures that are often used to represent relationships. For example in the field of Biology, a tree can be used to represent evolutionary relationships, with the leafs corresponding to species that exist today, and internal nodes to ancestor species that existed in the past. For a fixed set of $n$ species, different data (e.g. DNA, morphological) or construction methods (e.g. Q*~\cite{Qmethod}, neighbor joining~\cite{NJmethod}) can lead to trees that look structurally different. An interesting question that arises then is, given two trees $T_{1}$ and $T_{2}$ over $n$ species, how different are they? An answer to this question could potentially be used to determine whether the difference is statistically significant or not, which in turn could help with evolutionary inferences. 

Several distance measures have been proposed in the past to compare two trees. A class of them includes distance measures that are based on how often certain features are different in the two trees. Common distance measures of this kind are the Robinson-Foulds distance~\cite{RobinsonFoulds}, the triplet distance~\cite{Dobson} and the quartet distance~\cite{Quartet}. The Robinson-Foulds distance counts how many leaf bipartitions are different, where a bipartition in a given tree is generated by removing a single edge from the tree. The triplet distance is only defined for rooted trees, and counts how many leaf triples induce different topologies in the two trees. The counterpart of the triplet distance for unrooted trees, is the quartet distance, which counts how many leaf quadruples induce different topologies in the two trees. 

Algorithms exist that can efficiently compute these distance measures. The Robinson-Foulds distance can be optimally computed in $\mathrm{O}(n)$ time~\cite{DaysAlgorithm}. The triplet distance can be computed in $\mathrm{O}(n \log n)$ time~\cite{TripletQuartetSoda13}. The quartet distance can be computed in $\mathrm{O}(dn\log n)$ time~\cite{TripletQuartetSoda13}, where $d$ is the maximal degree of any node in the two input trees. 

Note that the above bounds are in the RAM model. Previous work did not consider any other models, for example external memory models like the I/O model~\cite{IOmodel} and the cache oblivious model~\cite{COmodel}. Typically when talking about algorithms for external memory models, one might (sometimes incorrectly) think of algorithms that have to deal with large amounts of data. So any practical improvement that comes from an algorithm that scales to external memory compared to an equivalent that does not, can only be noticed if the inputs are large. However, this is not necessarily the case for cache oblivious algorithms. A cache oblivious algorithm, if built and implemented correctly, can take advantage of the L1, L2 and L3 caches that exist in the vast majority of computers these days and give a significant performance boost even for small inputs.

A trivial modification of the algorithm in~\cite{DaysAlgorithm}, can give a cache oblivious algorithm for computing the Robinson-Foulds distance that achieves the sorting bound by requiring $\mathrm{O}(\frac{n}{B}\log_{\frac{M}{B}} \frac{n}{M})$ I/Os.  For the triplet and quartet distance measures, no such trivial modifications exist. 

In this paper we focus on the triplet distance computation and present the first non trivial algorithms for computing the triplet distance between two rooted trees, that for the first time for this problem, also scale to external memory.

\subparagraph*{Problem Definition.}

\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.22\textwidth}
        \centering
        \includegraphics{tripletTopologies/2.pdf}
        \caption{$xy|z$}
    \end{subfigure}%
    ~ 
    \begin{subfigure}{0.22\textwidth}
        \centering
        \includegraphics{tripletTopologies/3.pdf}
        \caption{$xz|y$}
    \end{subfigure}%
    ~ 
    \begin{subfigure}{0.22\textwidth}
        \centering
        \includegraphics{tripletTopologies/4.pdf}
        \caption{$yz|x$}
    \end{subfigure}%
    ~
     \begin{subfigure}{0.22\textwidth}
        \centering
        \includegraphics{tripletTopologies/1.pdf}
        \caption{$xyz$}
    \end{subfigure}%
    \caption{Triplet topologies.}
    \label{fig:tripletTopologies}
\end{figure}

For a given rooted unordered tree $T$ where each leaf has a unique label, a \emph{triplet} is defined by a set of three leaf labels $x$, $y$ and $z$ and their induced topology in $T$. The four possible topologies are illustrated in Figure~\ref{fig:tripletTopologies}. The notation $xy|z$ is used to describe a triplet where the lowest common ancestor of $x$ and $y$ is at a lower depth than the lowest common ancestor of $z$ with either $x$ or $y$. Note that the triplet $xy|z$ is the same as the triplet $yx|z$ because $T$ is considered to be unordered. Similarly, notation $xyz$ is used to describe a triplet for which every pair of leafs has the same lowest common ancestor. This triplet can only appear if we allow nodes with degree three or larger in $T$.

For two such trees $T_{1}$ and $T_{2}$ that are built on $n$ identical leaf labels, the \emph{triplet distance} $D(T_{1}, T_{2})$ is the number of triplets that are different in $T_{1}$ and $T_{2}$. Let $S(T_{1}, T_{2})$ be the number of \emph{shared} triplets in the two trees, i.e. leaf triples with identical topologies in the two trees.  We have the relationship that $D(T_{1}, T_{2}) + S(T_{1},T_{2}) = \binom{n}{3}$.

Previous and new results for computing the triplet distance are shown in the table below. Note that the papers \cite{Critchlow,Bansal2008,bmc13,TripletQuartetSoda13,Jansson2015} do not provide an analysis of the algorithms in the cache oblivious model, so here we provide an upper bound.

\begin{table}[ht]
  \centering
  \begin{tabular}{clcccc}
  	\toprule
  	Year & Reference & Time & IOs & Space & Non-Binary Trees\\
  	\midrule
    1996&Critchlow \ietal\cite{Critchlow} & $\mathrm{O}(n^2)$&$\mathrm{O}(n^2)$ &$\mathrm{O}(n^2)$&no\\
   2011& Bansal \ietal\cite{Bansal2008} & $\mathrm{O}(n^2)$ &$\mathrm{O}(n^2)$&$\mathrm{O}(n^2)$&yes\\
    2013&Brodal \ietal\cite{bmc13} & $\mathrm{O}(n \log^{2} n)$ &$\mathrm{O}(n \log^2 n$)&$\mathrm{O}(n)$&no\\
    2013&Brodal \ietal\cite{TripletQuartetSoda13} & $\mathrm{O}(n \log n$) &$\mathrm{O}(n \log n$)&$\mathrm{O}(n \log n)$&yes\\
    2015&Jansson \ietal\cite{Jansson2015} & $\mathrm{O}(n \log^{3} n)$ &$\mathrm{O}(n \log^{3} n)$ &$\mathrm{O}(n \log n)$&yes\\
    2017&new & $\mathrm{O}(n \log n)$ & $\mathrm{O}(\frac{n}{B} \log_{2} \frac{n}{M})$ &$\mathrm{O}(n)$&yes\\\hline
  \end{tabular}
  \label{tab:table1}
\end{table}

\subparagraph*{Related Work.}

The triplet distance was first suggested as a method of comparing the shapes of trees by Dobson in 1975 \cite{Dobson}. The first non-trivial algorithmic result dates back to 1996, when Critchlow \ietal\cite{Critchlow} proposed an $\mathrm{O}(n^2)$ algorithm that however works only for binary trees. Bansal \ietal\cite{Bansal2008} introduced an $\mathrm{O}(n^2)$ algorithm that works for general (binary and non-binary) trees. Both of these algorithms use $\mathrm{O}(n^2)$ space. Brodal \ietal\cite{bmc13} introduced a new $\mathrm{O}(n^2)$ algorithm using only $\mathrm{O}(n)$ space for the case of binary trees that they showed how to optimize to reduce the time to $\mathrm{O}(n\log^2 n)$. This algorithm was also implemented and shown to be the most efficient in practice. Soon after, Brodal \ietal\cite{TripletQuartetSoda13} managed to extend the $\mathrm{O}(n\log^2 n)$ algorithm to work for general trees, and at the same time brought the time down to $\mathrm{O}(n\log n)$ but now with the space increased to $\mathrm{O}(n \log n)$. Note that the space for binary trees was still $\mathrm{O}(n)$. The algorithms from \cite{bmc13} and \cite{TripletQuartetSoda13} were implemented and added to the library tqDist \cite{bioinformatics14}. Interestingly, it was shown in \cite{alenex14} that for binary trees the $\mathrm{O}(n\log^2 n)$ algorithm had a better practical performance than the $\mathrm{O}(n\log n)$ algorithm. Jansson \ietal\cite{Jansson2015,pmid27983874} showed that an even slower theoretically algorithm requiring worst case $\mathrm{O}(n\log^3 n)$ time and $\mathrm{O}(n \log n)$ space could give the best practical performance, both for binary and non-binary trees. A detailed survey over previous results until 2013 can be found in \cite{biology13}.

\subparagraph*{Contribution.}

The common main bottleneck with all previous approaches is that the data structures used rely intensively on $\Omega(n \log n)$ random memory accesses. This means that all algorithms are penalized by cache performance and thus do not scale to external memory. %(see Figures in Appendix \ref{appendix:figures}).%
We address this limitation by proposing new algorithms for computing the triplet distance on binary and non-binary trees, that match the previous best $\mathrm{O}(n\log n)$ time and $\mathrm{O}(n)$ space bounds in the RAM model,  but for the first time also scale to external memory. More specifically, in the cache oblivious model, the total number of I/Os required is $\mathrm{O}(\frac{n}{B} \log_{2} \frac{n}{M})$. The basic idea is to essentially replace the dependency of random access to data structures by scanning contracted versions of the input trees. A careful implementation of the algorithms is shown to achieve the best practical performance, thus essentially documenting that the theoretical results carry over to practice.

\subparagraph*{Outline of the Paper.}
In Section \ref{sec:prevApproaches} we provide an overview of previous approaches. In Section \ref{sec:binary} we describe the new algorithm for the case where $T_{1}$ and $T_{2}$ are binary trees. In Section \ref{sec:general} we extend it to work for general trees. In Section \ref{sec:implementation} we provide some details related to our implementation of the new algorithms. In Section \ref{sec:experiments} we include some experiments illustrating the efficiency of the new algorithms compared to the previously known algorithms for this problem. Finally, in Section \ref{sec:conclusion} we provide some concluding remarks.

\section{Previous Approaches}
\label{sec:prevApproaches}
A naive algorithm would enumerate over all $\binom{n}{3}$ sets of 3 labels and find for each set whether the induced topologies in $T_{1}$ and $T_{2}$ differ or not, giving an $\mathrm{O}(n^3)$ algorithm. This naive approach does not exploit the fact that the triplets are not completely independent. For example the triplets $xy|z$ and $yx|u$ share the leafs $x$ and $y$ and the fact that the lowest common ancestor of $x$ and $y$ is at a lower depth than the lowest common ancestor of $z$ with either $x$ or $y$ and the lowest common ancestor of $u$ with either $x$ or $y$. Dependencies like this can be exploited to count the number of shared triplets faster.

Critchlow \ietal\cite{Critchlow} exploit the depth of the leafs' ancestors to achieve the first improvement over the naive approach. Bansal \ietal\cite{Bansal2008} exploit the shared leafs between subtrees and essentially reduce the problem to computing the intersection size (number of shared leafs) of all pairs of subtrees, one from $T_{1}$ and one from $T_{2}$, which has an elegant dynamic programming solution.

\subparagraph*{The $\mathrm{O}(n^2)$ Algorithm for Binary Trees in \cite{bmc13}.}

The algorithm for binary trees in \cite{bmc13} is the basis for all subsequent improvements \cite{bmc13, TripletQuartetSoda13, Jansson2015}, including ours as well, so we will describe it in more detail here. The dependency that was exploited is the same as in \cite{Bansal2008} but the procedure for counting the shared triplets is completely different. More specifically, each triplet in $T_{1}$ and $T_{2}$, defined by the leafs $i$, $j$ and $k$, is implicitly \emph{anchored} in the lowest common ancestor of $i$, $j$ and $k$. For a node $u$ in $T_{1}$ and $v$ in $T_{2}$, let $s(u)$ and $s(v)$ be the set of triplets that are anchored in $u$ and $v$ respectively. For the number of shared triplets $S(T_{1},T_{2})$ we then have that $$S(T_{1},T_{2}) = \sum_{u \in T_{1}}\sum_{v \in T_{2}} \lvert {s(u) \cap s(v)} \rvert \;.$$ 

For the algorithm to be $\mathrm{O}(n^2)$ the value  $\lvert {s(u) \cap s(v)} \rvert$ must be computed in $\mathrm{O}(1)$ time. This is achieved by a leaf colouring procedure as follows: Fix a node $u$ in $T_{1}$ and color the leafs in the left subtree of $u$ \textit{red}, the leafs in the right subtree of $u$  \textit{blue}, let every other leaf have no color and then transfer this coloring to the leafs in $T_{2}$, i.e. identically labelled leafs get the same color. To compute $\lvert {s(u) \cap s(v)} \rvert$ we do as follows: let $l$ and $r$ be the left and right children of $v$, and let $w_{\texttt{red}}$ and $w_{\texttt{blue}}$ be the number of red and blue leafs in a subtree rooted at a node $w$ in $T_{2}$. We then have that
%
\begin{equation}
\label{eq:eq1}
\lvert {s(u) \cap s(v)} \rvert = \binom{l_{\texttt{red}}}{2}r_{\texttt{blue}} + \binom{l_{\texttt{blue}}}{2}r_{\texttt{red}} + \binom{r_{\texttt{red}}}{2}l_{\texttt{blue}} + \binom{r_{\texttt{blue}}}{2}l_{\texttt{red}} \;.
\end{equation}
%
\subparagraph*{Subquadratic Algorithms.}
To reduce the time, Brodal \ietal\cite{bmc13} applied the \textit{smaller half trick}, which specifies a depth first search order to visit the nodes $u$ of $T_{1}$, so that each leaf in $T_{1}$  changes color at most $\mathrm{O}(\log n)$ times. To count shared triplets efficiently without scanning $T_{2}$ completely for each node $u$ in $T_{1}$, $T_{2}$ is stored in a data structure denoted a \emph{hierarchical decomposition tree} (\emph{HDT}). This HDT maintains for the current visited node $u$ in $T_{1}$ the sum $\sum_{v \in T_{2}} \lvert {s(u) \cap s(v) \rvert}$, so that each color change in $T_{1}$ can be updated efficiently in $T_{2}$. In  \cite{bmc13} the HDT is a binary tree of height $\mathrm{O}(\log n)$ and every update can be done in a leaf to root path traversal in the HDT, which in total gives $\mathrm{O}(n \log^2 n)$ time. In \cite{TripletQuartetSoda13} the HDT is generalized to also handle non-binary trees, each query operates the same, and now due to a contraction scheme of the HDT the total time is reduced to $\mathrm{O}(n \log n)$. Finally, in \cite{Jansson2015} as an HDT the so called \emph{heavy-light tree decomposition} is used. Note that the only difference in all $\mathrm{O}(n \operatorname{polylog} n)$ results that are available until now is the type of HDT used.

In terms of external memory efficiency, every $\mathrm{O}(n \operatorname{polylog} n)$ algorithm performs $\mathrm{\Theta}(n \log n)$ updates to an HDT data structure, which means that for sufficiently large input trees every algorithm requires $\mathrm{\Omega}(n \log n)$ I/Os.

\section{The New Algorithm for Binary Trees}
\label{sec:binary}

\subparagraph{Overview.}

We will use the $\mathrm{O}(n^2)$ algorithm described in Section \ref{sec:prevApproaches} as a basis. The main difference lies in the order that we visit the nodes of $T_{1}$ and how we process $T_{2}$ when we count. We propose a new order of visiting the nodes of $T_{1}$, which we find by applying a hierarchical decomposition on $T_{1}$. Every component in this decomposition corresponds to a connected part of $T_{1}$ and a contracted version of $T_{2}$. In simple terms, if $\Lambda$ is the set of leafs in a component of $T_{1}$, the contracted version of $T_{2}$ is a binary tree on $\Lambda$ that preserves the topologies induced by $\Lambda$ in $T_{2}$ and has size $\mathrm{O}(|\Lambda|)$. To count shared triplets, every component of $T_{1}$ has a representative node $u$ that we use to scan the corresponding contracted version of $T_{2}$ in order to find  $\sum_{v \in T_{2}} \lvert {s(u) \cap s(v) \rvert}$. Unlike previous algorithms, we do not store $T_{2}$ in a data structure. We process $T_{2}$ by contracting and counting, both of which can be done by scanning. At the same time, even though we apply a hierarchical decomposition on $T_{1}$, the only reason why we do so, is so we can find the order in which to visit the nodes of $T_{1}$. This means that we do not need to store $T_{1}$ in a data structure either. Thus, we completely remove the need of data structures and scanning becomes the basic primitive in the algorithm. To make our algorithm I/O efficient, all that remains to be done is use a proper layout to store the trees in memory, so that every time we scan a tree of size $s$ we spend $\mathrm{O}(s/B)$ I/Os.

\subparagraph*{Preprocessing.}

As a preprocessing step, first we make $T_{1}$ left heavy, by swapping children so that for every node $u$ in $T_{1}$ the left subtree is larger than the right subtree, by a depth first traversal. Second, we change the leaf labels of $T_{1}$, which can also be done by a depth first traversal of $T_{1}$, so that the leafs are numbered $1$ to $n$ from left to right. This step takes $\mathrm{O}(n)$ time in the RAM model.

We do this to simplify the process of transferring the leaf colors between $T_{1}$ and $T_{2}$. The coloring of a subtree in $T_{1}$ will correspond to assigning the same color to a contiguous range of leaf labels. Determining the color of a leaf in $T_{2}$ will then require one if statement to find in what range (red or blue) its label belongs to.

\subparagraph*{Centroid Decomposition.}

For a given rooted binary tree $T$ we let $|T|$ denote the number of nodes in $T$ (internal nodes and leafs). For a node $u$ in $T$ we let $l$ and $r$ be the left and right children of $u$, and $p$ the parent. Removing $u$ from $T$ partitions $T$ into three (possibly empty) \emph{connected components} $T_{l}$, $T_{r}$ and $T_{p}$ containing $l$, $r$ and $p$, respectively. The \emph{centroid} is a node $u$ in $T$ such that $\max\{|T_{l}|, |T_{r}|, |T_{p}|\} \leq |T|/2$. A centroid always exists and can be found by starting from the root of $T$ and iteratively visiting the child with a largest subtree, eventually we will reach a centroid. Finding the size of every subtree and identifying $u$ takes $O(|T|)$ time in the RAM model. By recursively finding centroids in each of the three components, we will in the end get a ternary tree of centroids, which is called the \emph{centroid decomposition} of $T$, denoted $\mathit{CD}(T)$. We generate a level of $\mathit{CD}(T)$ in $\mathrm{O}(|T|)$ time. Since we have to generate at most $1+\log_{2}(|T|)$ levels, the total time required to build $\mathit{CD}(T)$ is $\mathrm{O}(|T|\log |T|)$, hence we get Lemma \ref{lemma1}.

\begin{lemma}
For any rooted binary tree $T$ with $n$ leafs, building $\mathit{CD}(T)$ takes $\mathrm{O}(n \log n)$ time in the RAM model.
\label{lemma1}
\end{lemma}

\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.50\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{cdExample/11.pdf}
       \caption{Picking the first centroid in $T_{1}$.}
    \end{subfigure}%
    \begin{subfigure}{0.50\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{cdExample/12.pdf}
        \caption{Recursing to component defined by $c_{1r}$.}
    \end{subfigure}\\
     \begin{subfigure}{0.50\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{cdExample/13.pdf}
        \caption{Recursing to component defined by $c_{2p}$.}
    \end{subfigure}%
     \begin{subfigure}{0.50\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{cdExample/14.pdf}
        \caption{Recursing to component defined by $c_{3p}$.}
    \end{subfigure}\\
 
    \caption{Generating a component in $\mathit{CD}(T_{1})$ that has two edges from below. The black polygon is the component.}
    \label{fig:CDExample}
\end{figure}

A component in a centroid decomposition $\mathit{CD}(T)$, might have many edges crossing its boundaries (connecting nodes inside and outside the component). An example with a component that has two edges from below can be seen in Figure \ref{fig:CDExample}. It is trivial to see that by following the same pattern of generating components we can have a component with an arbitrary number of edges from below.
The below \emph{modified centroid decomposition}, denoted $\mathit{MCD}(T)$, generates components with at most two edges crossing the boundary, one going towards the root and one down to exactly one subtree.

\subparagraph*{Modified Centroid Decomposition.}

A $\mathit{MCD}(T)$ is built recursively as follows: If a component $C$ has no edge from below, we select the centroid $c$ of $C$ as described above. Otherwise, let $(x,y)$ be the edge that crosses the boundary from below, where $x$ is in $C$ and let $c$ be centroid of $C$. As a splitting node choose the lowest common ancestor of $x$ and $c$. By induction every component has at most one edge from below and one edge from above. A useful property of $\mathit{MCD}(T)$ is presented by the following lemma:

\begin{lemma}
For any rooted binary tree $T$ with $n$ leafs, we have that $h(\mathit{MCD}(T))$ $\leq$ $2 + 2 \log_{2}n$, where $h(\mathit{MCD}(T))$ denotes the height of $\mathit{MCD}(T)$.
\label{lemma2}
\end{lemma}

\begin{proof}
In $\mathit{MCD}(T)$ if a component $C$ does not have an edge from below then the centroid of $C$ is used as a splitting node, thus generating three components $C_{l}$, $C_{r}$ and $C_{p}$ such that $|C_{l}| \leq \frac{|C|}{2}$, $|C_{r}|\leq \frac{|C|}{2}$ and $|C_{p}|\leq \frac{|C|}{2}$. Otherwise, $C$ has one edge $(x,y)$ from below, with $x$ being the node that is part of $C$. Let $c$ be the centroid of $C$. We have to consider the following two cases: if $c$ happens to be the lowest common ancestor of $c$ and $x$, then our algorithm will split $C$ according to the actual centroid, so we will have that $|C_{l}| \leq \frac{|C|}{2}$, $|C_{r}|\leq \frac{|C|}{2}$ and $|C_{p}|\leq \frac{|C|}{2}$. Otherwise, the splitting node will produce components $C_{l}$, $C_{r}$, and $C_{p}$ such that $|C_{l}| + |C_{p}| \leq \frac{|C|}{2}$ and $|C_{r}| \geq \frac{|C|}{2}$. From the first inequality we have that $|C_{l}| \leq \frac{|C|}{2}$ and $|C_{p}| \leq \frac{|C|}{2}$.  Notice however that $C_{r}$ is going to be a component corresponding to a complete subtree of $T$, so it will have no edges from below. This means that in the next recursion level when working with $C_{r}$ the actual centroid of $C_{r}$ will be chosen as a splitting node, then in the next level the three components produced from $C_{r}$ will be such that their sizes are at most half the size of $C$. From the analysis given so far, it becomes clear that when we have a component of size $|C|$ with one edge from below, then we will need at most 2 levels in $\mathit{MCD}(T)$ before producing components all of which will have a guaranteed size at most $\frac{|C|}{2}$. Hence we proved the given statement.
\end{proof}

It follows,

\begin{theorem}
For any rooted binary tree $T$ with $n$ leafs, building $\mathit{MCD}$($T$) takes $\mathrm{O}(n \log n)$ time in the RAM model.
\label{theorem1}
\end{theorem}

To return to our original problem, we visit the nodes of $T_{1}$, given by the depth first traversal of the ternary tree $\mathit{MCD}(T_{1})$, where the children of every node $u$ in $\mathit{MCD}(T_{1})$ are visited from left to right. 

For every such node $u$ we process $T_{2}$ in two phases, the \emph{contraction} phase and the \emph{counting} phase. 


\subparagraph*{Contraction.}

Let $L(T_{2})$ denote the set of leafs in $T_{2}$ and  $\Lambda \subseteq L(T_{2})$. In the contraction phase, $T_{2}$ is compressed into a binary tree of size $\mathrm{O}(|\Lambda|)$ whose leaf set is $\Lambda$.  The contraction is done in a way so that all the topologies induced by $\Lambda$ in $T_{2}$ are preserved in the compressed binary tree. This is achieved by the following three sequential steps:

\begin{itemize}
\item Prune all leafs of $T_{2}$ that are not in $\Lambda$,
\item Repeatedly prune all internal nodes of $T_{2}$ with no children, and
\item Repeatedly contract unary internal nodes, i.e. nodes having exactly one child.
\end{itemize}

Let $u$ be a node of $\mathit{MCD}(T_{1})$ and $C_{u}$ the corresponding component of $T_{1}$. For every such node $u$ we have a contracted version of $T_{2}$, from now on referred to as $T_{2}(u)$, where $L(T_{2}(u)) = L(C_{u})$. The goal is to augment $T_{2}(u)$ with counters (see counting phase), so that we can find  $\sum_{v \in T_{2}} \lvert {s(u) \cap s(v) \rvert}$ by scanning $T_{2}(u)$. One can imagine $\mathit{MCD}(T_{1})$ as being a tree where each node $u$ is augmented with $T_{2}(u)$. To generate all contractions of level $i$, which correspond to a set of disjoint connected components in $T_{1}$, we can reuse the contractions of level $i-1$. This means that we have to spend $\mathrm{O}(n)$ time to generate the contractions of level $i$, so to generate all contractions we need $\mathrm{O}(n \log n)$ time. Note that by explicitly storing all contractions, we will also need to use $\mathrm{O}(n \log n)$ space. For our problem, we traverse $\mathit{MCD}(T_{1})$ in a depth first search manner, so we only have to store a stack of contractions corresponding to the stack of nodes of $\mathit{MCD}(T_{1})$ that we have to remember during our traversal. Lemma \ref{lemma3} states that the size of this stack is always $\mathrm{O}(n)$.

\begin{lemma}
Let $T_{1}$ and $T_{2}$ be two rooted binary trees with $n$ leafs and $u_{1}, u_{2},...,u_{k}$ a root to leaf path of $\mathit{MCD}(T_{1})$. For the corresponding contracted versions $T_{2}(u_{1}), T_{2}(u_{2}), ..., T_{2}(u_{k})$ we have that $\sum_{i=1}^{k} \lvert T_{2}(u_{i}) \rvert = \mathrm{O}(n)$.
\label{lemma3}
\end{lemma}

\begin{proof}
For the root $u_{1}$ we have that $T_{2}(u_{1}) = T_{2}$ so $|T_{2}(u_{1})| \leq 2n$. From the proof of Lemma \ref{lemma2} we have that for every component of size $x$, we need at most two levels in $\mathit{MCD}(T_{1})$ before producing components all of which will have a guaranteed size of at most $\frac{x}{2}$. This means that $\sum_{i=1}^{k} \lvert T_{2}(u_{i}) \rvert \leq 2n + 2n + \frac{2n}{2} + \frac{2n}{2} + \frac{2n}{4} + \frac{2n}{4} + ... + \frac{2n}{2^{i}} + \frac{2n}{2^{i}} + ... = 2\sum_{j=0}^{\infty}\frac{2n}{2^{j}} \leq 8n = \mathrm{O}(n)$.
\end{proof}

\subparagraph*{Counting.}

In the counting phase, we find $\sum_{v \in T_{2}} \lvert {s(u) \cap s(v) \rvert}$ by scanning $T_{2}(u)$ instead of $T_{2}$. This makes the total time of the algorithm in the RAM model $\mathrm{O}(n \log n)$. We consider the following two cases:

\begin{itemize}[listparindent=\parindent]

\item \textit{$C_{u}$ has no edges from below.}

In this case $C_{u}$ corresponds to a complete subtree of $T_{1}$. We act exactly like in the $\mathrm{O}(n^2)$ algorithm (Section \ref{sec:prevApproaches}) but now instead of scanning $T_{2}$ we scan $T_{2}(u)$.
\item \textit{$C_{u}$ has one edge from below.}

\begin{figure}[ht]
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.33\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{tripletsAnchored/contractExample11.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.33\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{tripletsAnchored/contractExample22.pdf}
    \end{subfigure}%
     \begin{subfigure}{0.33\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{tripletsAnchored/contractExample33.pdf}
    \end{subfigure}
    \caption{$\mathit{MCD}(T_{1}$): Triplets that can be anchored in $u$ with the leafs not being in the component $C_{u}$.}
    \label{fig:tripletsWithXu}
\end{figure}

\begin{figure}[ht]
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.50\textwidth}
        \centering
        \includegraphics[width=1\textwidth]{tripletT2edges/1.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.50\textwidth}
        \centering
        \includegraphics[width=1\textwidth]{tripletT2edges/2.pdf}
    \end{subfigure}%
    \caption{Contracted subtrees on edges in $T_{2}(u)$ and shared triplets rooted on contracted nodes.}
    \label{fig:T2uEdges}
\end{figure}

In this case $C_{u}$ does not correspond to a complete subtree of $T_{1}$, since the edge from below $C_{u}$, will point to a subtree $X_{u}$, that is located outside of $C_{u}$. Note that because in the preprocessing step $T_{1}$ was made to be left heavy, $X_{u}$ is always rooted in a node on the left most path from $u$. The leafs in $X_{u}$ are important because they can be used to form triplets that are anchored in $u$. An illustration can be found in Figure \ref{fig:tripletsWithXu}. Acting in the exact same manner as in the previous case is not sufficient because we need to count these triplets as well.

To address this problem, every edge $(p_{v}, v)$ in $T_{2}(u)$ between a node $v$ and its parent $p_{v}$, is augmented with some counters about the leafs from $X_{u}$ that were contracted away in $T_{2}$. For every such edge $(p_{v},v)$, let $s_{1},s_{2},...,s_{k}$ be the contracted subtrees rooted on the edge (see Figure \ref{fig:T2uEdges}). Every such subtree contains either leafs with no color or leafs from $X_{u}$ that have the color red (the color can not be blue because $T_{1}$ was made to be left heavy). For every node $v$ in $T_{2}(u)$ the counters that we have are the following:

\begin{itemize}
\item $v_{\texttt{red}}$: the total number of red leafs in the subtree of $v$ (including those coming from $X_{u}$).
\item $v_{\texttt{blue}}$: the total number of blue leafs in the subtree of $v$.
\item $v_{ts}$: the total number of red leafs in $s_{1},s_{2},...,s_{k}$.
\item $v_{ps}$: the total number of pairs of red leafs from $s_{1},s_{2},...,s_{k}$ such that each pair comes from the same contracted subtree, i.e. $\sum_{i=1}^{k} \binom{r_{i}}{2}$ where $r_{i}$ is the number of red leafs in $s_{i}$.
\end{itemize}

The number of shared triplets that are anchored in a non-contracted node $v$ of $T_{2}(v)$ can be found like in the $\mathrm{O}(n^2)$ algorithm using the counters $v_{\texttt{red}}$ and $v_{\texttt{blue}}$ in (\ref{eq:eq1}).

As for the number of shared triplets that are anchored in a contracted node on edge $(p_{v}, v)$, this value is exactly $\binom{v.\texttt{blue}}{2} \cdot v_{ts} + v_{\texttt{blue}} \cdot v_{ps}$.

\end{itemize}

\subparagraph*{Scaling to External Memory.} 

If we store $T_{1}$ in an array of size $2n-1$ by using a preorder layout, i.e. if a node $v$ is stored in position $p$, the left child of $v$ is stored in position $p+1$ and if $x$ is the size of the left subtree of $v$ the right child of $v$ is stored in position $p + x + 1$, making $T_{1}$ left heavy requires $\mathrm{O}(n/B)$ I/Os. The preprocessing step that changes the labels of the leafs in $T_{1}$ and $T_{2}$ can be done in $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$ I/Os with a cache oblivious sorting routine, e.g. using merge sort. Finding the splitting node of a component $C_{u}$ requires $\mathrm{O}(|C_{u}|/B)$ I/Os, so in total the number of I/Os spent processing $T_{1}$ becomes $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$.

We use the proof of Lemma \ref{lemma3} to initialize a large enough array that can fit the contractions that we need to remember while traversing $\mathit{MCD}(T_{1})$. This array is used as a stack that we use to push and pop the contractions of $T_{2}$. Each contraction of $T_{2}$ is stored in memory using a post order layout, i.e. if a node $v$ is stored in position $p$ and $y$ is the size of the right subtree of $v$, the left child of $v$ is stored in position $v-y-1$ and the right child of $v$ is stored in position $v-1$. By using a stack, counting and contracting $T_{2}(u)$ requires $\mathrm{O}(|T_{2}(u)|/B)$ I/Os, so the total number of I/Os spent processing $T_{2}$ becomes $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$ as well.

Overall, our algorithm requires $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$ I/Os in the cache oblivious model.

\section{The New Algorithm for General Trees}
\label{sec:general}

Unlike a binary tree, a general tree allows having internal nodes with an arbitrary number of children. By introducing four colors, we will show that we can count the shared triplets in the two trees. We start by describing the new $\mathrm{O}(n^2)$ algorithm for general trees. We will then show how we can use the same ideas presented in the previous section to reduce the time to $\mathrm{O}(n \log n)$.

\subsection{\texorpdfstring{$\mathrm{O}(n^2)$}~ Algorithm for General Trees.}

Let $t$ be a triplet with leafs $i$, $j$ and $k$ that is either a resolved triplet $ij|k$  or an unresolved triplet $ijk$, where $i$ is to the left of $j$ and for the triplet $ijk$, $k$ is also to the right of $j$. Let $w$ be the lowest common ancestor of $i$ and $j$ and $(w,c)$ the edge from $w$ to the child $c$ whose subtree contains $j$. We anchor $t$ in the edge $(w,c)$. Let $s'(w,c)$ to be a set containing all triplets anchored in edge $(w,c)$. For the number of triplets $S(T_{1},T_{2})$ we have 
%
$$S(T_{1},T_{2}) = \sum_{(u,c) \in T_{1}}\sum_{(v,c') \in T_{2}} \lvert {s'(u,c) \cap s'(v,c')} \rvert \;.$$ 
%
To count shared triplets efficiently we use the following coloring procedure: in $T_{1}$ color the leafs of every child subtree of $u$ to the left of $c$ red, the leafs of the child subtree of $c$ blue, the leafs of every child subtree to the right of $c$ green and give the color black to every other leaf. We then transfer this coloring to the leafs in $T_{2}$. For the resolved triplet $ij|k$, $i$ corresponds to the red color, $j$ corresponds to the blue color and $k$ corresponds to the black color. For the unresolved triplet $ijk$, $i$ corresponds to the red color, $j$ corresponds to the blue color and $k$ corresponds to the green color.

Suppose that the node $v$ in $T_{2}$ has $k$ children. We will compute all shared triplets that are anchored in the $k$ edges that point to the children of $v$ in $\mathrm{O}(k)$ time, giving us the $\mathrm{O}(n^2)$ total time that we want. In $v$ we have the following counters:

\begin{itemize}
\item $v_{\texttt{red}}$: total number of red leafs in the subtree of $v$.
\item $v_{\texttt{blue}}$: total number of blue leafs in the subtree of $v$.
\item $v_{\texttt{green}}$: total number of green leafs in the subtree of $v$.
\item $\overline{v}_{\texttt{black}}$: total number of black leafs not in the subtree of $v$.
\end{itemize}

While scanning the $k$ children edges of $v$ from left to right, for a child $c'$ that is the $m-th$ child of $v$, we want to maintain the following counters:

\begin{itemize}
\item $a_{\texttt{red}}$: total number of red leafs from the first $m-1$ children subtrees.
\item $a_{\texttt{blue}}$: total number of blue leafs from the first $m-1$ children subtrees.
\item $a_{\texttt{green}}$: total number of green leafs from the first $m-1$ children subtrees.
\item $p_{\texttt{red,green}}$: total number of pairs of leafs from the first $m-1$ children subtrees, where one is red, the other is green and they both come from different children subtrees.
\item $p_{\texttt{red,blue}}$ : total number of pairs of leafs from the first $m-1$ children subtrees, where one is red, the other is blue and they both come from different children subtrees.
\item $p_{\texttt{blue,green}}$ : total number of pairs of leafs from the first $m-1$ children subtrees, where one is blue, the other is green and they both come from different children subtrees.
\item $t_{\texttt{red,blue,green}}$: the total number of leaf triples from the first $m-1$ children subtrees, where one is red, one is blue, one is green, and all three leafs come from different children subtrees.
\end{itemize}

Before beginning the scanning process of the children edges of $v$, every variable is initialized to 0. Then for the child $c'$ every variable can then be updated in $\mathrm{O}(1)$ time as follows:

\begin{itemize}
\item $a_{\texttt{red}} = a_{\texttt{red}} + c'_{\texttt{red}}$
\item $a_{\texttt{blue}} = a_{\texttt{blue}} + c'_{\texttt{blue}}$
\item $a_{\texttt{green}} = a_{\texttt{green}} + c'_{\texttt{green}}$
\item $p_{\texttt{red,green}} = p_{\texttt{red,green}} +  a_{\texttt{green}} \cdot c'_{\texttt{red}} +  a_{\texttt{red}} \cdot c'_{\texttt{green}}$
\item $p_{\texttt{red,blue}} = p_{\texttt{red,blue}} +  a_{\texttt{blue}} \cdot c'_{\texttt{red}} +  a_{\texttt{red}} \cdot c'_{\texttt{blue}}$
\item $p_{\texttt{blue,green}} = p_{\texttt{blue,green}} +  a_{\texttt{green}} \cdot c'_{\texttt{blue}} +  a_{\texttt{blue}} \cdot c'_{\texttt{green}}$
\item $t_{\texttt{red,blue,green}} = t_{\texttt{red,blue,green}} + p_{\texttt{red,green}} \cdot c'_{\texttt{blue}} + p_{\texttt{red,blue}} \cdot c'_{\texttt{green}} + p_{\texttt{blue,green}} \cdot c'_{\texttt{red}} $
\end{itemize}

After finishing scanning the $k$ children edges of $v$, we can compute the shared triplets that are anchored in every child edge of $v$ as follows: for the total number of shared resolved triplets, denoted $\texttt{tot}_{\texttt{res}}$, we have that $\texttt{tot}_{\texttt{res}} = p_{\texttt{red,blue}} \cdot \overline{v}_{\texttt{black}}$ and for the total number of shared unresolved triplets, denoted $\texttt{tot}_{\texttt{unres}}$, we have that $\texttt{tot}_{\texttt{unres}} = t_{\texttt{red,blue,green}}$. Now we are ready to describe the $\mathrm{O}(n \log n)$ algorithm.

\subsection{\texorpdfstring{$\mathrm{O}(n\log n)$}~ Algorithm for General Trees.}

\subparagraph*{Preprocessing.}

\begin{figure}
    \centering
        \includegraphics[width=0.6\textwidth]{buildingBOFT1/2.pdf}
    \caption{Transforming $T_{1}$ into $b(T_{1})$.}
    \label{fig:buildingboft1}
\end{figure}

\begin{figure}[ht]
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}[b]{0.50\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{MCDinBOFT1/1.pdf}
    \end{subfigure}%
    \begin{subfigure}[b]{0.50\textwidth}
        \centering
        \includegraphics[width=0.95\textwidth]{MCDinBOFT1/2.pdf}
    \end{subfigure}\\
    \caption{How a component in $b(T_{1})$ translates to a component in $T_{1}$.}
    \label{fig:componentsGeneral}
\end{figure}

In the preprocessing step of the algorithm, we start by transforming $T_{1}$ into a binary tree, denoted $b(T_{1})$. Let $w$ be a node of $T_{1}$ that has exactly $k$ children. The $k$ edges that connect $w$ to its children will be replaced by a binary tree,  the leafs of which will be the $k$ children of $w$, and each internal node will have the color \emph{orange}. Every such binary tree will be built in a way so that every orange node is on the left most path that starts from $w$, and the left most path is part of the heavy path starting from node $w$ in $b(T_{1})$, i.e. $b(T_{1})$ is left heavy (see Figure \ref{fig:buildingboft1}). Note that the number of orange nodes added because of node $w$ is exactly $k-2 $ which is less than the number of $w$'s children in $T_{1}$, so the total number of nodes in $b(T_{1})$ is still $\mathrm{O}(n)$.

We can generate the colorings of the $\mathrm{O}(n^2)$ algorithm by using $b(T_{1})$ instead of $T_{1}$. Let $u$ be some internal node of $b(T_{1})$ and $c$ be the right child of $u$. The leafs of $b(T_{1})$ are colored according to edge $(u,c)$ as follows: the left subtree of $u$ is colored red, the right subtree of $u$ is colored blue. If $u$ is an orange node, then the remaining leafs in the orange binary tree that $u$ is part of are colored green. All other leafs of $b(T_{1})$ are colored black. 

Let $u_{\texttt{root}}$ be the root node of the orange tree that $u$ is part of. If $u$ is not part of an orange tree, then $u_{\texttt{root}}$ is the parent of $u$. The nodes $u_{\texttt{root}}$ and $c$ exist in $T_{1}$ and $(u_{\texttt{root}}, c)$ is an edge in $T_{1}$. So coloring the leafs of $b(T_{1})$ according to edge $(u,c)$ is equivalent to coloring the leafs of $T_{1}$ according to edge $(u_{\texttt{root}}, c)$. Every right edge of $b(T_{1})$ corresponds to exactly one edge of $T_{1}$, coloring the leafs of $b(T_{1})$ according to this right edge as described above will produce exactly the same coloring as the coloring according to the corresponding edge in $T_{1}$ as described in the $\mathrm{O}(n^2)$ algorithm. This implies correctness, i.e. in the $\mathrm{O}(n^2)$ algorithm using $b(T_{1})$ instead of $T_{1}$ will yield $S(T_{1},T_{2})$. 

For notation purposes, for the edge $(u,c)$ in $b(T_{1})$, we let $s''(u,c)$ be the set of triplets that are anchored in the edge $(u_{root},c)$ in $T_{1}$, i.e. $s''(u,c) = s'(u_{\texttt{root}},c)$. For the number of shared triplets we then have:
%
$$S(T_{1},T_{2}) = \sum_{(u,c) \in b(T_{1})}\sum_{(v,c') \in T_{2}} \lvert {s''(u,c) \cap s'(v,c')} \rvert \;.$$ 
%
The reason behind this transformation is because now we can use exactly the same core ideas described in Section \ref{sec:binary}. The tree $b(T_{1})$ is a binary tree, so we apply the same preprocessing step, except we do not make it left heavy because by construction it already is. However, we do change the labels of the leafs in $b(T_{1})$ and $T_{2}$, so that the leafs in $b(T_{1})$ are numbered $1$ to $n$ from left to right. 

\subparagraph*{Modified Centroid Decomposition.}

After the preprocessing step, we build $\mathit{MCD}(b(T_{1}))$ as described in Section \ref{sec:binary}. Then we traverse the nodes of $b(T_{1})$, given by a depth first traversal of $\mathit{MCD}(b(T_{1}))$, where the children of every node $u$ in $\mathit{MCD}(b(T_{1}))$ are visited from left to right.


In Figure \ref{fig:componentsGeneral} we can see that a component $C_{u}$ of $b(T_{1})$ structurally looks like a component of $T_{1}$ in the binary algorithm of Section \ref{sec:binary}. However, the edges crossing its boundary can now be orange edges as well, which in $T_{1}$ that now is a general tree, translates to more than one consecutive subtrees.

Like in the binary algorithm, while traversing  $\mathit{MCD}(b(T_{1}))$ we process $T_{2}$ in two phases, the contraction phase and the counting phase. 

\subparagraph*{Contraction.} 

The contraction of $T_{2}$ with respect to a set of leafs $\Lambda \subset L(T_{2})$, happens in the exact same way as described in Section \ref{sec:binary}, i.e. we start by pruning all leafs of $T_{2}$ that are not in $\Lambda$, then we prune all internal nodes of $T_{2}$ with no children, and finally, we contract the nodes that have exactly one child. 

Let $u$ be a node of $\mathit{MCD}(b(T_{1}))$ and $C_{u}$ the corresponding component of $b(T_{1})$. For every such node $u$ we have a contracted version of $T_{2}$, denoted $T_{2}(u)$, where $L(T_{2}(u)) = L(C_{u})$. Just like in the binary algorithm, the goal is to augment $T_{2}(u)$ with counters, so that we can find  $\sum_{(v,c') \in T_{2}} \lvert {s''(u,c) \cap s'(v,c')} \rvert$ by scanning $T_{2}(u)$ instead of $T_{2}$. %From the proprocessing step, we have that $\sum_{(v,c') \in T_{2}} \lvert {s''(u,c) \cap s'(v,c')} \rvert = \sum_{(v,c') \in T_{2}} \lvert {s'(u_{\texttt{root},c}) \cap s'(v,c')} \rvert$.

Because of the location where the triplets are anchored, in $T_{2}(u)$ every leaf that was contracted away, must have a color and be stored in some way. The color of each leaf depends on the type of the component that we have in $b(T_{1})$ and the splitting node that is used for that component. For example, in Figure \ref{fig:componentsGeneral} the contracted leafs from $X_{u}$ will have the red color because like in the binary algorithm $b(T_{1})$ is left heavy. The contracted leafs from the children subtrees of $u_{p}$ in $T_{1}$ can either have the color green or black. If $u$ in $b(T_{1})$ happens to be part of the orange tree that $u_{p}$ is the root of, then the color must be green, otherwise black. Finally, every leaf that is not in the subtree defined by $u_{p}$, and thus is in $Y_{u}$, must have the color black. The way we store this information is described in the counting phase.

\subparagraph*{Counting.}
\begin{figure}[ht]
    \centering
    \includegraphics[width=0.7\textwidth]{generalContractedNodesEdges/1.pdf}
    \caption{$T_{2}(u)$: Contracted children subtrees on node $v$ and contracted subtrees on contracted nodes (gray color) in edge $(p_{v},v)$.}
    \label{fig:T2contractedNodesEdgesGeneral}
\end{figure}

In Figure \ref{fig:T2contractedNodesEdgesGeneral} we illustrate how a node $v$ in $T_{2}(u)$ can look like. The contracted subtrees are illustrated with the dark gray color. Every such subtree contains some number of red, green and black leafs. The counters that we need to maintain should be so that if $v$ has $k$ children in $T_{2}(u)$, then we can count all shared triplets that are anchored in every child edge (including those of the contracted children subtrees) of $v$ in $\mathrm{O}(k)$ time. At the same time, in $\mathrm{O}(1)$ time we should be able to count all shared triplets that are anchored in every child edge of every contracted node that lies on edge $(p_{v}, v)$. In this way, the counting phase will require $\mathrm{O}(|T_{2}(u)|)$ time, hence we will get the same bounds like in the binary algorithm.

In $v$ we have the following counters:


\begin{itemize}
\item $v_{\texttt{red}}$: total number of red leafs (including the contracted leafs) in the subtree of $v$.
\item $v_{\texttt{blue}}$: total number of blue leafs in the subtree of $v$.
\item $v_{\texttt{green}}$: total number of green leafs (including the contracted leafs) in the subtree of $v$.
\item $\overline{v}_{\texttt{black}}$: total number of black leafs (including the contracted leafs) not in the subtree of $v$.
\end{itemize}

We divide the rest of the counters into two categories. The first category corresponds to the leafs in the contracted children subtrees of $v$ and each counter will be stored in a variable of the form $v_{A.\{\}}$. The second category corresponds to the leafs in the subtrees in edge $(p_{v}, v)$ and each counter will be stored in a variable of the form $v_{B.\{\}}$. 

For the first category we have the following counters:
\begin{itemize}
\item $v_{A.\texttt{red}}$: total number of red leafs in the contracted children subtrees of $v$.
\item $v_{A.\texttt{green}}$: total number of green leafs in the contracted children subtrees of $v$.
\item $v_{A.\texttt{black}}$: total number of black leafs in the contracted children subtrees of $v$.
\item $v_{A.\texttt{red,green}}$: total number of pairs of leafs where one is red, the other is green and one leaf comes from one contracted child subtree of $v$ and the other leaf comes from a different contracted child subtree of $v$.
\end{itemize}

While scanning the $k$ edges of the children of $v$ from left to right, for a child $c'$ that is the $m-th$ child of $v$, we also maintain the following:

\begin{itemize}
\item $a_{\texttt{red}}$: total number of red leafs from the first $m-1$ children subtrees, including the contracted children subtrees.
\item $a_{\texttt{blue}}$: total number of blue leafs from the first $m-1$ children subtrees.
\item $a_{\texttt{green}}$: total number of green leafs from the first $m-1$ children subtrees, including the contracted children subtrees.
\item $p_{\texttt{red,green}}$: total number of pairs of leafs from the first $m-1$ children subtrees and all contracted children subtrees, where one is red, the other is green and they both come from different subtrees (one might be contracted and the other non-contracted).
\item $p_{\texttt{red,blue}}$ : total number of pairs of leafs from the first $m-1$ children subtrees and all contracted children subtrees, where one is red, the other is blue and they both come from different subtrees (one might be contracted and the other non-contracted).
\item $p_{\texttt{blue,green}}$ : total number of pairs of leafs from the first $m-1$ children subtrees and all contracted children subtrees, where one is blue, the other is green and they both come from different subtrees (one might be contracted and the other non-contracted). 
\item $t_{\texttt{red,blue,green}}$: the total number of leaf triples from the first $m-1$ children subtrees and all contracted children subtrees, where one is red, one is blue and one is green, and all three leafs come from different subtrees (some might be contracted, some might be non-contracted).
\end{itemize}

For every child $c'$, each variable can be updated in $\mathrm{O}(1)$ in exactly the same manner like in the $\mathrm{O}(n^2)$ algorithm. The main difference is the values of the variables before we begin scanning the children edges of $v$. Each variable is initialized as follows:

\begin{itemize}
\item $a_{\texttt{red}} = v_{A.\texttt{red}}$
\item $a_{\texttt{blue}} = 0$
\item $a_{\texttt{green}} = v_{A.\texttt{green}}$
\item $p_{\texttt{red,green}} = v_{A.\texttt{red,green}}$
\item $p_{\texttt{red,blue}} = p_{\texttt{blue,green}} = t_{\texttt{red,blue,green}} = 0$
\end{itemize}

After finishing scanning the $k$ children edges of $v$, we can compute the shared triplets that are anchored in every child edge of $v$ (including the children edges pointing to contracted subtrees) as follows: for the total number of shared resolved triplets, denoted $\texttt{tot}_{\texttt{A.res}}$, we have that $\texttt{tot}_{\texttt{A.res}} = p_{\texttt{red,blue}} \cdot \overline{v}_{\texttt{black}}$ and for the total number of shared unresolved triplets, denoted $\texttt{tot}_{\texttt{A.unres}}$, we have that $\texttt{tot}_{\texttt{unres}} = t_{\texttt{red,blue,green}}$.

The second category of counters will help us count triplets involving leafs (contracted and non-contracted) from the subtree of $v$ and leafs from the contracted subtrees in edge $(p_{v},v)$. We maintain the following:

\begin{itemize}
\item $v_{B.\texttt{red}}$: total number of red leafs in all contracted subtrees in edge $(p_{v},v)$.
\item $v_{B.\texttt{green}}$: total number of green leafs in all contracted subtrees in edge $(p_{v},v)$.
\item $v_{B.\texttt{black}}$: total number of black leafs in all contracted subtrees in edge $(p_{v},v)$.
\item $v_{B.\texttt{red,green}}$: total number of pairs of leafs where one is red and the other is green such that one leaf comes from one contracted child subtree of a contracted node $v'$ and the other leaf comes from a different contracted child subtree of the same contracted node $v'$.
\item $v_{B.\texttt{red,black}}$: total number of pairs of leafs where one is red and the other is black such that the red leaf comes from one contracted child subtree of a contracted node $v'$ and the black leaf comes from one contracted child subtree of a contracted node $v''$. For $v'$ and $v''$ we have that $v''$ is closer to $v_{p}$ than $v'$.
\end{itemize}

For the total number of shared unresolved triplets, denoted $\texttt{tot}_{B.\texttt{unres}}$, that are anchored in the children edges of every contracted node that exists in edge $(v_{p},v)$, we have that $\texttt{tot}_{B.\texttt{unres}} = v_{\texttt{blue}} \cdot v_{B.\texttt{red,green}}$. For the total number of shared resolved triplets, denoted $\texttt{tot}_{B.\texttt{res}}$, that are anchored in the children edges of every contracted node that exists in edge $(v_{p},v)$, we have that  $\texttt{tot}_{B.\texttt{res}} = v_{\texttt{blue}} \cdot v_{B.\texttt{red,black}} + v_{\texttt{blue}} \cdot v_{B.\texttt{red}} \cdot (\overline{v}_{\texttt{black}} - v_{B.\texttt{black}})$.

\subparagraph*{Scaling to External Memory.} 

Just like in Section \ref{sec:binary}, we store $T_{1}$ in an array of size $2n-1$ by using a preorder layout, that is if a node $v$ is stored in position $p$, the left child of $v$ is stored in position $p+1$ and if $x$ is the size of the left subtree of $v$ the right child of $v$ is stored in position $p + x + 1$. We build $b(T_{1})$ in $\mathrm{O}(n/B)$ I/Os and store it in a different array following the same layout as $T_{1}$.  Making $b(T_{1})$ left heavy requires $\mathrm{O}(n/B)$ I/Os and the preprocessing step that changes the labels of the leafs in $b(T_{1})$ and $T_{2}$ can be done in $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$ I/Os with a cache oblivious sorting routine, e.g. using merge sort. Finding the splitting node of a component $C_{u}$ requires $\mathrm{O}(|C_{u}|/B)$ I/Os, so in total the number of I/Os spent from processing $T_{1}$ and $b(T_{1})$ becomes $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$.

The proof of Lemma \ref{lemma3} can be trivially modified to apply to general trees as well, so we use it to initialize a large enough array that can fit the contractions that we need to remember while traversing $\mathit{MCD}(b(T_{1}))$. This array is used as a stack that we use to push and pop the contractions of $T_{2}$. Each contraction of $T_{2}$ is stored in memory using a post order layout, i.e. if a node $v$ is stored in position $p$ and $y$ is the size of the right subtree of $v$, the left child of $v$ is stored in position $v-y-1$ and the right child of $v$ is stored in position $v-1$. By using a stack, counting and contracting $T_{2}(u)$ requires $\mathrm{O}(|T_{2}(u)|/B)$ I/Os. 

Overall our algorithm for general trees requires $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$ I/Os in the cache oblivious model.

\section{Implementation}
\label{sec:implementation}

\begin{figure}[ht]
    \centering
        \includegraphics[width=0.6\textwidth]{impl/binaryOverview/1.pdf}
    \caption{Implementation overview.}
    \label{fig:implOverview}
\end{figure}

The implementation of both algorithms was made using the C++ programming language. A high level overview is illustrated in Figure \ref{fig:implOverview}. The source code can be found in \url{https://github.com/kmampent/CacheTD}.


\subsection{Input}

The two input trees $T_{1}$ and $T_{2}$ are stored in two separate text files following the Newick format. Both trees have $n$ leafs and the label of each leaf is assumed to be a number in \{1,2,...,n\}. Two leafs can not have the same label.

\subsection{Parser}

The parser receives the files that store $T_{1}$ and $T_{2}$ in Newick format, and returns $T_{1}$ and $T_{2}$ but now with $T_{1}$ stored inside an array following the preorder layout and $T_{2}$ inside an array following the postorder layout. The parser requires $\mathrm{O}(n)$ time and $\mathrm{O}{(n/B)}$ I/Os.

\subsection{Algorithm}

Having $T_{1}$ and $T_{2}$ stored in the right layouts, we can now proceed with the main part of the algorithm. Both implementations (binary, general) follow the same approach. There exists an \emph{initialization} phase and a \emph{distance computation} phase.

\subparagraph*{Initialization.}

In the initialization phase the first component of $T_{1}$ is built, and the corresponding contracted version of $T_{2}$, from now on referred to as \emph{corresponding component} of $T_{2}$, is built as well. At the end of this phase, the first component of $T_{1}$ will be stored in the preorder layout in a new array (different than the one produced by the parser) and similarly the corresponding component of $T_{2}$, which however will be stored in the postorder layout instead. 

Note that the first component of $T_{1}$ is $T_{1}$ and the component of $T_{2}$ is $T_{2}$. The main difference lies in the representation of the nodes of $T_{1}$ and $T_{2}$ that now have to store different variables/counters. The reason why we do not immediately build these components while parsing the input, is because we want to make the experiments as fair as possible. The main focus in the experiments is comparing the algorithms and not the different parser implementations. All previous implementations after parsing the input, stored the trees in a layout format and then proceeded with the main part of the algorithm. We do exactly the same, so that only the time of the initialization phase contributes to the total execution time of the algorithm, and not the parser.

\subparagraph*{Distance Computation.}

At the end of the initialization phase, the first component of $T_{1}$, denoted $\texttt{comp}(T_{1})$, and the corresponding component of $T_{2}$, denoted $\texttt{comp}(T_{2})$, are available. We are ready to start counting shared triplets to compute $S(T_{1}, T_{2})$. The following steps are recursively applied:

\begin{itemize}
\item Starting from the root of $\texttt{comp}(T_{1})$, scan the left most path of $\texttt{comp}(T_{1})$ to find the splitting node $u$.
\item Scan $\texttt{comp}(T_{2})$ to compute for the binary algorithm $\sum_{v \in T_{2}} \lvert {s(u) \cap s(v) \rvert}$, or for the general algorithm $\sum_{(v,c') \in T_{2}} \lvert {s''(u,c) \cap s'(v,c')} \rvert$.
\item Using the splitting node $u$, generate the next three components of $\texttt{comp}(T_{1})$.  Call the component determined by the left child of $u$ $\texttt{comp}(T_{1}(u_{l}))$, the component determined by the right child of $u$ $\texttt{comp}(T_{1}(u_{r}))$ and the component determined by the parent of $\texttt{comp}(T_{1})$ $\texttt{comp}(T_{1}(u_{p}))$.
\item Scan and contract $\texttt{comp}(T_{2})$ to generate the three corresponding components, denoted $\texttt{comp}(T_{2}(u_{l}))$, $\texttt{comp}(T_{2}(u_{r}))$ and $\texttt{comp}(T_{2}(u_{p}))$, which are contracted versions of $T_{2}$ with all the necessary counters (see Section \ref{sec:binary} and Section \ref{sec:general}) properly maintained.
\item Recurse on the three pairs of components, one from $T_{1}$ and one from $T_{2}$, and for each pair repeat the steps above.
\end{itemize}

As a last step, print $\binom{n}{3} - S(T_{1}, T_{2})$, which is the distance $D(T_{1},T_{2})$ that we wanted to find.

\subparagraph*{Correctness.}

The correctness was extensively tested by generating hundreds of thousands of random trees of varying size and varying degree and comparing the output of our implementations against the output of the implementations of the $\mathrm{O}(n \log^{3} n)$ time algorithm in \cite{Jansson2015} and the $\mathrm{O}(n \log n)$ time algorithm in \cite{bioinformatics14}. 


\subparagraph*{Changing the Leaf Labels.} To get the right theory bounds, changing the leaf labels of $T_{1}$ and $T_{2}$ must be done with a cache oblivious sorting routine, e.g. merge sort. In the RAM model this approach requires $\mathrm{O}(n \log n)$ time and in the cache oblivious model $\mathrm{O}(\frac{n}{B}\log_{2} \frac{n}{M})$ I/Os. A second approach is to exploit the fact that each label is between $1$ and $n$ and use an auxiliary array that stores the new labels of the leafs in $T_{1}$, which we then use to update the leaf labels of $T_{2}$. In the RAM model the second approach requires $\mathrm{O}(n)$ time but in the cache oblivious model $\mathrm{O}(n)$ I/Os. In practice the problem with the first approach is that the number of instructions that it incurs eliminates any advantage that we would expect to get due to its cache related efficiency for $L_{1}$, $L_{2}$ and $L_{3}$ cache. For the input sizes tested, the array of labels easily fits into RAM, so in our implementation of both algorithms we use the second approach. 

\section{Experiments}
\label{sec:experiments}

In this section we present experiments illustrating the practical performance of the algorithms described in Sections \ref{sec:binary} and \ref{sec:general}.

\subsection{The Setup}

The experiments were performed on a machine with 8GB RAM, Intel(R) Core(TM) i5-3470 CPU @ 3.20GHz, 32K L1 cache, 256K L2 cache and 6144K L3 cache. The operating system was Ubuntu, version 16.04.2 LTS. The compiler that was used was g++ version 5.4 and cmake version 3.5.1.

\subparagraph*{Generating Random Trees.}

We use two different models for generating input trees. The first model is called the \emph{random model}. A tree $T$ with $n$ leafs in this model is generated as follows: 

\begin{itemize}
\item Create a binary tree with $n$ leafs by the following procedure:
\begin{itemize}
\item Start with a binary tree $T'$ with two leafs.
\item Iteratively pick a leaf $l$ uniformly at random.
\item Make $l$ an internal node by appending a left child node and a right child node to $l$, thus increasing the number of leafs in $T'$ by exactly 1.
\end{itemize}
\item With probability $p$ contract every internal node $u$ of $T'$, i.e make the children of $u$ be the children of $u$'s parent and remove $u$.
\end{itemize}

The second model is called \emph{alpha} model. In this model, we can control more directly the shape of the input trees. A tree $T$ with $n$ leafs in this model is generated as follows: 

\begin{itemize}
\item Create a binary tree $T'$ with $n$ leafs as follows: let $0 \leq \alpha \leq 1$ be a parameter, $u$ some internal node in $T$, $l$ and $r$ the left and right children of $u$, and $T(u)$, $T(l)$ and $T(r)$ the subtrees rooted on $u$, $l$ and $r$ respectively. Create $T'$ so that for every internal node $u$ the fraction $\frac{|T(l)|}{|T(u)|}$ approaches $\alpha$, more precisely $|T_{l}| = \max(1, \min(\floor{a\cdot n},n-1))$ and $|T_{r}| = 1 - |T_{l}|$ where $|T_{l}|$ and $|T_{r}|$ are the number of leafs in $T(l)$ and $T(r)$ respectively.
\item With probability $p$ contract every internal node $u$ of $T'$, i.e make the children of $u$ be the children of $u$'s parent and remove $u$.
\end{itemize}

Note that in both models, after creating a tree, the leaf labels are shuffled by using \texttt{std::shuffle}\footnote{\url{http://www.cplusplus.com/reference/algorithm/shuffle/}} together with \texttt{std::default\_random\_engine}\footnote{\url{http://www.cplusplus.com/reference/random/default_random_engine/}}.

\subparagraph*{Implementations Tested.}

Let $p_{1}$ and $p_{2}$ denote the contraction probability of $T_{1}$ and $T_{2}$ respectively. When $p_{1} = p_{2} = 0$, the trees $T_{1}$ and $T_{2}$ are binary trees, so in our experiments we use the algorithm from Section \ref{sec:binary}. In all other cases, the algorithm from Section \ref{sec:general} is used. Note that the algorithm from Section \ref{sec:general} can handle binary trees just fine, however there is a natural overhead compared to the algorithm from Section \ref{sec:binary} that comes due to the additional counters that we have to maintain for the contractions of $T_{2}$.

We compared our implementation with previous implementations of \cite{Jansson2015} and \cite{bmc13, TripletQuartetSoda13} available at \url{http://sunflower.kuicr.kyoto-u.ac.jp/~jj/Software/} and \url{http://users-cs.au.dk/cstorm/software/tqdist/} respectively. The implementation of the $\mathrm{O}(n \log ^3 n)$ algorithm in \cite{Jansson2015} has two versions, one that uses \texttt{unordered\_map}\footnote{\url{http://en.cppreference.com/w/cpp/container/unordered_map}}, which we refer to as \texttt{CPDT}, and another that uses sparsehash\footnote{\url{https://github.com/sparsehash/sparsehash}}, which we refer to as \texttt{CPDTg}. The tqdist library \cite{bioinformatics14}, which we will refer to as \texttt{tqDist}, has an implementation of the binary $\mathrm{O}(n \log^2 n)$ algorithm from \cite{bmc13} and the general  $\mathrm{O}(n \log n)$ algorithm from \cite{TripletQuartetSoda13}. If the two input trees are binary the $\mathrm{O}(n \log^2 n)$ algorithm is used. We will refer to our new algorithm as \texttt{CacheTD}.

\subparagraph*{Statistics.}

We measured the execution time of the algorithms with the \texttt{clock\_gettime} function in C++. We used the PAPI library\footnote{\url{http://icl.utk.edu/papi/}} for statistics related to L1, L2 and L3 cache accesses and misses. Finally, we count the space of the algorithms by considering the \emph{Maximum resident set size} returned by \texttt{/usr/bin/time -v}.

\subsection{Results}

The experiments are divided into two parts. In the first part we look at how the different algorithms perform when the memory requirements do not exceed the available main memory (8G RAM).  In the second part we look at how the different algorithms perform when the memory requirements exceed the available main memory (by limiting the available RAM to the operating system to be 1GB), thus forcing the operating system to use the swap space which in turn yields the very expensive disk I/Os. All figures can be found in Appendix \ref{appendix:experimentFigures}.

\subparagraph*{RAM experiments.}

For the random model, in Figure \ref{fig:expTimeSmall} we illustrate a time comparison of all implementations for trees of up to $2^{21}$ leafs ($\sim 2$ million) with varying contraction probabilities. Every data point is the average of 10 runs. In all cases the new algorithms achieve the best performance. We note that for the case where $p_{1} = 0.95$ and $p_{2} = 0.2$, \texttt{CPDT} behaves in a different way compared to the experiments in \cite{Jansson2015}. The reason is because of the differences in the implementation of \texttt{unordered\_map} that exist between the different versions of the g++ compilers. For example, if we compile with GCC 4.7, \texttt{CPDT} will produce the same performance as in \cite{Jansson2015} (see Figure \ref{fig:gpp}). In Figure \ref{fig:expSpaceSmall} we compare the space consumption of the algorithms. \texttt{CacheTD} is the only algorithm that uses $\mathrm{O}(n)$ space for both binary and general trees. In theory we expect that the space consumption is better and this is also what we get in practice. In Figures \ref{fig:expPtimeSmall} and \ref{fig:expPspaceSmall} we can see how the contraction parameter affects the running time and the space consumption of the algorithms respectively. Finally, in Figures \ref{fig:l1misses}, \ref{fig:l2misses} and \ref{fig:l3misses} we compare the cache performance of the algorithms, i.e. how many cache misses (L1, L2 and L3 respectively) the algorithms perform for increasing input sizes and varying contraction parameters. As expected, the new algorithm achieves a significant improvement over all previous algorithms.

For the alpha model, the main interesting experimental results are illustrated in Figure \ref{fig:alphaFigure}, where we plot the alpha parameter against the execution time of the algorithms, when $n = 2^{21}$. The alpha parameter has the least effect on \texttt{CacheTD}, with the maximum running time being only a factor of 1.1 larger than the minimum. As mentioned in Section \ref{sec:prevApproaches}, \texttt{CPDT} and \texttt{CPDTg} use the heavy light decomposition for $T_{2}$. When $\alpha$ approaches 0 or 1, the number of heavy paths that will be updated because of a leaf color change decreases, thus the total number of operations of the algorithm decreases as well.

\subparagraph*{I/O experiments.}

\begin{table}[ht]
  \centering
  \caption{Random model: Time performance when limiting the available RAM to 1GB. For the left table we have $p_{1} = p_{2} = 0$ and for the right table $p_{1} = p_{2} = 0.5$.}
  \begin{tabular}{ccccc}
  	\toprule
  	$n$ &\texttt{CPDT/CPDTg}&\texttt{tqDist}& \texttt{CacheTD}\\
  	\midrule
    $2^{15}$ & 0m:01s & 0m:02s & 0m:01s\\
    $2^{16}$ & 0m:01s& 0m:02s & 0m:01s\\
    $2^{17}$ & 0m:01s & 0m:05s & 0m:01s\\
    $2^{18}$ & 0m:03s  & 2m:38s & 0m:01s\\
    $2^{19}$ & 0m:06s  & 1h:32m & 0m:01s\\
    $2^{20}$ & 0m:36s & - & 0m:02s\\
    $2^{21}$ & 19m:40s  & - & 0m:08s\\
    $2^{22}$ & $>$9h & - & 1m:13s\\
    $2^{23}$ & -  & - & 8m:12s\\
    $2^{24}$ & - & - & 31m:31s\\\hline
  \end{tabular}
  \begin{tabular}{ccccc}
  	\toprule
  	$n$ &\texttt{CPDT}&\texttt{CPDTg}&\texttt{tqDist}& \texttt{CacheTD}\\
  	\midrule
     $2^{15}$ &  0m:01s & 0m:01s &  0m:01s & 0m:01s\\
     $2^{16}$ & 0m:01s & 0m:01s &  0m:01s & 0m:01s\\
    $2^{17}$ &  0m:01s & 0m:01s &  0m:04s & 0m:01s\\
    $2^{18}$ & 0m:03s & 0m:03s &  0m:22s & 0m:01s\\
     $2^{19}$ & 0m:07s & 0m:08s &  16m:49s & 0m:01s\\
     $2^{20}$ &  29m:06s  &4h:36m& $>$10h & 0m:02s\\
    $2^{21}$ & >5h  &-& - & 1m:04s\\
     $2^{22}$ &  -  &-& - & 4m:58s\\
      $2^{23}$ & - &-& - & 27m:12s\\
     $2^{24}$ &  -  &-& - & 2h:17m\\\hline
  \end{tabular}
  \label{tab:randomIOtime}
\end{table}

In Table \ref{tab:randomIOtime} we include the results of the I/O experiments in the random model. Each cell in the table corresponds to the execution time (including the waiting time due to disk I/Os when the trees are large and can not fit in RAM). For this experiment we used the \texttt{time} function of Ubuntu and thus also took into account the time taken to parse the input trees. Every cell contains the result of 1 run and for the input trees of size $2^{23}$ and $2^{24}$ we used the 128 bit implementation of the new algorithm in order to avoid overflows. Unlike \texttt{CacheTD}, the performance of \texttt{CPDT} and \texttt{tqDist} deteriorates significantly the moment they start performing disk I/Os.


\begin{table}
  \centering
  \caption{Alpha model: Time performance when limiting the available RAM to 1GB. For both tables we have $\alpha = 0.5$. For the left table we have $p_{1} = p_{2} = 0$ and for the right table $p_{1} = p_{2} = 0.5$.}
  \begin{tabular}{ccccc}
  	\toprule
  	$n$ &\texttt{CPDT/CPDTg}&\texttt{tqDist}& \texttt{CacheTD}\\
  	\midrule
    $2^{15}$ & 0m:01s & 0m:01s & 0m:01s\\
    $2^{16}$ & 0m:01s& 0m:02s & 0m:01s\\
    $2^{17}$ & 0m:01s & 0m:08s & 0m:01s\\
    $2^{18}$ & 0m:02s  & 3m:10s & 0m:01s\\
    $2^{19}$ & 0m:05s  & 2h:16m & 0m:01s\\
    $2^{20}$ & 0m:34s & - & 0m:01s\\
    $2^{21}$ & 11h:20m  & - & 0m:03s\\
    $2^{22}$ & - & - & 0m:35s\\
    $2^{23}$ & -  & - & 10m:09s\\
    $2^{24}$ & - & - & 43m:52s\\\hline
  \end{tabular}
  \begin{tabular}{ccccc}
  	\toprule
  	$n$ &\texttt{CPDT}&\texttt{CPDTg}&\texttt{tqDist}& \texttt{CacheTD}\\
  	\midrule
     $2^{15}$ &  0m:01s & 0m:01s &  0m:01s & 0m:01s\\
     $2^{16}$ & 0m:01s & 0m:01s &  0m:01s & 0m:01s\\
    $2^{17}$ &  0m:01s & 0m:01s &  0m:03s & 0m:01s\\
    $2^{18}$ & 0m:03s & 0m:03s &  1m:18s & 0m:01s\\
     $2^{19}$ & 0m:10s & 0m:07s &  19m:02s & 0m:01s\\
     $2^{20}$ &  1h:58m  &6h:32m& $>$10h & 0m:02s\\
    $2^{21}$ & -  &-& - & 0m:56s\\
     $2^{22}$ &  -  &-& - & 4m:11s\\
      $2^{23}$ & - &-& - & 24m:44s\\
     $2^{24}$ &  -  &-& - & 2h:13m\\\hline
  \end{tabular}
  \label{tab:alphaIOtime}
\end{table}

In Table \ref{tab:alphaIOtime} we have the same experiments but now for the alpha model, where we choose $\alpha = 0.5$. Here we can clearly see how much \texttt{CPDT} and \texttt{CPDTg} are affected by the alpha parameter, e.g. for binary trees, when we go from $n=2^{20}$ tp $n=2^{21}$ the performance in the alpha model degrades by a factor of 1200, while in the random model, the performance degrades only by a factor of 32. We do want to emphasize that if we want to be more precise about these factors, we have to run each experiment several times and not just once. However, the overall outcome of the experiments is that the new algorithm, \texttt{CacheTD}, is clearly scaling to external memory much better than all previous algorithms.

\section{Conclusion}
\label{sec:conclusion}%

In this paper we presented two cache oblivious algorithms for computing the triplet distance between two rooted unordered trees, one that works for binary trees and one that works for arbitrary degree trees. Both require $\mathrm{O}(n \log n)$ time in the RAM model and $\mathrm{O}(\frac{n}{B} \log_{2} \frac{n}{M})$ I/Os in the cache oblivious model. We implemented the algorithms in C++ and showed with experiments that their performance surpasses the performance of previous implementations for this problem. In particular, our algorithms are the first to scale to external memory.

\subparagraph*{Future Work/Open Problems.}

\begin{itemize}

\item In the cache oblivious model could we change our algorithm so that the base of the logarithm becomes $M/B$, i.e. achieve the sorting bound? Would that algorithm be even more efficient practice?

\item Is it possible to compute the triplet distance in $O(n)$ time? 

\item For the quartet distance computation, could we apply similar techniques to those described in this paper in order to get an algorithm with better time bounds in the RAM model that also scales to external memory?

\end{itemize}
%%
%% Bibliography
%%

%% Either use bibtex (recommended), 

\bibliography{lipics-v2016-sample-article}

%% .. or use the thebibliography environment explicitely

\appendix
\section{Experiment Figures}
\label{appendix:experimentFigures}

\begin{figure}[ht]
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/10iterSMALL/time.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/10iterSMALL/time.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/10iterSMALL/time.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/10iterSMALL/time.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=5,width=1\textwidth]{figures/experiments/10iterSMALL/time.pdf}
    \end{subfigure}%
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=6,width=1\textwidth]{figures/experiments/10iterSMALL/time.pdf}
    \end{subfigure}
    \caption{Time performance in the random model.}
    \label{fig:expTimeSmall}
\end{figure}


\begin{figure}
    \centering
        \includegraphics[page=2,width=0.8\textwidth]{figures/experiments/10iterSMALL/gpp.pdf}
    \caption{Comparing the time performance (random model) of \texttt{CPDT} when compiling with g++ version 4.7 and g++ version 5.4.}
    \label{fig:gpp}
\end{figure}


\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/10iterSMALL/space.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/10iterSMALL/space.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/10iterSMALL/space.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/10iterSMALL/space.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=5,width=1\textwidth]{figures/experiments/10iterSMALL/space.pdf}
    \end{subfigure}%
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=6,width=1\textwidth]{figures/experiments/10iterSMALL/space.pdf}
    \end{subfigure}
    \caption{Space performance in the random model.}
    \label{fig:expSpaceSmall}
\end{figure}


\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/10iterSMALL/pTime.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/10iterSMALL/pTime.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/10iterSMALL/pTime.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/10iterSMALL/pTime.pdf}
    \end{subfigure}
    \caption{How the contraction parameter affects execution time (random model).}
    \label{fig:expPtimeSmall}
\end{figure}

\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/10iterSMALL/pSpace.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/10iterSMALL/pSpace.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/10iterSMALL/pSpace.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/10iterSMALL/pSpace.pdf}
    \end{subfigure}
    \caption{How the contraction parameter affects space (random model).}
    \label{fig:expPspaceSmall}
\end{figure}

\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/10iterSMALL/L1misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/10iterSMALL/L1misses.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/10iterSMALL/L1misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/10iterSMALL/L1misses.pdf}
    \end{subfigure}    
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=5,width=1\textwidth]{figures/experiments/10iterSMALL/L1misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=6,width=1\textwidth]{figures/experiments/10iterSMALL/L1misses.pdf}
    \end{subfigure}
    \caption{L1 misses in the random model.}
    \label{fig:l1misses}
\end{figure}

\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/10iterSMALL/L2misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/10iterSMALL/L2misses.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/10iterSMALL/L2misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/10iterSMALL/L2misses.pdf}
    \end{subfigure}
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=5,width=1\textwidth]{figures/experiments/10iterSMALL/L2misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=6,width=1\textwidth]{figures/experiments/10iterSMALL/L2misses.pdf}
    \end{subfigure}
    \caption{L2 misses in the random model.}
    \label{fig:l2misses}
\end{figure}

\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/10iterSMALL/L3misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/10iterSMALL/L3misses.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/10iterSMALL/L3misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/10iterSMALL/L3misses.pdf}
    \end{subfigure}
     \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=5,width=1\textwidth]{figures/experiments/10iterSMALL/L3misses.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=6,width=1\textwidth]{figures/experiments/10iterSMALL/L3misses.pdf}
    \end{subfigure}
    \caption{L3 misses in the random model.}
    \label{fig:l3misses}
\end{figure}



\begin{figure}
\captionsetup[subfigure]{justification=centering}
    \centering
    \begin{subfigure}{0.5\textwidth}
        \centering
        \includegraphics[page=1,width=1\textwidth]{figures/experiments/alpha/1.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=2,width=1\textwidth]{figures/experiments/alpha/1.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=3,width=1\textwidth]{figures/experiments/alpha/1.pdf}
    \end{subfigure}%
    \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=4,width=1\textwidth]{figures/experiments/alpha/1.pdf}
    \end{subfigure}%
    \\
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=5,width=1\textwidth]{figures/experiments/alpha/1.pdf}
    \end{subfigure}%
     \begin{subfigure}{0.5\textwidth}
        \centering
         \includegraphics[page=6,width=1\textwidth]{figures/experiments/alpha/1.pdf}
    \end{subfigure}
    \caption{Alpha parameter and running time.}
    \label{fig:alphaFigure}
\end{figure}

\end{document}

lipics-v2016-sample-article.bbl
\begin{thebibliography}{10}

\bibitem{IOmodel}
A.~Aggarwal and J.S. Vitter.
\newblock The input/output complexity of sorting and related problems.
\newblock {\em Communications of the ACM}, 31(9):1116--1127, 1988.

\bibitem{Bansal2008}
M.S. Bansal, J.~Dong, and D.~Fernández-Baca.
\newblock Comparing and aggregating partially resolved trees.
\newblock {\em Theoretical Computer Science}, 412(48):6634 -- 6652, 2011.

\bibitem{Qmethod}
V.~Berry and O.~Gascuel.
\newblock Inferring evolutionary trees with strong combinatorial evidence.
\newblock {\em Theoretical Computer Science}, 240(2):271--298, 2000.

\bibitem{TripletQuartetSoda13}
G.S. Brodal, R.~Fagerberg, C.N.S. Pedersen, T.~Mailund, and A.~Sand.
\newblock Efficient algorithms for computing the triplet and quartet distance
  between trees of arbitrary degree.
\newblock {\em 24th Annual ACM-SIAM Symposium on Discrete Algorithms}, pages
  1814--1832, 2013.

\bibitem{Critchlow}
D.E. Critchlow, D.K. Pearl, and C.L. Qian.
\newblock {T}he {T}riples {D}istance for {R}ooted {B}ifurcating {P}hylogenetic
  {T}rees.
\newblock {\em Systematic Biology}, 45(3):323, 1996.

\bibitem{DaysAlgorithm}
William~H.E. Day.
\newblock Optimal algorithms for comparing trees with labeled leaves.
\newblock {\em Journal of Classification}, 2(1):7--28, 1985.

\bibitem{Dobson}
A.J. Dobson.
\newblock Comparing the shapes of trees.
\newblock {\em Combinatorial Mathematics III. Lecture Notes in Mathematics},
  pages 95--100, 1975.

\bibitem{Quartet}
G.F. Estabrook, F.R. McMorris, and C.A. Meacham.
\newblock Comparison of undirected phylogenetic trees based on subtrees of four
  evolutionary units.
\newblock {\em Systematic Zoology}, 34(2):193--200, 1985.

\bibitem{COmodel}
M.~Frigo, C.E. Leiserson, H.~Prokop, and S.~Ramachandran.
\newblock Cache-oblivious algorithms.
\newblock {\em 40th Annual IEEE Symposium on Foundations of Computer Science},
  pages 285--297, 1999.

\bibitem{alenex14}
M.K. Holt, J.~Johansen, and G.S. Brodal.
\newblock On the scalability of computing triplet and quartet distances.
\newblock {\em 16th Workshop on Algorithm Engineering and Experiments}, pages
  9--19, 2014.

\bibitem{Jansson2015}
J.~Jansson and R.~Rajaby.
\newblock A more practical algorithm for the rooted triplet distance.
\newblock {\em International Conference on Algorithms for Computational
  Biology}, pages 109--125, 2015.

\bibitem{pmid27983874}
J.~Jansson and R.~Rajaby.
\newblock {{A} {M}ore {P}ractical {A}lgorithm for the {R}ooted {T}riplet
  {D}istance}.
\newblock {\em Journal of Computational Biology}, 24(2):106--126, 2017.

\bibitem{RobinsonFoulds}
D.F. Robinson and L.R. Foulds.
\newblock Comparison of phylogenetic trees.
\newblock {\em Mathematical Biosciences}, 53(1):131--147, 1981.

\bibitem{NJmethod}
N.~Saitou and M.~Nei.
\newblock The neighbor-joining method: a new method for reconstructing
  phylogenetic trees.
\newblock {\em Molecular Biology and Evolution}, 4(4):406, 1987.

\bibitem{bmc13}
A.~Sand, G.S. Brodal, R.~Fagerberg, C.N.S. Pedersen, and T.~Mailund.
\newblock A practical {O}($n \log^2 n$) time algorithm for computing the
  triplet distance on binary trees.
\newblock {\em BMC Bioinformatics}, 14(2):S18, 2013.

\bibitem{bioinformatics14}
A.~Sand, M.K. Holt, J.~Johansen, G.S. Brodal, T.~Mailund, and C.N.S. Pedersen.
\newblock tqdist: A library for computing the quartet and triplet distances
  between binary or general trees.
\newblock {\em Bioinformatics}, 30(14):2079, 2014.

\bibitem{biology13}
A.~Sand, M.K. Holt, J.~Johansen, R.~Fagerberg, G.S. Brodal, C.N.S. Pedersen,
  and T.~Mailund.
\newblock Algorithms for computing the triplet and quartet distances for binary
  and general trees.
\newblock {\em Biology - Special Issue on Developments in Bioinformatic
  Algorithms}, 2(4):1189--1209, 2013.

\end{thebibliography}
