text:
\section{Introduction}

Probabilistic graphical models play a central role in modern statistical data analysis.
Their compact and elegant way to visualise and represent complex dependence structures in multivariate probability distributions have shown to be successfully applicable in many scientific domains, ranging from disciplines such as social sciences and image analysis to biology, medical diagnosis and epidemiology \citep[see \emph{e.g.}][]{elwert2013graphical, friedman2000using, friedman2004inferring,10.1093/schbul/sbx013, kuipers2018mutational,  kuipers_moffa_kuipers_freeman_bebbington_2019}.
\felix{Intro to and application of GMs.}

One of the main advantages of graphical models is that they provide a tool for experts and researchers from non-statistical fields to easily specify their assumptions in a specific problem domain, in a way that is very transparent.
However, in many realistic situations the number of variables is either too large to build networks by hand, as it is for example the case in genomic applications, or simply no prior knowledge is available about the relationship between different variables.
As a consequence, there has been a growing interest in automated strategies to infer the graph component of a probabilistic graphical model from data, so-called \emph{structure learning}.
For recent reviews of the wide flora of algorithms that have emerged in the last two decades, the reader is referred to \cite{koski2012review, constantinou2021large, scutari2019learns}.

\felix{Intro to structure learning.}
Since the scientific principle of reproducibility is being recognised as instrumental to achieve high quality standards in research, there is an increasing demand to present new algorithms with publicly available source code and datasets \citep[see \emph{e.g.}][]{Munafo:2017tp, Lamprecht:2020uj}.
Ideally with publicly available software we should be able to easily compare different methods with respect to their performance in different settings. The reality is that the practical implementation of new approaches display high levels of heterogeneity on many aspects, which renders comparative studies both challenging and time consuming.
One of the main difficulties is that different algorithms may be implemented in different programming languages, with different dependencies.
Different packages may differ with respect to the representation of graphical models, the data formats and the way they interface to the user, \emph{e.g.} through the command-line arguments, a configuration file or a function in a specific programming language.
Another complication is that the number of computations increase rapidly in comparison studies, in particular when input parameters are altered, which requires parallel computation capabilities and careful bookkeeping to record and organise the results for reporting.
As a result, researchers may invest considerable time on tailored benchmarking scripts that intend to address these issues on a subset of algorithms.
To complicate matters beyond the implementational issues, there is no single well established metrics to evaluate the performance of structure learning algorithms.
Quite the contrary there is a wide range of different metrics to choose from \citep[see \emph{e.g.}][]{tsamardinos2006max,constantinou2021large} and researchers tend to use only a small
subset to evaluate a new algorithm, typically with the aim of highlighting the problem their algorithm addresses.

\felix{Identify problem of comparing }
With the problems described above in mind, the objective of this work is to develop a unified framework to facilitate execution and the benchmarking procedure of different algorithms.
We present a novel workflow called \hmm{\emph{b}}\emph{enchpress}\xspace, which is implemented in \emph{snakemake} \citep{10.1093/bioinformatics/bts480} and designed for making reproducible, platform-independent, and fully scalable benchmarks of structure learning algorithms.
\hmm{\emph{b}}\emph{enchpress}\xspace is interfaced via an easy-to-use configuration file in JSON-format \citep{pezoa2016foundations},
where a benchmark setup is specified in a module based \emph{probabilistic programming} style, separating model specification, in terms of graph and corresponding parameters, from algorithm execution and evaluation.
\hmm{\emph{s}}\emph{nakemake}\xspace enables \hmm{\emph{b}}\emph{enchpress}\xspace to seamlessly scale the computations to server, cluster, grid and cloud environments, without any changes in the configuration file.
The support for \emph{Docker} \citep{merkel2014docker}, together with platform-independent representations of data sets and graphs enables \hmm{\emph{b}}\emph{enchpress}\xspace to compare algorithms from different libraries, possibly implemented in different programming languages with different dependencies.

\hmm{\emph{b}}\emph{enchpress}\xspace is implemented in a modular coding style which facilitated for researchers to contribute with additional modules for generating models, structure learning algorithms and evaluating results.
In the first publicly released version we already included 19 algorithms from some of the most popular libraries such as \emph{BiDAG} \citep{suter2021bayesian}, \emph{bnlearn} \citep{JSSv035i03},   \emph{GOBNILP} \citep{cussens2020gobnilp}, \emph{pcalg} \citep{kalisch2012causal}, \emph{r.blip} \citep{scanagatta2015learning}, \emph{scikit-learn} \citep{scikit-learn}, \emph{TETRAD} \citep{glymour1986causal}, and \emph{trilearn} \citep{olsson2019} along with some of the models and data sets from the standard literature.

\felix{The solution}

Recently \cite{constantinou2020bayesys} developed \emph{Bayesys}, a software with similar intentions, albeit exclusively designed for \emph{Bayesian networks}, a specific type of probabilistic graphical models \citep{pearl1997probabilistic}.
Furthermore \emph{Bayesys} is a \emph{Java} implementation only including one algorithm for structure learning, called \emph{SaiyanH} \citep{constantinou2020learning}, suited for high dimensional settings.
\felix{Relate to other attempts to solve the same problem.}
However, the functionality of causal-compare is restricted to the algorithm implementations which are part of the \emph{TETRAD} project, which represents only a subset of available algorithms \citep{JMLR:v21:19-773}.
This is a huge restriction compared to \hmm{\emph{b}}\emph{enchpress}\xspace, which only requires a Linux system while there is no limitation for the individual algorithms requirements.

\hmm{\emph{b}}\emph{enchpress}\xspace is distributed under the Apache 2.0 Licence, encouraging its use in both academic
and commercial settings.
Source code and documentation are available from \url{https://github.com/felixleopoldo/benchpress}.

The rest of the paper is structured as follows.
Section~\ref{sec:pgms} presents an introduction to graphical models along with some notational conventions.
Section~\ref{sec:structure_learning} provides background on current strategies for structure learning.
Section~\ref{sec:json} describes the structure of the JSON interface and the modules already available for benchmarking.
Installation and usage guidelines are provided in Section~\ref{sec:installation}.
In Section~\ref{sec:case_study}, we show how to use \hmm{\emph{b}}\emph{enchpress}\xspace in four typical benchmarking scenarios.
\felix{Paper structure}

\section{Probabilistic graphical models}\label{sec:pgms}
This section provides a short introduction to graphical models and graph theory, with definitions as in  \emph{e.g.} \cite{diestel2005graph, lauritzen1996}.

Let \(\mathbb P\) be a \(p\)-dimensional distribution and let \(Y := (Y_i)_{i=1}^p\) be a random vector such that \(Y \sim \mathbb P\).
Further let  $\graph=(V, E)$ be a graph, where  \(V=\{Y_i\}_{i=1}^p\) is the node set and \(E\) is the edge set, consisting of tuples of distinct elements of \(V\).
Then \(\mathbb P\), or alternatively \(Y\), is said to be Markov with respect to \(\graph\) if separation, based on some graphical criteria, of two subsets of nodes,  $Y_A$ and $Y_B$, by a third set $Y_C$ in \(\graph\) implies their conditional independence in \(\mathbb P\), where $Y_A := (Y_i)_{i\in A}$.
The term probabilistic graphical model or graphical model has been used interchangeably in the literature to refer to either the tuple \((\graph, \mathbb P)\) or a collection of distributions \( \mathcal P_{\graph}\), usually restricted to some specific parametric family, where every \(\mathbb P' \in \mathcal P_{\graph} \) is Markov with respect to \(\graph{}\).
In this text we will use the former meaning.
In the case where \(\mathbb P\) is restricted to a specific parametric family, it is usually uniquely determined by some parameter vector \(\Theta\), thus either notation may be used.
\felix{Define graphical models}

Assuming \(f \) is a density for \(\mathbb P\), conditional independence between the random variables  $Y_A$ and $Y_B$ given $Y_C$ is well-defined as \(f(y_A|y_B, y_C) = f(y_A|y_C),\) where \( (y_A,y_B, y_C) \in \mathcalY_A
\times \mathcalY_B \times \mathcalY_C \) and \(\mathcalY_A\) denotes the range of \(Y_A\).
In contrast, node separation may have different definitions depending on the type of the graph considered.
Typical classes of graphs are: \emph{undirected graphs}, where edges are specified as unordered tuples and separation of  $Y_A$ and $Y_B$  given $Y_C$ means that every path between $Y_A$ and $Y_B$ must pass $Y_C$; \emph{directed acyclic graphs} (DAGs), where edges are represented as ordered tuples and separation is defined by means of a concept known as \emph{d-separation}, see \emph{e.g.}  \cite{pearl1997probabilistic}.
Next follows a short introduction to the type of graphical models currently implemented in \hmm{\emph{b}}\emph{enchpress}\xspace.

\felix{Examples of graphs: Bayesian networks}
\emph{Bayesian networks} constitute one of the most popular classes of graphical models, referring to distributions that are Markov with respect to a DAG.
The conditional independence encoded by the DAG implies that the density can be factorised into local conditional densities, where each node \(Y_i\) is directly dependent on its \emph{parents}, \( pa(Y_i)\), as
\begin{align*}
    f(y) = \prod_{i=1}^p f(y_i | pa(y_i)),
\end{align*}
where \(f(y_i | pa(y_i))\) denotes the conditional density of \(Y_i\) given \(pa(Y_i)=pa(y_i)\).
This property is appealing for several reasons.
Firstly, the partition into local conditional probability distribution provides an intuitive way to express and visualise knowledge about a specific problem domain.
In particular, an expert could express a causal mechanism in terms of a DAG, which in turn provides a mathematical language to answer causal queries from non-experimental data \citep{pearl1995causal}.
Further, it enables fast computations of conditional and marginal distribution using message-passing algorithms by exploiting the factorisation property in the distribution.
However, when the DAG is not used for expressing causal knowledge, the direction of the edges may be misleading since the conditional independencies of a distribution are not uniquely encoded by a single DAG in general but rather by a class of \emph{Markov equivalent} DAGs, all of which encode the same conditional independence statements.
One way to represent a Markov equivalence class is by a \emph{completed partially directed acyclic graph} \citep[CPDAG, see][]{chickering2013transformational}.
A CPDAG is a partially directed DAG which can be obtained from a DAG by regarding an edge \((Y_i,Y_j)\) as undirected if and only if both directed edges \((Y_i,Y_j)\) and \((Y_j,Y_i)\) are present in some of the DAGs in the same equivalence class.
The CPDAG is also sometimes called the \emph{essential graph}.

\felix{Introduce Markov networks}

\emph{Markov networks} or \emph{Markov random fields} are another popular class of graphical models, referring to distributions that are Markov with respect to an undirected graph.
There is in general no decomposition into local densities for the models in this class.
However, there exist functions \(\{\phi_{C}\}_{C\in \mathcal{C}}\) such that
\begin{align*}
        f(y) = \prod_{C \in \mathcal{C}} \phi_{C}(y_C),
\end{align*}
where \(\mathcal{C}\) is the set of  maximal subgraphs where all nodes are directly connected with each other, usually called \emph{cliques}, which  may be leveraged to calculate marginal distributions efficiently.
There is a sub-class of undirected graphs called \emph{decomposable} (DG) or \emph{chordal} or \emph{triangulated}, where each cycle with more than three nodes has a \emph{chord}, an edge which breaks it.
Distributions which are Markov with respect to decomposable graphs are called \emph{decomposable} and their densities are built by local, clique-specific densities, \(\{f_C\}_{C \in \mathcal{C}}\) and their marginals so that the joint density factorises as
\begin{align*}
    f(y) = \frac{\prod_{C \in \mathcal{C}} f_C(y_{C})}{\prod_{S \in \mathcal{S}} f_S(y_S)},
\end{align*}
where \(\mathcal{S}\) is the multi-set of separators in \(\graph\), found \emph{e.g.} as the intersection of neighboring cliques in a so-called \emph{junction tree}, built from the cliques in the graph; see \emph{e.g.} \cite{lauritzen1996} for details.
In general Markov networks and Bayesian networks encode different sets of conditional independence statements, and there are sets which we can represent with one model but not the other.
Decomposable graphs represent the intersection of those sets of statements which we can represent  by both undirected graphs and DAGs \citep[see \emph{e.g.}][]{lauritzen1996, cowell2003probabilistic, koller2009probabilistic}.
The DAGs representing the intersection set are called \emph{perfect}, meaning that all the parents of each node are connected by an edge.
Further characterisations of decomposable graphs are found in \emph{e.g.} \cite{lauritzen1996, cowell2003probabilistic} and \cite{duarte2021new}.

\felix{Introduce Markov networks}

\felix{Summarize the above discussion}

\section{Structure learning}
Structure learning algorithms for Bayesian networks was proved to be NP-hard, even when the number of parents of each node is restricted to two \citep{chickering2004large}.
Also the number of undirected graphs grows as \( \mathcal O(2^{p^2})\) even when restricted to decomposable graphs, making an exhaustive search for a specific graph infeasible \citep[see \emph{e.g.}][]{countchordal, cta}.
Therefore, in practice we must rely on approximation methods which we can divide into three main categories: \emph{constraint-based}, \emph{score-based}, and \emph{hybrid} algorithms, discussed briefly below.

\felix{Intro to structure learning}\label{sec:structure_learning}

The seed to the first constraint-based algorithm was planted by \cite{verma1991equivalence} in the context of causal Bayesian networks.
In line with their work, constraint-based methods employ conditional independence tests among the variables to first estimate an undirected graph by, \emph{e.g.}, starting with the complete graph, and excluding an edge \((Y_i,Y_j)\) if there exists some set \(Y_A\) such that the hypothesis \(Y_i  \ci_{\mathbb P} Y_j \mid Y_A\) can be not rejected according to a suitable test procedure.
Such methods tend to be very fast but testing can suffer from large numbers of false negatives, \emph{i.e.} falsely not including edges when they should be present.
Relaxing the significance level of the tests to reduce false negatives can however rapidly inflate false positives and runtimes.
\felix{Intro to constraint based methods}

Score-based methods on the other hand aim to optimise a global score function defined on the graph space.
As a work around to cope with the immense number of graphs, the search is often limited to certain types of graphs.
For example in  DAG models, one may restrict the number of parents or combine DAGs into bigger and more easily represented classes \citep[\emph{e.g.} orders and partitions][]{Friedman2003, doi:10.1080/01621459.2015.1133426}.
For decomposable graphs a possibility is to limit the maximal clique size.
Score functions may rely on penalised log-likelihoods or on Bayesian posterior probabilities of graphical structures conditional on the observed data \citep[see \emph{e.g.}][]{koller2009probabilistic}.
Besides mere optimisation, Bayesian structure learning also focuses on sampling procedure to provide a finer characterisation of the posterior graph distribution.
When an expression for the posterior distribution is available at least up to a normalising constant, we can implement  Markov chain Monte Carlo (MCMC) schemes to sample from it, see \emph{e.g.} \cite{madigan1995graphicalmodels, giudici2003improving, Friedman2003} for early works and \cite{ doi:10.1080/01621459.2015.1133426, jennings2018birth} for more recent strategies.
In Gaussian Markov networks, Monte Carlo sampling turns out to be infeasible for larger network since there is in general no analytic expression for the un-normalised posterior \citep{Atay-Kayis01062005}.
On the other hand, score-based methods such as \emph{graphical lasso}, relies on the fact that the edges in the graph are revealed from the non-zero pattern of the elements in the precision (inverse covariance) matrix.
Optimisation techniques can successfully estimate the \(L_1\)-penalised log-likelihood of the precision matrix, and the graph as a by-product; see \cite{citeulike:2134265, Meinshausen2008880, doi:10.1198/jcgs.2011.11051a}.
Decomposable models, do possess a closed form expression for the un-normalised graph posterior and most inferential methods built on the representation via undirected graphs focus on Monte Carlo sampling \cite[see \emph{e.g.}][]{Giudici01121999, Green01032013, elmasri2017decomposable, olsson2019} with a few focusing on optimisation \cite[see \emph{e.g.}][]{carvalho2006structure, studeny2017towards, rantanen2017learning}.

\felix{Add something.. One example of decomposable scores is  hyper directed Markov property .}

\felix{Intro to score methods}

Hybrid algorithms typically use constraint-based algorithms as a preliminary step called \emph{restrict phase} in order to reduce the search space for a score-based method.
A score based-algorithm is then used in the \emph{maximise phase}  \citep[see \emph{e.g.}][]{tsamardinos2003algorithms, scanagatta2015learning,  kuipers2018efficient}.

\felix{Intro to hybrids}

\felix{Summarize methods}
\section{The JSON configuration file and supported modules}\label{sec:json}

In this section we describe the structure of the JSON \citep{pezoa2016foundations} configuration file which serves as interface to \hmm{\emph{b}}\emph{enchpress}\xspace.
The JSON file provides a simple \emph{key-value} format which is easy to interpret and modify.
We will use the \texttt{typewriter} style to indicate keys and \vals{bold italic} style for values.
When the value of a key is another JSON object or a list, we will sometimes refer to it as a section.
At the highest level there are two main sections, \texttt{benchmark\_setup} and  \texttt{resources}.

The \texttt{resources} section contains separate subsections of the available \emph{modules} for generating or defining graphs (\texttt{graph}), parameters (\texttt{parameters}), data (\texttt{data}), and algorithms for structure learning (\texttt{structure\_learning\_algorithms}).
Each module in turn has a list, where each of its elements are objects defining a parameter setting, identified with a unique  \texttt{id}.
See Listings~\ref{alg:pc} for an example object in the \texttt{pcalg\_pc} module.
The \texttt{benchmark\_setup} section specifies the data models (\texttt{data}) and evaluation methods (\texttt{evaluation}) to be considered for analysis.
The module objects used in \texttt{benchmark\_setup} are those defined in \texttt{resources} and referenced by their corresponding \texttt{id}'s.
The output of each module is saved systematically under the \emph{results/} directory based on the corresponding objects values.
A more detailed description of the structure of each section and how to use the modules supported so far follows in the next section.

It is important to note that one of the main strengths of \hmm{\emph{b}}\emph{enchpress}\xspace is that the modular design provides a flexible framework to integrate new modules into any of the sections.
To do so, one essentially needs either a \emph{Docker} image or a local installation, where the new module is installed and obeys the input and output formats specified in Appendix~\ref{appendix}.
The \hmm{\emph{s}}\emph{nakemake}\xspace workflow should then be adapted for the new module following the same structure as the existing modules in its section.
The reader is referred to \cite{boettiger2017introduction, 10.1371/journal.pcbi.1008316} for a comprehensive introduction to \emph{Docker} in data science.

\subsection{Benchmark setup (\lowercase{\texttt{benchmark\_setup}})}
\subsubsection{Data setup section (\lowercase{\texttt{data}})}
\hmm{\emph{b}}\emph{enchpress}\xspace supports different strategies for defining models and data.
One possibility is to generate data, graphs and parameters using the existing modules, referenced by  their objects \texttt{id}, as described above.
An alternative is to  provide user-specified  graphs, parameters and data sets, referenced by their corresponding filenames.
Further we can combine the two approaches as needed.

The \texttt{data} section consists of a list of objects,
 each containing the \texttt{id} of a graph (\texttt{graph\_id}), parameters (\texttt{parameters\_id}), data (\texttt{data\_id}) module object, and a tuple (\texttt{seed\_range}), specifying the range of the random seeds used in the modules.
In case of fixed data sets, parameters, or graphs, we simply use the corresponding filenames in place of  \texttt{id}'s

For each seed number \(i\) in the range defined by \texttt{seed\_range}, a graph \(\graph_i\) is obtained as specified according to \(\texttt{graph\_id}\).
Given \(\graph_{i}\), the parameters in the model $\Theta_i$ are obtained according to the specifics in \texttt{parameters\_id}.
A data matrix, $\mathbf Y_i^T=(Y_{1:p}^j)_{j=1}^n$, is then sampled from $(\graph_i,\Theta_i)$ as specified by the data model in \texttt{data\_id}.
More details and examples are given later in Section~\ref{sec:resources}.

\subsubsection{Evaluation section (\texttt{evaluation})} \label{sec:evaluation}
There is typically no algorithm which performs well with respect to all performance metrics, and no single metric which is generally preferred.
Therefore, to get an overall picture of the performance of an algorithm \hmm{\emph{b}}\emph{enchpress}\xspace supports different metrics.
This section describes the metrics implemented in \hmm{\emph{b}}\emph{enchpress}\xspace, where we let \(\graph=(V,E)\) and \(\graph'=(V',E')\) denote the true and the estimated graph respectively.

The first choice is the space on which to compare the true and inferred structures. The currently supported Markov-equivalent types are CPDAGs, \emph{pattern graphs} and the skeleton of the graphs (undirected version of a graph). The pattern graph of a DAG is the graph obtained by keeping the directions of the \emph{v-structures}, and regarding the remaining edges as undirected. v-structures are characterised by triples of nodes where the edges of two of the nodes, which should not be neighbors, are pointing to the third. Within a given comparison space, the \texttt{evaluation} section below contains a list of the currently available evaluation modules.
\\\\
\noindent \emph{Structural Hamming distance} \vspace{0.2cm}\\
The \emph{structural Hamming distance} (SHD) is one of the most commonly used metrics to compare graphs.
It describes the number of changes, in terms of adding, removing or reversing edges, needed so that \(E=E'\).
\\\\
\noindent \emph{Receiver operating characteristic} (\texttt{roc}) \vspace{0.2cm}\\
The \texttt{roc} module defines a metric for comparing graphs which may also be applied to \emph{mixed graphs}, having a combination of edges oriented in different ways (directed or undirected), thus we distinguish between the orientations.

We assign to every edge \(e \in E'\) the true positive score \(TP(e)=1\) if \(e\) is contained in \(E\) with the same orientation (\(e\) being undirected in both \(E\) and \(E'\) or \(e\) having the same direction in both \(E\) and \(E'\)), \(TP(e)=1/2\) if \(e\) is contained in \(E\) with a different orientation, otherwise \(TP(e)=0\).
The false positive score \(FP(e)\) is defined analogously.
First let \( \bar E\) denote the complementary set of \(E\) in the sense that directed edges are reversed and undirected edges (non-edge) are removed (added).
For every edge \(e \in \barE'\), \(FP(e)=1\) if \(e\) is contained in \( \bar E\) with the same orientation
and \(FP(e)=1/2\) if \(e\) is contained in \( \bar E\)  with a different orientation, otherwise \(FP(e)=0\).

The total true and false positive edges (TP and FP) are obtained  as the sum of the individual edge scores, \emph{i.e.}
\begin{align*}
    TP := \sum_{e\in E} TP(e), \quad FP := \sum_{e\in E} FP(e).
\end{align*}
We define the true and false positive rates (FPRp and TPR) as
\begin{align*}
FPRp := \frac{FP}{P},\quad TPR := \frac{TP}{P},
\end{align*}
where \(P:=|E|\) denotes the total number of edges in \(\graph\).
Note that TP and FP reduces to the ordinary counting true and false positives when all edges are undirected in both \(\graph\) and \(\graph'\).

The \texttt{roc} module evaluates the TPR and FPRp for each seed in \texttt{seed\_range} and the median (mean) FPRp is plotted against the median (mean) TPR in ROC type sub figures using \emph{ggplot2} \citep{tidyverse} based on the graph, parameters and data modules used.
The top three title rows of each sub figure describe, from top to bottom, the objects from the graph, parameters and data modules.
The varying parameter forming the curve should be given as a list of values in the corresponding algorithm object.
A parameter could for example be the significance level in a test procedure of a constraint-based algorithm, such as \texttt{alpha} in Listing~\ref{alg:pc}, or a parameter of the score function in a score-based method.

The \texttt{id}'s for the objects that should be part of the evaluation should be placed in a list called \texttt{ids}.
The \texttt{filename\_prefix} field specifies a prefix of the produced files.
The user can choose to show points (\texttt{points}), parameter values (\texttt{text}), lines between the points or values (\texttt{path}), and error bars showing the 5\% and 95\% quantiles of the TPR (\texttt{errorbar}).

There is an interesting connection between the SHD metric and the ROC curve when TPR and FPRp are defined as above. \felix{move to graphical models section on CPDAGS, skeleton+vstructures.. = ME}
In such (possibly mixed) graphs, the SHD can easily be determined in a ROC figure as the scaled \emph{Manhattan distance} between the points (0,1) and (TPR, FPRp), \emph{i.e.}
\begin{align*}
    SHD/P= 1-TPR + FPRp.
\end{align*}
Consequently, the ROC curves enable the simultaneous analysis of SHD, TPR, and FPRp in the same figure.

The benchmarks are saved in \emph{.csv} files in the directory \emph{results/benchmarks/}.
A summary of the benchmarks for each algorithm is saved directly in \emph{results/output} according to the algorithm's module name and the value of the \texttt{filename\_prefix} field.
For example, if  \texttt{filename\_prefix} is set to \vals{myplots/}  the results of the \texttt{tetrad\_gfci} module are stored in \emph{results/output/myplots/tetrad\_gfci.csv}.
The plots are saved under \emph{results/output/myplots/} along with the source file \emph{results/output/myplots/ROC\_data.csv}.
\\\\
\noindent \emph{Plot true adjacency matrix} (\texttt{adjmat\_true\_plots}) \vspace{0.2cm}\\
The \texttt{adjmat\_true\_plots} module plots the adjacency matrices of the true graphs.
The plots are saved  in sub directories of \emph{results/adjmat/} and copied to \emph{results/output/adjmat\_true\_plots/} for easy reference.
\\\\
\noindent \emph{Plot adjacency matrix} (\texttt{adjmat\_plots}) \vspace{0.2cm}\\
The \texttt{adjmat\_plots} module plots the adjacency matrices of the estimated  graphs.
The figures are saved in \emph{results/adjmat} and copied to \emph{results/output/adjmat\_plots/}.
\felix{Not clear if I should couple MCMC and Baysian methods as one..}
\\\\
\noindent \emph{Plot estimated graphs} (\texttt{graph\_plots}) \vspace{0.2cm}\\
The \texttt{graph\_plots} module plots and saves the estimated graphs in \emph{.dot} format.
The figures are saved in \emph{results/adjmat} and copied to \emph{results/output/graph\_plots/}.
\\\\
\noindent \emph{Plot true graphs} (\texttt{graph\_true\_plots}) \vspace{0.2cm}\\
The \texttt{graph\_true\_plots} module plots the true underlying  graphs.
The figures are saved in \emph{results/adjmat} and copied to \emph{results/output/graph\_true\_plots/}.
\\\\
\noindent \emph{MCMC mean graphs} (\texttt{mcmc\_heatmaps}) \vspace{0.2cm}\\
As mentioned in Section~\ref{sec:structure_learning}, for Bayesian inference it is custom to use MCMC methods to simulate a Markov chain of graphs  \(\{\graph^l\}_{l=0}^\infty\) having the graph posterior as stationary distribution.
Suppose we have a realisation of length \(M+1\) of such chain \(\{\graph^l\}_{l=0}^M\), then the posterior edge probability of an edge \(e\) is estimated by \(\frac{1}{M+1-b} \sum_{l=b}^{M} \mathbf{1}_{e}(e^l)\), where the first samples \(b\)  are disregarded as a burn-in period.

The \texttt{mcmc\_heatmaps} module has a list of objects, where each object has an \texttt{id} field for the algorithm object \texttt{id} and a burn-in field (\texttt{burn\_in}) for specifying the burn-in period.
The estimated probabilities are plotted in heatmaps using \emph{seaborn} \citep{Waskom2021} which are saved in \emph{results/mcmc\_heatmaps/} and copied to \emph{results/output/mcmc\_heatmaps/} for easy reference.
\\\\
\noindent \emph{MCMC trajectory plots} (\texttt{mcmc\_traj\_plots}) \vspace{0.2cm}\\
The \texttt{mcmc\_traj\_plots} module plots the so called score values in the trajectory or the value of a given functional.
The currently supported functionals are the number of edges for the graphs (\vals{size}) and the graph \vals{score}.
The \texttt{mcmc\_traj\_plots} module has a list of objects, where  each object has an \texttt{id} field for the algorithm object \texttt{id}, a burn-in field (\texttt{burn\_in}) and a field specifying the functional to be considered (\texttt{functional}).
Since the trajectories tend to be very long, the user may choose to thin out the trajectory by only considering every graph at a given interval length specified by the \texttt{thinning} field.
The plots are saved in \emph{results/mcmc\_traj\_plots/} and copied to \emph{results/output/mcmc\_traj\_plots/}.
\\\\
\noindent \emph{MCMC auto-correlation plots} (\texttt{mcmc\_autocorr})\vspace{0.2cm}\\
The \texttt{mcmc\_autocorr} module plots the auto-correlation of a functional of the graphs in a MCMC trajectory.
Similar to the \texttt{mcmc\_traj\_plots} module, the  \texttt{mcmc\_autocorr} module has a list of objects, where each object has an \texttt{id}, \texttt{burn\_in}, \texttt{thinning}, and a \texttt{functional} field.
The maximum number of lags after thinning, is specified by the \texttt{lags} field.
The plots are saved in  \emph{results/mcmc\_autocorr/} and copied to \emph{results/output/mcmc\_autocorr/}.
\\\\

\subsection{Resources section (\texttt{resources})}\label{sec:resources}

\subsubsection{Graph, \(\graph\) (\texttt{graph})} \label{sec:graph}
Several possibilities exist for defining the graph component, depending on whether we wish to draw random graphs from a given distribution, or maybe provide a given structure on which we wish to perform a benchmarking study.
Below is the list of currently available modules.

The adjacency matrices of the generated graphs are saved in the directory \emph{results/adjmat/}, where each generation method has its own sub directory, named by their module name and structured by their objects values.
\\\\
\noindent \emph{Random DAG with constrained number of parents per node} (\texttt{pcalg\_randdag}) \vspace{0.2cm}\\
An object of the \texttt{pcalg\_randdag}\, module specifies a graph generated using the \emph{randDAG} function from the \emph{pcalg} package \citep{kalisch2012causal}, with the extra functionality of restricting the maximal number of parents per node (\texttt{max\_parents}).
It samples a random graph with \texttt{n} nodes and \texttt{d} neighbours per node on average, using the algorithm specified by \texttt{method} with parameters \texttt{par1} and \texttt{par2}, for further details see \cite{kalisch2012causal}.
\\\\
\noindent \emph{Band matrix} (\texttt{bandmat}) \vspace{0.2cm}\\
An object of the \texttt{bandmat} module specifies a decomposable graphs with adjacency matrix with given bandwidth (\texttt{bandwith}).
\\\\
\noindent \emph{Band matrix with random width} (\texttt{rand\_bandmat}) \vspace{0.2cm}\\
An object of the \texttt{max\_bandwith}\, module specifies a decomposable graph by randomly sampling a band structured adjacency matrix of a given maximum width (\texttt{max\_bandwith}).
\\\\
\noindent \emph{Christmas tree algorithm} (\texttt{trilearn\_cta}) \vspace{0.2cm}\\
An object of the  \texttt{trilearn\_cta}\, module defines decomposable graph randomly generated using the  \emph{Christmas tree algorithm} \citep{cta}, with the sparsity parameters \texttt{alpha} and \texttt{beta}.
\\\\
\noindent \emph{Fixed adjacency matrix} \vspace{0.2cm}\\
A fixed graphs is represented as an adjacency matrix and should be formatted as specified in Appendix~\ref{appendix}.
The file should be stored with the \emph{.csv} extension in the directory \emph{resources/adjmat/myadjmats/}.

\subsubsection{Model parameters, \(\Theta\) (\texttt{parameters})}\label{sec:parameters}
Depending on how the user chooses to provide data and graph there are different possibilities for setting the model parameters.
As for the graph module, one possibility is to define parameters is by a filename.
Another is to use one of the modules listed below.
In any case, it is assumed that the underlying graph \(\graph\) is the one specified by \texttt{graph\_id} in the \texttt{data} section.

The model parameters are saved in the directory \emph{results/parameters/}, where each generation method has its own sub directory named by their module and structured by the object values and corresponding graph.
\\\\
\noindent \emph{Random binary Bayesian network} (\texttt{bin\_bn}) \vspace{0.2cm}\\
Provided that \(\graph\) is a DAG, an object of the \texttt{bin\_bn}\, module defines a binary Bayesian network (where the nodes are binary variables) by sampling its conditional probability tables.
More specifically, the conditional probability tables for each variable \(Y_i\) and parent configuration \(pa(y_i)\) are distributed as
\begin{align*}
    pr(Y_i=0 | pa(y_i) ) \sim \mathrm{Unif}([a, b]),
\end{align*}
where  \((a,b) \in [0,1]^2, a<b\) and \( \mathrm{Unif}(c)\) denotes the uniform distribution on the set \(c\).

Using an \texttt{id} of this  module in the \texttt{parameters\_id} field of the \texttt{data} section requires that \texttt{graph\_id} represents a DAG (DAG module object or DAG adjacency matrix).
\\\\
\noindent \emph{Random linear Gaussian Bayesian network} (\texttt{sem\_params}) \vspace{0.2cm}\\
Provided that \(\graph\) is a DAG, an object of the
\texttt{sem\_params} module defines a Bayesian network by generating the weight matrix \(W\) of a Gaussian linear structural equation model (SEM) of the form
\begin{align}
Y_i=\sum_{j:Y_j\in pa(Y_i)} W_{ij}Y_j + Z_i,
\label{eq:sem}
\end{align}
where \(Z_i\sim \mathcal N(\mu, \sigma^2)\) and elements of \(W\) are distributed as
\begin{align*}
W_{ij} \sim
\begin{cases}
\mathrm{Unif}([a, b])\mathrm{Unif}(\{-1,1\}) & \text{ if }(i, j) \in E\\
0 & \text{ otherwise.}
\end{cases}
\end{align*}

Using an object \texttt{id} of this module  in the \texttt{parameters\_id} field of the \texttt{data} section requires that \texttt{graph\_id} represents a DAG.
\\\\
\noindent \emph{Fixed Bayesian network} \vspace{0.2cm}\\
A Bayesian network object in the \emph{R}-package \emph{bnlearn} \citep{JSSv035i03} is an instance of the \emph{bn.fit} class and contains both the DAG and corresponding parameters.
Motivated by the modular design of \hmm{\emph{b}}\emph{enchpress}\xspace, a \emph{bn.fit} object may be used to specify only the parameters  of a Bayesian network, by storing it in \emph{.rds} format in \emph{resources/parameters/myparams/bn.fit\_networks/}.
The adjacency matrix of the DAG should be stored in the folder for graphs, as specified in Section~\ref{sec:graph} and \texttt{graph\_id} should be set to the filename.
\hmm{\emph{b}}\emph{enchpress}\xspace already contains a number of standard Bayesian network models, \emph{e.g.} those available at the website of \emph{bnlearn}, including the so-called \emph{Bayesian network repository} (UAI)  \citep{friedman1997impact}.

A matrix defining the weights \(W\) in a SEM as defined in \eqref{eq:sem} can be stored  in \nolinebreak \emph{resources/parameters/myparams/sem\_params/} with the \emph{.csv} extension, formatted similarly to the adjacency matrix specified in Appendix~\ref{appendix}.
This requires that \texttt{graph\_id} represents a DAG.
\\\\
\noindent \emph{Gaussian graph intra-class model} (\texttt{intraclass}) \vspace{0.2cm}\\
Provided that \(\graph\) is decomposable, an object of the \texttt{intraclass} module defines
a zero mean multivariate Gaussian distribution by its covariance matrix $\Sigma$ as
\begin{align*}
        \Sigma_{ij} = \begin{cases}
        \sigma^2, &\text{ if } i=j\\
        \rho\sigma^2, &\text{ if } (i,j) \in E \\
    \end{cases}
\end{align*}
and $(\Sigma^{-1})_{ij} =0$ if $(i,j) \notin E$, where \(\sigma^2>0\) (\texttt{sigma2}) and \(\rho \in [0,1]\) (\texttt{rho}) denote the variance and correlation coefficient, respectively.

Using an object \texttt{id} of this module in the \texttt{parameters\_id} field of the \texttt{data} section requires that \texttt{graph\_id} represents a decomposable graph.
\\\\
\noindent \emph{Hyper Dirichlet} (\texttt{trilearn\_hyper-dir}) \vspace{0.2cm}\\
Provided that \(\graph\) is decomposable, an object of the \texttt{trilearn\_hyper-dir} module defines a categorical decomposable model by sampling the parameters from the \emph{hyper Dirichlet} distribution \citep{dawid1993}, with specified concentration parameter (\texttt{pseudo\_obs}) and number of levels per variable (\texttt{n\_levels}).\\\\

\noindent \emph{Hyper inverse-Wishart} (\texttt{trilearn\_hyper\_inv\_wishart}) \vspace{0.2cm}\\
Provided that \(\graph\) is decomposable, an object of the \texttt{trilearn\_hyper\_inv\_wishart} module defines a zero-mean Gaussian decomposable model by sampling the covariance matrix from the \emph{Hyper inverse-Wishart} distribution \citep{dawid1993} according to the procedure outlined by \cite{Carvalho01082007}.
The clique-wise scale matrices are fixed to the identity matrix and the degrees of freedom subtracted by \(p\) (\texttt{dof}) should be specified by the user.
\\\\
\subsubsection{Data, \(\mathbf Y\) (\texttt{data})}
There are two possibilities for the data section.
The user may provide their own data in \emph{.csv} format for benchmarking on real data.
Alternatively \hmm{\emph{b}}\emph{enchpress}\xspace can generate data according to a data model for benchmarking on simulated data.
The data sets are saved in \emph{results/data/}, structured according to the generation modules of the graph and parameters and their object values or filenames.
\\\\
\noindent \emph{Fixed data sets} \vspace{0.2cm}\\
The are two ways of providing fixed data sets.
The first option is to place data files directly in the directory \emph{resources/data/mydatasets/} with the \emph{.csv} extension, formatted as specified in Appendix~\ref{appendix}.
The second option is to place the files in a sub directory of \emph{resources/data/mydatasets/}.
In the latter case, the \texttt{data\_id} field should be the name of the directory and all the files in it will be considered for evaluation.

Note that, when a filename is used as value for \texttt{data\_id}, since the data is fixed, \texttt{seed\_range} should be set to \vals{null}.
However, depending on the evaluation method used, the \texttt{graph\_id} could be either \vals{null} or the filename of the true graph underlying the data.
For example, the \texttt{roc} module requires that \texttt{graph\_id} set to a proper \texttt{id} for a graph object, however it is not a requirement for the \texttt{adjmat\_plots} module.
When the value of \texttt{data\_id} is a directory, the other fields should be set to \vals{null}.
\\\\
\noindent \emph{Independent identically distributed (i.i.d.) samples} (\texttt{iid}) \vspace{0.2cm}\\
An object of the \texttt{iid}\, module will draw a specified (\texttt{sample\_sizes}) number of independent samples according to the module \texttt{id}'s or filenames specified in \texttt{graph\_id} and \texttt{parameters\_id}.
Continuous data may be standardized by setting \texttt{standardized} to \vals{true}.
See \cite{reisach2021beware} for a discussion about standardising data in a structure learning context.

\subsubsection{Structure learning algorithms (\texttt{structure\_learning\_algorithms})}
This section contains a short summary, sorted in a temporal ordering based on publication year,  of the currently 19 structure learning modules supported by \hmm{\emph{b}}\emph{enchpress}\xspace that are also used in the simulation studies in Section~\ref{sec:case_study}.
See Table~\ref{tab:algorithms} for a summary of their properties and implementations.
The full list of available algorithms is found in the documentation of \hmm{\emph{b}}\emph{enchpress}\xspace.
\\\\
\noindent \emph{PC} (\texttt{pcalg\_pc})\vspace{0.2cm}\\
The PC algorithm \citep{spirtes1991algorithm}, is a constraint based method consisting of two main steps.
The first step is called the \emph{adjacency search} and amounts to finding the undirected skeleton of the DAG.
The second step amounts to estimating a CPDAG.
See Listing~\ref{alg:pc} for an example.
\begin{listing}
    \begin{minted}[frame=single,
                   framesep=3mm,
                   linenos=true,
                   xleftmargin=21pt,
                   tabsize=4]{js}
{
    "id": "pc-binCItest",
    "alpha": [0.01, 0.05, 0.1],
    "NAdelete": true,
    "mmax": "Inf",
    "u2pd": "relaxed",
    "skelmethod": "stable",
    "conservative": false,
    "majrule": false,
    "solveconfl": false,
    "numCores": 1,
    "verbose": false,
    "indepTest": "binCItest"
}
\end{minted}
\caption{Example object in the \texttt{pcalg\_pc} module.}
\label{alg:pc}
\end{listing}
\\\\
\noindent \emph{Fast causal inference} (\texttt{tetrad\_fci}) \vspace{0.2cm}\\
 \emph{Fast causal inference} (FCI) is a two-stage algorithm similar to the PC algorithm, where the first step finds the skeleton and the second step orients the edges in a DAG \citep{spirtes2000causation}.
The FCI is designed for models with latent variables.
\\\\
\noindent \emph{Hill climbing} (\texttt{bnlearn\_hc})\vspace{0.2cm}\\
\emph{Hill climbing} (HC) is a score-based algorithm which starts with a DAG with no edges and adds, deletes or reverses edges in a greedy manner until an optimum is reached \citep{russell2002artificial,scutari2019learning}.
\\\\
\noindent \emph{Tabu} (\texttt{bnlearn\_tabu})\vspace{0.2cm}\\
\emph{Tabu} is a less greedy version of the HC algorithm allowing for non-optimal moves that might be beneficial from a global perspective to avoid local maxima \citep{russell2002artificial, scutari2019learning}.
\\\\
\noindent \emph{Really fast causal inference} (\texttt{tetrad\_rfci})\vspace{0.2cm}\\
\emph{Really fast causal inference} (RFCI) is a modified version of the the FCI algorithm, where some additional tests before orienting \emph{v-structures} and discriminating paths are introduced, making it much faster, in particular for sparse graphs in high-dimension \citep{colombo2012learning}.
\\\\
\noindent \emph{Fast greedy equivalent search} (\texttt{tetrad\_fges})\vspace{0.2cm}\\
\emph{Fast greedy equivalent search} (FGES) is a score based method based on the  the greedy equivalence search (GES) of \cite{meek1997graphical}.
This algorithm operates on the space of CPDAG's \citep{chickering2002optimal}.
Its complexity is polynomial in the number of nodes.
\\\\
\noindent \emph{Greedy fast causal inference} (\texttt{tetrad\_gfci})\vspace{0.2cm}\\
The FGES is asymptotically correct under the unrealistic assumption that there are no unmeasured confounders \citep{ogarrio2016hybrid}.
However, FCI and RFCI often perform poorly in scenarios with small sample sizes.
 \emph{Greedy fast causal inference} (GFCI) is a hybrid method combining elements from the FCI and FGES algorithms, while making no assumptions about unmeasured confounders.
It is expected to perform better on small samples.\\
\\\\
\noindent \emph{Order Markov chain Monte Carlo} (\texttt{bidag\_order\_mcmc})\vspace{0.2cm}\\
This technique relies on a Bayesian perspective on structure learning, where the score of a DAG is defined as its posterior distribution.
To overcome the limitation of simple structure-based MCMC schemes, \cite{Friedman2003} turned to a score defined as the sum of the posterior scores of all DAG which are consistent with a given topological ordering of the nodes. One can then run a Metropolis-Hasting algorithm to sample from the distribution induced by the order score, and later draw a DAG consistent with the order.
This strategy substantially improves convergence with respect to earlier structure MCMC scheme, though it unfortunately produces a biased sample on the space of DAGs.
The implementation considered in \hmm{\emph{b}}\emph{enchpress}\xspace is a hybrid version with the sampling performed on a restricted search space initialised with constraint-based testing and improved with score-based search \citep{kuipers2018efficient}.
\\\\
\noindent \emph{Grow-shrink} (\texttt{bnlearn\_gs})\vspace{0.2cm}\\
The \emph{grow-shrink} (GS) algorithm is based on the Markov blanket of the nodes in a DAG. For a specific node, the Markov blanket it the set of nodes which conditioning upon renders it conditionally independent from all other variables \citep{margaritis2003learning}.
It is a constraint-based method which estimates the Markov blanket of a node in a two-stage forward-backward procedure using conditional independence tests.
The Markov blankets are used to first estimate an undirected graph and then estimate a DAG in a four-step procedure.
\\\\
\noindent \emph{Incremental association Markov blanket} (\texttt{bnlearn\_interiamb})\vspace{0.2cm}\\
As in GS, this algorithm is also based on the Markov blanket method to first determine the undirected skeleton.
However, in \emph{incremental association Markov blanket} (inter-IAMB) the variable to be included in the Markov blankets are not considered in static order as in GS and the forward-backward stages are combined into a single procedure, which has the effect of reducing the size of the blankets \citep{tsamardinos2003algorithms}.
\\\\
\noindent \emph{Max-min hill-climbing} (\texttt{bnlearn\_mmhc})\vspace{0.2cm}\\
 \emph{Max-min hill-climbing} (MMHC) is a hybrid method which first estimates the skeleton of a DAG using an algorithm called \emph{Max-Min Parents and Children} and then performs a  greedy hill-climbing search to orient the edges with respect to a Bayesian score \citep{tsamardinos2006max}.
It is a popular approach used as standard benchmark and also well suited for high-dimensional domains.
\\\\
\noindent \emph{Globally optimal Bayesian network learning using integer linear programming} (\texttt{gobnilp})\vspace{0.2cm}\\
A score based method using integer linear programming (ILP) for learning an optimal DAG for a Bayesian network with limit on the maximal number of parents for each node \citep{cussens2012bayesian}.
It is a two-stage approach where candidate parent sets for each node are discovered in the first phase and the optimal sets are determined in a second phase.
\\\\
\noindent \emph{Acyclic selection ordering-based search} (\texttt{rblip\_asobs})\vspace{0.2cm}\\
A score-based two-phase algorithm where the first phase aims to identify the possible parent sets, \cite{scanagatta2015learning, scanagatta2018approximate}.
The second phase performs an optimisation on a modification of the space of node orders introduced in \cite{teyssier2012ordering}, allowing edges from nodes of higher to lower orders, provided that no cycles are introduced.
\\\\
\noindent \emph{Iterative search} (\texttt{bidag\_itsearch})\vspace{0.2cm}\\
This is a hybrid score-based optimisation technique based on Markov chain Monte Carlo schemes \citep{kuipers2018efficient}. The algorithm starts from a skeleton obtained through a fast method (e.g. a constraint based method, or GES). Then it performs score and search on the DAGs belonging to the space defined by the starting skeleton. To correct for edges which may be missed, the search space is iteratively expanded to include one additional parent for each variable from outside the current search space. The score and search phase relies on an MCMC scheme producing a chain of DAGs from their posterior probability given the data.
\\\\
\noindent \emph{No tears} (\texttt{notears})\vspace{0.2cm}\\
This score-based method recasts the combinatorial problem of estimating a DAG into a purely continuous non-convex optimization problem over real matrices with a smooth constraint to ensure acyclicity \citep{zheng2018dags}.

\begin{center}
    \begin{table}
        \small
        \begin{center}
            \begin{tabular}{ l l l l l }
                Algorithm & Space & Language & Package & Type \\
                \hline
                Iterative search & DAG&R & \emph{BiDAG} & H \\
                Order MCMC &DAG & R & \emph{BiDAG} & S \\
                GS & DAG &R& \emph{bnlearn} & C\\
                MMHC & DAG &R& \emph{bnlearn} & H \\
                HC & DAG &R& \emph{bnlearn} & S\\
                Inter-IAMB &CPDAG & R&\emph{bnlearn} & C\\
                Tabu & DAG &R& \emph{bnlearn} & S \\
                GOBNILP &DAG &C& \emph{GOBNILP} & S \\
                No tears & DAG &Python& \emph{notears} & S \\
                PC &CPDAG &R& \emph{pcalg} & C \\
                ASOBS & DAG &R/Java& \emph{r.blip} & S \\
                FCI& DAG  &Java& \emph{TETRAD} & C \\
                GFCI& DAG  &Java& \emph{TETRAD} & H \\
                FGES& CPDAG &Java& \emph{TETRAD} & S \\
                RFCI& CPDAG &Java& \emph{TETRAD} & C \\
            \end{tabular}
        \end{center}
        \caption{Structure learning algorithms for Bayesian networks along with the space on which the algorithm operates (second column), the language the algorithms is implemented in (third column) the package where it is implemented (fourth column), and the type of algorithm (C)onstraint-based, (S)core-based or (H)ybrid (last column).}
        \label{tab:algorithms}
    \end{table}
\end{center}

\section{Installation and usage}\label{sec:installation}
\hmm{\emph{b}}\emph{enchpress}\xspace runs under Linux operating systems and requires that the following softwares, with corresponding version numbers, are installed: \hmm{\emph{s}}\emph{nakemake}\xspace \(\geq 5.2\),  \emph{Docker} \(\geq 20.0\),  and \emph{Singularity} \( \geq 3.5\) \citep{kurtzer2017singularity}, which is required by \hmm{\emph{s}}\emph{nakemake}\xspace to run \emph{Docker} containers.
Following \hmm{\emph{s}}\emph{nakemake}\xspace standards, configuration files should be placed in the \emph{config/} directory.

The easiest way to get started is to download or clone the \emph{git} repository of \hmm{\emph{b}}\emph{enchpress}\xspace from \url{https://github.com/felixleopoldo/benchpress} and use the example config file \emph{config/ex.json}.
\hmm{\emph{b}}\emph{enchpress}\xspace can then be run through standard \hmm{\emph{s}}\emph{nakemake}\xspace commands and we refer to \citep{10.1093/bioinformatics/bts480} for further details.
The steps are summarised as follows:
\begin{verbatim}
    $ git clone https://github.com/felixleopoldo/benchpress.git
    $ cd benchpress
    $ snakemake --cores 4 --use-singularity --configfile config/ex.json
\end{verbatim}
If \hmm{\emph{s}}\emph{nakemake}\xspace is installed with \emph{conda} \citep{anaconda} in an environment named \emph{e.g.} \emph{snakemake}, then it should first be activated by the following command:
\begin{verbatim}
    $ conda activate snakemake
\end{verbatim}
The configuration files for the case studies in Section~\ref{sec:case_study} have filename \emph{config/sec6.1.json} - \emph{config/sec6.4.json}.
However, these studies is rather large so laptop users are recommended to use a shorter \texttt{seed\_range} and only a subset of the algorithms.
\section{Simulation study} \label{sec:case_study}

The following subsections demonstrate four data scenarios to benchmark structure learning algorithms for Bayesian networks; where we consider binary and continuous data types.
For both types of data we consider two situations: one where we simulate both the graphical model and the data, and one where the graphical model is given and we only simulate data.
For performance evaluation, we demonstrate the use of the SHD, TPR, and the FPRp metrics in a ROC type curve on the pattern graph as described in Section~\ref{sec:evaluation}.
Each colored curve corresponds to a specific object (\texttt{id}) of an algorithm and the ROC shape comes from different values of a parameter specific to each algorithms, typically the parameter for a score or a test procedure.
Since not all methods support the same test procedures or scores, we have chosen the default settings for some of the algorithms.
The parameter values are selected so as to indicate how the performance changes.
However, some curves are shorter than would be desired due to range limits for the parameter in the programs that were used.

For binary data we use the \emph{Bayesian Dirichlet likelihood-equivalence} (BDe) score \citep{heckerman1995learning} in the score-based methods, while varying the equivalent sample size parameter \citep{silander2012sensitivity}.
The benchmarking exercise considers the following score-based algorithms:  HC (\texttt{hc-bde}), TABU (\texttt{tabu-bde}), ASOBS (\texttt{asobs-bdeu}), GOBNILP (\texttt{gobnilp-bde}), iterative search with \emph{maximum a posteriori} (MAP) estimation (\texttt{itsearch\_map-bde}), iterative search with Monte Carlo sampling estimation (\texttt{itsearch\_sample-bde}), and order MCMC with the end space of \texttt{itsearch\_sample-bde} as start space  (\texttt{omcmc\_itsample-bde}).
Furthermore, we include a number of constraint-based methods: MMHC (\texttt{mmhc-bde-mi}), Inter-IAMB (\texttt{interiamb-mi}), and GS (\texttt{gs-mi}) with the \emph{mutual information} test; the PC algorithm with the \(G^2\)-test (\texttt{pc-binCItest}); and the FCI (\texttt{fci-chi-square}), RFCI (\texttt{rfci-chi-square}), and  GFCI (\texttt{gfci-bdeu-chi-square}) with the \(\chi^2\)-test.
Since GFCI is a hybrid method it also uses the BDe score.
The varying parameter to draw ROC type curves for all constraint based methods is the significance level for the independence tests.

For the continuous datasets, we use two different scores in the score-based methods. The benchmarking exercise in this case includes Tabu (\texttt{tabu-bge}), GOBNILP (\texttt{gobnilp-bge}), HC (\texttt{hc-bge}), iterative search with MAP estimation (\texttt{itsearch\_map-bge}), and order MCMC using the end space of \texttt{itsearch\_sample-bge} as start space (\texttt{omcmc\_itsample-bge}) with the \emph{Bayesian Gaussian likelihood-equivalence} (BGe) score \citep{geiger1994learning,geiger2002parameter, kuipers2014addendum}, and varying levels of equivalent sample size.
Furthermore we include the FGES (\texttt{fges-sem-bic}) and GFCI (\texttt{gfci-sem-bic-fisher-z}) with the \emph{structural equation model Bayesian information criterion} (SEM-BIC) score \citep{glymour1986causal}, and varying levels of the penalty discount.
For constraint-based methods we include the PC algorithm (\texttt{pc-gaussCItest}), GFCI (\texttt{gfci-sem-bic-fisher-z}), FCI (\texttt{fci-fisher-z}), GFCI (\texttt{gfci-fisher-z}),  RFCI (\texttt{rfci-fisher-z}),  MMHC (\texttt{mmhc-bge-zf}), and Inter-IAMB (\texttt{interiamb-zf}) with the \emph{Fisher's z-transformation} test of partial correlations \citep{edwards2012introduction} with varying significance levels.

\subsection{Binary valued Bayesian network with random parameters and random DAG} \label{sec:randbin}
In this example we study a binary valued Bayesian network, where both the graph \(\graph\) and the parameters \(\Theta\) are regarded as random variables.
More specifically, we consider 50 models \(\{(\graph_i,\Theta_i)\}_{i=1}^{50}\), where each \(\graph_i\) is sampled according to the  Erdős–Rényi random DAG model using the \texttt{pcalg\_randdag}\, module (Section~\ref{sec:graph}), where the number of nodes is \(p=80\), the average number of neighbours (parents) per node is 4 (2) and the maximal number of parents per node is 5, see Listing~\ref{json:randomtridag}.
The parameters \(\Theta_i\) are sampled using the \texttt{bin\_bn}\, module (Section~\ref{sec:parameters}) and restricting the conditional probabilities within the range \([0.1, 0.9]\), see Listing~\ref{json:randombinbn}.
From each model, we draw two datasets \(\mathbf Y_i^{320},\) and \(\mathbf Y_i^{640}\) of sizes \(n=320\) and \(n=640\) using the \texttt{iid}\, module, see Listing~\ref{json:datasetup}.
The overall data setup is shown in Listing~\ref{json:bsetup}.

Figure~\ref{fig:binarybn} shows the ROC curves for the algorithms considered for the discrete data as described above, and Listing~\ref{json:roc} displays the \texttt{evaluation} section corresponding to this example.
We note that for the smaller number of observations (\(n=320\)), the results are more spread out in terms of FPRp than for the larger data sets (\(n=640\)).
For the smaller data sets (\(n=320\)), the algorithms standing out in terms of low SHD in combination with low best median FPRp \((< 0.12)\) and higher best median TPR \((>0.5)\), are FGES (\texttt{fges-bdeu}), GFCI (\texttt{gfci-bdeu-chi-square}) and iterative order MCMC (\texttt{omcmc\_itsample-bde}).
For the case with larger sample size the same algorithms still achieve higher performance, together with iterative search (\texttt{itsearch\_sample-bde})
\begin{listing}
    \begin{minted}[frame=single,
                   framesep=3mm,
                   linenos=true,
                   xleftmargin=21pt,
                   tabsize=4]{js}
{
    "id": "DAG80",
    "d": 4,
    "n": 80,
    "max_parents": 5,
    "par1": null,
    "par2": null,
    "method": "er"
}
    \end{minted}
\caption{Object in the \texttt{pcalg\_randdag}\, module for the example in Section~\ref{sec:randbin}.}
\label{json:randomtridag}
\end{listing}
\begin{listing}
    \begin{minted}[frame=single,
                   framesep=3mm,
                   linenos=true,
                   xleftmargin=21pt,
                   tabsize=4]{js}
{
    "id": "binarybn",
    "min": 0.1,
    "max": 0.9
}
    \end{minted}
\caption{Object in the \texttt{bin\_bn}\, module for the example in Section~\ref{sec:randbin}.}
\label{json:randombinbn}
\end{listing}

\begin{listing}
    \begin{minted}[frame=single,
                   framesep=3mm,
                   linenos=true,
                   xleftmargin=21pt,
                   tabsize=4]{js}
{
    "id": "iid",
    "standardized": false,
    "sample_sizes": [320, 640]
}
    \end{minted}
\caption{Object in the \texttt{data}\, module for the example in Section~\ref{sec:randbin}.}
\label{json:datasetup}
\end{listing}

\begin{listing}
    \begin{minted}[frame=single,
                   framesep=3mm,
                   linenos=true,
                   xleftmargin=21pt,
                   tabsize=4]{js}
{
    "graph_id": "DAG80",
    "parameters_id": "binarybn",
    "data_id": "iid",
    "seed_range": [1, 50]
}
    \end{minted}
\caption{Object in the \texttt{data} section in the \texttt{benchmark\_setup} section for the example in Section~\ref{sec:randbin}.}
\label{json:bsetup}
\end{listing}

\begin{listing}
    \begin{minted}[frame=single,
                   framesep=3mm,
                   linenos=true,
                   xleftmargin=21pt,
                   tabsize=4]{js}
{
    "roc":
    {
        "ids": ["asobs-bdeu", "gobnilp-bde", "tabu-bde",
                "mmhc-bde-mi", "hc-bde", "gs-mi", "interiamb-mi",
                "fci-chi-square", "rfci-chi-square",
                "gfci-bdeu-chi-square", "fges-bdeu",
                "itsearch_sample-bde", "pc-binCItest",
                "omcmc_itsample-bde"],
        "filename_prefix": "examples/section_6.1_",
        "point": true,
        "errorbar": true,
        "path": true,
        "text": false
    },
    "adjmat_true_plots": false,
    "graph_true_plots": false,
    "adjmat_plots": [],
    "graph_plots": [],
    "mcmc_heatmaps": [],
    "mcmc_traj_plots": [],
    "mcmc_autocorr": []
 }
 \end{minted}
\caption{The \texttt{evaluation} section for the example in Section~\ref{sec:randbin}.}
\label{json:roc}
\end{listing}

\begin{figure}[ht]

\begin{subfigure}[b]{.45\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{{figures/ROC_randbinarybnreps50.eps}}
  \caption{Binary valued Bayesian network with random parameters and random DAG.}
    \label{fig:binarybn}
\end{subfigure}
~
\begin{subfigure}[b]{.45\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{{figures/ROC_gaussreps50.eps}}
\caption{Linear Gaussian SEM with random weights and random DAG.}
  \label{fig:gaussbn}
\end{subfigure}

\begin{subfigure}[b]{.45\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{{figures/ROC_hepar2binreps50.eps}}
  \caption{Binary valued Bayesian network with random parameters and HEPAR II DAG.}
    \label{fig:hepar2bin}
\end{subfigure}
~
\begin{subfigure}[b]{.45\textwidth}
  \centering
  \includegraphics[width=1.0\linewidth]{{figures/ROC_hepar2gaussreps50.eps}}
\caption{Linear Gaussian SEM with random weights and HEPAR II DAG.}
\label{fig:hepar2gauss}
\end{subfigure}
\caption{ROC curves showing the FPRp and TPR evaluated on 50 datasets from 4 different Bayesian network models.
Each color corresponds to a structure learning algorithm with object \texttt{id} shown to the right.}
\label{fig:rocs}
\end{figure}

\subsection{Linear Gaussian SEM with random weights and random DAG}
\label{sec:randsem}
In this example we study a continuous random Bayesian network.
More specifically,
we consider 50  models \(\{(\graph_i,\Theta_i)\}_{i=1}^{50}\), where each graph \(\graph_i\) has \(p=80\) nodes and is sampled using the \texttt{pcalg\_randdag}\, module (Section~\ref{sec:graph}) as described in Section~\ref{sec:randbin}.
The parameters \(\Theta_i\) are sampled from the random linear Gaussian SEM using the \texttt{sem\_params}\, module (Section~\ref{sec:parameters}) with \(a=0.25, b=1, \mu=0\) and \(\sigma=1\).
We draw two standardised datasets \(\mathbf Y_i^{320}\) and \(\mathbf Y_i^{640}\), of sizes \(n=320\) and \(n=640\) (as for the discrete data scenario) from each of the models using the \texttt{iid}\, module.

Figure~\ref{fig:gaussbn} shows ROC results for all the algorithms considered for continuous data as described above.
For both sample sizes, the constraint based methods FCI (\texttt{fci-sem-bic-fisher-z}), RFCI (\texttt{rfci-fisher-z}), PC (\texttt{pc-gaussCItest}), and the score based GOBNILP (\texttt{gobnilp-bge}) have comparable and lower best median TPR (\(<0.7\)) then the remaining algorithms.
In terms of achieving high TPR (\(>0.7\)) iterative order MCMC (\texttt{omcmc\_itsample-bge}) and iterative search MCMC (\texttt{itsearch-bge}) stand out with near perfect performance, \emph{i.e.} SHD \(\approx 0\).
Among the other algorithms GFCI performs next best for both sample sizes with TPR \(\approx 0.85\) and FPRp \(\approx 0.08\), followed by HC (\texttt{hc-bge}), MMHC (\texttt{mmhc-bge-zf}), FGES (\texttt{fges-sem-bic}), and GS (\texttt{gs-zf}).
All algorithms completed except for GOBNILP (\texttt{gobnilp-bge}),  which had to be terminated by hand, even with maximal number of parents set to 2.

\subsection{Bayesian network with random parameters and fixed DAG}

In this example we study a random binary Bayesian network where the graph \(\graph\) is fixed and the associated parameters \(\Theta\) are regarded as random.
More specifically, we consider 50 models \(\{(\graph,\Theta_i)\}_{i=1}^{50}\), where
the graph structure \(\graph\) is that of the well known Bayesian network \emph{HEPAR II} (\vals{hepar2.csv}), first introduced in \cite{onisko2003probabilistic}.
This graph has 70 nodes and 123 edges and has become one of the standard benchmarks for evaluating the performance of structure learning algorithms.
The maximum number of parents per node is 6
and we sample the parameters \(\Theta_i\) using the \texttt{bin\_bn}\, module (Section~\ref{sec:parameters}), in the same manner as described in Section~\ref{sec:randbin}.
From each model we draw, as before, two datasets \(\mathbf Y_i^{320},\) and \(\mathbf Y_i^{640}\) of sample sizes \(n=320\) and \(n=640\), respectively, using the \texttt{iid}\, module.

Figure~\ref{fig:hepar2bin} shows the ROC curves for this scenario.
The algorithms appear to be divided in two groups with respect to their performance in terms of TPR. Constraint based methods including PC (\texttt{pc-binCItest}), FCI (\texttt{fci-chi-square}), RFCI (\texttt{rfci-chi-square}), inter-IAMB (\texttt{interiamb-mi}), MMHC (\texttt{mmhc-bde-mi}), and GS (\texttt{gs-mi}) appear to cluster in the lower scoring region (TPR \(<0.5\)), with GS achieving the lowest FPRp.
Score based methods on the other seem to concentrate in the higher scoring region (TPR \(>0.5\)).
For the smaller dataset (\(n=320\)), the group of better performing algorithms includes FGES (\texttt{fges-bdeu}), iteratve order MCMC (\texttt{omcmc\_itsample-bde}), iterative search MCMC (\texttt{itsearch\_sample-bde}), GFCI (\texttt{gfci-bdeu-chi-square}), all with FPRp \(\approx0.15\). For the larger data set (\(n=640\)) the higher performing group in terms of FPRp \(\approx0.11\) includes again iterative order MCMC (\texttt{omcmc\_itsample-bde}) and iterative search MCMC (\texttt{itsearch\_sample-bde}), this time joined by GOBNILP (\texttt{gobnilp-bde}).

\subsection{Linear Gaussian SEM with random weights and fixed DAG} \label{sec:fixeddagsem}
In this example we draw again 50  models \(\{(\graph,\Theta_i)\}_{i=1}^{50}\),
where \(\graph\) corresponds to the HEPARII network (\vals{hepar2.csv}), and \(\Theta_i\) are the parameters of a linear Gaussian SEM sampled using the \texttt{sem\_params}\, module (Section~\ref{sec:parameters}) with the same settings as in Section~\ref{sec:randsem}.
From each of the models \((\graph,\Theta_i)\), we draw once more independent data sets \(\mathbf Y_i^{320}\) and \(\mathbf Y_i^{640}\) of sizes \(n=320\) and \(n=640\)  (\texttt{iid}).

The ROC results of Figure \ref{fig:hepar2gauss} highlight that order MCMC (\texttt{omcmc\_itsample-bge}), iterative search MCMC (\texttt{itsearch\_sample\_bge}) visibly separate themselves from the rest in terms of SHD (combining both low FPRp and high TPR) for both sample sizes. In terms of best median TPR (\(>0.8\)) only HC (\texttt{hc-bge}) and Tabu (\texttt{tabu-bge}) display similar performaces, while performing considerably worse with respect to median FPRp (\(\approx0.1\) vs \(\approx0.17\)).
There is in general not a big difference in performance of the algorithms between the two sample sizes.
All algorithm completed except for GOBNILP (\texttt{gobnilp-bge}), which again had to be terminated by hand, even with maximal number of parents set to 2.

\section{Conclusions}
\hmm{\emph{b}}\emph{enchpress}\xspace provides a novel \emph{snakemake} workflow for scalable and re-producible execution and benchmarking of structure learning algorithms.
\emph{Snakemake}'s support for \emph{Docker} together with the simple data and graph representation enables \hmm{\emph{b}}\emph{enchpress}\xspace to compare algorithms implemented in different programming languages.
In its current version \hmm{\emph{b}}\emph{enchpress}\xspace already implements 19 of the state-of-the-art learning algorithms, as well as a number of modules for generating data models and evaluating performance.
As \hmm{\emph{b}}\emph{enchpress}\xspace is built in a completely modular form, \emph{snakemake} allows for seamless scaling of computations over multiple cores, grids or servers without any extra additional effort.
In addition, it is straightforward to integrate new modules into the software.
Even though the \hmm{\emph{b}}\emph{enchpress}\xspace project so far has a focus on structure learning for graphical models, we also see the potential in extending \hmm{\emph{b}}\emph{enchpress}\xspace to evaluate other statistical models.

\acks{We would like to thank all the developers and researchers that have made their software available in the public domain.}