\documentclass[nojss]{jss}
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{Efficient computing with BART}
\usepackage{verbatim}
\usepackage{statex2}
\usepackage[authoryear,round]{natbib}
\usepackage{rotating}

%% -- Article metainformation (author, title, ...) -----------------------------

%% - \author{} with primary affiliation
%% - \Plainauthor{} without affiliations
%% - Separate authors by \And or \AND (in \author) or by comma (in \Plainauthor).
%% - \AND starts a new line, \And does not.
%\author{Achim Zeileis\\Universit\"at Innsbruck
%   \And Second Author\\Plus Affiliation}
%\Plainauthor{Achim Zeileis, Second Author}
\author{Rodney Sparapani\\Medical College of Wisconsin
\And Robert McCulloch\\Arizona State University}
\Plainauthor{Rodney Sparapani, Robert McCulloch}

%% - \title{} in title case
%% - \Plaintitle{} without LaTeX markup (if any)
%% - \Shorttitle{} with LaTeX markup (if any), used as running title
%\title{A Short Demo Article: Regression Models for Count Data in \proglang{R}}
%\Plaintitle{A Short Demo Article: Regression Models for Count Data in R}
%\Shorttitle{A Short Demo Article in \proglang{R}}

\title{Efficient computing with BART}
\Plaintitle{Efficient computing with BART}
\Shorttitle{Efficient computing with BART}

%% - \Abstract{} almost as usual
\Abstract{ This short article illustrates how to perform efficient
  computing with the {\bf BART} package.  }

%% - \Keywords{} with LaTeX markup, at least one required
%% - \Plainkeywords{} without LaTeX markup (if necessary)
%% - Should be comma-separated and in sentence case.
\Keywords{Bayesian Additive Regression Trees, multi-threading,
  \proglang{R}, \proglang{C++}, forking, OpenMP} 
\Plainkeywords{Bayesian Additive Regression Trees, 
  multi-threading, R, C++, forking, OpenMP}

%% - \Address{} of at least one author
%% - May contain multiple affiliations for each author
%%   (in extra lines, separated by \emph{and}\\).
%% - May contain multiple authors for the same affiliation
%%   (in the same first line, separated by comma).

\Address{
  Rodney Sparapani {rsparapa@mcw.edu}\\
  Division of Biostatistics, Institute for Health and Equity\\
  Medical College of Wisconsin, Milwaukee campus
%  8701 Watertown Plank Road\\
%  Milwaukee, WI\ \ 53226, USA\\
%  E-mail: {rsparapa@mcw.edu}
%URL: \url{https://www.mcw.edu/Biostatistics/Faculty-Staff/Rodney-Sparapani-PhD.htm}
}

\begin{document}

\maketitle

\section{Efficient computing with BART}

If you had the task of creating an efficient implementation for a
black-box model such as BART, which tools would you use?
Surprisingly, linear algebra routines which are a traditional building
block of scientific computing will be of little use for a tree-based
method such as BART.  So what is needed?  Restricting ourselves to
widely available off-the-shelf hardware and open-source sofware,
we believe there are four key technologies necessary for a successful
BART implementation.
\begin{itemize}
\item an object-oriented language to facilitate working with trees
\item a parallel (or distributed) CPU computing framework for faster processing
\item a high-quality parallel random number generator
%\item an interpreted shell for high-level data handling
\end{itemize}
In our implementation of BART, we pair the objected-oriented languages
of \proglang{R} and \proglang{C++} to satisfy these requirements.  In
this Section, we give a brief introduction to the concepts and
technologies harnessed for efficient computing by our {\bf BART}
package.

\subsection{A brief history of multi-threading}

We now present a short synopsis of multi-threading.  This is not meant
to be exhaustive; rather, we only provide enough detail to explain the
capability and popularity of multi-threading today.  Multi-threading
emerged rather early in the digital computer age with pioneers laying the
research groundwork in the 1950s and 60s.  In 1961, Burroughs released
the B5000 which was the first commercial hardware capable of
multi-threading \citep{Lync65}.  The B5000 performed asymmetric
multiprocessing which is commonly employed in modern hardware via
numerical co-processors and/or graphical processors today.  In 1962,
Burroughs released the D825 which was the first commercial hardware
capable of symmetric multiprocessing (SMP) with CPUs
\citep{AndeHoff62}.  In 1967, Gene Amdahl derived the theoretical
limits for multi-threading which came to be known as Amdahl's law
\citep{Amda67}.  If $C$ is the number of CPUs and $b$ is the fraction
of work that can't be parallelized, then the gain due to
multi-threading is $((1-b)/C+b)^{-1}$.

Let's fast-forward to the modern era of multi-threading.  Although,
not directly related to multi-threading, in 2000, Advanced Micro
Devices (AMD) released the AMD64 specification that created a new
64-bit x86 instruction set which was capable of co-existing with
16-bit and 32-bit x86 legacy instructions.  This was an important
advance since 64-bit math is capable of addressing vastly more memory
than 16-bit or 32-bit ($2^{64}$ vs.\ $2^{16}$ or $2^{32}$) and
multi-threading inherently requires more memory resources.  In 2003,
version 2.6 of the Linux kernel incorporated full SMP support; prior
Linux kernels had either no support or very limited/crippled support.
From 2005 to 2011, AMD released a series of Opteron chips with
multiple cores for multi-threading: 2 cores in 2005, 4 cores in 2007,
6 cores in 2009, 12 cores in 2010 and 16 cores in 2011.  From 2008 to
2010, Intel entered the market with Xeon chips and their
hyperthreading technology that allows each core to issue two
instructions per clock cycle: 4 cores (8 threads) in 2008 and 8 cores
(16 threads) in 2010.  In this era including today, most off-the-shelf
hardware available features 1 to 4 CPUs.  Therefore, in the span of
only a few years, multi-threading rapidly trickled down from servers
at large firms to mass-market products such as desktops/laptops.  For
example, the consumer machine that {\bf BART} is developed on,
purchased in 2016, is capable of 8 threads (and hence many of the
examples default to 8 threads).

\subsection{Modern multi-threading software frameworks}

In the late 1990s, the Message Passing Interface (MPI) \citep{WalkDong96}
was introduced which is the dominant distributed computing framework
in use today \citep{GabrFagg04}, i.e., distributed meaning tasks which
span multiple computers called nodes.  \proglang{R} has some support
for MPI provided in the {\bf parallel} \citep{UrbaRipl17} package and
additional support is provided by several other CRAN packages such as
{\bf snow} \citep{TierRoss16}, {\bf Rmpi} \citep{Yu17} and {\bf
  pbdMPI} \citep{ChenOstr18}.  To support MPI, BART sofware was
re-written with a simple, readable C++ object schema.  Although, not
an \proglang{R} package, this project was documented by
\cite{PratChip14}.  The current {\bf BART} package source code is a
descendent of the MPI BART project which is deprecated.

Furthermore, the current {\bf BART} package no longer supports MPI;
rather, the multi-threading available is now based on the OpenMP
standard \citep{DaguMeno98} and the {\bf parallel} package.  OpenMP
takes advantage of modern hardware by performing multi-threading on
single machines which often have multiple CPUs each with multiple
cores.  Currently, the {\bf BART} package only uses OpenMP for
parallelizing \code{predict} function calculations.  The challenge
with OpenMP, besides the programming per se, is that it is not widely
available on all platforms.  Operating system support can be detected
by the GNU autotools \citep{Calc10} which define a C pre-processor
macro if it is available, \code{_OPENMP}, or not.  There are numerous
exceptions for operating systems so it is difficult to generalize.
But, generally, Microsoft Windows lacks OpenMP detection since the GNU
autotools do not natively exist on this platform.  And, Apple macOS
lacks OpenMP support since the standard Xcode toolkit does not provide
it.  Thankfully, most Linux and UNIX distributions do provide OpenMP
(although, macOS is technically a UNIX distribution, yet it is a
notable exception in this regard).  We provide the function
\code{mc.cores.openmp} which returns $>0$ (0) if the \code{predict}
function is (not) capable of utilizing OpenMP.

The {\bf parallel} package provides multi-threading via forking.
Forking is available on Unix platforms, but not Windows (we use the
term Unix to refer to UNIX, Linux and macOS since they are all in
the UNIX family tree).  The {\bf BART} package uses forking for
posterior sampling of the $f$ function, and also for the
\code{predict} function when OpenMP is not available.  Except for
\code{predict}, all functions that use forking start with \code{mc}.
And, regardless of whether OpenMP or forking is employed, these
functions except the argument \code{mc.cores} which controls the
number of threads to be used.  The {\bf parallel} package provides the
function \code{detectCores} which returns the number of threads that
your hardware can support.

\subsection{BART implementations on CRAN}

Currently, there are four BART implementations on the Comprehensive R
Archive Network (CRAN); see the Appendix for a tabulated comparative
summary of their features.

{\bf BayesTree} was the first released in 2006 \citep{ChipMcCu16}.
Reported bugs will be fixed, but no future improvements are planned;
so, we suggest choosing one of the newer packages such as {\bf BART}.
The basic interface and workflow of {\bf BayesTree} has strongly
influenced the other packages which followed.  However, the {\bf
  BayesTree} source code is difficult to maintain and, therefore,
improvements were limited leaving it with relatively fewer
features than the other entries.

The next entrant is {\bf bartMachine} which is written in
\proglang{java} and was first released in 2013 \citep{KapeBlei16}.  It
provides advanced features like multi-threading, variable selection
\citep{BleiKape14}, a \code{predict} function, convergence diagnostics
and missing data handling.  However, the \proglang{R} to
\proglang{java} interface can be challenging to deal with.
\proglang{R} is written in \proglang{C} and \proglang{Fortran},
consequentally, functions written in \proglang{java} do not have a
natural interface to \proglang{R}.  This interface is provided by the
{\bf rJava} \citep{Urba17} package which requires the Java Development
Kit (JDK).  Therefore, we recommend {\bf bartMachine} only for those
users who have a firm grounding in the \proglang{java} language and
its tools in order to install/upgrade the package and get the best
performance out of it.

The next entrant is {\bf dbarts} which is written in \proglang{C++}
and was first released in 2014 \citep{DoriChip16}.  It is a clone of
the {\bf BayesTree} interface, but it does not share the source code;
{\bf dbarts} source has been re-written from scratch for efficiency
and maintainability.  {\bf dbarts} is a drop-in replacement for 
{\bf BayesTree}.  However, {\bf dbarts} has
relatively fewer features than the other entries.

The {\bf BART} package which is written in \proglang{C++} was first
released in 2017 \citep{McCuSpar18}.  It provides advanced features
like multi-threading, variable selection \citep{Line16}, a
\code{predict} function and convergence diagnostics.  The source code
is a descendent of the MPI BART project.  Although, \proglang{R} is
mainly written in \proglang{C} and \proglang{Fortran} (at the time of
this writing, 39.2\% and 26.8\% lines of source code respectively),
\proglang{C++} is a natural choice for creating \proglang{R} functions
since they are both object-oriented languages.  The \proglang{C++}
interface to \proglang{R} has been seamlessly provided by the {\bf
  Rcpp} package \citep{EddeFran11} which efficiently passes object
references from \proglang{R} to \proglang{C++} (and vice versa) as
well as providing direct accesss to the \proglang{R} random number
generator.  The source code can also be called from \proglang{C++}
alone without an \proglang{R} instance where the random number
generation is provided by either the standalone Rmath library
\citep{Rcore17} or the \proglang{C++} \code{random} Standard Template
Library.  Also, it is the only BART package to support categorical and
time-to-event outcomes \citep{SparLoga16}.  It does not provide
missing data imputation; rather, we recommend the {\bf sbart} package
for this niche which is performed by the so-called Sequential BART
algorithm \citep{DaniSing17,XuDani16} ({\bf sbart} is also a
descendent of MPI BART).

\subsection{MCMC is embarrassingly parallel}

In general, Bayesian Markov chain Monte Carlo (MCMC) posterior
sampling is considered to be embarrassingly parallel
\citep{RossTier07}, i.e., since the chains only share the data and
don't have to communicate with each other, parallel implementations
are considered to be trivial.  BART MCMC also falls into this class.
Typical practice for Bayesian MCMC is to start in some initial state,
perform a limited number of samples to generate a new random starting
position and throw away the preceding samples which we call burn-in
(the amount of burn-in in the {\bf BART} package is controlled by the
argument \code{nskip} which defaults to either 100 or 250).  The total
length of the chain returned is controlled by the argument
\code{ndpost} which defaults to 1000.  The theoretical gain due to
multi-threading can be calculated by what we call the MCMC Corollary
to Amdahl's Law.  Let $b$ be the burn-in fraction and $C$ be the
number of threads, then the gain limit is $((1-b)/C+b)^{-1}$.  (As an
aside, note that we can derive Amdahl's Law as follows where the
amount of work done is in the numerator and elapsed time is in the
denominator: $\frac{1-b+b}{(1-b)/C+b}=\frac{1}{(1-b)/C+b}$).  For
example, see the diagram in Figure~\ref{MCMC} where the burn-in
fraction, $b=\frac{100}{1100}=0.09$, and the number of CPUs, $C=5$,
results in an elapsed time of only $((1-b)/C+b)=0.27$ or a
$((1-b)/C+b)^{-1}=3.67$ fold reduction which is the gain in
efficiency.  In Figure~\ref{Amdahl}, we plot theoretical gains
on the y-axis and the number of CPUs on the x-axis
for two settings: $b \in \{0.025, 0.1\}$.

\subsection{Multi-threading and random access memory (RAM)}

The IEEE standard 754-2008 \citep{IEEE08} specifies that every
double-precision number consumes 8 bytes (64 bits).  Therefore, it is
quite simple to estimate the amount of random access memory (RAM)
required to store a matrix.  If $A$ is $m \times n$, then the amount
of RAM needed is $8 \times m \times n$ bytes.  Large matrices held in
RAM can present a challenge to system performance.  If you consume all
of the physical RAM, the system will ``swap'' segments out to virtual
RAM which are disk files and this can degrade performance and possibly
even crash the system.  On Unix, you can monitor memory and swap usage
with the \code{top} command-line utility.  And, within \proglang{R},
you can determine the size of an object with the \code{object.size}
function.

Mathematically, a matrix is represented as follows.
\begin{align*}
A=\wrap{\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots&\vdots & \vdots & \vdots \\
a_{m1} & a_{m2} & \cdots & a_{mn} \\
\end{array}}
\end{align*}
\proglang{R} is a column-major language, i.e., matrices are laid out
in consecutive memory locations by traversing the columns:
$\wrap{a_{11}, a_{21}, \., a_{12}, a_{22}, \.}$.  \proglang{R} is
written in \proglang{C} and \proglang{Fortran} where
\proglang{Fortran} is a column-major language as well.  However,
\proglang{C} and \proglang{C++} are row-major languages, i.e.,
matrices are laid out in consecutive memory locations by traversing
the rows: $\wrap{a_{11}, a_{12}, \., a_{21}, a_{22}, \.}$.  So, if you
have written an \proglang{R} function in \proglang{C/C++}, then you
need to be cognizant of the clash in paradigms (also note that
\proglang{R/Fortran} array indexing goes from 1 to $m$ while
\proglang{C/C++} indexing goes from 0 to $m-1$).  As you might
surmise, this is easily addressed with a transpose, i.e., instead of
passing $A$ from \proglang{R} to \proglang{C/C++}, pass $A^t$.

\proglang{R} is very efficient in passing objects; rather, than
passing an object (along with all of its memory consumption) on the
stack, it passes objects merely by a pointer referencing the original
memory location.  However, \proglang{R} follows copy-on-write memory
allocation, i.e., all objects present in the parent thread can be read
by a child thread without a copy, but when an object is
altered/written by the child, then a new copy is created in memory.
Therefore, if we pass $A$ from \proglang{R} to \proglang{C/C++}, and then
transpose, we will create multiple copies of $A$ consuming
$8 \times m \times n \times C$ where $C$ is the number of children.
If $A$ is a large matrix, then you may stress the system's limits.
The simple solution is for the parent to create the transpose before
passing $A$ and avoiding the multiple copies, i.e., \code{A <- t(A)}.
And this is the philosophy that the {\bf BART} package follows.

\subsection{Multi-threading: interactive and batch processing}

Interactive jobs must take precedence over batch jobs to prevent the
user experience from suffering high latency.  For example, have you
ever experienced a system slowdown while you are typing and the
display of your keystrokes can not keep up; this should never happen
and is the sign of something amiss.  With large multi-threaded jobs,
it is surprisingly easy to naively degrade system performance.  But,
this can easily be avoided by operating system support provided by
\proglang{R}.  In the {\bf tools} package \citep{HornLeis17}, there is
the \code{psnice} function.  Paraphrased from the \code{?psnice} help
page.
\begin{quote}
  Unix has a concept of process priority.  Priority is assigned values from 0
  to 39 with 20 being the normal priority and (counter-intuitively)
  larger numeric values denoting lower priority.  Adding to the
  complexity, there is a ``nice'' value, the amount by which the
  priority exceeds 20.  Processes with higher nice values
  will receive less CPU time than those with normal priority.
  Generally, processes with nice 19 are only run when the system would
  otherwise be idle.
\end{quote}
Therefore, by default, the {\bf BART} package children have their nice
value set to 19.

\subsection{Continuous BART: serial and parallel implementations}

Here we present snippets of \proglang{R} code to run BART in serial
and parallel for a continuous outcome which we call continuous BART.
While we only demonstrate continuous outcomes, the other outcomes
are as similarly handled by the {\bf BART} package as possible to
present a consistent interface.
The serial function is \code{wbart} and the parallel, \code{mc.wbart}.
The 'w' in the name stands for weighted since you can provide known
weights (with the \code{w} argument) for the following model:
$y_i ~ \N{f(\bm{x}_i}{w_i^2\sd^2} \where (f, \sd) \prior \mathrm{BART}$
and $i=1, \., N$ indexes subjects.  Now,
we can perform the calculations in serial,\\ 
\code{set.seed(99); post <-
  wbart(x.train, y.train, ..., ndpost=M)}\\ % keepevery=1)}\\
 or in parallel  (when said support is available),\\
\code{post <- {mc.wbart}(x.train, y.train, ..., ndpost=M, %keepevery=1, 
mc.cores=8, seed=99)}.\\
  Notice the difference in how the seed is set;
we will return to this detail later on.  The {\bf BART} package
allows \code{x.train} (and \code{x.test}) to be provided as matrices
or data frames, but for simplicity we present them as matrices.
\begin{align*}
\mbox{Input: \code{x.train} and, optionally, \code{x.test}:\ } &
\wrap{\begin{array}{c}
\bm{x}_{1} \\
\bm{x}_{2} \\
\vdots \\
\bm{x}_{N} \\
\end{array}} \where \bm{x}_{i} \mbox{\ is the $i^{th}$ row} \\
\mbox{Output: \code{post\$yhat.train} and \code{post\$yhat.test}:\ } &
\wrap{\begin{array}{ccc}
\hat{y}_{11}& \. & \hat{y}_{N1} \\
\vdots & \vdots & \vdots \\
\hat{y}_{1M}& \. & \hat{y}_{NM} \\
\end{array}} \where \hat{y}_{im}=f_m(\bm{x}_i) \\
\end{align*}
The \code{post} object returned is of type \code{wbart} which is
essentially a list.  There are other items returned in the list, but
here we only focus on \code{post\$yhat.train} and
\code{post\$yhat.test}; the latter only being returned if
\code{x.test} is provided.  In the above display, $m=1, \., M$ are the
MCMC samples which are the rows of \code{post\$yhat.train} and
\code{post\$yhat.test}.
Note that each outcome has a different return type, i.e., \code{post}
object of type \code{wbart} (continuous), \code{pbart} (binary
probit), \code{lbart} (binary logistic), \code{survbart} (survival
analysis), \code{criskbart} (competing risks) or \code{recurbart}
(recurrent events).

\subsection{Continuous BART: predicting with a previous fit}

Often when we are fitting a BART model, we have not specified an
\code{x.test} matrix of hypothetical values for the evaluation of $f$.
For fire-and-forget packages like {\bf BayesTree} and {\bf dbarts}, we
would have to re-fit the model every time we want to evaluate {\bf
  x.test} which can be very time-consuming.  Therefore, the {\bf BART}
package takes a unique approach: it returns the ensemble of trees in
the \code{post} object for later use; specifically, they are encoded
in an ASCII character string, \code{post\$treedraws\$trees}.  This allows us to
construct \code{x.test} after the fact which is often convenient when
it is large since we can partition it into smaller chunks.  Then we
can evaluate predictions via the S3 method \code{predict.wbart}.  The
predictions are generated in serial by default,\\
\code{pred <- predict({post}, {x.test}, ...)} \\
but can be parallelized (when said support is available),\\
\code{pred <- predict({post}, {x.test}, {mc.cores=8}, ...)}.\\

\begin{align*}
\mbox{Input: \code{x.test}:\ } &
\wrap{\begin{array}{c}
x_{1} \\
x_{2} \\
\vdots \\
x_{Q} \\
\end{array}} \where \bm{x}_{i} \mbox{\ is the $i^{th}$ row and\ } i=1, \., Q \\
\mbox{Output matrix:\ } &
\wrap{\begin{array}{ccc}
\hat{y}_{11}& \. & \hat{y}_{Q1} \\
\vdots & \vdots & \vdots \\
\hat{y}_{1M}& \. & \hat{y}_{QM} \\
\end{array}} \where \hat{y}_{im}=f_m(\bm{x}_i) \\
\end{align*}
In the above display, $m=1, \., M$ are the
MCMC samples which are the rows of the output matrix.

\subsection{Binary Trees}

BART relies on an ensemble of $H$ binary trees.  We exploit the tree
metaphor to its fullest.  Each of these trees grows from the ground up
starting out as a root node.  The root node is generally a branch
decision rule, but it doesn't have to be; occasionally there is only a
root terminal node and its corresponding leaf output value.  If the
root is a branch decision rule, then it spawns a left and a right node
which each can be either a branch decision rule or a terminal leaf
value and so on.  
In tree, ${T}$, there are $C$ nodes which are made of $B$
branches and $L$ leaves: $C=B+L$.  There is an algebraic relationship
between the number of branches and leaves which we express as $C= 2L-1$.

The ensemble of trees is encoded in an ASCII string
which is returned in the \code{treedraws\$trees} list item.  This
string can be easily exported and imported by \proglang{R} with the
following:\\
\code{write(post\$treedraws\$trees, 'trees.txt')\\
  tc <- textConnection(post\$treedraws\$tree)\\
  trees <- read.table(file=tc, fill=TRUE, row.names=NULL, header=FALSE,\\
  col.names=c('node', 'var', 'cut', 'leaf'))\\
  close(tc)}\\

The string is encoded via the following binary tree notation.  The
first line is an exception which has the number of MCMC samples, $M$,
in the field \code{node}; the number of trees, $H$, in the field
\code{var}; and the number of variables, $P$, in the field \code{cut}.
For the rest of the file, the field \code{node} is used for the number
of nodes in the tree when all other fields are \code{NA}; or for a
specific node when the other fields are present.  The nodes are
numbered in relation to the tier level, \code{t=floor(log2(node))}, as
follows.
\begin{table}[!h]
\begin{tabular}{r|ccccccc}
Tier & \\
$t$ & \multicolumn{3}{c}{$2^t+0$} & $\dots$ &
\multicolumn{3}{c}{$2^t+2^t-1$} \\ \hline
$\vdots$ & \\
2 & 4 &   & 5 &   & 6 &   & 7 \\
1 &   & 2 &   &   &   & 3 &   \\
0 &   &   &   & 1 &   &   &   \\
\end{tabular}
\end{table}

The \code{var} field is the variable in the branch decision rule which
is encoded $0, \dots, P-1$ as a \proglang{C/C++} array index (rather
than an R index).  Similarly, the \code{cut} field is the cutpoint of
the variable in the branch decision rule which is encoded
$0, \dots, c_j-1$ for variable $j$; note that the cutpoints are
returned in the \code{treedraws\$cutpoints} list item.  The terminal
leaf output value is contained in the field \code{leaf}.  It is not
immediately obvious which nodes are branches vs.\ leaves since it
appears that the \code{leaf} field is given for both branches and
leaves.  Leaves are always associated with \code{var=0} and
\code{cut=0}; however, note that this is also a valid branch 
variable/cutpoint since these are \proglang{C/C++} indices.
% The key insight is that the 
% first $B$ rows of each node are branches and the rest are leaves.
The key to discriminating between branches and leaves is the algebraic
relationship between a branch, $v$, at tier $t$ leading to its left,
$l$, and right, $r$, nodes at tier $t+1$.  If $v=2^t+k$, then
$l=2^{t+1}+2k$ and $r=2^{t+1}+2k+1$, i.e., for each node, besides root,
you can determine from which branch it arose and those nodes that are
not a branch (since they have no leaves) are necessarily leaves.
  
\subsection{Creating a BART executable}

Occasionally, you may need to create a BART executable that you can
run without an \proglang{R} instance.  This is especially useful if
you need to include BART in another \proglang{C++} program.  Or, when
you need to debug the {\bf BART} package \proglang{C++} source code
which is more difficult to do when you are calling the function from
\proglang{R}.  Several examples of these are provided with the {\bf
  BART} package.  With \proglang{R}, you can find the \code{Makefile}
and the weighted BART example with
\code{system.file('cxx-ex/Makefile', package='BART')} and\\
\code{system.file('cxx-ex/wmain.cpp', package='BART')} respectively.
Note that these examples require the installation of the standalone
Rmath library \citep{Rcore17} which is contained in the \proglang{R}
source code distribution.  Rmath provides common \proglang{R}
functions and random number generation, e.g., \code{pnorm} and
\code{rnorm}.  You will likely need to copy the \code{cxx-ex}
directory to your workspace.  Once done, you can build and run the
weighted BART executable example from the command line as
follows.\\
\code{> make wmain.out \#\# to build\\> ./wmain.out \#\# to run}\\
By default, these examples are based on the Rmath random number
generator.  However, you can specify the \proglang{C++} Standard
Template Library random number generator (contained in the
STL \code{random} header file) by uncommenting the following line in
the \code{Makefile} (by removing the pound, \#, symbols):\\
\code{\#\# CPPFLAGS = -I. -I/usr/local/include -DMATHLIB_STANDALONE -DRNG_random}\\
 (which still requires Rmath for other purposes).
These examples were developed on Linux and macOS, but they should be
readily adaptable to UNIX and Windows as well.

\bibliography{references}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.48]{figures/parallel.pdf}
\end{center}
\caption{The theoretical gain due to multi-threading can be calculated
  by Amdahl's Law.  Let $b$ be the burn-in fraction and $C$ be the
  number of threads, then the gain limit is $((1-b)/C+b)^{-1}$.  In
  this diagram, the burn-in fraction, $b=\frac{100}{1100}=0.09$, and
  the number of CPUs, $C=5$, results in an elapsed time of only
  $((1-b)/C+b)=0.27$ or a $((1-b)/C+b)^{-1}=3.67$ fold reduction which
  is the gain in efficiency. \label{MCMC} }
\end{figure}

\begin{figure}[t!]
\begin{center}
\includegraphics[scale=0.42]{figures/amdahl.pdf}
\end{center}
\caption{The theoretical gain due to multi-threading can be calculated
  by Amdahl's Law.  Let $b$ be the burn-in fraction and $C$ be the
  number of threads, then the gain limit is $((1-b)/C+b)^{-1}$.  In
  this figure, the theoretical gains are on the y-axis and the number
  of CPUs, the x-axis, for two settings: $b \in \{0.025,
  0.1\}$. \label{Amdahl} }
\end{figure}

\appendix

\begin{sidewaystable}
\begin{center}%\label{BARTcompare}
\begin{tabular}{l|llll}
Category                 & {\bf BayesTree}      & {\bf bartMachine}  & {\bf dbarts}     & {\bf BART}           \\ \hline
First release            & 2006                 & 2013               & 2014             & 2017                 \\
Authors                  & Chipman              & Kapelner           & Dorie,           & McCulloch, Sparapani \\
                         & \& McCulloch         & \& Bleich          & Chipman          & Gramacy, Spanbauer   \\
                         &                      &                    & \& McCulloch     & \& Pratola           \\
Source code              & C++                  & java               & C++              & C++                  \\
CRAN dependencies        & nnet                 & rJava, car,        &                  & Rcpp                 \\
                         &                      & randomForest,      &                  &                      \\
                         &                      & missForest         &                  &                      \\
Tree transition proposals& 4                    & 3                  & 4                & 3                    \\
Multi-threaded           & No                   & Yes                & No               & Yes                  \\
\code{predict} function  & No                   & Yes                & No               & Yes                  \\
Variable selection       & No                   & Yes                & No               & Yes                  \\
Continuous outcomes      & Yes                  & Yes                & Yes              & Yes                  \\
Dichotomous outcomes with Normal latents      
                         & Yes                  & Yes                & Yes              & Yes                  \\
Dichotomous outcomes with Logistic latents    
                         & No                   & No                 & No               & Yes                  \\
Categorical outcomes     & No                   & No                 & No               & Yes                  \\
Time-to-event outcomes   & No                   & No                 & No               & Yes                  \\
Convergence diagnostics  & No                   & Yes                & No               & Yes                  \\
Thinning                 & Yes                  & No                 & Yes              & Yes                  \\
Cross-validation         & No                   & Yes                & Yes              & No                   \\
Missing data handling    & No                   & Yes                & No               & No                   \\
Partial dependence plots & Yes                  & Yes                & Yes              & No                   \\ 
Citations                & \multicolumn{4}{l}{\cite{ChipMcCu16}}                                               \\ 
                         &                      & \multicolumn{3}{l}{\cite{KapeBlei16}}                        \\ 
                         &                      &                    & \multicolumn{2}{l}{\cite{DoriChip16}}   \\ 
                         &                      &                    &                  & \cite{McCuSpar18}    \\ \hline
\end{tabular}
\end{center}
\end{sidewaystable}

\end{document}

%Local Variables:
%TeX-command-default: "w"
%End:
