Tifenn Rault

Abdelmadjid Bouabdallah

Yacine Challal

Frédéric Marin

Context-aware energy-efficient wireless sensor architecture for body activity recognition

In this paper, we present EEWAA, a new Energy-Efficient Wireless sensor Architecture for Activity recognition Applications. Contrary to existing approaches, our solution aims at optimizing the lifetime of every battery-operated devices that may contribute to the activity recognition, including sensors but also smart phones or tablets. Depending on the environment opportunities, our approach balance the load among multiple base stations (mobile phone, tablet, laptop) and when beneficial, enables the energy-limited base stations and/or sensors to offload a part of their tasks to back-end servers. In order to do this, we propose a new decision metric to select the best execution configuration, considering the latency and accuracy requirements, as well as the energy of every device involved in the computation.

I. INTRODUCTION

The recognition of human activities has received a significant attention in the medical field, because it is a task of interest for many applications, including assisted living of elderlies, physical assessment, and supervision of patient with chronic diseases such as diabetes, obesity, cognitive disorder or arrhythmia [START_REF] Lara | A Survey on Human Activity Recognition based on Wearable Sensors[END_REF]. In a typical scenario, wearable sensors like accelerometers or ECG are attached to the patient to measure inertial and physiological data. These data are further analyzed to provide a feedback to caregivers who can assess the efficiency of a new treatment, adjust medication, better study a disease or supervise a patient behavior [START_REF] Lorincz | Mercury: a wearable sensor network platform for high-fidelity motion analysis[END_REF]. In practice, such systems face a major challenge that limit their widespread adoption : the significant energy consumption of batteryoperated devices for continuous sensing and communications. Indeed, radio transmissions, data processing and sensing tasks actively stress the nodes reserves, whereas they are expected to operate during the course of the day without maintenance.

To tackle this problem, many solutions have been proposed to save energy in activity recognition applications. Among them we can cite sensor set selection [START_REF] Noshadi | Energy optimization in wireless medical systems using physiological behavior[END_REF]- [START_REF] Gao | Activity recognition using dynamic multiple sensor fusion in body sensor networks[END_REF], deactivation of power-hungry sensors [START_REF] Lorincz | Mercury: a wearable sensor network platform for high-fidelity motion analysis[END_REF], [START_REF] Sun | PEAR: Power efficiency through activity recognition (for ECG-based sensing)[END_REF], [START_REF] Park | E-Gesture: a collaborative architecture for energy-efficient gesture recognition with hand-worn sensor and mobile devices[END_REF], adaptive sampling rate [START_REF] Qi | AdaSense: Adapting Sampling Rates for Activity Recognition in Body Sensor Networks[END_REF], [START_REF] Yan | Energy-Efficient Continuous Activity Recognition on Mobile Phones: An Activity-Adaptive Approach[END_REF], communications reduction [START_REF] Wang | A hierarchical approach to realtime activity recognition in body sensor networks[END_REF], duty-cycling [START_REF] Au | Episodic sampling: Towards energy-efficient patient monitoring with wearable sensors[END_REF] and resource sharing [START_REF] Keally | Remora: Sensing resource sharing among smartphone-based body sensor networks[END_REF], [START_REF] Rachuri | METIS: Exploring mobile phone sensing offloading for efficiently supporting social sensing applications[END_REF]. However, these works are usually interested in the energy consumption of the sensors but never consider that the base station can also be energy-constrained. Traditionally, powerful laptop acted as base station, but for mobility consideration most of the works now consider mobile phone as base station. Although phones are more powerful devices compared to wearable sensors, the truth is that they also face energy consumption problems [START_REF] Balasubramanian | Energy consumption in mobile phones: a measurement study and implications for network applications[END_REF].

We believe that activity recognition architectures can take advantage of the proliferation of smart devices (phones, tablets, laptops) and their daily usage to balance the load among base stations. For example, when at home and at the office, body sensors can communicate with a laptop or an access point. When outdoor, a mobile phone or a tablet can alternatively take the relay depending on their availability and battery state. Furthermore, additional energy can be saved at a base station through computation offloading, that is, a base station can execute (a part of) its tasks on a backend server. This approach has proven to be energy-efficient under some conditions [START_REF] Rachuri | METIS: Exploring mobile phone sensing offloading for efficiently supporting social sensing applications[END_REF], [START_REF] Kumar | Cloud Computing for Mobile Users: Can Offloading Computation Save Energy[END_REF], [START_REF] Rachuri | Sociable-Sense: exploring the trade-offs of adaptive sampling and computation offloading for social sensing[END_REF], and the usage of devices with higher computing, storage and energy capabilities allows to incorporate more complex algorithms and models. As computation offloading is a trade-off between energy and performances, we propose a new metric to determine when offloading is beneficial, based on the energy consumption of every devices involved in activity recognition, the latency and the classification accuracy. Besides, our solution is selfadaptive to the environment, that is the balancing and the offloading strategies are highly dependent on the available connections and devices, as well as their internal state. Finally, our architecture is complementary to the energy-efficient mechanisms already proposed in the literature, as we are still able to implement them to further enhance the system lifetime.

In what follows, in Section II we study traditional architectures for activity recognition applications. In section III, we present the general idea of our EEWAA architecture. In section IV, we propose a new metric for offloading decisions. Section V concludes the paper and highlights challenges.

II. ACTIVITY RECOGNITION

The main stages for activity recognition include sensing, processing and analysis, as illustrated in Figure 1. More specifically, sensing corresponds to the data acquisition at sensor nodes. Data processing generally consists in filters and features extraction to obtain relevant characteristics of the signal. The analysis refers to the activity recognition which is done through classification. Depending on where these tasks are executed, the energy/performance profile of the body sensor network varies. For example, the sensors have the possibility to directly store the raw data, leaving all the processing for later. The drawbacks of this method is that it prevents real-time feedback to the patient and it requires a lot of memory, even if to save space, the sensor can store extracted features. Another solution consists in sending the raw data to the base station, but it results in a huge amount of energy consumption due to transmissions. In the meantime, using a more powerful base station allows to have access to more complex algorithms and models. Furthermore, to save communications, the sensor can send only the extracted features, the compressed signal or the results, which decreases communication cost but increases the computation load. From the energy point of view it may be optimal, but sensor scarce resources limit the implementation to lightweight algorithms. Actually, some classification algorithms are very expensive in terms of memory and computation requirements which makes them not suitable for sensor platforms [START_REF] Lara | A Survey on Human Activity Recognition based on Wearable Sensors[END_REF].

Wang et al. [START_REF] Wang | Energy and Accuracy Trade-Offs in Accelerometry-Based Activity Recognition[END_REF] compare the energy/performance trade off between an on-node scheme where the classification is performed on the accelerometer, and an off-node scheme where the accelerometer sends its raw data to the base station. They show that with the on-node scheme 40% energy can be saved with 13% reduction in classification accuracy. So, their results suggest that it is more energy-efficient to perform all the tasks on the sensors. However, we can make several observations to nuance this statement. First, as mentioned by the authors, if the user stays in active state or if he/she changes activity often, on-node scheme may consume more energy due to additional processing and data communications. Second, the study only consider accelerometers, but other kind of sensors may require more complex processing like ECG. In this case, on-node signal processing may not be possible, or at the cost of accuracy loss. Third, computation offloading is more energy-efficient than local computation in some scenarios. For example, in SociableSense [START_REF] Rachuri | Sociable-Sense: exploring the trade-offs of adaptive sampling and computation offloading for social sensing[END_REF] it is more energy-efficient to offload the speaker identification task when using a microphone. In METIS [START_REF] Rachuri | Sociable-Sense: exploring the trade-offs of adaptive sampling and computation offloading for social sensing[END_REF], sensing offloading can save up to 35% energy over pure phone sensing. Fourth, energy is not the unique metric of interest and some applications may not accept a loss in accuracy. In these cases, leveraging energy consumption for more accurate classification may be necessary to maintain the desired performance.

As we can see, a static architecture can not satisfy the requirements of every body sensor networks. Indeed, activity recognition applications are constrained by the types of sensors and processing to perform, the patient behavior and the resource availability. That's why we propose a more general architecture that optimizes the energy/performance trade-off depending on the applications requirements and the environment opportunities. The idea is to switch from an operating mode to another when there is a change in the context. For example, if a patient leaves its office to go at a meeting, the sensors can stop sending their data to the PC and redirect their transmissions towards the patient's phone. If the phone has enough energy to last until the end of the day, it can process the received signals. Otherwise, the phone will only carry out the classification task. So, EEWAA supports many execution scenarios thanks to context-aware adaptation, which enables us to leverage energy consumption for accuracy and vice versa.

III. EEWAA ARCHITECTURE

Our architecture is composed of:

• A set of body sensors (e.g. accelerometer, EMG, gyroscope, ECG). • A set of ambient sensors (e.g. a microphone in a room, an RFID tag on objects, a GPS in a phone). • A set of candidate base stations, mobile or not, with different energy, storage and computing capabilities (e.g. laptop, mobile phone, tablet, access point). • A back-end server. At a given time t, only a subset of base stations may be available (possibly none), for two reasons : a base station may not be in the patient environment or it is in the patient environment but there is no wireless connections between the sensors and the base station. Consequently, the connection to the back-end server does not always exists. Similarly, there may not be ambient sensors in the environment. We denote by D t = {d 1 , d 2 , .., d q } the set of devices available at time t. Each device (sensor/base station/server) can execute a set of tasks, like sensing, signal processing, classification or result display. For a given environment, there exist several execution configurations, depending on where each task is executed and resulting in different energy/performance profile. Through the load balancing and offloading strategies, the objective of EEWAA is to select the configuration that achieves a satisfactory tradeoff. We give in section IV a new metric for configuration selection.

In Figure 2 we illustrate an execution configuration. In this scenario, the phone is elected as base station (BS). The phone is equipped with a microphone, but it offloads its sensing task to a microphone embedded in the environment. The BS receives extracted features from the accelerometer and the microphone, and it receives raw data from the ECG and the gyroscope. The BS computes the features for the gyroscope but it offloads the task for the ECG. So, the back-end server perform the feature extraction of ECG data, send the result back to the phone and store them in the data base. Finally, the phone is able to perform the classification using the features extracted from every signal and display the result. We can imagine that after a while, if the phone battery has decreased below a threshold, the tablet can take the relay as base station.

We first construct a matrix of constraints, where (i, j) represents the costs of computing a task j ∈ [1, n] on a device i ∈ [1, q], as illustrated in Table I. If the device i cannot execute the task j, then the cost is set to infinity. The cost may be a tuple of different costs, e.g. (e ij , l ij ) which are respectively the energy and latency costs of executing the task j on device i. We then periodically select a configuration using the metric defined in section IV, so that every task is assigned to one and only one device.

IV. OFFLOADING STRATEGIES

In this section, we present a new metric for computation offloading decision in EEWAA. Our objective is to decide where to execute each task, based on energy-efficiency, latency and accuracy requirements. It is often assumed a client/server architecture, where a mobile device can offload a part of its tasks to a back-end server. In this scenario, only the energy consumption of the mobile is of interest. This is the case of SociableSense [START_REF] Rachuri | Sociable-Sense: exploring the trade-offs of adaptive sampling and computation offloading for social sensing[END_REF], an architecture dedicated to activity recognition for sociable interactions, where a phone with integrated sensors can offload its tasks to the cloud. SociableSense selects the configuration c i that maximizes a weighted additive utility function [START_REF] Lara | A Survey on Human Activity Recognition based on Wearable Sensors[END_REF], which takes into consideration the energy consumption, the latency and the amount of data sent over the network.

u ci = w e u ei + w l u li + w d u di (1)
with

u xi = x min -x i x i x ∈ (e, l, d) (2)
where e i is the energy consumption of the architecture when choosing the configuration c i , l i is the resulting latency and d i is the amount of data transmitted over the network. Depending on the phone battery state, more prevalence can be given to the energy metric by adjusting the weights w i dynamically. Few works have considered that mobile devices can offload their data between them [START_REF] Shi | Serendipity: enabling remote computing among intermittently connected mobile devices[END_REF], [START_REF] Mtibaa | Towards resource sharing in mobile device clouds: power balancing across mobile devices[END_REF]. In this case, the energy of every device is of interest. Mtibaa et al. [START_REF] Mtibaa | Towards resource sharing in mobile device clouds: power balancing across mobile devices[END_REF] maximize the lifetime of the whole network, which is define as the time until the first mobile runs out of energy. They take into consideration the cost for computing locally E comp and the costs for transmitting and receiving E comm . A node u offloads a task to a neighboring node v if the response time of migrating the task to v is within the task deadline and the estimated remaining energy is maximized as in Eq. ( 3) [START_REF] Mtibaa | Towards resource sharing in mobile device clouds: power balancing across mobile devices[END_REF].

max v∈G {(E u t -E comm(u,v) )+(E v t -E comp(v) -E comm(v,u) )} (3
) where E u t is the remaining energy of node u. After selecting v, if v == u the node executes its task locally, otherwise it offloads the task. The results show that this approach manages to fairly distribute the load among nodes.

Shi et al. [START_REF] Shi | Serendipity: enabling remote computing among intermittently connected mobile devices[END_REF] maximizes the lifetime of the first depleted node under the constraint that tasks complete before their deadline. They minimize the following utility function to take into consideration the energy consumption of all nodes involved in the remote computing and the residual energy available to these nodes:

u(T ) = i∈N T e T i R i (4)
where N T is the set of nodes involved in the remote computing of task T , e T i is the energy consumption of node i for task T , and R i is the residual energy of node i [START_REF] Shi | Serendipity: enabling remote computing among intermittently connected mobile devices[END_REF]. We propose a new metric for decision offloading in EE-WAA, that takes into account the application accuracy, the latency, and the energy spent by every battery-operated devices involved in the activity recognition. We aim at maximizing a weight additive utility function as in [START_REF] Rachuri | Sociable-Sense: exploring the trade-offs of adaptive sampling and computation offloading for social sensing[END_REF], while considering :

• δ i , the error accuracy of the classifier used in configuration c i . • l i , the latency of the application in configuration c i .

• e i , the energy consumption of the architecture in c i .

u ci = w e u ei + w l u li + w δ u δi (5)
We can define e i as the sum of the energy spent in every battery-operated devices involved in the computation over their residual energy as in [START_REF] Shi | Serendipity: enabling remote computing among intermittently connected mobile devices[END_REF]. However, it is not fair because phones are not only dedicated to the health supervision application, so the fraction of their energy dedicated to the activity recognition application must not be as high as the one of the sensors. In other words, a part of phones energy must be preserved to be dedicated to other tasks, like calling or text messaging. That's why we propose to introduce a factor α j that represents the desired implication of the device j in the application. We redefine e i as in Eq [START_REF] Sun | PEAR: Power efficiency through activity recognition (for ECG-based sensing)[END_REF]:

e i = j∈Ni α j E ij R j (6)
where N i is the set of devices involved in the computation in configuration c i , E ij is the energy spent by the device j in configuration c i , R j is its residual energy and α j is a weight that represents how much a device can dedicate its energy to the activity recognition application. An α j equal to 1 means that the device is exclusively dedicated to the healthcare application, while an higher value means that the device will not spend all its energy in the application. The energy consumption E ij of a device j in c i includes the energy spent for sensing, computing and communications. Notably, the energy spent in communications will depend on the size of data to be transferred [START_REF] Kumar | Cloud Computing for Mobile Users: Can Offloading Computation Save Energy[END_REF], the bandwidth availability [START_REF] Wolski | Using bandwidth data to make computation offloading decisions[END_REF] and the communication technology (Wi-Fi, 3G) [START_REF] Cuervo | MAUI: making smartphones last longer with code offload[END_REF].

V. CONCLUSION

This paper describes the work in progress towards the development of EEWAA, an Energy-Efficient Architecture for Activity recognition Applications. EEWAA adapts its operating mode to context changes in order to take advantage of the opportunities provided by the environment. Many challenges need to be addressed. First, it will be necessary to manage data fusion and synchronization when signal processing are distributed over multiple devices, like in the scenario of Figure 2. Moreover, the presence of a remote server can be exploited to implement more complex classification algorithms based on prediction. It would be possible to progressively construct of a data base of the patient habits and infer additional knowledge.

Fig. 1 .

 1 Fig. 1. The main stages of activity recognition applications.

Fig. 2 .

 2 Fig. 2. An example of an execution configuration.

ACKNOWLEDGMENTS

This work was carried out and funded in the framework of the Labex MS2T. It was supported by the French Government, through the program "Investments for the future" managed by the National Agency for Research (ANR-11-IDEX-0004-02).


Keywords: persistent homology, simplicial complex, software library, computational topology, generic programming Bibliothèque de Calcul Gudhi : Complexes Simpliciaux et Homologie Persistante homologie persistante, complexes simpliciaux, bibliothèque de calcul, topologie algorithmique, programmation générique

Nous présentons le design de la bibliothèque de calcul Gudhi. Gudhi est une bibliothèque

Introduction

The principle of algebraic topology is to attach algebraic invariants to topological spaces in order to classify them up to homeomorphism. One can consequently study the property of a discrete algebraic structure (a sequence of homology groups in our case) instead of studying a continuous domain directly, which would be hard to handle algorithmically. Persistent homology [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF][START_REF] Zomorodian | Computing persistent homology[END_REF] may be considered as a "dynamic version" of this principle: given a sequence of topological spaces connected by continuous maps, we study the corresponding sequence of homology groups connected by group homomorphisms, induced by the topological space maps. The whole sequence of groups together with their homomorphisms form an algebraic structure (specifically a module) that we study. Very efficient methods have been developped for computing persistent homology [START_REF] Bauer | Clear and compress: Computing persistent homology in chunks[END_REF][START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF] and its dual, persistent cohomology [START_REF] Boissonnat | The compressed annotation matrix: An efficient data structure for computing persistent cohomology[END_REF][START_REF] Vin De Silva | Persistent cohomology and circular coordinates[END_REF][START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF]. The generality and stability with regard to noise [START_REF] Cohen-Steiner | Stability of persistence diagrams[END_REF] of persistence have made it a widely used tool in practice.

An application of interest for computational topology is topological data analysis, where one is interested in learning topological invariants of a shape, sampled by a point cloud. A popular approach is to construct, at different scales, an approximation of the shape using complexes built on top of the points, and then compute the persistent homology of these complexes. This approach has been successfully used in various areas of science and engineering, as for example in sensor networks [START_REF] Silva | Coverage in sensor network via persistent homology[END_REF], image analysis [START_REF] Carlsson | On the local behavior of spaces of natural images[END_REF], and data analysis [START_REF] Chazal | Towards persistence-based reconstruction in euclidean spaces[END_REF], where one typically needs to deal with big data sets in high dimensions and with general metrics. The simplicial complex and persistent homology packages in Gudhi provide all software components for this approach.

The challenge is twofold. On the one hand we need to design a generic library in computational topology, in order to adapt to the various configurations of the problem: nature of the complexes (simplicial, cubical, etc) and their representation, nature of the maps between them (inclusions, edge contractions, etc), ordering of the maps (linear, zigzag, etc) and types of algorithm for persistence (homology, cohomology). On the other hand, we need to implement a high-performance library to handle complex practical examples.

We recall in Section 2 the definition of homology and persistent homology constructed from simplicial complexes. In Section 3, we describe the design of the Gudhi library. In Section 4, we discuss the implementation choices and the user interface. Specifically, simplicial complexes are implemented with a simplex tree data structure [START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF]. The simplex tree is an efficient and flexible data structure for representing general (filtered) simplicial complexes. The persistent homology of a filtered simplicial complex is computed by means of the persistent cohomology algorithm [START_REF] Vin De Silva | Persistent cohomology and circular coordinates[END_REF][START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF], implemented with a compressed annotation matrix [START_REF] Boissonnat | The compressed annotation matrix: An efficient data structure for computing persistent cohomology[END_REF]. The persistent homology package provides the computation of persistence with different coefficient fields, including the implementation of the multi-field persistence algorithm of [START_REF] Boissonnat | Computing Persistent Homology with Various Coefficient Fields in a Single Pass[END_REF], i.e. the simultaneous computation of persistence with various coefficient fields. Finally, in Section 5 we discuss the future components of the library and their integration in the design.

Theoretical Foundation of Persistent Homology

The theory of homology consists in attaching to a topological space a sequence of (homology) groups, capturing global topological features like connected components, holes, cavities, etc. Persistent homology studies the evolution -birth, life and death -of these features when the topological space is changing. Consequently, the theory is essentially composed of three elements: topological spaces, their homology groups and an evolution scheme.

Indexing

Complex

Homology 0 0 0 0

1 0 3 4 2 5 6 7 [c] [c] [c] [c], [c ′ ]
Figure 1: Indexing of eight simplicial complexes and corresponding sequence of homology groups in dimension 1.

Simplicial Complexes:

In computer science, topological spaces are represented by their discrete counterpart: (cell) complexes. On the following, we focus on simplicial complexes, but our approach applies to all kinds of cell complexes.

Let V = {1, • • • , |V |} be a set of vertices. A simplex σ is a subset of vertices σ ⊆ V . A simplicial complex K on V is a collection of simplices {σ}, σ ⊆ V , such that τ ⊆ σ ∈ K ⇒ τ ∈ K. The dimension n = |σ| -1 of σ is its number of elements minus 1.
A simplicial map f : K → K ′ between simplicial complexes K and K ′ , with respective vertex sets V and V ′ , is a map that sends every vertex v ∈ V to a vertex f (v) ∈ V ′ , and every simplex

[v 0 , • • • , v n ] ∈ K to a simplex [f (v 0 ), • • • , f (v n )] ∈ K ′ . Note that they may be redundancy in the set {f (v 0 ), • • • , f (v n )}, in
which case the simplex image has lower dimension that its pre-image. In the following, we focus on inclusions, which are a particular case of simplicial maps, and discuss the case of general simplicial maps in Section 5.

Homology:

For a ring R, the group of n-chains, denoted C n (K, R), of K is the group of formal sums of n-simplices with R coefficients. The boundary operator is a linear operator

∂ n : C n (K, R) → C n-1 (K, R) such that ∂ n σ = ∂ n [v 0 , • • • , v n ] = n i=0 (-1) i [v 0 , • • • , v i , • • • , v n ]
, where v i means v i is omitted from the list. The chain groups form a sequence:

• • • C n (K, R) ∂n ---→ C n-1 (K, R) ∂n-1 ---→ • • • ∂2 ---→ C 1 (K, R) ∂1 ---→ C 0 (K, R)
of finitely many groups C n (K, R) and homomorphisms ∂ n , indexed by the dimension n ≥ 0. The boundary operators satisfy the property ∂ n • ∂ n+1 = 0 for every n > 0 and we define the homology groups:

H n (K, R) = ker ∂ n /im ∂ n+1
We refer to [START_REF] Munkres | Elements of algebraic topology[END_REF] for an introduction to homology theory and to [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF] for an introduction to persistent homology.

Indexing Scheme:

"Changing" a simplicial complex consists in applying a simplicial map. An indexing scheme is a directed graph together with a traversal order, such that two consecutive nodes in the graph are connected by an arrow (either forward or backward). The nodes represent simplicial complexes and the directed edges simplicial maps. From the computational point of view, there are two types of indexing schemes of interest in persistent homology 1 : [START_REF] Zomorodian | Computing persistent homology[END_REF], and [START_REF] Carlsson | Zigzag persistence[END_REF]. These indexing schemes have a natural left-to-right traversal order, and we describe them with ranges and iterators. We focus in the following on the linear case, and discuss the zigzag case in Section 5.

Inria

linear ones • -→ • -→ • • • -→ • -→ • in persistent homology
zigzag ones • -→ • ←-• • • -→ • ←-• in zigzag persistent homology
In the following, we consider the case where the indexing scheme is induced by a filtration. A filtration of a simplicial complex is a function f :

K → R satisfying f (τ ) ≤ f (σ) whenever τ ⊆ σ.
Ordering the simplices by increasing filtration values (breaking ties so as a simplex appears after its subsimplices of same filtration value) provides an indexing scheme.

We refer to Figure 1 for an illustration of the three components of the theory and their connections. The figure pictures the linear indexing of eight simplicial complexes connected by inclusions, and the corresponding sequence of homology groups in dimension 1. Every inclusion induces a group homomorphism at the homology level. Persistent homology studies this sequence of homology groups connected by homomorphisms. Specifically, computing persistent homology consists in computing a primary decomposition of this sequence of homology groups (forming a module); the decomposition is usually represented by means of a persistence diagram [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF].

Remark:

The reader may have found a category theory taste to this presentation of persistent homology. In particular, the vertical arrows in Figure 1 represent functors of categories. We refer to [START_REF] Bubenik | Categorification of persistent homology[END_REF] for more details on the categorification of persistent homology.

Design of the Library

A concept is a set of requirements (valid expression, associated types, etc) for a type. If a type satisfies these requirements, it is a model of the concept. The general idea under our design is to associate a concept per component presented in Section 2: the three components of the theory (indexing, complex and homology) are illustrated in Figure 1. Given two components related by a vertical arrow in Figure 1, and two models A and B of their respective associated concepts, we connect B with A through a template argument B<A>.

IndexingTag Concept:

In order to describe the indexing scheme, we use a tag IndexingTag that is either void compute_persistent_homology( FilteredComplex cpx ) { for( Simplex_handle sh : cpx.filtration_simplex_range() ) { int dim = cpx.dimension(sh); update_cohomology_groups( dim, sh, cpx ); //inside update_cohomology_groups for( Simplex_handle b_sh : cpx.boundary_simplex_range(sh) ) {...} //out } } } linear_indexing_tag or zigzag_indexing_tag, corresponding to the two indexing schemes of interest mentioned above. The tag is passed as template argument to a model of the concept FilteredComplex (described below and representing filtered cell complexes).

FilteredComplex Concept:

We define the concept FilteredComplex that describes the requirement for a type to implement a filtered cell complex. We use the vocabulary of simplicial complexes, but the concept is valid for any type of cell complex. The main requirements are the definition of: 1. type Indexing_tag, which is a model of the concept IndexingTag, describing the nature of the indexing scheme, 2. type Simplex_handle to manipulate simplices, 3. method int dimension(Simplex_handle) returning the dimension of a simplex, 4. type and method Boundary_simplex_range boundary_simplex_range( Simplex_handle) that returns a range giving access to the codimension 1 subsimplices of the input simplex, as-well-as the coefficients (-1) i in the definition of the operator ∂. The iterators have value type Simplex_handle, 5. type and method Filtration_simplex_range filtration_simplex_range() that returns a range giving access to all the simplices of the complex read in the order assigned by the indexing scheme, 6. type and method Filtration_value filtration(Simplex_handle) that returns the value of the filtration on the simplex represented by the handle.

Figure 3 illustrates the use of a model of the concept FilteredComplex. It sketches the algorithm used for computing persistent homology via the approach of [START_REF] Vin De Silva | Persistent cohomology and circular coordinates[END_REF][START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF].

PersistentHomology Concept:

The concept PersistentHomology describes the requirement for a type to compute the persistent homology of a filtered complex. The requirement are the definition of:

1. a type Filtered_complex, which is a model of FilteredComplex and provides the type of complex on which persistence is computed, Inria 2. a type Coefficient_field, which is a model of CoefficientField and provides the coefficient field on which homology is computed.

The requirements of the concept CoefficientField are essentially the definition of field operations (addition, multiplication, inversion, etc).

We refer to Figure 2 for a presentation of the concepts and their connections.

Implementation

In this section we describe how these concepts are implemented. The code will be available soon at project.inria.fr/gudhi/software/.

Simplicial Complex:

We use a Simplex Tree [START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF] to represent simplicial complexes. The class Simplex_tree is a model of FilteredComplex and hence furnishes all requirements of the concept. Moreover, it furnishes algorithms to construct efficiently simplicial complexes, and in particular flag complexes [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF]. Details on the implementation of the algorithms may be found in [START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF].

Persistent Homology:

We use the Compressed Annotation Matrix [START_REF] Boissonnat | The compressed annotation matrix: An efficient data structure for computing persistent cohomology[END_REF] to implement the persistent cohomology algorithm [START_REF] Vin De Silva | Persistent cohomology and circular coordinates[END_REF][START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF] for persistence. This leads to the class Persistent_cohomology, which is a model of PersistentHomology. The class Persistent_cohomology allows the computation of the persistence diagram of a filtered complex, using the method compute_persistent_homology (see Figure 3). The coefficient fields available as models of CoefficientField are Field_Zp for Z p (for any prime p) and Multi_field for the multi-field persistence algorithm -computing persistence simultaneously in various coefficient fields -described in [START_REF] Boissonnat | Computing Persistent Homology with Various Coefficient Fields in a Single Pass[END_REF].

Example of Use of the Library:

Figure 4 illustrates the user interface for constructing a flag complex [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF] from a graph and computing its persistent homology with various coefficient fields.

Graph g; ... //compute the graph Simplex_tree< linear_indexing_tag > st; //linear ordering st.insert(g); //insert the graph as 1-skeleton of the complex st.expand(5); //construct the 5-skeleton of the associated flag complex Persistent_cohomology< Simplex_tree<linear_indexing_tag>, Multi_field > pcoh;

//persistence with "multi field coefficients" defined on a simplex tree pcoh.compute_persistent_homology(st,2,1223); //compute persistent homology of st in all fields Zp for p prime between 2 and 1223

Figure 4: User interface for the construction of a filtered flag complex with a simplex tree and the computation of its persistent homology.

Data

Experiments:

Figure 5 presents timings T st for the construction of flag complexes with a simplex tree using the algorithm of [START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF], T ph Z2 and T ph Z1223 for the computation of persistent homology with coefficient is Z 2 and Z 1223 respectively, using the implementation of [START_REF] Boissonnat | The compressed annotation matrix: An efficient data structure for computing persistent cohomology[END_REF], and T ph

Z 2 1223
for the simultaneous computation of persistent homology in the 200 coefficient fields Z p with p prime, for 2 ≤ p ≤ 1223, using the multi-field persistent homology algorithm described in [START_REF] Boissonnat | Computing Persistent Homology with Various Coefficient Fields in a Single Pass[END_REF]. Experiments have been realized on a Linux machine with 3.00 GHz processor and 32 GB RAM, for Rips complexes [START_REF] Edelsbrunner | Computational Topology -an Introduction[END_REF] built on a variety of data points. Datasets are listed in Figure 5 with the size of points sets |P|, the ambient dimension D and intrinsic dimension d of the sample points ("?" if unknown), the parameter r for the Rips complex and the size of the complex |K|. More details about the implementation, the experimental protocol, the data sets as-well-as additional experiments can be found in [START_REF] Boissonnat | The compressed annotation matrix: An efficient data structure for computing persistent cohomology[END_REF][START_REF] Boissonnat | Computing Persistent Homology with Various Coefficient Fields in a Single Pass[END_REF][START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF].

The average timings per simplex of the various algorithms are ranging between 4.08•10 -8 and 7.28 • 10 -8 seconds per simplex for the construction of the simplex tree, between 1.27 • 10 -6 and 2.00 • 10 -6 seconds per simplex for the computation of persistent homology with coefficient field Z 2 or Z 1223 , and between 1.68 • 10 -6 and 3.52 • 10 -6 seconds per simplex for the computation of multi-field persistent homology in all fields Z p for p prime, 2 ≤ p ≤ 1223. Note that most of the time for the computation of persistent homology is spent computing boundaries in the simplex tree.

Future Components

The library may be extended in various directions that fit naturally in the design. The first direction is to allow zigzag indexing schemes, by the creation of a tag zigzag_indexing_tag. In this case, the method filtration_simplex_range must indicate the direction of the arrows.

New implementations and models for FilteredComplex may be added. For example, the construction of witness complexes [START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF] will be added to the class Simplex_tree. Additionnaly, new types of complexes (like cubical complexes) and new data structures to represent them may be added to the library: in order to compute their persistent homology, they only need to satisfy the requirements of the concept FilteredComplex.

So far, only inclusions have been considered for simplicial maps between simplicial complexes. As explained in [START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF], any simplicial map may be implemented with a sequence of inclusions and edge contractions. We will consequently add edge contractions as updates in the class Simplex_tree and implement the induced updates in the class Persistent_cohomology (algorithms exist for edge contractions in a simplex tree [START_REF] Boissonnat | The simplex tree: An efficient data structure for general simplicial complexes[END_REF] and for the corresponding updates at the cohomology level [START_REF] Tamal | Computing topological persistence for simplicial maps[END_REF]). This way, we will be able to compute persistent homology of simplicial maps. In this case, the range provided by filtration_simplex_range must indicate the nature of the map between complexes.

Inria

Future works include also the implementation of a class Field_Q, model of concept CoefficientField, for homology with Q coefficients. Finally the interface between complexes and persistent homology allows us to implement more persistent homology algorithms.

FilteredComplexFigure 2 :

 2 Figure 2: Overview of the design of the library.

Figure 3 :

 3 Figure 3: Sample code for the computation of persistence, illustrating the use of a model of concept FilteredComplex.

Figure 5 :

 5 Figure 5: Timings in seconds for the various algorithms.

RR n°8548

i.e. from which an interval decomposition of the persistence module exists: Gabriel's theorem[START_REF] Derksen | Quiver representations[END_REF] in quiver theory classifies these graphs.RR n°8548

InriaRESEARCH CENTRE SOPHIA ANTIPOLIS -MÉDITERRANÉE 2004 route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgment

This research has been partially supported by the European Research Council under Advanced Grant 339025 GUDHI (Algorithmic Foundations of Geometric Understanding in Higher Dimensions).


Alexis Fouilhe
email: alexis.fouilhe@imag.fr

Sylvain Boulmé
email: sylvain.boulme@imag.fr

A certifying frontend for (sub)polyhedral abstract domains ⋆

Keywords: abstract interpretation, abstract domain of polyhedra, program verification in Coq

Convex polyhedra provide a relational abstraction of numerical properties for static analysis of programs by abstract interpretation. We describe a lightweight certification of polyhedral abstract domains using the Coq proof assistant. Our approach consists in delegating most computations to an untrusted backend and in checking its outputs with a certified frontend. The backend is free to implement relaxations of domain operators (i.e. a subpolyhedral abstract domain) in order to trade some precision for more efficiency, but must produce hints about the soundness of its results. Previously published experimental results show that the certification overhead with a full-precision backend is small and that the resulting certified abstract domain has comparable performance to non-certifying state-of-the-art implementations.

Introduction

Astrée [START_REF] Blanchet | A static analyzer for large safety-critical software[END_REF] is a major success of semantics-based static analysis of programs: it is capable proving the absence of runtime undefined behaviours in large scale real world C programs from avionics. Abstract interpretation [START_REF] Cousot | Abstract interpretation: a unified lattice model for static analysis of programs by construction or approximation of fixpoints[END_REF], on which Astrée is based, formalizes the state analysis of programs and guarantees that the analyzer soundly over-approximates the behaviours of the program under analysis. However, Astrée is itself a complex piece of software. Despite the care put in its development, it may contain bugs. One possible solution consists in proving that the analyzer implementation is sound and having this proof mechanically checked by a proof assistant. Trusting the result of the analyzer is thereby reduced to trusting the proof checker and answering the question: Is what has been proved what we want to prove?

This question is especially relevant in the context of automated C program manipulation, where the semantics of the C programming language are part of the specification. The CompCert C compiler [START_REF] Leroy | Formal verification of a realistic compiler[END_REF] is a successful project built with the Coq proof assistant. The Verasco project aims at building an abstract interpreter in a similar manner, reusing the CompCert infrastructure: mainly the formalized semantics of the C programming language and the frontend. Once the semantics of the program under analysis are defined, the correctness proof is composed of two components: the abstract domain must be shown to soundly over-approximate manipulations of sets of reachable states, and the link between the program semantics and the abstract domain must be proved correct.

Our work addresses the problem of proving correct in Coq an implementation of the abstract domain of polyhedra [START_REF] Cousot | Automatic discovery of linear restraints among variables of a program[END_REF], which capture linear relationships between program variables. The abstract domain we built is similar both in features and performance to the core of the polyhedra library in the PPL [START_REF] Bagnara | The Parma Polyhedra Library: Toward a complete set of numerical abstractions for the analysis and verification of hardware and software systems[END_REF] and Apron [START_REF] Jeannet | Apron: A library of numerical abstract domains for static analysis[END_REF]. We adopted the same two tier architecture as Besson et al. [START_REF] Besson | Result certification for relational program analysis[END_REF]: an untrusted Ocaml backend performs most of the computations and outputs proof hints for the results it produces, which are used by a frontend developed in Coq to build trustworthy results.

The efficient generation of proof hints, which we call certificates, is described elsewhere [START_REF] Fouilhe | Efficient Generation of Correctness Certificates for the Abstract Domain of Polyhedra[END_REF], along with an experimental evaluation of the overall abstract domain. The main contribution of the work described here is the design of the link between the Coq frontend and the untrusted backend. It avoids the conversion and transfer of polyhedra. This makes the coupling between the frontend and the backend very loose. As a result, building other certificate-producing backends is easy and has no impact on the Coq frontend code. Complete freedom is given on the choice of data structures: a backend could use constraint or double representation for polyhedra. Furthermore, since the backend does not give formal precision guarantees, a backend could implement relaxations of domain operators [START_REF] Laviron | Subpolyhedra: a (more) scalable approach to infer linear inequalities[END_REF][START_REF] Sankaranarayanan | Efficient strongly relational polyhedral analysis[END_REF], trading precision for efficiency.

We also present a lightweight method to ensure the soundness of Ocaml code extracted from the Coq frontend, even when the backend has an internal state, or when the functional purity of the backend is not trusted. Last, we describe the architecture of the frontend as a collection of functors which extends a bare metal abstract domain interface. This approach makes the proof modular: it is simpler and more robust to change.

A certified interface of polyhedral abstract domains

Let us introduce an small imperative programming language, named PL. The syntax of PL programs is described on figure 1. Letter t stands for an affine term and c is a condition over numerical variables with the following syntax: Let us now sketch how to build a "sound-by-construction" static analyzer performing a value analysis for this toy language. This will also introduce our logical interface of abstract domains. For simplicity, we assume here that PL programs are annotated, by an untrusted analyzer, with candidate invariants p i where they are hard to infer: at the loop headers. Hence, we only have to prove the soundness of a postcondition computation, described on figure 1, which checks whether the candidate invariants are inductive. Given a precondition p, the postcondition p s of a statement s is computed using recursion on the syntax of s. If a candidate invariant cannot be shown to be inductive, it is replaced by ⊤, which is always safe. This happens when the candidate invariant does not include the postcondition of the loop body and may have two causes. Either the candidate invariant is not inductive, or the abtract domain used for checking inductiveness is not precise enough.

c ::= t 1 ⊲⊳ t 2 | ¬c | c 1 ∧ c 2 | c 1 ∨ c 2 with ⊲⊳ ∈ {=, =, ≤, ≥, <, >}. All numbers are rationals. s x := t s1; s2 if (c){s1}else{s2} while(c){s : pi} p s p[x := t] p s1 s2 p ⊓ c s1 ⊔ p ⊓ ¬c s2 pi ⊓ ¬c if p ⊑ pi ∧ pi ⊓ c s ⊑ pi ⊤ ⊓ ¬c otherwise
The postcondition computation relies on the operators of the abstract domain. Let us introduce them on the example of the abstract domain of polyhedra. Their Coq formal specifications are presented on figure 2. A polyhedron p encodes a formula i a i .x ≤ b i , where a i is a row vector of rational constants, b i is a rational constant and x is a column vector of numerical variables of the program. Its semantics, or concretization, is the predicate p defined as λm. i j a ij .m(x j ) ≤ b i , where m is a total map from variables to rationals representing a memory state. We omit the definitions of the semantics t of t and c of c, as they are standard.

-Polyhedron ⊤ corresponds to the predicate True.

-Polyhedron ⊥ corresponds to the predicate False.

-Polyhedron p ⊓ c over-approximates the conjunction of p and c (hence, the forward predicate transformer for guard). -Polyhedron p 1 ⊔p 2 over-approximates the disjunction of p 1 and p 2 (hence, the forward predicate transformer for join). -Given a term t and a variable x, polyhedron p[x := t] over-approximates the result of applying the forward predicate transformer for x := t on p . -Boolean p 1 ⊑ p 2 over-approximates the inclusion of predicates: if it is true, then p 2 is a logical consequence of p 1 .

Although we have omitted them here, the Coq code of the invariant checker needs to formalize the semantics of PL and prove that the reachable states of a PL program are soundly captured by the postcondition computation defined on figure 1. This relies on the abstract domain operations satisfying the specifications on figure 2.

These specifications are weak: they only enforce that the operators of the abstract domain perform safe over-approximations. They give no information on the precision of the results. Building an abstract domain satisfying these specifications is our focus in this paper.

⊤ m ¬ ⊥ m p m ∧ c m ⇒ p ⊓ c m p1 m ∨ p2 m ⇒ p1 ⊔ p2 m p m ⇒ p[x := t] (m[x := t m]) p1 ⊑ p2 ∧ p1 m ⇒ p2 m

Result certification of polyhedral abstract domains

While using Coq enhances the reliability of software, it sets a number of restrictions on the programs which can be reasoned about. First, the Coq programming language is restricted to pure functions that must be shown to terminate. The algorithms used by the abstract domain of polyhedra are complex to implement and these Coq requirements would have made their implementation even harder.

The most representative example is the simplex algorithm. Furthermore, Coq programs cannot use native machine arithmetic for computing. Instead, numbers are represented as lists of bits. The algorithms operating on polyhedra being arithmetic intensive, this suggests carrying out as much computation outside Coq as possible. Again, the simplex algorithm is the most representative example.

The arguments required to prove the correctness of the operators of the abstract domain of polyhedra make it convenient to offload much computation to an untrusted oracle and keep only a small amount of code to be proved correct in Coq. We back up this claim with some background on polyhedra.

Representing polyhedra for certification

A polyhedron can be represented in two ways: as a conjunction of constraints (i.e. affine inequalities) or as a set of generators, as illustrated on figure 3.

When working with generator representation, proving correctness of the polyhedral operations specified on figure 2 requires proving completeness results. Indeed, forgetting one vertex of the hypercube yields an under-approximation of this hypercube. Correctness of static forward analysis is not preserved through under-approximation, but through over-approximation.

Proving correctness of polyhedral operations in constraint representation is easier, as forgetting one constraint of the result produces a safe over-approxi- mation. The proof can be built incrementally by proving that each produced constraint includes the exact result.

p m ⇒ ∀a, p\x (m[x := a]) p (m[x1 := m(x2)]) ⇒ p[x1 ← x2] m

Expressing correctness as inclusions of polyhedra

The correctness of each operation reduces to inclusions of polyhedra. However, this reduction requires to first break the complex operations given in figure 2 into simpler ones, which compose a low level interface of polyhedral abstract domains. This low level interface has the same inclusion test and join operators as before. The guard is restricted to one affine constraint. There is no forward predicate transformer for assignment, but projection and renaming operators are provided, from which it can be built (see §6.3). The specifications for projection and renaming are given in figure 4. Polyhedron p\x results from the projection of p on the space of variables where dimension x has been removed. Renaming p[x 1 ← x 2 ] over-approximates the renaming of x 1 as x 2 in p. Variable x 2 is required to be fresh, but this precondition is not formalized in the Coq specification as it is not needed for our correctness proofs. A violation of this precondition may result in a precision bug, see §4.1.

Correctness of the operators of the low level interface can now be reduced to inclusions of polyhedra, with the exception of renaming. However, under the freshness precondition, renaming is a purely syntactic transformation. Three operations remain: the guard of a polyhedron p with an affine constraint c, the projection of a variable x from a polyhedron p and the join of two polyhedra p 1 and p 2 . Each constraint c ′ of their result must be shown to satisfy the inclusion properties specified below.

guard. i c i ∧ c ⊑ c ′ , with p i c i projection. p ⊑ c ′ (and x should have a nil coefficient in c ′ ) join. p 1 ⊑ c ′ and p 2 ⊑ c ′

Checking inclusion of polyhedra

The correctness of the complex low level operations reduces to inclusions of polyhedra. Farkas's lemma further reduces polyhedra inclusion to a linear programming problem on constraint representation of polyhedra. Below, we say that "constraint

a 1 .x ≤ b 1 syntactically entails a 2 .x ≤ b 2 " if and only if a 1 = a 2 and b 1 ≤ b 2 .
Farkas's lemma. A polyhedron p c i is included in a one-constraint polyhedron c ′ if and only if there exists λ i ≥ 0, such that i λ i .c i syntactically entails c ′ . Given λ, a vector of λ i , checking that p is included in c ′ is straightforward: build the linear combination c λ.p and check that c syntactically entails c ′ . This generalizes to a polyhedron p ′ c ′ j by supplying a vector of coefficients λ j for each constraint c ′ j of p ′ . The vectors λ j form a matrix Λ such that p s Λ.p and the constraints of p s entail those of p ′ syntactically. The matrix Λ can be used by a checker to validate the result of an operator: we call Λ an inclusion certificate. Nothing is proved when the check fails, however.

Core architecture of the abstract domain

Farkas's lemma makes result verification cheap. Moreover, it guarantees that producing a certificate to justify an inclusion property is always possible. This motivates the two tier architecture we have chosen for our polyhedral abstract domain.

The abstract domain is split in an untrusted Ocaml backend and a frontend which is developed in Coq. The backend performs most complex computations of the low level interface. Its interface is given on figure 5. The backend provides certificates of type cert that allow the frontend to produce certified results. Type poly is the internal representation of polyhedra used by the backend: it remains opaque for the frontend. The functions isEmpty and isIncl produce a certificate only when inclusion in ⊥ or in another polyhedron holds. Other operations produce both polyhedra and certificates, except for renaming where a certificate is not needed.

The communication protocol between the backend and the frontend is detailed in next section. Section 5 describes the formalization in Coq of the backend functions. Last, section 6 describes how complex polyhedra operations are built from the low level interface.

Using certificates as build instructions

Three polyhedra operators use a certificate from the backend and produce a polyhedron: the guard, join and projection operators. What we have presented Definition project ( pF , pB ) x := let ( pB ' , ce ) := Backend . project pB x in let pF ' := projectUsing ce pF x in ( pF ' , pB ') Fig. 6. The implementation of the projection in the frontend leads naturally to a pattern of algorithms for the frontend, which we illustrate for the projection operator. First, polyhedron p s Λ.p is built, using the certificate Λ provided by the backend. Syntactic entailment is then checked with the result p\x actually provided by the backend. An extra check is specific to the projection: verifying that x is free in p.

Checking syntactic entailment is actually unnecessary: p ′ can be used as a result of the projection operator. It satisfies the inclusion property, by construction. On top of sparing the entailment check, this approach removes the need for the backend to communicate its result to the frontend. The certificate is sufficient. This remark applies to the projection operator, as well as to the guard and join operators.

As a result, the operators follow a simpler pattern, illustrated for the projection operator on figure 6. The backend and the frontend both have their own representation of a polyhedron, which we call pB and pF, respectively.

That pB and pF represent the same polyhedron is an invariant property. An operator of the abstract domain consists in invoking the corresponding operator of the backend, thereby obtaining the backend representation pB' for the resulting polyhedron. The backend also produces a certificate ce, from which the frontend computes its representation pF' of the result of the operator, along with a proof that it is correct. This restores the synchronisation between the frontend and backend: pF' and pB' represent the same polyhedron.

The impact of bugs

Previous discussion makes the assumption that all goes well: the certificate is well-formed and yields a representation of the result computed by the backend. However, bugs might lurk in the backend, leading to incorrect results or erroneous certificates. Two possible effects can be observed by the user of the abstract domain.

-If the certificate is well-formed but yields a result different from that of the backend, synchronization is lost and the results built by the abstract domain are likely to be wildly over-approximated, yet correct. -If an ill-formed certificate (e.g. refers to nonexistent constraints) is output by the backend, the frontend will report a failure. Two failure modes are supported: abort or return a correct ⊤ result. Unless the backend aborts, the frontend returns correct results in all cases: soundness bugs in the backend induce precision bugs of the abstract domain. These bugs are uncovered using standard software engineering methods.

The certificate language

Inductive cert := | Implies : list ( positive * consCert ) → cert | Empty : linComb → cert | Bind : positive → consCert → cert → cert .

Fig. 7. Coq definition of polyhedron build instructions

The frontend builds correct by construction results using certificates provided by the backend. The type cert of the certificates is given in figure 7. We will describe the design of the certificates from the ground up on the example of a projection p\x for which the backend has produced a certificate Implies l.

From a high level of abstraction, Implies l is the sparse representation of a matrix Λ which defines the result p ′ Λ.p of the projection. In order to make the certificate compact, the constraints of p are identified by positive numbers and the descriptions of linear combinations, the type linComb, refer to constraints by their identifier. Identifier generation is handled by the backend: the frontend requests freshId pB when it needs a constraint identifier that does not appear in polyhedron pB. The frontend does not check the freshness of identifiers: as described in §4.1, invalid identifiers may result in precision bugs.

Type consCert describes the various ways to build one constraint of p ′ 1 . Its definition appears in figure 8. The Direct construct is the standard application of Farkas's lemma. For efficiency reasons, a backend may handle equality constraints specially, instead of representing them as pairs of inequalities. Two applications of Farkas's lemma are necessary to build an equality a.x = b from p i c i . One builds a.x ≤ b and the other builds a.x ≥ b. The equality follows from their conjunction and we introduced the SplitEq construct to handle this case.

The join operator requires a special construct, JoinCons. For each constraint c of the result of p 1 ⊔ p 2 , it must be shown that p 1 ⊑ c and p 2 ⊑ c. To this end, a JoinCons certificate contains one linear combination to build c 1 a 1 .x ≤ b 1 such that p 1 ⊑ c 1 and another for c 2 a 2 .x ≤ b 2 such that p 2 ⊑ c 2 . The frontend checks that a 1 = a 2 and then chooses c 1 as the resulting constraint if b 1 ≥ b 2 , or c 2 otherwise. If a 1 = a 2 , the certificate is considered ill-formed.

Type cert also provides a construct to build ⊥, as the result of a guard for example. An Empty l certificate is used for this purpose, where the linear combination l yields a trivially contradictory constraint, like 0 ≤ -1.

Let us motivate the last construct of type cert through a glimpse of the redundancy elimination behind a backend implementation of the guard p ⊓ c, with p i c i . Constraint c is rewritten using the equality constraints in p, so as to lower the number of variables involved. The result c ′ could then be involved in proving that the system of inequalities i c i ∧ c ′ hides an implicit equality e.

The new equality e could then used for further rewriting. Building a certificate in that setting is hard. The construct Bind j cc ce helps by introducing a new constraint resulting from the linear combination cc and giving it identifier j. The remainder of the inclusion certificate, ce, may then use it.

Formalizing the backend in Coq

Our abstract domain is split in two components: the frontend, which is developped in Coq, and the backend, which is written in Ocaml. In order to execute the code, the Coq frontend must be extracted to Ocaml code through Coq extraction mechanism. Extraction roughly consists in removing all the proofrelated information from a Coq development, as Ocaml type system is not powerful enough to represent it.

Once extracted, the frontend calls to the backend appear as function calls in the operators code. For the extraction to generate these calls, the backend functions must be declared to Coq as axioms. Let f be an external function of Ocaml type A → B. It is declared to Coq as a function f , of Coq type A → B and the extractor is instructed to replace calls to f with calls to f . Types A and B must be the extracted versions of A and B. The Ocaml compiler will report an error otherwise.

These declarations prevent the execution of the Coq development in Coq virtual machine: the body of the backend functions is not available to Coq. Furthermore, this process of linking certified Coq code to untrusted Ocaml code may lead to a number of serious pitfalls.

Inconsistency. An axiom like failwith : ∀B, string → B introduces inconsistency as it builds a proof of any B from a string. In particular, failwith False "" gives a proof of False. This pitfall is avoided by providing a model of axioms in Coq: a proof that their Coq type is inhabited.

Implicit axioms. Inductive type B (e.g. {x : Z | x < 5}) may be extracted into a strictly larger extracted type B (e.g. Z). This introduces an implicit requirement on f (i.e. its results are lower than 5) that Ocaml typechecker cannot ensure. For our frontend, we have thus carefully checked that Coq inductive types involved in backend functions are identical to their Ocaml extraction. Memory corruption. Our backend uses the GMP [START_REF]Free Software Foundation: The GNU Multiple Precision Arithmetic Library[END_REF] C library. A bug in GMP or its Ocaml frontend, zarith [START_REF] Miné | ZArith[END_REF], may corrupt arbitrary memory locations. However, it seems unlikely that such a bug breaks soundness silently.

Inductive consCert := | Direct : linComb → consCert | SplitEq : linComb → linComb → consCert | JoinCons : linComb → linComb → consCert .
Implicit purity axiom. Semantics of → are different in Coq and in Ocaml. In Coq, f is implicitly a pure function: hence ∀x, f x = f x is provable. On the contrary, f in Ocaml may use an implicit state such that, for a given x two distinct calls f x give different results. In other word, axiomatizing f as A → B in Coq introduces an implicit functional requirement: f is observationally pure.

Having an implicit state is allowed only if the effect of this implicit state remains hidden (e.g. for memoization). See [START_REF] Pottier | Syntactic soundness proof of a type-and-capability system with hidden state[END_REF] for details. However, it may be difficult to ensure that a backend has no observable side effects. In ours, a bug in GMP or zarith may break this requirement. Furthermore, our proofs do not rely on the purity of backend functions. The following describes the theory of impure computations we have formalized in Coq in order to declare the backend functions as potentially impure. This theory is inspired by simulable monads [START_REF] Claret | Lightweight proof by reflection using a posteriori simulation of effectful computation[END_REF], but from which we drop the notion of prophecy, because we are not interested in generating the backend from Coq.

May-return monads: a simple theory of impure computations

Impure computations are Coq computations that may use external computations in Ocaml. For any Coq type A, we assume a type ?A to denote impure computations returning values of type A. Type transformer " ?. " is equipped with a monad:

-Operator bind A,B : ?A → (A → ?B) → ?B encodes Ocaml "let x = k 1 in k 2 " as "bind k 1 λx.k 2 ". -Operator unit A : A → ?A lifts a pure computation as an impure one.

-Relation ≡ A : ?A → ?A → Prop is a congruence (w.r.t. bind) which represents equivalence of semantics between Ocaml computations. Moreover, operator bind is associative and admits unit as neutral element. Last, we assume a relation A : ?A → A → Prop and write "k a" to denote the property that "computation k may return a". This relation is assumed to be compatible with ≡ A and to satisfy the following axioms:

unit a 1 a 2 ⇒ a 1 = a 2 bind k 1 k 2 b ⇒ ∃a, k 1 a ∧ k 2 a b ?A A k1≡k2 k1=k2 k a k = a unit a a bind k1 k2 k2 k1
Fig. 10. A trivial implementation of the may-return monad

The theory of may-return monads is a very abstract axiomatization of impurity: it does not provide any information about effects of impure computations. However, as our frontend only cares about results of backend functions, this theory suffices to our needs. Hence, backend functions f of type A → B are simply axiomatized in Coq as f : A → ?B.

Our frontend is parameterized by an implementation of may-return monads: it does not depend on a particular model. Simple transformers over a global state have a denotation in the state monad defined in figure 9, using S as type of states. Even if building a model where any Ocaml computation is denoted is complex [START_REF] Birkedal | Step-indexed kripke models over recursive worlds[END_REF] and beyond the scope of this work, this gives us confidence in our frontend being correct when used with a side-effecting backend.

Extraction of impure computations

The may-return monad is useful to control Coq assumptions that would otherwise be left implicit. However, it is of no other practical interest and is removed at extraction time by providing the trivial implementation given on figure 10. The extractor inlines these definitions so that the monad has no runtime overhead.

The trivial implementation of the may-return monad is also used to provide a pure Coq interface to our abstract domain, by exposing that is equality. Although this partly puts the backend in the trusted computing base (TCB), this was actually required to plug our library as an abstract domain of the analyzer developed as part of the Verasco project [START_REF] Blazy | Formal Verification of a C Value Analysis Based on Abstract Interpretation[END_REF].

Backward reasoning on impure computations

Having introduced axioms for impure computations in section 5.1, we sketch below how we automate Coq reasonings about such computations, by using a weakest-precondition calculus programmed as a Ltac tactic.

First, we define in Coq an operator wlp A : ?A → (A → Prop) → Prop such that wlp k P ∀a, k a ⇒ (P a) expresses the weakest precondition ensuring that any result returned by computation k satisfies postcondition P .

For example, let us consider a Coq function g that first calls an external f returning a natural number of N and second, increments its result. We define g x bind (f x) λn.(unit n+1) and express the property that "g returns only strictly positive naturals" as the goal "∀x, wlp (g x) λn.n = 0". Our Ltac tactic simplifies this goal into a trivial consequence of "∀n : N, n+1 = 0". This tactic proceeds backward on wlp-goals, by applying repeatedly lemmas which are represented below as rules. It first tries to apply backward a decomposition rule: one for unit or bind below, or one for pattern-matching over some usual types (booleans, option types, product types, etc.). When no decomposition applies, the tactic applies Cut&Unfold. Actually, it tries to discharge the left premise using existing lemmas; if this fails, the definition of wlp is simply unfolded; otherwise, the goal is replaced using the right premise: the unfolding is thus performed with a lemma injection in hypothesis. In our Coq development, this tactic automates most of the bureaucratic reasoning on first-order impure computations. For higher-order impure computations (e.g. invoking a list iterator), equational reasoning is also needed.

Modular construction of the abstract domain

We have described in the last two sections a basic interface to the abstract domain of polyhedra. It is a restricted version of the interface described in section 2: the forward predicate transformer for assignment is missing, for example. The gap between the low level interface, closer to what the backend provides, and the fully-fledged interface, that our abstract domain offers to the user, is bridged entirely in the Coq frontend. The extra functionality is provided through the use of functors. Each functor takes an abstract domain and builds a richer one while lifting the proofs as necessary. This decomposition makes the proofs more manageable and modular. The overall architecture of the abstract domain is pictured on figure 11. The shaded left-hand side is the Coq frontend. Each of the pictured layers represents a functor. The untrusted backend stands on the right-hand side. While communication between the two is represented by arrows, it reduces to function calls in the extracted frontend code.

Building the guard operator

A first example of our modular construction of some abstract domain features is the guard operator presented in section 3. This operator p⊓c accepts an arbitrary propositional formula as constraint c. However, the backend guard operator takes only constraints of the form t ⊲⊳ 0, with ⊲⊳∈ {=, ≤, <} and t an affine term. The transformation from the more expressive guards to the basic guards is performed by the frontend by the following steps:

1. Negations are eliminated using De Morgan's laws on binary operators to push negation inwards, eliminating double negations, and taking the dual comparison on atomic formula. 2. Comparison = is rewritten as a disjunction of strict inequalities. 3. On Z, t 1 < t 2 is rewritten as t 1 -t 2 +1 ≤ 0. This increases precision of our polyhedra computations where all variables are in Q. 4. Disjunctions are over-approximated by joins. For a guard p ⊓ c, our algorithm performs a number of polyhedra operations that is linear in the number of operations in c. The functor which provides this extended guard operator to an abstract domain featuring only a basic one also contains the proof that the algorithm described above is sound.

Framing potentially constrained variables

Generating fresh variables has many applications for program verification: handling local variables, parameter passing during function inlining, implementing the forward predicate transformer for assignment, etc. However, our Coq specifications of abstract domains do not provide any information about the set of variables constrained by a polyhedron. Indeed, these loose specifications allow a modular management of fresh variables: in particular, certification of low-level operations presented at section 3 is not intricated with fresh variables handling.

Hence, this section and the next one introduce generic abstract domains (functors) gluing additional data about constrained variables to the value of an underlying abstract domain. We certify these functors by expressing the invariant of theses additional data through the concretization function of the newly introduced domain. We now illustrate the necessity of this trick on a functor that simply frames the variables constrained by an abstract value (e.g. a polyhedron).

Formally, if F is a set of variables, we note m ≡ F m ′ if and only if memories m and m ′ coincide on F . Then, given a polyhedron p, we say that F frames p if and only if ∀m 1 ∀m 2 , m 1 ≡ F m 2 ⇒ ( p m 1 ⇔ p m 2 ), and that x is free in p if and only if {x ′ | x ′ = x} frames p. An operator bnd(p) can then be defined such that {x | x ≤ bnd(p)} frames p (variables are represented by positive integers). This operator returns an upper bound β on the variables constrained by p: we say that β bounds p. These definitions also apply to conditions and terms.

Operator bnd is provided by a new abstract domain P bnd wrapping each element p of the underlying domain into a pair (p, β) such that β bounds (p, β). Operations of P bnd are given figure 12.

Naive definitions of P bnd fail to make provable the correctness of operations: property "∀(p, β) ∈ P bnd , β bounds p" may not be preserved by the operations of figure 12. For example, let us consider (p 1 , β 1 ) ⊔ bnd (p 2 , β 2 ). We expect (p 1 ⊔ p 2 , max(β 1 , β 2 )) to be a good candidate implementation, since if β 1 bounds p 1 and β 2 bounds p 2 , then max(β 1 , β 2 ) bounds both p 1 and p 2 . However, max(β 1 , β 2 ) may not bound p 1 ⊔p 2 , as this somewhat contorted, yet correct, implementation of p 1 ⊔ p 2 shows when x is chosen above the bound 1:

p 1 ⊔ p 2 x ≤ 0 if p 1 = p 2 = ⊥ ⊤ otherwise
Our solution consists in keeping the definitions given in figure 12, but changing that of (p, β) so that it implies the property "β bounds (p, β)". Given a concrete memory m, we impose that variables above β are free in (p, β) m by quantifying over any abstract memory m ♯ that results from the arbitrary update of m on these variables:

(p, β) m ∀m ♯ , m ♯ ≡ { x | x≤β } m ⇒ p m ♯

Assignment with buffered renaming

The P bnd functor can be used to over-approximate the forward predicate transformer of assignment. Indeed, it allows to introduce an auxiliary fresh variable x 0 which names the value of variable x before the assignment:

p[x := t] (p[x ← x 0 ] ⊓ x = t[x ← x 0 ]) \x 0 where x 0 max(bnd(t), bnd(p))+1
However, our abstract domain uses the P := functor described below instead, because it performs a lower amortized number of polyhedra renamings.

Functor P := makes it possible to express relations between memory states in the intermediary computations of the operators. This achieved by duplicating the set of variable names: each variable x can be represented as x@0 or x@1. Of these two representatives of x, the concretization imposes that exactly one refers to a concrete memory cell. Similarly to what is done in P bnd , the other representative is arbitrarily updated in abstract memory m ♯ . The concretization involves a function σ that associates its current representative to each variable. It also bnd involves a function π that associates concrete x to both abstract variables x@0 and x@1, for all x.

((p, β)) β ⊤ bnd (⊤, 1) ⊥ bnd (⊥, 1) (p1, β1) ⊔ bnd (p2, β2) (p1 ⊔ p2, max(β1, β2)) (p, β) ⊓ bnd c (p ⊓ c, max(β, bnd(c))) (p1, β1) ⊑ bnd (p2, β2) p1 ⊑ p2 ∧ β1 ≤ β2
(p, σ) m ∀m ♯ , m

♯ ≡ { x ♯ | x ♯ =σ(π(x ♯ )) } m • π ⇒ p m ♯
In the P := functor, assignment to x switches the representative of x, instead of renaming the variable in the underlying polyhedron as with P bnd . Renamings from assignments are buffered until joins or inclusions, where they may eventually be performed (only when representatives of identical variables need to be unified). Furthermore, two successive renamings on the same variable in the buffer annihilate (by involution of representative switch).

This functor could be extended so as to buffer projections, which can then be reordered to get smaller intermediate results (in terms of size of representation). The decision to apply projections is delegated to the backend. In this version, the functor introduces a unique representative at each assignment: a kind of SSA form is thus computed on-the-fly in the abstract domain. This extension is not implemented yet.

In conclusion, our modular treatments of assignment depart from [START_REF] Besson | Result certification for relational program analysis[END_REF]: this results in more manageable proofs. In [START_REF] Besson | Result certification for relational program analysis[END_REF], projections are systematically delayed until inclusion tests: we believe that the choice of when to apply projections should be delegated to the backend.

Conclusion

We presented one solution to prove the correctness of an implementation of the abstract domain of polyhedra using the Coq proof assistant. In this setting, correctness reduces to inclusions of polyhedra which, through Farkas's lemma, makes a posteriori verification of results a convenient approach. As a result, our domain is composed of an untrusted backend, to which most of the complex computations are offloaded, and a Coq frontend which validates the results produced by the backend. This work makes two main contributions.

On one hand, we consider the implicit requirements set when linking certified code to untrusted external code in the Coq proof assistant. This delicate issue might be carefully considered by certification authorities for e.g. avionics. We partly address it through a lightweight method for declaring the backend functions to Coq in such a way that the proofs remain trustworthy even when the backend is not functionally pure.

On the other hand, we show how communication between the frontend and the backend can be reduced to certificates, which serve as build instructions for the frontend. The certificate language induces a low coupling between the frontend and the backend: the latter could implement relaxations of some operators [START_REF] Laviron | Subpolyhedra: a (more) scalable approach to infer linear inequalities[END_REF][START_REF] Sankaranarayanan | Efficient strongly relational polyhedral analysis[END_REF] or use entirely different data structures without requiring changes to the frontend. Although it does not make abstract domain development easier, our approach reduces the impact of bugs.

The complete domain further distinguishes itself from previous work by integrating certificate generation to the backend and by a more modular proof architecture. Experiments shows that it has comparable performance to noncertifying state-of-the-art implementations [START_REF] Fouilhe | Efficient Generation of Correctness Certificates for the Abstract Domain of Polyhedra[END_REF].

The complete code is available on the Web, along with a demonstration application, from www-verimag.imag.fr/~boulme/vstte2014.html.

Fig. 1 .

 1 Fig. 1. Syntax and postcondition computation of PL

Fig. 2 .Fig. 3 .

 23 Fig. 2. Correctness specifications of our main polyhedral operations

Fig. 4 .

 4 Fig. 4. Correctness specifications of the projection and renaming operators

Fig. 5 .

 5 Fig. 5. Ocaml interface of the backend

Fig. 8 .Fig. 9 .

 89 Fig. 8. Coq definition of constraint build instructions

  wlp k1 λa.(wlp (k2 a) P ) wlp (bind k1 k2) P Cut&Unfold wlp k P1 ∀a, k a ∧ P1 a ⇒ P2 a wlp k P2

Fig. 11 .

 11 Fig. 11. Overview of the domain architecture

Fig. 12 .

 12 Fig. 12. Main operators of P bnd

⋆ This work was partially supported by ANR project VERASCO (INS 2011).

Acknowledgements. We would like to thank Michaël Périn and David Monniaux for their continuous feedback all along this work. We also thank the members of the Verasco project for their motivating interaction.


Nawfal Blal

Loïc Daridon

Yann Monerie

Stéphane Pagano

A micromechanical-based model for the ductile cohesive-volumetric damage

For numerical purposes, the cohesive-volumetric finite element method has gained much popularity in crack and damage simulations. A micromechanical-based damage model is here proposed for such ductile "cohesive-volumetric" media. The proposed model defines theoretical and practical criteria for the calibration of cohesive zone models (CZMs) in these simulations.

The studied medium is made of a hardening matrix (Hencky plasticity) containing penny-shaped cohesive inclusions (tractionseparation laws). The spatial distribution of the cohesive inclusions can fit a prescribed finite element discretization. The overall elastoplastic and damageable behavior is derived using a non-linear homogenization technique (variational approach of P. P. Castañeda, [START_REF] Ponte Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF]). The proposed model: i/ can be applied whatever the shape of the cohesive law contrary to what is proposed in the literature, and ii/ has the capacity to exhibit the influence of the triaxiality loading rate on the overall ductile damage properties.

For the case of a perfect-plastic bulk medium, a closed-form expression of the macroscopic potential is obtained. And whatever the volumetric plastic behavior, direct relationships between the local cohesive parameters and the overall material properties are developed. These relationships depend on: 1/ the mesh size and morphology, 2/ the applied triaxiality loading rate and 3/ the macroscopic material properties (maximal stress, failure energy, etc). More particularly, it is shown that the cohesive parameters are triaxiality-dependent. This dependency is consistent with previous results available in the literature and based on numerical or experimental studies (e.g. [START_REF] Siegmund | Prediction of the work of separation and implications to modeling[END_REF]).

Figure 1 .

 1 Figure 1. Effect of the triaxiality loading rate, T, on the overall behavior: normalized overall stress vs overall strain

3rd INTERNATIONAL WORKSHOP ON PHYSICS BASED MATERIAL MODELS AND EXPERIMENTAL OBSERVATIONS 2-4 JUNE, 2014, Çeşme/Izmir -TURKEY A micromechanical-based model for the ductile cohesive-volumetric damage

	Laboratoire de Mécanique et Génie Civil, LMGC, Université Montpellier 2, CNRS, CC 048 Place
	Eugène Bataillon, 34095 Montpellier
	Laboratoire de Micromécanique et d'Intégrité des Structures, MIST Laboratory, IRSN-CNRS-
	Université Montpellier 2
	{nawfal.blal,loic.daridon,yann.monerie,stephane.pagano}@um2.fr
	Key Words: Micromechanics, Plasticity, Cohesive Zone Models, Homogenization

N. BLAL, L. DARIDON, Y. MONERIE and S. PAGANO