11institutetext: Universidad Nacional de Rosario and CIFASIS, Rosario,
Argentina 22institutetext: ENSIIE, lab. Samovar, Évry-Courcouronnes, France

# Comparing EventB, $\\{log\\}$ and Why3 Models of Sparse Sets

Maximiliano Cristiá 11 Catherine Dubois 22

###### Abstract

Many representations for sets are available in programming languages
libraries. The paper focuses on sparse sets used, e.g., in some constraint
solvers for representing integer variable domains which are finite sets of
values, as an alternative to range sequence. We propose in this paper verified
implementations of sparse sets, in three deductive formal verification tools,
namely EventB, $\\{log\\}$ and Why3. Furthermore, we draw some comparisons
regarding specifications and proofs.

##  1 Introduction

Sets are widely used in programs. They are sometimes first-class objects of
programming languages, e.g. SETL [23] or $\\{log\\}$ [11], but more frequently
they are data structures provided in libraries. Many different representations
are available, depending on the targeted set operations. In this paper, we
focus on sparse sets, introduced by Briggs and Torczon in [3], used in
different contexts and freely available for different programming languages
(Rust, C++ and many others). In particular, sparse sets are used in constraint
solvers as an alternative to range sequences or bit vectors for implementing
domains of integer variables [18] which are nothing else than mathematical
finite sets of integers. Their use in solvers implementations is motivated by
-at least- the two following properties: searching and removing an element are
constant-time operations—removing requires only two swapping operations on
arrays; sparse sets are cheap to trail and restore, which is a key point when
backtracking.

Confidence on constraint solvers using sparse sets can be improved if the
algorithms implementing the main operations are formally verified, as it has
been done by Ledein and Dubois in [19] for the traditional implementation of
domains as range sequences. Hence, the main contribution of this paper is a
verified implementation of sparse sets for representing finite sets of
integers in EventB, $\\{log\\}$ and Why3. We prove that the implemented
operations preserve the invariants and we also prove properties that can be
seen as formal foundations of trailing and restoring. As far as we know, this
is the first formally verified implementation of sparse sets, whereas it has
been done for other representations e.g. [16, 19]. All the specifications and
proofs can be found here: https://gitlab.com/cdubois/sets2023.git.

It has been known for decades that there is no silver bullet for software
engineering or software development. The best we can do as software engineers
is to increase our toolbox as much as possible and use the best available tool
in it for the problem at hand. This software engineer practical principle
still applies when it comes to formal development, formal methods and formal
verification. In our opinion the Formal Methods (FM for short) community
should have as much information as possible about the relative advantages and
disadvantages of different FM methods and tools. With the intention to shed
some light on the ups and downs of different FM, we specified and verified
sparse sets with three different FM techniques. Then, a second contribution of
this paper is a comparison of these FM w.r.t. aspects such as expressiveness,
specification analysis and automated proof.

##  2 Sparse sets

We deal here with sets as subsets of natural numbers up to $N-1$, where $N$ is
any non null natural number. A sparse set $S$ is represented by two arrays of
length $N$ called $mapD$ and $domD$ (as in [18]), and a natural number
$sizeD$. The array $mapD$ maps any value $v\in[0,N-1]$ to its index $ind_{v}$
in $domD$, the value indexed by $ind_{v}$ in $domD$ is $v$. The main idea that
brings efficiency when removing an element or testing membership is to split
$domD$ into two sub-arrays, $domD[0,sizeD-1]$ and $domD[sizeD,N-1]$,
containing resp. the elements of $S$ and the elements of $[0,N-1]$ not in $S$.
Then, if $S$ is empty, $sizeD$ is equal to 0, if $S$ is the full set, then
$sizeD$ is $N$. Checking if an element $i$ belongs to the sparse set $S$
simply consists in the evaluation of the expression $mapD[i]<sizeD$. Removing
an element from the set consists in moving this element to $domD[sizeD,N-1]$
(with 2 swaps in $mapD$ and $domD$ and decreasing $sizeD$). Binding $S$ to the
singleton set $\\{v\\}$ follows the same idea: moving this element at the
first place in $domD$ and assigning the value 1 to $sizeD$.

In our formalizations, we only deal with two operations consisting in removing
an element in a sparse set and bind a sparse set to a singleton set since
these two operations are fundamental when solving constraints. In this
context, we may also need to walk through all the elements of a variable
domain, it means exploring $domD[0..sizeD-1]$. If minimal and maximal values
are required, then they have to be maintained in parallel. This is outside the
scope of this work.

##  3 EventB formal development

In this section we succinctly introduce the EventB formal specification
language and with more detail the EventB models for sparse sets.

###  3.1 EventB

EventB [1] is a deductive formal method based on set theory and first order
logic allowing users to design correct-by-construction systems. It relies on a
state-based modeling language in which a model, called a machine, is made of a
state and a collection of events allowing for state changes. The state
consists of variables constrained by invariants. Proof obligations are
generated to verify the preservation of invariants by events. A machine may
use a -mathematical- context which introduces abstract sets, constants, axioms
or theorems. A formal design in EventB starts with an abstract machine which
is usually refined several times. Proof obligations are generated to verify
the correctness of a refinement step.

An event may have parameters. When its guards are satisfied, its actions, if
any, are executed, updating state variables. Actions may be -multiple-
deterministic assignments, $x,y:=e,f$, or -multiple- nondeterministic ones,
$x,y:\mid BAP(x,x^{\prime},y,y^{\prime})$ where $BAP$ is called a Before-After
Predicate relating current ($x$, $y$) and next ($x^{\prime}$, $y^{\prime}$)
values of state variables $x$ and $y$. In the latter case, $x$ and $y$ are
assigned arbitrary values satisfying the BAP predicate. When using such a non-
deterministic form of assignment, a feasibility proof obligation is generated
in order to check that there exist values for $x^{\prime}$ and $y^{\prime}$
such that $BAP(x,x^{\prime},y,y^{\prime})$ holds when the invariants and
guards hold. Furthermore when this kind of action is used and refined, the
concrete action updating $x$ and $y$ is required to assign them values which
satisfy the BAP predicate.

In the following, we use Rodin, an Eclipse based IDE for EventB project
management, model edition, refinement and proof, automatic proof obligations
generation, model animation and code generation. Rodin supports automatic and
interactive provers [14]. In this work we used the standard provers (AtelierB
provers) and also the SMT solvers VeriT, CVC3 and CVC4. More details about
EventB and Rodin can be found in [1] and [2].

###  3.2 EventB formalization

The formalization is made of six components, i.e. two contexts, a machine and
three refinements. Context $Ctx$ introduces the bound $N$ as a non-zero
natural number and context $Ctx1$ extends the latter with helper theorems. The
high level machine gives the abstract specification. This model contains a
state composed of a finite set $D$, constrained to be a subset of the
(integer) range $0..N-1$, and two events, to remove an element from $D$ or set
$D$ as a singleton set (see Fig. 1 in which $bind$ is removed for lack of
space).

MACHINE Domain
SEES Ctx
VARIABLES D
INVARIANTS inv1: $D\subseteq 0\mathbin{\ldotp\ldotp}N-1$
EVENTS
Initialisation begin act1: $D:=0\mathbin{\ldotp\ldotp}{}N-1$ end
Event remove $\langle$ordinary$\rangle$$\mathrel{\widehat{=}}$
any v where grd1: $v\in{}D$ then act1: $D:=D\setminus\\{v\\}$ end
END

Figure 1: EventB abstract specification, the Domain machine

The first refinement (see Fig.2) introduces the representation of the domain
as a sparse set, i.e. two arrays $mapD$ and $domD$ modeled as total functions
and also the variable $sizeD$ which is a natural number in the range $0..N$.
Invariants $inv4$ and $inv5$ constrain $mapD$ and $domD$ to be inverse
functions of each other. The gluing invariant $inv6$ relates the states
between the concrete and former abstract machines. So the set
$domD[0..sizeD-1]$ containing the elements of the subarray from 0 to $sizeD-1$
is exactly the set $D$.

Theorem $inv7$ is introduced to ease some interactive proofs, it is proved as
a consequence of the previous formulas ($inv1$ to $inv6$). It follows directly
from a theorem of $Ctx1$ whose statement is $inv7$ where $domD$ and $mapD$ are
universally quantified. Theorem $inv8$, also used in an interactive proof, and
automatically proved by CVC3, states that $domD$ is an injective function.

Variables $mapD$ and $domD$ are both set initially to the identity function on
$0..N-1$ and $sizeD$ to $N$. So invariants are satisfied at the initial state.
Machine SparseSets  ref1 refines the events of the initial machine by non
deterministic events. So here the $remove$ event assigns the three state
variables with values that satisfy invariants and also such that $sizeD$
strictly decreases and removed elements in $domD$ are kept at the same place
(properties in bold font). Event $bind$ follows the same pattern (again not
shown here).

MACHINE SparseSets  ref1
REFINES Domain
SEES Ctx1
VARIABLES domD mapD sizeD
INVARIANTS
inv1: $domD\in 0\mathbin{\ldotp\ldotp}N-1\rightarrow
0\mathbin{\ldotp\ldotp}N-1$
inv2: $mapD\in 0\mathbin{\ldotp\ldotp}N-1\rightarrow
0\mathbin{\ldotp\ldotp}N-1$
inv3: $sizeD\in 0\mathbin{\ldotp\ldotp}N$
inv4: $domD\mathchar
24635\relax\;mapD=\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1}$
inv5: $mapD\mathchar
24635\relax\;domD=\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1}$
inv6: $domD[0\mathbin{\ldotp\ldotp}sizeD-1]=D$
inv7: $\langle$ theorem $\rangle$
$\mathop{\mathstrut{\forall}}\nolimits x,v\cdot x\in
0\mathbin{\ldotp\ldotp}N-1\mathrel{\wedge}v\in
0\mathbin{\ldotp\ldotp}N-1\Rightarrow(mapD(v)=x\Leftrightarrow domD(x)=v)$
inv8: $\langle$ theorem $\rangle$ $domD\in
0\mathbin{\ldotp\ldotp}N-1\rightarrowtail 0\mathbin{\ldotp\ldotp}N-1$
EVENTS
Initialisation
act1:
$mapD,domD:=\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1},\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1}$
act2: $sizeD:=N$
Event remove $\langle$ordinary$\rangle$$\mathrel{\widehat{=}}$ refines remove
any v
where grd1: $v\in 0\mathbin{\ldotp\ldotp}N-1\mathrel{\wedge}{}$ grd2:
$0<sizeD\mathrel{\wedge}{}$ grd3: $mapD(v)<sizeD$
then act1: $mapD,domD,sizeD:\mid$
$(domD^{\prime}\in 0\mathbin{\ldotp\ldotp}N-1\rightarrow
0\mathbin{\ldotp\ldotp}N-1\mathrel{\wedge}mapD^{\prime}\in
0\mathbin{\ldotp\ldotp}N-1\rightarrow 0\mathbin{\ldotp\ldotp}N-1$
${}\mathrel{\wedge}domD^{\prime}\mathchar
24635\relax\;mapD^{\prime}=\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1}\mathrel{\wedge}mapD^{\prime}\mathchar
24635\relax\;domD^{\prime}=\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1}$
${}\mathrel{\wedge}domD^{\prime}[0\mathbin{\ldotp\ldotp}sizeD^{\prime}-1]=domD[0\mathbin{\ldotp\ldotp}sizeD-1]\setminus\\{v\\}\mathrel{\wedge}{\bf
sizeD^{\prime}<sizeD}$
${}\mathrel{\wedge}{\bf(sizeD\mathbin{\ldotp\ldotp}N-1)\dres
domD^{\prime}=(sizeD\mathbin{\ldotp\ldotp}N-1)\dres DomD}$
end

Figure 2: EventB first refinement

The second refinement has the same state than the previous refinement (see
Fig. 3). Its events implement the operations using the new state variables. It
is a straightforward translation of the algorithms described in [18].

The only reason to have introduced the intermediate model SparseSets  ref1 is
to express the properties written in bold font and thus generate, in the next
refinement, proof obligations which, when discharged, will not only ensure
that the events refined in Fig. 3 preserve the invariants $inv1$, $inv2$
…$inv6$ but also the local properties regarding $sizeD$ and $domD[sizeD..N-1]$
(SIM proof obligations).

MACHINE SparseSets  ref2
REFINES SparseSets  ref1
SEES Ctx1
VARIABLES domD mapD sizeD
EVENTS
Initialisation
act1:
$mapD,domD:=\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1},\mathop{\mathrm{id}}_{0\mathbin{\ldotp\ldotp}N-1}$
act2: $sizeD:=N$
Event remove $\langle$ordinary$\rangle$$\mathrel{\widehat{=}}$ refines remove
any v
where grd1: $v\in 0\mathbin{\ldotp\ldotp}N-1\mathrel{\wedge}{}$ grd2:
$0<sizeD\mathrel{\wedge}{}$ grd3: $mapD(v)<sizeD$
then
act1: $domD:=domD\mathbin{{\dres}\hbox
to0.0pt{\hss+\\!\\!}\,}\\{mapD(v)\mapsto domD(sizeD-1),sizeD-1\mapsto v\\}$
act2: $mapD:=mapD\mathbin{{\dres}\hbox to0.0pt{\hss+\\!\\!}\,}\\{v\mapsto
sizeD-1,domD(sizeD-1)\mapsto mapD(v)\\}$
act3: $sizeD:=sizeD-1$
end

Figure 3: EventB second refinement

The feasibility (FIS) proof obligations generated by the non-deterministic
events of SparseSets  ref1 require to prove that there exist values such that
the BAP predicate holds. We can prove it using the new values of $domD$,
$mapD$ and $sizeD$ specified in the last refinement as witnesses. The
simulation (SIM) proof obligations generated by events of SparseSets  ref2
require to prove that the latter values again satisfy the BAP predicate used
in SparseSets  ref1. In order not to do these -interactive- proofs twice, we
generalize them and prove them as theorems of the context. Thus to discharge
the FIS and SIM proof obligations, we only have to instanciate these theorems
to provide a proof.

A last algorithmic refinement, omitted here, refines the $remove$ event in two
events, $removeLastButOne$ and $removeLast$. The former differs from $remove$
only by its more restrictive guard; the latter is dedicated to the case where
the element with index $sizeD-1$ in $domD$ is removed thus avoiding the
unnecessary swapping.

##  4 $\\{log\\}$ formal development

In this section we briefly present the $\\{log\\}$ tool and how we used it to
encode the EventB model of sparse sets.

###  4.1 $\\{log\\}$

$\\{log\\}$ is a constraint logic programming (CLP) language and
satisfiability solver where sets and binary relations are first-class citizens
[21, 15, 6]. The tool implements several decision procedures for expressive
fragments of set theory and set relation algebra including cardinality
constraints [13], restricted universal quantifiers [12], set-builder notation
[8] and integer intervals [10]. In previous works $\\{log\\}$ has been
satisfactory tested against some known case studies [7, 9, 5].

$\\{log\\}$ code enjoys the _formula-program duality_. This means that
$\\{log\\}$ code can behave as both a formula and a program. When seen as a
formula, it can be used as a specification on which verification conditions
can be (sometimes automatically) proved. When seen as a program, it can be
used as a (less efficient) regular program. Due to the formula-program
duality, a piece of $\\{log\\}$ code is sometimes called _forgram_ —a
portmanteau word resulting from combining _for_ mula with prog _gram_.

###  4.2 $\\{log\\}$ formalization

The $\\{log\\}$ formalization presented in this paper is the result of
translating the EventB abstract specification (i.e., Fig. 1) and the second
refinement (i.e. Fig. 3). Both EventB models can be easily translated into
$\\{log\\}$ by using the (still under development) state machine specification
language (SMSL) defined on top of $\\{log\\}$ (see Fig. 4 and 5) [22]. The
notions of context and refinement are not available in SMSL. For this reason,
refinements introduced in the EventB model have to be manually encoded in
$\\{log\\}$. The context is encoded simply as an axiom. In order to ensure
that the $\\{log\\}$ code verifies the properties highlighted in bold in Fig.
2 as well as the gluing invariant (i.e., $inv6$), a few user-defined
verification conditions are introduced as theorems. Since the first EventB
refinement is introduced to express the properties written in bold, its events
have not been encoded in $\\{log\\}$.

    parameters([N,I]).
    variables([D,DomD,MapD,SizeD]).

    axiom(axm1).
    axm1(N) :- 1 =< N.

    axiom(axm2).
    axm2(N,I) :- M is N - 1 & id(int(0,M),I).

    invariant(inv11).
    inv11(DomD) :- pfun(DomD).

    n_inv11(DomD) :- neg(  pfun(DomD)  ).

    invariant(inv12).
    inv12(N,DomD) :- N1 is N - 1 & dom(DomD,int(0,N1)).

    invariant(inv13).
    inv13(N,DomD) :- N1 is N - 1 & ran(DomD,R) & subset(R,int(0,N1)).

    invariant(inv4).
    inv4(N,I,DomD,MapD) :- axm2(N,I) & comppf(DomD,MapD,I).

    inv6(D,DomD,SizeD) :-
      S is SizeD - 1 &
      foreach([X,Y] in DomD, X in int(0,S) implies Y in D) &
      foreach(X in D, exists([A,B] in DomD, A in int(0,S) & B = X)).

    inv7(MapD,DomD) :-
      foreach([[V,Y1] in MapD, [X,Y2] in DomD],
        (Y1 = X implies Y2 = V) & (Y2 = V implies Y1 = X)  ).

    theorem(inv7_th).
    inv7_th(N,MapD,DomD) :-
      neg(inv4(N,I,DomD,MapD) & inv5(N,I,DomD,MapD) implies inv7(MapD,DomD)).

Figure 4: Some representative axioms and invariants of the $\\{log\\}$ forgram

Figures 4 and 5 list only representative parts of the $\\{log\\}$ forgram. We
tried to use the same identifiers as for the EventB models as much as
possible. In this way, for example, the invariant labeled as $inv6$ in the
SparseSets  ref1 machine (Fig. 2), is named `inv6` in the $\\{log\\}$ forgram.
The name of variables in $\\{log\\}$ cannot fully complain with those used in
the EventB models because $\\{log\\}$ requires all variables to begin with a
capital letter. So, for example, $domD$ in the SparseSets  ref1 machine
becomes `DomD` in $\\{log\\}$.

As can be seen in Fig. 4, the state machine specification language defined on
top of $\\{log\\}$ allows for the declaration of parameters (similar to EventB
context constants), state variables, axioms (similar to EventB context axioms)
and invariants. Parameter `I` is used to compute the identity relation on the
integer interval $[0,N-1]$ as shown in axiom `axm2`, which in turn is used in
invariant `inv4`. As $\\{log\\}$ is a CLP language implemented on top of
Prolog, it inherits many of Prolog’s features. In particular, integer
expressions are evaluated by means of the `is` predicate. Along the same
lines, all set operators are implemented in $\\{log\\}$ as constraints. For
example, `id(A,R)` is true when `R` is the identity relation on the set `A`.
The term `int(0,M)` corresponds to the integer interval $[0,M]$.

Invariants named `inv11`, `inv12` and `inv13` correspond to invariant $inv1$
of the SparseSets  ref1 machine. Splitting invariants in smaller pieces, is a
good practice when using $\\{log\\}$ as a prover because it increases the
chances of automated proofs. `n_inv11` implements the negation of invariant
`inv11`. $\\{log\\}$ does not automatically compute the negation of user-
defined predicates. As a user-defined predicate can contain existential
variables, its negation could involve introducing universal quantifiers which
fall outside $\\{log\\}$’s decision procedures. Then, users are responsible
for ensuring that all predicates are safe.

In invariant `inv6` we can see the `foreach` constraint. This constraint
implements the notion of _restricted universal quantifier_ (RUQ). That is, for
some $\\{log\\}$ formula $\phi$ and set `A`, foreach(X in A, $\phi($X$)$)
corresponds to $\mathop{\mathstrut{\forall}}\nolimits X.(X\in
A\Rightarrow\phi(X))$. In a `foreach` constraint it is possible to quantify
over binary relations, as is the case of `inv6`. Hence, we have a quantified
ordered pair (`[X,Y]`), rather than just a variable. Likewise, $\\{log\\}$
offers the `exists` constraint implementing the notion of _restricted
existential quantifier_ (REQ). The important point about REQ and RUQ is not
only their expressiveness but the fact that there is a decision procedure
involving them [12]. In `inv6` these constraints are used to state a double
set inclusion equivalent to the EventB formula $domD[0..sizeD-1]=D$. If the
user is not convinced or unsure about the validity of this equivalence (s)he
can use $\\{log\\}$ itself to prove it.

Note that `inv7` is not declared as an invariant because in Fig. 2 it is a
theorem that can be deduced from previous invariants. Therefore, we introduce
it as a simple predicate but then we declare a theorem whose conclusion is
`inv7`. Later, $\\{log\\}$ will include `inv7_th` as a proof obligation and
will attempt to discharge it. Given that $\\{log\\}$ is a satisfiability
solver, if $\Phi$ is intended to be a theorem then we ask it to prove the
unsatisfiability of $\neg\;\Phi$.

    operation(remove).
    remove(N,SizeD,MapD,DomD,V,SizeD_,MapD_,DomD_) :-
      M is N - 1 & V in int(0,M) & 0 < SizeD & S is SizeD - 1 &
      MapD = {[V,Y1],[Y2,Y4] / MapD1} & disj({[V,Y1],[Y2,Y4]},MapD1) &
      Y1 < SizeD  &
      DomD = {[S,Y2],[Y1,Y5] / DomD1} & disj({[S,Y2],[Y1,Y5]},DomD1) &
      DomD_ = {[S,V],[Y1,Y2] / DomD1} &
      MapD_ = {[V,S],[Y2,Y1] / MapD1} &
      SizeD_ = S.
    theorem(remove_pi_inv6).
    remove_pi_invr6(N,SizeD,MapD,DomD,V,SizeD_,MapD_,DomD_) :-
      inv7(MapD,DomD) &
      neg(    inv6(D,DomD,SizeD) &
              remove(V,D,D_) &
              remove(N,SizeD,MapD,DomD,V,SizeD_,MapD_,DomD_)
                  implies inv6(D_,DomD_,SizeD_)    ).
    theorem(remove_b2).
    remove_b2(N,SizeD,MapD,DomD,V,SizeD_,MapD_,DomD_) :-
      neg(    N1 is N - 1 &
              remove(N,SizeD,MapD,DomD,V,SizeD_,MapD_,DomD_) &
              fimg(int(SizeD,N1),DomD_,D)
                  implies fimg(int(SizeD,N1),DomD,D)  ).

Figure 5: The remove operation and some user-defined proof obligations

Moving into in Fig. 5 we can see the encoding of the $remove$ operation
specified in the SparseSets  ref2 machine of Fig. 3, along with two user-
defined proof obligations. In $\\{log\\}$, there is no global state so state
variables have to be included as explicit arguments of clauses representing
operations. Next-state variables are denoted by decorating the base name with
an underscore character (e.g., `SizeD_` corresponds to the value of `SizeD` in
the next state). Another important difference between the EventB and the
$\\{log\\}$ specifications is that in the latter we can use _set unification_
to implement function application. For instance, `DomD = {[S,Y2],[Y1,Y5] /
DomD1}` is equivalent to the EventB predicate:
$\mathop{\mathstrut{\exists}}\nolimits
y_{2},y_{5},domD_{1}.(domD=\\{sizeD-1\mapsto y_{2},y_{1}\mapsto y_{5}\\}\cup
domD_{1})$, where $y_{1}=mapD(v)$ (due to the previous set unification). The
not-membership constraints following the equality constraint prevent
$\\{log\\}$ to generate repeated solutions. Hence, when `remove` is called
with some set term in its fourth argument, this term is unified with
`{[S,Y2],[Y1,Y5] / DomD1}`. If the unification succeeds, then the images of
`S` and `Y1` are available.

As said before, some user-defined proof obligations are introduced as theorems
to ensure that the $\\{log\\}$ forgram verifies the gluing invariant (i.e.,
$inv6$) and the properties written in bold in machine SparseSets  ref1.
Precisely, theorem `remove_pi_inv6` states that if `inv6` holds and `remove`
and its abstract version (not shown in the paper) are executed, then `inv6`
holds in the next state.111remove and its abstract version can be
distinguished by their arities.

Likewise, theorem `remove_b2` ensures that the second property written in bold
in machine SparseSets  ref1 is indeed a property of the $\\{log\\}$ forgram.
As can be seen, the theorem states that if `remove` is executed and the
functional image222fimg is a user-defined $\\{log\\}$ predicate computing the
relational image through a function—fimg stands for _functional image_. of the
interval from `SizeD` up to `N-1` through `DomD_` is `D`, then it must
coincide with the functional image of the same interval but through `DomD`.

Once the specification is ready, we can call the verification condition
generator (VCG) and run the verification conditions (VC) so generated:

    {log}=> vcg(’sp.pl’) & consult(’sp-vc.pl’) & check_vcs_sp.

VCs include the satisfiability of the conjunction of all axioms, the
satisfiability of each operation and preservation lemmas for each and every
operation and invariant. The last command above will attempt to automatically
discharge every VC. Part of the output is as follows:

    Checking remove_is_sat ... OK
    Checking remove_pi_inv11 ... ERROR

An `ERROR` answer means that, for some reason, $\\{log\\}$ is unable to
discharge the VC. Most of the times this is due to some missing hypothesis
which, in turn, is due to the way the VCG generates the VCs. Briefly, when it
comes to invariance lemmas, the VCG generates them with the minimum number of
hypothesis. So, for instance, the invariance lemma named `remove_pi_inv11` is
as follows:

    neg(  inv11(DomD) &
          remove(N,SizeD,MapD,DomD,V,SizeD_,MapD_,DomD_) implies
          inv11(DomD_)  ).

By including minimum hypothesis, $\\{log\\}$ will have to solve a simpler goal
which reduces the possibilities to have a complexity explosion. If the
hypothesis is not enough, the `findh` command can be used to find potential
missing hypothesis. In this way, users can edit the VC file, add the missing
hypothesis and run the VC again. If more hypotheses are still missing, the
process can be executed until the proof is done—or the complexity explosion
cannot be avoided.

$\\{log\\}$ discharges all the VC generated by the VCG for the present
forgram.

##  5 Why3 formal development

In this section we briefly introduce the Why3 platform and describe with some
details our specification of sparse sets.

###  5.1 Why3

Why3 [17] is a platform for deductive program verification providing a
language for specification and programming, called WhyML, and relies on
external automated and interactive theorem provers, to discharge verification
conditions. In the context of this paper, we used Why3 with the SMT provers
CVC4 and Z3.

Proof tactics are also provided, making Why3 a proof environment close to the
one of Rodin for interactive proofs. Why3 supports modular verification.

WhyML allows the user to write functional or imperative programs featuring
polymorphism, algebraic data types, pattern-matching, exceptions, references,
arrays, etc. These programs can be annotated by contracts and assertions and
thus verified. User-defined types with invariants can be introduced, the
invariants are verified at the function call boundaries. Furthermore to
prevent logical inconsistencies, Why3 generates a verification condition to
show the existence of at least one value satisfying the invariant. To help the
verification, a witness is explicitly given by the user (see the by clause in
Fig. 6). The old and at operators can be used inside post-conditions and
assertions to refer to the value of a mutable program variable at some past
moment of execution. In particular old t in a function post-condition refers
to the value of term t when the function is called. Programs may also contain
ghost variables and ghost code to facilitate specification and verification.
From verified WhyML programs, correct-by-construction OCaml programs (and
recently C programs) can be automatically extracted.

###  5.2 Why3 formalization

From the Why3 library, we use pre-defined theories for integer arithmetic,
polymorphic finite sets and arrays. In the latter, we use in particular the
swap operation that exchanges two elements in an array and its specification
using the exchange predicate.

We first define a record type, sparse, whose mutable fields are a record of
type sparse  data containing the computational elements of a sparse set
representation and a ghost finite set of integer numbers which is the abstract
model of the data structure. The type invariant of sparse relates the abstract
model with the concrete representation. It is used to enforce consistency
between them. Invariants enforcing consistency between the two arrays mapD and
domD and the bound sizeD are attached to the sparse  data type: lengths of the
arrays is n, contents are belonging to $0..\mbox{{n}}-1$ and the two arrays
are inverse of each other, sized is in the interval $0..\mbox{{n}}$. These
type definitions and related predicates are shown in Fig. 6.

⬇

predicate dom_ran (a: array int) (n: int) =

0 <= n && a.length = n && forall i. 0<=i<n -> 0<=a[i]< n

predicate comp_is_id (a: array int) (b: array int) (n: int) =

forall v,i. (0<=i<n && 0<=v<n) -> (a[i]=v <-> b[v]=i)

type sparse_data = {n: int; mutable domD: array int;

mutable mapD: array int; mutable sizeD: int; }

invariant {dom_ran domD n && dom_ran mapD n &&

comp_is_id domD mapD n && 0 <= sizeD <= n }

by {n = 0 ; domD = make 0 0 ; mapD = make 0 0; sizeD=0}

type sparse = {mutable data: sparse_data; mutable ghost setD: fset int;}

invariant {subset setD (interval 0 data.n) &&

forall i:int. 0 <= i < data.n ->

(i < data.sizeD <-> mem data.domD[i] setD)}

by {data = {n = 0 ; domD = make 0 0 ; mapD = make 0 0; sizeD=0} ;

setD=FsetInt.empty }

Figure 6: WhyML types for sparse sets

Our formalization (see Fig. 7, where, again, $bind$ is removed for lack of
place) contains three functions, swap  sparse  data, remove  sparse and bind
sparse, which update their arguments. They are the straightforward translation
of the algorithms in [18] in WhyML, except for the supplementary ghost code
(the last statement in both remove  sparse and bind  sparse) which updates the
abstract model contained in a.setD. Function swap  sparse  data is a helper
function which is called in the other ones. The contract of swap  sparse  data
makes explicit the modifications of both arrays a.mapD and a.domD, using the
exchange predicate defined in the library. Verification conditions for this
function concern the conformance of the code to the two post-conditions
(trivial as it is ensured by swap) and also the preservation of the invariant
attached to the sparse  data type—i.e. mainly that a.mapD and a.DomD after
swapping elements remain inverse from each other. Both remove  sparse and bind
sparse act not only on the two arrays and the bound but also on the ghost
part, i.e. the corresponding mathematical set a.setD. Thus the verification
conditions here not only concern the structural invariants related to mapD,
domD and sizeD but also the ones deriving from the use of the sparse type,
proving the link between the abstract logical view (using finite sets) and the
computational one implemented through arrays.

⬇

predicate same_end (a : array int) (b : array int) (s : int) (n : int) =

forall i. s <= i < n -> a[i] = b[i]

let swap_sparse_data (a : sparse_data) (i : int) (j : int)

requires {0<=i<a.n}

requires {0<=j<a.n}

ensures {exchange (old a.domD) a.domD i j}

ensures {exchange (old a.mapD) a.mapD a.domD[i] a.domD[j]} =

swap a.domD i j;

a.mapD[a.domD[i]] <- i;

a.mapD[a.domD[j]] <- j;

let remove_sparse (v : int) (a : sparse)

requires {0<=v<a.data.n}

requires {a.data.mapD[v] < a.data.sizeD}

requires {a.data.sizeD > 0}

ensures {old a.data.sizeD > a.data.sizeD}

ensures {same_end a.data.domD (old a.data.domD) (old a.data.sizeD) a.data.n} =

swap_sparse_data a.data a.data.mapD[v] (a.data.sizeD - 1);

a.data.sizeD <- a.data.sizeD - 1;

a.setD <- remove v a.setD

Figure 7: WhyML functions for sparse sets

Observe that types sparse  data and sparse correspond to the state and
invariants of the EventB refinements. The abstract specification presented in
the first EventB machine becomes a ghost field in WhyML. The invariant of the
sparse type corresponds to the EventB gluing invariant ($inv6$). A similar
transposition happens for the operations. Actions in the EventB abstract
events, i.e. updating the abstract set, appear as ghost code in WhyML.

All proofs are discovered by the automatic provers except for some proof
obligations related to the remove function. Nevertheless these proofs are
simplified thanks to some Why3 tactics that inject some hints that can be used
by the external provers to finish the proofs.

##  6 Comparison and discussion

Set theory is primitive in EventB and $\\{log\\}$ whereas Why3 which permits
to express other theories, provides a theory for it. Rodin uses provers where
set theory is primitive but can also call external provers such as VeriT, Z3
and CVC4—where set theory is not primitive. However a big effort has been done
to process set theory in VeriT, which is often recognized as allowing
significant improvements in proofs [20]. Why3 relies entirely on external
provers where set theory is not primitive. Conversely, $\\{log\\}$ is a
satisfiability solver that can only work with set theory—and linear integer
algebra. It is the only of the three tools implementing advanced decision
procedures for set theory. Likely, this proved to be crucial for $\\{log\\}$
being able to be the only tool that automatically discharged all the VC,
although it required a simple hypothesis discovery procedure. It should be a
concern the time $\\{log\\}$ needs to discharge all the VC because with more
complex models the resolution time might be prohibitive. It worth to be
studied ways of avoiding the algorithmic complexity of the decision procedures
implemented in $\\{log\\}$. Results on Computable Set Theory should be
revisited (eg. [4]). Why3 and Rodin interactive proofs are not numerous and
remain quite simple.

In EventB, 51 proof obligations were generated for the whole development,
around half of them coming from the first refinement. 37 were proven
automatically by the standard provers (AtelierB provers), 18 automatically by
SMT provers, mainly VeriT, either directly or after applying the Rodin lasso
allowing for adding additional, backup hypotheses having identifiers in common
with the goal. Only two proof obligations required real human intervention,
mainly instantiations of the general theorems introduced in $Ctx1$ or explicit
witnesses introduction in the case of feasibility proof obligations.

After working in the way described in Sect. 4, $\\{log\\}$ discharges all the
38 VC generated by the VCG in around 7 minutes.

Why3 makes it possible to apply transformations (e.g. split conjunctions) on a
proof goal instead of calling an automatic prover on it. Some of these
transformations are very simple, e.g. splitting conjunctions, and can then
been applied systematically and automatically. Most of the generated VC in our
formalization were proven automatically thanks to the split transformation.
Only two of them about pieces of type invariants, required human interaction
to insert some more complex transformations, e.g a case analysis on indexes in
mapD (case (i=a  data.mapD[v]). At the end, 55 VC were proved by CVC4, except
two of them discharged by Z3, in a total amount of time of 30 seconds.

Clearly, all three tools are expressive enough for the problem at hand.
However, the EventB specification is probably the most readable. The tools
permit to express axioms, invariants and automatically generate similar VC.
$\\{log\\}$ still needs work to express how two models are linked in terms of
abstraction/refinement relations. Writing some key properties proved to be
complex in EventB. Indeed, it was necessary to add a somewhat artificial
refinement level for Rodin being able to generate the desired VC linking.
These properties can be easily defined by the user in $\\{log\\}$. However, in
Why3 and EventB, proof obligations are automatically generated from the
specifications, in particular the abstract and concrete models can be
naturally linked and the tool automatically generates the corresponding VC. In
that regard, Why3 and EventB are safer than $\\{log\\}$.

The possibility to count with executable code without much effort enables many
lightweight analysis that can be put into practice before attempting complex
proofs. In $\\{log\\}$ tool where specification and implementation are
described by only one piece of code (cf. forgrams). This tool is not the
integration of an interpreter and a prover; the same set of rewrite rules are
used to compute and prove. In EventB/Rodin there is only a specification—later
it can be converted into an executable representation if tools such as ProB
are used. Why3 can execute WhyML programs natively thanks to its interpreter
and the execute command. Furthermore, once the the program is proved to verify
the specification, correct-by-construction OCaml and C programs can be
automatically extracted. These programs will be orders of magnitude more
efficient than the equivalent $\\{log\\}$ forgrams.

##  7 Conclusion

We formally verified the implementation of sparse sets using three formal
languages and associated tools, focusing on the operations and correctness
properties required by a constraint solver when domains of integer variables
are implemented with sparse sets. We compared in particular the several
statements of invariants and pre-post properties and their proofs. As future
work, two directions can be investigated. The first one is to complete the
formal developments with other set operations. A second one is to implement
and verify, in Why3 or EventB, a labeling procedure such as the ones used in
constraint solvers, it would need to backtrack on the values of some domains,
and thus make use of the theorems proven in this paper. Labeling is native in
$\\{log\\}$when the CLP(FD) solver is active.

## References

  * [1] J. Abrial.  Modeling in Event-B - System and Software Engineering.  Cambridge University Press, 2010.
  * [2] J. Abrial, M. J. Butler, S. arxlerstede, T. S. Hoang, F. Mehta, and L. Voisin.  Rodin: an open toolset for modelling and reasoning in event-b.  Int. J. Softw. Tools Technol. Transf., 12(6):447–466, 2010.
  * [3] P. Briggs and L. Torczon.  An efficient representation for sparse sets.  LOPLAS, 2(1-4):59–69, 1993.
  * [4] D. Cantone, E. G. Omodeo, and A. Policriti.  Set Theory for Computing - From Decision Procedures to Declarative Programming with Sets.  Monographs in Computer Science. Springer, 2001.
  * [5] M. Cristiá, G. De Luca, and C. Luna.  An automatically verified prototype of the android permissions system.  Journal of Automated Reasoning, 67(2):17, May 2023.
  * [6] M. Cristiá and G. Rossi.  Solving quantifier-free first-order constraints over finite sets and binary relations.  J. Autom. Reason., 64(2):295–330, 2020.
  * [7] M. Cristiá and G. Rossi.  Automated proof of Bell-LaPadula security properties.  J. Autom. Reason., 65(4):463–478, 2021.
  * [8] M. Cristiá and G. Rossi.  Automated reasoning with restricted intensional sets.  J. Autom. Reason., 65(6):809–890, 2021.
  * [9] M. Cristiá and G. Rossi.  An automatically verified prototype of the Tokeneer ID station specification.  J. Autom. Reason., 65(8):1125–1151, 2021.
  * [10] M. Cristiá and G. Rossi.  A decision procedure for a theory of finite sets with finite integer intervals.  CoRR, abs/2105.03005, 2021.
  * [11] M. Cristiá and G. Rossi.  $\\{\mathit{log}\\}$: set formulas as programs.  Rend. Ist. Mat. Univ. Trieste, 53:24, 2021.  Id/No 23.
  * [12] M. Cristiá and G. Rossi.  A set-theoretic decision procedure for quantifier-free, decidable languages extended with restricted quantifiers.  CoRR, abs/2208.03518, 2022.  Under consideration in Journal of Automated Reasoning.
  * [13] M. Cristiá and G. Rossi.  Integrating cardinality constraints into constraint logic programming with sets.  Theory Pract. Log. Program., 23(2):468–502, 2023.
  * [14] D. Déharbe, P. Fontaine, Y. Guyot, and L. Voisin.  Integrating SMT solvers in rodin.  Sci. Comput. Program., 94:130–143, 2014.
  * [15] A. Dovier, C. Piazza, E. Pontelli, and G. Rossi.  Sets and constraint logic programming.  ACM Trans. Program. Lang. Syst., 22(5):861–931, 2000.
  * [16] J. Filliâtre and P. Letouzey.  Functors for proofs and programs.  In D. A. Schmidt, editor, 13th European Symposium on Programming, ESOP 2004, Held as Part of ETAPS 2004, Barcelona, Spain, Proceedings, volume 2986 of LNCS, pages 370–384. Springer, 2004.
  * [17] J. Filliâtre and A. Paskevich.  Why3 - where programs meet provers.  In M. Felleisen and P. Gardner, editors, 22nd European Symposium on Programming, ESOP 2013, Held as Part of ETAPS 2013, Rome, Italy, Proceedings, volume 7792 of LNCS, pages 125–128. Springer, 2013.
  * [18] V. Le clément de saint-Marcq, P. Schaus, C. Solnon, and C. Lecoutre.  Sparse-Sets for Domain Implementation.  In CP workshop on Techniques foR Implementing Constraint programming Systems (TRICS), pages 1–10, 2013.
  * [19] A. Ledein and C. Dubois.  Facile en coq : vérification formelle des listes d’intervalles.  In 31ème Journées Francophones des Langages Applicatifs, 2019\.
  * [20] D. Mentré, C. Marché, J. Filliâtre, and M. Asuka.  Discharging proof obligations from atelier B using multiple automated provers.  In J. Derrick, J. S. Fitzgerald, S. Gnesi, S. Khurshid, M. Leuschel, S. Reeves, and E. Riccobene, editors, Abstract State Machines, Alloy, B, VDM, and Z - Third Int. Conf., ABZ 2012, Pisa, Italy. Proceedings, volume 7316 of LNCS, pages 238–251. Springer, 2012.
  * [21] G. Rossi.  $\\{log\\}$.  http://www.clpset.unipr.it/setlog.Home.html, 2008.  Last access 2023.
  * [22] G. Rossi and M. Cristiá.  $\\{log\\}$ user’s manual.  Technical report, Dipartimento di Matematica, Universitá di Parma, 2020\.  https://www.clpset.unipr.it/SETLOG/setlog-man.pdf.
  * [23] J. T. Schwartz, R. B. K. Dewar, E. Dubinsky, and E. Schonberg.  Programming with Sets - An Introduction to SETL.  Texts and Monographs in Computer Science. Springer, 1986.


2021

[1]Dillon S. Fitzgerald

[1]Department of Physics, University of Michigan, 450 Church St, Ann Arbor,
48109, Michigan, USA

2]European Organization for Nuclear Research (CERN), Geneva, Switzerland

3]Universität Bonn – Helmholtz-Institut für Strahlen und Kernphysik, Bonn,
Germany

4]RWTH Aachen University, Aachen, Germany

# Ntuple Wizard: An Application to Access Large-Scale Open Data from LHCb

Christine A. Aidala  Christopher Burr  Marco Cattaneo<EMAIL_ADDRESS>Adam
Morris  Sebastian Neubert  Donijor Tropmann  *  [  [  [

###### Abstract

Making the large data sets collected at the Large Hadron Collider (LHC)
accessible to the world is a considerable challenge because of both the
complexity and the volume of data. This paper presents the Ntuple Wizard, an
application that leverages the existing computing infrastructure available to
the LHCb collaboration in order to enable third-party users to request
specific data. An intuitive web interface allows the discovery of accessible
data sets and guides the user through the process of specifying a
configuration-based request. The application allows for fine-grained control
of the level of access granted to the public.

###### keywords:

Open Data, Open Access, LHCb, LHC, CERN, HEP, Outreach

##  1 Introduction

In an increasingly diverse research landscape, management and curation of
public data are becoming critical components of transdisciplinary science.
Keys to the realization of an open research ecosystem that adds scientific
value have been identified in the FAIR principles of scientific data
management and stewardship fair_2016 . Making data Findable, Accessible,
Interoperable, and Reusable, however, requires a considerable amount of
tooling and infrastructure.

A common problem, which is acute for data in high-energy physics but
increasingly an issue in other fields as well, is the sheer size of data sets
stored in custom file formats. For large-scale experimental facilities, such
as the Large Hadron Collider (LHC) at the European Organization for Nuclear
Research (CERN), the data sets are so large that even access by the directly
involved scientists has to be centrally managed. LHCb is one of the four major
experiments at the LHC LHCb_det , and as an example, the data collected by the
LHCb experiment in the years 2011-12 (corresponding to $\sim 3$ fb-1 of
proton-proton collisions) amount to a volume of 900 TB of reconstructed data.
For the purpose of processing these data, extensive computing infrastructure
has been set up by the countries participating in this type of research
Bird:2014ctt . Making use of this infrastructure to enable external users FAIR
access to LHCb data in a resource-effective way is the goal of the application
presented in this paper.

In the following, an application is presented that exposes a data query
service to allow the public to request sub-samples of data collected and
published by the LHCb experiment. The samples are delivered as ROOT Ntuples
ROOT  a data format that requires no special LHCb-specific software to read
and for which converters to other standard file formats exist. We call the
application the Ntuple Wizard.

The Ntuple Wizard interface guides users with basic knowledge in particle
physics through the process of discovering the available data and formulating
a useful query. The queries can be processed by the existing data production
infrastructure, and results will be delivered through the CERN Open Data
Portal ODP . By splitting the data request into the construction of a data
query and subsequent processing of the query on the internal infrastructure,
the LHCb collaboration retains fine-grained control over access to the data.
Crucially this system protects the compute infrastructure from attacks by
malicious code injection.

###  1.1 Accessible open data

In 2020, the LHC experiments at CERN adopted a new Open Data Policy CERN-
OPEN-2020-013 , the scope of which expanded in 2022 to an Open Science Policy
CERN-OPEN-2022-013 . These documents define the commitments of CERN to make
the data collected at the LHC, at several levels of complexity, publicly
available DPHEPStudyGroup:2012dsv :

Level 1

Published results — this can include tables and figures but also preprocessed
Ntuples or binned and unbinned fit likelihood functions.

Level 2

Outreach and education — usually in the form of highly preprocessed Ntuples.

Level 3

Reconstructed data — these data have been preprocessed to derive physics
objects, such as charged particle candidates, photons, or particle jets.
Reconstructed data may or may not be corrected for detector effects, such as
efficiency and resolution.

Level 4

Raw data – the basic quantities recorded by the experimental instruments.

Both Level 1 and 2 data are considered to be highly processed, abstracted, and
manageable using commonly available computers. Level 4 raw data will not be
made available due to practical reasons concerning data size but also
detector-specific information needed for the interpretation of these data.
This leaves Level 3 data as the most versatile and basic data set which will
be publicly accessible.

All LHC experiments have long and intricate data reconstruction pipelines,
which yield several intermediate output data formats. During a pipeline, the
raw data are converted to physical objects such as charged particle
trajectories, jets, and vertices. Furthermore, the raw data are classified and
filtered to obtain samples enriched in interesting signatures.

Figure 1: LHCb data flow in Runs 1 and 2. The output of the stripping step
will be made public through the CERN Open Data Portal ODP .

Figure 1 shows an overview of the data processing pipeline in LHCb as it has
been used during LHC data-taking Runs 1 and 2 (2011–18). The various steps of
the pipeline are outlined further in the LHCb computing technical design
reports TDR1 ; TDR2 . Level 3 data have been defined as the output of the
stripping step. The stripping consists of a large number of selection
algorithms called lines, which are designed to filter the data and sort events
into several collections, which are called streams. Streams are defined
according to common physics signatures and aim to collect selections with
significant overlaps into a common set of files, to reduce duplication of the
data. The LHCb data organization is discussed in more detail in Appendix A,
including a list of streams available in Runs 1 and 2.

The stripping selections are based on the concept of physics candidates. A
candidate refers to a set of data matching a particular physics signature. In
most cases, this signature will be a particular particle decay, such as for
example $B^{+}\rightarrow\bar{D}^{0}\pi^{+}$ with the subsequent decay
$\bar{D}^{0}\rightarrow K^{+}\pi^{-}$, where $B,D,K,$ and $\pi$ mesons are the
lightest hadrons containing $b,c,s$, and $u/d$ quarks respectively. Such
cascading decays are represented as tree-like data structures, where the nodes
represent (intermediate) particles and the edges indicate a parent-child
relationship in the decay. These data structures are referred to as decay
trees (see Fig. 5 for an example). The root particle of the decay tree (the
$B^{+}$ in our example) is called its head. Stripping selections attempt to
find sets of physics objects in the reconstructed LHCb data, which match the
desired decay tree and any additional criteria that might be applied to
distinguish the intended signal process from background. Typical selection
criteria include kinematic variables, vertex and track reconstruction
qualities, and particle identification variables. Some stripping lines rely on
multivariate classifiers to combine several observables into a single powerful
selection criterion. The output of this procedure is collections of decay
candidates specified by their particular decay trees in a custom LHCb-specific
data format.

It is important to note that candidates are distinct from the concept of
events in the LHCb data processing. An event is defined during the data
acquisition and refers to a particular time window in which collisions can
occur. Several proton-proton collisions can happen during this time window,
and in principle, it can happen that several candidates for a particular decay
are identified for a single collision. In such cases, relevant quantities
(related to vertex reconstruction and flight distances) can be computed for
every primary vertex (e.g. collision point) in the event.

In order to convert these data into a framework-independent format a useful
concept is the aforementioned Ntuples. The idea of a Ntuple is simple: each
candidate is described by a tuple of variables, i.e. physical observables of
interest measured on the particular candidate, or referring to the global
event in which the candidate was found. A data set consists of $N$ such
tuples, much like a simple CSV file. Ntuples are saved in ROOT files ROOT  and
only basic data types are allowed for the variables. As a small complication
in some instances, the variables can actually be arrays of basic data types.
In such cases, the Ntuple Wizard provides the necessary documentation for the
interpretation.

###  1.2 Principle of Ntuple creation and the Ntuple Wizard

Both the stripping as well as the Ntuple-making step in Fig. 1 are handled by
DaVinci DaVinci ; TDR1 ; TDR2 , an LHCb application for event selection and
data analysis using the Gaudi framework Gaudi ; TDR1 ; TDR2 . DaVinci is
configured via Python scripts and used to process entire data sets with batch
processing. Both the Python configuration as well as the batch production
system are intentionally hidden from users of the Ntuple Wizard for security
reasons.

The DaVinci application provides access to a number of algorithms that can be
combined in sequence for event selection and processing. In order to produce a
Ntuple the user has to specify which variables should appear in the output
data. This Ntuple configuration is handled by an algorithm named
DecayTreeTuple, in which variables are registered through the use of so-called
TupleTools and LoKi functors. A large collection of those tools and functors
are available for the user to choose from. In general, a TupleTool will add a
set of variables to the Ntuple, while a LoKi functor usually computes a single
number. The LoKi::Hybrid::TupleTool can be used to write the output of
functors into the tuple. Functors can be combined with standard arithmetic and
logic operations, providing a flexible and powerful system to compute derived
quantities. A list of important available tools is presented in Appendix B.

Figure 2: Architecture of the Ntuple Wizard.

Figure 2 shows an overview of the Ntuple Wizard architecture, the core
functionality of which is the configuration of DaVinci. The metadata and
documentation describing the available data, preselections, as well as
available selection operations are generated from the original provenance
traces of the data and the stripping selection code. The web interface
presents this metadata and documentation to the user in a pedagogical way that
facilitates data discovery and formulation of the query. The query to the data
has two principal parts: Data set discovery and Ntuple configuration. First,
the Ntuple Wizard allows the user to select from the available predefined
stripping selections, data-taking periods, and conditions. In the second step,
the user defines what quantities should be computed and written to the output
Ntuple. Standard tools for the computation of typical quantities, such as
kinematic variables, particle identification (PID) variables, etc., are
available. The query formulated by the user is stored in a set of
configuration files. These files can be converted into a Python configuration
compatible with the internal LHCb Analysis Productions system
AnalysisProductions . This conversion and the final submission of the query to
the computing infrastructure are handled through an LHCb Analysis Productions
manager.

###  1.3 Security considerations

Accepting arbitrary external code to run on the LHCb computing resources has
obvious unacceptable security risks. Therefore, the Ntuple Wizard is designed
to generate the configuration in a pure data-structure format. As shown in
Figure 2, the configuration of the query is captured in YAML files, which can
be submitted to an LHCb Analysis Productions manager for further processing.

##  2 Metadata and documentation acquisition

In order to facilitate the core functionality of the Ntuple Wizard — namely
data set discovery (Sec. 4) and algorithm configuration (Sec. 5), metadata and
documentation from several sources are required. In particular, the Ntuple
Wizard needs to know what types of decays can be queried and what tools are
available to compute derived quantities of interest about these candidates.

Since these metadata are unchanging, and providing direct access to the
various sources requires authentication and introduces more points of failure,
the metadata are collated and served as static files over HTTP. No additional
access to the LHCb code or database is needed by the Ntuple Wizard once it has
been deployed.

The sources of metadata can be grouped into two coarse categories: the LHCb
software stack and the LHCb database. Metadata are acquired from the LHCb
software stack in two ways. The first is from the Gaudi Python interface;
particularly under the DaVinci application environment. Metadata about the
configuration interface of each TupleTool are extracted from DaVinci. Details
of the stripping lines, including the chain of selection algorithms that
define them, are extracted from the DaVinci versions used in the corresponding
stripping campaigns.

The process of building decay candidates in a stripping line often involves a
combination of many algorithms from the LHCb selection framework, which
combine particles, impose selection requirements, perform PID substitution,
and build final-state particle candidates from trajectories of charged
particles and calorimeter clusters. The algorithms can be related to each
other through their input and output locations. The full list of decays
(including all sub-decays) must be inferred by traversing the ‘dependency
tree’ of the selection algorithms. This is performed using custom code during
metadata acquisition.

The second, more indirect way is from the LHCb Doxygen pages, which themselves
are generated from the source code of the LHCb software stack. The latest
Doxygen pages for Run 1 or Run 2 DaVinci versions are used to extract the
documentation for each TupleTool and LoKi functor. A campaign to improve the
Doxygen documentation at its source was undertaken during the development of
the Ntuple Wizard.

The LHCb database provides metadata about the centrally managed data sets,
which is necessary to configure the Ntupling productions as explained above.
In order not to duplicate effort, a common code base is employed to extract
metadata from the LHCb database for both the Ntuple Wizard and the CERN Open
Data Portal.

##  3 User interface

The user interface consists of a sequence of dialogues that guide the user
through the configuration steps. This is designed as a client-side dynamic web
page that reads metadata acquired at deployment time to serve as static files
(see Sec. 2).

Since users of LHCb open data do not, in general, have access to the same
support network of experienced users and developers enjoyed by LHCb
collaboration members, a key design element of the Wizard is to provide the
necessary documentation for a novice user to complete each step of the
configuration process.

The existing documentation of DaVinci DaVinci ; TDR1 ; TDR2  is fragmented
across several sources (Twiki twiki , the Starterkit starterkit , Doxygen
doxygen  and the source code itself), so where possible, the Wizard pulls text
from each of these disparate sources and renders it in the relevant context
within the user interface.

There are two main steps to formulate a query using the Ntuple Wizard: Data
set discovery and Ntuple configuration. These steps are explained in the
following.

##  4 Data set discovery and production configuration

The available data contain a wide range of stripping selections, data-taking
periods, and running conditions. The Production configuration dialogue of the
Ntuple Wizard guides the user through the selection of the desired subsets.
The interface allows the selection of several decays to be processed
simultaneously as part of one query. For each decay, a separate Ntuple will be
produced.

###  4.1 Discovering available candidate decays

In the Decay search dialogue, the Ntuple Wizard presents a list of all decays
selected by the stripping, accompanied by decay descriptors in LoKi and LaTeX
formats, information about which stripping lines build them, as well as ‘tags’
that can be used to filter different types of decays. Decays are searchable
through various filters, including the identity or properties of the parent
particle and decay products, whether the candidates are built by a specific
stripping line, and the aforementioned tags. An example of the decay search is
shown in Figure 3. The selected candidate of interest is highlighted in blue,
and the collection was narrowed down from the list of all possible decays by
using the filters and tags at the top of the page. The ‘none of’ option of the
‘tags’ drop-down menu is chosen by default, indicating that decays with the
displayed tags are hidden from the list of selectable decays. The tags
‘charge-violating‘ and ‘undefined-unstable‘ corresponding to decays that
violate charge conservation and that contain unstable particles without
defined decays respectively are hidden by default. If the user wishes to
instead isolate decays that meet the criteria of any or all of the selected
tags, the options ‘any’ or ‘all’ can be chosen from the ‘tags’ drop-down menu.
It is possible to select several decays for further processing at this stage.

Figure 3: Example of the decay candidate search function of the Ntuple Wizard.

###  4.2 Stripping line and data set selection

Once a decay is selected by the user, all corresponding stripping lines and
data sets from the various running periods are listed, and the desired
combination(s) can be selected. The case can arise where the same stripping
line shows up in multiple stripping versions within the same data set (stream,
running year, and magnet polarity). These are rendered as separate options in
the data set selection drop-down menu of the Ntuple Wizard. For a given decay,
it is recommended to choose only one data set for each magnet polarity within
a given running year, and to use the most recent stripping version in the case
of duplicates. The data organization of LHCb is elaborated on in Appendix A,
including a table of running years, as well as corresponding collision
energies and stripping versions.

Links to documentation about each stripping line including selection
algorithms that went into building the decay candidates are displayed to the
user to guide them in choosing the most suitable stripping line and data
stream for their physics interest. Figure 4 shows an example of the production
configuration page, where an available stripping line and data set have been
chosen from lists of all possibilities corresponding to the selected decay
channel. The blue question mark button contains links to the aforementioned
stripping documentation. At this point, the query is specified up to deciding
what information to write into the Ntuple.

Figure 4: Example of the data set selection and production configuration step
of the Ntuple Wizard.

##  5 Ntuple configuration

The DecayTreeTuple configuration dialogue is designed to guide the user
through customization of the quantities written to the Ntuple for the selected
candidates. For each decay, a separate DecayTreeTuple has to be configured.
Care should be taken to name the Ntuples appropriately. The Ntuple Wizard
requires a unique name for each Ntuple.

Figure 5: Example of an interactive graph used to configure DecayTreeTuple,
with selected TupleTools displayed for both the entire candidate (top) and
selected nodes (bottom).

Selected decay trees are visually represented as graphs, where each physics
object (e.g. particle) is represented by a node as shown in the screenshots in
Figure 5. The user can interact with this graph by selecting one or multiple
nodes at a time and determining which TupleTools will be added to each node,
which in turn determines which quantities are saved to the Ntuple. A list is
rendered on screen depending on the selected node(s), each element of which
corresponds to a selected TupleTool, with buttons for configuring and removing
the tool. The TupleTool configuration interface includes links to relevant
documentation about the tool, including lists of quantities written by the
tool where available. Each node in the graph comes with the standard set of
TupleTools for LHCb analyses, but more will often be needed depending on the
particular physics interests of the user. Furthermore, any added tool will
come with the standard configuration, which can be further modified if the
user desires. A custom set of standard LoKi variables and functions of these
variables can also be saved to the Ntuple for each node, using the
Loki::Hybrid::TupleTool. Appendix B contains a brief description of the
standard set of TupleTools included with each node on the graph, as well as
other useful TupleTools for physics analysis. Figure 5 shows an example of the
configurable graph corresponding to the selected candidate shown in Figures 3
and 4, as well as a list of TupleTools corresponding to the entire decay
candidate (top), and particular nodes selected on the graph (bottom). It can
be seen from the figure that nodes can also be selected through the categories
shown below the graph and that TupleTools can be added, removed, or configured
for each node or grouping of nodes.

Figure 6 shows an example of the user interface for configuring TupleTools,
with the particular example showing TupleToolTISTOS, which saves trigger
information to the Ntuple. It can be seen at the bottom how relevant
information is provided.

Figure 6: Example of the configuration interface of a TupleTool within the
Ntuple Wizard, (in particular, TupleToolTISTOS for saving trigger
information), including links to relevant documentation at the bottom of the
modal.

###  5.1 Configuration output

Figure 7 shows an example of the output YAML file used to configure the
DecayTreeTuple algorithm that was populated via configurations captured in
Figs. 5 – 6, where the tools, groups and branches keys are shown specifying
which TupleTools and therefore which information will be saved to the Ntuple.
The top-level key tools contains a list of TupleTool configurations, from
which the parsing functions create and configure TupleTool algorithms attached
to the DecayTreeTuple itself, which will thus write either particle-level
information about the decay or event-level information, depending on the class
of the TupleTool. The keys branches and groups themselves contain lists of
dictionaries whose keys specify particles and have their own tools lists which
are used similarly to attach TupleTool algorithms to the specified particle(s)
in the decay tree. Note that groups differs from branches in that it specifies
multiple particles to be looped over and have identically configured TupleTool
algorithms attached.

Figure 7: Output of Btree.yaml, the data file used to configure the
DecayTreeTuple algorithm.

⬇

inputs:

- /Event/BhadronCompleteEvent/Phys/B2D0PiD2HHBeauty2CharmLine/Particles

descriptorTemplate: ${Bplus}[B+ -> ${D_0}(D~0 -> ${Kplus}K+
${piminus}pi-)${piplus}pi+]CC

tools:

- TupleToolKinematic:

ExtraName: ’’

Verbose: false

MaxPV: 100

Transporter: ParticleTransporter:PUBLIC

- TupleToolPid:

ExtraName: ’’

Verbose: false

MaxPV: 100

- TupleToolANNPID:

ExtraName: ’’

Verbose: false

MaxPV: 100

ANNPIDTunes:

- MC12TuneV2

- MC12TuneV3

- MC12TuneV4

- MC15TuneV1

PIDTypes:

- Electron

- Muon

- Pion

- Kaon

- Proton

- Ghost

- TupleToolGeometry:

ExtraName: ’’

Verbose: false

MaxPV: 100

RefitPVs: false

PVReFitter: LoKi::PVReFitter:PUBLIC

FillMultiPV: false

- TupleToolEventInfo:

ExtraName: ’’

Verbose: false

MaxPV: 100

branches:

Bplus:

particle: B+

tools: []

⬇

D_0:

particle: D~0

tools: []

Kplus:

particle: K+

tools: []

piminus:

particle: pi-

tools: []

piplus:

particle: pi+

tools: []

groups:

Kplus,piminus:

particles:

- K+

- pi-

tools:

- TupleToolTISTOS:

ExtraName: ’’

Verbose: false

MaxPV: 100

VerboseL0: false

VerboseHlt1: false

VerboseHlt2: false

VerboseStripping: false

FillL0: true

FillHlt1: true

FillHlt2: true

FillStripping: false

TriggerList: []

Hlt1TriggerTisTosName: Hlt1TriggerTisTos

Hlt2TriggerTisTosName: Hlt2TriggerTisTos

L0TriggerTisTosName: L0TriggerTisTos

PIDList: []

TopParticleOnly: false

Hlt1Phys: >-

Hlt1(?!ODIN)(?!L0)(?!Lumi)(?!Tell1)(?!MB)(?!NZS)(?!Velo)(?!BeamGas)(?!Incident).*Decision

Hlt2Phys: >-

Hlt2(?!Forward)(?!DebugEvent)(?!Express)(?!Lumi)(?!Transparent)(?!PassThrough).*Decision

TIS: true

TOS: true

TUS: false

TPS: false

name: DecayTreeTuple/Btree

###  5.2 Future developments

It is planned to extend the current functionality of the Ntuple Wizard by
including the ability to create custom candidates from standard collections of
LHCb particles. Another important planned addition is the ability to configure
custom jet reconstruction. Ideally, support will be included for the full set
of algorithms available in DaVinci for data analysis and event/candidate
selection as resources allow.

As of Run 3, which started in 2022, the majority of the filtering and
preselection of the data will be done in real time within the LHCb high-level
trigger (HLT). In this architecture, the data will be fully reconstructed
online and the final preselection algorithms will run in the HLT. Offline
preselections will be feasible for a subset of the events. In both cases the
output will have the same level of abstraction as the output of the stripping,
allowing for a relatively simple adaptation of the Ntuple Wizard once the Run
3 data are made public.

##  6 Request submission and execution

Once the candidate(s) of interest, data set(s), and information to be saved in
the Ntuple(s) are specified, and a name and email address have been provided
for the production, the configuration files can be submitted to an LHCb
Analysis Productions manager.

The Ntuple Wizard will be integrated with the CERN Open Data Portal, including
an interface that allows users to make requests, check the status of their
productions, and retrieve their output Ntuples all in one place. This will
also provide a means to keep track of users through the CERN Open Data Portal
database system.

Requests for Ntuple creation are handled using the Analysis Productions
package. The files describing a new request are committed to a repository
hosted on the CERN GitLab git , and a merge request is created once they are
ready for review. The Continuous Integration feature of GitLab is used to
submit test productions to LHCbDIRAC TDR1 ; TDR2 , which automatically
processes a small fraction of the data when the remote repository is updated.

Once the request is submitted, it is handled by the LHCbDIRAC production
system. A production defines how a data set is to be processed, and LHCbDIRAC
will launch and manage computing jobs until the data set is fully processed.
Productions are defined in ‘steps’ that specify which application to run and
which configuration files to read, and may be chained together such that the
output of the first step is the input to the second, etc. The info.yaml file
produced by the Ntuple Wizard defines one production per data set, each
consisting of a single DaVinci step.

Within the production jobs, DaVinci is configured by functions defined in an
external Python module according to the YAML files produced by the Ntuple
Wizard. The data structure configured in Section 5 and displayed in Figure 7
is traversed, and the configurable properties of the DecayTreeTuple algorithm
are assigned the corresponding values.

After the Analysis Production jobs are complete, the produced Ntuples will be
delivered to the CERN Open Data Portal for retrieval.

##  7 Summary

Providing public access to the large data sets at the LHC is a significant
technical challenge, but it is becoming increasingly important for the
longevity of high-energy physics in order to optimize acquired knowledge from
the collected data. The volume and complexity of the data collected at LHCb
make providing direct access to reconstructed (Level 3) data suitable for
physics research difficult, motivating the design of the Ntuple Wizard, where
users can submit queries to obtain skimmed data samples (Ntuples) of the
reconstructed data suitable for their physics interests. The Ntuple Wizard is
a web-based application that intuitively guides the user through specifying a
query, from discovering a data set from a physics candidate (e.g. decay) of
interest to configuring the information to be saved in the output Ntuple. The
output of the Ntuple Wizard is a pure data structure (YAML) format, which is
to be submitted to an LHCb Analysis Productions manager so it can be parsed
internally to provide the necessary Python scripts needed to configure the
DaVinci application. The Ntuples will ultimately be delivered to the CERN Open
Data Portal for retrieval.

## Appendices

##  Appendix A LHCb Data Organization

Table 1 shows a list of running years, including corresponding collision
energies and stripping versions available in the Ntuple Wizard.

Table 1: Table of running years, including collision energy ($\sqrt{s}$) and relevant stripping versions available in the Ntuple Wizard. Running Year |  $\sqrt{s}$ (TeV) | Stripping Versions
---|---|---
Run 1
2011 | 7 | s21r1, s21r1p1, s21r1p2
2012 | 8 | s21, s21r0p1, s21r0p2
Run 2
2015 | 13 | s24r2
2016 | 13 | s28r2, s28r2p1
2017 | 13 | s29r2, s29r2p1, s29r2p2
2018 | 13 | s34, s34r0p1, s34r0p2

LHCb data streams come in two formats, Data Summary Tape (DST) files, which
contain the full event information, and micro Data Summary Tape (MDST) files,
which only contain the physics objects directly matched in at least one
stripping selection. In MDST streams, the rest of the information in the
events, apart from a few global event properties, is discarded. Table 2 shows
a list of streams from Run 1 and Run 2, with the DST vs MDST format indicated
in the stream name.

Table 2: Table of data streams from Runs 1 and 2 available through the Ntuple
Wizard. Stream
---
BHADRON.MDST
BHADRONCOMPLETEEVENT.DST
CALIBRATION.DST
CHARM.MDST
CHARMCOMPLETEEVENT.DST
CHARMTOBESWUM.DST
DIMUON.DST
EW.DST
LEPTONIC.MDST
MINIBIAS.DST
PID.MDST
RADIATIVE.DST
SEMILEPTONIC.DST

##  Appendix B List of useful TupleTools

###  B.1 Default TupleTools

A set of TupleTools is included by default for all Ntuple configuration files
produced by the Ntuple Wizard. These tools can be removed by the user if
desired, but are standard tools used in LHCb analyses, and are recommended to
keep for physics analyses. Given the flexible data structure of the Ntuple, it
is easy to produce a reduced data structure with a subset of variables at a
later stage in data processing, while still maintaining the full set of
variables in the original Ntuple.

  * •

TupleToolANNPID — A tool used to add artificial neural network particle
identification information about the physics candidate to the Ntuple.

  * •

TupleToolEventInfo — A tool used to add event and run information to the
Ntuple.

  * •

TupleToolGeometry — A tool used to add information about physics candidate
geometry and event geometry to the Ntuple.

  * •

TupleToolKinematic — A tool used to add kinematic information about the
physics candidate to the Ntuple.

  * •

TupleToolPid — A tool used to add particle identification information about
the physics candidate to the Ntuple, with additional information than that in
TupleToolANNPID, including information about which PID detector subsystems
were used in the probability calculations.

###  B.2 Other useful TupleTools

  * •

TupleToolTISTOS — A tool that saves the trigger TIS/TOS (Trigger independent
of Signal/Trigger on Signal) decisions for selected particles to the Ntuple.

  * •

LoKi::Hybrid::TupleTool — A tool that allows the user to add LoKi variables or
expressions of these variables known as LoKi functors to the Ntuple.

## Acknowledgements

We thank our colleagues at LHCb for providing the necessary data and software,
and within the Data Processing & Analysis (DPA) project for their incredibly
valuable discussions. We additionally would like to thank Jose Marco Arias for
systematic testing of the web interface. All authors acknowledge support from
CERN. In addition, C.A.A. and D.S.F. acknowledge support from the National
Science Foundation under Award Number 2012926, and A.M. and S.N. acknowledge
support from the DFG Grant NE 2185/1-1.

## References

  * (1) M. D. Wilkinson et al.  The FAIR Guiding Principles for scientific data management and stewardship.  Scientific Data, 3(1):160018, March 2016.
  * (2) A. Augusto Alves, Jr. et al.  The LHCb Detector at the LHC.  JINST, 3:S08005, 2008.
  * (3) I. Bird et al.  Update of the Computing Models of the WLCG and the LHC Experiments.  (CERN-LHCC-2014-014, LCG-TDR-002), Apr 2014.
  * (4) I. Antcheva et al.  ROOT: A C++ framework for petabyte data storage, statistical analysis and visualization.  Comput. Phys. Commun., 180:2499–2512, 2009.
  * (5) CERN Open Data Portal.  https://opendata.cern.ch/.
  * (6) CERN Open Data Policy for the LHC Experiments.  (CERN-OPEN-2020-013), Nov 2020.  https://cds.cern.ch/record/2745133.
  * (7) CERN Open Science Policy.  (CERN-OPEN-2022-013), 2022.  https://cds.cern.ch/record/2835057.
  * (8) Z. Akopov et al.  Status Report of the DPHEP Study Group: Towards a Global Effort for Sustainable Data Preservation in High Energy Physics.  (DPHEP-2012-001, FERMILAB-PUB-12-878-PPD, arXiv:1205.4667), May 2012.
  * (9) The LHCb Collaboration.  LHCb computing: Technical Design Report.  (CERN-LHCC-2005-019; LHCb-TDR-11), 2005.  http://cds.cern.ch/record/835156.
  * (10) The LHCb Collaboration.  Upgrade Software and Computing.  (CERN-LHCC-2018-007, LHCB-TDR-017), 2018.  https://cds.cern.ch/record/2310827.
  * (11) DaVinci.  https://lhcbdoc.web.cern.ch/lhcbdoc/davinci/.
  * (12) Gaudi Project.  https://gaudi.web.cern.ch/gaudi/.
  * (13) M. Ferrillo.  New generation offline software for the LHCb upgrade I.  2022\.  https://cds.cern.ch/record/2806414.
  * (14) LHCb Twiki.  https://twiki.cern.ch/twiki/bin/view/LHCb/.
  * (15) LHCb Starterkit.  https://lhcb.github.io/starterkit-lessons/index.html.
  * (16) LHCb Doxygen.  https://lhcb.web.cern.ch/computing/Support/Doxygen/doxygen.htm.
  * (17) GitLab.  https://gitlab.cern.ch/.