extfieldsofstudy:
Mathematics
Computer Science

text:
Clear and Compress: Computing Persistent Homology in Chunks

We present a parallelizable algorithm for computing the persistent homology of a filtered chain complex. Our approach differs from the commonly used reduction algorithm by first computing persistence pairs within local chunks, then simplifying the unpaired columns, and finally applying standard reduction on the simplified matrix. The approach generalizes a technique by G\"unther et al., which uses discrete Morse Theory to compute persistence; we derive the same worst-case complexity bound in a more general context. The algorithm employs several practical optimization techniques which are of independent interest. Our sequential implementation of the algorithm is competitive with state-of-the-art methods, and we improve the performance through parallelized computation.

Introduction
Persistent homology has developed from a theoretical idea to an entire research area within the field of computational topology. One of its core features is its multi-scale approach to analyzing and quantifying topological features in data. Recent examples of application areas are shape classification [2], topological denoising [1], or developmental biology [7].
A second major feature of persistent homology is the existence of a simple yet efficient computation method: The standard reduction algorithm as described in [8,17] computes the persistence pairs by a simple sequence of column operations; Algorithm 1 gives a complete description in just 10 lines of pseudo-code. The worst-case complexity is cubic in the input size of the complex, but the practical behavior has been observed to be closer to linear on average. Various variants have been proposed in order to improve the theoretical bounds [14,3] or the practical behavior [5].
Our first contribution consists of two simple optimization techniques of the standard reduction algorithm, which we call clearing and compression. Both approaches exploit the special structure of a filtered chain complex in order to significantly reduce the number of operations on real-world instances. However, the two methods cannot be easily combined because they require the columns of the boundary matrix to be processed in different orders.
Our second contribution is a novel algorithm that incorporates both of the above optimization techniques, and is also suitable for parallelization. It proceeds in three steps: In the first step, the columns of the matrix are partitioned into consecutive chunks. Each chunk is reduced independently, applying the clearing optimization mentioned above. In this step, the algorithm finds at least all persistence pairs with (index) persistence less than the size of the smallest chunk; let g be the number of columns not paired within this first step. In the second step, the g unpaired columns are compressed using the method mentioned above. After compression, each column has at most g non-zero entries and the unpaired columns form a nested (g × g)-matrix. In the third and final step, this nested matrix is reduced, again applying the clearing optimization.
The chunk algorithm is closely related to two other methods for computing persistence. First of all, the spectral sequence algorithm [6,§VII.4] decomposes the matrix into blocks and proceeds in several phases, computing in phase r the persistence pairs lying r − 1 blocks apart. The first step of the chunk algorithm is equivalent to applying the first two phases of the spectral sequence approach. Furthermore, the three step chunk algorithm is inspired by the approach of Günther et al. [10], which combines persistence computation and discrete Morse theory for 3D image data. The first step of that algorithm consists in constructing a discrete gradient field consistent with the input function; such a gradient field can be interpreted as a set of persistence pairs that are incident in the complex and have persistence 0. We replace this method by local persistence computations, allowing us to find pairs with are not incident in the complex.
We analyze the chunk algorithm in terms of time complexity. Let n be the number of generators (simplices, cells) of the chain complex, m the number of chunks, the maximal size of a chunk, and g as above. We obtain a worst-case bound of where the three terms reflect the worst-case running times of the three steps. For the filtration of a cubical complex induced by a d-dimensional grayscale image (with d some fixed constant), this bound simplifies to if the chunks are given by the cells appearing simultaneously in the filtration. This bound improves on the previous general bound of O(g 2 n log n) from [4]. Moreover, it matches the bound in [10], but applies to arbitrary dimensions. Of course, the bound is still cubic if expressed only in terms of n because g ∈ O(n) in the worst case. We implemented a sequential and a parallelized version of the chunk algorithm; both are publicly available in our new PHAT library (http://phat.googlecode.com/). The sequential code already outperforms the standard reduction algorithm and is competitive to other variants with a good practical behavior. The parallelized version using 12 cores yields a speed-up factor between 3 and 11 (depending on the example) in our tests, making the implementation the fastest among the considered choices. This is the first result where the usefulness of parallelization is shown for the problem of persistence computation through practical experiments.

Background
This section summarizes the theoretical foundations of persistent homology as needed in this work. We limit our scope to simplicial homology over Z 2 just for the sake of simplicity in the description; our methods generalize to chain complexes over arbitrary fields.
Homology Homology is an algebraic tool for analyzing the connectivity of topological spaces. Let K be a simplicial complex of dimension d. In any dimension p, we call a pchain a formal sum of the p-simplices of K with Z 2 coefficients. The p-chains form a group called the pth chain group C p . The boundary of a p-simplex σ is the (p − 1)-chain formed by the sum of all faces of σ of codimension 1. This operation extends linearly to a boundary operator δ : C p → C p−1 . A p-chain γ is a p-cycle if δ (γ) = 0. The p-cycles form a subgroup of the p-chains, which we call the pth cycle group Z p . A p-chain γ is called a p-boundary if γ = δ (ξ ) for some (p + 1)-chain ξ . Again, the p-boundaries form a group B p , and since δ (δ (ξ )) = 0 for any chain ξ , p-boundaries are p-cycles, and so B p is a subgroup of Z p . The pth homology group H p is defined as the quotient group Z p /B p . The rank of H p is denoted by β p and is called the pth Betti number. In our case of Z 2 coefficients, the homology group is a vector space isomorphic to Z β p 2 , hence it is completely determined by the Betti number. Roughly speaking, the Betti numbers in dimension 0, 1, and 2 yield the number of connected components, tunnels, and voids of K, respectively.
Persistence Let {σ 1 , . . . , σ n } denote the simplices of K. We assume that for each i ≤ n, For every dimension p and every K i , we have a homology group H p (K i ); we usually consider all dimensions at once and write H(K i ) for the direct sum of the homology groups of K i in all dimensions. The inclusion K i → K i+1 induces a homomorphism g i+1 i : H(K i ) → H(K i+1 ) on the homology groups. These homomorphisms compose and we can define g j i : H(K i ) → H(K j ) for any i ≤ j. We say that a class α ∈ H(K ) is born at (index) i if α ∈ im g i but α / ∈ im g i−1 . A class α born at index i dies entering (index) j if g j i (α) ∈ im g j i−1 but g j−1 i (α) / ∈ im g j−1 i−1 . In this case, the index pair (i, j) is called a persistence pair, and the difference j − i is the (index) persistence of the pair. The transition from K i−1 to K i either causes the birth or the death of an homology class. We call the added simplex σ i positive if it causes a birth and negative if it causes a death. Note that homology classes of the full complex K do not die during the filtration. We call a simplex σ i that gives birth to such a class essential. All other simplices are called inessential.
Boundary matrix For a matrix M ∈ Z n×n 2 , we let M j denote its j-th column, M i its ith row, and M i j ∈ Z 2 its entry in row i and column j. For a non-zero column 0 = M j = (m 1 , . . . , m n ) ∈ Z n 2 , we set pivot(M j ) := max{i = 1, . . . , n | m i = 1} and call it the pivot index of that column.
The boundary matrix D ∈ (Z 2 ) n×n of a simplexwise filtration (K i ) i is a n × n matrix with D i j = 1 if and only if σ i is a face of σ j of codimension 1. In other words, the jth column of D encodes the boundary of σ j . D is an upper-triangular matrix because any face of σ j must precede σ j in the filtration. Since the jth row and column of D corresponds to the jth simplex σ j of the filtration, we can talk about positive columns, negative columns, and essential columns in a natural way, and similarly for rows.
The reduction algorithm A column operation of the form M j ← M j + M k is called leftto-right if k < j. We call a matrix M derived from M if M can be transformed into M by left-to-right operations. Note that in a derivation M of M, the jth column can be expressed as a linear combination of the columns 1, . . . , j of M, and this linear combination includes M j . We call a matrix R reduced if no two non-zero columns have the same pivot index. If R is derived from M, we call it a reduction of M. In this case, we define Although the reduction matrix R is not unique, the sets P R and E R are the same for any choice of reduction; therefore, we can define P M and E M to be equal to P R and E R for any reduction R of M. We call the set P the persistence pairs of M. When obvious from the context, we omit the subscripts and simply write P for the persistence pairs. For the boundary matrix D of K, the pairs (i, j) ∈ P are the persistence pairs of the filtration (K i ) 0≤i≤n , and the indices in E correspond to the essential simplices of the complex. Note that E is uniquely determined by P and n as the indices between 1 and n that do not appear in any pair of P.
The simplest way of reducing D is to process columns from left to right; for every column, other columns are added from the left until the pivot index is unique (Algorithm 1). A lookup table can be used to identify the next column to be added in constant time. A flag is used for every column denoting whether a persistence pair with the column index has already been found. After termination, the unpaired columns correspond to the essential columns. The running time is at most cubic in n, and this bound is actually tight for certain input filtrations, as demonstrated in [16].
Mark columns i and j as paired and add (i, j) to P 10: With this definition, a matrix M is a reduction of D if and only if every column is reduced.

Speed-ups
Algorithm 1 describes the simplest way of reducing the boundary matrix, but it performs more operations than actually necessary to compute the persistence pairs. We now present two simple techniques which both lead to a significant decrease in the number of required operations.
Clearing positive columns The key insight behind our first optimization is the following fact: if i appears as the pivot in a reduced column of M, the index i is positive and hence there exists a sequence of left-to-right operations on M i that turn it to zero. Instead of explicitly executing this sequence of operations, we define the clear operation by setting column M i to zero directly. Informally speaking, a clear is a shortcut to avoid some column operations in the reduction when the result is evident.
In order to apply this optimization, we change the traversal order in the reduction by first reducing the columns corresponding to simplices with dimension d (from left to right), then all columns with dimension d − 1, and so on. After having reduced all columns with dimension δ , we have found all positive inessential columns with dimension δ − 1 and clear them before continuing with δ − 1. This way all positive inessential columns of the complex are cleared without performing any column additions on them. See [4] for a more detailed description.
Compression Alternatively, we can try to save arithmetic operations by reducing the number of non-zero rows among the unpaired columns. A useful observation in this context is given next. Proof. The statement is clearly true if M j is reduced, because in this case (i, j) is a persistence pair. If M j is not reduced, this means that after applying some sequence of left-to-right column operations, some reduced column has i as pivot index. As a consequence, whenever a negative column with index j has been reduced, row j can be set to zero before further reducing. We can even do more: let i be the pivot index of the reduced column M j and assume that the submatrix of M with column indices {1, . . . , j} and row indices {i, . . . , n} is reduced, i.e., the pivot indices are unique in this submatrix. By adding column j to each unreduced column in the matrix that has a non-zero entry at row i, we can eliminate all non-zero entries in row i from the unreduced columns. Note that if k < j and M i k = 0, then pivot(M k ) ≥ i and thus, by assumption, M k must be a reduced negative column. Therefore, for each unreduced column M k , the operation M k ← M k + M j is a left-to-right addition and thus does not affect the pairs.

Reduction in chunks
The two optimization techniques from Section 3 both yield significant speed-ups, but they are not easily combinable, because clearing requires to process a simplex before its faces, whereas compression works in the opposite direction. In this section, we present an algorithm which combines both optimization techniques.
Let m ∈ N. Fix m + 1 numbers 0 = t 0 < t 1 < . . . < t m−1 < t m = n and define the ith chunk of D to be the columns of D with indices {t i−1 + 1, . . . ,t i }. We call a column D j local if it forms a persistence pair with another column in the same chunk or in one of the adjacent chunks. In this case, we also call the persistence pair local. Non-local columns (and pairs) are called global. If is a lower bound on the size of each chuck, then every global persistence pair has index persistence at least . We also call an index j local if the jth column of D is local, and the same for global. We denote the number of global columns in D by g. The high-level description of our new algorithm consists of three steps: 1. Partially reduce every chunk independently, applying the clearing optimization, so that all local columns are completely reduced.
2. Independently compress every global column such that all its non-zero entries are global.
3. Reduce the submatrix consisting only of the global rows and columns.
We give details about the three steps in the rest of this section. The first two steps can be performed in parallel, whereas the third step only needs to reduce a matrix of size g × g instead of n × n. In many situations, g is significantly smaller than n.
Local chunk reduction The first step of our algorithm computes the local pairs by performing two phases of the spectral sequence algorithm [6,§VII.4]. Concretely, we apply left-to-right operations as usual, but in the first phase we only add columns from the same chunk, and in the second phase we only add columns from both the same chunk and its left neighbor. After phase r, for each b ∈ {r, . . . , m} the submatrix with column indices {1, . . . ,t b } and row indices {t b−2 +1, . . . , n} is reduced. If the reduction of column j stops at a pivot index i > t b−r , row j cannot be reduced any further by adding any column, so we identify (i, j) as a local persistence pair. Conversely, any local pair (i, j) is detected by this method after two phases. We incorporate the clearing operation for efficiency, that is, we proceed in decreasing dimension and set detected local positive columns to zero; see Algorithm 2. After its execution, L[i] contains the index of the local negative column with pivot index i for any local positive column i, and the resulting matrix R is a derivation of D in which all local columns are reduced. for r = 1, 2 do Perform two phases of the spectral sequence algorithm 5: for b = r, . . . , m do Loop is parallelizable 6: for if j is not marked as paired then 8: 10: if R j = 0 then 11: i ← pivot(R j ) Mark i and j as paired and add (i, j) to P 16: return (R, L, P) if k is marked as active/inactive then return true/false 6: for each non-zero entry index i of R k do 7: if is unpaired then 8: mark k as active and return true 9: else if i is positive then 10: (i, j) is persistence pair 11: if j = k and MARK COLUMN(R, j) then 12: mark k as active and return true 13: mark k as inactive and return false Global column compression Let R be the matrix returned by Algorithm 2. Before computing the global persistence pairs, we first compress the global columns, using the ideas from Section 3; recall that negative rows can simply be set to zero, while entries in positive rows can be eliminated by an appropriate column addition. Note, however, that a full column addition might actually be unnecessary: for instance, if all non-zero row indices in the added column belong to negative columns (except for the pivot), the entry in the local positive row could just have been zeroed out in the same way as in Corollary 3. Speaking more generally, it is more efficient to avoid column additions that have no consequences for global indices, neither directly nor indirectly.
In the spirit of this observation, we call an index i inactive if either it is a local negative index or if (i, j) is a local pair and all indices of non-zero entries in column R j apart from i are inactive. Otherwise, the index is called active. By induction and Corollary 3, we can show: The compression proceeds in two steps: first, every non-zero entry of a global column is classified as active or inactive (using depth-first search; see Algorithm 3). Then, we iterate over the global columns, set all entries with inactive index to zero, and eliminate any non-zero entry with a local positive index by column addition with L[ ] (see Algorithm 4). After this process, we obtain a matrix R with the same persistence pairs as R, such that the global columns of R have non-zero entries only in the global rows.
Submatrix reduction After having compressed all global columns, these form a g × g matrix "nested" in R (recall that g is the number of global columns). To complete the computation of the persistence pairs, we simply perform standard reduction on the remaining matrix. For efficiency, we perform steps 2 and 3 alternatingly for all dimensions in decreasing order and apply the clearing optimization; this way, we avoid the compression of positive global columns. Algorithm 5 summarizes the whole method. for each non-zero entry index of R k in decreasing order do 4: if is paired then 5: if is inactive then for j = 1, . . . , n with dim σ j = δ do Loop is parallelizable 7: if column j is not paired then for j = 1, . . . , n with dim σ j = δ do step 3: reduce global columns 10:

12:
if R j = 0 then Mark i and j as paired and add (i, j) to P 17: return P

Analysis
Algorithm 5 permits a complexity analysis depending on the following parameters: n, the number of simplices; m, the number of chunks; , the maximal size of a chunk; and g, the number of global columns. We assume that for any simplex, the number of non-zero faces of codimension 1 is bounded by a constant (this is equivalent to assuming that the dimension of the complex is a constant).

General complexity We show that the complexity of Algorithm 5 is bounded by
The three summands correspond to the running times of the three steps 1 . Note that g ∈ O(n) in the worst case. For the complexity of Algorithm 2, we consider the complexity of reducing one chunk, which consists of up to columns. Within the local chunk reduction, every column is only added with columns of the same or the previous chunk, so there are only up to 2 column additions per column. Moreover, since the number of non-zero entries per column in D is assumed to be constant, there are only O( ) many entries that can possibly become non-zero during the local chunk reduction. It follows that the local chunk reduction can be considered as a reduction on a matrix with columns and O( ) rows. If we represent columns by linked lists (containing the non-zero indices in sorted order), one column operation can be done in O( ) primitive operations, which leads to a total complexity of O( 3 ) per chunk.
The computation of active columns in Algorithm 3 is done by depth-first search on a graph whose vertices are given by the columns and whose edges correspond to their non-zero entries. The number of edges is O(n ), so we obtain a running time of O(n ).
Next, we consider the cost of compressing a global column with index j. After the previous step, the column has at most O( ) non-zero entries. We transform the presentation of the column from a linked list into a bit vector of size n. In this representation, adding another column in list representation with v entries to column j takes time proportional to v. In the worst case, we need to add all columns with indices 1, . . . , j − 1 to j. Each such column has O( ) entries. At the end, we transform the bit vector back into a linked list representation. The total cost is O(n + ( j − 1) + n) = O(n ) per global column.
Finally, the complexity of the global reduction is O(g 3 ), as in the standard reduction.
Choosing chunks We discuss different choices of chunk size and their complexities. A generic choice for an arbitrary complex is to choose O( √ n) chunks of size O( √ n) each. With that, the complexity of (1) becomes O(n 2 + g 1 n √ n + g 3 1 ).
Alternatively, choosing O( n log n ) chunks of size O(log n), the complexity becomes O(n log 2 n + g 2 n log n + g 3 2 ).
We replaced g by g 1 and g 2 to express that the number of global columns is different in both variants. In general, choosing larger chunks is likely to produce less global columns, since every global persistence pair has index persistence at least (the size of the smallest chunk) .
Cubical complexes We consider an important special case of boundary matrices: consider a d-dimensional image with p hypercubes, where each vertex contains a grayscale value. We assume that the cubes are triangulated conformally in order to get simplicial input -the argument also works, however, for the case of cubical cells. We assign function values inductively, assigning to each simplex the maximal value of its faces. Assuming that all vertex values are distinct, the lower star of vertex v is the set of all simplices which have the same function value as v. Filtering the simplices in a way that respects the order of the function values, we get a lower star filtration of the image. Now choose the lower stars as the chunks in our reduction algorithm. Note that the lower star is a subset of the star of the corresponding vertex, which is of constant size (assuming that the dimension d is constant). Therefore, the complexity bound (1) reduces to O(n + gn + g 3 ) = O(gn + g 3 ).
Note that global columns with large index persistence might still have very small, or even zero, persistence with respect to the function values, for instance in the presence of a flat region in the image where many vertices have similar values.

Experiments
We implemented two versions of the algorithm presented in Section 4: a sequential and a parallel version (using OPENMP), in which the first two steps of the algorithm are performed simultaneously on each chunk and on each global column, respectively. In both cases, we use √ n as the chunk size. For a fair comparison, we also re-implemented the algorithms introduced in [4,8] in the same framework, that means, using the same data representations and low-level operations such as column additions. Our implementation is publicly available in our new PHAT library for computing persistence homology, available at http://phat.googlecode.com/. Additionally, we compare to the memory efficient algorithm [10] based on discrete Morse theory [9] and to the implementation of the persistent cohomology algorithm [5] found in the DIONYSUS library [15].
To find out how these algorithms behave in practice, we apply them to five representative data sets. The first three are 3D image data sets with a resolution of 128 3 . The first of these is given by a Fourier sum with random coefficients and is representative of smooth data. The second is uniform noise. The third is the sum of the first two and represents largescale structures with some small-scale noise. These data sets are illustrated in Figure 1 by an isosurface. a) b) c) Figure 1: A single isosurface of the representative data sets used to benchmark and compare our algorithm: a) a smooth data set, b) uniformly distributed noise, c) the sum of a) and b).
In addition to the lower star filtrations of these image data sets, we also consider an alpha shape filtration defined by 10000 samples of a torus embedded in R 3 , and the 4-skeleton of the Rips filtration given by 50 points randomly chosen from the Mumford data set [11].
As pointed out in [5], the pairs of persistent cohomology are the same as those of persistent homology. We therefore also applied all algorithms to the corresponding cochain filtration, given in matrix form by transposing the boundary matrix and reversing the order of the columns; this operation is denoted by (·) ⊥ . When reducing such a coboundary matrix with the clearing optimization, columns are processed in order of increasing instead of decreasing dimension. Table 1 contains the running times of the above algorithms applied to filtrations of these five data sets run on a PC with two Intel Xeon E5645 CPUs. We can observe a huge speed-up caused by the clearing optimization as already reported in [4]. We can see that our chunk algorithm performs slightly worse than the one of [4] when executed sequentially, but faster when parallelized. We also observe that standard, twist, and chunk generally behave worse when computing persistent cohomology, except for the Rips filtration of the Mumford data set, where the converse is true.

Dataset
n · 10 −6 std. [8] twist [4] cohom. [ Table 1: Running time comparison of various persistent homology algorithms applied to the data sets described in Section 6. The last three columns contain information of the algorithm presented in this paper: the fraction of global columns g/n, and the running times using one and twelve cores, respectively.

Conclusion and Outlook
We have presented an algorithm for persistent homology that includes two simple optimization techniques into the reduction algorithm. It can be fully parallelized, except for the reduction of compressed global columns, whose number is often small. Besides our asymptotic complexity bounds, which give a detailed dependence on the parameters of the algorithm, our experiments show that significant speed-ups can be achieved through parallelized persistence computation. Similar observations have been made recently by Lewis and Zomorodian [12] for the computation of (non-persistent) homology; see also [13]. We plan a more extensive discussion of the practical effects of our optimizations and parallelization in an extended version of this paper.


extfieldsofstudy:
Engineering
Computer Science

text:
Who watches the watchers: Validating the ProB Validation Tool

Over the years, ProB has moved from a tool that complemented proving, to a development environment that is now sometimes used instead of proving for applications, such as exhaustive model checking or data validation. This has led to much more stringent requirements on the integrity of ProB. In this paper we present a summary of our validation efforts for ProB, in particular within the context of the norm EN 50128 and safety critical applications in the railway domain.

Introduction
Over the years, PROB has moved from being a tool that complemented proving, to being a tool that is now sometimes used instead of proving for certain applications. PROB can be used as a plug-in for proof-centric development environmants like Rodin or Atelier-B, but it can also be used as a standalone development environment that focuses on validation or verification using animation and model checking. For example, PROB is sometimes used to exhaustively model check B specifications which are not easily amenable to proving (e.g., where it is not easy to find inductive invariants). In 2009 PROB started to be used for data validation within Siemens [8], replacing custom proof procedures. In this application, complicated properties had to be validated on concrete data values for safety critical railway applications. As this turned out to be highly successful, PROB is now being used for data validation by a variety of other companies (Alstom, ClearSy, Systerel) and for a variety of industrial end users (e.g., RATP).
All this means that we can now longer rely on the integrity of the B provers in case a bug in PROB prevents the detection of an error. As such, validation of the PROB tool became very important and even mandatory. Indeed, according to the railway norm EN 50128 [4] [Sect. 3.1], PROB is now being used as a tool of class T2, which mandates a certain amount of documentation and validation. Our long term goal is to enable PROB to also be used as a tool of class T3, i.e., moving from a tool of class T2 that "supports the test or verification of the design or executable code, where errors in the tool can fail to reveal defects but cannot directly create errors in the executable software" [4] to a tool that "generates outputs which can directly or indirectly contribute to the executable code (including data) of the safety related system" [4]. In addition to the validation efforts described in this paper, T3 would require a secondary toolchain as described in "Checking Computations of Formal Method Tools -A Secondary Toolchain for PROB". 1 In this paper we outline the various validation efforts that have been conducted in that light, and present a summary of our experiences.

Source Code
The kernel of PROB is developed in Prolog, but many other programming languages come into play. More precisely, the source code of PROB contains more than 45,000 lines of Prolog, more than 15,000 lines of Tcl/Tk, more than 5,000 lines of C (for LTL and symmetry reduction), 1,951 lines of grammar and more than 10,000 lines of Java that are expanded by the SableCC compiler generator [6] into about 91,000 lines of Java for our parsers. In addition, there are slightly more than 5,000 lines of Haskell code for the CSP parser and about 70,000 lines of Java code for the Rodin [2] plugin. These are the statistics for version 1.3.7-beta9 of PROB.
We use Git for our source code repositories, managing a history that traces back to 2005. The repository has been ported to different SCMs multiple times, starting from an initial CVS repository, to SVN and to GIT. In addition to our main Git repositories we manage several SVN repositories for examples, test cases and documentation. These have not been moved to Git for practical reasons.
We have two versions of the PROB binary: • the ProB Tcl/Tk version, which has a GUI and • the probcli binary. It does not depend on Tcl/Tk and can be used from the command-line, from within Makefiles and via a socket interface which listens to requests. The latter is used for integrating PROB into the Rodin Eclipse platform.
Both binaries are used for industrial applications such as data validation; the Tcl/Tk version has the additional benefit of providing a predicate debugger and graphical visualization of formulas.

Important Components
A graphical representation of the PROB architecture can be found in Figure 1. The most important components are the following: • The PROB kernel deals with B's datatypes (integers, booleans, pairs, sets, ...) and the various operations on them (arithmetic, set operations, ...).
• The B interpreter interprets the structure of a B machine to compute the initial states, the enabled operations, the effect of performing an operation, etc. The B interpreter calls the PROB kernel (via the kernel mappings interface file).
• The B Type Checker reads in the parsed abstract syntax tree (AST), type checks it and performs certain pre-processing.
• the B Parser reads in a series of B files and produces an abstract syntax tree (AST) in Prolog form (written into a file ending with ".prob").
Any error in those components is likely to result in a potential error during data validation or model checking, and as such these components need to be tested very thoroughly.

ProB kernel
The kernel supports the underlying data structures and operators of the B language: arithmetic, sets, relations, functions and sequences. As such, the specification for the ProB kernel is the reference book  by Abrial [1], in particular Chapters 2 and 3. As additional reference, the Atelier-B handbook [5] has been used, for those operators not covered in [1].
Many of the laws from [1] are tested below, in various forms, e.g., instances are used as unit tests or more systematically, in Section 3.4.2 we ensure that ProB can find no counter-example to those laws. For example, on page 124 we can find the law union(SS) = {x| ∃ y . (y ∈ SS ∧ x∈y)}, which is added as one of the laws in the machine checking set laws. The particular example from that page u = {{0,5,2,4}, {2,4,5}, {2,1,7,5}} with union(u) = {0,1,2,4,5,7} -is added as a regression test to the machine checking explicit computations of laws.
Several mistakes in the B-Book were found in this process; e.g., the first definition of override on page 80 uses p instead of r. In the second definition, it should be (a,b) instead of (x,y). On page 102, the property involving ran(p; q) = ran(p) is false, it should be ran(q).

B Interpreter
The B interpreter evaluates B predicates, expressions 2 and substitutions. As such, again the book by Abrial [1] is the main specification.
In particular, Chapter 1.2 is used for predicates, and Chapter 4 is used for substitutions. As additional reference, the Atelier-B handbook [5] has been used, for those aspects not covered in [1].
For substitutions, the definitions in [1] are in form of rules for the weakest pre-condition. PROB uses an interpreter which, given a predecessor state, computes all possible successor states. (An alternate view is that PROB finds solutions for the before-after predicate of a substitution.) As such, the weakest pre-condition definitions do not apply directly, and PROB is validated in other ways.

Parser, Type Checker and Syntax Tree Manipulations
For the parser, the syntax of Atelier-B [5] is the reference specification. The type checking basics are specified in [1]. For the visibility rules, the tables in [5] have been used as reference.

Building and Testing through Continuous Integration
PROB contains unit tests, integration and regression tests as well as self-model check tests for mathematical laws (see below). All of these tests are run automatically on our continuous integration platform "Jenkins". 3 The tests are run on all platforms for which we ship PROB. When a test fails, an email is sent automatically to the PROB development team. The status of the tests can be monitored on the web pages of our Jenkins instance (see Figure 2): http://cobra.cs.uni-duesseldorf.de/jenkins/. Over the years many tests have been developed. They can be classified into unit tests for the PROB kernel, parser tests, type checker tests, and system tests exercising the whole PROB tooling infrastructure. There are two situation in which we manually add tests. The first is to add unit tests following the test-driven development process. The second is to add tests aiming to increase code coverage as we will explain in Section 4.

Unit Tests for Kernel
In addition, we have developed an automatic unit test generator, which tests the PROB kernel predicates for many different scenarios and different set representations. For example, starting from the initial Prolog call union([int (1)],[int (2)],[int(1),int(2)]), the test generator will derive 1358 individual unit tests. The various tests will use different set representations (AVL-tree or symbolic representations), will swap the order of the first two arguments as union is declared commutative, and will check various orderings in which the sets are instantiated (e.g., it could be that first the result of the union is known, then the second argument). The latter point is particularly important for the PROB kernel, which relies on co-routines. Indeed, we have to check that the kernel predicates behave correctly no matter in which order (partial) information is propagated.

Validation of the parser
The PROB parser is written in Java and has been developed using SableCC. We use about 1800 JUnit based tests to check various parts of the parser for errors. While most of them verify the correct behavior on small examples or single expressions, some employ full scale B machines. The unit tests are run independently on both Windows and Linux. Several of these tests are automatically adopted to different newline characters and encodings.
Additionally, we execute our parser on a large number of our regression test machines and pretty print the internal representation. We then parse the internal representation and pretty print it again, verifying (with diff) that we get exactly the same result (see Figure 3). This type of validation can easily be applied to a large number of B machines, and will detect if the parser omits, reorders or modifies expressions, provided that the pretty printer does not compensate errors of the parser. On the downside, the validation will only detect those errors in machines generated by the pretty printer, which may prevent us from catching errors which only appear in non-pretty printed machines, e.g., when parentheses in expressions are set incorrectly.

Validation of the type checker
For the moment we also read in a large number of our regression test machines and pretty print the internal representation, this time with explicit typing information inserted. We now run this automatically generated file through the Atelier B parser and type checker (see Figure 4). With this we test whether the typing information inferred by our tool is compatible with the Atelier B type checker. (Of course, we cannot use this approach in cases where our type checker detects a type error.) Also, as the pretty printer only uses the minimal number of parentheses, we also ensure to some extent that our parser is compatible with the Atelier B parser (see below). Again, this validation can easily be applied to a large number of B machines. More importantly, it can be systematically applied to the particular B machines that PROB validates: provided the parser and pretty printer are correct, this gives us a guarantee that the typing information for those machines is correct. The latest version of PROB Tcl/Tk has a command to cross check the typing of the internal representation with Atelier B in this manner.
Errors Detected During validation of the type checker we found several issues. (We did not detect errors in the type checker, but rather differences between our parser and the Atelier B parser bcomp.) For example, PROB accepts identifiers with single letter without warning and PROB accepts enumerated set elements that are used as identifiers while bcomp does not. We also detected that bcomp reports a lexical error ("illegal token |-") if the vertical bar (|) of a lambda abstraction is followed directly by the minus sign. Most importantly, however, we found errors in the priority table of the english Atelier B "B Language Reference Manual 1.8.6" [5]. All in all there are 26 errors in the English reference manuals, upon which our pretty printer and BParser was based. These issues have been fixed both in our parser and in the handbook.

System Tests
These tests exercise the entire PROB tool chain, including the parser, the type checker, the PROB interpreter and the PROB kernel.

Regression Tests
PROB contains over 1100 regression tests which are made up of B models along with saved animation traces. These models are loaded, the saved animation traces are replayed, and the models are also run through the model checker. These tests have turned out to be extremely valuable in ensuring that a bug once fixed remains fixed. They are also very effective at uncovering errors in arbitrary parts of the system (e.g., the parser, type checker, the interpreter, the PROB kernel, etc.). The various features of PROB are law1 == ( dom ( f f ∪ gg ) = dom ( f f ) ∪ dom ( gg ) ) ; law2 == ( r a n ( f f ∪ gg ) = r a n ( f f ) ∪ r a n ( gg ) ) ; law3 == ( dom ( f f ∩ gg ) ⊆ dom ( f f ) ∩ dom ( gg ) ) ; law4 == ( r a n ( f f ∩ gg ) ⊆ r a n ( f f ) ∩ r a n ( gg ) ) ; . . . .

Self-Model Check with Mathematical Laws
With this approach we use PROB's model checker to check itself, in particular the PROB kernel and the B interpreter. The idea is to formulate a wide variety of mathematical laws and then use the model checker to ensure that no counterexample to these laws can be found.
Concretely, PROB now checks itself for over 500 mathematical laws. In version 1.3.3 there were laws for booleans (39 laws), arithmetic laws (40 laws), laws for sets (81 laws), relations (189 laws), functions (73 laws) and sequences (61 laws), as well as some specific laws about integer ranges (24 laws) and the various basic integer sets (7 laws). The number of laws is continuously being expanded. Figure 5 contains some of these laws about function. 4 This technique can be viewed as using the model checker to generate a very large number of tests from the mathematical laws. The self-model check has been very effective at uncovering errors in the PROB kernel and interpreter. Furthermore, the self-model checking tests rely on every component of the PROB execution environment working perfectly; otherwise a violation to a mathematical law could be found. In addition to the PROB main code, the parser, type checker, Prolog compiler, hardware and operating system all have to work perfectly. Indeed, we have identified a bug in our parser (FIN was treated the same as FIN1) using the self-model check. Furthermore, we have even uncovered two bugs in the underlying SICStus Prolog compiler using self-model check: • A bug in the AVL library (notably in the predicate avl max computing the maximum element of an AVL-tree) was found and reported to SICS. The bug was fixed in SICStus Prolog 4.0.5.
Note that these problems would not have been detected by validating or proving the code of PROB correct. It was essential to test the actual code of PROB along with the entire execution environment. The model checker together with the mathematical laws enabled this testing to be performed very effectively (see also Example 1 below).
We also applied some of these tests to other tools. For example, we found various bugs in the first implementation of the JEB animator [9], and the TLC [10] model checker when applied to TLA translations of our B laws.

Coverage
The validation techniques explained above are complemented by code coverage analysis techniques. In particular, we try to ensure that the unit tests (Section 3.1) and the self-model checks (Section 3.4.2) above cover all predicates and clauses of the PROB kernel. 5 We count a clause as covered if it has been executed successfully by at least one test case. A predicate is regarded as covered if at least one of the clauses belonging to it is covered.
Initially, we have developed our own code coverage tool for SICStus Prolog. The tool uses Prolog's term expansion facility to keep a record of which program points are covered. As of version 4.2, SICStus Prolog now provides code coverage profiling, along with visualizations in Eclipse and Emacs. The coverage is slightly less fine-grained than our original implementation, but it is much faster (as the profiling is directly supported by the Prolog abstract machine) and supports all of Prolog's features out of the box. We have thus decided to switch to SICStus' code coverage profiler, but have implemented various backends to output the information in HTML, LaTeX, or other formats (see http://nightly.cobra.cs.uni-duesseldorf.de/coverage/). One of the backends features a colorized version of our source code, highlighting the coverage status. With this backend, it is possible to inspect the source code itself, finding out exactly which parts of the code are covered. This can be used to identify poorly tested source code and helps in adding missing test cases.
The code coverage statistics are also systematically computed by our continuous integration server (see Section 3). As of version 1.3.5 of PROB, we now also keep track of clauses which should not be covered as they relate to internal errors: • the error manager module was extended to distinguish internal and normal errors. A normal error is, for example, a division by zero in a B machine. An internal error would be a typing error or unexpected failure within a piece of PROB Prolog code. Normal errors can be covered by tests; internal errors must not be coverable.
• a Prolog term expander keeps track of clauses and execution paths which relate to internal errors, and excludes those from the coverage statistics. So far, the code coverage analysis has been very useful in extending our unit tests and mathematical laws. As an example, 96.5 % of the clauses of the kernel mappings module are now covered by the unit tests and the mathematical laws. The only uncovered clauses relate to the Z compaction operator, and one error condition that was not triggered by the tests. In summary, the code coverage has helped us write better tests and has allowed us to uncover a few undetected errors in the kernel.
However, even fully covered code might still contain errors. Executing all parts of the source code does not necessarily mean that all paths were covered. There might, for instance, be bugs that only occur if certain combinations of branches are executed. Thus we can not rely solely on the coverage reports to judge the effectiveness of our test suite. We experimented with artificially introduced bugs to evaluate how certain classes of programming errors can be found. The following section will give a small overview.

Judging Effectiveness of Tests
In Section 3 we already described the effectiveness of the current testing procedure at uncovering various errors in PROB and in the underlying compiler. We now check the effectiveness for artificially introduced errors.

Validating ProB
The tests in this subsection were conducted on PROB version 1 Passed As you can see, all errors were detected and the mathematical laws checks detected most of the errors. However, as checking the mathematical laws is based on searching for counter-examples, there is a certain class of errors that can not be found. If due to a bug the search space is pruned, mathematical law checking will obviously not be able to find more counter-examples within the pruned area. Therefore the mathematical law checks do not make the unit tests obsolete.
In the future, we plan to conduct a more extensive study about the effectiveness of our tests.

Performance and Codespeed
We run several benchmarks for PROB on a daily basis. The results are stored using the codespeed web application 6 which allows us to compare selected revisions on selected benchmarks. Furthermore, an unusual increase or decrease in performance is reported right on the front page. Strictly speaking this does not add value to our test infrastructure, as the benchmarks do not report error or success. However, in several occasions the existence of a bug in PROB lead to strange speed-ups or slowdowns, for instance, if a bug causes PROB to check every possible value of an existential quantification instead of stopping after finding a solution. Therefore, monitoring the overall speed of our applications is a useful addition to testing. In order to identify the particular revision causing unexpected benchmark results, we have found git bisect to be useful. It allows us to mark certain revisions as good (i.e. the last benchmarked revision that returned the expected results) and certain other revisions as bad (i.e. the revision showing strange behavior). Git then performs an interactive binary search, checking out revisions in between until the one responsible has been found.
The process of comparing certain revisions of PROB has been partly automated by several scripts that control checkout, benchmarking and result reporting.

Static Analysis
Some bugs are hard to detect by testing but easy to detect with static analysis. For instance, we had some calls to error handling predicates that were not used correctly, e.g., they were called with the wrong arity. In most cases these calls happened in source code locations that are unreachable unless a previous bug occurs in the underlying execution engine. There were also bugs in the reachable parts of ProB that were not well tested, i.e., those parts that are of a more experimental nature.
Prolog is not statically typed and therefore we cannot rely on a typechecker to identify incorrect calls to predicates, however it is very easy to implement an analysis for these sort of problems in Prolog by augmenting the interpreter with a term expander. Term expansion is very similar to macro expansion in Lisp. A term expander hooks into the loading process of a Prolog file and allows modification of a term while it is being loaded. In our case, we did not modify the terms, but extract information from it. Hence, the expanded code is equal to the code used in production. This allows us to relate all detected bugs directly to our code. For instance, if we analyze the term :-dynamic foo/1 in a module bar, we store the information that a predicate foo of arity 1 exists in the module bar and that it is a dynamic predicate. We also extract information about predicate and fact definitions, mode declarations, multifile and meta annotations, exported and imported predicates, blocking declarations for co-routines and most importantly calls. Strictly speaking this is not traditional static analysis, because Prolog directives such as use module are executed while loading the Prolog program. But the program is not run either, so we think it is appropriate to classify the approach as a static analysis. To run the analyzer, we first load our term expander and then the entry point module of ProB. During loading the module, the term expander inspects all terms that are being loaded and stores the extracted information in a Prolog database. The analyzer automatically detects calls where no definition is known at compile time. These calls are potential bugs. Using the tool, we were able to find 17 missing import statements, 9 cases where predicates are used with the wrong arity, 5 cases where the wrong predicate was called and 5 missing dynamic declarations. The missing imports did not actually cause wrong behavior at runtime because of a Sicstus Prolog implementation detail. As long as a module is loaded (even if it was loaded by any other module ), a predicate can be called using a fully qualified form module:predicate(...) at runtime. This is rather fragile because if for some reasons the module is not loaded the calls will break. Each module should declare its dependencies explicitly. The calls to predicates with wrong arity were also not critical because they were in dead code (if the Prolog engine works properly).
Three cases where the wrong predicate was called were actual bugs in some experimental features of ProB. One feature was removed, the other two were fixed. The other two cases were dead code. The missing dynamic declarations are not dangerous in the current version of Sisctus Prolog, but there is no guarantee that this won't cause actual bugs in future versions. The problems we found using the tool are summarized in the following table:

Problem Type
Fragile code/Bad practice Bug Dead code Import missing 17 0 0 Wrong arity 0 0 9 Wrong predicate 0 3 2 Missing dynamic delaration 5 0 0 Because Prolog is a dynamic language, finding dead code is harder than in a static language such as Java. The compiler cannot know if a call is produced dynamically at runtime. This can easily lead to dead code, if we change our code. For instance, if we replace a call to predicate p by a call to predicate q, is p now dead code or is it used somewhere else? Using our analysis we get a list of predicates that are not being called in a static way. This clearly produces false positives, but it is possible to manually inspect the candidates. We think that we can further improve the analysis by detecting typical dynamic usage pattern such as using a predicate in a higher order predicate, e.g. maplist.
In addition to the static analysis, we are currently working on a tool that can be used to inspect the ProB source code. The tool should contain a code browser that allows to directly navigate beween caller and callee, search for predicates, filter according to specific criteria (e.g., parts of the name). This tool should also incorporate the documentation. We want to display certain properties of the source code (e.g., percentage of predicates that are documented or percentage of predicates with a mode declaration) to encourage the developers to increase the quality of the documentation and code.
Our custom tools are accompanied by Spider, an Eclipse based IDE specifically developed by SIC-Stus for Prolog development. It features several static analyses, including calling mode and checks for potential non-determinacy and potential non-termination. In addition, there is an option to check for common mistakes and bad practice.

Dynamic Runtime Analysis
In this section we examine various ways we ensure consistency of PROB's output at runtime.

Double chain
One way to dramatically increase the confidence in PROB's output is to double check the result with a separately developed tool. Thus far, this has been done with the following three tools: 1. PredicateB. Together with the company ClearSy we have performed Data Validation tasks for Alstom [7]. Here, the Java based tool PredicateB was used to double check PROB's output. Quite often, PredicateB took considerably more time than PROB to check the properties. Some properties needed to be rewritten to be suitable for PredicateB (whose constraint solving capabilities are much weaker). Overall, this approach was very successful and led to the development of the DTVT (Data Table Validation Tool) which integrates PROB and PredicateB.
2. Ovado. In this work [3], the company Systerel is using the tool Ovado as primary tool, and uses PROB to double check Ovado's output.
3. PyB. Finally, in recent work we are developing an independent tool written in Python to check B predicates. The idea is to feed PROB's output into PyB (in a more integrated fashion, e.g., passing witnesses for existential predicates) and to let PyB check the output. An independent companion paper about PyB has been submitted to this workshop.

Run Time Checking
The Prolog code contains a monitoring module which -when turned on -will check pre-and postconditions of certain predicate calls and also detect unexpected failures. Many kernel predicates also check for unexpected arguments. All of this overcomes the fact that Prolog has no static typing to some extent. So far we have not systematically relied on this feature; but we could turn the run time checking on during property validation. The only drawback is that the runtimes will be increased.

Positive and Negative Evaluation
As already mentioned, all properties and assertions were checked twice, both positively and negatively (see Figure 6). Indeed, PROB has two Prolog predicates to evaluate B predicates: one positive version which will succeed and enumerate solutions if the predicate is true and a negative version, which will succeed if the predicate is false and then enumerate solutions to the negation of the predicate. The reason for the existence of these two Prolog predicates is that Prolog's built-in negation is generally unsound and cannot be used to enumerate solutions in case of failure. For example, given the Prolog rule p(eq(X,Y)) :-X=Y. the query not(p(eq(X,0))) would fail, i.e., one could erroneously conclude that there is no value X which is different from 0. In PROB, we have a dedicated predicate for the negation, which will suspend until it can determine the result correctly. For the above example, one could write not p(eq(X,Y)) :-dif(X,Y).
timeout P is (probably) true false timeout P is false or undefined timeout true P is (probably) false timeout false P is true or undefined timeout timeout P is unknown With these two predicates we can uncover undefined predicates: if for a given B predicate both the positive and negative Prolog predicates fail then the formula is undefined. For example, the property x = 2/y & y = x-x over the constants x and y would be detected as being undefined, and would be visualized by our graphical formula viewer as in Figure 7 (yellow and orange parts are undefined).
In the context of validation, this approach has another advantage: for a formula to be classified as true the positive Prolog predicate must succeed and the negative Prolog predicate must fail, introducing a certain amount of redundancy (admittedly with common error modes). In fact, if both the positive and negative Prolog predicates would succeed for a particular B predicate then a bug in PROB would have been uncovered. If both fail, then either the B predicate is undefined or we have a bug in PROB. 7 This validation aspect can detect errors in the predicate evaluation parts of PROB, i.e., the treatment of the Boolean connectives ∨, ∧, ⇒, ¬, ⇔, quantification ∀, ∃, and the various predicate operators such as ∈, ∈, =, =, <, ... This redundancy can not detect bugs inside expressions (e.g., +, −, ...) or substitutions (but the other validation aspects mentioned above can). Example 1. We want to illustrate the error detection of double evaluation by a simple example. Take the following very simple B machine, with two assertions which should obviously be true.

MACHINE D o u b l e E v a l u a t i o n T e s t SETS ID={ aa , bb
Say that we introduce an error in the part of the source code that deals with testing membership of elements in a singleton set. An element is a member of a singleton set if it is equal to the single element. We replaced this test by its negation, and thus introduced an error into the PROB kernel. When analyzing the assertions, we now get the following result, which tells us that there is a bug in the ProB kernel, because the assertion stating that iv is not a member of the set {bb} is both true and false at the same time: r u e f a l s e Also note that the self-model checking described in Section 3.4.2 immediately provides a counter example (SS={el1}, TT=∅) for several of the mathematical laws for sets, for example: { xx | xx ∈ SS ∨ xx ∈ TT } = SS ∪ TT == f a l s e

Future Work and Discussion
In the longer term we want to further increase the confidence in PROB's output so that it can eventually be used as a tool of class T3 within the railway norm EN 50128.
In addition to the development of the double chain PyB, there are several techniques that might help to increase the confidence in our primary validation tool even more. One of these techniques is fuzzing or robustness testing. The idea behind fuzzing would be to generate random input data and feed it to PROB. As mentioned before, test coverage of error cases is hard to archive. Random inputs of data could help us to identify cases in which some part of our tool chain crashes instead of reaching a controlled error state.
We could also further refine the results of the coverage analysis mentioned in Section 4. Instead of just observing which predicates were executed we could strive for coverage of possible execution paths or the evaluation of each branch. However, this would again increase the complexity of the analysis.
Aside from improving the static checks strengthening the run time checks would be another reasonable approach. With the increased usage of our code documentation tools more annotations will be added to the source code. The correctness of several of these annotations, could be enforced at runtime. One example, is the Prolog mode declaration which specifies which arguments are input arguments or output arguments respectively.

Conclusion
In summary, the following points increase our confidence in the proper working for PROB in general and for data validation in particular: 1. repeated successful application on case studies (e.g., by Siemens on the lines of San Juan, Paris Line 1, Sao Paolo, Barcelona Line 9, ... [8]), discovering all known problems and in some cases even discovering previously unknown problems in the data properties.
2. extensive regression, unit-tests, automated unit-test generation, and validation via mathematical laws and model checking (see Section 3). These tests have even discovered several errors in the underlying Prolog compiler.
3. coverage analysis of the testing, ensuring that all critical parts are covered by the tests; this has enabled us to detect a few more issues (and has prompted us to develop the automated unit-test generator).
4. validation of the parser and type-checker via pretty-printing and cross-checking with the Atelier-B parser (bcomp) (see Section 3.2). The internal representation of PROB for a particular model can be fed at runtime to bcomp, for additional validation that the type inference was properly conducted.
5. double evaluation of predicates (see Section 7.3), to catch undefined predicates as well as errors in the PROB kernel and interpreter relating to predicates.