Alexandre Perrot

Romain Bourqui

Nicolas Hanusse

Fr Éd Éric Lalanne

David Auber

Large Interactive Visualization of Density Functions on Big Data Infrastructure

Keywords:

Point set visualization is required in lots of visualization techniques. Scatter plots as well as geographic heat-maps are straightforward examples. Data analysts are now well trained to use such visualization techniques. The availability of larger and larger datasets raises the need to make these techniques scale as fast as the data grows. The Big Data Infrastructure offers the possibility to scale horizontally. Designing point set visualization methods that fit into that new paradigm is thus a crucial challenge. In this paper, we present a complete architecture which fully fits into the Big Data paradigm and so enables interactive visualization of heatmaps at ultra-scale. A new distributed algorithm for multi-scale aggregation of point set is given and an adaptive GPU based method for kernel density estimation is proposed. A complete prototype working with Hadoop, HBase, Spark and WebGL has been implemented. We give a benchmark of our solution on a dataset having more than 2 billion points.

INTRODUCTION

Point set visualization is without any doubt a key weapon in the arsenal of a data visualizer. The John Snow Map of deaths due to the 1854 London cholera epidemic is typically used to demonstrate the efficiency of point set visualization. By visually detecting a higher density of points (i.e., concentration of death) near a pump in the Soho district, John Snow was able to find that the source of the epidemic was infected water from that pump. The Figure 1.a shows a point set visualization of the John Snow dataset. If we zoom out the Figure 1.a (Fig. 1.c), one can see that due to point occlusion, it becomes harder to detect regions with high density of points. In [START_REF] Kindlmann | An algebraic process for visualization design[END_REF], Kindlman et al. have proposed a method to evaluate the efficiency of a visualization technique. They formalize the fact that a change in the visualization must be proportional to a change in the data. According to the measures they introduce, standard point set visualization technique is not efficient. Indeed, even if the data does not change between Figure 1.a and Figure 1.c, the visualization varies significantly.

To solve this occlusion problem, many techniques have been set up. The core idea of all of them is first to count the number of points c i, j that are drawn at the same location i, j on the screen. Then one point p i, j is rendered for each unique location i, j and c i, j is used to determine the color of p i, j . In the literature, the term heatmap is often used to describe that technique. More elaborate heatmap visualizations construct a continuous function by applying a convolution on c i, j . This continuous function is called a density function. In the Figure 1.c, one can see that density function eases identification of dense regions on the map, even when there is plenty of point * emails: {aperrot,bourqui,hanusse,flalanne,auber}@labri.fr heatmap, even at low resolution, dense regions are easy to detect. occlusion (Fig. 1.d). Thus, as discussed in [START_REF] Kindlmann | An algebraic process for visualization design[END_REF], density function visualization helps to solve the data change vs visualization change problem illustrated by Figure 1. Density functions have been used in a large number of applications. For instance, in [START_REF] Špakov | Visualization of eye gaze data using heat maps[END_REF][START_REF] Cutrell | What are you looking for?: an eye-tracking study of information usage in web search[END_REF], it enables the analysis of eye-tracking data. Analysis of geolocalised data is also often done with density function visualization [START_REF] Hotta | Online geovisualization with fast kernel density estimator[END_REF][START_REF] Fisher | Hotmap: Looking at geographic attention[END_REF][START_REF] Liu | imMens: Real-time visual querying of big data[END_REF]. In this case, the heatmaps are displayed on top of a geographical map.

The main motivation of our research was to provide a fully working solution for standard Big Data Infrastructure (BDI). The Big Data ecosystem enables to process big datasets, but introduces some constraints. From an algorithmic perspective, it is preferable to scan the data a constant number of times. Moreover, the data has to be stored in a distributed way. In order to bound latency, the entire dataset must not be transferred to the location where the visual analysis is done. A pre-computation can be done in order to allow aggregate or part of the dataset to be transferred on the user's computer. For all these reasons, visualization of Big Data is best done in a web browser. This enables visualizations to be stored on a server and sent on demand.

The contribution of this article is a complete fully usable architecture for heatmap visualization at large scale. All the constraints mentioned above have been taken into account to allow its direct integration into modern Big Data Infrastructure. To sum up, for a point set of size n stored within a distributed system of k nodes, there is a preprocessing in time O( n×log(n) k × i) when the points are uniformly spread on the plane, where i is the number of levels of detail desired. It enables to guarantee the approximation of the density function and the quality of visualization. Such a preprocessing step guarantees a constant size message transfer for any interaction event (zoom in/out or pan). Then, the client computes the visualization in O(p) time, with p the number of pixels on the screen. In a theoretical model of computation, the interaction with our system is in constant time. In practice, our experiments show the relevance of our proposition.

The paper is organized as follows. First, we give an overview of the most related work. Second, we present an overview of the architecture we have set up. We then provide the algorithm details. We conclude by a discussion supported by the benchmark we have done with datasets having more than a billion points.

PREVIOUS WORK 2.1 Classical density visualization

The main technique used in the literature to compute the density function of a point set is called Kernel Density Estimation (KDE) [START_REF] Silverman | Density estimation for statistics and data analysis[END_REF]. Introduced independently by Rosenblatt [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF] and Parzen [START_REF] Parzen | On estimation of a probability density function and mode[END_REF], KDE was first used to estimate the density function of a random variable. Using a kernel function, one can generate a continuous function from a discrete point set. Follow the equation of the KDE and the equation of the Gaussian kernel, the most commonly used kernel:

KDE(x) = 1 n n ∑ i=1 w i f σ (dist(x, x i )) f σ (x) = 1 σ √ 2π e -x 2 2σ 2
Unformally, a point of weight w spreads a weight w f σ (d) for any location at distance d. The bigger σ is, the further a point will spread its weight. Visually, this represents the amount of "smoothing" of the density function. In the Figure 2 one can see the density function computed using the black points as the original point set. In this Figure, we used a gaussian kernel.

One can find lots of applications of kernel density estimation in the literature. The graph splatting technique [START_REF] Van Liere | Graphsplatting: Visualizing graphs as continuous fields[END_REF] introduced by van Liere and de Leeuw uses it to visualize node density in a graph. In that paper, the density function is called a "splat field" and is computed using the OpenGL texture rendering engine and the accumulation buffer. In [START_REF] Lampe | Interactive visualization of streaming data with kernel density estimation[END_REF], Lampe and Hauser have proposed another way to compute Kernel Density Estimations with OpenGL, using shaders instead of textures. In both of these works, the density function is only an estimation obtained by limiting the distance at which a point contributes to the density function.

Based on those works, Zinsmaier and Brandes proposed a levelof-detail visualization for large graphs [START_REF] Zinsmaier | Interactive level-of-detail rendering of large graphs[END_REF]. In that work, the KDE technique is used to show nodes density. The sources and targets of edges are then moved to the nearest local maximum of the node density function using a hill-climbing algorithm. A similar idea has been used for edge bundling with Kernel Density Estimation by Hurter et al. [START_REF] Hurter | Graph bundling by kernel density estimation[END_REF].

When using KDE to visualize density function, one needs to compute the density value for every pixel of the visualization. The complexity of this operation is O(n * p), with n being the number of points in the dataset and p the number of pixels in the visualization. For instance, computing the exact density function of 1 million points on a 2 mega pixels screen (1920x1080) requires 2 trillion operations. This complexity issue makes KDE a candidate for parallel computing. In addition to the GPU techniques and approximation method cited earlier, Michailidis and Margaritis [START_REF] Michailidis | Parallel computing of kernel density estimation with different multi-core programming models[END_REF][START_REF] Michailidis | Accelerating kernel density estimation on the gpu using the cuda framework[END_REF] compared KDE implementations with different parallel computing frameworks and Łukasik [START_REF] Łukasik | Parallel computing of kernel density estimates with mpi[END_REF] compared parallel implementations with MPI. Such methods are essential to enabling interactive visualization on large point sets.

Enabling density visualization for big datasets

In the following we present methods that enable to reduce the size of the original point set before starting the evaluation of the density. These techniques can be used as a pre process of techniques presented above.

Sampling: Sampling is a technique that removes data points from the dataset in order to reduce its size. It is successfully used to reduce overplot in scatterplots. Mayorga and Gleicher [START_REF] Mayorga | Splatterplots: Overcoming overdraw in scatter plots[END_REF] used it in conjunction with KDE to help visualize density in cluttered scatterplots. Chen et al. [START_REF] Chen | Visual abstraction and exploration of multi-class scatterplots[END_REF] also used an adaptive hierarchical resampling scheme to abstract multi-class scatterplots.

Clustering: Contrary to sampling, clustering merges points, instead of removing them, and represents the merged data by a new point. Clustering is often repeated to produce a hierarchical representation of the data called a cluster tree [START_REF] Abello | Ask-graphview: A large scale graph visualization system[END_REF][START_REF] Auber | Interactive refinement of multi-scale network clusterings[END_REF][START_REF] Delort | Vizualizing large spatial datasets in interactive maps[END_REF][START_REF] Elmqvist | Hierarchical aggregation for information visualization: Overview, techniques, and design guidelines[END_REF]. One can select an antichain on this tree to serve as dataset abstraction. The selection can depend on view parameters such as zoom and position, enabling interactive visualization.

Binned Aggregation : In ImMens [START_REF] Liu | imMens: Real-time visual querying of big data[END_REF], Liu et al. avoid the cost of the KDE computation by using binned aggregation. This process groups data points into predefined "bins" by partitioning the data space. The bin in which a data point is grouped does not depend on other points. It can thus be computed in constant time. This technique has the advantage to be fast to compute and can easily be computed in parallel, since the binning process is independent for each point. This technique can be considered as a case of geometrical clustering. The disadvantage of this technique is that it creates edge effects at the border of the bins. Li et al. [START_REF] Li | Interactive visualization of high density streaming points with heat-map[END_REF] applied KDE to binned points to produce a multilevel heatmap visualization.

Out of core visualization: In order to scale up to bigger datasets, it is essential to use out of core architectures. This kind of architecture is useful for multilevel visualization, as demonstrated by Hadwiger et al. [START_REF] Hadwiger | Interactive volume exploration of petascale microscopy data streams using a visualization-driven virtual memory approach[END_REF]. Moreover, the Map-reduce framework is useful for distributing existing algorithms [START_REF] Vo | Parallel visualization on large clusters using mapreduce[END_REF]. Using a specialized system, it enables to process spatial datasets and render images of heatmaps [START_REF] Eldawy | Shahed: A mapreduce-based system for querying and visualizing spatio-temporal satellite data[END_REF]. Similarly, Meier et al. [START_REF] Meier | Heattile, a new method for heatmap implementations for mobile web-based cartographic applications[END_REF] have designed an application where the data is stored on a server and delivered on the fly. There are no pre-rendered images and point density function are computed on the client side. Every datapoint is stored with a tile-based index for every zoom-level, enabling fast retrieval of a portion of the map. However, as the authors themselves emphasize, their method is not efficient for large datasets.

None of these methods provide guarantees both toward big data constraints and quality of visualization for large datasets. In the following we propose a technique that shares similar ideas and that enables to scale to datasets with billions of elements, while conserving guarantees toward the quality of the visualization.

PROCESS OVERVIEW

The interactive visual exploration of density functions at large scale requires using distributed storage and distributed computation as well as GPU intensive computation. In this section we present an overview of our technique (see Figure 3). The technical details for each step are given in the following sections. [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], we obtain a similarity score of 0.959, which corresponds to a "good" rating, according to Zinner et al. [START_REF] Zinner | Impact of frame rate and resolution on objective qoe metrics[END_REF]. Image (c) shows the result of another merging pass. Points B' through F' have been merged into A", whose weight is now 10. Differences in the density visualization become visible because the grouping distance is high in respect to the screen space and to the σ used for the KDE.

Figure 3: Overview of the system's different steps. Each step is presented with its role and the technology we used. The data stream is first stored in a distributed file system before aggregation. The aggregation is realized on a distributed computing cluster. The result of the aggregation step is stored in a distributed NoSQL database.

A front server provides access to the stored aggregation via a REST interface. The client web application requests tiles from the server and render them using an adaptive rendering algorithm. Thanks to the aggregation step, the amount of data transferred is bounded.

Our solution is based on the use of a multi-scale approximation of density functions. For that purpose, we compute several levels. The level max is the entire point set and level i-1 is obtained by merging points of level i . Thus, the more aggregation steps we do, the lowest the number of points in the final level, level 0 . This multiscale approach enables us to guarantee the quality and processing time of our approximate density function visualization.

To work at large scale all the data is stored into a distributed file system and we designed a distributed algorithm to automatically generate aggregated levels. Combining data distribution and distributed computing enables to run computations directly on the computers containing the data. Even if they are aggregated, the generated levels can be very large, thus we also store them into a distributed filesystem through a NoSQL database. Distributed storage as well as distributed computation are the key features that enable our solution to scale to large datasets.

Once the levels have been computed, the next step is to compute and display the approximate density function on the screen. In our solution that stage is achieved on the client side. That operation requires computation of the density functions for each pixel on the screen and mapping that density function to a color scheme. It necessitates to transfer data from the distributed storage to the client side. To prevent the transfer of an entire level to the client, every level is split into tiles, which are then stored into a highly accessible distributed database. Only requested tiles are transferred to the client. To get a similar data size per tile the size of tiles are adapted according to the level of aggregation (i.e. the more abstract a level is, the fewer tiles it has). That property guarantees the size of packets transferred on the network and thus smooth interaction according to network bandwidth.

On the client side, we use a Least Recently Used cache system to store tiles of point set transferred during the exploration. For the computation of the density function and the color mapping, we designed a multiscale approximation algorithm that can be fully implemented on the GPU. One of our requirements was to have a lightweight client interface. Therefore, all the design choices for our algorithm have been made to ensure the possibility to implement it on standard recent web browsers.

GENERATION OF LEVELS OF AGGREGATION

Starting from the original dataset, we compute several abstractions of the dataset using the Approximate Kernel Density Estimation introduced by Zheng et al. [START_REF] Zheng | Quality and efficiency for kernel density estimates in large data[END_REF].

Let D be a square subspace of R 2 . Given the kernel f σ and P the original point set, an Approximate Kernel Density Estimation (AKDE) Q is a weighted point set, such that the KDE on Q approximates the KDE on P. The distance between P and Q is defined as follows :

KDE P -KDE Q = max x∈D |KDE P (x) -KDE Q (x)|
Q is said to be an ε-approximation of P if the distance between P and Q is less than ε. We also want |Q| to be smaller than |P|, so that computing the density function on Q is less time consuming than on P.

This technique enables to reduce the number of points for a level of detail, while generating an almost identical density function. Figure 2 shows density function visualization for a point set and two approximations. We use the Grouping Selection (GS) technique, presented by Zheng et al. in [START_REF] Zheng | Quality and efficiency for kernel density estimates in large data[END_REF] to generate our AKDE. The GS technique groups points closer than a threshold distance d and generates a set of weighted points. According to the authors, the error using this technique is bounded by d/σ , with d the grouping distance and σ the bandwidth (representing the standard deviation) of the kernel used for KDE. Therefore, we can guarantee and control how our visualization approximates the original point set.

Finding points on a surface distant of at least d is related to the circle packing problem [START_REF] Goldberg | The packing of equal circles in a square[END_REF]. Given d and a finite square of area S, the number of circles with diameter d that can fit inside this square without intersections is bounded, naively by S d 2 . This is the same as determining points distant of at least d. Therefore, the number of points generated by the GS technique is bounded, independently of the size of the starting point set. Using results on the circle packing problem, it is possible to determine a grouping distance that ensures a specific upper bound on the size of the point set generated by GS.

We implement the GS technique by using the canopy clustering algorithm [START_REF] Mccallum | Efficient clustering of high-dimensional data sets with application to reference matching[END_REF]. This algorithm groups points that are at distance at most d from each other. The canopy clustering algorithm. The distance to the current point is checked for every canopy. This result in a complexity of O(n 2 ) for the pseudo code presented here. A better complexity can be achieved using a spatial index data structure, such as a Quad Tree or a KDTree. With those data structures, the complexity is O(n * log(|canopies|)).

Canopies of level i and the threshold distance d i

In the following, we build a hierarchical set of canopy centers such that level 0 ⊇ level 1 . . . ⊇ level i . d 0 is the threshold distance of level 0 defined so that n 0 = |level 0 | is the maximal number of points (and message size) to transmit to the lightweight client. Depending on the context, we aim at upper bounding n 0 by a small constant (from 100 to 1000). To provide this guarantee, it is easy to show that taking any

d 0 ≥ 2A n 0 √ 3
is enough assuming that A corresponds to the number of pixels of the screen. The levels of detail of the heatmap correspond to zoom levels in the visualization. Each zoom level corresponds to a 2 i zoom factor, that is the distance between pixels for the same pair of points seen on the screen is doubled each time we zoom in. The surface of a zoom level on screen is therefore A i = A 0 /4 i . Thus, taking d i = d i-1 /2 guarantees to have at most n 0 visible points among n i ≤ min (n, 4 i n 0 ).

In order to get an efficient and quick computation, we start by computing level i , the initial threshold d i can be chosen arbitrarily depending on the maximal precision we want. The set of canopy centers of level i will then be used as the entry of the computation of level i-1 taking d i-1 = 2d i . The smaller the input of the reducers, the quicker the global computation. Ideally, we aim at getting a theoretical time complexity of O(n/k) but this turns out to be difficult whenever a well-known algorithm is adapted into a MapReduce implementation.

The Mahout big data algorithm library1 provides a MapReduce implementation for the canopy clustering. In this implementation, the dataset is randomly partitioned between the mappers. Each mapper runs a local canopy clustering on its subset of the data. The resulting canopies are then transferred to a single reducer. This reducer applies another local canopy clustering. This implementation is well suited in cases where the set of canopies resulting from the mappers is small enough to fit on a single machine. However, due to the necessity to generate very abstract levels as well as detailed ones, our system must be able to cope with a very large (although bounded) number of canopies. Therefore, we could not use the Mahout implementation of the canopy clustering for our system. We thus designed a new implementation. Our implementation of the canopy clustering no longer uses a single reducer and thus, our approach can scale horizontally.

The random partitioning used in Mahout can lead to calculate the distance between points that are very far apart. We propose to use a more specialized data partition guaranteeing that the n i points are scanned at most twice. We start from a partition of data into strips that cuts the dataset along one of the dimensions. Let h = 4d be the height of a strip, if two points are in the same strip, their distance along one of the dimensions is at most h. This reduces the distance calculations to a more useful subset. In both of our jobs, the mappers will be responsible for the data partitioning and the reducers will execute the effective clustering process. One trick is that a reducer is in charge of a random set of strips. Even if some strips contain much less than n i /h points, the reducer will receive at When the strips are merged, as shown on image (b), it creates a zone called "conflict zone", whose height is 2d. This zone is colored in grey. In this zone, two points coming from the two strips can be closer than d. Image (c) shows the partitioning used in the second job. In order to be able to process the entirety of the conflict zone in a single strip, the strip limit is shifted downwards by d. As illustrated on image (d), the conflict zone is then totally contained in the upper strip. The conflict can be detected and solved. One of the points will be removed from the canopy set during the second job.

most O(n i /k) whenever the data distribution is close to be uniform.

First job: In the first job, the mappers receive points from the dataset. For each point, the index of the strip containing it is computed. Then the mapper emits a key value pair for each point, where the key is the strip index and the value is the point itself. Every record with the same key is then transmitted to the same reducer. Each reducer thus has every point inside a given strip. The canopy clustering is then executed inside each strip. We use a QuadTree to reduce the complexity, as suggested in algorithm 1. The pseudo code for the Map and Reduce phases of the job can be seen in algorithm 2. The Figure 4 shows how data is partitioned, transferred and assigned to reducers. After this process, the points are marked according to the result of the clustering. Second job: Once the first job is completed, the canopy clustering has been executed inside of each strip. However, as shown on Figure 5, when the strips are put back together, a point A close to the border of its strip can be closer than d from a point B in the neighbouring strip. In this case, both A and B are closer than d from the limit of the strips. We can define the zone which contains every such point. This zone is located at the border of adjacent strips and extends to a distance d from either side of the border. We call this zone "conflict zone".

The second Map Reduce job of our implementation aims at executing the canopy clustering on the conflict zone. We retain the idea of strips for data partitioning. However, the strip limit is shifted downwards by d, as shown on Figure 5. This way, the conflict zone is now contained in a single strip. This will enable to execute the canopy clustering on this area inside of each reducer.

This time, it is not necessary to execute a full canopy clustering Only canopies inside the conflict zone need to be treated. Using the QuadTree, we check if a canopy is already closer than d. If not, we add this canopy to the QuadTree. The data partition technique used in this implementation can induce load-balancing issues if the data points distribution is not uniform enough. It is then possible for a single strip to contain many more points than the others. However, when applying this algorithm recursively to generate the levels, only the first one can suffer from this issue. Once the algorithm has been executed a first time, the point set has a bounded size. This property is conserved when dividing it into strips, thus guaranteeing a maximum number of points per strip.

KDE RENDERING

The final step of our system is the rendering of the density function. This part happens on the client side after the transfer of the visible tiles. Due to the complexity of KDE calculation, approximations of the computation by limiting the kernel size have been developed in previous works [START_REF] Lampe | Interactive visualization of streaming data with kernel density estimation[END_REF][START_REF] Zinsmaier | Interactive level-of-detail rendering of large graphs[END_REF]. Here, we compute an exact KDE to be able to interactively change the kernel's size and type.

We use a simple technique that enables us to guarantee a constant rendering time on a given client at the expanse of KDE precision. The main idea of our technique is to adapt the rendering resolution according to the power of the GPU. The KDE is rendered on a smaller resolution when there are too many points, ensuring that the total number of operations needed is constant. The resulting smaller resolution KDE is then stretched to the original screen's resolution.

We introduce a threshold N of points that can be rendered at the original screen resolution in a small enough time to maintain satisfactory framerate. When we need to render more points than N, we will render the KDE on a smaller resolution.We define m = n N , with n the number of points to render. Instinctively, m measures the difference of magnitude between the data to render and the ac-ceptable amount. We want to render a smaller KDE on p = p m pixels. In order to keep the screen's aspect ratio, we take w = w √ m and h = h √ m , with w and h respectively the width and height of the screen. We thus have p = w * h . This interrelation is illustrated in Figure 6.

Rendering the n points costs p * n = p * n n N

. Using the properties of the ceiling function, it is easy to show that p * n n N ≤ p * N, so the complexity of the smaller KDE rendering is O(p * N). This only depends on the display resolution and the threshold constant. Considering that these values only depend on the client computer, one can consider that, for a given computer, our rendering algorithm runs in O(1).

IMPLEMENTATION DETAILS

The techniques and algorithms we describe here could be implemented using several technologies. We provide here an exhaustive list of the technologies we used. This constitutes a working set for a proof of concept. Figure 3 shows the complete system.

In the domain of Big Data, Hadoop is a well established standard and covers the majority of our needs. We use the Hadoop Distributed File System (HDFS) as raw data repository before processing. There are several implementations of the Map Reduce paradigm on top of Hadoop, via its resource manager YARN (Yet Another Resource Negotiator). The main one is Hadoop MapRe-duce2 (HMR2), which mainly uses disk based I/O between the computations. We used another Map Reduce implementation called Spark, which is also built on top of YARN and HDFS. Spark favors in-memory computation instead of the disk-based I/O used by HMR2. Spark is based on the principle of Resilient Distributed Datasets (RDDs) [START_REF] Zaharia | Resilient distributed datasets: A fault-tolerant abstraction for in-memory cluster computing[END_REF] to distribute computations closest to the data. In addition to the Map and Reduce operations, Spark provides several other operations optimized for distributed computing. This flexibility allows Spark to be easier to use than HMR2. In practice, we noted speed gains of order 100 from using Spark instead of HMR2.

The size of the data generated by our clustering algorithm demands a NoSQL distributed database. We chose to use the Hadoop database, HBase, which is an open source implementation of Google's BigTable [START_REF] Chang | Bigtable: A distributed storage system for structured data[END_REF] database system.

The link to the client application is made using a front server implemented using a JavaEE servlet. The task of this front server is to request data from HBase, and transform it into a format understandable by the client. This way, the client application only needs to know the front server location and does not directly contact the big data cluster.

We implemented the client application to run on modern browsers using HTML5 and WebGL for rendering. As a threshold for the KDE rendering, we use the variable MAX FRAGMENT UNIFORM V ECT ORS, which is the maximum number of variables that can be sent to the fragment shader.

BENCHMARKS

The visualization of density functions on top of geographical maps (Geographical heatmaps) is frequently done by datavisualizers. To run this benchmark, a complete prototype for that specific application has been implemented. It uses the technologies presented in the section 6, as well as Google Maps.

To enable the use of our system on low bandwidth network, the number of points per tile n 0 (cf. sec 4.1) has been set to 1000. This parameter enables to transfer only a few kilobytes (maximum 15KB in all our tests) on the network during interactions.

Infrastructure: The clustering benchmark has been run on our lab infrastructure composed of 16 computers (non virtualized). Each computer has 64GB of RAM, 2x6 hyperthreaded cores at 2.1GHz and 2 hard drives of 1 TB each. The computers are linked by a 1Gb/sec network infrastructure. During program execution, one computer serves as the master node and 15 as slave computing nodes. The rendering benchmark was conducted with an NVIDIA 970M graphics card for the GPU rendering and a Core i7-4710HQ for the CPU rendering. Datasets: We used four free datasets with sizes of different orders of magnitude, two small and two large. We give the datasets size in terms of records and in terms of unique positions, since the computational complexity of a density function depends on the number of different positions. The smallest one is the Brightkite dataset (4.7M records and 693K different positions). This dataset [START_REF] Cho | Friendship and mobility: user movement in location-based social networks[END_REF] has been built by recording the positions of each Brightkite user between April 2008 and October 2010. Heatmaps of this dataset have been made by Liu et al. [START_REF] Liu | imMens: Real-time visual querying of big data[END_REF] and Li et al. [START_REF] Li | Interactive visualization of high density streaming points with heat-map[END_REF]. A view of this dataset in our application can be seen on Figure 8. The Gowalla dataset [START_REF] Cho | Friendship and mobility: user movement in location-based social networks[END_REF] (6.4M records and 1.2M different positions) is similar to the Brightkite dataset but for Gowalla's users between february 2009 and october 2010. The OpenStreetMap (OSM) dataset2 (2.2B records and 1.7B different positions), this one represents every point of interest registered in the OpenStreetMap database. On Figure 8, one can grasp the scale of this dataset. Finally, the biggest dataset we used is a collection of GPS traces registered in the OpenStreetMap database (2.7B records and 2.2B different positions). The OSM and GPS datasets have very different point distributions, with OSM being more evenly spread than GPS. This made the OSM dataset more demanding in terms of memory needed for the clustering process, despite its lower number of unique positions.

Clustering benchmark: For each dataset, 21 levels of detail were generated. In HBase, the storage of these 21 levels necessitates 381MB for the brightkite dataset, 668MB for the gowalla dataset, 609GB for the OpenStreetMap dataset and 321GB for the GPS dataset. Unsurprisingly, there is an additional cost due to the storage of the entire level hierarchy. The difference in storage size between OSM and GPS is due to their different point distribution, leading to quicker aggregation. The figure 7 presents the results obtained on our lab infrastructure. We measured the scalability of our implementation using Spark executors. In Spark, executors are computation workers, enabling a finer granularity than computers. This enables to test finely the behaviour of the program 1468). On the GPS view, the high density in Europe is clearly highlighted. From all those view, it is clear that the density function is accurately represented using much fewer points. Furthermore, independently of the zoom level, the number of points displayed has the same magnitude. Figure 9: GPU rendering time on an NVIDIA 970M, up to 32000 points. In this case, N = 512. The sudden drops correspond to the changes of KDE resolution. One can see that the rendering is bounded by a constant. In practice, with n 0 = 1000, there is no more than 16000 points in a 1000*1000 visualization.

when the resources grow. Our benchmark demonstrates clearly that our approach is able to scale horizontally when treating really large datasets. For example, generating the levels for the GPS dataset took less than 45 minutes with 64 executors having 2 cores each.

Rendering benchmark: We now demonstrate that our technique is able to render KDEs with a large number of points in bounded time. Figure 9 shows the results of a benchmark with up to 32000 random points. The rendering time always stays below 11ms, which corresponds to more than 90fps. As expected, this technique maintains an acceptable rendering time by diminishing the KDE resolution when the number of points grow, as highlighted by the drops in Figure 9. Between the resolution changes, the rendering time scales linearly with the number of points.

Image Quality benchmark: To show the relevance and efficiency of the pipeline proposed here, we compare its results with a complete rendering of the entirety of the data. For this experiment, we chose the brightkite dataset because its small size enables local computations. We use a CPU rendering of this dataset as the ground truth and compare it to the levels of abstraction produced by our pipeline, rendered both on CPU and GPU. For this comparison, we use the SSIM image metric [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], which gives a score of similarity between two images. This metric is known to be closer to human perception than per-pixel error metrics. Figure 10 shows the resulting similarity scores. The quality is first limited by the clustering algorithm and then by the rendering algorithm. This creates a "bell" shape of optimal settings. From the angle of image quality, it is thus equivalent to render with more points on a smaller resolution or at full resolution with fewer points. This property demonstrates the validity of our approach, and opens the perspective to adapt the abstraction level to the device used.

CONCLUSION & FUTURE WORK

We presented a technique that enables interactive exploration of point sets at really large scale. We introduced a new distributed algorithm to compute a canopy clustering as well as a GPU based method to render density functions in constant time, supported by a theoretical justification. To validate our assumptions, a full implementation of the technique was made and a complete benchmark on large datasets has been run. All our results prove the efficiency of our approach in term of both horizontal scalability and quality of the produced visualization. Our technique is thus able to interactively explore sets of points of any size, the only limiting factor is the size of the Big Data Infrastructure.

In its current implementation, this system is not able to support incremental updates, however, the clustering algorithm used is fully suited for the task. Future works will be to extend that method for the visualization of dynamic set of points. New Big Data technologies, like Storm, must be investigated to produce efficient solutions.

Figure 1 :

 1 Figure 1: Point set visualization of the 1854 London cholera epidemic. (a) 512×512 plots of the original point set, (b) 512×512 visualization using a density function, (c) 256×256 rendering, due to occlusion it is difficult to detect dense regions, (d) 256×256 heatmap, even at low resolution, dense regions are easy to detect.

Figure 2 :

 2 Figure 2: On image (a), the density function visualization of the original point set. The black dots represent the data points. The density has been computed using Kernel Density Estimation with a Gaussian kernel. On image (b), the density function visualization of the corresponding point set for Approximate Kernel Density Estimation. The closest points have been merged into a single point whose weight is the sum of the weights of the merged points. For instance, points A, B and C have been merged into A'. The weight of A' is the sum of the weights of A, B and C, which is 3. The differences between the two generated pictures are almost invisible to the naked eye. Using the SSIM image comparison algorithm[START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF], we obtain a similarity score of 0.959, which corresponds to a "good" rating, according to Zinner et al.[START_REF] Zinner | Impact of frame rate and resolution on objective qoe metrics[END_REF]. Image (c) shows the result of another merging pass. Points B' through F' have been merged into A", whose weight is now 10. Differences in the density visualization become visible because the grouping distance is high in respect to the screen space and to the σ used for the KDE.

Figure 4 :

 4 Figure 4: Process of the first Map Reduce job. The map phase assigns a strip to each point, and the reduce phase executes the canopy clustering locally.

Figure 5 :

 5 Figure 5: An example of the data partition in the second job, focusing on the border of two strips. Image (a) shows the result of the first job.When the strips are merged, as shown on image (b), it creates a zone called "conflict zone", whose height is 2d. This zone is colored in grey. In this zone, two points coming from the two strips can be closer than d. Image (c) shows the partitioning used in the second job. In order to be able to process the entirety of the conflict zone in a single strip, the strip limit is shifted downwards by d. As illustrated on image (d), the conflict zone is then totally contained in the upper strip. The conflict can be detected and solved. One of the points will be removed from the canopy set during the second job.

  Map point: emit (strip(point), point); Reduce (strip, points): canopies = new QuadTree; for every point p in points do s = elements in canopies with dist(p) less than d; if s is empty then add p to canopies; for c in canopies do emit c; Algorithm 2: First job pseudo code for the Map and Reduce functions. The Map computes the strip for each point. The Reduce executes the canopy clustering on each strip.

Figure 6 :

 6 Figure 6: As the data grows, the resolution of the KDE diminishes to keep a constant rendering time.

Figure 7 :

 7 Figure 7: (a): Benchmark for the brightkite and gowalla datasets. Each executor has 4 cores. The optimal number of executors is attained around 10 because of the small size of those datasets. (b): Benchmark for the OSM and GPS datasets. Each executor has 2 cores. The benchmark for the OSM dataset was stopped at 30 executors due to memory constraints. This benchmark shows the scalability of our clustering algorithm.

Figure 8 :

 8 Figure 8: Views of the four datasets, and the number of points displayed. Top row, from left to right: brightkite (1671), gowalla (1334), and OSM (2873). Bottom row: GPS dataset, global view (3185), Europe view (3218) and Chicago view (1468). On the GPS view, the high density in Europe is clearly highlighted. From all those view, it is clear that the density function is accurately represented using much fewer points. Furthermore, independently of the zoom level, the number of points displayed has the same magnitude.

  The resulting group is called a canopy. Every canopy has a center point. A point can start a new canopy if it is further apart from the other centers than d. If a point is closer than d from a center, then it is in its canopy. Each canopy can be represented by its center with a weight equal to the sum of weights of points in the canopy. While the original algorithm uses two different distances, we use a simplified version where those distances are equal. This version of the algorithm is presented in algorithm 1.

	Data: A set of points, a threshold distance d
	Result: The set of canopies
	canopies = empty set;
	for every point p in the set do
	for every canopy c in canopies do
	if distance(p,c) ¡ d then
	go to next point;
	add p to canopies;
	return canopies
	Algorithm 1:

http://mahout.apache.org

http://spatialhadoop.cs.umn.edu/index.html

ACKNOWLEDGEMENTS

This work has been carried out as part of the "REQUEST" (PIAO18062-645401) and "SpeedData" (PIAO17298-398711) projects supported by the French "Investissement d'Avenir" Program (Big Data -Cloud Computing topic).

Figure 10: Comparison of the ground truth to the abstractions generated. The colored zones denote the user ratings proposed by Zinner et al. [35]: dark green for "Excellent", light green for "Good", yellow for "Fair" and orange for "Poor".


National INIS Centre: New Directions and Main Achievements in an Evolving Context

In a letter dated October 8, 1969 and addressed to the Director General of the IAEA, the French Governor of the IAEA Board of Governors gave his official consent for France to participate in the newborn INIS. At the same time, it was decided to establish the France National INIS Centre at the premises of the Central Documentation Department, located at the Saclay Nuclear Research Centre (French Atomic Energy Commission -CEA), near Paris. Forty-five years later, the Central Documentation has become the "Information Valorization Service" (SVI in French) which incorporates INIS, the CEA-Saclay Scientific Library and several other scientific and technical information-related activities (archives, scientific watch, collaborative work, bibliometric surveys, etc.).

INIS-France: workforce versus production

Until 2006, France was also a member of the Energy Technology Data Exchange (ETDE). During the pre-2000 era, the workforce dedicated to both INIS and ETDE input preparation was quite high (Fig. 1). The workflow process was organized so that all bibliographic inputs prepared for the INIS database were simultaneously submitted to ETDE's 'Energy' database. The 'pure' ETDE inputs (corresponding to non-INIS topics) represented about 30% of France's annual production of the period. In 2006, for both budgetary and workforce reduction reasons, France decided to withdraw from the ETDE and to focus its efforts on INIS. Additionally, France took this opportunity to cut down on subcontracting over the following years. However, the gradual decline of the workforce dedicated to this activity -which slowly began in the early 2000s -became more pronounced by 2008. Thus it became necessary for France to subcontract once again and to look for new means of production in order to maintain a good level of coverage and contributions to INIS, comparable to that of the preceding period.

Classical production versus voluntary inputting

A solution emerged during an informal discussion at the 11th INIS/ETDE Joint Technical Meeting (Vienna, November 2007). At that time, the INIS Secretariat was already making extensive use of its Computer-Assisted Indexing (CAI) system for the preparation of tens of thousands of journal article inputs (per year) to the database, using records purchased from big international editors. As indexing was no longer a problem for the Secretariat, their main concern was now access to specific documents, in particular the conference proceedings, and their bibliographic description. For France, the situation was the other way around: we had many conference proceedings on hand waiting to be processed, with the capability to prepare the bibliographic description of hundreds of articles, but not enough workforce or time to analyze and index each of them. France's annual production (Fig. 2). collection of CEA's technical notes (the 'CEA-N' collection). This second project will be completed before mid-2015 and all documents will be made available in full text through the INIS Collection. Last but not least, a third project concerns the collection, digitization and preservation of all historical and recent PhD theses prepared at the CEA. A 'CEA PhD Theses database' has been specifically created on the CEA-Saclay Intranet to make the full text of these PhD Theses accessible to all CEA researchers. Among all collected theses, about one third deal with nuclear topics and are consequently referenced in INIS in parallel. This project is bound to be developed since the number of PhD theses prepared at the CEA has dramatically increased in recent years with the number of thesis defenses exceeding 500 (Fig. 3). Such a project requires considerable time and effort, not only to process the documents, but also to retrieve the original manuscripts, since compulsory depositing of PhD Theses has never been implemented at the CEA. For this reason, even in recent years, manuscript collection is never 100% successful.

Grey literature digitization and knowledge preservation

PhD Theses and promotion of INIS

France participated in the Open Access Pilot Project launched in 2008 which led, the following year, to open access of the INIS database on the Internet. Since then, promoting INIS really started to make sense in order to inform the public at large about the richness of the database resources. Attracting new users is a task we have managed to carry out in two ways: internally, using flyers which, for instance, are made available to the users of the Scientific Library (Fig. 4), and externally, by targeting the university laboratories potentially interested in using INIS. These laboratories are identified thanks to the information recorded in the CEA PhD Theses database. In practice, for each recent CEA thesis referenced in INIS, we address a personalized E-mail to the PhD student and to both his university and CEA supervisors. The text of the message stresses the particular interest of the PhD work for the nuclear sciences and technologies, which is the reason for its referencing in INIS, and provides the direct link to the INIS record describing the document. In this way, for 145 CEA-theses published in 2012 and entered into INIS, more than 400 E-mails were prepared and sent to PhD students and laboratory researchers or managers, representing as many new potential INIS users.

Statistics of INIS usage in France

The indexing of the INIS Collection Search (ICS) by Google.com and Google Scholar search engines at the end of 2013 led to an impressive increase in the number of unique visitors worldwide. A similar result has been observed in France (Fig. 5) with the monthly number of unique visitors oscillating between 2500 and 6000 since January 2014 (monthly average of 4250 visitors compared with the previous years' monthly average of 150-200 visitors). Among those users identified in 2014, CEA is leading with 3726 unique users, followed by Areva (1605 users), the National Nuclear and Particle Physics Institutes (IN2P3 -1047 users), Electricite de France (EDF -1036 users) and followed by a huge list of universities and research institutes that highlight the pertinence and efficiency of our promotional actions.

Conclusion

This brief look back over the last 20 years of the France National INIS Centre allows us to measure how far we have come and what choices we have had to make in a constantly evolving context. We have particularly stressed what, in our view, represents the most important mission of INIS, which is the preservation and dissemination of nuclear knowledge, in particular through the coverage and full text archiving of conferences and grey literature documents, with special attention paid to PhD Theses.

FIG. 3 .

 3 FIG. 1. INIS France:1996-2014 staff distribution per category

FIG. 4 .

 4 FIG. 4. Two INIS promotional flyers made available for library users.


Abdelhakim Latoui
email: latoui.ahakim@yahoo.fr

Raphael Canals
email: raphael.canals@univ-orleans.fr

Remote Online Testing Of Embedded Systems Using Optical BILBO

Keywords: Embedded Systems, On-Line Testing, Test and Testability, VHDL-AMS, Optical Sensors

In this paper, an online testing approach that checks continuously at remote system the correctness of embedded systems in full operation is proposed. Our new method exploits optical beams produced by an Optical Built-In Logic-Block Observation (OBILBO) register based on the captured data presented on the outputs of all stages of the registers in the BILBO. Then it is possible to send the captured data to remote system equipped with optical sensors in real time way. The final captured response data can also be used for final comparison with stored data of a golden circuit. Preliminary simulation results showed that faults are concurrently detected without affecting normal system operation and without having recourse to any external or internal Automatic Test Pattern Generation (ATPG).

I. INTRODUCTION

Nowadays, embedded systems are omnipresent in many applications whose contributions are significant in terms of weight, size, power consumption, response time and flexibility. Although their sophistication continues to increase, embedded systems are increasingly used in critical applications (e.g. medical or transport applications) in which their misbehavior may generate unacceptable consequences that may go up to human life loss [START_REF] Ma | An Adaptive Risk Control and Security Management for Embedded Real-Time System[END_REF]. It is therefore vital to ensure immediate detection of errors caused by a failure in these systems by having recourse to online testing. These tests can check continuously the results produced by these systems in exploitation mode to confirm the validity of their functioning [START_REF] Ferranti | CAST: Automating Software Tests for Embedded Systems[END_REF]. However, this kind of online test constitutes nowadays a big challenge for test designers of such systems.

So, a large variety of on-line testing techniques has been proposed in the literature [START_REF] Nicolaids | On-Line Testing for VLSI-A Compendium of Approaches[END_REF][START_REF] Papa | Optimal On-Line Built-In Self-Test Structure for System-Reliability Improvement[END_REF]. Such techniques are for instance: self-checking design, signature monitoring, on-line monitoring of reliability relevant parameters such as current or temperature, implementation of BIST techniques specific to on line testing, etc.

However, in modern System-on-Chip (SoC) design, many cores are integrated into a single chip. Some of them are embedded and cannot be accessed directly from the outside of the chip. Such SoC designs make the test of these embedded cores a great challenge. Although the BIST technique is one of the most popular solution [START_REF] Abramovici | Digital Systems Testing and Testable Design[END_REF][START_REF] Mourad | Principles of Testing Electronic Systems[END_REF] adopted for the test of embedded cores, it is used to realize a discontinuous online test by activating periodically the necessary test sequences. Detection of the faults is then carried out with a delay which is related to the periodicity of the test.

Furthermore, the International Technology Roadmap for Semiconductors (ITRS) [START_REF]The International Technology Roadmap for Semiconductors[END_REF] projects that a Giga-Scale Integration (GSI) chip manufactured at the 18nm technology node will have a chip-to-board speed of 56GHz. To this end, it identifies optical interconnects as a potential enabling technology. Recent reports relating to the monolithic fabrication of photo detectors, modulators and emitters on Si have given rise to the exciting prospect of combining optoelectronics and GSI in a high-volume manufacturing environment [START_REF] Thacker | Wafer-Testing of Optoelectonic-Gigascale CMOS Integrated Circuits[END_REF]. In fact, the domain which involves measuring the electrical quantity such as voltage or current, using optical methods is adopted to perform testing modern integrated circuits [START_REF] Tsang | Time-resolved optical characterization of electrical activity in integrated circuits[END_REF]. The technology that uses electro-optic (EO) probing of electrical circuits with laser beam is developed by several research teams [START_REF] Yang | Electric field mapping system using an optical-fiber-based electrooptic probe[END_REF]. These techniques have for object the characterization of off-line circuits. EO probing can also be used, as detailed in [START_REF] Pannetier | External Optical Probe for Novel Online Testing Method of Digital Circuits[END_REF][START_REF] Aktouf | On-line Testing of embedded Systems Using Optical Probes: System Modeling and Probing technology[END_REF], to make real time measurement of the state (i.e. high or low) of digital signals present at the outputs or the inputs of a circuit board or integrated chip in order to detect any failure of the latter. So, the combination of BIST approach and EO probing can lead to viable solution for on-line testing of nowadays circuits. For example, in [START_REF] Latoui | An Optical BILBO for On-Line Testing of Embedded Systems[END_REF], an Optical Built-In Logic-Block Observation (OBILBO) was proposed to achieve on line detection of stuck at fault of embedded systems with very low error latency. Nevertheless, the solution has an area overhead due to the use of the hardware duplication.

This paper presents a more efficient solution in terms of area overhead. The novel proposed approach is based on the OBILBO which permits the emission at remote system of optical beams according to the logical states of its outputs. When operating in parallel load mode, the OBILBO acts as snapshot of outputs (inputs) data of the Circuit Under Test (CUT). Thus, the transmitted data, corresponding to the outputs (inputs) of the CUT, via the OBILBO, can be compared at remote system in a real time way to the stored data of a Reference Model (RM) of the CUT. The stored data replace and thus eliminate the complete duplication of the CUT.

The rest of the paper is organized as follows. Section 2 describes the new remote on-line testing scheme. The corresponding VHDL-AMS model is presented in Section 3. Analysis and experimental results are presented in Section 4. Section 5 concludes this paper.

II. REMOTE ONLINE TESTING SCHEME

An ideal on-line testing scheme would have 100% error coverage, error latency of 1 clock cycle, no space redundancy, and no time redundancy. Considering all four parameters in the design of an on-line testing scheme may create conflicting goals. High coverage, for example, requires high error latency, space redundancy and/or time redundancy. On the other hand, schemes with immediate detection (error latency equaling 1) minimize time redundancy but require more hardware [START_REF] Hayes | J "Online BIST for Embedded systems[END_REF]. However, reducing the error latency is often considered as a primary goal in on-line testing [START_REF] Al-Asaad | Efficient Techniques for Reducing Error Latency in Online Periodic Built-in Self-Test[END_REF]. As a matter of fact, techniques using COTS (Component-Off-The-Shelfsuch) as duplication or triplication are more often adopted even at the extra area overhead, since they have low error detection latency compared to others techniques described earlier [START_REF] Nicolaids | On-Line Testing for VLSI-A Compendium of Approaches[END_REF]. This paper presents a new on line testing scheme that allows errors detection with very low error latency and minimum area overhead.

The overall scheme of the proposed approach is illustrated in Figure 1 where the CUT can be any chip supporting Built-In Logic-Block Observation (BILBO) architecture. The BILBO is an embedded and off-line BIST architecture because it uses existing flip-flops from the CUT to construct the TPG (Test Pattern Generator) and ORA (Output Response Analyzer) functions [START_REF] Koenemann | Built-In Test for Complex Digital Integrated Circuits[END_REF]. Optical probes are then added to the BILBO register to form an OBILBO [START_REF] Latoui | An Optical BILBO for On-Line Testing of Embedded Systems[END_REF]. In addition to the standards functions [START_REF] Koenemann | Built-In Test for Complex Digital Integrated Circuits[END_REF] of the BILBO, the OBILBO continuously produces optical beams according to the logic states of its outputs thanks to optical probes integrated in its architecture. Figure 2 [START_REF] Latoui | An Optical BILBO for On-Line Testing of Embedded Systems[END_REF] shows the logic diagram of four stages OBILBO register.

To ensure real time checking of the CUT integrity without affecting its normal operation, the OBILBO1 and OBILBO2 (Figure 1) should be however configured in their parallel load mode. Thus, OBILBO1 and OBILBO2 produce optical signals corresponding to the logical states of the CUT inputs and outputs respectively. The transmitted optical signals can then be captured at remote system equipped with photo-detectors, such as PIN diodes, and Analog to Digital Converters (ADC). The final captured data are compared to the stored data of a RM, i.e. for each input data vector (corresponding to signals transmitted by OBILBO1), we establish a corresponding reference response data vector which is compared to the response data (corresponding to the signals transmitted by OBILBO2) of the CUT. It is worth noting that the constituting elements of the proposed scheme exchange signals which are of nature heterogeneous. In the next section, we give a possible modeling of the described scheme with the VHDL-AMS language.

III. VHDL-AMS MODELING

The VHDL-AMS language is an extension of the IEEE 1076 (VHDL) standard for describing digital, analog and mixedsignal systems. It provides mixed-disciplines modeling, wherein different domains can be described and simulated in a single entity. In fact, the system of the Figure 1 can be viewed as simple components connections which can be modeled using components instantiation.

Signal-processing chips are basic elements in today Embedded and Real-Time systems. A typical example of these is the TMS32010 signal-processing chip which supports BILBO architecture. Indeed, various BIST architectures were applied to this chip in order to achieve offline testing. Four BIST designs are proposed in [START_REF] Abramovici | Digital Systems Testing and Testable Design[END_REF]. Among these architectures a BILBO solution was proposed as shown in Figure 3 in which some registers are replaced by BILBO registers to test the multiplier and the ALU. In our case, all MISR and PRPG are considered as BILBO registers which are all replaced now by OBILBO registers. However, in order to test this chip in on-line mode, the CUT of Figure 1 is considered now as portion of the TMS32010 signalprocessing chip with its novel OBILBO registers, as depicted in Figure 4.

We assume now that the ALU of TMS32010 is considered as single adder cell as illustrated in Figure 4 (since more signals are required to represent the simulation results of the ALU). We suppose that the adder is faulty due the stuck-at-0 at the input of the grey OR-gate. Table 1 shows the results which should be produces in the faulty free case and in the faulty one. The corresponding VHDL-AMS model is given in Figure 5.

0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 1 0 1 1 0 0 0 1 1 0 0 1 0 1 0 0 1 1 0 0 1 0 1 0 0 1 0 1 1 0 0 0 1 1 1 1 1 1 1 1 1

IV. ANALYSIS AND EXPERIMENTAL RESULTS

Simulation results are reported in Figure 6 where b1 [START_REF] Ferranti | CAST: Automating Software Tests for Embedded Systems[END_REF] and b1 [START_REF] Ma | An Adaptive Risk Control and Security Management for Embedded Real-Time System[END_REF] are the control lines of OBILBO1 whereas a, b and cin are its inputs data. q1 [START_REF] Ferranti | CAST: Automating Software Tests for Embedded Systems[END_REF], . . . q1[0] are the outputs of OBILBO1 where their corresponding optical signals are ql21, . . . ql01 respectively which are also the inputs of the first set of the PIN diodes.

The ADCs placed after this first set deliver the digital signals d21, d11and d01 which correspond in fact to the inputs data a,b and cin respectively. However, d02 and d12 are the signals delivered by the second set of ADC placed after the second set of the PIN diodes which correspond to the sum and the carry produced by the faulty adder respectively. The optical signals for the second set are not shown in the figure. The signals s1 and c1 are the sum and the carry references delivered according to the inputs data d01, d11 and d12 which are in fact the converted data of the optical beams produced by OBILBO1 (corresponding to the inputs data a, b and cin).

The checker has delivered the signals er0 and er1 where er0 is the result of comparison between the s1 (reference sum) and d02 (corresponding to the faulty sum) and er1 is the result of comparison between c1 and d12 (corresponding to the faulty carry). We note that er0 is equal to 0 when s1 and d02 are identical and equal to 1 otherwise (Idem for er1). We can see from these results that the considered stuck-at-0 is detected after only one clock cycle (delay introduced by the OBILBO) with no having recourse to any ATPG. However, for a duplication technique errors are not detected until they propagate to the CUT outputs. For the sake of argument, considering the system of Figure 3, if we assume that the inputs of R1 and R2 are the inputs of the CUT and the outputs of R4 are its outputs, a fault at the input of the multiplier cannot be detected until it propagates to the outputs of R4. Therefore, the ALU and the R4 register will introduce a delay like for the multiplier and the R3 register. Thus, more clock cycles are required to detect this fault. As a result, the error latency will increase and it will increase more and more with the complexity of the CUT. It follows that it is clear the proposed solution is more efficient since it reduces considerably the error detection latency with minimum area overhead and also permits a comparison at remote system, contrary to the duplication technique.

V. CONCLUSION In this work, a new remote online testing scheme for embedded systems based on the Optical BILBO (OBILBO) is proposed. Simulation results show that the proposed solution not only ensures errors detection concurrently to the normal function of the CUT with very low error latency but also requires less area overhead, compared to the duplication technique. Then, the optical technology is promising for online testing of embedded systems but it may be nevertheless accompanied by other problems due mainly to interferences.

Figure 1 .

 1 Figure 1. Schematic representation of the proposed method.

Figure 2 .

 2 Figure 2. Logic diagram of OBILBO.

Figure 3 .

 3 Figure 3. Portion of the TMS32010 Signal-processing chip.

Figure 4 .

 4 Figure 4. Full Adder with OBILBO.

Figure 5 .

 5 Figure 5. VHDL-AMS Model of the proposed scheme.

Table I .

 I TRUTH TABLE

	A	B	CIN	SUM (Free)	SUM (Faulty)	Carry (Free)	Carry (Faulty)


Martin Monloubou

Arnaud Saint-Jalmes

Benjamin Dollet

I Cantat

Influence of bubble

Keywords: 47.57.Bc -Foams, Complex fluids 43.25.Ed -Acoustical absorption, Nonlinear acoustics

teaching and research institutions in France or abroad, or from public or private research centers.

Liquid foams have a surprisingly high effective viscosity, which can be several orders of magnitude larger than the foaming solution viscosity [START_REF] Cantat | Foams. Structure and Dynamics[END_REF]. Energy dissipation in foams may have very different local origins and their identification, as a function of the characteristic time scales and of the foam properties, is a long lasting fundamental and applied challenge [START_REF] Buzza | [END_REF]. For simple shear, dissipation has been shown to be directly coupled to the foam coarsening at very low frequency (10 -2 Hz) [3,4] and to arise from viscous dissipation between bubbles for frequencies around the Hz [5]. At higher frequencies, the foam loss modulus varies as the square root of the frequency and the local origin of this dissipation is still a matter of debate [6][7][8][9]. For isolated films vibrating at 1 kHz, the viscous dissipation in air is dominant [10].

Compressive waves in foams have been studied in the context of acoustic and blast waves. The need to optimize the blast wave absorption, in terms of liquid fraction, bubble size or chemical formulation, stimulated experimental studies at the film scale [11] and at the foam scale, both on field, using real detonations [12,13], and in the lab, using shock tubes [14,15]. Most theoretical approaches on blast wave absorption disregard the foam structure, and focus only on the role of the liquid fraction [16]. In the acoustic pressure range however, several studies clearly established the important role of the bubble size on the sound propagation, both theoretically and experimentally [17][18][19][20][21][22].

Compression waves modify the gas temperature and thus induce thermal effects. Thermal dissipation is exacerbated in the foam, as its structure makes especially efficient heat transfer between the liquid and gas phases. It has been computed for a bubbly liquid, in the limit of small gas fraction, evidencing two limiting cases, for bubbles of radius R much smaller or much larger than the thermal length ℓ T = √ χτ , with τ the excitation characteristic time and χ the heat diffusivity. The attenuation coefficient is predicted to scale as R 2 /τ 2 for small bubbles (isothermal regime), and as 1/(R √ τ ) for large ones (isentropic regime) [23,24]. A good agreement has been found with the first scaling for acoustic waves of frequencies in the range 40 -80 kHz propagating in commercial shaving foams [18].

This paper focuses on the attenuation of a compressive wave propagating in a liquid foam, with an overpressure duration τ of the order of one millisecond. We work on the moderate overpressure regime, below the critical pressure at which the foam is fully destroyed, and above the linear regime. Recording the pressure evolution in a foam sample, we determine the wave amplitude attenuation. We establish the scaling laws for the thermal attenuation in the context of dry foams and recover the two scaling laws proposed for bubbly liquids [24]. Experimentally, the attenuation varies as R 2 for small bubbles as in [18], but we evidence that it departs from this law for larger bubbles. We thus identify the cross-over between the isothermal and isentropic regimes, for R ≈ 5ℓ T , that coincides with an optimal energy absorption, distinct from the optimal energy absorption induced by the resonance [19,21,25].

The foam is produced by turbulent mixing of nitrogen and foaming solution at a flow rate of the order of 100 p-1 mL/s [26]. The solution is composed of SDS (sodium dodecyl sulfate) at a concentration 10 g/L in pure water. The foam is then injected in a cell made of two transparent square plates, through a hole in the middle of the top plate. The plates are 30 cm wide and the gap between them is H = 3 cm (see Fig. 1). The liquid fraction of the foam is controlled by the nitrogen flux and is measured a posteriori. Its value is φ l = (4.8 ± 0.3) % for all the results shown in the paper. As the foam evolves by coarsening, the average bubble size increases with time and bigger bubbles are obtained after some delay. During this delay, the cell is flipped upside down every 30 seconds to ensure a homogeneous liquid fraction in the cell. With this method, the accessible bubble radius is 100 -300 µm. Larger bubble radii, up to 500 µm, have been obtained by blowing nitrogen through four porous frits immersed in the foaming solution, below the cell, held vertically.

Just after the shock, a few hundreds of bubbles are extracted from the cell and spread on a thin layer of foaming solution, forming a monolayer of spherical bubbles. Their average radius R, determined by image processing, is used to quantify the bubble size in the foam [25]. The diameter distribution is roughly log-normal, with a normalized standard deviation of the order of 0.4.

A pressure wave is generated with a home-made shock tube inspired by [27], composed of a high pressure chamber separated from a low pressure one by a diaphragm. Nitrogen is injected into the high pressure chamber until the rupture of the diaphragm, which generates a shock wave (i. e. an overpressure step) that propagates in the tube. When reaching the exit of the tube, which is connected to the middle of the cell (see Fig. 1), the shock wave turns into a blast wave that propagates in the foam. The wave amplitude is varied by changing the diaphragm. We used aluminium foil (1 to 4 layers) and, for the weakest shocks, one parafilm layer. Four pressure sensors PCB 113 B28 are installed on the top plate at a distance r i from the centre, with i = [1 -4] and r i = [3.10, 5.15, 7.20, 9.25] cm. They are connected to an oscilloscope TDS 2004 B, and record the pressure p i (t) with a sampling time of 4µs. Fig. 2 shows typical pressure signals as a function of time, at different distances from the cell center. At sensor 1, a rapid pressure increase is followed by a slower pressure decrease and a small underpressure. This kind of profile is typical of a blast wave [12]. The duration of the overpressure is of the order of 2 ms and finishes before the first wave reflections come back from the external foam boundary, in contact with air. The results do not depend on the cell thickness H for H < 5 cm and the wave propagation is cylindrical.

The transmission coefficient from the shock tube to the foam depends on the bubble size and on the shock amplitude, but the rescaled pressure p 1 (t)/P 1 , with P 1 the maximal pressure at sensor 1, is independent on those parameters (see Fig. 2, inset). The initial signal p 1 (t) is thus fully characterized by the single parameter P 1 which is our control parameter for the shock amplitude in the whole study. The pressure ratio, α = P 4 /P 1 , with P 4 the maximum pressure value measured at sensor 4, is shown in Fig. 3 as a function of the amplitude, for 6 bubble sizes. α increases with the amplitude and can be fitted by a phenomenological affine law α f it (R, P 1 ) = α 0 (R) + KP 1 , with K = 6.5 Pa -1 the slope leading to the best agreement. The standard deviation is σ(R) = (α exp -α f it ) 2 1/2 . In the linear regime, the wave attenuation is classically quantified by an attenuation length ℓ a defined by the relation α = α ref exp[(r 1 -r 4 )/ℓ a ], α ref being the expected attenuation, in the non-dissipative case. In our cylindrical geometry, α ref = (r 1 /r 4 ) 1/2 = 0.58. The length ℓ exp a was computed for each bubble size using the value of α extrapp-2 The last two points have been obtained with the vertical cell. The liquid fraction is φ l = (4.8 ± 0.3) %. The solid lines are numerical results obtained for ν = 0.001, 0.17, 0.25 and 1.5 m 2 /s (from top to bottom) and γ = 1. The dashed line are phenomenological fits α f it = α 0 (R) + 6.5 P1 (in Pa). Inset: α 0 (R) ± σ(R) as a function of the bubble size. olated at small amplitude, α 0 (R) ± σ(R), shown in Fig. 3 (inset). The main result is that the pressure attenuation increases with bubble size for bubbles smaller than 300 µm (and the attenuation length thus decreases) and saturates for larger bubbles, at constant liquid fraction. The experimental pressure ratio, in the small amplitude limit, remains below the geometrical limit α ref as expected. However, at larger amplitudes the non-linearities reduce the attenuation, and α may cross this limit.

When the wave reaches a bubble, its pressure increases by typically δp ≈ 10 4 Pa during τ ≈ 1 ms. Considering an ideal gas behavior in the bubbles, this may lead to a temperature increase of δT s = δP/(ρ g c p ) ≈ 10 K if the process is isentropic, with ρ g = 1.2 kg/m 3 the gas density and c p = 10 3 J • kg -1 • K -1 its specific heat at constant pressure. The bubble is surrounded by the liquid phase, that plays the role of a thermostat at temperature T . Indeed, in a given volume of foam, the ratio of the thermal capacities of gas and liquid phases equals ρ g c p /(ρ s c p,s φ l ) = 5 • 10 -3 ≪ 1, with ρ s = 10 3 kg/m -3 the density of the solution and c p,s = 4 • 10 3 J•kg -1 •K -1 its specific heat. At a time τ after the bubble compression, the thickness of the thermal boundary layer is ℓ T = √ χτ ≈ 10 -4 m, with χ = 20 • 10 -6 m 2 /s the heat diffusivity. This length ℓ T is much smaller than the wave length λ = cτ ≈ 5 • 10 -2 m, with c ≈ 50 m/s the sound velocity (see inset of Fig. 4). It is, in contrast, just in the range of the investigated bubble size.

Thermal properties first influence the foam effective compressibility, and thus the wave velocity. A classical prediction for the linear sound velocity c is the Wood velocity. In this model, the foam is replaced by a continuous medium of density ρ = φ l ρ s related to the local pressure p tot by a constitutive relation [START_REF] Wood | A Textbook of Sound[END_REF][START_REF] Surov | [END_REF]

p tot = p 0 (1 -φ l,0 )ρ ρ 0 -φ l,0 ρ γ , (1)
with p 0 the initial atmospheric pressure, ρ 0 and φ 0 the foam density and its liquid fraction at p tot = p 0 . The gas can be modeled as a ideal diatomic gas and γ = 1.4 or γ = 1 for the isentropic and isothermal cases, respectively. The Wood velocity is

c 2 W = ∂p tot ∂ρ ρ=ρ0 = γp 0 (1 -φ l,0 )φ l,0 ρ s . (2)
For φ l = 0.048, c t W = 47 m/s for γ = 1, and c s W = 56 m/s for γ = 1.4. The small amplitude front velocity, measured optically for a better precision (see Appendix A), is better fitted with the isothermal Wood velocity for small bubbles and by the isentropic velocity for large bubbles (see inset of Fig. 4). Note that the front velocity we measured may be sligthly larger than the average velocity of the whole signal, because of the signal spreading. It thus should be considered as an upper boundary for the sound velocity.

More importantly, thermal gradients also lead to energy dissipation and significantly contribute to the wave attenuation. The isentropic and isothermal assumptions are valid if the bubble size is respectively much larger and much smaller than ℓ T [24]. We thus derive the scaling laws of the thermal dissipation in the linear regime for small bubbles (isothermal model) and for large bubbles (isentropic model). We then extrapolate these two limiting cases to cover our whole bubble size range.

The mechanical wave carries the energy per bubble E 0 ∼ ρ 0 u 2 R 3 ∼ ρ 0 [δp/(ρ 0 c)] 2 R 3 . This energy decreases with time as [START_REF] Landau | Physique théorique : mécanique des fluides[END_REF]

Ė ∼ κ T ∇T 2 Ω dis with κ = 2 • 10 -2 W•m -1 •K -1 the
thermal conductivity of nitrogen and Ω dis the volume of dissipation, in one bubble.

The attenuation length is larger than the wave extension, so we can use the properties of the inviscid wave propagation to determine the scaling law of the energy dissipation [START_REF] Landau | Physique théorique : mécanique des fluides[END_REF].

We first consider the isentropic limit. In that case, the temperature in a bubble during the overpressure is T +δT s and varies close to the liquid phase at T over the typical distance ℓ T = √ χτ < R, so that Ω dis ∼ R 2 ℓ T and ∇T ∼ δT s /ℓ T . Using χ = κ/(ρ g c p ) we get:

1 ℓ s a ∼ Ė cE 0 ∼ cρ 0 ρ g c p T R κ ρ g c p τ , (3)
giving an attenuation length ℓ s a ≈ 10 cm for τ ≈ 10 -3 s and R ≈ 10 -4 m.

The same kind of analysis can be performed for the isothermal wave propagation. In that case, the gas remains at constant temperature and thus at constant energy (ideal gas assumption). The mechanical work on a p-3 bubble δW ∼ pδV ∼ R 3 δp thus equals the heat exchange with the liquid phase and δQ ∼ R 3 δp for each bubble, during the characteristic time τ . This heat flux is ensured by thermal conduction and thus δQ ∼ κ∇T t R 2 τ ∼ κδT t Rτ . From these two relations, we deduce the small temperature variation δT t between the bubble center and the bubble boundary. The temperature gradients scale as δT t /R and Ω dis = R 3 . We finally get:

1 ℓ t a ∼ Ė cE 0 ∼ R 2 cρ 0 κT τ 2 , (4)
giving ℓ t a ≈ 10 cm, like ℓ s a , for the same parameters. The dissipation in the continuous phase may also be relevant in the problem. When the foam is compressed, the liquid phase remains incompressible, but the films and the meniscus are sheared or elongated in a complex way which depends on the local foam structure. This induces viscous dissipation in the liquid phase, as well as a specific interfacial dissipation related to the interfacial shear and dilatational viscosities. When the foam is simply sheared, the gas is not compressed and the thermal dissipation is thus probably irrelevant. In contrast, the film and meniscus deformations are a priori of similar nature as under foam compression. Even if the local processes associated to the global deformation of the foam are mostly unknown, the dissipation in the continuous phase can thus be qualitatively estimated from the effective viscosity of the foam measured in simple shear experiments. Kinematic shear viscosities at 100 Hz of the order of ν v = 0.02 m 2 /s have been reported [30,31]. The viscous attenuation length ℓ v a scales as 1/ℓ v a = ν v /(τ 2 c 3 ) [START_REF] Landau | Physique théorique : mécanique des fluides[END_REF], which would lead to ℓ v a ∼ 5 m. As the foam is known to be a shear thinning material, its effective viscosity should be even smaller at higher frequencies. Viscous dissipation is thus a priori negligible in comparison with thermal dissipation.

Whatever their physical origin, all dissipative processes play the mechanical role of an effective viscosity, with ν ∼ τ 2 c 3 /ℓ a [START_REF] Landau | Physique théorique : mécanique des fluides[END_REF]. From the order of magnitude determined above, ν is governed by the attenuation (3) for small bubbles and (4) for large ones. The cross-over is reached for R ∼ ℓ T , which corresponds to a maximum of viscosity. This predicts that, in a compressive regime, the effective kinematic viscocity of a foam increases as R 2 for R < ℓ T , reaches a maximum and then decreases as 1/R, as already proposed in [18].

We define the experimental viscosity by ν exp = bτ 2 c 3 /ℓ exp a , with b a numerical prefactor close to one. The parameter b = 1.2 leads to the best agreement with the numerical simulation discussed below. Using this value, we obtained the effective kinematic viscosity ν exp for all bubble sizes, as shown in Fig. 4. It lies in the range [0.2 -2] m 2 /s and increases with the bubble radius up to 300 µm. These experimental data are fitted by the expression ν t = 0.4 τ 2 c 3 /ℓ t a , which shows a good agreement for the bubble radii in the range [100-300] µm. In contrast, the last point at 500 µm is incompatible with the quadratic increase of the viscosity with the bubble size. This radius is much larger than the thermal layer thickness ℓ T = 100 µm, hence the adiabatic limit is reached. The effective viscosity ν s = 7 τ 2 c 3 /ℓ s a is plotted in Fig. 4 and reconciles the theoretical and experimental results. Unfortunately, we have not been able to produce foams with larger bubbles and still a 5% liquid fraction to evidence the viscosity decrease at large bubbles. We have performed experiments at larger R, albeit at smaler φ l , that suggest this decrease of attenuation.

Beyond the linear regime, the foam properties must be compared to numerical simulations. A minimal model consists in considering the foam as a compressible and dissipative continuum of density ρ 0 , undergoing weakly nonlinear sound propagation. Hence, we performed numerical simulations based on the second-order equation of sound propagation,

The wave propagation equation is obtained from the series expansion of the compressible, axisymmetric, Navier-Stokes equation at the second order in the overpressure p(r, t) = p tot -p 0 and in the radial foam velocity u(r, t) [START_REF] Morse | Theoretical Acoustics[END_REF] (eq. 14.4.12):

∂ rr + 1 r ∂ r p + ν c 2 ∂ t p - 1 c 2 ∂ tt p = -ρ 0 ∂ rr + 2 r ∂ r (u 2 ) + Γ ρ 0 c 4 ∂ tt p 2 , (5)
where ∂ X and ∂ XX respectively stand for partial derivatives of first and second order with respect to variable X. The parameter Γ is ρ 0 c 4 (∂ 2 ρ/∂p 2 ) 0 /2. For an ideal gas, Γ = (1 -γ)/2 < 0. The parameter ν is a dissipative parameter with the dimension of a kinematic viscosity. We only kept first order terms in ν, the validity domain thus being limited by the conditions p ≪ p tot and ν ≪ c 2 τ with τ ≈ 1 ms the typical duration of pressure variation. This last condition yields ν ≪ 3 m 2 •s -1 , with c ≈ 50 m/s. Practically, reliable results have been obtained for ν < 2 m 2 /s. The mass conservation closes the problem and only needs to be computed at order one:

1 c 2 ∂ t p + ρ 0 ∂ r + 1 r u = 0 (6)
We used the phenomenological fit of the pressure p 1 (t) at sensor 1 shown in Fig. 2 as a boundary condition at the inner boundary. The foam is assumed to be at rest at large r. This set of equations is solved numerically with an implicit scheme implemented in a Matlab code, which has been validated with a precision of 1% using known nonlinear analytical solutions (see Appendix B).

Qualitatively, the viscous term in (5) induces a signal attenuation and a signal spreading, thus leading to an enhanced pressure attenuation. In contrast, the non-linear term acts against the signal spreading, and thus slows down the maximal value decrease. These features, observed experimentally, are reproduced by the model (see p-4 Fig. 2). However, for large bubbles and large amplitudes, the model fails to predict the high attenuation experimentally observed. Moreover, the front pressure gradient is always underestimated. This could come from the fact that the numerical simulation disregards any acoustic dispersion [21].

The numerical value of the attenuation α num is plotted in Fig. 3. For each bubble size the dissipative coefficient has been chosen to get the best agreement with α f it . Its value ν num (R) is plotted in Fig. 4 as a function of the bubble size. The error bars have been obtained by fitting numerically α f it ± σ. (o) viscosities deduced directly from the experimental attenuation, as defined in the text, (x) viscosities obtained as fitting parameter in the numerical simulations, when possible. The error bars are deduced from the experimental dispersion σ(R). The solid and dashed lines are the predictions of eq. ( 4) and eq. ( 3), rescaled vertically by a factor 0.4 and 7, respectivelly. Inset: Sound velocity as a function of the bubble size, in the small amplitude limit. Rectangles represent the isothermal (bottom) and isentropic (top) Wood velocities, with the error bar arising from the liquid fraction error bar.

In conclusion, using a foam with a well controlled chemical composition, bubble size and liquid fraction, we show that the main dissipative processes that mitigates a compressive wave, in a time scale domain of the order of τ = 1 ms, are of thermal origin. The associated effective viscosity increases as R 2 /τ 2 at small bubbles and departs from this law for larger bubbles. Our theoretical analysis is based on the Wood model. This has been validated by checking that the propagation velocity does not exceed the Wood velocity. It has recently been shown that above a certain radius, Wood's model breaks down and a resonant behavior is expected. For our characteristic time of 1 ms, the resonant bubble radius is of the order of 1 mm [21]. However strong variations of the sound velocity are expected at the resonance, with a sound velocity of order 200 m/s for R > R r , which is not compatible with our y r observations. We thus believe that the wave propagation we observe is below the resonance for our whole parameter range (i.e. at a radius smaller than R r ). Moreover, the local process responsible for the attenuation during the resonance was not identified, and the thermal dissipation may be the relevant one. These results seem limited to compressive waves. However, it is known that plastic events triggered in foams by simple shear induce pressure variations, and thus temperature variations. In this context, the associated thermal dissipation may thus also be important in some parameter range. . Numerically, setting the inner boundary at the position r 0 = 3000 δ to be in a plane wave geometry and using the eq. ( 7) as pressure condition at this point, we got an agreement between numerical and analytical solution better than 1 % on the velocity and on the thickness of the pressure step, for the same time and space discretisation and the same global time and space scales as used in the next section. Second, we tested the energy conservation for an inviscid case at very small amplitude. Using the low amplitude experimental data as boundary condition at the inner boundary (at its experimental position) we found that ∞ r=0 rp 2 dr, which is proportionnal to the wave energy, varies by less than 1% between t max 1 and t max 4 , for a vanishing viscosity.

FoamFig. 1 :

 1 Fig. 1: Sketch of the experimental set-up.

1 Fig. 2 :

 12 Fig. 2: Pressure at sensors 1 to 4 (solid lines, from left to right) as a function of time for R = 210 µm. Dash-dotted line: signal used as numerical boundary condition for the simulations. Dashed lines: numerical simulations of the pressure at sensor 2 to 4, with ν = 1.1 m 2 /s and γ = 1. Inset: rescaled signals p1(t)/P1 at sensor 1 for R = 115 (black), 210 (blue) and 265 (red) µm for 4 aluminium foils (solid line), 2 foils (dash-dotted) and parafilm (dashed) (color online).

Fig. 3 :

 3 Fig. 3: Pressure ratio in the foam α = P4/P1 as a function of the initial wave amplitude P1, for different bubble sizes R: 117 ( ), 182 (⋆), 210 ( ), 265 ( ), 303 (•), and 535 ( )µm.The last two points have been obtained with the vertical cell. The liquid fraction is φ l = (4.8 ± 0.3) %. The solid lines are numerical results obtained for ν = 0.001, 0.17, 0.25 and 1.5 m 2 /s (from top to bottom) and γ = 1. The dashed line are phenomenological fits α f it = α 0 (R) + 6.5 P1 (in Pa). Inset: α 0 (R) ± σ(R) as a function of the bubble size.

Fig. 4 :

 4 Fig.4: Foam effective viscosity as a function of the bubble size; (o) viscosities deduced directly from the experimental attenuation, as defined in the text, (x) viscosities obtained as fitting parameter in the numerical simulations, when possible. The error bars are deduced from the experimental dispersion σ(R). The solid and dashed lines are the predictions of eq. (4) and eq. (3), rescaled vertically by a factor 0.4 and 7, respectivelly. Inset: Sound velocity as a function of the bubble size, in the small amplitude limit. Rectangles represent the isothermal (bottom) and isentropic (top) Wood velocities, with the error bar arising from the liquid fraction error bar.

Fig. 5 :

 5 Fig. 5: Image of the foam few ms after the shock. The rectangle is the reduced field of view used for the measurements.

Fig. 6 :

 6 Fig. 6: Space-time image of the foam. The position of the wave front deduced from the image is plotted in black solid circles.

= c 1 +

 1 lytical cases have been used to test the code. First, the nonlinear steady solution in a plane wave geometry, p = p 0 + δp {1 + tanh [(x -vt)/δ] /2} , with v δp(1 -Γ)/(2 ρ 0 c 2 ) and δ = 2ν ρ 0 c/(δp(1 -Γ))

We thank G. Jourdan and L. Houas for their advice on shock tubes and S. Kerampran, C. Clanet, V. Leroy and all the members of ANR SAMOUSSE for enlightening discussions. M. M. acknowledges a financial support by the DGA.

Appendix A. Wave velocity measurements. -The foam sample is lightened from below with an intense white light. The foam dynamics is then recorded from top using a high speed camera Photron FastCam Sa3. The camera is triggered by the same signal as the pressure sensor and both measures are therefore synchronised.

Images are recorded with the field of view shown on Fig. 5 (bottom). Each image is averaged in the y direction, and the resulting line is used to produce a space-time plot as shown on Fig. 6. At each position r, the first gray level variation is automatically detected by image processing. This first foam motion occurs at time t fr (r) and coincides with the wave arrival at position r. This provides the front wave position r fr (t).

The wave velocity is deduced from the slope obtained in the linear regime, far from the center. Appendix B. Code validation. -The set of equations ( 5) and ( 6) is solved numerically with an implicit scheme implemented in a Matlab code. Two known ana-