Teng Ma

George Bosilca
email: bosilca@eecs.utk.edu

Aurelien Bouteiller
email: bouteill@eecs.utk.edu

Brice Goglin
email: brice.goglin@inria.fr

Jeffrey M Squyres
email: jsquyres@cisco.com

Jack J Dongarra
email: dongarra@eecs.utk.edu

Kernel Assisted Collective Intra-node MPI Communication Among Multi-core and Many-core CPUs

Keywords: MPI, multi-core, many-core, shared memory, NUMA, kernel, collective communication

Shared memory is among the most common approaches to implementing message passing within multicore nodes. However, current shared memory techniques do not scale with increasing numbers of cores and expanding memory hierarchies -most notably when handling large data transfers and collective communication. Neglecting the underlying hardware topology, using copy-in/copy-out memory transfer operations, and overloading the memory subsystem using one-to-many types of operations are some of the most common mistakes in today's shared memory implementations. Unfortunately, they all negatively impact the performance and scalability of MPI libraries -and therefore applications.

In this paper, we present several kernel-assisted intra-node collective communication techniques that address these three issues on many-core systems. We also present a new Open MPI collective communication component that uses the KNEM Linux module for direct inter-process memory copying. Our Open MPI component implements several novel strategies to decrease the number of intermediate memory copies and improve data locality in order to diminish both cache pollution and memory pressure. Experimental results show that our KNEM-enabled Open MPI collective component can outperform state-of-art MPI libraries (Open MPI and MPICH2) on synthetic benchmarks, resulting in a significant improvement for a typical graph application.

I. INTRODUCTION

In recent years, thermic issues have prevented the straightforward performance improvement of increasing processor operating frequency. Maintaining an increase in performance therefore had to be achieved a different way. Parallelism on a chip -in the form of multi-core processors -has been widely adopted as the solution. This trend is even more pronounced when considering the systems of the TOP500 list of supercomputers [START_REF] Meuer | Top500 supercomputing sites[END_REF]. These systems are expected to feature, in the near future, fat many-core nodes composed of one hundred (or more) cores. While increasing the number of cores raises the theoretical peak performance, keeping that many processing units busy requires a significant amount of data to be transferred to and from main memory. The flat memory bus, as featured in many legacy Symmetric Multi Processors (SMP) north-bridge chipsets, is not able to sustain such a bandwidth and request throughput, practically limiting the achievable performance to a tiny fraction of the computing peak -a problem known as the memory wall. To avoid this problem, most recent multi-core designs embrace Non Uniform Memory Access (NUMA) and hierarchical memory interconnects to enable core count scalability and adequate bandwidth between the cores and the memory banks, but at the expense of an excruciating programming complexity.

Message passing has been the dominant programming model in High Performance Computing (HPC) applications for over a decade. However, MPI applications that are oblivious of the NUMA topologies and the associated performance traps, are bound to suffer from unacceptable performance because they generate a load pattern on the memory subsystem that crashes into the memory wall. While the MPI programming model is expressive enough to enable a mapping between the underlying shared memory topology and the application communication pattern, such an approach requires the modification of every code on every platform, defeating one key feature that has empowered the prevalence of MPI: performance portability. Indeed, application developers have come to expect performance portability not just in point-to-point MPI communications, but in MPI's collective communications. Although researchers have been advancing the state of the art in collective communication performance for years, the investigative focus has typically been on the complexity and performance difficulties posed by hierarchical network interconnects. We believe that the collective communication techniques used with hierarchical network interconnects should now also be applied to distributed machines featuring many cores and NUMA architectures.

The usual approach for transporting message payloads between MPI processes across shared memory is based on the copy-in/copy-out algorithm (as illustrated by the SM component in Open MPI [START_REF] Gabriel | Open MPI: Goals, concept, and design of a next generation MPI implementation[END_REF] and the Nemesis [START_REF] Buntinas | Design and evaluation of Nemesis, a scalable, low-latency, message-passing communication subsystem[END_REF] device in MPICH2 [START_REF] Gropp | A highperformance, portable implementation of the MPI message passing interface standard[END_REF]). This algorithm uses an extra pre-allocated shared memory buffer as an exchange zone between processes. Each message is copied to this intermediate zone by the sender process and then copied to the destination buffer by the receiver process. With the rise of multi-core processors, alternative methods of transferring data between processes have emerged in the form of kernel assisted memory copy (such as KNEM [START_REF] Buntinas | Cache-Efficient, Intranode Large-Message MPI Communication with MPICH2-Nemesis[END_REF] and Limic [START_REF] Jin | LiMIC: support for high-performance MPI intra-node communication on linux cluster[END_REF]). These approaches utilize one-copy memory transfers, and have proved beneficial in the context of point-to-point communications [START_REF] Ma | Locality and topology aware intra-node communication among multicore CPUs[END_REF].

As the premise of this work, we have identified three critical issues specific to collective communications between cores in NUMA architectures that prevent fully utilizing the benefits of kernel-assisted copy methods. The first problem arises from the contentions imposed on the root process in one-to-all or all-to-one collective operations. Every process has to wait for the progression of the core hosting the root process in the copy of data to or from the intermediate buffer.

This effect actually prevents any potential acceleration from having multiple cores available to undertake multiple copies simultaneously. Second, many algorithms do not take into account temporal locality, which results in cache-ready data being discarded and then reloaded multiple times. More cache pollution, in turn, leads to a plummeting memory bandwidth as more data lines are reclaimed from the slow and contention-prone memory banks. Last, many implementations ignore topological characteristics such as NUMA and network-style processor interconnects. The blind application of the one size fits all collective algorithm can lead to unnecessary traffic between sockets, potentially overwhelming some memory links while under-utilizing others. These issues can be tackled by 1) extending the feature set of kernelassisted memory copy mechanisms to allow the specification of copy direction and granularity, and 2) reworking the collective algorithms themselves to detect and exploit data locality in NUMA architectures.

In this paper, we propose new MPI collective communication algorithms that take advantage of the NUMA memory subsystem to avoid memory contention and to maximize the overall sustained bandwidth. Our approach is based on the KNEM Linux module -a software mechanism that enables direct memory copying between processes. We investigate several different optimizations to collective algorithms that maximize both parallelism and pipelining, all of which are NUMA topology-aware. A key point in the design is that multiple processes can access the same buffer -or different parts of the same buffer -simultaneously, without the need for more than a single memory copy between processes. Moreover, stream direction control enables the collective algorithm to select sender-writing or receiverreading according to the communication pattern (all-to-one or one-to-all) to avoid the root process bottleneck. Last, our collective algorithms can detect distance between hardware units to build an optimized communication topology that minimizes inter-socket traffic. Each of these approaches are evaluated experimentally on a variety of different hardware setups, exhibiting a better scalability when increasing the number of cores, leading to substantial performance gains on many-core hardware.

The rest of this paper is organized as follows: Section II introduces related work on intra-node collective operations and kernel-assisted memory copy. Section III introduces some key concepts of the KNEM kernel copy framework. Next, Section IV discusses the simultaneous use of KNEMbased collective algorithms with a NUMA topology aware hierarchical layout. Then, Section V describes the linear KNEM collective algorithms: one-to-all (Broadcast and Scatter), all-to-one (Gather), all-to-all (Alltoall and Allgather), and their corresponding implementations in a new collective component for Open MPI. All those algorithms are compared experimentally with state-of-the-art MPI implementations: Open MPI and MPICH2 in Section VI. Finally, Section VII concludes the paper with a discussion of the results.

II. RELATED WORK

The work presented in this paper is the meeting point between two different trends focusing at optimizing MPI communications within shared memory machines. It therefore relates both the optimization of MPI collective operations on shared memory machines and the advances in kernelassisted process-to-process copies.

Several optimizations have been used to maximize throughput of collective communications on shared memory nodes. Most of them have been based on adopting different communication topologies (linear, chain, split binary tree, binomial tree, etc.), [START_REF] Huse | Collective communication on dedicated clusters of workstations[END_REF] and by enabling parallel treatment of the message through pipelining. Both MPICH2 and Open MPI feature many of those algorithms and select among them with highly tuned and optimized switch-points based on the message size. Consequently, they deliver good performance on SMP nodes, even though all those algorithms rely on the double memory copy shared memory transport device. However, this last aspect is greatly challenged by the multiplication of the number of cores in currently deployed fat supercomputer nodes.

One of the most recent of those efforts is due to Richard Graham et al. [START_REF] Graham | MPI support for multi-core architectures: Optimized shared memory collectives[END_REF], who proposed a shared memory-based fanin/fan-out implementation for multi-core MPI collectives, implemented in the Open MPI SM collective component. Their optimization focuses on lightweight synchronization, reducing memory copy times, increasing parallelism by copying messages in a pipeline way, and controlling working set size to fit into caches by building a logical fixed degree tree. This shared memory based method simply takes multicore/many-core as a SMP system and ignores other architecture characteristics, such as NUMA, memory hierarchy, and core distance. The fixed degree tree is built following the logical ranks layout, which cannot always reflect architecture characteristics. With more heterogeneity in modern NUMA multi-core designs, it is hard to optimally tune a shared memory based implementation for different platforms.

Open MPI also features another interesting intra-node collective component named tuned. The fundamental idea of the tuned collectives is to make available several different algorithms, and use a runtime decision to select the best algorithm according to message size, communicator size, and other parameters [START_REF] Fagg | Tuned: A flexible high performance collective communication component developed for Open MPI[END_REF]. As an example, for a Broadcast in the tuned collective component, a binomial algorithm is used to deliver small messages, a split binary tree algorithm is selected for intermediate messages, and large messages are transferred by a pipeline algorithm in which the pipeline size varies with the message and communicator size. However, it is still hard to tune for an unknown platform, even for expert developers. Indeed, there are too many parameters such as pipeline size, thresholds, etc, and any wrong selection might ruin the overall performance of the tuned collectives.

These previous approaches are orthogonal works to the proposed ideas of this paper. They try to maximize the throughput of the collective operation by developing new collective topologies, or selecting, among the available algorithms, the most suitable one. For those approaches to reach their full potential, there is a need to cooperate with another approach to alleviate the penalty due to heterogeneous NUMA architectures, a feature that kernel assisted approaches are able to deliver, if used properly. Unlike those collective algorithms, our component is fully aware of the existence of kernel assisted collective, and uses strategies specifically tailored to maximize the resulting benefits.

Several platform-specific efforts offered single-copy large message communication. For instance BIP-SMP implemented such an optimization for Myrinet based clusters [START_REF] Geoffray | BIP-SMP: high performance message passing over a cluster of commodity SMPs[END_REF]. This idea has spread into most vendor specific HPC stacks, such as Myricom MX, Qlogic IPath, and Bull MPI. Some lightweight kernels enabled an even bigger rework of the model on Cray platforms, thanks to the ability of the operating system to make all processes address spaces accessible to any of them. This unusual feature enabled single-copy RMA-based communication (SMARTMAP [START_REF] Brightwell | Smartmap: operating system support for efficient data sharing among processes on a multi-core processor[END_REF]) which greatly reduces memory copies needed by intra-node message passing, especially for collectives. Recent Linux kernels support the remapping of others' processes memory thanks to the XPMEM module which enables similar optimizations but is restricted to SGI machines.

Lei Chai et al. [START_REF] Jin | LiMIC: support for high-performance MPI intra-node communication on linux cluster[END_REF] introduced a kernel module interface called LiMIC. This kernel-based approach can reduce the number of necessary memory copies to one. KNEM [START_REF]KNEM: High-Performance Intra-Node MPI Communication[END_REF] is another similar kernel module used in MPICH2 and Open MPI. KNEM offers additional features such as an asynchronous copy model, vectorial buffer support, and copy offload on dedicated hardware. This approach has already proved to be valuable to increase point-to-point bandwidth between processes communicating over shared memory [START_REF] Buntinas | Cache-Efficient, Intranode Large-Message MPI Communication with MPICH2-Nemesis[END_REF]. However, beyond the free performance upgrade offered by using more efficient point-to-point operations [START_REF] Moreaud | Optimizing MPI Communication within large Multicore nodes with Kernel assistance[END_REF], hierarchy aware MPI collective components need more control over the underlying memory copy mechanism to reach their full potential. In this paper, we present new collective algorithms that take into account the specificities of kernel assisted memory copies, and require a new feature compared to state of the art kernel assisted copies: directional control of transfers. The details of those novel KNEM features are introduced in the following section.

III. KERNEL ASSISTED MEMORY COPIES

Copy-in/copy-out in a shared memory segment is still the most common approach for transferring data between processes on a shared memory platform. The most prominent drawback is the necessity to copy every data twice. Kernelassisted memory copy alleviates this issue by using system calls to offload the copy to the kernel. Because the kernel has complete access to the memory space of both processes, it can perform the copy from the source buffer in the sender's address space directly to the target buffer in the receiver's address space without the need for an intermediate buffer.

KNEM is an example of such a Linux kernel module that enables high-performance, inter-process, one-copy memory copies. It offers support for asynchronous and vector data transfers. KNEM can also offload memory copies to a hardware DMA engine (such as the Intel I/O AT hardware), if available.

It was not long before the above-cited features raised interest for improving intra-node MPI communications. Open MPI v1.5 includes KNEM support in its shared memory point-to-point communications component. MPICH2 v1.1.1 uses KNEM in the DMA LMT to improve large message performance within a single node. The general operating principles of the integration between KNEM and MPI point-to-point messages are the following (more details can be found in [START_REF] Buntinas | Cache-Efficient, Intranode Large-Message MPI Communication with MPICH2-Nemesis[END_REF]): [START_REF] Meuer | Top500 supercomputing sites[END_REF] The sender process declares a send buffer to KNEM. The kernel module saves the list of virtual segments contained in the buffer and associates them with a unique cookie. [START_REF] Gabriel | Open MPI: Goals, concept, and design of a next generation MPI implementation[END_REF] The sender passes the cookie to the process which is interested in this buffer by an out-ofband transfer. (3) The receiver gives this incoming cookie to KNEM along with a receive buffer. (4) The KNEM module copies the data from the send buffer to the receive buffer within the kernel.

The security model of this strategy is similar to System V IPC shared memory segments. Declaring a memory buffer to the KNEM driver makes it available to any other process. However, other memory regions cannot be accessed unless explicitly declared to KNEM as well. A malicious user modifying a cookie value would either get an invalid parameter error or get a valid access to a previously properly declared buffer.

A. Issues with MPI Collective Operations

While the beneficial effect of KNEM on point-topoint performance also translates into collective improvements [START_REF] Moreaud | Optimizing MPI Communication within large Multicore nodes with Kernel assistance[END_REF], we have identified a series of additional optimizations that further boost collective communication performance. They require that the collective component has more control over the movement of data (vs. simply using MPI point-to-point primitives): (1 Because control of the kernel module is delegated to the point-to-point MPI message passing engine, using inter-process kernel-assist memory copies results in the same data region being registered multiple times when sent to different destination processes. The overhead of synchronizing and exchanging cookies therefore cannot be amortized. The collective component knows when the same buffer is used with multiple recipients, and can therefore eliminate redundant registrations. (2 Many collective algorithms exhibit a one-to-all communication pattern. Such a communication pattern stresses the root process resources. Although in kernel execution space when using KNEM, only the core hosting the root process is performing all the data movements. This actually serializes all memory copies and limits the progression of the collective to the speed of only a single core (i.e., where the root process is executing), even though multiple cores are available (and probably waiting for the collective operation to progress). Attempting to alleviate this issue by inverting the point-to-point algorithm -i.e., having the receiver(s) make the copies -simply results in degrading all-to-one communication patterns, instead. Thus, it is desirable for the collective algorithm to be able to express the direction of data transfer, independently of concurrent point-to-point strategies.

B. KNEM Direction and Granularity Control

Implementing collective operations directly on top of the original KNEM interface wastes system resources. Collective patterns such as one-to-all or all-to-one would declare the same multiply-accessed memory buffer multiple times and pass multiple KNEM cookies between processes for synchronization. To solve this issue, we introduced, (available since KNEM 0.7), an extended programming interface designed to address the needs of collective operations. Instead of only offering a point-to-point send-receive interface, KNEM now offers the ability to register persistent memory regions and access them multiple times from different processes. Such accesses may either touch all or only part of the region, enabling the actual copy a single message at once or as multiple chunks simultaneously. This model avoids wasting system resources and reduces the overhead of synchronizing processes when creating and passing region identifiers.

Another extension added to KNEM is the ability to read (receiver-reading) or write (sender-writing) to each region, enabling effective direction control of the data transfer.

Direction control is a very important feature to unleash the performance of collective operations. The effective direction control of the copies can be decided by the collective component to match the communication pattern (one-to-all or all-to-one), with the goal of maximizing the number of cores participating to the progression of the data transfers, and parallelize the progression of the collective algorithm as much as possible.

Those two novel features introduced into the KNEM kernel module are used by our KNEM collective component in Open MPI. Unlike previous components that only used kernel-assisted copies simply through the point-to-point interface, this component takes advantage of directional control and persistent registrations to further increase the performance achieved on collective communications. A significant amount of work has been done over the past decade toward improving the collective communication performance by taking into account the network features. Some algorithms take advantage of the low latency of some specific high-performance, while others capitalize on the bandwidth capabilities. Additionally, the network topology (butterfly, torus) has been another important factor in redesigning collective communications algorithms. While this work is significant, and has clearly influenced our approach, one has to keep in mind the tight conditions required for collective algorithms in shared memory. The memory access latency and bandwidth are important, but not more important than the topology of the links to and from the memory banks. Therefore, the algorithms designed in this context are significantly more complex than the usual two-level hierarchical collective algorithms [START_REF] Kielmann | MagPIe: MPI's collective communication operations for clustered wide area systems[END_REF], and take into account the topology as well as the memory accesses performance, that usually exhibit several order of magnitude performance differences for different levels in the hierarchy.

IV. TOPOLOGY AND NUMA AWARE COLLECTIVES

In order to explain how the collective algorithm will be affected by the hardware architectural features, let us take our large NUMA node (IG) as an example. This machine consists of 8 NUMA nodes, each of them containing a sixcore AMD Opteron processor and 16GB of local memory. 8 NUMA nodes are interconnected by AMD HyperTransport (HT). The Figure 1 shows how a 48 processes broadcast will unfold on this large NUMA node. Processes are split into 8 sets according to their NUMA locality information, which means processes within the same socket and NUMA node are in the same set. A two levels' tree is then built accordingly, one process per NUMA node will belong to the first tree level (the green background circles), while all the remaining processes per NUMA node will belong to the second level (the blue background circles), and behave as leafs to the tree. This tree structure reflects the architecture topology such as core distances and relationships which can be retrieved thanks to the Hardware Locality [START_REF] Broquedis | HWLOC: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF] software.

Dividing the processes based on the topology and NUMA information has the advantage of reducing inter-socket traffic since a single data transfer is performed towards each NUMA node. Moreover, since a cache is shared between all processes inside the same NUMA node, multiple copies between processes in the same set benefit from cache hits. However, one of the disadvantages of such an algorithm is the reduction in the degree of parallelism between the memory copies toward NUMA nodes, and between the leaf nodes and their corresponding intermediary node (a leaf process cannot start a memory copy until the intermediary process received the entire data). To alleviate this strong synchronization requirement, and therefore minimize the unnecessary waiting on the leaf nodes, we divide the data in several segments and pipeline the operations corresponding to each one of these segments (pink arrows).

V. IMPLEMENTATION

A. Framework

We implemented the KNEM based collectives as a new component named KNEM coll in Open MPI. Open MPI is based on a flexible component architecture, collective algorithms are placed under the COLL framework (as depicted in Figure 2). Multiple collective components are available (Tuned, Basic, SM, and our proposed KNEM coll), and can be selected at runtime. In Open MPI ' [START_REF] Squyres | The component architecture of Open MPI: Enabling Third-party Collective Algorithms[END_REF]. By default, the Tuned collective component and the SM BTL component are teamed to provide for MPI collective operations, resulting in the default setup to use the best collective algorithms with the copy-in/copy-out point-to-point transport. In previous works the SM/KNEM BTL have been introduced. This BTL enhances shared-memory point-to-point operations by using the KNEM driver for large messages. While this is not a default setup, the Tuned collective component can be teamed with this BTL to benefit from KNEM speedup to some extent for large messages.

Unlike other collective components, the KNEM collective component uses the shared memory BTL only as an out of band channel for synchronization or delivering "cookie". All data movements are directly handled inside the collective component by resorting to direct calls to the KNEM kernel driver. While the KNEM kernel module is used directly, the collective algorithms themselves remain implemented in user-space. Direct calls to KNEM within the collective component are required to have enough flexibility to express the new algorithms intended to avoid unnecessary buffer registrations, handle directional copies, and fragment the messages to establish a pipeline suitable for the most complex memory hierarchies. The KNEM RMA API that we introduced in KNEM (see Section III) is used to create memory regions dynamically and read or write all or part of them multiple times whenever needed. However, trapping into kernel mode has a non-negligible overhead (about 100 ns on modern processors) when delivering small messages. So we only consider KNEM for optimizing collectives for intermediate and large messages (larger than 16KBytes). We started by adapting the most useful algorithms first (Gather(v), Scatter(v), AlltoAll(v), AllGather(v) and Broadcast, v stands for vector variants of corresponding operations). For smaller messages, or unimplemented collective calls, the operation is delegated to the regular Open MPI component.

B. Rooted Operations: Broadcast, Scatter and Gather

The implementation of KNEM Broadcast is a straightforward adaptation of the KNEM point-to-point model: 1)

The root process declares a send buffer to KNEM and gets the corresponding cookie in return. 2) It passes the cookie to all non-root processes in the communicator through an out-of-band transfer (Open MPI SM BTL). 3) Each receiver process passes the incoming cookie to KNEM along with a receive buffer. 4) KNEM triggers a memory copy from the send buffer to receive buffer within the kernel. The copy is performed by each receiver core in parallel. 5) Each receiver process sends back a synchronization message to root process after the completetion of KNEM copy. 6) After the root process receives all synchronization messages from non-root processes, it deregisters the send buffer from the KNEM driver. The KNEM Broadcast also features a hierarchical pipelining algorithm (as sketched in Section IV), that can be turned on or off depending on the properties of the hardware. The topology mapping is static for now, but will be dynamic in future works.

The KNEM Scatter implementation is overall similar to the broadcast, except that each non-root process reads only parts of the root buffer. Their starting offset is calculated from their ranks and the data type.

The KNEM Gather consists of the opposite communication pattern of Scatter. Therefore it can benefit from the direction control of the new KNEM version to declare the root process memory as a write buffer. Unlike the regular Gather, all non-root processes can write simultaneously to that buffer.

C. AllGather and Alltoall

The KNEM AllGather is an assembly of a gather followed by a broadcast. During the first step, all processes do a KNEM Gather operation to rank 0. In the second stage, the root process performs a KNEM Broadcast. This method is far from being optimal, because it puts the memory controller of the the core hosting the root process under a high pressure. Although it could be perfected, it is a simple and straightforward way to capitalize on the improvements achieved in Broadcast and Gather. Unlike one-to-all operations, in the AlltoAll communication pattern, cache reuse plays little to no role. Therefore, our AlltoAll algorithm focuses on avoiding bandwidth sharing, by rotating the access pattern so that at any instant, a core is sending and receiving exactly one fragment of data. First, each process declares its send buffer to KNEM and gets back the corresponding cookie. Those cookies are exchanged by doing an out of band AllGather operation based on shared memory(not KNEM AllGather); it is necessary to pre-allocate an integer array the size of the communicator to store cookies from all other processes. A loop of KNEM copies is performed to fetch the corresponding messages from other nodes (receiver-reading). Each process offsets the starting point of this loop in a round-robin manner. Finally, each process deregisters memory buffer from the KNEM driver after a barrier operation. Figure 3 gives an example of this communication scheduling; as one can see, for every step of the algorithm (marked between parenthesis on the right buffer), the memory belonging to a particular sender is accessed only once, and the workload is evenly spread on the entire duration of the collective.

VI. EXPERIMENTAL EVALUATION

A. Experimental Conditions

Our experimental platform is composed of four multicore/many-core machines that cover the spectrum of typical current commodity high-performance computing nodes, but also expected designs for the upcoming years. They feature Intel and AMD processors, and represent a wide variety of setups, from SMP to massively NUMA machines. As our approach does not impact the performance of internode communications, our experiment focuses on intra-node communications performance. Zoot is a 16 core machine with 32GB of memory. The system has four sockets with a quad-core 2.40 GHz Intel Xeon Tigerton E7340 featuring 4 MB L2 caches shared between pairs of cores. A single SMP memory controller in the north-bridge chipset connects all the sockets with the global shared memory. Dancer is an 8 core machine with 4GB of memory. The system is composed of two sockets populated with a quadcore 2.27 GHz Intel Xeon Nehalem-EP E5520 with 8 MB L3 caches and 2 GB of memory on each NUMA socket. Hyper-threading is disabled in the configuration. IG is a 48 cores machine with 128GB of memory. The system is composed of 8 sockets with a six-core 2.8 GHz AMD Opteron 8439 SE, 5 MB L3 caches and 16 GB memory per NUMA node. The sockets are further divided as two sets of 4 sockets on two separate boards connected by a low performance interlink. Saturn is a 16 core machine with 64GB of memory. It is composed of two sockets with an octo-core 2.00 GHz Intel Xeon Nehalem-EX X7550, 18 MB L3 caches and 32 GB memory on each NUMA socket. Hyper-threading is enabled but not used.

Software setup includes KNEM version 0.9.2 [START_REF]KNEM: High-Performance Intra-Node MPI Communication[END_REF]. The Intel MPI benchmark suite IMB-3.2 was used to assert the difference between the collective components with the offcache option in enabled to avoid cache reuse. Open MPI version 1.7a1 and MPICH2-1.3.1, both properly tuned, are compared to our approach. Tuning parameters between all components based on Open MPI are identical, unless stated. On a particular machine, the mapping between physical cores and MPI processes is identical, regardless of the MPI implementation used. Due to the large number of combinations, we restrict the discussion only to the most meaningful algorithms, but still covering the entire spectrum of collective patterns.

Because our own component (referred to as KNEM-Coll from now on) is inside Open MPI, we undergo a deeper comparison with its default collective component: the Tuned component. It does not use kernel assisted copies by default, but uses the SM device that relies on copy-in/copy-out; results obtained with this setup are called Tuned-SM. To assert clearly what are the extra benefits of our approach, we also present the results obtained when simply using the Tuned component on top of KNEM point-to-point messages (Tuned-KNEM). MPICH2-SM and MPICH2-KNEM are similar to Tuned-SM and Tuned-KNEM, with respectively shared-memory or KNEM as the underneath point-to-point communication transport. The MPICH2 KNEM implementation is broken for messages of size beyond 64 KB on the Saturn platform; we could not obtain a revised version from the author at the date of publication of this paper.

To ease the comparison, our KNEM collective component is taken as the reference point for performance, and the execution time of other implementations are normalized against it. The smaller these normalized values, the better the performance of the corresponding collective component. hierarchical pipelined KNEM Broadcast on the large NUMA machine IG. The pipeline sizes range from 4KB to 2MB. The execution time of different pipeline sizes is normalized to the execution time of hierarchical KNEM Broadcast without pipeline. Compared with the linear KNEM broadcast, the hierarchical approach itself contributes a 2.2× to 2.4× speedup on this large NUMA node. And the pipelining provides hierarchical KNEM Broadcast extra up to 1.25× speedup. One can see that except using too small pipeline size (4KB), hierarchical KNEM Broadcast always can get significant benefits from pipelining. Selecting a suitable pipeline size can be a challenging problem for the hierarchical pipelined strategies. A too small pipeline size induces more synchronization between each segment and makes the transfer efficiency suffer due to KNEM copy startup overhead when delivering small messages (4KB and 8KB in Fig 4), while a too large pipeline size (2MB in Fig 4) leads to a long initialization time for the pipeline algorithm to take effect. One can see that the best pipeline size is 16KB for intermediate message size (smaller than 2MB), and 512KB for large message sizes. In the rest of this paper, we settled the pipeline size according to this tuning on IG: 16KB for intermediate message size and 512KB for large message size.

B. Hierarchical effect and tuning the pipeline size

C. One-to-all and All-to-One Operations

Figure 5 shows the performance comparison of the Broadcast implementations on all platforms. The KNEM Broadcast outperforms other collective components in all cases. Compared with Open MPI's best collective component, the KNEM Broadcast can provide a speedup of about 1-2.5× on Zoot, 1.2-2.8× on Dancer, 1-1.8× on Saturn and 1.5-2.1× on IG.

Figure 6 presents the performance comparison of the Gather operation. The linear KNEM Gather tremendously outperforms all other components in all cases. Compared with the best Gather in Open MPI and MPICH2, the maximum speedup, thanks to KNEM collective component, is 3.1× on Zoot, 2.2× on Dancer, 2.6× on Saturn, and 3.2× on IG. KNEM Gather and Scatter operations are very similar, except from the different copy direction: sender-writing for Gather and receiver-reading for Scatter. Consequently, those two algorithms exhibit very similar performance profiles and the Scatter results are not presented here. Compared with Open MPI's best Tuned Scatter implementation, the maximum speedup of KNEM Scatter is about 3× on Zoot, 2× on Dancer, 4× on Saturn, and 4× on IG.

The KNEM collective component has a huge performance speedup in these "rooted" collective operations, thanks to unleashing parallel access to the buffer of the root process. Compared with the approach of indirectly using KNEM copy as underneath communication (Tuned-KNEM), the KNEM collective component can provide more reliable improvement in all cases, benefiting from KNEM drivers' new features we mentioned at section III.

D. All-to-all Operations

Figure 7 shows the performance comparison for the AlltoAllv collective operation. The performance of these all-to-all collective operations is restricted by the press over memory buses due to delivering too many messages simultaneously. As a consequence, the relative performance benefits are smaller when compared to Tuned-KNEM than for the one-to-all or all-to-one rooted operations presented in the previous paragraphs. However, compared with Tuned-SM based on the shared memory approach, the KNEM AlltoAllv can still show significant improvement, with a maximum speedup of 2× on Zoot, 1.9× on Dancer, 1.25× on Saturn, and 2.7× on IG. Finally, Figure 8 presents the performance of the collective components on the AllGather operation. On the SMP machine (Zoot), and the medium size NUMA machines (Dancer, Saturn), the KNEM AllGather has the best performance among all collective components except some medium size messages on Dancer and Saturn. On the large NUMA node (IG), Tuned-KNEM performs better than KNEM AllGather by up to 25%. The loss of KNEM AllGather's performance on large NUMA nodes lies in the KNEM AllGather operation not being optimized as explained in Section V-C. The operation is split into two separated stages. Although the KNEM Gather and Broadcast are optimized in each stage, overlapping between these two stages is eliminated by this simple concatenation of the KNEM Gather and Broadcast algorithms. And the selected root process forms a single point in KNEM AllGather implementation, forcing the throughput of AllGather operation restricted by the limited memory bandwidth of the NUMA node owning the root process. This is also the reason why the KNEM AllGather's performance suffers more on the large NUMA nodes than on SMP or small NUMA nodes. However, even on large NUMA node IG, the KNEM All-Gather still performs better than Tuned-SM and MPICH2-SM, which are based on shared memory approach. The implementation of these two all-to-all collective operations (e.g. AlltoAllv and AllGather) benefits greatly from the adoption of KNEM copy, thanks to reducing memory copies and cache pollution in these communication-intensive operations. In the next release of the KNEM collective component, we will borrow some ideas in Tuned collective components such as adopting a ring-style algorithm to distribute memory accesses evenly across memory controllers in the KNEM AllGather, especially on large NUMA nodes.

E. Application Performance

To evaluate the impact of the improvement due to using KNEM collective operations on real application performance, we use the ASP [START_REF] Plaat | Sensitivity of parallel applications to large differences in bandwidth and latency in two-layer interconnects[END_REF] application, a parallel implementation of the Floyd-Warshall algorithm used to solve the all-pairs-shortest-path problem. The main MPI collective operation used in this application is MPI Bcast. We tested this application on two machines: Zoot and IG, the two extreme platforms regarding the degree of complexity of the core hierarchy. The problem is scaled to match the available memory; the matrix size is 16384 2 on Zoot and 32768 2 on IG (32bits integers). Matrices are distributed by rows across all the available cores. The MPI Bcast operation is called 16384 times (64 KB message) on Zoot and 32768 times (128 KB message) on IG. The KNEM collective component uses the linear KNEM algorithm on Zoot and the hierarchical pipelined algorithm on IG. All tests use the same mapping between cores and processes. Table I presents the application execution time of ASP when using different collective components. The improvement is the relative difference between the best performing MPI library (between Open MPI and MPICH2) and our KNEM collective components. By using the KNEM collective components, the application can see a significant improvement in the time it spends doing Broadcast operations, with the improvement of 82% on Zoot and 33% on IG. One can notice that the performance improvement of the Broadcast only on the SMP machine is even more pronounced than for the synthetic benchmark, because unlike the synthetic benchmark, the application does not systematically invalidate the cache before performing the operation. As a consequence of the shorter time spent in the collective operations, the overall application runtime is improved when compared to other optimized collective components, with an improvement of 5% on Zoot and 2% on IG.

VII. CONCLUSION

In this paper, we showed that an MPI implementation can successfully take advantage of new features in the kernel, namely kernel assisted memory copies between process spaces, to greatly improve the performance of shared memory collective communication while successfully hiding the complexity of the NUMA characteristics of modern manycore nodes. Three challenges specific to many-core systems have been addressed to meet that goal: NUMA complexity, cache pollution, and collective algorithm sequentialization.

The main contributions of this paper are: (1) extend kernel-assisted memory copy operations to include direction and granularity control, thereby giving more flexibility and data movement control to the collective component implementation, and (2) embed the use of kernel-assisted copies directly into collective algorithms rather than relying on kernel-assisted speedups from MPI point-to-point-based primitives. (3) Our collective component also maps the inherent parallelism of MPI collective communications on to the specific characteristics of NUMA multi-core systems. For example, the root process can offload memory copies to non-root processes in order to parallelize operations with multiple receivers. Finally, the pipelining algorithm strategically overlaps the latency incurred by transferring first to the upper levels of the NUMA hierarchies with data transfers to the leaf nodes in the collective topology. Experimental results show that our approach outperforms all other types of optimizations on many-core systems, including the Tuned Open MPI component (even though Tuned benefits from KNEM-enabled point-to-point operations). These performance improvements translate into a substantial gain for our showcase application.

Figure 1 .

 1 Figure 1. Progression of the hierarchical pipeline KNEM Broadcast

Figure 3 .

 3 Figure 3. An example of copy sequence for KNEM AlltoAll on four processors.

Figure 4 .

 4 Figure 4.Performance comparison between linear KNEM Broadcast, hierarchical KNEM without pipeline, and different pipeline sizes in the hierarchical pipelined KNEM Broadcast on the IG platform. Results are normalized to the runtime of hierarchical KNEM Broadcast without pipeline (lower is better).

Figure 4 Figure 5 .

 45 Figure 4 presents the effect of the pipeline size on the

Figure 6 .

 6 Figure 6. Performance comparison of Gather operations between Open MPI Tuned components, MPICH2 and the KNEM collective component. Results are normalized to the runtime of the KNEM Gather (lower is better).

Figure 8 .

 8 Figure 8. Performance comparison of AllGather operations between Open MPI Tuned components, MPICH2 and the KNEM collective component. Results are normalized to the runtime of the KNEM AllGather (lower is better).

  s component architecture, different collective components can use different point-to-point communication components underneath (the

		MPI API
		COLL Framework
	KNEM coll	Basic coll	Tuned coll	SM coll	...
		Point-to-point
		Communication (BTL)
	SM/KNEM BTL	SM BTL
	KNEM driver		Shared Memory
		Operating System
	Figure 2. Open MPI collective communication framework
	byte transfer layer, BTL)


Hugues Garnier

Data-based continuous-time modelling of dynamic systems

Data-based continuous-time model identification of continuous-time dynamic systems is a mature subject. In this contribution, we focus first on a refined instrumental variable method that yields parameter estimates with optimal statistical properties for hybrid continuous-time Box-Jenkins transfer function models. The second part of the paper describes further recent developments of this reliable estimation technique, including its extension to handle non-uniformly sampled data situation, closed-loop and nonlinear model identification. It also discusses how the recently developed methods are implemented in the CONTSID toolbox for Matlab and the advantages of these direct schemes to continuous-time model identification.

I. INTRODUCTION

The identification of continuous-time (CT) models is a problem of considerable importance that has applications in virtually all disciplines of science. Early research on this topic focussed on identification of CT models from CT data (see e.g. [START_REF] Young | In flight dynamic checkout -a discussion[END_REF], [START_REF] Young | The determination of the parameters of a dynamic process[END_REF]). Subsequently, however, rapid developments in digital data acquisition and computers have resulted in attention being shifted to the identification of discretetime (DT) models from sampled data, as documented in many books (see e.g [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF], [START_REF] Söderström | System Identification, Series in Systems and Control Engineering[END_REF] and [START_REF] Ljung | System Identification. Theory for the User[END_REF]). Much less attention has been devoted to CT modelling from DT data and many practitioners appear unaware that such alternative methods not only exist but may be better suited to their modelling problems.

In order to identify a continuous-time model from timedomain sampled data, two main time-domain approaches are possible. In the first, 'indirect' approach, a DT model is identified first using DT model identification methods, and this is then converted into a CT model using a standard algorithm for discrete to continuous-time conversion. In the second, 'direct' approach the CT model is identified directly from DT data. Direct data-based CT modelling is often incorrectly presented as being too complicated but, as we will see, the approaches are straightforward, reliable and have proven useful in many practical applications. These approaches have recently regained interest showing better performance than indirect approaches for both linear and nonlinear models, see e.g. [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF], [START_REF] Rao | Identification of continuous-time systems: direct or indirect?[END_REF], [START_REF] Rao | Identification of continuous-time systems[END_REF], [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF], [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF]. The main motivations for identifying CT models directly from sampled data have been recently discussed in [START_REF] Garnier | Direct Identification of Continuous-time Models from Sampled Data: Issues, Basic Solutions and Relevance. In Identification of continuous-time models from sampled data[END_REF] (see also the Conclusions Section in this paper). Exhaustive reviews of direct estimation methods can be found in [START_REF] Young | Parameter estimation for continuous-time models -a survey[END_REF], [START_REF] Garnier | Continuous-time model identification from sampled data. Implementation issues and performance evaluation[END_REF], [START_REF]Identification of Continuous-time Models from Sampled Data[END_REF] and [START_REF] Rao | Identification of continuous-time systems[END_REF]. Amongst the available identification approaches for CT input-output models, the interest for instrumental variable (IV) methods has been growing in the last years [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF], H. Garnier is with Centre de Recherche en Automatique de Nancy (CRAN), Nancy-University, CNRS, BP 70239, 54506 Vandoeuvre-les-Nancy Cedex, France hugues.garnier@uhp-nancy.fr [START_REF] Young | The refined instrumental variable method: unified estimation of discrete and continuous-time transfer function models[END_REF], [START_REF]System Identification, Environmetric Modelling and Control System Design[END_REF]. The main reason of this increasing interest is that IV methods offer similar performance as extended least square (LS) methods or other prediction-error-minimization (PEM) methods (see [START_REF] Rao | Identification of continuous-time systems: direct or indirect?[END_REF], [START_REF] Ljung | Experiments with identification of continuous-time models[END_REF]) and provide consistent results even for an imperfect noise structure which is the case in most practical applications. The IV schemes considered here present the major advantages over PEM methods to be much less sensitive to the initialisation stage (see [START_REF] Rao | Numerical illustrations of the relevance of direct continuous-time model identification[END_REF], [START_REF] Rao | Identification of continuous-time systems: direct or indirect?[END_REF], [START_REF] Ljung | Initialisation aspects for subspace and output-error identification methods[END_REF]) . These IV approaches lead to optimal estimates in the linear time-invariant case if the system belongs to the model set defined. This paper concentrates on a reliable Instrumental Variable (IV)-based estimation method in particular, and presents the latest developments, including its use for closedloop and nonlinear model identification.

II. RIVC FOR CT LINEAR MODELS

We focus on a statistically optimal method for the identification of continuous-time hybrid Box-Jenkins (BJ) transfer function models from discrete-time data [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models. In Identification of continuous-time models from sampled data[END_REF]. Here, the model of the dynamic system is estimated in continuoustime, differential equation form, while the associated additive noise model is estimated as a discrete-time, autoregressive moving average (ARMA) process. This refined instrumental variable method for continuous-time systems (RIVC) was first developed in 1980 by Young and Jakeman [START_REF] Young | Refined instrumental variable methods of time-series analysis: Part III, extensions[END_REF] and its simplest embodiment, the simplified RIVC (SRIVC) method, has been used successfully for many years, demonstrating the advantages that this stochastic formulation of the continuoustime estimation problem provides in practical applications (see, e.g., some recent such examples in [START_REF] Young | Data-based mechanistic modeling of engineering systems[END_REF], [START_REF] Young | Identification and estimation of continuous-time, data-based mechanistic models for environmental systems[END_REF]).

However, the 'simplification' that characterises the name of the SRIVC method is the assumption, for the purposes of simplicity and algorithmic development, that the additive noise is purely white in form. Such an approach is optimal under this assumption and the inherent instrumental variable aspects of the resulting algorithm ensure that the parameter estimates are consistent and asymptotically unbiased in statistical terms, even if the noise happens to be coloured. However, the SRIVC estimates are not, in general, statistically efficient (minimum variance) in this situation because the prefilters are not designed to account for the colour in the noise process.

The hybrid RIVC estimation procedure follows logically from the refined instrumental variable (RIV) method for discrete-time models, first developed within a maximum likelihood (ML) context by Young in 1976 [START_REF] Young | Some observations on instrumental variable methods of time-series analysis[END_REF] and comprehensively evaluated by Young and Jakeman [START_REF] Young | Refined instrumental variable methods of time-series analysis: Part I, SISO systems[END_REF], [START_REF] Jakeman | Refined instrumental variable methods of time-series analysis: Part II, multivariable systems[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF].

The RIV algorithm involves concurrent DT noise model estimation and uses this estimated noise model in the iterative-adaptive design of statistically optimal prefilters that effectively attenuate noise outside the passband of the system and prewhiten the noise remaining within the bandpass. Similarly motivated prefilters are utilised in the RIVC algorithm but they also provide a very convenient way of generating the prefiltered derivatives of the input and output variables, as required for CT model estimation.

The alternative hybrid form of the continuous-time transfer function model is considered here for two reasons. First, the approach is simple and straightforward: the theoretical and practical problems associated with the estimation of purely stochastic, continuous-time CAR or CARMA models are avoided by formulating the problem in this manner. Second, as pointed out above, one of the main functions of the noise estimation is to improve the statistical efficiency of the parameter estimation by introducing appropriately defined prefilters into the estimation procedure. And, as we shall see, this can be achieved adequately on the basis of hybrid prefilters defined by reference to discrete-time AR or ARMA noise models.

A. Problem Formulation

For simplicity of presentation, the formulation and solution of the CT estimation problem will be restricted to the case of a linear, single-input, single-output system. It is assumed that the input u(t) and the noise-free output x(t) are related by the following constant coefficient, differentialdelay equation,

x (n) (t) + a o 1 x (n-1) (t) + • • • + a o n x(t) = b o 0 u (m) (t -τ ) + • • • + b o m u(t -τ ) (1)
where x (i) (t) denotes the ith time derivative of the continuous-time signal x(t) and τ is a pure time delay in time units. This is often assumed to be an integer number related to the sampling time: i.e., τ = n k T s but this is not essential: in this CT environment, 'fractional' time delays can be introduced if required (e.g., see [START_REF] Mahata | Identification of continuous-time Box-Jenkins models with arbitrary time-delay[END_REF], [START_REF] Young | The data-based mechanistic approach to the modelling, forecasting and control of environmental systems[END_REF]). For simplicity, the time delay will not be considered in the following analysis but it can be accommodated straightforwardly if identified from the data. Equation (1) can also be written in the following compact transfer function (TF) form,

x(t) = G o (p)u(t) = B o (p) A o (p) u(t) (2)
with

B o (p) =b o 0 p m + b o 1 p m-1 + • • • + b o m , (2a) A o (p) =p n + a o 1 p n-1 + • • • + a o n , n ≥ m (2b)
where x(t) is the deterministic output of the system; p is the differential operator, i.e., p i x(t) = d i x(t) dt i ; B o (p) and A o (p) are assumed to be coprime; and the system is asymptotically stable. It is assumed that the input signal {u(t), t 1 < t < t N } is applied to the system and this gives rise to an output signal {x(t), t 1 < t < t N }.

In order to obtain high-quality statistical estimation results, it is vital to consider the inevitable errors that will affect the measured output signal. It is assumed here that x(t) is corrupted by an additive, coloured measurement noise ξ(t), so that the complete equation for the data-generating system, denoted by S, can be written in the form,

S : y(t) = G o (p)u(t) + H o (p)e o (t) (3)
or, in the alternative decomposed form that is more appropriate in the present context

S      x(t) = G o (p)u(t) ξ(t) = H o (p)e o (t) y(t) = x(t) + ξ(t) (4)
where H o (p) is stable and invertible, while e o (t) is a zeromean, continuous-time white noise source, which is assumed to be uncorrelated with the input u(t). Finally, if the additive coloured noise ξ(t) has rational spectral density, then a suitable parametric representation is the following continuoustime, autoregressive moving average (CARMA) model

ξ(t) = H o (p)e o (t) = C o (p) D o (p) e o (t) (5)
where C o (p) and D o (p) are suitably defined polynomials in the p operator. Of course, in most practical situations, the input and output signals u(t) and y(t) will be sampled in discrete time. In the case of uniform sampling, at a constant sampling interval T s , these sampled signals will be denoted by u(t k ) and y(t k ) and the output observation equation then takes the form,

y(t k ) = x(t k ) + ξ(t k ) k = 1, • • • N (6)
where x(t k ) is the sampled value of the unobserved, noisefree output x(t). The objective is then to identify a suitable model structure for (4) and estimate the parameters that characterise this structure, based on these sampled input and output data

Z N = {u(t k ); y(t k )} N k=1 .
Given the discrete-time, sampled nature of the data, an obvious assumption is that the discrete-time, coloured noise associated with the sampled output measurement y(t k ) has rational spectral density and so can be represented by a discrete-time ARMA model. The model set to be identified and estimated, as denoted by M with system (G) and noise (H) models parameterised independently, then takes the form,

M : {G(p, ρ), H(q, η)} (7)
where ρ and η are parameter vectors that characterise the system and noise models, respectively. In particular, the system model is formulated in continuous-time terms

G : G(p, ρ) = B(p, ρ) A(p, ρ) = b 0 p m + b 1 p m-1 + • • • + b m p n + a 1 p n-1 • • • + a n (8)
and the associated model parameters are stacked columnwise in the parameter vector,

ρ = a 1 • • • a n b 0 • • • b m T ∈ R n+m+1 (9)
while the noise model is in discrete-time form

H : H(q, η) = C(q -1 , η) D(q -1 , η) = 1 + c 1 q -1 + • • • + c q q -q 1 + d 1 q -1 + • • • + d p q -p (10
) where q -r is the backward shift operator, i.e., q -r y(t k ) = y(t k-r ) and the associated model parameters are stacked columnwise in the parameter vector,

η = c 1 • • • c q d 1 • • • d p T ∈ R p+q (11)
Consequently, the noise TF takes the usual ARMA model form

ξ(t k ) = C(q -1 , η) D(q -1 , η) e(t k ) e(t k ) ∼ N (0, σ 2 ) (12)
where, as shown, e(t k ) is a zero-mean, normally distributed, discrete-time white noise sequence. The structure S does not specify any common factors in the plant (G o ) and noise (H o ) components, so that these models can be parameterised independently. More formally, there exists the following decomposition of the parameter vector θ for the whole hybrid model,

θ = ρ η (13)
such that the model equations can be written in the form

M      x(t) = G(p, ρ)u(t) ξ(t k ) = H(q, η)e(t k ) y(t k ) = x(t k ) + ξ(t k ) (14)
This model is considered as a hybrid Box-Jenkins model because of its close relationship to the DT model considered in great detail by Box and Jenkins in their seminal book on time-series analysis, forecasting and control [START_REF] Box | Time Series Analysis Forecasting and Control[END_REF] and used as the basis for the development of the original RIVC algorithm [START_REF] Young | Refined instrumental variable methods of time-series analysis: Part III, extensions[END_REF]. Alternatively, the model can be written in the following vector terms

M      x (n) (t) = ϕ T (t)ρ ξ(t k ) = ψ T (t k )η + e(t k ) y(t k ) = x(t k ) + ξ(t k ) (15)
where,

ϕ T (t) = -x (n-1) (t) • • • -x(t) u (m) (t) • • • u(t) (15a) ψ T (t k ) = [-ξ(t k-1 ) • • • -ξ(t k-p ) e(t k-1 ) • • • e(t k-q )] (15b)
For the purposes of identification, the order of this singleinput model (with the pure time delay τ now added for completeness) is denoted by [n m τ p q] and the complete identification problem can now be stated as follows:

Based on N uniformly sampled measurements of the input and output,

Z N = {u(t k ); y(t k )} N k=1
, identify the orders n, m, p and q of the polynomials in the system and noise TF models, as well as any pure time delay τ , and estimate the parameter vector θ in (13) whose parameters characterise these polynomials.

B. Optimal RIVC Estimation: Theoretical Motivation

The RIVC algorithm derives from the RIV algorithm for DT systems. This was evolved by converting the maximum likelihood (ML) estimation equations to a pseudo-linear form [START_REF] Solo | Time Series Recursions and Stochastic Approximation[END_REF] involving optimal prefilters [START_REF] Young | Some observations on instrumental variable methods of time-series analysis[END_REF], [START_REF] Young | Refined instrumental variable methods of time-series analysis: Part III, extensions[END_REF], [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF]. A similar analysis can be utilised in the present situation because the problem is very similar, in both algebraic and statistical terms. However, to conserve space, the discussion here will be restricted to a simpler development of the RIVC algorithm and we leave the interested reader to consult with these earlier references for details of the ML analysis.

1) The Hybrid Box-Jenkins Estimation Model: Following the usual prediction error minimisation (PEM) approach in the present hybrid situation (which is ML estimation because of the Gaussian assumptions on e(t k )), a suitable error function ε(t k ), at the kth sampling instant, is given by,

ε(t k ) = D(q -1 , η) C(q -1 , η) y(t k ) - B(p, ρ) A(p, ρ) u(t k )
which can be written as,

ε(t k ) = D(q -1 , η) C(q -1 , η)  1 A(p, ρ) [A(p, ρ)y(t k ) -B(p, ρ)u(t k )] ff (16)
where the discrete-time prefilter D(q -1 , η)/C(q -1 , η) will be recognised as the inverse of the ARMA(p,q) noise model. Note that in these equations, we are mixing discrete and continuous-time operators somewhat informally in order to indicate the hydrid computational nature of the estimation problem being considered here. Thus, operations such as, B(p, ρ)

A(p, ρ) u(t k )
imply that the input variable u(t k ) is interpolated in some manner. This is to allow for the inter-sample behaviour that is not available from the sampled data and so has to be inferred in order to allow for the continuous-time numerical integration of the associated differential equations. Minimisation of a least squares criterion function in ε(t k ), measured at the sampling instants, provides the basis for stochastic estimation. However, since the polynomial operators commute in this linear case, ( 16) can be considered in the alternative form,

ε(t k ) = A(p, ρ)y f (t k ) -B(p, ρ)u f (t k ) (17)
where y f (t k ) and u f (t k ) represent the sampled outputs of the complete hybrid prefiltering operation involving the continuous-time filtering operations using the filter

f c (p, ρ) = 1 A(p, ρ) (18)
as well as discrete-time filtering operations, using the inverse noise model filter

f d (q -1 , η) = D(q -1 , η) C(q -1 , η) (19)
The associated, linear-in-the-parameters estimation model then takes the form

y (n) f (t k ) = ϕ T f (t k )ρ + η(t k ) (20)
where,

ϕ T f (t k ) = [-y (n-1) f (t k ) • • • -y f (t) u (m) f (t k ) • • • u f (t k )] (21
) and η(t k ) is the continuous-time noise signal η(t) = A(p, ρ)ξ(t) sampled at the kth sampling instant.

2) RIVC Estimation: Optimal methods of IV estimation (see, e.g., [START_REF] Young | Some observations on instrumental variable methods of time-series analysis[END_REF], [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]) normally involve an iterative (or relaxation) algorithm in which, at each iteration, the 'auxiliary model' used to generate the instrumental variables, as well as the associated prefilters, are updated, based on the parameter estimates obtained at the previous iteration. Let us consider, therefore, the jth iteration where we have access to the estimate,

θj-1 = ρj-1 ηj-1 (22)
obtained previously at iteration j -1. The most important aspect of optimal IV estimation is the definition of an optimal instrumental variable. In the present context, this is generated from the output of the continuous-time auxiliary model,

x(t, ρj-1 ) = G(p, ρj-1 )u(t) (23)
which is prefiltered in the same hybrid manner as the other variables. The associated optimal IV vector φf (t k ), is then an estimate of the noise-free version of the vector ϕ f (t k ) in [START_REF] Rao | Identification of continuous-time systems: direct or indirect?[END_REF] and is defined as follows

φf (t k ) = -x (n-1) f (t k ) • • • -xf (t k ) u (m) f (t k ) • • • u f (t k ) T (24)
where it should be noted that

φf (t k ) = φf (t k , ρj-1 , ηj ) (25)
because the instrumental variables are now prefiltered and so are a function of both the system parameter estimates at the previous iteration and the most recent noise model parameter estimates. For simplicity, however, these additional arguments will be omitted in the subsequent analysis. Note also that the noise-free version of the vector ϕ f (t k ) in ( 21), which we will define as follows,

φT f (t k ) = -x (n-1) f (t k ) • • • -x f (t k ) u (m) f (t k ) • • • u f (t k ) (26) where x(t) = G o (p)u(t)
, is referred to in Section II-D when considering the statistical properties of the optimal IV parameter estimates.

The IV optimisation problem can now be stated in the form

ρj (N ) = arg min ρ 1 N N k=1 φf (t k )ϕ T f (t k ) ρ - 1 N N k=1 φf (t k )y (n) f (t k ) 2 Q ( 27
)
where x 2 = x T Qx and Q = I. This results in the wellknown solution of the IV estimation (IV normal) equations

ρj (N ) = N k=1 φf (t k )ϕ T f (t k ) -1 N k=1 φf (t k )y (n) f (t k ) ( 28
)
where the ρj (N ) is the IV estimate of the system model parameter vector at the jth iteration based on the appropriately prefiltered input/output data Z N = {u(t k ); y(t k )} N k=1 . As regards the hybrid prefiltering, it will be noted from [START_REF] Solo | Time Series Recursions and Stochastic Approximation[END_REF] that this involves the inverse noise model parameters ηj obtained at the current jth iteration. This is because, given ρj-1 , an estimate of the sampled noise signal ξ(t k ), at the jth iteration, is obtained by subtracting the sampled output of the auxiliary model equation ( 23) from the measured output y(t k ), i.e.,

ξ(t k ) = y(t k ) -x(t, ρj-1 ) (29)
This estimate provides the basis for the estimation of the noise model parameter vector η j , using whatever ARMA model estimation algorithm is selected for this task.

C. The RIVC and SRIVC Algorithms

The iterative RIVC and SRIVC algorithms follow directly from the RIV and SRIV algorithms for DT systems (e.g., [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF]). This section summarises both algorithms.

1) The RIVC Algorithm: Bearing the analysis of the previous subsection II-B.2 in mind, the main steps in the RIVC algorithm are as follows:

Step 1. Initialisation: generate an initial estimate of the TF model parameter vector ρo using the simplified RIVC (SRIVC) algorithm (see subsection II-C.2) and use this to define the initial CT prefilter f c (p, ρo ).

Step 2. Iterative estimation.

for j = 1 : convergence (i) Generate the IV series x(t, ρj-1 ) using the auxiliary model built up from the estimated polynomials A(p, ρj-1 ) and B(p, ρj-1 ) based on ρj-1 at the previous (j -1)th iteration. (ii) Prefilter the input u(t k ), output y(t k ) and instrumental variable x(t, ρj-1 ) by the continuous-time filter fc(p, ρj-1 ) in order to generate the filtered derivatives of these variables. (iii) Obtain an optimal estimate of the noise model parameter vector ηj based on the estimated noise sequence ξ(t k ) from (29), using a selected ARMA estimation algorithm. (iv) Sample the filtered derivative signals at the discretetime sampling interval Ts and prefilter these by the discrete-time filter f d (q -1 , ηj ), in order to define all the required elements in the data vector ϕ f (t k ), the IV vector φf (t k ) and the nth-order filtered derivative y

(n) f (t k ). (v)
Based on these prefiltered data, generate the latest estimate ρj of the system model parameter vector using the en bloc IV solution [START_REF]System Identification, Environmetric Modelling and Control System Design[END_REF], or its recursive equivalent. Together with the estimate ηj of the noise model parameter estimate from (iii), this provides the estimate θj of the composite parameter vector at the jth iteration.

end

Step 3. After the convergence of the iterations is complete, compute the estimated parametric error covariance matrix Pρ , associated with the converged estimate ρ of the system model parameter vector, from the expression (see Section II-D),

Pρ = σ2 N k=1 φf (t k ) φT f (t k ) -1 (30)
where φf (t k ) is the IV vector obtained at convergence and σ2 is the estimated residual variance.

2) The SRIVC Algorithm: It will be noted that the above formulation of the RIVC estimation problem is considerably simplified if it is assumed that the additive noise is white, i.e., C(q -1 , η) = D(q -1 , η) = 1. In this case, simplified RIVC (SRIVC) estimation involves only the parameters in the A(p, ρ) and B(p, ρ) polynomials and the prefiltering only involves the continuous-time prefilter f c (p, ρ) = 1/A(p, ρ). Consequently, the main steps in the SRIVC algorithm are the same as those in the RIVC algorithm, except that the noise model estimation and subsequent discrete-time prefiltering in steps (ii) and (iii) of the iterative procedure are no longer required and are omitted.

It is worth noting that the RIVC algorithm has a much longer computation time than the SRIVC algorithm. As a result, it is advantageous to use the SRIVC algorithm for initial model order identification and only employ the full RIVC algorithm in those situations where the theoretical assumptions are satisfied and it is essential to have the most efficient parameter estimates and better estimates of the uncertainty on the parameters. For day-to-day usage, the SRIVC algorithm provides a quick and reliable approach to continuous-time model identification and estimation.

D. Theoretical Background and Statistical Properties of the RIVC Estimates

The motivational arguments presented in Section II-B suggest that, upon convergence, the RIVC parameter estimates will possess the optimal statistical properties of consistency and asymptotic efficiency when the additive noise has a Gaussian normal probability distribution and rational spectral density. This section presents more formal analysis to verify further the optimality of the estimates and confirm the asymptotic independence of the system and noise model parameter estimates.

1) Optimality of RIVC Estimation: In the control and systems literature, optimal IV estimation is usually considered in relation to the so-called 'extended IV' approach to estimation, as developed for the DT case [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF]. A similar approach can be applied in the present CT case by re-writing the IV optimisation equation [START_REF] Van Den Hof | Closed-loop issues in system identification[END_REF] in the following alternative form that explicitly reveals a continuous-time prefilter f (p)

ρ(N ) = arg min ρ [ 1 N N X k=1 ζ f (t k )f (p)ϕ T (t k )ρ - " 1 N N X k=1 ζ f (t k )f (p)y (n) (t k ) # 2 Q ( 31
)
where f (p) is the stable prefilter,

ζ f (t k ) is the prefiltered instrumental vector ζ f (t k ) = f (p)ζ(t k )
and Q is a positivedefinite matrix. By definition, when G o ∈ G, the extended IV estimate provides a consistent estimate under the following two conditions

Ē{ζ f (t k )f (p)ϕ T (t k )} is non-singular, Ē{ζ f (t k )f (p)ξ(t k )} = 0 (32)
Clearly, the selection of the instrumental variable vector ζ f (t k ), the weighting matrix Q and the prefilter f (p) may have a considerable effect on the covariance matrix P θ produced by the IV estimation algorithm.

In the open-loop situation, the Cramér-Rao lower bound on P θ for any unbiased identification method (e.g., [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Ljung | System Identification. Theory for the User[END_REF]) defines the optimal solution. In this regard, it has been shown that the minimum value of the covariance matrix P θ , as a function of the design variables ζ f (t k ), f (p) and Q, exists and is given by

P θ ≥ P opt θ with P opt θ = [ Ē{ ζf (t k ) ζT f (t k )}] -1 (33)
where ζf (t k ) is the optimally prefiltered IV vector, with the associated design variables defined as

Q = I, (34a)
f (p) = 1 H o (p)A o (p) = D o (p) C o (p)A o (p) , (34b)
ζ(t k ) = -x (n-1) (t k ) • • • -x(t k ) u (m) (t k ) • • • u(t k ) T (34c)
so that,

ζf (t k ) = f (p) ζ(t k ) (35)
which will be recognised as the noise-free, prefiltered vector φT f (t k ) defined earlier in [START_REF] Thil | Third-order cumulants based methods for continuous-time errors-in-variables model identification[END_REF].

2) Comments:

• Not surprisingly, the above analysis justifies the RIVC algorithmic design that iteratively updates those aspects of the theoretical solution that are not known a priori: in this case, the unknown model polynomials and the noise-free output of the system that is, of course, the source of the instrumental variables. If it is assumed that, in all identifiable situations, the RIVC algorithm converges in the sense that ρ ⇒ ρ and η ⇒ η, then the RIVC estimates will be consistent and asymptotically efficient.

• The optimal filter f (p) in (34b) is formulated in CT terms. In the proposed RIVC algorithm, this filter takes a hybrid form, as discussed in the previous sections.

One very important aspect of TF modelling is the identification of the model structure: i.e., the degrees n, m, p, and q of the model polynomials and any associated pure time delay τ . A model order selection method associated to the SRIVC model estimation method allows the user to automatically search over a whole range of different model orders. Two statistical measures are then used to help to user choose the best model structure (see Subsection II.E and [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models. In Identification of continuous-time models from sampled data[END_REF].

• Both RIVC/SRIVC routines are avalaible in the CON-TSID (see Section IV below) and CAPTAIN 1 Toolboxes for MATLAB .

E. Model Order Identification

One very important aspect of TF modelling is the identification of the model structure: i.e., the degrees n, m, p, and q of the model polynomials and any associated pure time delay τ . One statistical measure that is useful in this regard is the coefficient of determination R 2 T , defined as follows

R 2 T = 1 - σ 2 ξ σ 2 y ( 36
)
where σ 2 ξ is the variance of the estimated noise ξ(t k ) and σ 2 y is the variance of the measured output y(t k ). R 2 T is clearly a normalised measure of how much of the output variance is explained by the deterministic system part of the estimated model. However, it is well known that this measure, on its own, is not sufficient to avoid over-parametrisation and identify a parsimonious model, so that other model order identification statistics are required. In this regard, because the SRIVC and RIVC methods exploit optimal instrumental variable methodology, they are able to utilise the special properties of the instrumental product matrix (IPM) [START_REF] Young | An instrumental variable method for model order identification[END_REF]; in particular, the YIC statistic [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF] which is defined as follows

YIC = log e σ2 σ 2 y + log e {NEVN}; NEVN = 1 n θ n θ i=1 pii θ2 i (37)
Here, n θ = n + m + p + q + 1 is the number of estimated parameters; pii is the ith diagonal element of the blockdiagonal covariance matrix P θ , where,

P θ = P ρ 0 0 P η ( 38
)
and so is an estimate of the variance of the estimated uncertainty on the ith parameter estimate. θ2 i is the square of the ith parameter estimate in the θ vector, so that ratio pii / θ2 i is a normalised measure of the uncertainty on the ith parameter estimate.

From the definition of R 2 T , we see that the first term in the YIC is simply a relative measure of how well the model explains the data: the smaller the model residuals the more negative the term becomes. The normalised error variance norm (NEVN) term, on the other hand, provides a measure of the conditioning of the IPM, which needs to be inverted when the IV normal equations are solved (see, e.g., [START_REF] Young | Recursive Estimation and Time-Series Analysis[END_REF]): if the model is overparameterised, then it can be shown that the IPM will tend to singularity and, because of its illconditioning, the elements of its inverse (in the form here of the covariance matrix P θ ) will increase in value, often by several orders of magnitude. When this happens, the second term in the YIC tends to dominate the criterion function, indicating over-parametrisation.

1 http://www.es.lancs.ac.uk/cres/captain/ It is important to note that, based on practical experience, the YIC is normally best considered during SRIVC identification, which is much less computationally intensive than RIVC identification, so allowing for much faster investigation of the model order range selected by the user. In this situation, n θ is replaced by n ρ = n + m + 1 the pii are obtained by reference to the covariance matrix P ρ .

Although heuristic, the YIC has proven very useful in practical identification terms. It should not, however, be used as a sole arbiter of model order: rather the combination of R 2

T and YIC provides an indication of the best parsimonious models that can be evaluated by other standard statistical measures (e.g., the auto and partial autocorrelation of the model residuals, the cross-correlation of the residuals with the input signal u(t k ), etc.). Also, the physical interpretation of the model can often provide valuable information on the model adequacy: for instance, a model with complex eigenvalues caused by overparametrisation may prove incompatible with the non-oscillatory nature of the physical system under study.

III. LATEST DEVELOPMENTS FOR THE RIVC METHOD

Recent developments aimed at extending the RIVC method to handle wider practical situations in order to enhance the application field of direct CT model identification.

A. Multiple-input Systems

It is clearly straightforward to extend the RIVC/SRIVC methods to the multiple-input situation if the TF denominator is common to all input channels. The situation is not so straightforward in the case where there are different denominator polynomials for each input channel. However, following the RIV approach for DT systems [START_REF] Jakeman | Instrumental variable algorithms for multiple input systems described by multiple transfer functions[END_REF], the algorithms can be extended to handle this situation [START_REF] Garnier | An optimal IV technique for identifying continuous-time transfer function model of multiple input systems[END_REF]: indeed, the current version of RIVC in the CONTSID Toolbox provides this option.

B. Non-uniformly Sampled Data

One advantage of the SRIVC approach to continuoustime modelling is that it can be based on irregularly sampled data and can handle 'fractional' pure time delays. The current implementation of the SRIVC algorithm in the CONTSID Toolbox can handle irregularly sampled data. However, the RIVC algorithm has not yet been upgraded in this regard because it requires additional interpolation and re-sampling in order to generate a regularly sampled series for the ARMA noise model estimation parts of the algorithm.

C. Closed-loop Model Identification

Provided there is an external command input signal, the identification and estimation of a system within a closed automatic control loop has always been straightforward when using IV estimation methodology [START_REF] Young | An instrumental variable method for real-time identification of a noisy process[END_REF], [START_REF] Gilson | Instrumental variable methods for continuous-time closed-loop model identification[END_REF]. In the case of the RIVC/SRIVC algorithms, a two-stage approach, such as that used by Van den Hof [START_REF] Van Den Hof | Closed-loop issues in system identification[END_REF] for discrete-time systems, is the most effective, since it does not require prior knowledge of the control system. Recent research [START_REF] Young | Simple refined IV methods of closed-loop system identification[END_REF] has shown that a modification of this approach employing the SRIVC algorithm (rather than the FIR model estimation used by Van den Hof) for estimating the control input signal, followed by full RIVC estimation of the system, based on this estimated control input, works extremely well.

D. Hammerstein and LPV Model Identification

Direct identification of CT nonlinear models is still an unmature subject. This section discusses briefly the extension of the RIVC method for the identification of Hammerstein and LPV CT Box-Jenkins models. In the case of Hammerstein hybrid BJ model, the nonlinear function f (.) is assumed to be a sum of known basis functions γ 1 , γ 2 , . . . , γ l given as:

ū(t) = l i=1 α i γ i (u(t)) with α 1 = 1. (39)
The hybrid CT BJ Hammerstein model is described by the following input-output relationship:

     x(t) = G(p)ū(t) ξ(t k ) = H(q -1 )e(t k ), y(t k ) = x(t k ) + ξ(t k ), (40)
where

G(p) = B(p) A(p) . (41)
where the coloured noise associated with the sampled output measurement y(t k ) has rational spectral density and can be represented by a discrete-time autoregressive moving average ARMA model:

ξ(t k ) = H(q -1 )e(t k ) = C(q -1 ) D(q -1 ) e(t k ) (42)
The RIVC method has very recently been extended to estimate the parameters of such CT hybrid BJ Hammerstein models [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF], [START_REF] Laurain | Refined instrumental variable methods for Hammerstein Box-Jenkins models[END_REF]. So called Linear Parameter Varying (LPV) models have been the subject of recent interest. The RIVC approach has recently been extended to estimate CT LPV input/output models [START_REF] Laurain | Direct identification of continuous-time LPV input/ouput models[END_REF].

IV. SOFTWARE ASPECTS -THE CONTSID TOOLBOX

The field of system identification is an extensive and versatile area. It is easy to get confused by the vast number of approaches and variants of methods available. We have seen so far that direct continuous-time model identification from sampled data is now a mature subject and it is important to package the identification tools in a user-friendly way. An attempt to do that was carried out with the CONTinuous-time System IDentification (CONTSID) toolbox for MATLAB . The CONTSID toolbox was first released in 1999 [START_REF] Garnier | CONTSID: a continuous-time system identification toolbox for Matlab[END_REF]. It has gone through several updates. The key features of the CONTSID toolbox are [START_REF] Garnier | CONTSID toolbox: a software support for continuous-time data-based modelling[END_REF]:

• it supports most of the time-domain methods developed over the last thirty years [START_REF] Garnier | Continuous-time model identification from sampled data. Implementation issues and performance evaluation[END_REF] for identifying linear dynamic continuous-time parametric models from measured input/output sampled data;

• it provides transfer function and state-space model identification methods for single-input single-output (SISO) and multiple-input multiple-output (MIMO) systems, including both traditional and more recent approaches;

• it can handle mild irregularly sampled data in a straightforward way;

• it may be seen as an add-on to the system identification (SID) toolbox for MATLAB . To facilitate its use, it has been given a similar setup to the SID toolbox;

• it provides a flexible graphical user interface (GUI) that lets the user analyse the experimental data, identify and evaluate models in an easy way;

• It can be freely downloaded from http://www.cran.uhp-nancy.fr/contsid/

The latest version of the CONTSID toolbox has the following three major additions:

• it supports errors-in-variables CT transfer function model identification [START_REF] Mahata | Identification of continuous-time errorsin-variables models[END_REF], [START_REF] Thil | Third-order cumulants based methods for continuous-time errors-in-variables model identification[END_REF];

• it provides routines to estimate linear CT transfer function model in closed loop [START_REF] Gilson | Instrumental variable methods for continuous-time closed-loop model identification[END_REF], [START_REF] Young | Simple refined IV methods of closed-loop system identification[END_REF];

• it includes methods to identify nonlinear CT Hammerstein models [START_REF] Laurain | Refined instrumental variable methods for identification of Hammerstein continuous-time Box-Jenkins models[END_REF], [START_REF] Laurain | Refined instrumental variable methods for Hammerstein Box-Jenkins models[END_REF].

V. CONCLUSIONS This paper has first described the full RIVC algorithm for identifying hybrid Box-Jenkins transfer function models for linear, continuous-time systems from discrete-time, sampled data. The latest developments of the RIVC approach for nonuniformly sampled data, closed-loop identification as well as for nonlinear Hammerstein and LPV model identification have also been briefly discussed.

It is felt that continuous-time model identification, based on a stochastic formulation of the transfer function estimation problem, provides a theoretically elegant and practically useful approach to the modelling of stochastic dynamic systems from sampled data.

It is an approach that has many advantages in scientific terms since it provides differential equation models that conform with models used in most scientific research, where conservation equations are normally formulated in terms of differential equations. It is also a model defined by a unique set of parameter values that are not dependent on the sampling interval, so eliminating the need for conversion from discrete to continuous time that is an essential element of indirect approaches to estimation based on discrete-time model estimation. These direct continuous-time model identification methods have proven to be particularly well suited in the case of mild non-uniformly sampled data, dominant system modes with widely different natural frequencies (stiff systems), fast sampled data, or when the input does not respect the zeroorder hold assumption. Finally but not the least, these direct data-based CT modelling methods have proven successful in many practical applications and are available as user-friendly

and computationally efficient algorithms in the CONTSID toolbox for Matlab TM .

VI. ACKNOWLEDGMENTS

I am extremely grateful to Professor Peter Young for teaching me the basics of the Refined Instrumental Variable estimation concept. I would like also gratefully acknowledge the contributions of all of my students in the developments presented in the paper: Michel Mensler, Marion Gilson, Eric Huselstein, Stephane Thil, Damien Kuss and Vincent Laurain.