Kahina Ghidouche
email: kahina.ghidouche@univ-bejaia.dz

Abderrahmane Sider
email: ar.sider@univ-bejaia.dz

Lilia Ziane Khodja
email: zianekhodja.lilia@gmail.com

Raphael Couturier
email: raphael.couturier@univ-fcomte.fr

Raphaël Couturier

Two parallel implementations of

Keywords: finding method, Ehrlich-Aberth method, GPU, MPI, OpenMP

come

I. INTRODUCTION

Finding the roots of polynomials of very high degrees arises in many complex problems of various domains such as algebra, biology or physics. A polynomial p(x) in C in one variable x is an algebraic expression in x of the form:

p(x) = n i=0 α i x i , α n = 0, (1)
where {α i } 0≤i≤n are complex coefficients and n is a high integer number. If α n = 0 then n is called the degree of the polynomial. The root-finding problem consists in finding the n different values of the unknown variable x for which p(x) = 0. Such values are called roots of p(x). Let {z i } 1≤i≤n be the roots of polynomial p(x), then p(x) can be written as :

p(x) = α n n i=1 (x -z i ), α n = 0. (2)
Most of the numerical methods that deal with the polynomial root-finding problems are simultaneous methods, i.e. the iterative methods to find simultaneous approximations of the n polynomial roots. These methods start from the initial approximation of all n polynomial roots and give a sequence of approximations that converge to the roots of the polynomial. Two examples of well-known simultaneous methods for root-finding problem of polynomials are the Durand-Kerner method [START_REF] Durand | Solutions numériques des équations algébriques. Tome I: Équations du type F (x) = 0; racines d'un polynôme[END_REF], [START_REF] Kerner | Ein Gesamtschrittverfahren zur Berechnung der Nullstellen von Polynomen. (German) [A complete step method for the computation of zeros of polynomials[END_REF] and the Ehrlich-Aberth method [START_REF] Ehrlich | A modified newton method for polynomials[END_REF], [START_REF] Aberth | Iteration methods for finding all zeros of a polynomial simultaneously[END_REF].

The convergence time of simultaneous methods drastically increases with the increasing of the polynomial's degree. The great challenge with simultaneous methods is to parallelize them and to improve their convergence. Many authors have proposed parallel simultaneous methods [START_REF] Freeman | Calculating polynomial zeros on a local memory parallel computer[END_REF], [START_REF] Loizou | Higher-order iteration functions for simultaneously approximating polynomial zeros[END_REF], [START_REF] Freeman | Asynchronous polynomial zero-finding algorithms[END_REF], [START_REF] Bini | Numerical computation of polynomial zeros by means of aberth's method[END_REF], [START_REF] Couturier | Extraction de racines dans des polynômes creux de degré élevé[END_REF], [START_REF] Couturier | Iterative algorithms on heterogeneous network computing: Parallel polynomial root extracting[END_REF], using several paradigms of parallelization (synchronous or asynchronous computations, mechanism of shared or distributed memory, etc). However, so far until now, only polynomials not exceeding degrees of less than 100,000 have been solved.

The recent advent of the Compute Unified Device Architecture (CUDA) [START_REF]CUDA C programming guide[END_REF], a programming model and a parallel computing architecture developed by NVIDIA, has revived parallel programming interest in this problem. Indeed, the computing power of GPUs (Graphics Processing Units) has exceeded that of traditional CPUs processors, which makes it very appealing to the research community to investigate new parallel implementations for a whole set of scientific problems in the reasonable hope to solve bigger instances of well known computationally demanding issues such as the one beforehand. However, CUDA provides an efficient massive data computing model which is suited to GPU architectures. Ghidouche et al. [START_REF] Ghidouche | Parallel implementation of the Durand-Kerner algorithm for polynomial root-finding on GPU[END_REF] proposed an implementation of the Durand-Kerner method on a single GPU Tesla 2070. Their main results showed that a parallel CUDA implementation is about 10 times faster than the sequential implementation on a single CPU Intel(R) Xeon(R) CPU E5620@2.40GHz for sparse polynomials of degree 48,000.

In this paper we propose the parallelization of the Ehrlich-Aberth (EA) method which has a much better cubic convergence rate than the quadratic rate of the Durand-Kerner method that has already been investigated in [START_REF] Ghidouche | Parallel implementation of the Durand-Kerner algorithm for polynomial root-finding on GPU[END_REF]. In the other hand, EA is suitable to be implemented in parallel computers according to the data-parallel paradigm. In this model, computing elements carry computations on the data they are assigned and communicate with other computing elements in order to get fresh data or to synchronize. Classi-cally, two parallel programming paradigms OpenMP and MPI are used to code such solutions. But in our case, computing elements are CUDA multi-GPU platforms. This architectural setting poses new programming challenges but offers also new opportunities to efficiently solve huge problems, otherwise considered intractable until recently. To the best of our knowledge, our CUDA-MPI and CUDA-OpenMP codes are the first implementations of EA method with multiple GPUs for finding roots of polynomials. Our major contributions include:

• The parallel implementation of EA algorithm on a multi-GPU platform with a shared memory using OpenMP API. It is based on threads created from the same system process, such that each thread is attached to one GPU.

In this case the communications between GPUs are done by OpenMP threads through shared memory.

• The parallel implementation of EA algorithm on a multi-GPU platform with a distributed memory using MPI API, such that each GPU is attached and managed by a MPI process. The GPUs exchange their data by messagepassing communications. This approach is more used on clusters to solve very complex problems that are too large for traditional supercomputers, which are very expensive to build and run. • Our method is efficient to compute the roots of sparse and full polynomials of degree up to 5 million. The paper is organized as follows. In Section II we present three different parallel programming models OpenMP, MPI and CUDA. In Section III we present the implementation of the Ehrlich-Aberth algorithm on a single GPU. In Section IV we present the parallel implementations of the Ehrlich-Aberth algorithm on multiple GPUs using the OpenMP and MPI approaches. In section V we present our experiments and discuss them. Finally, Section VI concludes this paper and gives some hints for future research directions in this topic.

II. PARALLEL PROGRAMMING MODELS

Our objective consists in implementing a root-finding algorithm of polynomials on multiple GPUs. To this end, it is essential to know how to manage the CUDA contexts of different GPUs. A direct method to control the various GPUs is to use as many threads or processes as GPU devices. We investigate two parallel paradigms: OpenMP and MPI. In this case, the GPU indices are defined according to the identifiers of the OpenMP threads or the ranks of the MPI processes. In this section we present the parallel programming models: OpenMP, MPI and CUDA.

A. OpenMP

OpenMP (Open Multi-processing) is an application programming interface for parallel programming [START_REF]OpenMP application program interface[END_REF]. It is a portable approach based on the multithreading designed for shared memory computers, where a master thread forks a number of slave threads which execute blocks of code in parallel. An OpenMP program alternates sequential regions and parallel regions of code, where the sequential regions are executed by the master thread and the parallel ones may be executed by multiple threads. During the execution of an OpenMP program the threads communicate their data (read and modified) in the shared memory. One advantage of OpenMP is the global view of the memory address space of an application. This allows a relatively fast development of parallel applications with easier maintenance. However, it is often difficult to get high rates of performances in large scaleapplications.

B. MPI

MPI (Message Passing Interface

) is a portable message passing style of the parallel programming designed specifically for distributed memory architectures [START_REF]Parallel Programming with MPI[END_REF]. In most MPI implementations, a computation contains a fixed set of processes created at the initialization of the program in such a way that one process is created per processor. The processes synchronize their computations and communicate by sending/receiving messages to/from other processes. In this case, the data are explicitly exchanged by message passing while the data exchanges are implicit in a multithread programming model like OpenMP and Pthreads. However in the MPI programming model, the processes may either execute different programs referred to as multiple program multiple data (MPMD) or every process executes the same program (SPMD). The MPI approach is one of the most used HPC programming model to solve large scale and complex applications.

C. CUDA

CUDA (Compute Unified Device Architecture) is a parallel computing architecture developed by NVIDIA [START_REF]CUDA C programming guide[END_REF] for GPUs. It provides a high level GPGPU-based programming model to program GPUs for general purpose computations. The GPU is viewed as an accelerator such that data-parallel operations of a CUDA program running on a CPU are off-loaded onto GPU and executed by this latter. The data-parallel operations executed by GPUs are called kernels. The same kernel is executed in parallel by a large number of threads organized in grids of thread blocks, such that each GPU multiprocessor executes one or more thread blocks in SIMD fashion (Single Instruction, Multiple Data) and in turn each core of the multiprocessor executes one or more threads within a block. Threads within a block can cooperate by sharing data through a fast shared memory and coordinate their execution through synchronization points. In contrast, within a grid of thread blocks, there is no synchronization at all between blocks. The GPU only works on data filled in the global memory and the final results of the kernel executions must be transferred out of the GPU. In the GPU, the global memory has lower bandwidth than the shared memory associated to each multiprocessor. Thus with CUDA programming, it is necessary to design carefully the arrangement of the thread blocks in order to ensure a low latency and a proper use of the shared memory. As for the global memory accesses, it should also be minimized.

III. THE EHRLICH-ABERTH ALGORITHM ON A GPU

A. The Ehrlich-Aberth method

The Ehrlich-Aberth method is a simultaneous method [START_REF] Aberth | Iteration methods for finding all zeros of a polynomial simultaneously[END_REF] using the following iteration

z k+1 i = z k i - p(z k i ) p (z k i ) 1 - p(z k i ) p (z k i ) j=n j=1,j =i 1 (z k i -z k j ) , i = 1, . . . , n (3)
This method contains 4 steps. The first step consists in the initializing the polynomial. The second step initializes the solution vector Z using the Guggenheimer method [START_REF] Guggenheimer | Initial approximations in Durand-Kerner's root finding method[END_REF] to ensure that initial roots are all distinct from each other. In step 3, the iterative function based on the Newton's method [START_REF] Newton | Tractatus de methodis serierum et fluxionum[END_REF] and Weiestrass operator [START_REF] Weierstrass | Neuer beweis des satzes, dass jede ganze rationale function einer veranderlichen dagestellt werden kann als ein product aus linearen functionen derselben veranderlichen[END_REF] is applied. In our case, the Ehrlich-Aberth is applied as in [START_REF] Ehrlich | A modified newton method for polynomials[END_REF]. Iterations of the EA method will converge to the roots of the considered polynomial. In order to stop the iterative function, a stop condition is applied, this is the 4th step. This condition checks that all the root modules are lower than a fixed value .

∀i ∈ [1, n], | z k i -z k-1 i z k i | < (4)

B. Improving Ehrlich-Aberth method

With high degree polynomials, the EA method suffers from floating point overflows due to the mantissa of floating points representations. This induces errors in the computation of p(z) when z is large.

In order to solve this problem, we propose to modify the iterative function by using the logarithm and the exponential of a complex and we propose a new version of the EA method. This method allows us to exceed the computation of the polynomials of degree 100,000 and to reach a degree up to more than 1,000,000. The reformulation of the iteration (3) of the EA method with exponential and logarithm operators is defined as follows, for i = 1, . . . , n:

z k+1 i = z k i -exp(ln(p(z k i ))-ln(p (z k i ))-ln(1-Q(z k i ))), (5)
where:

Q(z k i ) = exp(ln(p(z k i ))-ln(p (z k i ))+ln( n i =j 1 z k i -z k j )). ( 6
)
Using the logarithm and the exponential operators, we can replace any multiplications and divisions with additions and subtractions. Consequently, computations manipulate lower values in absolute values [START_REF] Rhofir | Perfectionnements de la méthode asynchrone de Durand-Kerner pour les polynômes complexes[END_REF]. In practice, the exponential and logarithm mode is used when a root is outside the circle unit represented by the radius R evaluated in C language with: R = exp(log(DBL_M AX)/(2 * n)); [START_REF] Freeman | Asynchronous polynomial zero-finding algorithms[END_REF] where DBL_MAX stands for the maximum representable double value and n is the degree of the polynomial.

C. The Ehrlich-Aberth parallel implementation on CUDA

The algorithm 1 shows sketch of the Ehrlich-Aberth method using CUDA. The first step consists in the initialization of the input data, for exemple the polynomial P, the derivative of P and the vector solution Z. Then, all data of the root finding problem are transfered from the CPU memory to the GPU global memory, because the GPUs only work on the data filled in their memories. Next, all the data-parallel arithmetic operations inside the main loop (while(...)) are executed as kernels by the GPU. The first kernel named Kernelsave in line 5 of Algorithm 1 consists in saving the vector of polynomial roots found at the previous time-step in GPU memory, in order to check the convergence of the roots after each iteration (line 7, Algorithm 1). Then the new roots with the new iterations are computed using the EA method with a Gauss-Seidel iteration mode in order to use the latest updated roots (line 6). This improves the convergence compared to the Jacobi method. This kernel is, in practice, very long since it performs all the operations with complex numbers with the normal mode of the EA method as in Eq. 3 but also with the logarithm-exponential one as in Eq. [START_REF] Loizou | Higher-order iteration functions for simultaneously approximating polynomial zeros[END_REF][START_REF] Freeman | Calculating polynomial zeros on a local memory parallel computer[END_REF]. The last kernel checks the convergence of the roots after each update of Z k , according to formula Eq. 4 line [START_REF] Freeman | Asynchronous polynomial zero-finding algorithms[END_REF]. We used the functions of the CUBLAS Library (CUDA Basic Linear Algebra Subroutines) to implement this kernel. The algorithm terminates its computations when all the roots have converged. Listing 1 shows the a simplified version of second kernel code (some parameters in the kernels have been simplified in order to increase the readability). As can be seen this kernel calls multiple kernels, all the kernels for complex numbers and kernels for the evaluation of a polynomial are not detailed.

Listing 1. Kernels to update the roots / / Normal v e r s i o n o f t h e E h r l i c h -A b e r t h method _ _ d e v i c e _ _ cuDoubleComplex F i r s t H _ E A ( i n t i , cuDoubleComplex * Z ) { cuDoubleComplex r e s u l t ; cuDoubleComplex C , F , Fp ; i n t j ; cuDoubleComplex sum ; cuDoubleComplex un ; GPUs. Both approaches allow us to solve sparse and full polynomials of very high degrees. Using 4 GPUs allows us to achieve a quasi-linear speedup.

B. Comparison between the CUDA-OpenMP and the CUDA-MPI approaches

In the previous section we saw that both approaches are very efficient to reduce the execution times to solve sparse and full polynomials. In this section we try to compare these two approaches. In this experiment three sparse polynomials and three full polynomials of degrees 200,000, 800,000 and 1,400,000 are investigated. Figures 5 and6 show the comparison between CUDA-OpenMP and CUDA-MPI algorithms of the EA method to solve sparse and full polynomials, respectively.

In Figure 5 there is one curve for CUDA-OpenMP and another one for CUDA-MPI for each polynomial investigated. We can see that the results are quite similar between OpenMP and MPI for the polynomial degree of 200K. For the degree of 800K, the MPI version is a little bit slower than the OpenMP version but for the degree of 1,4 million, there is a slight advantage for the MPI version. In Figure 6, we can see that when it comes to full polynomials, both approaches are almost equivalent.

C. Solving sparse and full polynomials of the same degree on multiple GPUs

In this experiment we compare the execution times of the EA algorithm according to the number of GPUs to solve sparse and full polynomials on multiple GPUs using OpenMP or MPI approaches. We chose three sparse and three full polynomials of degrees 200,000, 800,000 and 1,400,000. Figures 7 and8 show the execution times to solve sparse and full polynomials of the same degrees with the CUDA-OpenMP version and the CUDA-MPI version, respectively. those to solve full polynomials. With sparse polynomials the number of monomials is reduced, consequently the number of operations is reduced and the execution time decreases. Figure 8 shows the impact of sparsity on the efficiency of the CUDA-MPI approach. We can see that the impact follows the same pattern, a difference in execution times in favor of the sparse polynomials.

D. Scalability of the EA method on multiple GPUs to solve very high degree polynomials

These experiments report the execution times of the EA method for sparse and full polynomials of high degrees ranging from 1,000,000 to 5,000,000. In Figure 9 we can see that both approaches (CUDA-OpenMP and CUDA-MPI) are scalable and can solve very high degree polynomials. In addition, with full polynomial as well as sparse ones, both approaches give very similar results. VI. CONCLUSION In this paper, we have presented parallel implementations of the Ehrlich-Aberth algorithm to solve full and sparse polynomials, on a single GPU with CUDA and on multiple GPUs using two parallel paradigms: shared memory with OpenMP and distributed memory with MPI. These architectures were addressed by a CUDA-OpenMP approach and CUDA-MPI approach, respectively. Experiments show that, using parallel programming model like OpenMP or MPI, we can efficiently manage multiple graphics cards to solve the same problem and accelerate the parallel execution with 4 GPUs and solve a polynomial of degree up-to 5,000,000 four times faster than on a single GPU.

Our next objective is to extend the model presented here to clusters of GPU nodes, with a three-level scheme: internode communications via MPI processes (distributed memory), management of multi-GPU nodes by OpenMP threads (shared memory). Actual platforms may probably also contain purely multi-core nodes without any GPU. This heterogeneous setting may lead to the integration of load balancing algorithms so as to allow an optimal use of hardware resources.

Algorithm 1 : 5 Z 6 Z

 156 Finding roots of polynomials with the Ehrlich-Aberth method on a GPU Input: (tolerance threshold) Output: Z (solution vector of roots) 1 Initialize the polynomial P and its derivative P ; 2 Set the initial values of vector Z; 3 Copy P , P and Z from CPU to GPU; 4 while error > do prev = KernelSave(Z); = KernelUpdate(P, P , Z); 7 error = KernelComputeError(Z, Z prev ); 8 Copy Z from GPU to CPU;

Fig. 1 .Fig. 2 .

 12 Fig. 1. Execution times in seconds of the Ehrlich-Aberth method to solve sparse polynomials on multiple GPUs with CUDA-OpenMP.

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 3456 Fig. 3. Execution times in seconds of the Ehrlich-Aberth method to solve sparse polynomials on multiple GPUs with CUDA-MPI.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Execution times to solve sparse and full polynomials of three distinct degrees on multiple GPUs using OpenMP. In Figure 7 the execution times of the CUDA-OpenMP version to solve sparse polynomials are very low compared to

Fig. 9 .

 9 Fig.9. Execution times in seconds of the Ehrlich-Aberth method to solve sparse and full polynomials of high degree on 4 GPUs for degrees ranging from 1M to 5M

ACKNOWLEDGMENT

This paper is partially funded by the Labex ACTION program (contract ANR-11-LABX-01-01) and the Franche-Comte regional council. Computations have been performed on the supercomputer facilities of the Mésocentre de calcul de Franche-Comté. We also would like to thank Nvidia for hardware donation under CUDA Research Center 2016.

/ / e v a l u a t e t h e p o l y n o m i a l F = F o n c t i o n ( Z [ i ] ) ; / / e v a l u a t e t h e d e r i v a t i v e o f t h e p o l y . The development of this code is a rather long task due to the development of all the kernels that compute the parts ported on the GPU. This comes in particular from the fact that it is very difficult to debug CUDA running threads like threads on a CPU host. In the following section the GPU parallel implementation of the Ehrlich-Aberth method with OpenMP and MPI is presented.

IV. THE EHRLICH-ABERTH ALGORITHM ON MULTIPLE GPUS

In order to manage the CUDA contexts of different GPUs, two parallel paradigms are investigated: OpenMP and MPI. In this section we present the OpenMP-CUDA and the MPI-CUDA approaches used to implement the Ehrlich-Aberth algorithm on multiple GPUs.

A. An OpenMP-CUDA approach Our OpenMP-CUDA implementation of EA algorithm is based on the hybrid OpenMP and CUDA programming model. This algorithm is presented in Algorithm 2. All the data are shared with OpenMP among all the OpenMP threads. The shared data are the solution vector Z, the polynomial to solve P , its derivative P , and the error vector error. The number of OpenMP threads is equal to the number of GPUs, each OpenMP thread binds to one GPU, and it controls a part of the shared memory. More precisely each OpenMP thread will be responsible for updating its own part of the vector Z. This part is called Z loc in the following. Then all GPUs will have a grid of computation organized according to the device performance and the size of data on which it runs the computation kernels.

To compute one iteration of the EA method each GPU performs the followings steps. First, roots are shared with OpenMP and the computation of the local size for each GPU is performed (line 4). Each thread starts by copying all the previous roots inside its GPU (line 5). At each iteration, the following operations are performed. First the vector Z is transferred from the CPU to the GPU (line 7). Each GPU copies the previous roots (line 8) and it computes an iteration of the EA method on its own roots (line 9). For that all the other roots are used. The local error is computed on the new roots (line 10) and the maximum of the local errors is computed on all OpenMP threads (line 11). At the end of an iteration, the updated roots are copied from the GPU to the CPU (line 12) and each CPU directly updates its own roots in the shared memory arrays containing all the roots. Algorithm 2: Finding roots of polynomials with the Ehrlich-Aberth method on multiple GPUs using OpenMP Input: (tolerance threshold) Output: Z (solution vector of roots) Copy Z loc from GPU to Z in CPU;

B. A MPI-CUDA approach

Our parallel implementation of EA to find the roots of polynomials using a CUDA-MPI approach follows a similar approach to the one used in CUDA-OpenMP. Each processor is responsible for computing its own part of roots using all the roots computed by other processors at the previous iteration. The difference between both approaches lies in the way processors communicate and exchange data. With MPI, processors need to send and receive data explicitly. So in Algorithm 3, after the initialization phase all the processors have the same Z vector. Then they need to compute the parameters used by the M P I_AlltoAll routines (line 4). In practice, each processor needs to compute its offset and its local size. Processors need to allocate memory on their GPU and need to copy their data on the GPU (line 5). At the beginning of each iteration, a processor starts by transferring the whole vector Z from the CPU to the GPU (line 7). Only the local part of Z prev is saved (line 8). After that, a processor is able to compute an updated version of its own roots (line 9) with the EA method. The local error is computed (line 10) and the global error is also computed using M P I_Reduce (line 11). Then the local roots are transferred from the GPU memory to the CPU memory (line 12) before being exchanged between all processors (line 13) in order to give to all processors the last version of the roots (with the M P I_AlltoAll routine). If the convergence is not satisfied, a new iteration is executed.

V. EXPERIMENTS

We study two categories of polynomials: sparse polynomials and full polynomials. A sparse polynomial is a polynomial for which only some coefficients are not null. In this paper, we consider sparse polynomials for which the roots are distributed on 2 distinct circles:

Algorithm 3: Finding roots of polynomials with the Ehrlich-Aberth method on multiple GPUs using MPI Input: (tolerance threshold) Output: Z (solution vector of roots) Copy Z loc from GPU to CPU;

13 Z = M P I_AlltoAll(Z loc );

A full polynomial is, in contrast, a polynomial for which all the coefficients are not null. A full polynomial is defined by:

In the following experiments, all the reported results have been obtained on a machine equipped with a CPU Intel(R) Xeon(R) CPU X5620@2.40GHz node with 64GB of ram and 4 Tesla Kepler K40 GPUs (2880 cores) with CUDA 7.5, OpenMP and the OpenMPI 1.7.5.

In order to evaluate both the GPU and Multi-GPU approaches, we performed a set of experiments on a single GPU and multiple GPUs using OpenMP or MPI with the EA algorithm, for both sparse and full polynomials of different degrees. All experimental results obtained are performed with double precision floating-point data and the convergence threshold of the EA method is set to 10 -7 . The initialization values of the vector solution of the methods are given by the Guggenheimer method [START_REF] Guggenheimer | Initial approximations in Durand-Kerner's root finding method[END_REF].

A. Evaluation of the multi-GPUs approaches

In this part, we evaluate the performances of the CUDA-OpenMP and CUDA-MPI approaches of the EA algorithm on different GPU platforms composed each of 1, 2, 3 or 4 GPUs. In this experiments we report the experimental results of the EA algorithms to find the roots of different sparse and full polynomials of high degrees ranging from 100,000 to 1,400,000. Figures 1 and2 show the execution times to solve, respectively, sparse and full polynomials with the CUDA-OpenMP algorithm, and Figures 3 and4 show those to solve, respectively, sparse and full polynomials with the CUDA-MPI algorithm.

All these figures show that the CUDA-OpenMP and the CUDA-MPI approaches of the EA algorithm, compared to the single GPU version, are efficient and scale well with multiple


Carol Habib
email: carol.habib@univ-fcomte.fr

Abdallah Makhoul
email: abdallah.makhoul@univ-fcomte.fr

Rony Darazi
email: rony.darazi@ua.edu.lb

Raphael Couturier
email: raphael.couturier@univ-fcomte.fr

Raphaël Couturier

Real-time Sampling Rate Adaptation based on Continuous Risk Level Evaluation in Wireless Body Sensor Networks

Keywords: Adaptive sampling, risk level evaluation, behavior function, health assessment, wireless body sensor networks

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

In the past decade, Wireless Body Sensor Networks (WB-SNs) emerged as a low cost solution for remote and continuous patient monitoring, thus reducing unnecessary hospitalization as well as keeping doctors and the healthcare center up-to-date with health-related matters [START_REF] Carmen | Body sensor networks: In the era of big data and beyond[END_REF]. Many populations can benefit from such a technology: elderly in nursing homes, patients at home after a surgical intervention, patients susceptible to cancer or heart disease, athletes etc. Some of the most popular applications in the literature include fall detection, stress evaluation systems, gait analysis, monitoring of vital signs and physilogical signals [START_REF] Miguel Pires | From data acquisition to data fusion: a comprehensive review and a roadmap for the identification of activities of daily living using mobile devices[END_REF], [START_REF] Fortino | A framework for collaborative computing and multi-sensor data fusion in body sensor networks[END_REF], [START_REF] Begum | Physiological sensor signals classification for healthcare using sensor data fusion and case-based reasoning[END_REF]. A WBSN is composed of sensor nodes and a coordinator. The former are placed on the body of the patient. They sense physiological signals, extract vital signs and wirelessly send these measurements to the coordinator for fusion in a periodic fashion. We suppose that each sensor node monitors only one vital sign (e.g. temperature, respiration rate, etc.). The coordinator usually is the smartphone, PDA or any other portable device which is carried by the patient. Many challenges exist in such a network such as: the heterogeneity of the collected data (heart rate, respiration rate, blood pressure etc.) and its huge amount, the energy consumption due to periodic transmission as well as privacy and security issues [START_REF] Carmen | Body sensor networks: In the era of big data and beyond[END_REF], [START_REF] Chen | Body area networks: A survey[END_REF], [START_REF] Salim | Adaptive sampling algorithms with local emergency detection for energy saving in wireless body sensor networks[END_REF]. In this paper we address the energy consumption and data reduction issues which are directly related to one another. In [START_REF] Rault | Energy efficiency in wireless sensor networks: A top-down survey[END_REF], the authors have classified energy-efficient mechanisms in Wireless Sensor Networks (WSNs) into five categories: data reduction approaches, sleep/wakeup schemes, radio optimization techniques, energyefficient routing methods and battery repletion. These approaches can be split up into software and hardware strategies [START_REF] Magno | Energy-efficient context aware power management with asynchronous protocol for body sensor network[END_REF]. In addition, some of these mechanisms are suitable for large scale networks such as environmental monitoring, industry, public safety or military systems applications. Thus, they cannot be applied in WBSNs where the network characteristics are different. For example, energy-efficient routing methods as well as transmission power control and topology control approaches cannot be directly used in WBSNs [START_REF] Chen | An adaptive sensor data segments selection method for wearable health care services[END_REF]. Many existing approaches in WBSNs have proposed energy-efficient data collection approaches. Some have used context-awareness based on activity recognition to perform adaptive sampling or adaptive sensing [START_REF] Magno | Energy-efficient context aware power management with asynchronous protocol for body sensor network[END_REF]. Some of them applied these approaches only on WBSNs composed of inertial detection sensor nodes such as accelerometers and gyroscopes [START_REF] Chen | An adaptive sensor data segments selection method for wearable health care services[END_REF]. Others suggested adapting the transmission of sensor nodes according to the relevance between symptoms and diseases [START_REF] Wang | Information-based energy efficient sensor selection in wireless body area networks[END_REF]. Others have used compressive sensing (CS) theory, which has emerged as an energy-efficient approach for wireless communication, due to the fact that biosignals such as the electrocardiogram (ECG) are sparse [START_REF] Mamaghanian | Compressed sensing for real-time energy-efficient ecg compression on wireless body sensor nodes[END_REF], [START_REF] Faust | Compressed sampling for heart rate monitoring[END_REF], [START_REF] Li | A continuous biomedical signal acquisition system based on compressed sensing in body sensor networks[END_REF], [START_REF] Wang | A configurable energy-efficient compressed sensing architecture with its application on body sensor networks[END_REF]. This technique reduces the amount of wirelessly transmitted data by sampling the signal of interest at a much lower rate than the traditional Shannon-Nyquist theorem. Thus, it has the potential of dramatically reducing the power consumption due to transmission. However, this technique has not been applied at the analog level since it requires extensive work on the analog sensor read-out electronics prior to the analog-to-digital converter (ADC). Its application on the digital level needs processing and energy resources and depends on the sensor's microcontroller [START_REF] Mamaghanian | Compressed sensing for real-time energy-efficient ecg compression on wireless body sensor nodes[END_REF]. In addition, it requires the sparsity of the signals of interest, thus in case of non-sparstity it does not perform as desired [START_REF] Faust | Compressed sampling for heart rate monitoring[END_REF]. In [START_REF] Jiang | A self-adaptively evolutionary screening approach for sepsis patient[END_REF], the authors have proposed AdapDBN a selfadaptively sepsis screening system allowing to shorten the time of syndrome detection. In a previous work ( [START_REF] Habib | Self-adaptive data collection and fusion for health monitoring based on body sensor networks[END_REF]), a distributed self-adaptive data collection approach has been proposed in the context of WBSNs. An adaptive sampling rate scheme, having a direct impact on the sensing, processing and especially the transmission tasks of the biosensor node, has been proposed. Using a Quadratic Bezier Curve as a BehaVior (BV) Function, it takes into account two parameters : the evolution of the monitored vital sign over time and its monitoring importance, based on a medical judgment, regarding the patient's health condition. These parameters are, respectively, determined by the Fisher Test and one-way Analysis of Variance (ANOVA) which study the variances of the sensed measurements over time and by a static risk level variable. However, the overall health condition of a patient, being continuously and remotely monitored on a long-term basis, changes over time. It is subject to many health events which can be acute or even chronic. Thus, it can vary from day to day as well as from an improvement state into a deterioration state and vice versa especially that acute disease go through many different stages. As a consequence, the monitoring importance given for each vital sign should be adapted with these changing conditions. This matter, has a direct influence on the collection of data, therefore on the energy consumption of the WBSN and the early detection of critical events. In this paper, we propose to dynamically adapt the risk level of a vital sign according to the changing health condition of the patient. The remainder of the paper is organized as follows. Section II presents our work's background. The risk level of the sensor node is defined and its adaptation is explained in Section III. Experimental results are shown and discussed in Section IV. Finally Section V concludes the paper with some directions and future work.

II. BACKGROUND

In this section, a brief overview about early warning score systems is given and the BehaVior (BV) function is extensively explained by adapting its equations based on the requirements of vital signs monitoring. The former is used to evaluate the patient's overall health condition as well as the severity of each vital sign, whereas, the latter is used for the sampling rate adaptation.

A. Early Warning Score Systems

An early warning score system (EWS) is a chart used by emergency medical services staff in hospitals to evaluate patients' heath status. The vital signs are weighed and aggregated in order to allow an early recognition of patients who are subject to an acute illness or those whose health condition is deteriorating [START_REF] Baskaran Ragesh | A survey on futuristic health care system: Wban[END_REF]. For each vital sign, a normal healthy range is defined. Values outside of this range are allocated a score according to the magnitude of the deviation from the normal range. The score weighing reflects the severity of the physiological disturbance. Such scoring systems can give the biosensor node the ability to locally evaluate the severity level of the vital sign being monitored and to assign to it a score. We have used the National EWS (NEWS) in our work in order to evaluate the severity level of each of the following vital signs: temperature, heart rate, respiration rate, oxygen saturation and systolic blood pressure. A score of 0 is assigned to the measurement if it is in the normal range and a score of 1, 2 or 3, according to its magnitude of deviation from the normal range, where 3 indicates the highest level of criticality. NEWS is standarized and employed in hospitals in the United Kingdom (UK) for the assessment of accute-illeness [START_REF]National early warning score[END_REF].

B. Behavior Function

Bezier curves are flexible parametric curves allowing the definition of shapes given the knowledge about some points of interest [START_REF] Makhoul | An adaptive scheme for data collection and aggregation in periodic sensor networks[END_REF]. The quadratic bezier curves, which curvature is defined following three points of interest, are used in the proposed sampling adaptation model. A quadratic bezier curve is limited by the two points P 0 (l x ; l y ) and P 2 (h x ; h y ) and its curvature is controlled by P 1 (b x ; b y ). It can be written under the following general form:

ax 2 + bxy + cy 2 + ux + vy + w = 0 (1)
where the values of the coefficients depend on P 0 , P 1 and P 2 .

In our approach, the x-axis represents the Fisher Test result F and the y-axis represents the sampling rate SR. Thus, the three points of the BV function are defined as follows: P 0 (0; l y ), P 1 (b x ; b y ), P 2 (h x ; h y ) where P 0 corresponds to F = 0 and SR = SR min , P 1 is the behavior point and P 2 corresponds to F = F t and SR = SR max . F t (h, N -h) corresponds to the critical value given by the Fisher Test for N collected samples during h periods. The BV function curve can be drawn using the quadratic bezier curve as follows:

BV (F ) =          hx+ly-2by 4b 2 x F 2 + hy-ly hx F + l y , h x = 2b x (h y + l y -2b y )α(F ) 2 + 2(b y -l y )α(F ) + l y , otherwise with α(F ) = -b x + √ b 2 x -2bxF +hxF hx-2bx
, such as 0 ≤ b x ≤ h x , 0 ≤ F ≤ h x and h x > 0. P 1 moves on the diagonal [AB] of the behavior rectangle, such as A(0; h y ) and B(h x ; l y ). Thus, the coordinates b x and b y of P 1 satisfy the equation of [AB] which is defined as follows:

y = l y -h y h x × x + h y (2)
In our approach, the position of P 1 should be determined according to the value of the risk level r. The coordinates b x and b y of P 1 are found using the following equation:

B r : [0; 1] -→ [0; h x ] × [l y ; h y ] r -→ (b x ; b y )
Thus, the following equations can be derived:

B r (r) = b x = (1 -r) × h x b y = l y + r × (h y -l y ) (3)
The closer the value of r is to 1 the more the vital sign is judged to be critical and its monitoring importance is increased. Figure 1 shows the BV function which is used to adapt the sampling rate according to the Fisher test result F and the risk level r of the vital sign given the patient's health condition. Having a maximum sampling rate SR max and a minimum sampling rate SR min , depending on the application's requirements, and the critical value F t given by the Fisher Test, the BV function is defined as follows:

BV (SR max , SR min , r, F, F t ) = SR (4)
The closer F is to F t , the more the collected measurements in h periods present variations. Thus, the higher the sampling rate is given to the sensor node in order not to miss any important variations. Whereas, the closer F is to 0, the less the collected measurements in h periods present variations. Thus, a lower sampling rate is given to the sensor node in order to preserve their energy level. A curve is associated for each risk level value. The higher the risk level value, the greater the sampling rate values. Indeed, for high risk levels the sampling rate of the sensor node should be increased in order not to miss any important changes and to keep track with the variations of the corresponding vital sign. As shown in Figure 2, for the same value of F , if the vital sign has a low risk level r, the sensor's measurement sampling rate will take a lower value than if the vital sign has a high risk level r. Thus, in the first case, the biosensor node will preserve its energy and will sense data and process measurements at a lower rate. However, in the latter case, the biosensor node, which is monitoring a vital sign being medically judged as essential regarding the patient's health condition, will be assigned a higher sampling rate in order not to miss important measurements and events.

SR

In the following section, we describe the risk level adaptation of a given vital sign according to the dynamic change of Fig. 2: Sampling rate adaptation according to the vital sign's risk level the patient's health condition over time. The sensor node will be able to adapt its sampling rate based on the BV curve that corresponds to its dynamically adapted risk level.

III. RISK LEVEL ADAPTATION

A. Scenario

Tom is an elderly person living in a nursing home. In order to keep track of his health condition, he is being remotely and continuously monitored by a WBSN. We suppose that the WBSN is composed of n biosensor nodes, where each node monitors one vital sign. At first, the monitoring importance for each vital sign is medically judged by the medical team. Accordingly, the risk levels r for all the vital sign are chosen. Since Tom's health condition can change with time, he can become sick or be subject to dangerous health events. Thus, the risk levels r should be adapted by the WBSN with Tom's changing health condition. Higher values should be given to the nodes in order to increase their sampling rate and capture important measurements and when the patient's health condition is at a lower risk and when it is normal, lower risk values should be given to the nodes in order to preserve their energy level. Next, we give a proper definition to the risk level r and discuss its initial setup.

B. Definition of the risk level r

The variable r represents the monitoring importance given to a specific vital sign. At first, healthcare experts chose the values for all the monitored vital signs. r is assigned a low value (< 0.5) if the following conditions are met:

• If the vital sign generally does not present many variations and is usually stable. • If the medical expert thinks that a low sampling rate is sufficient to capture any important variation. • If the variations do not have a drastic impact on the patient's health. For example, the temperature, the galvanic skin response and the oxygen saturation meet these requirements.

r is assigned a high value (≥ 0.5) if the following conditions are met:

• If the vital sign usually presents many variations.

• If the medical expert thinks that a high sampling rate is essential to capture any important variation. • If any variation can have a drastic impact on the patient's health. For example, the heart rate, the blood pressure and the respiration rate meet these requirements.

C. Risk Level Evaluation Function

We propose to adapt the risk level r of a vital sign, throughout the life of the WBSN, according to the patient's changing health condition. Let r global represent the patient's overall health condition where r global ∈ [0; 1] and r global ∈ R. r global is referred as the patient's global risk level and it indicates the severity level of his/her health condition. r global is evaluated by the coordinator based on the multi-sensor data fusion [START_REF] Habib | Multisensor data fusion for patient risk level determination and decisionsupport in wireless body sensor networks[END_REF]. Let S be the average score of the sensed measurements for a given vital sign during one round R, such as R = h×p where p is a time period and h ∈ N, and S max be the maximum score that a measurement can have according to the used EWS (cf. Section II-A). Then, the risk level evaluation function Eval : (r global , S) -→ [0; 1] is defined as follows:

Eval(r global , S) = α × r global + β × S S max (5)
where α and β are weight coefficients such as : 0 ≤ α ≤ 1 and 0 ≤ β ≤ 1 and α+β = 1. The Score ratio S Smax is equal to the mean score of the sensed measurements over the maximum score that can be given to a measurement. It represents the situation of a vital sign (its criticality) compared to the worst case (highest criticality level). S is calculated as follows:

S = 1 m × m i=1 s i ( 6
)
where m is the total number of sensed measurements and s i is the score of th i th measurement. Figure 3 shows the output of the proposed risk level evaluation function for two different parameter setups. The blue surface corresponds to a parameterization of α = β = 0.5, whereas the pink surface corresponds to a parameterization of α = 0.2 and β = 0.8. In the former case, both the score ratio and the overall health condition of the patient are given equal weights. Thus, both are equally important. In the latter case, the score ratio is given a greater importance than the overall health condition. For example, for both cases, the lower the global risk and the higher the score ratio, then the higher the risk level. However, in the latter case the risk level will have higher values than the former case for the same global risk value, since the impact of the score ratio is greater than the impact of the global risk. The α and β values are to be judged and parametrized by the healthcare experts, depending on whether the overall health condition of the patient or the status of the vital sign itself is more essential given the monitoring needs. For instance, some viruses such as flu or infections are accompanied by fever. Therefore, the temperature should be given a higher monitoring importance given its impact on the these types of sickeness and given that healthcare experts are interested in monitoring its variations over time. Thus, β better have a higher value than α in order to give a higher importance to the score ratio of the temperature rather than the overall health condition regrouping all the vital signs being monitored. for each round do 4:

for each period do end for 22: end while In this section, the algorithm for adapting the sampling rate in real-time is presented. The proposed algorithm is to be locally performed on each sensor node of the WBSN. The network is managed by one coordinator. Each sensor node locally adapts its sampling rate after adapting its risk level r using the proposed risk level evaluation function. On one side, for each round, the coordinator sends the value of r global , which is evaluated based on multisensor data fusion, to all the sensor nodes. On the other side, for each round, each sensor nodes computes the score ratio and evaluate its risk level based on these two values. Then, the sampling rate adaption can be made based on the result of Fisher test F and r using the BV function. The sampling rate adaptation is performed each round k such as Round k = h+k-1 i=k p i where k ∈ N,h ∈ N and p is the time period. In other words, it is performed at end of each period based on the measurements collected during the last h periods. The global risk r global , which is sent by the coordinator to the sensor nodes, is based on the data fusion of the collected measurements during these last h periods.

IV. EXPERIMENTAL RESULTS

Experiments are conducted using a cutom-based Java simulator. In order to evaluate the performance of the proposed approach, patient vital signs datasets are collected from Multiple Intelligent Monitoring in Intensive Care (MIMIC) I and II databases of PhysioNet [START_REF] Physionet | [END_REF]. We have tested our approach on different patient records as well as different vital signs such as the heart rate (HR), the respiration rate (RESP), the systolic blood pressure (ABPsys), the blood temperature (BLOODT) and the oxygen saturation (SpO2). We suppose that each vital sign is monitored by only one sensor node: if N vital signs are required to be monitored, then N sensors nodes are deployed on the patient's body. The proposed algorithm (cf. Algorithm 1) is implemented and tested at the level of each node. After running multiple experiments, noting that a minimum of 30 samples is required for the Fisher Test and ANOVA, the parameters settings are chosen as follows:

• Period p = 100 sec and Round R = 2 × p.

• Minimum sampling rate SR min = 1 samples/3.33 sec (corresponding to 30 samples per period) and Maximum sampling rate SR max = 1 sample/1.42 sec (corresponding to 70 samples per period). • Fisher Risk α f isher = 0.05.

• Initial risk level on all sensor nodes is randomnly set to r 0 = 0.4, indicating that all vital signs are slightly critical and have the same impact on the patient's health. The global risk level which represents the overall health condition of the patient is determined by the coordinator by running the algorithm Health-RAD [START_REF] Habib | Health risk assessment and decision-making for patient monitoring and decision-support using wireless body sensor networks[END_REF] for the multi-sensor data fusion and health assessment [START_REF] Habib | Multisensor data fusion for patient risk level determination and decisionsupport in wireless body sensor networks[END_REF].

First, we show the risk adaptation over time for different α and β setups and discuss their impact on the risk evaluation. Then we compare the results obtained in terms of data reduction, loss of information and energy consumption for 3 different scenarios : static risk level r = 0.4, static risk level r = 0.9 and our proposed approach.

A. Risk Adaptation over time

Figure 4 shows the ABPSys sensor node's risk level adaptation over 70 periods for 3 different parameter setups. In our proposed approach, α and β respectively denote the weights accorded for the patient's overall health condition and for the vital sign's severity level. Initially, the sensor node's risk level is set to 0.4, it is clear to see that in the 3 cases the risk evaluation function has yielded most of the time a risk value different than 0.4 and whose values range from low to high risk levels. Hence, pointing out the importance of dynamically adapting the risk level of each sensor node according to the changes in the patient's health conditon. When comparing Figures 4a and4b, we can clearly see the influence of both parameters on the risk level evaluation. For example, between period 45 and period 60, the risk level is mainly between 0.4 and 0.6 for α = 0.7 and β = 0.3, however it ranges between 0.6 and 0.8 for α = 0.3 and β = 0.7. This is due to the fact that between these periods the severity level of the ABPSys represented by the score ratio is more critical than the patient's overall health condition. Hence, giving the score ratio a higher impact than the overall health condition has yeilded a higher risk level (see Figure 4c) than in the opposite case (see Figure 4b). Whereas, between periods 8 and 18, the results show that the overall health condition is more critical than the severity level of the ABPSys. The risk level ranges between 0.2 and 0.4, in the case of α = 0.7 and β = 0.3 (see Figure 4b), however it ranges between 0.1 and 0.2 in the case of α = 0.3 and β = 0.7 (see Figure 4c). Hence, the parameters α and β of the risk level evaluation function should be medically judged by the healthcare experts on the basis of the importance given to each of the vital sign's severity level and the patient's overall health condition. In the rest of this section, the results correspond to a setup of α = β = 0.5.

B. Sampling Rate Adaptation and Energy Consumption

In this section, we compare the sampling rate adaptation and the energy consumption between two cases: static (constant) risk level r and our proposal of dynamic risk level r. We study the sampling rate adaptation and the energy consumption of the HR sensor node over 70 periods for 2 different patients. Patient 1 corresponds to record s01840-3454-10-24-18-46nm from MIMIC II database whereas patient 2 corresponds to record 276n from MIMIC database . According to the available datasets, all of the 5 vital signs are monitored for patient 1 and only the HR and the ABPSys are monitored for patient 2.

1) Sampling Rate Adaptation: Figure 5 shows the sampling rate adaptation of the HR sensor node for both patients for 3 different cases: static risk r=0.4, static risk r=0.9 and dynamic risk. For both patients, when comparing the two static cases (see Figures 5b,5c, 5e, 5f), we can notice that the sampling rate adaptation results in higher sampling rates when r = 0.9 than when r = 0.4. For example periods 2 and 10 in Figures 5b,5c and from period 50 till period 70 in Figures 5e and5f. This is due to the definition of the BV function which results in higher sampling rates for higher risk level values. However, if we compare the results obtained in the dynamic risk level case with the results of the two static cases for patient 1, we can notice that the total of sampled measurements during the 70 periods are less than the total of sampled measurements in the case of the static risk r = 0.9 and slightly less than the case of static risk r = 0.4. This indicates that our proposed approach has evaluated the HR risk level to less than 0.9. Thus, the amount of sampled data is reduced during the 70 periods which will have an impact on the energy consumption. While, if we compare the results obtained for patient 2, we can notice that the amount of sampled measurements in the case of dynamic risk is greater than the amount of samples in the case of static risk r = 0.4 and less than the amount of samples in the case of static risk r = 0.9. This indicates that our proposed approach has evaluated the HR risk level less than 0.9. Thus, the amount of sampled data is reduced during the 70 periods compared to a static risk of 0.9 and slighlty increased compared to a static risk of 0.4. Therefore, adapting the risk level value over time has an impact on the Fisher test result since the amount of sampled data changes from one period to another, as well as on the BV function and thus on the output of the sampling rate adaptation scheme. This will have an impact on the energy consumption.

2) Energy Consumption: Figure 6 shows the HR sensor node's remaining energy at the end of 70 periods for both patient 1 and 2. We suppose that 1 unit of energy is equal to 152 Joules : the sensing task consumes 6 Joules, the processing task consumes 24 Joules, the transmission task (TX) consumes 60 Joules and the receiving task (RX) consumes 62 Joules [START_REF] Magno | Energy-efficient context aware power management with asynchronous protocol for body sensor network[END_REF]. Thus, we suppose that each sampled measurement needs 0.6 units of energy to be sent to the coordinator. It includes 0.04 units of energy for the sensing task, 0.16 units of energy for the processing task and 0.4 units of energy for the transmission task. Having the initial energy randomnly set to 4000 units, the dynamic approach has consumed the least energy among the 3 cases in the case of patient 1. The remaining energy is 1463.8 units when setting the risk level to 0.9, is 1687 units when setting the risk level to 0.4 and 1753.6 units when adapting the risk level over the 70 periods. However, the results of the HR sensor node's sampling rate adaptation of patient 2 shows that adapting the risk level over time has reduced the energy consumption compared to setting the risk level value to 0.9 but it has increased the energy consumption compared to setting the risk level value to 0.4. The remaining energy is 1703.2 units when setting the risk level to 0.9, 2031.4 units when setting the risk level to 0.4 and 1888.6 units when adapting the risk level over the 70 periods. Thus, adapting the risk level of a vital sign over time adapts the lifetime of a sensor node to the patient's changing health condition.

C. Data Reduction vs Loss of information

In this section, we study the data reduction performed at the sensor node level following the sampling rate adaptation and discuss its impact on the loss of information. We compare the results obtained for the adaptive sampling rate scheme with static risk level and with dynamic risk level for patient s01840-3454-10-24-18-46nm. Table I shows the data reduction performed at the level of each sensor node for the three scenarios. Theoretically, if the Fisher Test null hypothesis is rejected for all 70 periods then a data reduction of only 30% is acheived given that the maximum sampling rate is SR max = 1 sample/1.42 sec (corresponding to 70 samples out of 100 per period). Whereas, if the Fisher Test null hypothesis is accepted for all 70 periods with F = 0 then a maximum data reduction of 70% is acheived given that the minimum sampling rate is SR max = 1 sample/3.33 sec (corresponding to 30 samples out of 100 per period) (cf. Algorithm 1). The results show that for both patients and for all the monitored vital signs the percentage of data reduction is bounded by 30% and 70% such as the data reduction is stricly higher than 30% demonstrating that the F-test is a suitable statistical test for this type of medical data. As shown in Table I, a total data reduction of about 48% compared to the original dataset has been acheived in our approach, of 43% in the case of a static risk r = 0.9 and of 47% in the case of a static risk r = 0.4.

We compare the scores of the measurements of the sampled datasets to the ones of the original dataset in each of the three scenarios, while taking the time granularity of the original dataset (1 sec). We choose the mean-squared error to measure the average of the squares of the errors or deviations given that

M SE = 1 M ×N M i=1 N j=1 (a(i, j) -b(i, j)) 2
where M is the total of monitored vital signs, N is the total of measurements for each vital sign during the simulation, a is the score of the j th measurement of the i th vital sign in the original dataset and b is the score of the j th measurement of the i th vital sign in the sampled dataset. The results (cf. Table II) show that the MSE of the three scenarios are very close to 1, thus the loss of information is negligible. Our approach has reduced the amount of sampled data 5% more than the case of static risk r = 0.9 with a MSE difference of 0.0017.

Table III shows the percentage of critical events detected over 70 periods compared to the critical events recorded in the original dataset. As we can see, the percentage is greater than about 70% for all vital signs for the three scenarios. However, it is important to note that the percentage of critical events in the original dataset for each of the HR, SpO2, BlOODT, RESP and ABPSys which do not last more than 2 seconds are respectively 38%, 1%, 20%, 5% and 17%. Thus, in three scenarios and for all vital signs, all of the persistant critical events have been detected. Finally, Figure 7 shows the HR sensor node's totals for each score for the 3 different scenarios in comparison with the totals of the original dataset. As we can see, although our proposed approach has reduced the amount of HR sampled data to 47% (cf. Table I), the totals of all scores were conserved. Consequently, our proposed approach allows a sensor node to dynamically adapt its risk level over time based on the state of the vital sign it is monitoring as well as the overall health condition of the patient. Fixing the risk level to a high value when it is not needed increases the energy consumption and the amount of sampled data over time whereas in some scenarios, fixing the risk level to a low value can increase the loss of information and some critical events can be passed unseen. V. CONCLUSION

In this paper, an adaptive sampling rate approach using dynamically adapted risk level has been proposed. Since the patient's health condition, remotly and continously monitored using a WBSN, changes over time, we have proposed to adapt the risk level of each sensor node based on the overall health condition of the patient and the state of their corresponding vital sign. The risk level represents the degree of monitoring importance given to a vital sign regarding the severity level of the patient's health condition. The results show that our approach, compared to an adaptive sampling rate approach with static risk levels, adapts the amount of sampled data and energy consumption according to the patient's status while preserving information about critical events as well as about the changing condition of the vital sign whether it is deteriorating or improving. As a future work, we intend to study other statistical tests which can be applied on medical data to study the variations of the measurements, propose other types of behavior function and compare them with our proposed adaptive sampling rate scheme. Moreover, a real implementation of the proposed approach is intended to be achieved in order to validate its performance on real-case scenarios.

Fig. 1 :

 1 Fig. 1: Behavior Function: Quadratic Bezier Curves

Fig. 3 :

 3 Fig. 3: Risk Level Evaluation Function

7 Fig. 4 : 9 Fig. 5 :

 7495 Fig. 4: ABPSys sensor node's risk level adaptation over 70 periods with different α and β setups

1 : 2 :Fig. 6 :

 126 Fig. 6: Comparison of the HR sensor node's remaining energy at the end of 70 periods between 3 different cases.

Fig. 7 :

 7 Fig. 7: The HR sensor node's totals for each score for three different scenarios in comparison with the original dataset.

TABLE I :

 I Data reduction performed for each monitored vital sign of record s01840-3454-10-24-18-46nm from MIMIC II compared to original dataset in 3 different cases.

	Risk Level Setup	HR	SpO2	BLOODT Resp	ABPSys
	Static Risk r = 0.9	40%	44%	65%	34%	34%
	Static Risk r = 0.4	45%	48%	65%	37%	37%
	Dynamic Risk	47%	48%	65%	38%	38%

TABLE II :

 II Mean squared error between original dataset and sampled dataset of patient s01840-3454-10-24-18-46nm from MIMIC II for 3 different scenarios.

	Risk Level Setup	MSE
	Static Risk r = 0.9	0.0128
	Static Risk r = 0.4	0.0144
	Dynamic Risk	0.0145

TABLE III :

 III Percentage of detected critical changes compared to original dataset of patient s01840-3454-10-24-18-46nm from MIMIC II for 3 different scenarios.

	Risk Level Setup	HR	SpO2	BLOODT Resp	ABPSys
	Static Risk r = 0.9	76%	97%	80%	97%	86%
	Static Risk r = 0.4	72%	96%	80%	97%	85%
	Dynamic Risk	69%	95%	80%	97%	86%

ACKNOWLEDGMENT

This project has been performed in cooperation with the Labex ACTION program (contract ANR-11-LABX-0001-01). This project has been funded with support from the National Council for Scientific Research in Lebanon.