Alejandro Reyes Amaro
email: alejandro.reyes@univ-nantes.fr

Eric Monfroy
email: eric.monfroy@univ-nantes.fr

Florian Richoux
email: florian.richoux@univ-nantes.fr

A Parallel-Oriented Language for Modeling Constraint-Based Solvers

This paper presents the Parallel-Oriented Solver Language (POSL, pronounced "puzzle"): a new framework to build interconnected meta-heuristic based solvers working in parallel. The goal of this work is to obtain a framework to easily build solvers and reduce the developing effort, by providing a mechanism for reusing code from other solvers. The novelty of this approach lies in looking at solver as a set of components with specific goals, written in a parallel-oriented language based on operators.

An interesting advantage of POSL is the possibility to share not only information, but also behaviors, allowing solvers modifications. POSL allows solver's components to be transferred and executed by other solvers. It provides an additional layer to dynamically define the connectivity between the solvers.

The implementation of POSL remains a work in progress, therefore this paper will focus on POSL concepts only.

Introduction

Combinatorial Optimization has important applications in several fields, including machine learning, artificial intelligence, and software engineering. In some cases, the main goal is only to find a solution, like for Constraint Satisfaction Problems (CSP). A solution will be an assignment of variables satisfying the constraints set. In other words: finding one feasible solution.

CSPs find a lot of application in the industry. For that reason, many techniques and methods are applied to the resolution of these problems. Although many of these techniques, like meta-heuristics, have shown themselves to be effective, sometimes problems we want to solve are too large, i.e. the search space is huge, and in most cases intractable.

However, the development of computer architecture is leading us toward massively multi/many-core computers. This evolution, closely related to the development of the super-computing, has opened a new way to find solutions for these problems in a more feasible manner, reducing the search times. Adaptive Search [START_REF] Diaz | Constraint-Based Local Search for the Costas Array Problem[END_REF] is one of the most efficient algorithm showing very good performances scaling to hundreds or even thousands of cores. It is an example of multi-walk local search method (i.e. algorithms that explore the search space by using independent search processes). For this algorithm, an implementation of a cooperative multi-walks strategy has been published in [START_REF] Munera | A Parametric Framework for Cooperative Parallel Local Search[END_REF]. These works have shown the efficiency of independent multi-walk, that is why we have oriented POSL towards this parallel scheme.

In the last years, a lot of efforts have been made in parallel constraint programming. In this field, the inter-process communication, useful to share information between solvers, is one of the more critical issues. In [START_REF] Hamadi | Control-Based Clause Sharing in Parallel SAT Solving[END_REF], an idea to include low-level reasoning components in the SAT problems resolution is proposed, dynamically adjusting the size of shared clauses to reduce the possible blow up in communication. The interaction between solvers is called solver cooperation and it is very popular in this field due to their good results. In [START_REF] Pajot | Separating Search and Strategy in Solver Cooperations[END_REF] is presented a paradigm that enables the user to properly separate strategies to follow to combine solvers applications in order to find the desired result, from the way that the search space is explored. Meta-S is a practical implementation of a theoretical framework proposed in [START_REF] Frank | Meta-S: A Strategy-Oriented Meta-Solver Framework[END_REF], which allows to tackle problems, through the cooperation of arbitrary domain-specific constraint solvers. POSL provides a mechanism of creating solver-independent communication strategies, giving the possibility of easily studying solving processes and results.

During the development of constraint programming field, much research focuses on fitting existing algorithms , trying to improve some metrics such as convergence speed and solution quality. However, it is required a deep study to find the right algorithm for the right problem. Hyperion [START_REF] Alexander | Hyperion 2. A toolkit for {meta[END_REF] is a Java framework for meta-and hyper-heuristics built with the principles of interoperability, generality by providing generic templates for a variety of local search and evolutionary computation algorithms, and efficiency allowing rapid prototyping, with the potential use of the same source code. POSL aims to offer these advantages, but provides also a mechanism to define communication protocols between solvers.

In this paper we explain a methodology to use POSL for easily build many and different cooperating solvers based on coupling four fundamental and independent components: operation modules, open channels, the computation strategy and communication channels or subscriptions. Recently, the hybridization approach, leads to very good results in constraint satisfaction [START_REF] El-Ghazali | Combining metaheuristics with mathematical programming, constraint programming and machine learning[END_REF]. For that reason, since the solver's component can be combined, POSL is designed to give the possibility to obtain sets of different solvers to be executed in parallel.

POSL provides, through a simple operator-based language, a way to create a computation strategy, combining already defined components (operation modules and open channels). A similar idea was proposed in [START_REF] Fukunaga | Automated discovery of local search heuristics for satisfiability testing[END_REF] without communication, where it is investigated an evolutionary approach that uses a simple composition operator, to automatically discover new local search heuristics for SAT, visualizing them as combinations of a set of building blocks. In the last phase of the coding process using POSL, solvers can be connected each others, depending on the structure of their open channels, and this way, they can share not only information, but also their behavior, giving the possibility of send/receive their operation modules. This approach gives to the solvers the possibility to evolve during the execution.

The work that we present in this paper focuses on the concepts of POSL only. It is part of an ongoing investigation that foresees a complete implementation.

Target problems

POSL is a framework to tackle Constraint Satisfaction Problems (CSP). A CSP is defined by a triple xX, D, Cy where X " tx 1 , x 2 , . . . , x n u is a finite set of variables, D " tD 1 , D 2 , . . . , D n u, is the set of associated domains to each variable in X, and C " tc 1 , c 2 , . . . , c m u, is a set of constraints. Each constraint is defined involving a set of variables, and specifies the possible combinations of values for these variables. We denote by P a CSP.

A configuration s P D 1 ˆD2 ˆ¨¨¨ˆD n , is a combinations of values for the variables of X. We say that s is a solution of P if and only if s mets all the constraints c i P C. exchange, but also to execute the components in parallel. In this stage, the information that is interesting to share with other solvers, is sent using operators for that purpose. This will be the solver's backbone.

4. Solvers are defined by instantiating and connecting the strategy, operation modules and open channels, and by connecting them each others (Figure 1d).

In this section we will explain in details each step.

Operation module

An operation module is the most basic and abstract way to define a piece of computation. It receives an input, executes an internal algorithm and gives an output. In this paper, we use this concept to describe and define one of the basic components of a solver, and they will be joined through computation strategies.

An operation module will represent a piece of the algorithm that can be changed or modified during the execution. They can be dynamically replaced by or combined with others operation modules, since the operation modules can be also information to be shared between solvers working in parallel. This way the solver can change/mutate its behavior during the execution, by combining its operation modules with operation modules from other solvers. They are represented by blue blocks in the Figure 4.

Definition 1 An operation module

Om is an mapping defined as follow:

Om : D Ñ I (1)
In ( 1), the nature D and I will depend on the type of operation module. They can be either a configuration, or a set of configurations, or a set of values of some data type, etc.

Open channels

The open channels are the solver's components in charge of the information reception in the communication between solvers. They can interact with the operation modules, depending on the computation strategy. The open channels will play the role of outlets, allowing solvers to be connected to them and receiving information. They are represented by red shapes in the Figure 4.

Definition 2 An open channel

Ch is an mapping defined as follow:

Ch : D Ñ D (2)
An open channel is an open door to outside to receive information. Inside this component will not exist any processing or data transformation. That is way in (2), the domain and the image are the same.

Computation strategy

The computation strategy is the solver's backbone. It joins the operation modules and the open channels coherently. It is independent from the operation modules and open channels used in the solver. It means they can be changed or modified during the execution, without altering the general algorithm, but still respecting the main structure.

Definition 3

We call module to (and it is denoted by the letter M):

1. An operation module or

2. An open channel or 3. rM 1 OP M 2 s: The composition of two modules M 1 and M 2 to be executed sequentially, one after the other, and returning an output depending on the nature of the operator OP; or

4. M 1 OP M 2 :
The composition of two modules M 1 and M 2 to be executed and returning an output depending on the nature of the operator OP. These two modules will be executed in parallel if and only if OP supports parallelism, (i.e. some modules will be executed sequentially although they were grouped this way); or sequentially otherwise.

We denote the space of the modules by M. We will call compound modules to the composition of modules described in 3 and 4.

Following we will define some operators included in POSL framework. In order to group modules, like in Definition 3(3) and 3(4), we will use |.| as generic grouper.

The following operator allows us to execute two modules sequentially one after the other.

Definition 4 Suppose two different modules: i) M 1 : D 1 Ñ I 1 and ii) M 2 : D 2 Ñ I 2 . Then, we can obtain a compound module through the following parametric operator:

Þ ÝÑ pM 1 , M 2 q : D 1 Ñ I 2 The operation |M 1 Þ ÝÑ M 2 | can be performed if and only if I 1 Ď D 2 . The operation |M 1 Þ ÝÑ M 2 |
returns a compound module as result of the sequential execution of M 1 followed by M 2 .

The following operator is useful to execute modules sequentially creating bifurcations, subject to some boolean condition:

Definition 5 Suppose three different modules: i) M 1 : D 1 Ñ I 1 , ii) M 2 : D 2 Ñ I 2 and iii) M 3 : D 3 Ñ I 3 . Then,
we can obtain a compound module through the following parametric operator:

Þ Ý Ý Ñ ă.ą pM 1 , M 2 , M 3 q : t0, 1u ˆD1 Ñ I o The operation ˇˇM 1 Þ Ý Ý Ñ ă.ą tM 2 , M 3 u ˇˇcan be performed if and only if: I 1 Ď D 2 X D 3 and I 2 Y I 3 Ď I o .
The operation ˇˇˇM 1 Þ ÝÝÝÝÑ ăcondą tM 2 , M 3 u ˇˇˇr eturns a compound module as result of the sequential execution of M 1 followed by M 2 if ă cond ą is true or by M 3 otherwise.

We can execute modules sequentially creating also cycles, by defining compound modules with another conditional operator: Definition 6 Suppose a module M 1 : D 1 Ñ I 1 . Then, we can obtain a cyclic compound module through the following parametric operator: pM 1 q : t0, 1u ˆD1 Ñ I 1

The operation pă cond ąq tM 1 u can be performed if and only if I 1 Ď D 1 . The operation pă cond ąq tM 1 u returns a compound module as result of the sequential execution of M 1 while ă cond ą remains true.

In Figure 2 we present a simple example of how combining modules using POSL operators introduced above. Algorithm 1 shows the corresponding code.

r ă stop cond ą t " M 1 Þ ÝÝÝÝÑ ăcondą tM 2 ; M 3 u  Þ ÝÑ M 4 us
So far, we have assumed that we can only execute two modules, by using the output of the first one as input of the second one; but sometimes some modules need some other information as input coming from module(s) executed not right before them. In that case, it is necessary to give additional inputs:

Definition 7 Suppose two different modules: i) M 1 : D 1 Ñ I 1 and ii) M 2 : D 2 Ñ I 2 .
Then, we can obtain a compound module through the following parametric operator:

ˆpM 1 , M 2 q : D o Ñ I 1 ˆI2
The operation ˇˇM 1 ˆM2 ˇˇcan be performed if and only if:

D o Ď D 1 X D 2 .
The operation ˇˇM 1 ˆM2 ˇˇreturns a compound module as result of the execution of M 1 and M 2 , and its image will be the cartesian product of its operand's images.

Figure 3 shows a scenario where it is necessary to use this operator. POSL offers the possibility to give variability to the solvers. Depending on the operator, one or both operand module will be executed, but only the output of one of them will be returned by the compound module. We present these operators in two definitions, grouping those which execute only one module operand (Definition 8) and those which execute both (Definition 9). Definition 8 Suppose two different modules: M 1 , M 2 : D Ñ I and a probability ρ. Then, we can obtain compound modules through the following parametric operators:

1. ρ pρ, M 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 ρ M 2 ˇˇexecutes " M 1 with probability ρ M 2 with probability p1 ´ρq 2. _ pM 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 _ M 2 ˇˇexecutes " M 1 if M 1 is not null M 2 otherwise
In the following Definition, the concepts of cooperative parallelism and competitive parallelism are implicitly included. We say that cooperative parallelism exists when two or more processes are running separately, they are independent, and the general result will be some combination of the results of all the involved processes (e.g. Definitions 9.1 and 9.2). On the other hand, competitive parallelism arise when the general result is the solution of the first process terminates (e.g. Definition 9.3). Definition 9 Suppose two different modules: M 1 , M 2 : D Ñ I. Suppose also that o 1 and o 2 are the outputs of M 1 and M 2 respectively, and that there exist some order criteria between them. Then, we can obtain a compound module through the following parametric operators:

1. M pM 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 M M 2 ˇˇreturns max to 1 , o 2 u 2. m pM 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 m M 2 ˇˇreturns min to 1 , o 2 u 3. pM 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 M 2 ˇˇreturns " o 1 if M 1 ends first o 2 otherwise
Algorithm 2: POSL code for Figure 3

M 1 Þ ÝÑ "" M 4 M 5 ı ˆrM 2 Þ ÝÑ M 3 s ı Þ ÝÑ M 6
The operators presented in Definitions 8 and 9 are very useful in terms of sharing not only information between solvers, but also behaviors. If one of the operands is an open channel it can receive an external solver's operation module, providing the opportunity to instantiate it in the current solver. The operator will either instantiate that module if it is not null, and execute it; or it will execute the other operand module otherwise.

Some others operators can be useful dealing with sets, especially for population-based solvers:

Definition 10 Suppose two different modules: M 1 , M 2 : D Ñ 2 I . Where 2 I is to denote sets of some data type. Suppose also that the sets V 1 and V 2 are the outputs of the execution of M 1 and M 2 . Then we can obtain a compound module through the following parametric operators:

1. Y pM 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 Y M 2 ˇˇreturns V 1 Y V 2 2. X pM 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 X M 2 ˇˇreturns V 1 X V 2 3. ´pM 1 , M 2 q : D o Ñ I o Where: ˇˇM 1 ´M2 ˇˇreturns V 1 zV 2
Now, we define the operators that allow us to send information to outside, i.e. others solvers. We can send two types of information: i) we can execute the operation module and send its output, or ii) we can send the operation module itself. . This utility is very useful in terms of sharing behaviors between solvers. In this stage, and using these operators, we can create the algorithm responsible for managing the components to find the solution for a given problem. These algorithms are fixed, but generic w.r.t. their components (operation modules and open channels). It means that we can build different solvers using the same strategy, but instantiating it with different components, as long as they have the right structure (i.e., input/output signature).

To define a computation strategy we will use the following environment:

oModule: < list of operation modules types > ; Before coding the computation strategy, it is necessary to declare which types of operation modules and open channels will be used. They are placed in < list of operation modules types > and < list of open channels types > respectively.

Inside the brackets, the field < ...computation strategy... > corresponds to the Parallel-Oriented Solver Language based on operators combining the modules. A clear example is provided in Algorithm 3.

Solver definition

With operation modules, open channels and computation strategy already defined, we can create the solvers by instantiating the mentioned components. POSL provides an environment to do that:

Communication definition

Once we have defined the strategy to our solver, the next step is to declare the communication channels, i.e. connecting the solvers each others. Up to here, the solvers are disconnected, but they have everything to establish the communication (Figure 1c). In this last stage, POSL provides to the user a platform to easily define the cooperation meta-strategy to be followed by the set of solvers declared before.

The communication is established by following the next rules guideline:

1. Each time a solver sends any kind of information by using the operator . o or . m , it creates a communication jack

Each time a solver places an open channel into its definition, it creates a communication outlet

3. Solvers can be connected each others by creating subscriptions, connecting communication jacks with communication outlet (see Figure 4).

With the operator p¨q we can have access to the operation modules sending information and to the open channel's names in a solver. For example: Solver 1 ¨M1 provides access to the operation module M 1 in Solver 1 if and only if it is affected by the operator . o (or . m ), and Solver 2 ¨Ch 2 provides access to the open channel Ch 2 in Solver 2 . Tacking this into account, we can define the subscriptions.

Definition 12 Suppose two different solvers: Solver 1 and Solver 2 . Then, we can connect them through the following operation:

Solver 1 ¨M1 Solver 2 ¨Ch 2
The connection can be defined if and only if:

An POSL solver

In this section we explain the structure of a solver created by POSL through an example. We choose one of the more classic solution methods for combinatorial problems: local search meta-heuristics algorithms. These algorithms have a common structure: they start by initializing some data structures (e.g. a tabu list for Tabu Search [START_REF] Gendreau | Tabu Search[END_REF], a temperature for Simulated Annealing [START_REF] Nikolaev | Simulated Annealing[END_REF], etc.). Then, an initial configuration s is generated (either randomly or by using an heuristic). After that, a new configuration s ˚is selected from the neighborhood V psq . If s ˚is a solution for the problem P, then the process stops, and s ˚is returned.

If not, the data structured are updated, and s ˚is accepted or not for the next iteration, depending on some criterion (e.g. penalizing features of local optimums, like in Guided Local Search [START_REF] Boussaïd | A survey on optimization metaheuristics[END_REF]). We can combine modules to create more complex ones. For example, if we want a operation module generating a random configuration but able to give sometimes a fixed configuration (one passed by parameter) following a given probability, our operation module 1 would be the combination of this two modules through the operator ( ρ ).

Let's make this solver a little bit more complex: suppose that we have a neighborhood function which is exploitation-oriented, but sometimes we will ask other solvers to send their neighborhood operation modules to be executed in the current solver, as a mechanism to scape from local minimum. This behavior is modeled by the following open channel: Open Channel ´1 : Asking for V psq operation module.

Let's suppose also that we would like to always broadcast our current solution. It is possible by applying the operator . o to the operation module 4. Figure 4 presents the example above. The Solver-1's open channel is represented by a red shapes. The Solver-1 has an open channel asking for additional neighborhood function, so it means that it can be connected with Solver-2, because Solver-2 has applied the operator . o to send (somehow/at some point) its neighborhood function operation module. In the solver definition layer it is created a subscription (represented by a green block) to define the communication channel. As we can see in the diagram, there is nobody attending to the operation module sending the best configuration. In that case, there are no links outside the channel, and no information will be sent.

Algorithm 3 shows POSL code for the computation strategy of the solver described above (Solver-1), using predefined operators. POSL provides information regarding the execution process, such as number of Iterations, solver execution Time, Best found solution, among others.

: M a 1 , M b 1 , M 2 , M 2 , M 4 ; oChannel : Ch 1 ; { r pIterations % 10000 != 0q t " M a 1 ρ M b 1 ı Þ ÝÑ r pIterations % 1000 != 0q t " Ch 1 _ M 2 ı Þ ÝÑ M 3 Þ ÝÑ M 4 o us us } Algorithm 4: Local search meta- heuristic solver definition Σ 1 Ð solver { strategy : St oModule : m a 1 , m b 1 , m 2 , m 3 , m 4 oChannel : ch 1 }
Algorithm 4 shows the solver definition. We place here the computation strategy, followed by instances of the modules (operation modules and open channels). The instances have to respect the type of the modules declared in the computation strategy. Supposing that there exist another solver Σ 2 with an operation module called M V sharing its neighborhood function, we can connect it with the solver Σ 1 as the Algorithm 5 shows.

Conclusions

In this paper we have presented POSL, a framework for building cooperating solvers. It provides an effective way to build solvers ables to exchange any kind of information, even other solver's behavior, sharing their operation modules. Using POSL, many different solvers can be created and ran in parallel, using only one generic strategy, but instantiating different operation modules and open channels for each other.

It is possible to implement different communication strategies, since POSL provides a layer to define communication channels connecting solvers dynamically using subscriptions.

At this point, the implementation of POSL remains in progress, in which our principal task is creating a design as general as possible, allowing to add new features in the near future. Our goal is obtaining a rich library of operation modules and open channels to be used by the user, based on a deep study of the classical meta-heuristics algorithms for solving combinatorial problems, in order to cover them as much as possible. In such a way, building new algorithms by using POSL will be easier.

At the same time we pretend to develop new operators, depending on the new needs and requirements. It is necessary, for example, to improve the solver definition language, allowing the process to build sets of many new solvers to be faster and easier. Furthermore, we are aiming to expand the communication definition language, in order to create versatile and more complex communication strategies, useful to study the solvers behavior.

As a medium term future work, it would be interesting also to include learning techniques, allowing the solver to change automatically, depending for instance on the results of their neighbor solvers.

Figure 1 :

 1 Figure 1: Building parallel solvers with POSL

Figure 2 Algorithm 1 :

 21 Figure 2

Definition 11

 11 Given an operation module M : D Ñ I, we can define the following parametric operators: 1. . o pMq : D Ñ I Where: M o executes M and sends its output to outside 2. . m pMq : D Ñ I Where: M m executes M and also sends the operation module itself to outside

3 oChannel: 5 <

 35 < list of open channels types > ; 4 { ...computation strategy... > 6 }

1

  solver_k := solver 2 { 3 cStrategy: < strategy > ; 4 oModule: < list of operation modules instances > ; 5 oChannel: < list of open channels instances > ; 6 }

1 .

 1 Solver 1 has an operation module called M 1 encapsulated into an operator . o or . m . 2. Solver 2 has an open channel called Ch 2 receiving the same type of information sent by M 1 .

  Restarts are classic mechanisms to avoid becoming trapped in local minimum. They are trigged by reaching no improvements or a timeout. The operation modules composing local search meta-heuristics are described bellow: Operation Module ´1 : Generating a configuration s Operation Module ´2 : Defining the neighborhood V pSq Operation Module ´3 : Selecting s ˚P V psq Operation Module ´4 : Evaluating an acceptance criteria for s

Figure 4 : 3 :

 43 Figure 4: Diagram of a solver cooperation strategy

Algorithm 5 :

 5 Inter-solvers communication definitionΣ 2 ¨MV Σ 1 ¨Ch 1

POSL parallel solversPOSL proposes a solver construction platform, following different stages:1. The solver algorithm is modeled by decomposing it into small pieces/modules of computation. After that, they are implemented as separated functions. We will give the name of operation module to these pieces/modules (Figure1a).2. The next step is to decide what information is interesting to receive from other solvers. This information is encapsulated into other objects called open channels, allowing data transmission among solvers (Figure1a).3. A generic strategy is coded through the Parallel-Oriented Solver Language based on operators it provides, using the mentioned components in the stages 1 and 2, allowing not only the information

St := cStrategy


Mohamed Boussaa
email: mohamed.boussaa@inria.fr

Olivier Barais
email: olivier.barais@inria.fr

Gerson Sunye
email: gerson.sunye@inria.fr

Benoit Baudry
email: benoit.baudry@inria.fr

A Novelty Search-based Test Data Generator for Object-oriented Programs

Keywords: Search-based Software Testing, Structural Coverage, Automated Test Data Generation, Novelty Search, Genetic Algorithm

In search-based structural testing, meta-heuristic search techniques have been frequently used to automate test data generation. In this paper, we introduce the use of novelty search algorithm to the test data generation problem based on statement-covered criterion. In this approach, we seek to explore the search space by considering diversity as the unique objective function to be optimized. In fact, instead of having a fitness-based selection, we select test cases based on a novelty score showing how different they are compared to all other solutions evaluated so far.

INTRODUCTION

In general, manually creating test cases for testing software systems is time consuming and error-prone. Thus, the automation of this process becomes necessary. In fact, meta-heuristic search techniques such as Genetic Algorithms (GAs) are frequently used to automate the test data generation process and gather relevant test cases through the wide search space [START_REF] Harman | A theoretical and empirical study of search-based testing: Local, global, and hybrid search[END_REF][START_REF] Ali | A systematic review of the application and empirical investigation of search-based test case generation[END_REF]. These techniques are especially applied for structural white-box testing. For coverage-oriented approaches, applying Evolutionary Algorithms (EAs) to test data generation [START_REF] Lakhotia | Automated test data generation for coverage: Haven't we solved this problem yet?[END_REF] has been focused on finding input data for a specific path of the program in accordance with a coverage criterion (e.g., longest path executed). The problem with coverage-oriented approaches is that search-based techniques cannot exploit the huge space of possible test data. In fact, some structures of the system may not be reached since they are executed only by a small portion of the input domain. The use of a fitness function as a coverage criterion to guide the search to detect relevant test data usually create one or more local optima to which the search may converge.

Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for third-party components of this work must be honored. For all other uses, contact the owner/author(s). This process is known as diversity loss. Therefore, diversity maintenance in the population level is key for avoiding premature convergence. This issue is a common problem in GAs and many techniques are proposed to escape from local optima in multi-modal search space [START_REF] Das | Real-parameter evolutionary multimodal optimizationâ ȂŤ a survey of the state-of-the-art[END_REF][START_REF] Pandey | A comparative review of approaches to prevent premature convergence in ga[END_REF]. However, all these alternatives use a fitness-based selection to guide the search.

GECCO '15

In this paper, we propose the use of Novelty Search (NS) algorithm [START_REF] Lehman | Exploiting open-endedness to solve problems through the search for novelty[END_REF] to the test data generation problem. In this approach, we seek to explore the search space of possible test input values by considering diversity as the unique objective function to be optimized. In fact, instead of having a fitnessbased selection, we rather select test cases based on a novelty score showing how different they are compared to all other test data evaluated so far.

The paper is organized as follows: section 2 describes the approach overview and our novelty search adaptation. Concluding remarks and future work are provided in section 3.

NOVELTY SEARCH FOR TEST DATA GENERATION 2.1 Approach Overview

Our test data generation framework aims to fully automate the test data generation process without requiring user intervention. In fact, we have considered our System Under Test (SUT) as a gray or semi transparent box in where the internal structure is partially known. Thus, we can design test cases through the exposed interfaces and conduct a code coverage analysis from the general structure of our target SUT. For example, within apache.commons.math library, Methods Under Test (MUTs) are accessed through some specific Java interfaces. Each interface belongs to a subpackage of the whole library and exhibits a set of methods. Our test data generation framework will rely on this concept to generate automatically test cases. In fact, as shown in Figure 1, starting from an input interface and a code source package, the testing framework is able to: (1) generate automatically sequences of method invocation through the input interface, (2) generate relative test data using NS, (3) execute test cases on target Classes Under Test (CUTs) and then (4) analyze code coverage within the source package. The process is iterated until a termination criterion is met (e.g., number of iterations)

In our proposal, we use the same logic as GAs to implement the NS algorithm. However, we make some changes and we add new settings to optimize the test data generation. In fact, instead of using a fitness function to evaluate generated test cases, we define a new measure of novelty to maximize. We replace, as well, the fitness-based selection (for fittest test cases) with a novelty-based one. This may favor the diversity of generated test data. Finally, we add an archive that acts as a memory of previously generated test cases. This archive is used to calculate the novelty metric. Otherwise, we keep the same settings as GAs, namely the crossover and mutation operators.

For our NS adaptation, we have to define the number of iterations N, the population size and the minimum coverage value. This latter defines the threshold of covered statements that should be reached. Test cases that exceed this threshold are automatically added to the set of relevant test cases. As well, we have to define a novelty threshold T that defines the threshold for how novel a test suite has to be before it is added to the archive. In addition, we define a maximum size limit for the archive L and a k number that will be used in calculating the novelty metric.

Novelty Metric

The Novelty metric expresses the sparseness of a test suite. It measures its distance to all other test cases in the current population and to all test cases that were discovered in the past (i.e., test cases in the archive). This measure expresses how unique the test suite is. We can quantify the sparseness of a set of test cases as the average distance to the k-nearest neighbors. The distance between two test suites is computed as a Manhattan distance between the input parameter values of all methods tested in the test suite. A candidate solution represents the set of methods signatures declared in the interface. Thereby, we represent this solution as a vector where each dimension has a method name, a list of parameters types and a list of test data (the genotype) Formally, we define this distance between two solutions as follows :

distance(t1, t2) = m i=1 p j=1 |t1(Mi, Pj) -t2(Mi, Pj)| (1)
where t1 et t2 are two selected test suites (solutions), m is the number of methods composing the test suite, p is the number of parameters composing a method. The couple (Mi, Pj) returns the j th parameter value of the i th method M relative to a test suite (t1 or t2). Since we are using Java primitive types for test data such as floats, integers, doubles, etc, it is easy to calculate this distance for numerical parameters values. However, for string data types we use the Levenshtein Algorithm1 to measure the strings distance.

To measure the sparseness of a test suite, we will use the previously defined distance to compute the average distance of a test suite to its k-nearest neighbors. In this context, we define the novelty metric of a particular solution as follows:

N M (S) = 1 k k i=1 distance(S, µi) (2)
where µi is the i th nearest neighbor of the solution S within the population and the archive of novel individuals. Finally, we have normalized the resulted novelty distance in the range [0-100].

CONCLUSION AND FUTURE WORK

In this paper, we proposed a new adaptation of the NS algorithm to test data generation problem. Using the NS approach is clearly a divergent evolutionary technique, inspired by natural evolution's drive to novelty that directly rewards novel behaviors instead of progress towards a fixed objective.

As a future work, we aim to conduct an empirical evaluation of our NS approach by comparing it to fitness-based and random approaches. In addition, we can optimize our approach by adding diversity as an addition goal to a traditional objective driven approach to form a multi-objective optimization problem. Finally, since we are testing graybox systems, we can apply this approach, as well, for blackbox testing. In this case, we will be able to measure some non-functional properties such as memory usage and CPU consumption.

Figure 1 :

 1 Figure 1: Approach overview

http://www.levenshtein.net/


Christophe Steiner

Michel Mehrenberger

Nicolas Crouseilles

Philippe Helluy

Quasi

come

Introduction

In strongly magnetized plasma, when collisions are not the dominant collective effect, one has to deal with kinetic models. Fluid models, which assume that the distribution function is close to an equilibrium, are no more suitable. However, the numerical solution of Vlasov type models is challenging since this model involves six dimensions in the phase space. Moreover, multi-scaled phenomena make the problem very difficult. Gyrokinetic theory enables to get rid of one of these constraints since the explicit dependence on the phase angle of the Vlasov equation is removed through gyrophase averaging while gyroradius effects are retained. The so-obtained five-dimensional function should be self-consistently coupled with Maxwell equations. In the following we consider the electrostatic approximation, where Maxwell equations are reduced to a Poisson equation (or its asymptotic counterpart, the so-called quasi-neutrality equation). The work presents here will still be applicable taking into account electromagnetic effects. The present work is devoted to the numerical resolution of the quasi-neutrality equation in a polar mesh which is the framework in the gyrokinetic code Gysela [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF]. We propose two alternatives to the classic Padé approximation by using (i) a Padé approximation of higher order or (ii) a direct method in the space configuration using interpolation operator.

The paper is organized as follows. In Section 2, we introduce the quasi-neutrality equation in the gyrokinetic framework. The numerical methods are described in Section 3 for the classic and new Padé approximations and in Section 4 for the method based on interpolation. These methods will be applied to gyrokinetic simulations in Section 5.

Quasi-neutrality equation in the gyrokinetic framework

The computational effort to numerically solve the 6 dimensional Vlasov-Maxwell systems which describes plasma turbulence in tokamak plasmas still remains out of reach for present day supercomputers. All the numerical simulations performed until now in this domain take care of the gyrokinetic ordering to reduce this problem of one dimension. This ordering takes into account the fact that (i) electromagnetic fluctuations occur on time scales much longer than charged particle gyration period (ω/Ω c 1 with ω the fluctuation frequency and Ω c the cyclotron frequency), and (ii) the wavelength of these fluctuations is much smaller than the characteristic scale length of the gradients of magnetic field, density and temperature. See [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF] for a detailed review on the gyrokinetic framework and simulations to compute turbulent transport in fusion plasmas. Within this gyro-ordering, the so-called gyrokinetic model can be derived (see [START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF]) by averaging on the fast gyration of charged particles around the magnetic field lines. The magnetic toroidal configuration considered in this paper is simplified. Indeed, magnetic flux surfaces are assumed concentric tori with circular cross-sections. The gyroaverage operator occurs in this reduction from 6 to 5 dimensions. The new 5D set of coordinates corresponds to: (i) 3D toroidal spatial coordinates (r, θ, ϕ) (with r the radial direction, θ and ϕ the poloidal (resp. toroidal) angle), and (ii) 2D in velocity space with v the velocity parallel to the magnetic field line and µ = m|v ⊥ | 2 /(2B) the magnetic moment where v ⊥ represents the velocity in the plane orthogonal to the magnetic field of amplitude |B|. It is important to note that in this ordering µ is an adiabatic invariant, so it plays the role of a parameter in the 5D gyrokinetic Vlasov equation.

In the following, the 4D problem which is treated corresponds to the case where we consider a unique value of µ, i.e. the same Larmor radius is taken for all particles.

T i , T e and n 0 refer to ion temperature, electron temperature and density profiles and will be defined at the beginning of Section 5. Numerical solutions are computed using normalized equations. The temperature is normalized to T e0 , where T e0 is defined by the initial temperature profile such that T e (r p )/T e0 = 1 where r p ∈ [r min , r max ] is the radial coordinate of the peak of the initial distribution function. The time is normalized to the inverse of the ion cyclotron frequency ω c = e i B 0 /m i . Velocities, including the parallel velocity, are expressed in units of the ion speed v T 0 = T e0 /m i , the electric potential is normalized to T e0 /e i and the magnetic field is normalized to B 0 . Consequently, lengths are normalized to the Larmor radius ρ = m i v T 0 /e i B 0 and the magnetic moment µ to T e0 /B 0 .

The time evolution of the gyrocenter distribution function f is given by the gyrokinetic conservative equation (see also Eqs ( 17)- [START_REF] Hatzky | Energy conservation in a nonlinear gyrokinetic particle-in-cell code for ion-temperature-gradientdriven modes in theta-pinch geometry[END_REF] in [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]):

B * ∂f ∂t + ∇ ∇ ∇ • B * dx G dt f + ∂ ∂v G B * dv G dt f = 0 (2.1)
where x G and v G are respectively the space coordinates and the parallel velocity of the guiding centers. In the electrostatic limit, for a particle of mass m and charge q the motion equations of the guiding centers are given by

dx G dt = v G b * + E × B + v D (2.2) m dv G dt = -µ∇ ∇ ∇ * B -q∇ ∇ ∇ * Φ + mv G E × B • ∇ ∇ ∇B B (2.3)
where ∇ ∇ ∇ * ≡ b * • ∇ ∇ ∇, while b * and B * are defined by:

b * ≡ B B * + mv G qB * B ∇ ∇ ∇ × B (2.4) B * ≡ B + mv G qB b • (∇ ∇ ∇ × B) . (2.5)
The drift reads

E × B = (1/B * )b × ∇ ∇ ∇ Φ while curvature drift is defined as v D = mv 2 G +µB qB * b × ∇ ∇ ∇B B .
The obtained function in five dimensions must be self-consistent coupled with Maxwell's equations. In the following, we consider the electrostatic approximation, where Maxwell's equations are reduced to a quasi-neutrality equation which is asymptotically equivalent to the Poisson equation. Since the Poisson equation is defined on the particles coordinates, the resolution of the gyrokinetic Vlasov-Poisson system requires an operator which transforms the gyrocentered phase space to the particles phase space. Let ρ be the gyro-radius which is transverse to b and which depends on the gyrophase α ∈ [0, 2π], i.e. ρ = ρ(cos(α) e ⊥1 + sin(α) e ⊥2 ).

Here e ⊥1 and e ⊥2 are the unit vectors of a cartesian basis in the plane perpendicular to the magnetic field direction b. Let x G be the guiding-center radial coordinate and x the position of the particle in the real space. These two quantities differ by a Larmor radius ρ, i.e

x = x G + ρ. Let f : (r, θ) ∈ R + × R → f (r, θ) be a polar function and g : (x 1 , x 2 ) ∈ R 2 → g(x 1 , x 2
) the function defined by g(r cos(θ), r sin(θ)) = f (r, θ) for all (r, θ). The function f (resp. g) stands for field quantity defined on x with polar (resp. cartesian) coordinates. The gyroaverage J ρ (f ) of f depending on the spatial coordinates is defined as

J ρ (f )(r, θ) = 1 2π 2π 0 g( x G + ρ)dα
where x G = r(cos(θ), sin(θ)). This gyroaverage process consists in computing an average on the Larmor circle. It tends to damp any fluctuation which develops at sub-Larmor scales. This operator is studied in [START_REF] Steiner | Gyroaverage operator for a polar mesh[END_REF] and we also refer to an abundant literature on this subject (see [START_REF] Brizard | Nonlinear gyrokinetic Maxwell-Vlasov equations using magnetic coordinates[END_REF][START_REF] Hahm | Nonlinear Gyrokinetic Equations for Tokamaks Microturbulence[END_REF][START_REF] Lin | Method for solving the gyrokinetic Poisson equation in general geometry[END_REF] and references in these articles).

The gyroaveraged electrostatic potential is solution of the self-consistently coupled 3D quasi-neutrality equation:

n 0 (r) T i (r) R + (Φ -J 2 √ 2µT i (Φ)) exp(-µ)dµ + n 0 (r) T e (r) (Φ -λ Φ ) = ni -n 0 (r). (2.6)
λ is an integer which is expand to 0 or 1 according to the fact that zonal flows are token into account or not.

We consider a simplified model of the system of equations (2.1)-(2.6). A periodic cylindrical plasma of radius a and length 2πR (with R the major radius) is considered as a limit case of a stretched torus. The plasma is confined by a strong magnetic which is uniform B = Be z where e z stand for the unit vector in the toroidal direction z. With theses assumptions the velocity drifts are reduced to the E × B drift. This SLAB 4D case is equivalent to the one treated in [START_REF] Klein | Bertrand Gyrokinetic waterbag modeling of a plasma column: Magnetic moment distribution and finite Larmor radius effects[END_REF] or [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. The equation satisfied by the distribution function of ions f (t, r, θ, z, v) following the guiding center movement reads:

∂ t f - ∂ θ J √ 2µ Φ r ∂ r f + ∂ r J √ 2µ Φ r ∂ θ f + v∂ z f -∂ z J √ 2µ Φ ∂ v f = 0 (2.7) for (r, θ, z, v) ∈ [r min , r max ] × [0, 2π] × [0, L] × [-v max , v max ].

Method based on Padé approximations

Padé approximation is a classic way to solve numerically the quasi-neutrality equation (2.6). More precisely, we can use such an approximation in order to compute the term

Φ = R + J 2 √ 2µT i (Φ)e -µ dµ.
In Fourier space associated to the variable θ, the gyroaverage is reduced to the multiplication by the Bessel function J 0 of argument k ⊥ ρ. Indeed, the Fourier transform of J ρ (f ) can be written as

J ρ (f )( k) = 1 2π R 2 2π 0 f ( x G + ρ)dα e -i x G • k d x G = 1 2π 2π 0 R 2 f ( x G + ρ)e -i( x G + ρ)• k d x G e i ρ• k dα = 1 2π 2π 0 e ikρ cos(α-θ) dα f ( k) which leads to J ρ (f )( k) = J 0 (kρ) f ( k).
where J 0 is the Bessel function of the first kind and of order 0. This relation implies

Φ(k) = Γ 0 (|k| 2 T i ) Φ(k) where the function Γ 0 is defined by Γ 0 (k) := R + exp(-x 2 /2)J 2 0 (x √ k)xdx.
Then, the left part of the quasi-neutrality equation (2.6) becomes

R + (Φ -J 2 √ 2µT i (Φ)) exp(-µ)dµ = (1 -Γ 0 (|k| 2 T i )) Φ.
Let I 0 be the modified Bessel function of first kind, one can prove that Γ 0 (k) = I 0 (k)e -k (see [START_REF] Lin | Method for solving the gyrokinetic Poisson equation in general geometry[END_REF]).

classic Padé approximation

The classic Padé approximation

1 1 -Γ 0 (x) ≈ 1 x
gives Γ 0 (|k| 2 T i ) ≈ 1 -|k| 2 T i and using this approximation, we obtain (see [START_REF] Hahm | Nonlinear Gyrokinetic Equations for Tokamaks Microturbulence[END_REF][START_REF] Dubin | Nonlinear gyrokinetic equations[END_REF][START_REF] Lee | On higher-order corrections to gyrokinetic Vlasov-Poisson equations in the long wavelength limit[END_REF]):

-∇ ⊥ • (n i ∇ ⊥ Φ) + n 0 T e (Φ -Φ ) = ni -n 0 .
A linearization of the diffusion term by assuming that n i (x) ≈ n 0 (r) (see [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Hahm | Nonlinear Gyrokinetic Equations for Tokamaks Microturbulence[END_REF]) leads to:

-∇ ⊥ • (n 0 ∇ ⊥ Φ) + n 0 T e (Φ -Φ ) = ni -n 0
and finally:

-∂ 2 r Φ + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r Φ + 1 r 2 ∂ 2 θ Φ + 1 T e (Φ -Φ ) = 1 n 0 (n i -n 0 ).
This Poisson equation can be solved using FFT in θ and finite differences in r. In the numerical results, we use centered differences of order 2 which leads to a tridiagonal system. As boundaries conditions, we consider homogeneous Neumann condition for the mode 0 and homogeneous Dirichlet condition for the other modes in r min and homogeneous Dirichlet condition for all modes in r max .

New Padé approximation

In the same spirit, we can increase the degree of the Padé approximation by considering a denominator of degree 2:

1 1 -Γ 0 (x) ≈    -4 3x 2 -4x -20x-36 7x 2 -36x .
These developments are not satisfactory since the denominator vanishes for a positive value (x = 4 3 for the first approximation and x = 36 7 for the second approximation). To overcome this problem, we are looking for an approximation of the form

1 1 -Γ 0 (x) ≈ ax + b x 2 + cx + d with a, b, c, d ∈ R. Since 1 1-Γ 0 (x) -→ x→0
+∞, we assume that d = 0. The following expansions:

ax + b x 2 + cx = b/c x + a -b/c c + O(x) 1 1 -Γ 0 (x) = 1 x + 3 4 + O(x)
imply b = c and a = 3 4 c + 1. Finally, by setting ε = 1/c, we obtain the following approximation:

1 1 -Γ 0 (x) ≈ 1 + ( 3 4 + ε)x x(1 + εx) .
So that the denominator does not vanish in R * + , we impose the condition ε ≥ 0 and moreover ε > 0 in order to damp the high frequencies. Note that this new approximation corresponds to the previous Padé approximation x → -20x-36 7x 2 -36x with ε = -7 36 . In Fig. 1, the function 1 1-Γ 0 is plotted, its classic approximation 1/x and the new approximation (with ε = 0.1 and 0.001) over the domain [0.2, 2]. We can see that the new Padé approximation fits better the function than the classic approximation. This new approximation gives

Γ 0 (|k| 2 T i ) ≈ 1 + (ε -1 4 )T i |k| 2 -εT 2 i |k| 4 1 + ( 3 4 + ε)T i |k| 2
.

Using this approximation of Γ 0 and assuming that Φ = 0, we obtain:

-εT 2 i ∆ 2 ⊥ Φ + T i 1 + T i T e 3 4 + ε ∆ ⊥ Φ - T i T e Φ = 3 4 + ε T i ∆ ⊥ -1 T i (n i -n 0 ) n 0
and finally:

-εT 2 i ∂ 4 r + 2 r 2 ∂ 2 r ∂ 2 θ + 1 r 4 ∂ 4 θ + 2 r ∂ 3 r - 2 r 3 ∂ r ∂ 2 θ - 1 r 2 ∂ 2 r + 4 r 4 ∂ 2 θ + 1 r 3 ∂ r Φ +T i 1 + T i T e 3 4 + ε ∂ 2 r + 1 r ∂ r + 1 r 2 ∂ 2 θ Φ - T i T e Φ = 3 4 + ε T i ∆ ⊥ -1 T i (n i -n 0 ) n 0 .
Again, this Poisson equation can be solved using FFT in θ and finite differences in r. In the numerical results, we use centered differences of order 2 which leads to a pentadiagonal system. As boundaries conditions, we consider homogeneous Dirichlet and homogeneous Neumann conditions in r min and r max .

Numerical method based on interpolation

We present in this part a method in the space configuration which consists in integrating on the gyrocircles using interpolation for solving the quasi-neutrality equation (2.6). This equation reads by making the change of variable μ = µT i :

R + 1 + T i T e Φ -J 2 √ 2μ (Φ) exp - μ T i dμ = T 2 i ni -n 0 n 0 .
This method is directly based on the computation of the gyroaverage described in [START_REF] Steiner | Gyroaverage operator for a polar mesh[END_REF] which consists to interpolate the distribution function in N uniformly distributed points on the Larmor circle. To do this, we consider a uniform polar mesh on the domain [r min , r max ] × [0, 2π] including N r × N θ cells:

C ij = [r i , r i+1 ] × [θ j , θ j+1 ], i = 0, . . . , N r , j = 0, . . . , N θ -1
where The computation of the gyroaverage writes, for (r i , θ j ) a point of the polar mesh:

r i = r min + i r max -r min N r , i = 0, ..., N r , θ j = j 2π N θ , j = 0, ..., N θ .
J ρ (Φ)(r i , θ j ) 1 N N -1 =0 P(Φ) r i cos (θ j ) + ρ cos 2 π N , r i sin (θ j ) + ρ sin 2 π N ,
where P is an interpolation operator. In the following, we will use the cubic splines interpolation. We make a radial projection on the boundaries for the points outside the domain and we consider 2π-periodic conditions in θ.

In order to detail the steps of the solver, we note Φ i,j := Φ(r i , θ j ), i = 0..N r , j = 0..N θ -1

J ρ Φ i,j := J ρ (Φ)(r i , θ j ), i = 0..N r , j = 0..N θ -1 φ := t (Φ 0,0 , ..., Φ 0,N θ -1 , Φ 1,0 , ..., Φ 1,N θ -1 , ..., Φ Nr,1 , ..., Φ Nr,N θ -1 )
J ρ (φ) := t (J ρ Φ 0,0 , ..., J ρ Φ 0,N θ -1 , J ρ Φ 1,0 , ..., J ρ Φ 1,N θ -1 , ..., J ρ Φ Nr,0 , ..., J ρ Φ Nr,N θ -1 ).

Construction of the matrix

A spl ∈ M (Nr+3)×N θ ,(Nr+1)×N θ (R) such that S = A spl φ
is the vector of splines coefficients. We consider Hermite splines with homoge-neous Neumann boundary condition in r and periodic conditions in θ. The matrix A spl in independent of the Larmor radius.

2. For each Larmor radius ρ j = 2μ j , construction of the matrix A contr ρ j ∈ M (Nr+1)×N θ ,(Nr+3)×N θ (R) giving the contribution of the gyroaverage of radius ρ j in each point in function of the splines coefficients. We have also J ρ j (φ) = A contr ρ j S and the matrix of the gyroaverage for the Larmor radius ρ j is given by G ρ j = A contr ρ j A spl . The matrix given the double gyroaverage is obtained by B ρ j = G 2 ρ j . We underline that for a given mesh in ρ j , these matrices can be computed once for all.

3. Evaluation of the integral by quadrature in μ:

R + (Φ -J 2 √ 2μ (Φ))e -μ/T i dμ ≈ μmax 0 (Φ -J 2 √ 2μ (Φ))e -μ/T i dμ ≈ Nµ j=1 c j (Φ -J 2 √ 2µ j (Φ))e -μ j i ≈   Nµ j=1 c j (Id -B ρ j )e -μ j /T i   Φ.
In numerical results, we use the Simpson quadrature.

4. Inversion of the matrix of the quasi-neutrality operator using LU decomposition.

Remark 1.

1. Note that all the previous steps are performed in precomputation.

2. An integration method is also presented in [START_REF] Lin | Method for solving the gyrokinetic Poisson equation in general geometry[END_REF] based on the approximation of the function Γ 0 .

We construct the matrix of the quasi-neutrality operator in the Fourier basis. in fact, the computation is more efficient in this basis since the matrix products, used for the construction of the double gyroaverage matrix, involve block diagonal matrix. This method was already used, for instance in [START_REF] Mehrenberger | Vlasov on GPU[END_REF]. More precisely, the periodicity in θ ensures that the matrices A contr ρ j and A spl are block circulant:

A spl =        A spl 0 A spl 1 . . . A spl N θ -1 A spl N θ -1 . . . . . . . . . . . . . . . . . . A spl 1 A spl 1 . . . A spl N θ -1 A spl 0        ∈ M (Nr+3)×N θ ,(Nr+1)×N θ (R) A contr ρ j =        A contr ρ j ,0 A contr ρ j ,1 . . . A contr ρ j ,N θ -1 A contr ρ j ,N θ -1 . . . . . . . . . . . . . . . . . . A contr ρ j ,1 A contr ρ j ,1 . . . A contr ρ j ,N θ -1 A contr ρ j ,0        ∈ M (Nr+1)×N θ ,(Nr+3)×N θ (R)
where A spl i ∈ M Nr+3,Nr+1 (R) and A contr ρ j ,i ∈ M Nr+1,Nr+3 (R) for i = 0, ..., N θ -1. These matrices are diagonalisable in the Fourier basis: where

D spl =    D spl 0 . . . D spl N θ -1    D contr ρ j =    D contr ρ j ,0 . . . D contr ρ j ,N θ -1    with D spl m = N θ -1 k=0 A spl k e -2iπkm N θ , D contr ρ j ,m = N θ -1 k=0 A contr ρ j ,k e -2iπkm N θ
and

U n =    U n,0,0 . . . U n,0,N θ -1 . . . . . . . . . U n,N θ -1,0 . . . U n,N θ -1,N θ -1    , U n,k, = 1 √ N θ e 2iπk N θ I n
where I n is the identity matrix of size n × n.

The advantage of the diagonalisation in the Fourier basis is that the computation of the gyroaverage of Φ by the product G G ρ j Φ can be performed using a fast algorithm:

1. Change to Fourier basis by FFT(Φ).

2. Computation of the product G ρ j Φ in the Fourier basis.

3. Change to real space by FFT -1 of the previous result.

The use of polar mesh and FFT allow to make more quickly computations and provides a base for a future work in more complex geometry.

Numerical results

To deal with the equation system (2.6)-(2.7), we have used the SELALIB platform [34] with a classic semi-Lagrangian method with cubic splines interpolation, predictor corrector method and Verlet algorithm for the characteristics. In our case the MPI parallelization is based on transpositions between (r, θ, v) domain decomposition and z domain decomposition.

In the simulations, we will take λ = 0 (no zonal flow case). The initial distribution function reads:

f (0, r, θ, z, v) = f eq (r, v) × 1 + ε exp - (r -r p ) 2 δr cos 2πn L z + mθ
where the equilibrium function f eq is

f eq (r, v) = n 0 (r) exp -v 2 2T i (r) (2πT i (r)) 1/2 .
The profiles T i , T e and n 0 are given by: P(r) = C P exp -κ P δr P tanh r -r p δr P where P ∈ {T i , T e , n 0 }, C T i = C Te = 1 and

C n 0 = r max -r min r min rmax exp -κ n 0 δr n 0 tanh r-rp δrn 0 dr .
We consider the parameters of [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF] [Medium case]:

r min = 0.1, r max = 14.5, v max = 7.32, κ n 0 = 0.055,

κ T i = κ Te = 0.27586, δr T i = δr Te = δr n 0 2 = 1.45, ε = 10 -6 , n = 1, m = 5, L = 1506.759067, r p = r min + r max 2 , δ r = 4δr n 0 δr T i .
We compare the performance of the different numerical methods for the quasineutrality operator on the machines:

• irma-hpc2 with the parameters: 64 × 64 × 32 × 64, ∆t = 5, 1000 iterations.

• Curie with the parameters: 128 × 256 × 32 × 128, ∆t = 2, 4000 iterations.

For the right side of the quasi-neutrality equation, we use the Hermite interpolation method in order to compute the gyroaverage (see [START_REF] Steiner | Gyroaverage operator for a polar mesh[END_REF]) using 1024 points on the Larmor circle with µ = 0.7 on irma-hpc2 and µ = 1 on Curie. The numerical and performance results are given in Table 1, Fig. 2 and Fig. 3. In Table 1, we see that the new Padé approximation is slower than the classic one but faster than the interpolation method. The performance results of the interpolation method seem to be only slightly affected by the number of quadrature points in µ. Indeed, the steps involving the quadrature in µ and the construction of the double gyroaverage matrix are performed once for all in precomputation. The Fig. 2 and3 show that the instability rate obtained with the new Padé approximation is higher than the instability rate obtained with the classical Padé, this is in agreement with Fig. 1. Moreover, the instability rate obtained by the interpolation method is also higher than the instability rate obtained with the classical Padé, this is the opposite behavior of the simple gyroaverage case but there is no compensation in general. The rate of the new Padé method converges to the rate obtained with the interpolation method and we observe that these 2 methods give simular results in Fig. 2 and reveal small structures in the non-linear phase which are not present with the classic Padé method. For memory limitation reasons, we cannot take more than 33 quadrature points in µ on Curie. An outlook would be to make a parallelization in µ in order to increase this number of quadrature points.

Interpolation

Conclusion

We have validated two new methods on polar geometry for solving the quasi-neutrality equation. Comparisons are made with classic Padé approximation, considering a 4D drift-kinetic model with one fixed Larmor radius for the distribution function and an integration in µ for the potential. We find that the instability rate obtained by the new Padé approximation is closer to the rate obtained by the interpolation method (which is accurate when the number of quadrature points increases) than for the classic Padé approximation which is smaller. This is a preliminary study; further study in GYSELA is envisioned.

Figure 1 :

 1 Figure 1: The function 11-Γ 0 (blue curve) compare to its classic Padé approximation 1/x (green curve) and the new approximations with ε = 0.1 (red curve) and ε = 0.001 (cyan curve).

Figure 2 :

 2 Figure 2: Poloidal cut f (r, θ, 0, 0) at time T = 8000 for 128 × 256 × 32 × 128, ∆t = 2, µ = 1. From top to bottom: classic Padé, new Padé (ε = 0.001), interpolation (N µ = 33) . On Curie.

Figure 3 :

 3 Figure 3: Time evolution of rmax r min 2π 0 Φ(r, θ, 0)rdrdθ on irma-hpc2 with µ = 0.7 (top) and on Curie with µ = 1 (bottom).

Table 1 :

 1 Classic Padé New Padé N µ = 9 N µ = 17 N Execution time (in s.) on Curie with 32 processors. Parameters: 128 × 256 × 32 × 128, ∆t = 2, 4000 iterations.

	µ = 33

 [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement N • 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.


Davide Giraudo

INTEGRABILITY CONDITIONS ON COBOUNDARY AND TRANSFER FUNCTION FOR LIMIT THEOREMS

Keywords: Mathematics Subject Classification, 60F17; 60G10 Invariance principle, strictly stationary process, 1

est destinée au dépôt   Copyright Integrability conditions on coboundary and transfer function for limit theorems

Introduction and notations

Let (Ω, F , µ) be a probability space and T : Ω → Ω be a bijective bi-measurable and measure preserving map. We assume that the dynamical system is ergodic (that is, if T -1 A = A for some A ∈ F , then µ(A) ∈ {0, 1}). If n 1 is an integer and f : Ω → R, we denote S n (f ) := n-1 j=0 f • T j and for a fixed t, define (1.1)

S pl n (f, t) := S [nt] (f ) + (nt -[nt])f • T [nt] , t ∈ [0, 1]
, where [x] denote the integer part of the real number x. Then for each ω ∈ Ω and each integer n 1, the map t → S pl n (f, t) is an element of the space of continuous functions on [0, 1], denoted by C[0, 1].

Let us state the limit theorems we are interested in.

Definition 1.1. Let f : Ω → R be a measurable function.

• We say that the function f satisfies the invariance principle if the sequence (n -1/2 S pl n (f, •)) n 1 weakly converges in the space C[0, 1] endowed with the topology of uniform convergence to a scalar multiple of a standard Brownian motion.

• We say that the function f satisfies the law of the iterated logarithm if there exists a constant C(f ) such that for almost every ω ∈ Ω,

(1.2) lim sup n→+∞ S n (f )(ω) √ n log log n = C(f ) and lim inf n→+∞ S n (f )(ω) √ n log log n = -C(f ).
• We say that the function f satisfies the functional law of the iterated logarithm if the sequence ( √ n log log n) -1 S pl n (f, •) n 1 is relatively compact and the set of its limit points coincides with the set of all absolutely continuous functions x ∈ C[0, 1] such that x(0) = 0 and 1 0 (x ′ (t)) 2 dt 1, where x ′ denotes the derivative of x defined almost everywhere with respect to the Lebesgue measure.

• Let 1 < p < 2. We say that the function f satisfies the p-strong law of large numbers if for any α ∈ [1/p, 1], the following holds:

(1.3) ∀ε > 0, +∞ n=1 n αp-2 µ max 1 k n |S k (f )| εn α < +∞.
If it is possible to find a decomposition of the function f

(1.4) f = m + g -g • T,
where g : Ω → R is a measurable function and m satisfies one of the previous limit theorems, then one can wonder if we can deduce the result for f . A known situation is when the sequence (m • T i ) i 0 is a square-integrable martingale differences sequence. A necessary and sufficient condition to have (1.4) with such an m and a square integrable g is known (see Theorem 2 in [START_REF] Volný | Approximating martingales and the central limit theorem for strictly stationary processes[END_REF]). If (m • T i ) i 0 is a squareintegrable martingale difference sequence, then the functional law of the iterated logarithm and the invariance principle take place.

If 1 < p < 2 and m ∈ L p , then Theorem 5 by Dedecker and Merlevède [START_REF] Dedecker | Convergence rates in the law of large numbers for Banachvalued dependent variables[END_REF] implies that m satisfies the p-strong law of large numbers. Actually, their results holds in a more general setting than strictly stationary sequences, as they only require a stochastic domination on the considered martingale differences sequence. A similar result as (1.3) takes place for α = 1 if we require a conditional stochastic domination (see [START_REF] Benoist | Central limit theorem for linear groups[END_REF], Theorem 2.2). A necessary and sufficient condition for (1.4) to hold with m, g ∈ L p , 1 < p < 2, is given by Volný in Theorem 1 of [START_REF]Martingale approximation of non-stationary stochastic processes[END_REF], and in this case, (1.3) is satisfied (see Theorem 6 of [START_REF] Dedecker | Convergence rates in the law of large numbers for Banachvalued dependent variables[END_REF]).

We call a coboundary a function of the form g -g • T , where g : Ω → R is a measurable function. The function g is called a transfer function. The following result is Theorem 1 of [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF]. It gives a necessary and sufficient condition on the transfer function to preserve the limit theorems mentioned in the previous definition. Sufficiency for the invariance principle and the law of the iterated logarithm was established in [START_REF] Hall | Martingale limit theory and its application[END_REF], pages 140-142).

Theorem 1.2 (The equivalence theorem, [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF]). Let us suppose that for the process (m • T i ) i∈Z the invariance principle, the law of the iterated logarithm (functional law of the iterated logarithm) respectively, holds true. Let g be a measurable function and

(1.5) f = m + g -g • T.
Then for the process (f

• T i ) i∈Z
• the invariance principle holds if and only if

(1.6) 1 √ n max 1 k n g • T k → n→∞ 0 in probability;
• the law of the iterated logarithm as well as the functional law of the iterated logarithm holds if and only if

(1.7) 1 √ n log log n g • T n → n→∞ 0 a.s.
Both conditions (1.6) and (1.7) take place when the function g is square-integrable. If 1 < p < 2, Theorem 6 in [START_REF] Dedecker | Convergence rates in the law of large numbers for Banachvalued dependent variables[END_REF] shows that (1.3) holds if g belongs to L p . However, it may happen that we obtain a decomposition (1.4) where m ∈ L 2 but the function g is only integrable (see [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF] for explicit counter-examples, and [START_REF] Esseen | On moment conditions for normed sums of independent variables and martingale differences[END_REF] for a condition which guarantees the square integrability of m) and in this case, the weak invariance principle does not need to hold. We investigate conditions on the functions t → µ {|g| > t} and t → µ {|g -g • T | > t} which guarantee (1.6) or (1.7). In order to state these conditions in a more concise way, we introduce the so-called weak L q -spaces.

Definition 1.3. Let q be a real number strictly greater than 1. We denote by L q,∞ the space of functions h : Ω → R such that

(1.8) h q q,∞ := sup t>0 t q µ {|h| > t} is finite.
The subspace of L q,∞ which consists of functions h such that lim t→+∞ t q µ {|h| > t} = 0 is denoted by L q,∞ 0 .

Main results

In this section, we state the main results of this note. In the first subsection, we give a sufficient condition on the functions t → µ {|g| > t} and t → µ {|g -g • T | > t} in order to preserve the weak invariance principle, the law of the iterated logarithm and the p-strong law of large numbers respectively. We provide projective conditions which guarantee these sufficient conditions an a martingale coboundary decomposition.

In the second subsection, we construct counter-examples which show that the found conditions are sharp when the considered dynamical system is aperiodic.

Finally, in the third subsection, we provide applications of the results of Subsection 2.1 to Bernoulli shifts.

Volný and Samek showed in [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF] that the conclusion of Theorems 2.1 and 2.4 (see the next subsections) holds when p (r + 2)/r and that of Theorem 2.10 when p < (r -1)/(r -3/2). In the case r > 2, we cannot conclude from their results if (r -1)/(r -3/2) p < (r + 2)/r, while it is the case with our conditions. 2.1. Sufficient conditions. Theorem 2.1. Let 1 < p < 2 and let g : Ω → R be a function such that g ∈ L p,∞ 0 and g -g • T ∈ L p/(p-1),∞ 0

. Then for any square integrable martingale differences sequence

(m • T i ) i 0 , the function f := m + g -g • T satisfies the weak invariance principle in C[0, 1].
A similar result has been obtained for the quenched functional central limit theorem (see [START_REF] Barrera | On the functional CLT for stationary Markov Chains started at a point[END_REF], Corollary 7).

Corollary 2.2. Let 1 < p < 2 and let M be a sub-σ-algebra of F such that T M ⊂ M. Assume that f is an M-measurable element of L p/(p-1),∞ 0 such that the following two conditions hold:

(2.1) the sequence (E [S n (f ) | M]) n 1 converges in L p,∞ ;
(2.2) the sequence

((I -U )E [S n (f ) | M]) n 1 converges in L p/(p-1),∞ .
Then the function f satisfies the invariance principle.

The conditions of Corollary 2.2 imply that f admit the martingale-coboundary representation (1.4) with m and g integrable and g satisfies (1.6). In [START_REF] Durieu | Comparison between criteria leading to the weak invariance principle[END_REF][START_REF] Durieu | Independence of four projective criteria for the weak invariance principle[END_REF], the later condition was compared with Dedecker and Rio projective criterion (cf. [START_REF] Dedecker | On the functional central limit theorem for stationary processes[END_REF]):

(2.3) the sequence f E[S n (f ) | M] converges in L 1 ,
which also implies that f satisfies the invariance principle. It was shown that there is an example of function f which satisfies the martingale-coboundary decomposition in L 1 and the invariance principle but not (2.3).

Here do not assume that f belongs to L p/(p-1) and that the convergence in (2.2) holds in L p , otherwise, the function f would satisfy Dedecker-Rio projective criterion. In this way, our condition is independent of (2.3).

Proposition 2.3. Let (Ω, F , µ) be a dynamical system of positive entropy. Let 1 < p < 2. There exists a sub-σ-algebra M such that T M ⊂ M and a function f ∈ L p/(p-1),∞ 0 satisfying (2.1) and (2.2) but not (2.3).

Theorem 2.4. Let 1 < p < 2 < r and let g : Ω → R be a function.

(i) If p > r/(r -1), g ∈ L p,∞ and g -g • T ∈ L r,∞ , then for any martingale differences sequence (m • T i ) i 0 , the function f := m + g -g • T satisfies the law of the iterated logarithm; (ii) if p = r/(r-1), g ∈ L p and g -g •T ∈ L r , then for any square integrable martingale differences sequence (m • T i ) i 0 , the function f := m + g -g • T satisfies the law of the iterated logarithm.

Corollary 2.5. Let 1 < p < 2 and let M be a sub-σ-algebra of F such that T M ⊂ M. Assume that f is an M-measurable element of L p/(p-1) such that the following two conditions hold:

(2.4) the sequence (E [S n (f ) | M]) n 1 converges in L p ;
(2.5) the sequence

((I -U )E [S n (f ) | M]) n 1 converges in L p/(p-1) .
Then the function f satisfies the functional law of the iterated logarithms.

Remark 2.6. In [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF], L q spaces are involved. It turns out that in the setting of Theorems 2.1 and 2.4, (i), we may work with weak L q -spaces. For the case (ii) in Theorem 2.4, it is an open question to determine whether strong moments are actually needed.

Theorem 2.7. Let 1 q < p < r < 2 be real numbers and let g : Ω → R be a function such that g ∈ L q and g -g • T ∈ L r . If q (p -1)r/(r -1), then for any martingale differences sequence (m

• T i ) i 0 with m ∈ L p , the function f := m + g -g • T satisfies (1.3).
Corollary 2.8. Let 1 < p < r < 2 and let M be a sub-σ-algebra of F such that T M ⊂ M. Assume that f is an M-measurable element of L r such that the following two conditions hold:

(2.6) the sequence (E [S n (f ) | M]) n 1 converges in L max{1,(p-1)r/(r-1)} ;
(2.7) the sequence

((I -U )E [S n (f ) | M]) n 1 converges in L r .
Then for each positive ε and any α ∈ [1/p, 1], we have

(2.8) +∞ n=1 n αp-2 µ max 1 k n |S k (f )| εn α < +∞.
Remark 2.9. Ergodicity of the dynamical system is required for the "only if" direction in the equivalence involving the law of the iterated logarithm of Theorem 1.2. Therefore, the results of this subsection remain valid in the non-ergodic setting.

2.2.

Counter-examples. In Theorems 2.1 and 2.4, we gave a sufficient condition on (p, r) for which g ∈ L p and g -g • T r ∈ L r guarantees the invariance principle, and the law of iterated logarithms for g -g • T . The next results show that when this condition is not satisfied, these limit theorems may fail, which shows its sharpness.

Theorem 2.10. Assume that the dynamical system (Ω, F , µ, T ) is aperiodic. Let 1 p < 2 r be real numbers such that p < r/(r -1). Then there exists a function g ∈ L p such that g -g • T ∈ L r and the function g -g • T satisfies neither the invariance principle nor the law of the iterated logarithm.

We also have a similar counter-example for the p-strong law of large numbers.

Theorem 2.11. Assume that the dynamical system (Ω, F , µ, T ) is aperiodic. Let 1 < p < 2 and let 1 q < p < r be real numbers such that q < (p -1)r/(r -1). Then there exists a function g ∈ L q such that g -g • T ∈ L r but the sequence (n -1/p S n (g -g • T )) n 1 does not converge almost surely to 0. In particular, g -g • T does not satisfy (1.3) 2.3. Applications. In Subsection 2.1, we provided sufficient conditions for the functional central limit theorem, the law of the iterated logarithm and the p-law of large numbers. It is natural to try to apply these conditions to some strictly stationary sequences.

In the case of a linear process f

• T k = i 0 a i ε k-i , where (ε i ) i∈Z is an i.i.d. centered sequence of random variables in L p/(p-1),∞ 0
, 1 < p < 2, condition (2.1) implies that f has a martingale coboundary decomposition in L q , 1 < q < p. Since q/2 < 1, the Marcinkievicz-Zygmund and Jensen's inequalities give the convergence of the series i 0

k i a k 1/2
, which in turn implies the martingale-coboundary decomposition in L 2 , from which the invariance principle could have been already deduced. This is why we shall focus on some functionals of a particular linear process: the socalled Bernoulli shifts. Let (ε k ) k∈Z by an i.i.d. sequence of random variables which take the values 0 and 1 with probability 1/2. We define for x

= +∞ k=1 2 -k-1 ε -k (2.9) T n x := +∞ k=1 2 -k-1 ε n-k
The map T preserves the Lebesgue measure λ on the unit interval endowed with the Borel σ-algebra. Given a function f : [0, 1] → R, we would like to give some sufficient conditions on the regularity of f which guarantee the invariance principle, the law of the iterated logarithm and the p-law of large numbers.

Proposition 2.12. Let f : Ω → R be a centered function such that for some p ∈ (1, 2), t p/(p-1) λ {|f | > t} = 0 and for some δ > 0,

(2.10)

1 0 1 0 |f (x) -f (y)| p |x -y| log 1 |x -y| p-1+δ dxdy < ∞ and
(2.11)

1 0 1 0 f (x) -f (y) p p-1 |x -y| log 1 |x -y| 1 p-1 +δ dxdy < ∞, where f (x) = f (x) -f (x/2)/2 -f ((x + 1)/2)/2.
Then f satisfies the invariance principle and the functional law of the iterated logarithm.

Proposition 2.13. Let p ∈ (1, 2). Let f : Ω → R be a centered function such that for some r ∈ (p, 1), f belongs to L r and for some δ > 0,

(2.12)

1 0 1 0 |f (x) -f (y)| q |x -y| log 1 |x -y| q-1+δ dxdy < ∞ and
(2.13)

1 0 1 0 f (x) -f (y) r |x -y| log 1 |x -y| r-1+δ dxdy < ∞,
where q := max {1, (p -1)r/(r -1)} and f (x) = f (x) -f (x/2)/2 -f ((x + 1)/2)/2. Then f satisfies the p-strong law of large numbers.

Proofs

If h : Ω → R is a measurable function, we define M * (h) := sup N 1 N -1 |S N (h)|.
The following lemma about Birkhoff averages will be used in the proof.

Lemma 3.1. Let q > 1 and let h : Ω → R be a measurable function.

(i) If h belongs to L q,∞ 0 then the function M * (h) belongs to L q,∞ 0 ; (ii) if h belongs to L q , then so does M * (h).

Proof. By the maximal ergodic theorem, we have for each positive t,

(3.1) t • µ {M * (h) t} E [|h| • 1 {M * (h) t}] .
(i) The expectation can be bounded by

(3.2) +∞ 0 min µ(A t ), s -q h • 1(A t ) q,∞ ds,
where

(3.3) A t := {M * (h) t} .
Therefore, we infer the bound

E [|h| • 1(A t )] µ(A t ) 1-1/q h1(A t ) q,∞ + (3.4) + h • 1(A t ) q q,∞ +∞ h•1(At) q,∞ µ(At) -1/q s -q ds = µ(A t ) 1-1/q h • 1(A t ) q,∞ + (3.5) + h • 1(A t ) q q,∞ ( h • 1(A t ) q,∞ µ(A t ) -1/q ) 1-q q -1 = µ(A t ) 1-1/q h • 1(A t ) q,∞ + h • 1(A t ) q,∞ µ(A t ) 1-1/q q -1 (3.6) = q q -1 µ(A t ) 1-1/q h • 1(A t ) q,∞ . (3.7)
Plugging this into (3.1), we obtain

(3.8) t • µ {M * (h) t} 1/q q q -1 h • 1(A t ) q,∞ ,
hence is is enough to prove that (3.9) lim t→+∞ sup s 0

s q µ ({|h| s} ∩ A t ) = 0.
To this aim, fix a positive ε; by assumption, there exists a positive real number s 0 such that for s s 0 , we have s q µ {|h| s} ε, hence

(3.10) sup s 0 s q µ ({|h| s} ∩ A t ) max {ε, s p 0 µ(A t )} ,
which is smaller than ε for t large enough. (ii) This follows by multiplying (3.1) by t q-2 , integrating over [0, +∞) with respect to the Lebesgue measure and switching the integrals. This concludes the proof of Lemma 3.1.

We now give the proofs of the main results, which combine Lemma 3.1 with the ideas of [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF].

Proof of sufficient conditions.

Proof of Theorem 2.1. In view of Theorem 1.2, we have to show that the sequence n -1/2 max 1 j n g • T j n 1 converges to 0 in probability.

Let ε be a positive fixed number. Let k, n be positive integers such that k < n. Denoting p n := µ max 1 j n g • T j > 2εn 1/2 , the following estimates take place:

p n µ max 1 i [ n k ]+1 max ik j<(i+1)k g • T ik + g • T j -g • T ik > 2εn 1/2 µ max 1 i [ n k ]+1 g • T ik > ε √ n + [ n k ]+1 i=1 µ max ik j<(i+1)k g • T j -g • T ik > εn 1/2 n k + 1 µ |g| > εn 1/2 + n k + 1 µ max 0 j<k g • T j -g > εn 1/2 n k + 1 µ |g| > εn 1/2 + n k + 1 µ max 0 j<k 1 j |S j (g -g • T )| > ε n 1/2 k .
This yields (3.11)

p n n k + 1 µ |g| > εn 1/2 + n k + 1 µ M * (g -g • T ) > ε n 1/2 k .
By the assumption on p and r, the inequality

(3.12) r/2 -1 r -1 = r -1 -r/2 r -1 = 1 - r 2(r -1) 1 - p 2
takes place, hence we may choose a number α such that

(3.13) 1 - p 2 α r/2 -1 r -1 .
We now use (3.11) with k := [n α ]. This yields, for some constant c depending only on p and r:

(3.14) p n c • n 1-α µ |g| > εn 1/2 + c • n 1-α µ M * (g -g • T ) > εn 1/2-α ,
and, by (3.13),

p n cn p/2 µ |g| > εn 1/2 + (3.15) + cn 1-α-r(1/2-α) n (1/2-α)r µ M * (g -g • T ) > εn 1/2-α = cn p/2 µ |g| > εn 1/2 + (3.16) + cn 1-r/2+α(r-1) n (1/2-α)r µ M * (g -g • T ) > εn 1/2-α cn p/2 µ |g| > εn 1/2 + cn (1/2-α)r µ M * (g -g • T ) > εn 1/2-α . (3.17) Since g ∈ L p,∞ 0 and g -g • T ∈ L r,∞
0 , we conclude by item (i) of Lemma 3.1 that the sequence (p n ) n 1 converges to 0.

This concludes the proof of Theorem 2.1.

Proof of Corollary 2.2. Condition (2.1) implies by [START_REF] Volný | Approximating martingales and the central limit theorem for strictly stationary processes[END_REF] that f may be written as f = m + g -g • T , where (m • T i ) i 0 is a martingale differences sequence with respect to the filtration (T -i M) i 0 . Since f is M-measurable, the function g may be written as i 0 E U i f | T M . Thus, by condition (2.1), we derive that g ∈ L p,∞ and since each term in the series defining g belongs to L p,∞ 0 , we derive that g belongs to L p,∞ 0 . Similarly, by condition (2.2), we infer that g -g • T belongs to L p/(p-1),∞ 0 . Since f ∈ L p/(p-1),∞ 0 , we have m ∈ L p/(p-1),∞ 0 , and accounting the inequality p/(p -1) > 2, we conclude that m is square integrable. The proof is complete since we showed that g satisfies the conditions of Theorem 2.1.

Proof of Proposition 2.3. We use the contruction given in [START_REF] Durieu | Comparison between criteria leading to the weak invariance principle[END_REF]. There exists two independent and T -invariant sub-σ-algebras B and C. We consider a B-measurable function e 0 : Ω → {-1, 1} such that µ({e 0 = 1}) = µ({e 0 = -1}) = 1/2, and define e i := e 0 • T i , i ∈ Z, and M := C ∨ σ {e i , i 0} (which satisfies T M ⊂ M). We introduce three sequences: (θ k ) k 1 ⊂ (0, +∞), (ρ k ) k 1 ⊂ (0, 1) a decreasing sequence such that k 1 ρ k < 1 and an increasing sequence of integers (N k ) k 1 . Once these sequences are fixed, we choose a decreasing sequence

(ε k ) k 1 ⊂ (0, 1) such that (3.18) k 1 θ k N k ε 1/p k .
We now consider for each fixed k 1 a set A k ∈ C such that (1) the sets A k are mutually disjoint;

(2)

1 -i 1 ρ i ρ k 2 µ(A k ) ρ k for all k 1;
(3) for each k 1 and all i, j ∈ {0, . . . ,

N k + 1}, µ(T -i A k ∆T -j A k ) ε k .
The existence of such a sequence of sets as well as that of B and C is explained in [START_REF] Durieu | Comparison between criteria leading to the weak invariance principle[END_REF].

The function f is defined by

(3.19) f = +∞ k=1 θ k e -N k 1(A k ).
Assume that the sequence (θ k ) k 1 is increasing and θ k → +∞. The function f belongs to

L p/(p-1),∞ 0 if (3.20) lim k→+∞ θ p/(p-1) k+1 i k ρ i = 0.
Now, we have

(3.21) E [S n (f ) | M] = +∞ k=1 θ k min{n,N k } i=1 e -N k +i 1(A k )+ +∞ k=1 θ k min{n,N k } i=1 e -N k +i 1(T -i A k \ A k ) -1(A k \ T -i A k ) ,
and since

(3.22) min{n,N k } i=1 e -N k +i 1(T -i A k \ A k ) -1(A k \ T -i A k ) p N k ε 1/p k , the sequence (E [S n (f ) | M]) n 1 converges in L p if (3.23) lim m,n→∞ E +∞ k=1 θ k min{n,N k } i=min{m,N k } e -N k +i 1(A k ) p = 0.
Since the family

(A k ) k 1 is disjoint, (3.23) is equivalent to (3.24) lim m,n→∞ +∞ k=1 θ p k E min{n,N k } i=min{m,N k } e -N k +i p ρ k = 0,
which is implied (by Marcinkievicz-Zygmund inequality) by

(3.25) +∞ k=1 θ p k N p/2 k k .
Notice also that

(I -U )E [S n (f ) | M] = +∞ k=1 θ k (I -U ) min{n,N k } i=1 e -N k +i 1(T -i A k ) = +∞ k=1 θ k e -N k +1 1(T -1 A k ) -e -N k +min{n,N k } 1(T -min{n,N k } A k ) , hence the sequence ((I -U )E [S n (f ) | T M]) n 1 converges in L p/(p-
+∞ k=1 θ 2 k √ N k ρ k < ∞. We thus have to takes sequences (θ k ) k 1 , (ρ k ) k 1 and (N k ) k 1 such that (θ k ) k 1 is increasing, (3.25) and (3.26) hold but +∞ k=1 θ 2 k √ N k ρ k = +∞.
Such a selection is possible; for example, take (3.27)

θ k := 2 k(p-1)/p (log k) 2/p , N k := 2 2k(2-p)/p k 2/p , ρ k := 2 -k .
Proof of Theorem 2.4. In view of Theorem 1.2, we have to prove that the convergence in (1.7) takes place. Let α ∈ (0, 1) be a number which will be specified later. We define (3.28)

m j := j-1 i=1 [i α ] , j 1.
Notice that for some constant κ depending only on α, we have (3.29) j α+1 κ m j κj α+1 , j 1.

By the Borel-Cantelli lemma, we have to prove the convergence of the series

(3.30) +∞ j=1 p j , with p j := µ max 0 i [j α ] 1 m j log log m j g • T mj+i > ε
for each positive ε. To this aim, we start from the inequalities

p j µ 1 m j log log m j |g • T mj | > ε/2 + + µ max 0 i [j α ] 1 m j log log m j g • T mj+i -T mj > ε/2 = µ 1 m j log log m j |g| > ε/2 + + µ max 0 i [j α ] 1 m j log log m j g • T i -g > ε/2 , from which we infer (3.31) p j µ 1 m j log log m j |g| > ε/2 + + µ 1 m j log log m j M * (g -g • T ) > ε 2 [j α ]
.

(i) Assume that p > r/(r -1). Using the definition of • p,∞ and inequality (3.29), we obtain

(3.32) µ 1 m j log log m j |g| > ε/2 c(p, ε, α)κ p/2 g p p,∞ 1 j (α+1)p/2 ,
where c(p, ε, α) is independent of j. Using (3.29) and (3.8), we derive

(3.33) µ 1 m j log log m j M * (g -g • T ) > ε 2 [j α ] c(r, ε, α)κ r/2 g -g • T r r,∞ j αr-(α+1)r/2
, where c(r, ε, α) is independent of j.

Combining (3.31), (3.32) and (3.33), we deduce the upper bound

(3.34) p j c(p, r, α, ε, g) 1 j (α+1)p/2 + 1 j (1-α)r/2 .
We have to take α such that

(3.35) (α + 1)p/2 > 1 and (1 -α)r/2 > 1.
This is equivalent to

(3.36) α > 2/p -1 and α < 1 -2/r.
Since 2/p -1 < 1 -2/r, inequalities (3.35) are satisfied and in view of (3.34), the series defined by (3.30) is convergent for any positive ε. This conclude part (i) of Theorem 2.4. (ii) Assume that p = r/(r -1). We pick α := 2/p -1 = 1 -2/r. In this case, for some constant c depending only on p and r, the inequality

p j µ |g| > cεj 1/p + µ M * (g -g • T ) > εcj -α+(α+1)/2 (3.37) = µ |g| > cεj 1/p + µ M * (g -g • T ) > εcj (1-α)/2 (3.38) = µ |g| > cεj 1/p + µ M * (g -g • T ) > εcj 1/r (3.39)
takes place. By item (ii) of Lemma 3.1, we conclude that the series defined by (3.30) is convergent for any positive ε and this concludes part (ii) of Theorem 2.4, hence the proof of Theorem 2.4.

Proof of Corollary 2.5. Like in the proof of Corollary 2.2, we derive that f admits a martingale-coboundary decomposition with g in L p and m, g -g • T ∈ L p/(p-1) . We thus may apply Theorem 2.4 to conclude that f satisfies the functional law of the iteratd logarithms.

Proof of Theorem 2.7. Let us fix a positive ε and α ∈ [1/p, 1]. Let 1 k < n be integers. By similar inequalities which leaded to (3.11) (we replace the exponent 1/2 by α), we derive (3.40) µ max

1 j n g -g • T j > εn α 2 n k + 1 µ {|g| > εn α /4} + + 2 n k + 1 µ M * (g -g • T ) > εn α 2k .
Let us choose k := n β , where (3.41) (p -q)α β α(r -p)/(r -1)

(the existence of such a β is guaranted by the assumptions on p, q and r). Then it suffices to check that for each positive constant c, the series

(3.42) Σ 1 := +∞ n=1 n pα-1-β µ {|g| cn α } and (3.43) Σ 2 := +∞ n=1 n pα-1-β µ M * (g -g • T ) cn α-β
are convergent. The convergence of Σ 1 is equivalent to the integrability of the function |g| p-β/α ; this holds since (3.41) implies q p -β/α. Note that the second series converges if

(3.44) E (M * (g -g • T )) pα-β α-β < +∞.
Notice that inequality (3.41) implies that (pα -β)/(α -β) r, hence we derive the convergence of S 2 by item (ii) of Lemma 3.1 (with the exponent (pα -β)/(α -β) > 1 since p > 1).

This concludes the proof of Theorem 2.7.

Proof of Corollary 2.8. Like in the proof of Corollaries 2.2 and 2.5, we derive that f admits a martingale-coboundary decomposition with g ∈ L (p-1)r/(r-1) and m, g -g • T ∈ L r . By Theorem 2.7, we derive that for each positive ε and any α ∈ [1/p, 1], (2.8) takes place.

Counter-examples.

Proof of Theorem 2.10. We recall the construction given in the proof of Theorem 3 of [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF]. We choose a real number α such that

(3.45) r -2 2(r -1) < α < 1 - p 2 .
This is possible because

(3.46) 1 - p 2 - r -2 2(r -1) = 1 2 2 -p - r -1 -1 r -1 = = 1 2 1 -p + 1 r -1 = 1 2 r r -1 -p > 0.
For each i 1, we define n i := 2 i and k i := 2 iα . By the Rokhlin lemma (see [START_REF] Kakutani | Induced measure preserving transformations[END_REF][START_REF] Rohlin | A "general" measure-preserving transformation is not mixing[END_REF]), one can find a set A i ∈ F such that (3.47) sets A i , T A i , . . . , T ni-1 A i are pairwise disjoint and

(3.48) µ   ni-1 j=0 T j A i   > 1/2.
In particular, the quantity µ(A i ) can be bounded as follows:

(3.49)

1 2n i µ(A i ) 1 n i .
We then define for i 1,

(3.50)

g i := √ n i log log n i k i   ki j=1 j1(T ni-j A i ) + 2ki-1 j=ki+1 (2k i -j)1(T ni-j A i )   ,
and g := +∞ i=i0 g i , where i 0 is such that 2k i < n i for each i i 0 . By the Borel-Cantelli lemma, since µ {g i = 0} 2k i /n i 2 1-(1-α)i , the series which defines g is almost surely convergent.

Since it has been shown in [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF] that the function f satisfies neither the invariance principle nor the law of the iterated logarithm, it remains to prove that the constructed function g belongs to L p and that the coboundary g -g • T belongs to L r . By (3.47) and (3.50), the equality

(3.51) |g i | p = √ n i log log n i k i p   ki j=1 j p 1(T ni-j A i ) + 2ki-1 j=ki +1 (2k i -j) p 1(T ni-j A i )  
takes place, hence integrating and accounting (3.49), we derive the estimates

E |g i | p √ n i log log n i k i p   ki j=1 j p + 2ki-1 j=ki +1 (2k i -j) p )   1 n i (3.52) 2 n p/2-1 i (log log n i ) p/2 k p i k p+1 i (3.53) = 2n p/2-1 i (log log n i ) p/2 k i , (3.54) hence (3.55) g i p 2 1/p n 1/2-1/p i (log log n i ) 1/2 k 1/p i .
By definition of n i and k i , one can find a constant c depending only on α (hence on p and r) such that for i large enough, (3.56)

g i p c • 2 i(1/2-1/p) (log i) 1/2 • 2 iα/p = c • (log i) 1/2 • 2 i(α-1+p/2)/p .
By (3.45), the series +∞ i=1 (log i) 1/2 • 2 i(α-1+p/2)/p is convergent, and we conclude by (3.56) that g belongs to L p .

It is proved in [START_REF] Volný | On the invariance principle and the law of iterated logarithm for stationary processes[END_REF] that by construction, the equality

(3.57) |g i -g i • T | = √ n i log log n i k i • 1   2ki j=1 T ni-j A i  
holds. By (3.47) and (3.50), we have

(3.58) g i -g i • T r √ n i log log n i k i 2k i n i 1/r
, hence by the definition of n i and k i , we have for i large enough and a constant c depending only on α,

(3.59) g i -g i • T r c • 2 i(1/2-1/r) 2 iα(1/r-1) (log i) 1/2 ,
from which we infer (by (3.45)) the convergence of the series

+∞ i=1 g i -g i • T r hence the fact that the function g -g • T belongs to L r .
This concludes the proof of Theorem 2.10.

Proof of Theorem 2.11. The construction is similar to that of the proof of Theorem 2.10. For each i 1, we define n i := 2 i and k i := 2 iβ , where β satisfies

(3.60) r -p p(r -1) < β < 1 - q p .
Such a choice is possible since

(3.61) p -q - r -p r -1 > p -(p -1) r r -1 - r -p r -1 = p(r -1) -(p -1)r -r + p r -1 = 0.
We take a set A i ∈ F such that (3.47) and (3.48) hold. We then define for i 1,

(3.62)

g i := n 1/p i k i   ki j=1 j1(T ni-j A i ) + 2ki-1 j=ki+1 (2k i -j)1(T ni-j A i )   ,
and g := +∞ i=i0 g i , where i 0 is such that 2k i < n i for each i i 0 . The proof will be complete if we show the following three assertions:

(1) the function g belongs to L q ;

(2) the function g -g • T belongs to L r ;

(3) the sequence (n -1/p g • T n ) n 1 does not converge almost surely to 0.

The first two items follow by completely similar computations as in the proof of Theorem 2.10. To show the last item, we notice that the sequence (2 -i/p max 2 i l 2 i+1 g • T l ) i 1 does not converge to 0 in probability. To see this, one can note that (3.63) µ 2 -i/p max

2 i l 2 i+1 g • T l 1 µ 2 -i/p max 2 i l 2 i+1 g i • T l 1 µ   ni-ki j=1 T j (A i )   1 2 .
Indeed, since g j is non-negative, we have g i g, which gives the first inequality of (3.63).

For the second, notice that since n i = 2 i ,

2 -i/p max 0 l 2 i g i • T l 1 =    max 0 l 2 i ki j=1 j k i 1 T ni-j A i • T l 1    (3.64) ⊃ 2 i l=0 ω | T l ω ∈ T ni-ki A i ⊃ ni-ki j=1 T j (A i ). (3.65)
This finishes the proof of Theorem 2.11.

Applications.

For the context, we refer the reader to Subsection 2.3. We define M := σ(ε i , i 0).

The proof of Propositions 2.12 and 2.13 will follow from Corollaries 2.2, 2.5 and 2.8 and the following intermediate step.

Lemma 3.2. Let q > 1. Then for each centered f : [0, 1] → R and each n 1, the following inequalities hold: < +∞.

(3.66) E [f • T n | M]

1111

  |x -y| 2 -n |f (x) -f (y)| p dxdy (3.67) (I -U )E [f • T n | M] |x -y| 2 -n f (x) -f (y) q dxdy < ∞, where f (x) = f (x) -f (x/2)/2 -f ((x + 1)/2)/2.Proof. Following [MW00], we have for each x ∈ [0, 1],(3.68) E [f • T n | M] (x) = 2 -n |x -y| 2 -n |f (x) -f (y)| q dxdy. Since (3.70) (I -U )E [f • T n | M] (x) = 2 -n 67) in a similar way.Lemma 3.3. Let q > 1. If for some positive δ, we havethen the series n 1 E [f • T n | M] q converges.A similar sufficient condition can be stated for the convergence of the series n 1 (I -U )E [f • T n | M] q .Proof of Lemma 3.3. In view of inequality (3.66), we have to prove that (3.72) |x -y| 2 -n |f (x) -f (y)| q dxdy 1/q

Acknowledgements. The author would like to thank anonymous referees for their comments, which led to an improvement of the presentation of the paper.

To this aim, we define β := (q -1 + δ)/q and bound, by Hölder's inequality,

1-1/q < +∞ and since βq/(q -1) = (q -1 + δ)/(q -1) > 1, it suffices to show that the series