Xavier Cavin

Olivier Demengeon

Shift-Based Parallel Image Compositing on InfiniBand TM Fat-Trees

Keywords: Descriptors (according to ACM CCS): I, 3, 1 [Computer Graphics]: Hardware Architecture-Parallel processing I, 3, 2 [Computer Graphics]: Graphics Systems-Distributed/network graphics I, 3, 8 [Computer Graphics]: Applications-Image compositing

Parallel image compositing has been widely studied over the past 20 years, as this is one, if not the most, crucial element in the implementation of a scalable parallel rendering system. Many algorithms have been proposed and implemented on a large variety of supercomputers. Among the existing supercomputers, InfiniBand TM (IB) PC clusters, and their associated fat-tree topology, are clearly becoming the dominant architecture, as they provide the scalability, high bandwidth and low latency required by the most demanding parallel applications. Surprisingly, very few efforts have been devoted to the implementation and performance evaluation of parallel image compositing algorithms on this kind of architecture. We propose in this paper a new parallel image compositing algorithm, called Shift-Based, relying on a well-known communication pattern called shift permutation. Indeed, shift permutation is one of the possible ways to get the maximum cross bisectional bandwidth provided by an IB fat-tree cluster. We show that our Shift-Based algorithm scales on any number of processing nodes (with peak performance on specific counts), allows overlapping communications with computations and exhibits contentionfree network communications. This is demonstrated with the image compositing of very high resolution images at interactive frame rates.

Introduction

Parallel Rendering

Both in the industry and in the academic worlds, larger and larger numerical datasets are being produced on a daily basis, whether they are large Computer Aided Design (CAD) assemblies, results of numerical simulations or 3D scanned information, just to cite a few. Unless a considerable preprocessing can be applied to the data, a single computer cannot accomplish real-time rendering of these data. In most cases, high-end parallel visualization solutions must be designed to tackle these ever-increasing amounts of data, sometimes on-the-fly as they are being produced.

Sort-last parallel rendering, as introduced by Molnar et al. [START_REF] Molnar | A sorting classification of parallel rendering[END_REF], is an efficient technique to render very large datasets in parallel. In this kind of parallel applications, the dataset to be visualized is first partitioned and distributed across multiple processors. Each processor renders its subset of the dataset independently using its local resources; this stage does not imply any communication between the processors and scales perfectly. Then the computed images are composed together by a parallel image compositing algorithm and the final image is delivered. The parallel image compositing stage requires a lot of inter-processor communication, and may become the bottleneck of a sort-last parallel pipeline if it is not designed with sufficient care.

Since the early 1990's, parallel image compositing has been widely studied, and numerous algorithms have been developed and experimented on many parallel computing architectures, ranging from supercomputers to PC clusters. These image compositing algorithms are commonly classified into three categories: Direct Send [START_REF] Neumann | Communication costs for parallel volume-rendering algorithms[END_REF] fies the Direct Send and Binary Swap and enables numerous other configurations.

The 2-3 Swap and Radix-k algorithms have been successfully deployed on very large supercomputer configurations. Yu et al. report compositing times of 2 seconds per image for a 16 Mega-pixels (16 bytes per pixel) resolution with 1024 processors of a Cray XT4 system connected to a SeaStar 2 router. Peterka et.al report compositing times of 0.7 seconds per image for an 8 Mega-pixels (16 bytes per pixel) resolution with up to 34,816 processes on an IBM Blue Gene/P.

Multistage Interconnection Networks

A large majority of high performance computing clusters today rely on the central switch network architecture. This network architecture considerably simplifies the work of the programmer, since each node of the cluster can communicate with every other node. In the particular case where the central switch is a true crossbar, the network has almost ideal properties: constant latency between all pairs of nodes and full bisection bandwidth (any half of the nodes can simultaneously communicate with the other half at full speed).

Most parallel applications and communication libraries typically treat central switch network architectures as if they were true crossbar switches, and do not take into account the actual internal characteristics of the switch. However, practical central switches are generally implemented as multistage interconnection network (MIN), and do not truly provide the latency and full bisection bandwidth characteristics of crossbars, as shown in [HSL08]. In particular, the routing policies, needed to route connections in a MIN, can greatly affect the performance of a given communication pattern. The fat-tree topology is a class of MIN and is one of the most commonly used topologies for high performance computing clusters. In particular for InfiniBand TM (IB) based clusters, this topology is the dominating one. This includes large installations such as CEA Tera-100, LANL Roadrunner, TACC Ranger or FZJ-JuRoPA (Source: http://www.top500.org/).

For IB fat-trees, as with most other topologies, the routing algorithm is crucial for efficient use of the underlying topology. IB provides a deterministic static-destination routing that is computed offline and loaded into the different switches of the MIN. Several options are offered to the system administrator to compute the routing tables, depending on the physical network architecture of the cluster. A given communication pattern can either achieve satisfactory performance or on the contrary can lead to heavy contention for different routing strategies.

Our Contributions

In this paper, we start in Section 2 with an overview of fattree IB clusters, which represent a significant part of today's high-performance supercomputers, and with a presentation of existing parallel image compositing algorithms. Surprisingly, we have not found any study of these algorithms on this ever-growing class of supercomputers. In Section 3, we introduce a new parallel image compositing algorithm called Shift-Based, relying on a very well-known parallel communication pattern called shift permutation. Indeed, a lot of efforts have been dedicated to efficiently support this communication pattern on fat-tree IB clusters. We continue in Section 4 with an analysis of the theoretical communication and computation costs of our new algorithm, and we compare them to those of the existing parallel image compositing algorithms, in an ideal environment. Then, we show how network contention, appearing in real-life environments, can negatively affect the expected performance of these algorithms. In Section 5, we experiment our algorithm on a 144node fat-tree IB cluster and we show how it maximizes the utilization of the bandwidth made available by the IB cluster. We demonstrate that our algorithm scales very well on any number of processors (with peak performance on specific counts) and performs at least twice better as compared to existing algorithms. We conclude and present future work in Section 6, anticipating similar performance results at larger scales.

Background

Fat-Tree IB Clusters

The dominant architecture of scalable high-performance supercomputers has clearly become the clustered architecture; the latest (36th) edition of the TOP500 list reports 414 clustered systems (82.8 %). For these clusters, as well as for other systems, a scalable, high bandwidth and low latency network is a key. The InfiniBand TM (IB) architecture standard [START_REF]IBTA -INFINIBAND TRADE ASSOCIATION: IBTA specification 1[END_REF] meets all these requirements and is the interconnection network of choice for such systems: the same TOP500 list reports 213 IB based systems (42.6 %) in total, including 178 IB clusters.

Large IB clusters are commonly built with multiple fixedarity crossbar switches (up to 36x 40 Gbps ports) and hostchannel adapters (HCAs) providing up to 40 Gbps links. The fat-tree topology is the most commonly used to interconnect the nodes and IB switches, providing varying ratios of crossbisectional bandwidth (CBB). Closely related to the chosen topology, the routing algorithm, used to route the data between switches from one node to the other, is of particular importance when one wants to maximize the utilization of the network.

We focus in this paper on large scale fat-tree IB clusters, as they are sufficiently representative of today's supercomputers, and we will review their main characteristics in the remaining of this Section. However, a similar analysis could be conducted on the other types of interconnection networks and supercomputers architectures.

IB Clusters

The basic building blocks of large IB clusters are integrated non-blocking switch elements (i.e. crossbars) with a relatively low number of ports. As an example, QDR (40 Gbps) IB switches come in 8, 18 and 36 ports configurations. Crossbar switches exhibit almost ideal network properties: constant latency between all pairs of endpoints and full bisection bandwidth (any half of the nodes can simultaneously communicate with the other half at full speed, i.e. without contention).

The IB standard does not impose any particular network topology. In other words, switches can be connected together in any arbitrary way. IB switches use oblivious destination based distributed routing, meaning that a switch routes a packet out a particular port based solely on the destination address of the packet. This decision is determined via a simple static routing table (Random or Linear Forwarding Table) that defines which port leads to which endpoint. On startup, or after a topology change, a central entity called the Subnet Manager (SM) discovers the network topology using special packets computes the routing table for each switch and uploads them into the switches. As this operation requires network downtime, it cannot be performed too frequently.

The current implementation of OpenSM offers five routing engines: Min Hop, UPDN Unicast, Fat-tree Unicast [START_REF] Zahavi | Optimized infiniband fat-tree routing for shift all-to-all communication patterns[END_REF], LASH Unicast, DOR Unicast. Additionally, OpenSM also supports a file method which can load routes from a table. The interested reader can refer to the OpenSM documentation for more details about these methods.

Fat-Tree Topology

The fat-tree topology is an indirect interconnection network based on a complete binary tree [Lei85]. The processors of a fat-tree are located at the leaves of a single-rooted complete binary tree, and the internal nodes are crossbar switches. Unlike traditional trees in computer science, fat-trees are more like real trees, because they get thicker near the root. Going up the fat-tree, the number of wires connecting a node with its father increases. The rate of growth influences the size and cost of the hardware, but also the communication bandwidth it can support. As the arity of the internal switches of the fat-tree increases as we get closer to the root, the physical implementation of this topology is impractical.

To overcome this limitation, the concept of a multi-rooted topology was introduced. The non-blocking nature of the single rooted tree is replaced by a weaker attribute: given a permutation of source and destination pairs, the routing on the tree can be made non-blocking (i.e. rearrangeably nonblocking). Several families of fat-trees have been introduced. The k-ary n-tree [START_REF] Petrini | k-ary n-trees: high performance networks for massively parallel architectures[END_REF] is the basic type of tree built out of switches with an equal number of connections going up or down the tree. A k-ary n-tree is composed of k n processing nodes and nk n-1 k × k constant arity commutation switches.

The k-ary n-tree topology has been extended with the introduction of Generalized Fat-Trees (GTF) which allow for a different number of up and down connections [START_REF] Ohring | On generalized fat trees[END_REF]. Extended Generalized Fat-Trees (XGFT) further extends the possible topologies allowing for different number of connections at each level. XGFTs cannot capture existing real life topologies, since they allow only a single connection between switches: the resulting family does not preserve cross bisectional bandwidth (the bandwidth towards the top of the tree is not equal to that of the leaves).

The latest topology introduced to represent fat-trees is the Parallel ports Generalized Fat Tree (PGFT) and its sub-classification into Real Life Fat-Tree (RLFT) [START_REF] Zahavi | D-Mod-K routing providing non-blocking traffic for shift permutations on real life fat trees[END_REF]. PGFT allows multiple connections between the switches and support a large variety of topologies, including some unpractical to build. RLFT is a subclass of PGFT with several restrictions (detailed in [START_REF] Zahavi | D-Mod-K routing providing non-blocking traffic for shift permutations on real life fat trees[END_REF]), and it encompasses the majority of real life configurations.

Hot Spots in Fat-Tree IB Clusters

Fat-tree IB clusters with full cross-bisectional bandwidth have been measured to provide an effective bisection bandwidth in the range of 55-60 % for MPI communication patterns [HSL08]. This performance degradation is due to cases where multiple sources congest a particular network link, creating a hot spot. Indeed, even if enough physical links are available to support any communication pattern, the routing mechanism might oversubscribe some physical links while others remain idle. The probability of interconnect hot spots increases with the size of cluster and might degrade the latency and effective bandwidth experienced by MPI operations.

One possible solution to the hot-spot problem is the IB Lid Mask Control (LMC) mechanism. The LMC mechanism assigns multiple Local Identifiers (LIDs) to hosts and enables multiple routes for every pair of peers. However, for multiple reasons discussed in [HSL08], LMC-based routing strategies can be seen as counter-productive at large scale.

In this paper, we will focus on single-path static routing. Consequently, the chosen routing engine must be taken into account when designing a parallel algorithm in order to avoid hot spots and benefit from most of the theoretical cross-bisectional bandwidth.

Non-blocking Shift Permutations

Given a set of p nodes, shift permutation is defined as the set of source/destination pairs such that processing node P i sends data to processing node P j :

{(P i , P j )| j = (i + S)mod p} i=p-1 i=0 (1)
for S ∈ {1 . . . p -1}. Shift-based communication involves p -1 communication steps: at each step S, processing node P i sends data to processing node P (i+S)mod p and receives data from processing node P (i-S+p)mod p . c The Eurographics Association 2012.

Shift-based communication enables all nodes to send and receive data simultaneously in all steps and ensures a regular and predictive communication pattern. This particular pattern has been shown to be representative for scientific applications [START_REF] Kerbyson | Infiniband routing table optimizations for scientific applications[END_REF] and has been widely studied, on multiple network topologies including IB fat-trees.

In particular, the Fat-tree Unicast routing engine available in the latest version of OpenSM provides routing tables for contention free shift communication pattern on var-ary ntrees [START_REF] Zahavi | Optimized infiniband fat-tree routing for shift all-to-all communication patterns[END_REF], a variant of k-ary n-trees. For more complex, real life fat-trees (i.e. RLFT), d-mod-k routing has recently been proven to provide non-blocking traffic for the same shift-permutations [START_REF] Zahavi | D-Mod-K routing providing non-blocking traffic for shift permutations on real life fat trees[END_REF].

In this paper, we will focus on var-ary n-trees IB clusters with Fat-tree Unicast routing enabled. However, we can anticipate similar results with d-mod-k routing on RLFT. Our goal will be to design a parallel image compositing algorithm based on the shift permutation pattern, in order to benefit from the full available cross-bisectional bandwidth.

Parallel Image Compositing Algorithms

In this Section, we will review the existing parallel image compositing algorithms. We will assume p processing nodes. Each processing node P i owns a vector x i (the local image) of n components (the pixels). The goal of a parallel image compositing algorithm is to blend all x i with all x j using a component wise linear binary operator.

Direct Send

The Direct Send algorithm [START_REF] Neumann | Communication costs for parallel volume-rendering algorithms[END_REF] is the most straightforward parallel image compositing algorithm to implement. Each node P i is responsible of compositing the (n/p) pixels ranging from i(n/p) to (i + 1)(n/p) -1. The Direct Send algorithm requires a single stage of communication: during this stage, in the worst case, each node P i sends (p -1) messages of (n/p) pixels, one to each other node, and receives (p -1) messages of (n/p) pixels, one from each other node.

The main advantage of the Direct Send is its flexibility, since it can accommodate any number of processing nodes. It is very easy to implement and allows overlapping computations (i.e. pixels blending) with communications. However, it involves multiple processors sending messages to the same processor at the same time in an unpredictable and non-deterministic communication pattern.

Binary Swap and 2-3 Swap

The Binary Swap algorithm [MPHK94] is based on a binary tree compositing strategy which keeps every node busy in all stages of the process. It takes log 2 (p) stages to complete, where p must be a power of two (i.e. p = 2 k ) to fully exploit parallelism.

At each stage i (0 < i < log 2 (p) + 1) of the algorithm, each node sends one message of n/2 i pixels to a pair node, and receives one message of n/2 i pixels from the same pair node.

The Binary Swap algorithm exhibits a regular and predictable communication pattern. However, it does only support powers of two numbers of processors and its synchronous nature does not permit to overlap computations with communications (except on a Gigabit Ethernet network with a strongly streamed communication approach [START_REF] Cavin | Cots cluster-based sort-last rendering: performance evaluation and pipelined implementation[END_REF]).

The 2-3 Swap algorithm [YWM08] extends the Binary Swap algorithm to non-powers of two numbers of processors, at the cost of increased latency and bandwidth requirements.

Radix-k

The Radix-k algorithm [PGR * 09] is a configurable image compositing algorithm which unifies the Binary Swap and Direct Send algorithms, and enables numerous other algorithms through appropriate choice of radices or k-values.

A Radix-k algorithm is completely specified by a vector k = [k 1 , k 2 , . . . , kr] of r values such that p = ∏ i=r i=1 k i . It consists in r communication/compositing rounds. During each round i, the p processing nodes are divided into p/k i groups of k i participants, which communicate only within a group in a Direct Send fashion. When all of the k-values are equal to 2, this gives the Binary Swap algorithm. When k = [p], there is only one single round and this gives the Direct Send algorithm.

The Radix-k algorithm combines the advantages and the drawbacks of the Direct Send and the Binary Swap algorithms. Indeed:

• It is not limited to powers of two numbers of processors.

• It allows to partially overlapping computations with communication within a round, but not between the rounds. • The communication pattern is irregular and unpredictable within a given round.

It has been shown to perform better than the Binary Swap and 2-3 Swap algorithms.

Parallel Pipeline

The Parallel Pipeline algorithm [LRN96] requires (p -1) communication/computation rounds. The processing nodes are organized on a circular ring. At each round, each processing node composes the (n/p) pixels it has received from its "previous" neighbor in the ring and sends these (n/p) composed pixels to its "next" neighbor in the ring. The subimages are accumulated in a pipelined fashion along the ring.

The Parallel Pipeline algorithm presents two main advantages. First, it can support any number of processors. Second, it can benefit from the full-bisection bandwidth provided by a fat-tree IB cluster, since each node always sends data to the same node: contention free communication can be achieved using Fat-tree Unicast routing if the processes are placed with sufficient care. However, similarly to the Binary Swap algorithm, its synchronous nature does not permit to overlap computations with communications.

Shift-Based Image Compositing

We introduce in this Section our new parallel image compositing algorithm, called Shift-Based. We rely on the analysis done in the previous Section to build a new algorithm that combines all the advantages of existing image compositing methods. In particular, we base it on a shift-permutation pattern to ensure full utilization of the cross-bisectional bandwidth made available by fat-tree IB clusters.

Algorithm Requirements

Our algorithm is based on the five following requirements:

Algorithm Description

We make no assumption on the number p of nodes, which can be any value greater than 1 (requirement 2). In order to fulfill requirement 4, we base our algorithm on a shift permutation pattern. Similarly to the Direct Send method, each processing node P i is responsible of compositing the (n/p) pixels ranging from i(n/p) to (i + 1)(n/p) -1. Our algorithm simply consists in (p -1) compositing stages, as illustrated by Figure 1. During each stage s (0 < s < p), each processing node P i :

• Sends (n/p) pixels to processing node P (i+s)mod p and receives (n/p) pixels from processing node P (i-s+p)mod p . • Composes the pixels received during previous stage s -1 (for s > 1) while it does the send and receive operations.

It is easy to understand that this simple algorithm (requirement 5) keeps all processors busy and load balanced at all stages (requirement 1) and allows to overlap computations with communications (requirement 3).

Similarly to the Direct Send and Radix-k algorithms, partial subimages can arrive in a non-continuous order (see the example of processing node P 0 in Figure 1). In these cases, the non-continuous subimage is simply buffered until the missing subimage has been received (see for example at the end of Stage 3 for processing node P 0 in Figure 1).

Communication and Computation costs

Similarly to Peterka et al. [PGR * 09], we assume a simple cost model as proposed by Chan et al. [START_REF] Heimlich | Collective communication: theory, practice, and experience[END_REF] to analyze the communication and computation costs of the different existing parallel image compositing algorithms. We recall here the main assumptions:

• The target architecture is a distributed memory parallel architecture composed of p computational nodes indexed from P 0 to P p-1 . We consider the worst case scenario, where all pixels are active, and we do not consider optimizations to take advantage of the sparseness of the partial images. The communication term of the Radix-k is the one reported in [PGR * 09], and can be easily verified. On the contrary, we disagree with the value given by the authors for the latency term. Indeed, during each round, there are k i -1 messages exchanged, and they should be summed up for the r rounds. We claim that the latency term is then: α ∑ i=r i=1 (k i -1). It can be easily verified that: • For r = 1 (i.e. Direct Send), the latency term reduces to:

α(p -1). • For all k i equal to 2 (i.e. Binary Swap), the latency term reduces to α log 2 (p). • Depending on the k i values, the latency term is between Binary Swap (best case) and Direct Send (worst case).

The latency term of the Parallel Pipeline and Shift-Based algorithms is α(p -1), since there are (p -1) stages of communication. The communication term is bounded by β(p -1)(n/p), since each message contains (n/p) pixels.

Table 1 summarizes the communication and computation costs of the different algorithms. We do not specify the computation term, but rather whether the computations can be overlapped with the communications for the given algorithm. Except for the 2-3 Swap algorithm, all algorithms are equivalent in terms of communications. They only differ on the latency term and on their capability to overlap computations with communications.

At that point of the analysis, we are interested in the impact of the latency term of the parallel image compositing algorithms for modern high performance, low latency interconnects. For instance, for a 4x DDR (Double Data Rate) IB network, the latency and transmission time can be measured as: α = 1 × 10 -6 s and β = 0.53 × 10 -9 s (corresponding to 1800 MB/s): Figure 2 The latency term of the image compositing algorithms does not have a big impact on the total network time, except for very high processor counts (several thousands).

However, we recall that only the Parallel Pipeline and Shift-Based algorithms can benefit from the full crossbisectional bandwidth provided by IB fat-trees. The Direct Send and Radix-k algorithms all suffer from hot spots on large IB fat-trees, and would not benefit from the full 1800 MB/s provided by DDR IB. Figure 3 plots the evolution of the communication time we may observe when running the different algorithms on real fat-tree IB clusters (we anticipate a performance degradation of 30%, which is a low estimate as compared to real runs). The main advantage of our Shift-Based algorithm over the Binary Swap algorithm is its ability to overlap computations with communications. This has indeed a huge impact on the overall performance of the image compositing algorithm. On the other hand, it could be demonstrated that a routing scheme (different from the Fat-tree Unicast routing) could be computed to provide contention-free communication for the Binary Swap communication pattern.

In terms of scalability of the communication time on large IB fat-trees, the Parallel Pipeline and the Shift-Based algorithms are completely equivalent. The last, but not least, difference between the two algorithms is that the Parallel Pipeline does not permit to overlap computations with communications (similarly to the Binary Swap algorithm), while our Shift-Based algorithm does. We can conclude from our analysis that the Shift-Based algorithm is the most scalable on IB fat-trees with Fat-tree Unicast routing enabled.

Results

Our experimentations have been conducted on a 144-node fat-tree DDR IB cluster, running Linux CentOS 64-bit. Each node is composed of a single quad-core Intel Xeon X3440, and the nodes are interconnected with 24-port DDR switches on a 12-ary 2-tree topology. The OpenSM has been configured to run the Fat-tree Unicast routing engine. We used the methodology and tools of Moody [Moo09] to measure the influence of contention. During the execution of the Shift-Based algorithm on the cluster, each process measures the time it takes to send its data to every other pair process, and computes the associated bandwidth. A 2D matrix is built from these collected values, where each element at row a and column b represents the send bandwidth process Pa measured when it sent to process Pb. The diagonal elements of this matrix represent a task sending to itself, which is not measured. The minimum, average and maximum values of the matrix are computed. The 2D matrix is then converted into a gray-scale bitmap image by representing each entry as a pixel. The maximum bandwidth of the matrix is set to a pixel value of 255 (white), while all other values are scales to a pixel value of 0 (black) depending on their percentage of the maximum bandwidth. A histogram of the pixel values is also computed and shown next to each gray-scale image. The horizontal axis plots the pixel value (from 0/black on the left to 255/white on the right), and the vertical axis plots the number of pixels of a given value in the image.

Performance on a Crossbar Switch

Figure 4 shows a run of our Shift-Based algorithm on 12 nodes connected on a single 24-port crossbar switch, with an image size of 4096x2048x16 bytes. The average bandwidth is 1827 MB/s and the final frame rate is 14.5 FPS. It can be easily concluded that there is very few contention, and that our algorithm benefits from the full cross-bisectional bandwidth of the switch. As a point of comparison, our implementation of the Direct Send algorithm on the same 12 nodes on the same single 24-port crossbar switch leads to a frame rate of 6.4 FPS.

Performance on IB Fat-Trees

Figure 5 shows a run of our Shift-Based algorithm on 144 nodes with Fat-tree Unicast routing enabled, with an image size of 4096x2048x16 bytes. The 144 nodes have been chosen in the exact same order given by the OpenSM routing engine. The average bandwidth is 1431 MB/s (that is 78% of the cross-bisectional bandwidth) and the final frame rate is 11.5 FPS.

In order to validate the effects of the routing mechanism on the performance of our Shift-Based algorithm, we performed the same run, this time with the 144 nodes listed in alphabetical order (that is not in the order given by the OpenSM routing engine). Figure 6 shows the result of this Final frame rate is 11.5 FPS. Nodes are taken in the order given by the OpenSM routing engine. run: the average bandwidth is 924 MB/s (that is 50% of the cross-bisectional bandwidth) and the final frame rate is 6.5 FPS. The gray-scale image and the associated histogram highlight the high level of link contention that occurred during this run.

As a point of comparison, our implementation of the Direct Send algorithm on the same 144 nodes leads to a frame rate of 4.7 FPS, whether we use the order given by the OpenSM routing engine or the alphabetical order. We also wanted to show the influence of the size of the data transferred over the network. Figure 7 shows a run of our Shift-Based algorithm on 144 nodes with Fat-tree Unicast routing enabled, with an image size of 4096x4096x16 bytes. The 144 nodes have been chosen in the exact same order given by the OpenSM routing engine. The average bandwidth is 1568 MB/s (that is 85% of the cross-bisectional bandwidth) and the final frame rate is 5.8 FPS. As compared to Figure 5 (for which we used a 4096x2048x16 image size), the gray-scale image is much clearer, showing even less contention on the network, and a better utilization of the crossbisectional bandwidth.

As a final experiment, we have run our Shift-Based algorithm using all nodes counts from 12 to 144 of the 144 nodes with Fat-tree Unicast routing enabled, with several image sizes: 2048x2048x16 bytes, 4096x2048x16 bytes and 4096x4096x16 bytes. Figure 8 shows the evolution of the final frame rate as a function of the number of nodes used. The nodes are taken in the order given by the OpenSM routing engine. We can see that the obtained final frame rate is consistent with the size of the image being composed (the 2048x2048x16 bytes image is four times faster to be composed as compared to the 4096x4096x16 bytes image) and stable with increasing numbers of processing nodes (we do not see any random and unexplained variation). It is important to note that the obtained results are fully reproducible from one run to another. This can be explained by the predictable shift communication pattern and the absence of hot spots made possible by the Fat-tree Unicast routing engine. The only unanticipated behavior is the presence of peak performance for counts of nodes that are multiple of 12, where 12 is half of the arity of the switches used in our fat-tree topology (i.e. the number of leaf processing nodes connected to the first level switches). These peaks can be explained by the fact that when we do not use a fully populated level-1 switch, some links are missing to provide contention-free routing for our shift permutation, leading to hot spots on the used links.

Figure 9 shows a run of our Shift-Based algorithm on 64 nodes with Fat-tree Unicast routing enabled, with an image size of 4096x2048x16 bytes. The 64 nodes have been chosen in the exact same order given by the OpenSM routing engine. The average bandwidth is 1159 MB/s (that is 63% of the cross-bisectional bandwidth) and the final frame rate is 8 FPS. The gray-scale image and the associated histogram highlight many hot spots in our shift permutation pattern.

Summary

Conclusion

In this paper, we have explored parallel image compositing algorithms on a growing class of supercomputers: the fattree IB clusters. We have built upon an in-depth study of the fat-tree IB topology on one hand and of existing parallel image compositing algorithms on the other hand, to introduce Shift-Based, a new algorithm based on the well-known shift permutation communication pattern. Our Shift-Based algorithm combines all the advantages of previous image compositing algorithms, and allows benefiting from a big part of the cross-bisectional bandwidth made available by IB fattrees. We have proven that it is the most scalable algorithm for this kind of topology, and we have demonstrated its scalability on a 144-node fat-tree IB cluster with Fat-tree Unicast routing enabled. This cluster is the biggest one we have been able to have access to so far. However, we are confident that our conclusions can be extended to the case of much larger IB fat-trees, as long as a routing mechanism, such as the d-mod-k routing [START_REF] Zahavi | D-Mod-K routing providing non-blocking traffic for shift permutations on real life fat trees[END_REF], can be found to provide nonblocking traffic for shit permutations. Furthermore, we claim that the difference will be even bigger at higher processor counts, due to the increasing number of network contention for the other algorithms not relying on a shift permutation.

Future Work

The research work presented in this paper tends to prove that it is definitely worth the effort to take into account the network topology and its associated routing mechanisms in order to build efficient and scalable parallel algorithms. In the future, we plan to make similar investigations for parallel image compositing on other types of supercomputer architectures and to extend this work to other kinds of parallel algorithms.

Acknowledgments

Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see http://www.grid5000.fr).

  , Binary Swap [MPHK94] and Parallel Pipeline [LRN96]. The two latest advances in the field include 2-3 Swap by Yu et al. [YWM08], which is an extension of the Binary Swap to non-power of two numbers of processors, and the Radix-k algorithm by Peterka et al. [PGR * 09, MKPH11] which unic The Eurographics Association 2012.

  1. It should keep all processors busy and load balanced at all stages. 2. It should support any number of processors, such as the Direct Send, the 2-3 Swap and the Radix-k algorithms. 3. Computations and communications should be overlapped, such as the Direct Send and partially the Radix-k algorithms. 4. It should provide contention free communication on IB fat-trees, such as the Parallel Pipeline algorithm. 5. Finally, it must be easy to understand and simple to implement.

Figure 1 :

 1 Figure 1: The three stages of our Shift-Based parallel image compositing algorithm with p = 4 processing nodes.

  The latency term of the Direct Send is bounded by α(p -1), since each node has to send at most (p -1) messages. The communication term is bounded by β(p -1)(n/p), since each message contains (n/p) pixels. The latency term of the Binary Swap is α log 2 (p), since there are log 2 (p) stages of communication. The communication term is equal to: β ∑ i=log 2 (p) i=1 n/2 i = βn(p -1)/p. The upper bounds of the latency and communication terms of the 2-3 Swap are the ones reported in [YWM08], and correspond to the non-power of two case.

  plots the evolution of the communication time for increasing numbers of processing nodes for Communication and computation costs of the parallel image compositing algorithms.the Binary Swap algorithm (best latency case) and the Direct Send, Parallel Pipeline and Shift-Based algorithms (worst latency case) for a 4096x4096x16 image.

Figure 2 :

 2 Figure 2: Theoretical evolution of the communication time for increasing numbers of processing nodes on an ideal DDR IB crossbar for compositing of a 4096x4096x16 image.

Figure 3 :

 3 Figure 3: Theoretical evolution of the communication time for increasing numbers of processing nodes on a real DDR IB fat-tree for compositing of a 4096x4096x16 image.

Figure 4 :

 4 Figure 4: Shift-Based run on 12 nodes on a single 24-port crossbar switch (image size: 4096x2048x16 bytes). Final frame rate is 14.5 FPS.

Figure 5 :

 5 Figure 5: Shift-Based run on 144 nodes with the Fat-tree Unicast routing engine (image size: 4096x2048x16 bytes).Final frame rate is 11.5 FPS. Nodes are taken in the order given by the OpenSM routing engine.

Figure 6 :

 6 Figure 6: Shift-Based run on 144 nodes with the Fat-tree Unicast routing engine (image size: 4096x2048x16 bytes).Final frame rate is 6.5 FPS. Nodes are taken in the alphabetical order.

Figure 7 :

 7 Figure 7: Shift-Based run on 144 nodes with the Fat-tree Unicast routing engine (image size: 4096x4096x16 bytes). Final frame rate is 5.8 FPS. Nodes are taken in the order given by the OpenSM routing engine.

  As a point of comparison, our implementation of the Direct Send algorithm on the same 64 nodes leads to a frame rate of 5.2 FPS. c The Eurographics Association 2012.

Figure 8 :

 8 Figure8: Evolution of the image compositing frame rate with increasing numbers of processing nodes (from 12 to 144, with a step of 1). Nodes are taken in the order given by the OpenSM routing engine.

Figure 9 :

 9 Figure 9: Shift-Based run on 64 nodes with the Fat-tree Unicast routing engine (image size: 4096x2048x16 bytes). Final frame rate is 8 FPS. Nodes are taken in the order given by the OpenSM routing engine.

•

  The interconnection network is fully connected: each node can send directly to any other node. • At any given time, a single node can send only one message to one other node. Similarly, it can only receive one message from one other node. It is also assumed that a node can send and receive simultaneously. • The communication cost of sending a message between two nodes is modeled by α + nβ in the absence of network conflict. Here, α represents the message startup time (i.e. latency) and β represents the per data item transmission c The Eurographics Association 2012. time (i.e. reciprocal of bandwidth). The cost of the message is not a function of the distance between nodes. • Messages traveling in opposite directions on a link do not conflict. • The computation time to reduce one data item is denoted by γ.

c The Eurographics Association 2012.


Xuan-Quy Dao

Christophe Collewet

Drag Reduction of the Plane Poiseuille Flow by Partitioned Visual Servo Control

We present in this paper a way to minimize the drag of the 2D plane Poiseuille flow. To do that, as in our previous work [1], we use a vision-based approach to estimate the state of the flow. Since visual measurements are used, we propose to exploit visual servoing techniques to derive an efficient control law. Therefore, contrary to the literature concerning drag reduction, we also simultaneously minimize the kinetic energy density of the flow. That is of great importance since the controlled flow may become turbulent when this kinetic energy density is growing. To cope with this problem we propose to design a control law based on partitioned visual servo control. This approach has been first proposed in the robotics community by [2]. Simulations results validate our control scheme. We also compare our approaches with the most relevant ones.

I. INTRODUCTION

The goal of this paper is to show the benefit of using advanced visual servoing techniques [START_REF] Chaumette | Visual servo control, Part II: Advanced approaches[END_REF] for fluid flow control. By controlling a flow we mean either to change its state to another state or to maintain its current state whatever external disturbances. This can be of great environmental and economical interests for industrial applications. For instance, fuel consumption of aircrafts can be highly reduced by decreasing the drag while enhancing the lift (up to 20 %, see [START_REF] Flaig | Eco-efficient by design[END_REF]).

Flow control can be achieved in two different ways: passive or active control. Passive control provides a permanent action on a system. Most often it consists in optimizing shapes or in choosing suitable surfacing. Conversely, in active control an external energy is required to act on the system like, for example, techniques based on blowing and suction [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]. This approach can be seen as an optimal problem where one has to apply an optimal control law based on a certain cost (minimization of the drag, minimization of the actuators power, etc.) [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF]. However, very often, only openloop control or even most often forcing is used [START_REF] Protas | Drag force in the open-loop control of the cylinder wake in the laminar regime[END_REF]. Indeed, using a closed-loop control law requires sensors that have to be at the same time non-intrusive, accurate and adapted to the time and space scale of the phenomenon that we have to control. Unfortunately, such sensors are hardly available in the context of real control applications. The most commonly used measurement, obtained from Micro Electro Mechanical Systems (MEMS), is the shear stress at a limited set of measurement points [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]. Thereafter, an observer is required to estimate the state of the flow that will be used in the X-Q. Dao and C. Collewet are with INRIA Rennes -Bretagne Atlantique, Fluminance team, Rennes, France. E-mail: firstname.name@inria.fr. control law. Indeed, the shear stress based LQG regulator is the standard effective approach for flow control [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF], [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF].

In our previous paper [START_REF] Collewet | Fluid flow control: a visionbased approach[END_REF] this estimation is obtained from visual measurements. More precisely, the optical flow [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF] is computed and used to build an observer-free closed-loop control law as detailed later. Our approach has been revealed to be much more robust than the classical approaches [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF], [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF]; it is no more concerned with the well known problem of initialization issue when an observer is used. It has also been shown to be highly robust to measurements noise since a large number of measurements is available. To compare our approach with the classical one we also used the LQG control law to control the flow. In this paper, since visual measurements have been used, we aim to fully exploit the capabilities of visual servoing techniques by designing a more efficient control law.

Visual servoing or visual servo control aims to control the motions of a robot by using data provided by a vision sensor [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. It is now a well established technique in the robotics community. More precisely, to achieve a visual servoing task, a set of visual features s(t) is selected from the image of the scene being observed. A control law is then designed so that this set of visual features s(t) reaches a desired value s * corresponding to a desired state of the system. The dynamic of the error vector e(t) = s(t)s * is then given by

ė(t) = ∂e(t) ∂t + L e (t)u(t), (1)
where u(t) is the system control inputs, L e (t) is the so-called interaction matrix [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] or the image jacobian that links the time variation of the visual features to the variation of the control signal acting on the system and ∂e(t)/∂t expresses the variation of the error vector due to the free motion of the visual features. The control principle is of course to regulate the error vector e to zero. The control law is built from (1) using the knowledge or an approximation of the interaction matrix and an approximation of the free motion of the visual features. Visual servoing has shown impressive results in numerous complex contexts such as underwater, medical and aerial (helicopters, blimps) robotic as shown in [START_REF] Bonin-Font | Visual navigation for mobile robots: a survey[END_REF].

The final goal of this paper is to reduce the skin friction drag. As already mentioned at the top of this introduction, reducing the drag is an important challenge. Drag reduction in turbulent boundary layers has been investigated in the recent years. An adjoint-based suboptimal control approach has been proposed in [START_REF] Bewley | Optimal feedback control of turbulent channel flow[END_REF], a drag reduction of approximately 17% has been obtained. However this way to proceed leads to a very high computation cost, not compatible within a real time control scheme. Moreover, it is not possible to implement this approach in practice, see [START_REF] Kim | Control of turbulent boundary layers[END_REF] for more details. More viable approaches has been proposed in [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Joshi | Finite dimensional optimal control of poiseuille flow[END_REF]- [START_REF] Cortelezzi | Skin-friction drag reduction via robust reduced-order linear feedback control[END_REF]. In [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], a simple proportional controller was used to suppress the wall shear stress leading consequently to a drag reduction, while in [START_REF] Joshi | Finite dimensional optimal control of poiseuille flow[END_REF]- [START_REF] Cortelezzi | Skin-friction drag reduction via robust reduced-order linear feedback control[END_REF] a LQR approach has been used in which the cost function to minimize includes the drag. However, all these works does not account the kinetic energy density of the flow. This is yet an important issue since it is well known that if this density of energy is growing, transition of the controlled flow to turbulence is possible. That is why lots of papers focused on this problem (see e.g. [START_REF] Whidborne | On the minimization of maximum transient energy growth[END_REF], [START_REF] Mckernan | Linear quadratic control of plane poiseuille flow -the transient behaviour[END_REF]).

The main contribution of this paper is to show that, by fully exploiting the controlled degrees of freedom, a control law can be derived to simultaneously reduce the drag and the kinetic energy density of the flow. We will validate this approach on the 2D plane Poiseuille flow (see a description of this flow in section II). Note that an interesting approach has been proposed in [START_REF] Aamo | Control of mixing by boundary feedback in 2d channel flow[END_REF], in this work both a drag reduction and a reduction of the kinetic energy density have been obtained. However, tangential actuation has been used. This way to acts has been revealed to lead to lower controllability than wall normal actuation (see [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]). This is why, the more classical wall normal actuation is preferred in this paper as for example in [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF], [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF], [START_REF] Joshi | Finite dimensional optimal control of poiseuille flow[END_REF]- [START_REF] Cortelezzi | Skin-friction drag reduction via robust reduced-order linear feedback control[END_REF].

This paper is organized as follows: in Section II we first present some basics of the 2D plane Poiseuille controlled flow and the modeling model of some physical quantities that will be involved in the control law. The control scheme is detailed in Section III while in Section IV the simulation results are presented, we compare our approach with the most relevant ones.

II. MODELING ISSUES FOR CONTROL DESIGN

In this section we first present some basics of the plane Poiseuille flow, then we introduce the control principle for this flow and present a reduced linearized model around the steady state. The modeling of skin friction drag and the kinetic energy density are next presented.

A. Basics of the plane Poiseuille flow

Poiseuille flow is a flow in an infinite length channel due to a pressure gradient. Fig. 1 illustrates the steady state velocities profile of the plane Poiseuille flow in a streamwise period of the infinite length channel according to the conceptual model proposed in [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF]. The x axis corresponds to the streamwise direction and the y axis to the wall normal direction. The non dimensionalized Navier-Stokes equation of this flow is given by

         ∂V ∂t + (∇V)V = -∇P + 1 R e ∇ 2 V ∇ • V = 0 V(x, y = ±1, t) = 0 (2)
High pressure

Low pressure where P is the pressure; V is the flow velocity; V(x, y = ±1, t) = 0 represents the no slip boundary condition; R e is a dimensionless number called the Reynolds number; ∇ = ( ∂ ∂x , ∂ ∂y , ∂ ∂z ) represents the gradient operator and ∇ 2 the Laplacian.

M V(M, t) x y x = L c y = -1 x = 0 y = 1 h c
Since the Poiseuille flow is simple, the analytical solution (V bx , V by , P b ) of ( 2) in the steady state case, i.e. dV dt = ∂V ∂t + V • ∇V = 0, can be found:

(V bx , V by , P b ) = (1 -y 2 , 0, - 2 R e x). (3)
This solution is illustrated in Fig. 1.

B. Boundary control of the Poiseuille flow

The typical way to control the 2D plane flow Poiseuille is boundary control.

It consists in modifying the boundary conditions of the system (2) either only on the lower boundary y = -1 [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], or on both the upper y = 1 and lower y = -1 boundaries [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF]. Boundary control on the upper and the lower channels can be theoretically represented by functions χ u and χ l that allows mass conservation in the controlled system as pictured in Fig. 2.

Practically, fluidic devices (mean to inject and to suck fluid) like synthetic jet actuators can be used. Note those actuators are only one possible type of actuators for flow control since there are numerous other actuator control mechanisms such as plasma or other type which involves moving the domain boundary as described in [START_REF] Cattafesta | Actuators for active flow control[END_REF].

C. Reduced linearized model of Poiseuille flow

The modeling, required to derive the control law, consists first of all in linearizing the Navier-Stokes equation around the steady state solution [START_REF] Chaumette | Visual servo control, Part II: Advanced approaches[END_REF]. Then the continuous linearized model of this equation is first represented in the Fourier domain at a specific wavenumber α n = 2π/L c containing the instability [START_REF] Orszag | Accurate solution of the Orr-Sommerfeld stability equation[END_REF]: V p (x, y, t) = V n p (y, t)e jαn . Then the spectral representation V n p (y, t) is discretised in the y direction by using a number of M Gauss-Lobatto collocation points y k = cos ((k -1)π/(M -1)) , 1 ≤ k ≤ M leading to the following reduced linearized model (see [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF] for all the details):

ẋn (t) = A n x n (t) + B n u(t) x n (0) = x n 0 (4)
where x n (t) is the state vector of dimension M , A n is the state matrix which dimensions are (u u (t), u l (t)) is the system input (i.e. blowing and suction actions at the channel boundaries), B n is the input matrix which dimensions are M ×2. Recall that in our case (see [START_REF] Collewet | Fluid flow control: a visionbased approach[END_REF]) the state vector is directly obtained from the computation of the optical flow through a visualization system as shown in Fig. 3. More precisely, optical flow is the apparent velocity vector field representing the motion of photometric pattern (pixels luminance) in successive image sequences. Optical flow techniques can be used to estimate instantaneous velocities of a fluid flow from any image sequences as detailed in [START_REF] Heitz | Variational fluid flow measurements from image sequences: synopsis and perspectives[END_REF].

M × M , u(t) = pressure High Low pressure M V(M, t) x y x = 0 x = L c y = 1 y = -1 χ l (x, y 1 , t) χ u (x, y M , t)
In the case of Reynolds number Re = 10 000 and wavenumber α n = 1, the reduced linearized model of Poiseuille flow ( 4) is unstable as shown in [START_REF] Orszag | Accurate solution of the Orr-Sommerfeld stability equation[END_REF]. Therefore, as in [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF] we focus on the stabilization of the flow, i.e. maintaining the flow in the steady state, whatever the external perturbations. In this case the flow is initially in the steady state but in an unstable equilibrium, i.e. a small velocity pertubation value V p (x, y, t) destabilizes the non-controlled fluid flow. This instability can be easy seen through the poles of the matrix A n , a pole has a real part positive (λ = 0.00373967 ± j0.23752649). We will consider this case in the remainder of this paper. Note also that we will omit now the upperscript n involved in (4) for the sake of clarity.

D. Modeling of the skin friction drag

The skin friction drag d(t) due to the perturbation is given by (see e.g. [START_REF] Aamo | Control of mixing by boundary feedback in 2d channel flow[END_REF])

d(t) = - ∂V px (x, +1, t) ∂y + ∂V px (x, -1, t) ∂y (5)
where V px (x, +1, t) and V px (x, -1, t) are the streamwise components of the perturbation velocity V p (x, y, t) at the walls. Note that ∂V px (x, -1, z, t)/∂y and ∂V px (x, +1, z, t)/∂y are the components of the wall shear stress s(t) due to the perturbation:

s(t) =    ∂V px (x, +1, z, t) ∂y ∂V px (x, -1, z, t) ∂y    . (6)
According to modeling process recalled in section II-C, it can be shown that s(t) writes simply as a linear function of the state through a matrix W [8] leading consequently from (5) to

s(t) = Wx(t) (7)
d(t) = D ⊤ x(t) (8)
with D ⊤ = [ -1 +1 ]W where D ⊤ denotes the conjugate transpose of D.

E. Modeling of the kinetic energy density

The kinetic energy density of flow perturbation is given by

ε(t) = 1 V 0 V0 V P (x, y, t) 2 2 dV 0 (9)
where V 0 is the volume of a period of the domain under consideration. It is possible to obtain a weighting matrix Q exactly as in [START_REF] Mckernan | Linear quadratic control of plane poiseuille flow -the transient behaviour[END_REF] such that the kinetic energy density reduces to ε(t) = x(t) ⊤ Qx(t).

All these physical quantities having been defined, the control law can now be derived.

III. CONTROL OF THE 2D PLANE POISEUILLE FLOW

As mentioned in the introduction of this paper, our goal is to simultaneously the drag and the kinetic energy density since the flow is in unstable state (see section II-C). To achieve this goal we propose to use the so-called partionned visual servo control [START_REF] Corke | A new partitioned approach to imagebased visual servo control[END_REF]. It is a well known approach in the robotics community used to decouple the rotational motions from the translational ones. To do that, we first have to compute the interaction matrix L d related to the drag.

A. Computation of the interaction matrix related to the drag

This computation is straightforward. Indeed, according to definition of the interaction matrix [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF], we have to express the total time variation of the drag

ḋ(t) = L d u(t) + ∂d(t) ∂t . (11)
L d encodes the variation of the drag due to the actions, while ∂d(t)/∂t expresses the variation of the drag due to the flow itself.

This computation can be done by derivating [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF] with respect to the time and next using (4). We obtain

L d = D ⊤ B ∂d(t) ∂t = D ⊤ Ax(t). (12)
Let us introduce the components of the interaction matrix related to the components of the input signal u(t) (see ( 4)).

L d = (L du , L dl ) . (13)
This expression will be useful in the next section to introduce the partitioned visual servo control.

B. Design of the control law

Partitioned visual servo control relies on the partition of the interaction matrix. From ( 13), the time variation of the drag [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] becomes

ḋ(t) = L du u u (t) + L dl u l (t) + ∂d(t) ∂t . ( 14
)
Since a decreasing of the drag is desired, we choose a particular behavior for the drag. An exponential decrease is imposed

ḋ(t) = -λd(t) ( 15
)
where λ is a positive gain able to tune the decrease rate. Thereafter, from ( 14) and ( 15), we have

-λd(t) = L du u u (t) + L dl u l (t) + ∂d(t) ∂t . (16)
Any of the two components of u(t) can be used to reduce the drag, we arbitrary choose u l (t). We thus have

u l (t) = - 1 L dl λd(t) + L du u u (t) + ∂d(t) ∂t . ( 17
)
Since the lower boundary control law u l (t) is known if the upper boundary control law u u (t) is known, the next step is to express u u (t).

In order to minimize the kinetic energy density [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] and the energy consumption of actuators, a LQR control scheme over an infinite time horizon is used by considering the following cost function

J = ∞ 0 x ⊤ (t)Qx(t) + ru 2 u (t) dt. ( 18
)
This can be done by expressing the time variation of the state vector with respect to the control signal u u . To do that, we first rewrite ( 17) from ( 12) to exhibit the state vector

u l (t) = - 1 L dl λD ⊤ + D ⊤ A x(t) + L du u u (t) (19)
that we have to plug in (4) leading to

ẋ(t) = A 1 x(t) + B 1 u u (t) (20)
where we have introduced the following matrices

A 1 = A -B l L dl (λD ⊤ + D ⊤ A) and B 1 = B u -L du L dl B l with B = (B u , B l ). Note that L dl is always non null.
Thereafter, it becomes straightforward to compute the LQR gain k involved in the optimal control u u = -k ⊤ x(t) by considering [START_REF] Aamo | Control of mixing by boundary feedback in 2d channel flow[END_REF] and solving the algebraic Ricatti equation.

Note that here, contrary to the works involved in flow control, we have fully exploited the controlled degrees of freedom.

The next section is dedicated to the study of the behavior of the closed loop system when measurement noise is considered.

C. Behavior of the closed loop system in presence of measurement noise

Of course, when measurement noise occurs, the main problem is the stability analysis of [START_REF] Aamo | Control of mixing by boundary feedback in 2d channel flow[END_REF]. In that case, the upper control signal becomes

ûu = -k ⊤ x(t). (21)
However, it has been shown in [START_REF] Collewet | Fluid flow control: a visionbased approach[END_REF] that in our vision-based approach the state estimation writes as

x(t) = x(t) + 1 N x e x (t) (22)
where N x is the number of pixels of the camera in the streamwise direction of the flow and e x (t) is related to the measurements noise (see [START_REF] Collewet | Fluid flow control: a visionbased approach[END_REF] for more details). It is clear from [START_REF] Cattafesta | Actuators for active flow control[END_REF], that the larger the value of N x the closer x(t) is from x(t). Consequently, from ( 21), [START_REF] Aamo | Control of mixing by boundary feedback in 2d channel flow[END_REF] becomes

ẋ(t) = (A 1 -B 1 k ⊤ )x(t) -B 1 kδ(t) (23)
where the estimation error δ(t) = x(t)x(t) has been introduced. When N x is large enough (it is always the case in practice), the estimation error tends toward 0, and, consequently the closed loop system (23) writes as

ẋ(t) = (A 1 -B 1 k ⊤ )x(t) (24)
which is stable. other issue concerns the drag. Indeed, we have to verify that it is still a decreasing function. To do that, we exam ( 14) by considering ûu , ûl and ∂ d(t)/∂t. From (19), we have

ûl (t) = - 1 L dl λD ⊤ + D ⊤ A x(t) + L du ûu (t) (25)
that writes simply in function of x(t) by considering ( 21) and ( 22)

ûl (t) = - 1 L dl λD ⊤ + D ⊤ A -L du k ⊤ x(t) + 1 N x e x (t) .
(26) The last step is the computation of ∂ d(t)/∂t, it can be easy done from [START_REF] Bonin-Font | Visual navigation for mobile robots: a survey[END_REF] and [START_REF] Cattafesta | Actuators for active flow control[END_REF]. All computations done, [START_REF] Kim | Control of turbulent boundary layers[END_REF] becomes

ḋ(t) = -λd(t) - λD ⊤ + D ⊤ A N x e x (t). (27)
Here again, since N x is a large value, we can consider than the drag follows the desired behavior ḋ(t) = -λd(t). These results will be validated in the next section.

IV. SIMULATION RESULTS

The simulation results are based on the code provided in [START_REF] Mckernan | Control of plane Poiseuille flow: a theoretical and computational investigation[END_REF]. As in [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Mckernan | A linear state-space representation of plane poiseuille flow for control design-a tutorial[END_REF], [START_REF] Bewley | Optimal and robust control and estimation of linear paths to transition[END_REF], we consider R e = 10000 and a channel length L c = 4π (see Fig. 1). The scalar r involved in [START_REF] Whidborne | On the minimization of maximum transient energy growth[END_REF] has been set to 2.5 while λ, involved in the desired behavior of the drag, has been set to 0.5. The dimension of the state vector M has been set to 40.

First, we present on Fig. 4 the behavior of the velocity field during the servoing. The flow is first not controlled and becomes thus turbulent due to the presence of a temporal perturbation and because its unstable state. After the application of the input signal, the flow tends towards a laminar behavior but, of course, still in an unstable state.

We compare our approach with the most relevant ones, that is [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], [START_REF] Cortelezzi | Skin-friction drag reduction via robust reduced-order linear feedback control[END_REF], [START_REF] Mckernan | Linear quadratic control of plane poiseuille flow -the transient behaviour[END_REF]; [START_REF] Joshi | Finite dimensional optimal control of poiseuille flow[END_REF], [START_REF] Lee | Application of robust reduced-order controller to turbulent flows for drag reduction[END_REF] are similar in mind than [START_REF] Cortelezzi | Skin-friction drag reduction via robust reduced-order linear feedback control[END_REF]. We respectively denote these approaches P , LQR 2 and LQR 1 , we denote ours by P V S. Let us briefly describe these works. In [START_REF] Joshi | A system theory approach to the feedback stabilization of infinitesimal and finite amplitude disturbances in plane Poiseuille flow[END_REF], a proportional output feedback control is used, the output is the wall shear stress. Therefore, vanishing the shear stress also vanishes the drag. In this approach the kinetic energy density is not considered. In [START_REF] Mckernan | Linear quadratic control of plane poiseuille flow -the transient behaviour[END_REF] the kinetic energy density is minimized by using a LQR approach. The matrix Q involved in [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] is used in the cost function to minimize. Note that is an indirect way to minimize the drag since the drag vanishes when the kinetic energy density does. Finally, in [START_REF] Cortelezzi | Skin-friction drag reduction via robust reduced-order linear feedback control[END_REF], the wall shear stress is minimized also through a LQR approach.

Note that we cannot exactly compared our approach with these works since none of them simultaneously minimizes the drag and the kinetic energy density. The results are depicted on Fig. 5. More precisely, Fig. 5a depicts the behavior of the drag versus time while Fig. 5b depicts the behavior of the kinetic energy density versus time. As expected, our approach provides better results than the other approaches. A nice decreasing of the drag and the kinetic energy density are observed. The LQR 1 approach leads also to a nice minimization of the kinetic energy density but the (indirect) minimization of the drag is not satisfactory. Both methods based on the minimization of the shear stress (P and LQR 2 ) lead to a worse behavior than P V S and LQR 1 . Note also that, except for our approach, the decreasing of the drag towards zero is very low.

However, the input signal provides both by the P and LQR 2 methods are very low in contrast to the P V S and LQR 1 approaches (see Fig. 6). Therefore, by tuning the matrix R involved in the LQR 2 approach or the gain for the P approach better results could be expected. Unfortunately, it is not the case as proved by the results depicted on Fig. 7. In particular, oscillations occur in the kinetic energy density (Fig. 7b).

Fig. 8 describes the variation of the kinetic energy density for various choices of λ, as can be seen this kinetic energy density does not depend, in practice, on this value. That means that the way the drag is reduced does not influence the way the kinetic energy density is reduced.

The next simulation concerns the behavior of the system when noisy measurement are considered. The results are depicted on Fig. 9, with or without noise. Respectively the drag, the kinetic energy density and the norm of the input signal are represented on Fig. 9a, Fig. 9b and Fig. 9c. A Gaussian noise with a standard deviation of 0.3 has been added to the velocity fields. The number of pixels N x has been set to 2048. As can be seen, despite this very low value, very good results are obtained. In particular, the convergence of the control law is ensured, as well as the decrease of the drag. Note that, contrary to our vision-based estimation scheme, all the works previously described are very sensitive to measurement noise. The control law may be diverged (see [START_REF] Collewet | Fluid flow control: a visionbased approach[END_REF] for more details).

V. CONCLUSION

We have shown in this paper that partitioned visual servo control is an efficient way to simultaneously minimize the drag and the kinetic energy density contrary to the existing methods. Moreover, we have shown that in practice, the way the drag is reduced does not influence the way the kinetic energy density is reduced. In addition, we have proved that our approach is very robust against measurement noise. This nice propriety is due to the large amount of flow velocities available through the computation of the optical flow. In contrast, the existing methods are sensitive to noise since the state estimation is based on shear stress measurements on only few points.

Fig. 1 .

 1 Fig. 1. Steady state velocities profile of the 2D plane Poiseuille flow: hc is the channel half height, Lc is the flow streamwise period length and M is a point in the flow.

Fig. 2 .

 2 Fig. 2. Boundaries control of the plane Poiseuille flow.

Fig. 3 .

 3 Fig.3. Visualization of the flow using a laser sheet which role is to enlighten the particles seeded in the fluid.

Fig. 4 .

 4 Fig. 4. Velocity field vs time. (a) due to an initial perturbation the flow becomes turbulent; (b) first iteration of the control law; (c) and (d) intermediate state of the controlled flow; (e) the flow becomes laminar; (f) the control law converged, the flow is laminar in an unstable state.

Fig. 5 .

 5 Fig. 5. Comparison of the different approaches (x axis in second). (a) drag vs time; (b) kinetic energy density vs time.

Fig. 6 .Fig. 7 .

 67 Fig.6. Comparison of the different approaches (x axis in second). Norm of the input signal u(t) vs time.

Fig. 8 .

 8 Fig. 8. Behavior of the kinetic energy density for different values of λ (x axis in second).

Fig. 9 .

 9 Fig. 9. Behavior when measurement noise is considered (x axis in second). (a) drag vs time; (b) kinetic energy density vs time; (c) Norm of the input signal u(t) vs time.


Denis Cousineau

Damien Doligez

Leslie Lamport

Stephan Merz

Daniel Ricketts

Hernán Vanzetto

TLA + Proofs

TLA + is a specification language based on standard set theory and temporal logic that has constructs for hierarchical proofs. We describe how to write TLA + proofs and check them with TLAPS, the TLA + Proof System. We use Peterson's mutual exclusion algorithm as a simple example and show how TLAPS and the Toolbox (an IDE for TLA + ) help users to manage large, complex proofs.

Introduction

TLA + [START_REF] Lamport | Specifying Systems: The TLA + Language and Tools for Hardware and Software Engineers[END_REF] is a specification language originally designed for specifying concurrent and distributed systems and their properties. It is based on Zermelo-Fraenkel set theory for modeling data structures and on the linear-time temporal logic TLA for specifying system executions and their properties. More recently, constructs for writing proofs have been added to TLA + , following a proposal for presenting rigorous hand proofs in a hierarchical style [START_REF] Lamport | How to write a 21st century proof[END_REF].

In this paper, we present the main ideas that guided the design of the proof language and its implementation in TLAPS, the TLA + proof system [START_REF] Chaudhuri | Verifying safety properties with the TLA + proof system[END_REF][START_REF]The TLAPS Project[END_REF]. The proof language and TLAPS have been designed to be independent of any particular theorem prover. All interaction takes place at the level of TLA + . Users need know only what sort of reasoning TLAPS's backend provers tend to be good at-for example, that SMT solvers excel at arithmetic. This knowledge is gained mostly by experience.

TLAPS has a Proof Manager (PM) that transforms a proof into individual proof obligations that it sends to backend provers. Currently, the main backend provers are Isabelle/TLA + , an encoding of TLA + as an object logic in Isabelle [START_REF] Wenzel | The Isabelle framework[END_REF], Zenon [START_REF] Bonichon | Zenon : An extensible automated theorem prover producing checkable proofs[END_REF], a tableau prover for classical first-order logic with equality, and a backend for SMT solvers. Isabelle serves as the most trusted backend prover, and when possible, we expect backend provers to produce a detailed proof that is checked by Isabelle. This is currently implemented for the Zenon backend.

TLAPS has been integrated into the TLA + Toolbox, an IDE (Integrated Development Environment) based on Eclipse for writing TLA + specifications and running the TLA + tools on them, including the TLC model checker. The Toolbox provides commands to hide and unhide parts of a proof, allowing a user to focus on a given proof step and its context. It is also invaluable to be able to run the model checker on the same formulas that one reasons about.

We explain how to write and check TLA + proofs, using a tiny well-known example: a proof that Peterson's algorithm [START_REF] Peterson | Myths about the mutual exclusion problem[END_REF] implements mutual exclusion. We start by writing the algorithm in PlusCal [START_REF] Lamport | The PlusCal algorithm language[END_REF], an algorithm language that is based on the expression language of TLA + . The PlusCal code is translated to a TLA + specification, which is what we reason about. Section 3 introduces the salient features of the proof language and of TLAPS with the proof of mutual exclusion. Liveness of Peterson's algorithm (processes eventually enter their critical section) can also be asserted and proved with TLA + . However, liveness reasoning makes full use of temporal logic, and TLAPS cannot yet check temporal logic proofs.

Section 4 indicates the features that make TLA + , TLAPS, and the Toolbox scale to realistic examples. A concluding section summarizes what we have done and our plans for future work.

Modeling Peterson's Algorithm In TLA +

Peterson's algorithm is a classic, very simple two-process mutual exclusion algorithm. We specify the algorithm in TLA + and prove that it satisfies mutual exclusion: no two processes are in their critical sections at the same time. 6A representation of Peterson's algorithm in the PlusCal algorithm language is shown on the left-hand side of Figure 1. The two processes are named 0 and 1; the PlusCal code is embedded in a TLA + module that defines an operator Not so that Not(0) = 1 and Not(1) = 0.

The variables statement declares the variables and their initial values. For example, the initial value of flag is an array such that flag[0] = flag [START_REF] Balser | Formal system development with KIV[END_REF] = false. (Mathematically, an array is a function; the TLA + notation [x ∈ S → e] for writing functions is similar to a lambda expression.) To specify a multiprocess algorithm, it is necessary to specify what its atomic actions are. In PlusCal, an atomic action consists of the execution from one label to the next. With this brief explanation, the reader should be able to figure out what the code means.

A translator, normally called from the Toolbox, generates a TLA + specification from the PlusCal code. We illustrate the structure of the TLA + translation in the right-hand part of Figure 1. The heart of the TLA + specification consists of the predicates Init describing the initial state and Next, which represents the next-state relation.

The PlusCal translator adds a variable pc to record the control state of each process. The meaning of formula Init in the figure is straightforward. The formula Next is the disjunction of the two formulas proc(0) and proc(1), which are in turn defined as disjunctions of formulas corresponding to the atomic steps of the process. In these formulas, unprimed variables refer to the old state and primed variables to the new state. The temporal formula Spec is the complete specification. It characterizes behaviors (ω-sequences of states) that start in a state satisfying Init and where every pair of successive states either satisfies Next or else leaves the values of the tuple vars unchanged. 7Before trying to prove that the algorithm is correct, we use TLC, the TLA + model checker, to check it for errors. The Toolbox runs TLC on a model of a TLA + specification. A model usually assigns particular values to specification constants, such as the number of processes. It can also restrict the set of states explored, which is useful if the specification allows an infinite number of reachable states. TLC easily verifies that the two processes can never both be at label cs by checking that the following formula is an invariant (true in all reachable states):

Init ∆ = ∧ flag = [i ∈ {0, 1} → false] ∧ turn = 0 ∧ pc = [self ∈ {0, 1} → "a0"] a3a(self ) ∆ = ∧ pc[self ] = "a3a" ∧ if flag[Not(self )] then pc = [pc except ![self ] = "a3b"] else pc = [pc except ![self ] = "cs"] ∧ unchanged flag, turn \ * remaining actions omitted proc(self ) ∆ = a0(self ) ∨ . . . ∨ a4(self ) Next ∆ = ∃ self ∈ {0, 1} : proc(self ) Spec ∆ = Init ∧ 2[Next]vars
MutualExclusion ∆ = (pc[0] = "cs") ∨ (pc[1] = "cs")
Peterson's algorithm is so simple that TLC can check all possible executions. For more interesting algorithms that have parameters (such as the number of processes) and perhaps an infinite set of reachable states, TLC cannot exhaustively verify all executions, and correctness can only be proved deductively. Still, TLC is invaluable for catching errors, and it is much easier to run TLC than to write a formal proof.

theorem Spec ⇒ 2MutualExclusion 1 1. Init ⇒ Inv 1 2. Inv ∧ [Next]vars ⇒ Inv 1 3. Inv ⇒ MutualExclusion 1 4. qed

Proving Mutual Exclusion For Peterson's Algorithm

The assertion that Peterson's algorithm implements mutual exclusion is formalized in TLA + as the theorem in Figure 2. The standard method of proving this invariance property is to find an inductive invariant Inv such that the steps 1 1-1 3 of Figure 2 are provable.

TLA + proofs are hierarchically structured and are generally written topdown. Each proof in the hierarchy ends with a qed step that asserts the proof's goal. We usually write the qed step's proof before the proofs of the intermediate steps. The qed step follows easily from steps 1 1-1 3 by standard proof rules of temporal logic. However, TLAPS does not yet handle temporal reasoning, so we omit that step's proof. When temporal reasoning is added to TLAPS, we expect it easily to check such a trivial proof.

Figure 3 defines the inductive invariant Inv as the conjunction of two formulas. (A definition must precede its use, so the definition of Inv appears in the module before the proof.) The first, TypeOK , asserts simply that the values of all variables are elements of the expected sets. (The expression [S → T ] is the set of all functions whose domain is S and whose range is a subset of T .) In an untyped logic like that of TLA + , almost any inductive invariant must assert type correctness. The second conjunct, I , is the interesting one that explains why Peterson's algorithm implements mutual exclusion. We again use TLC to check that Inv is indeed an invariant. In our simple example, TLC can even check that Inv is inductive, by checking that it is an (ordinary) invariant of the specification Inv ∧ 2[Next] vars , obtained from Spec by replacing the initial condition by Inv .

TypeOK ∆ = ∧ pc ∈ [ {0, 1} → {"
a0", "a1", "a2", "a3a", "a3b", "cs", "a4"} ] We now prove steps 1 1-1 3. We can prove them in any order; let us start with 1 1. This step follows easily from the definitions, and the following leaf proof is accepted by TLAPS: by def Init, Inv , TypeOK , I TLAPS will not expand definitions unless directed to so. In complex proofs, automatically expanding definitions often leads to formulas that are too big for provers to handle. Forgetting to expand some definition is a common mistake. If a proof does not succeed, the Toolbox displays the exact proof obligation that it passed to the prover. It is then usually easy to see which definitions need to be invoked.

∧ turn ∈ {0, 1} ∧ flag ∈ [ {0, 1} → boolean ] I ∆ = ∀i ∈ {0, 1} : ∧ pc[i] ∈ {"a2", "a3a", "a3b", "cs", "a4"} ⇒ flag[i] ∧ pc[i] ∈ {"cs", "a4"} ⇒ ∧ pc[Not(i)] / ∈ {"cs", "a4"} ∧ pc[Not(i)] ∈ {"a3a", "a3b"} ⇒ turn = i Inv ∆ = TypeOK ∧ I
Step 1 3 is proved the same way, by simply expanding the definitions of MutualExclusion, Inv , I , and Not. We next try the same technique on 1 2. A little thought shows that we have to tell TLAPS to expand all the definitions in the module up to and including the definition of Next, except for the definition of Init. Unfortunately, when we direct TLAPS to prove the step, it fails to do so, reporting a 65-line proof obligation.

TLAPS uses Zenon and Isabelle as its default backend provers. However, TLAPS also includes an SMT solver backend [START_REF] Merz | Automatic verification of TLA + proof obligations with SMT solvers[END_REF] that is capable of handling larger "shallow" proof obligations-in particular, ones that do not contain significant quantifier reasoning. We instruct TLAPS to use the SMT backend when proving the current step by writing by SMT def . . . The backend translates the proof obligation to the input language of SMT solvers. In this way, step 1 2 is proved in a few seconds. For sufficiently complicated algorithms, an SMT solver will not be able to prove inductive invariance as a single obligation. Instead, the proof will have to be hierarchically decomposed. We illustrate how this is done by writing a proof of 1 2 that can be checked using only the Zenon and Isabelle backend provers.

The outline of a hierarchical proof of step 1 2 appears in Figure 4. The proof introduces more elements of the TLA + proof language that we now explain.

A suffices step allows a user to introduce an auxiliary assertion, from which the current goal can be proved. For example, step 2 1 reduces the proof of the implication asserted in step 1 2 to assuming predicates Inv and Next, and proving Inv . In particular, this step establishes that the invariant is preserved by stuttering steps that leave the tuple vars unchanged. Steps 2 2 and 2 3 establish the two conjuncts in the definition of Inv . Whereas 2 2 can be proved directly by Isabelle, 2 3 needs some more interaction.

Following the definition of predicate I as a universally quantified formula, we introduce in step 3 1 a new variable j , assume that j ∈ {0, 1}, and prove I !(j ) , which denotes the body of the universally quantified formula, with j substituted for the bound variable, and with primed copies of all state variables. Similarly, step 3 2 introduces variable i to denote the process that makes a transition, following the definition of Next (which is assumed in step 2 1). Even after this elimination of two quantifiers, Isabelle and Zenon cannot prove the goal in a single step. The usual way of decomposing the proof is to reason separately about each atomic action a0(i ), . . . , a4(i ). However, Peterson's algorithm is simple enough that we can just split the proof into the two cases i = j and i = j with steps 3 3 and 3 4. Isabelle and Zenon can now prove all the steps.

Writing Real Proofs

Peterson's algorithm is a tiny example. Some larger case studies have been carried out using the system [START_REF] Lamport | Byzantizing Paxos by refinement[END_REF][START_REF] Lu | Towards verification of the Pastry protocol using TLA +[END_REF][START_REF] Parno | Memoir: Practical state continuity for protected modules[END_REF]. Several features of TLAPS and its Toolbox interface help in coping with the complexity of large proofs.

Hierarchical Proofs And The Proof Manager

Hierarchical structure is the key to managing complexity. TLA + 's hierarchical and declarative proof language enables a user to keep decomposing a complex proof into smaller steps until the steps become provable by one of the backend provers. In logical terms, proof steps correspond to natural-deduction sequents that must be proved in the current context. The Proof Manager tracks the context, which is modified by non-leaf proof steps. For leaf proof steps, it sends the corresponding sequent to the backend provers, and records the result of the step's proof that they report.

Proof obligations are independent of one another, so users can develop proofs in any order and work on different proof steps independently. The Toolbox makes it easy to instruct TLAPS to check the proof of everything in a file, of a single theorem, or of any step in the proof hierarchy. Its editor helps reading and writing large proofs, providing commands that show or hide subproofs. Although some other interactive proof systems offer hierarchical proofs, we do not know of other systems that provide the Toolbox's abilities to use that structure to aid in reading and writing proofs and to prove steps in any order.

Hierarchical proofs are much better than conventional lemmas for handling complexity. In a TLA + proof, each step with a non-leaf proof is effectively a lemma. One typical 1100-line invariance proof [START_REF] Lamport | Byzantizing Paxos by refinement[END_REF] contains 100 such steps. A conventional linear proof with 100 lemmas would be impossible to read.

Unlike most interactive proof assistants [START_REF] Wiedijk | The Seventeen Provers of the World[END_REF], TLAPS is independent of any specific backend prover. There is no way for a user to indicate how available facts should be used by backends. TLA + proofs are therefore less sensitive to changes in any prover's implementation.

Fingerprinting: Tracking The Status Of Proof Obligations

During proof development, a user repeatedly modifies the proof structure or changes details of the specification. By default, TLAPS does not re-prove an obligation that it has already proved-even if the proof has been reorganized. It can also show the user the impact of a change by indicating which parts of the existing proof must be re-proved.

The Proof Manager computes a fingerprint of every obligation, which it stores, along with the obligation's status, in a separate file. The fingerprint is a compact canonical representation of the obligation and the relevant part of its context. The Toolbox displays the proof status of each step, indicating by color whether the step has been proved or some obligation in its proof has failed or been omitted. The only other proof assistant that we know to offer a mechanism comparable to our fingerprinting facility is the KIV system [START_REF] Balser | Formal system development with KIV[END_REF].

Conclusion

The proof of Peterson's algorithm illustrates the main constructs of the hierarchical and declarative TLA + proof language. The algorithm is so simple that we had to eschew the use of the SMT solver backend so we could write a nontrivial proof. Section 4 explains why TLAPS, used with the TLA + Toolbox, can handle more complex algorithms and specifications.

A key feature of TLAPS is its use of multiple backend provers. Different proof techniques, such as resolution, tableau methods, rewriting, and SMT solving offer complementary strengths. Future versions of TLAPS will probably support additional backend provers. Because multiple backends raise concerns about soundness, TLAPS provides the option of having Isabelle certify proof traces produced by backend provers; and this has been implemented for Zenon. Still, it is much more likely that a proof is meaningless because of an error in the specification than that it is wrong because of an error in a backend. Soundness also depends on parts of the proof manager.

We cannot overstate the importance of having TLAPS integrated with the other TLA + tools-especially the TLC model checker. Finding errors by running TLC on finite instances of a specification is much faster and easier than discovering them when writing a proof. Also, verifying an algorithm or system may require standard mathematical results. For example, the correctness of a distributed algorithm might depend on known facts about graphs. Engineers want to assume such results, not prove them. However, it is easy to make a mistake when formalizing mathematics. TLC can check the exact TLA + formulas assumed in a proof (on finite instances), greatly reducing the chance of introducing an unsound assumption.

We are actively developing TLAPS. Our main short-term objective is to add support for temporal reasoning. We have designed a smooth extension of the existing proof language to sequents containing temporal formulas. We also plan to improve support for standard TLA + data structures such as sequences.

-

  -algorithm Peterson { variables flag = [i ∈ {0, 1} → false], turn = 0; process (proc ∈ {0, 1}) { a0: while (true) { a1: flag[self ] := true; a2: turn := Not(self ); a3a: if (flag[Not(self )]) {goto a3b} else {goto cs} ; a3b: if (turn = Not(self )) {goto a3a} else {goto cs} ; cs: skip; \ * critical section a4: flag[self ] := false; } \ * end while } \ * end process } \ * end algorithm variables flag, turn, pc vars ∆ = flag, turn, pc

Fig. 1 .

 1 Fig. 1. Peterson's algorithm in PlusCal (left) and in TLA + (excerpt, right).

Fig. 2 .

 2 Fig.2. The high-level proof.

Fig. 3 . 2 .Fig. 4 .

 324 Fig. 3. The inductive invariant.

The TLA + module containing the specification and proof as well as an extended version of this paper are accessible at the TLAPS Web page[START_REF]The TLAPS Project[END_REF].

"Stuttering steps" are allowed in order to make simple[START_REF] Lamport | What good is temporal logic?[END_REF].

This work was partially funded by Inria-Microsoft Research Joint Centre, France.