Salwa Kobeissi
email: salwa.kobeissi@inria.fr

Philippe Clauss
email: philippe.clauss@inria.fr

The Polyhedral Model Beyond Loops Recursion Optimization and Parallelization Through Polyhedral Modeling

Keywords: Compilation, dynamic analysis, loop optimization, parallelization, polyhedral model, recursion, recursive function ACM Reference

There may be a huge gap between the statements outlined by programmers in a program source code and instructions that are actually performed by a given processor architecture when running the executable code. This gap is due to the way the input code has been interpreted, translated and transformed by the compiler and the final processor hardware. Thus, there is an opportunity for efficient optimization strategies, that are dedicated to specific control structures and memory access patterns, to apply as soon as the actual runtime behavior has been discovered, even if they could not have been applied on the original source code.

In this paper, we develop this idea by identifying code extracts that behave as polyhedral-compliant loops at runtime, while not having been outlined at all as loops in the original source code. In particular, we are interested in recursive functions whose runtime behavior can be modeled as polyhedral loops. Therefore, the scope of this study exclusively includes recursive functions whose control flow and memory accesses exhibit an affine behavior, which means that there exists a semantically equivalent affine loop nest, candidate for polyhedral optimizations. Accordingly, our approach is based on analyzing early executions of a recursive program using a Nested Loop Recognition (NLR) algorithm, performing the affine loop modeling of the original program runtime behavior, which is then used to generate an equivalent iterative program, finally optimized using the polyhedral compiler Polly. We present some preliminary results showing that this approach brings recursion optimization techniques into a higher level in addition to widening the scope of the polyhedral model to include originally non-loop programs.

Introduction

In the area of automatic program optimization and compilation for imperative languages, loop structures are obviously legitimate targets, since they represent a large and computeintensive part of the whole execution time. Accordingly, numerous advanced loop analysis and optimization techniques have been developed and implemented in compilers and other tools dedicated to optimizing code transformations. Many of these sophisticated loop optimizers are based on the polyhedral model from which they acquire their powerful capabilities in this domain (e.g., Polly [START_REF] Grosser | Polly -performing polyhedral optimizations on a low-level intermediate representation[END_REF], Pluto [START_REF] Bondhugula | A practical automatic polyhedral parallelizer and locality optimizer[END_REF]).

The polyhedral model is a mathematical framework that contributes a substantial abstraction and representation for programs having affine loop nests and accessing multi-dimensional arrays through affine array references, i.e., programs with static control parts (SCoPs). This framework provides powerful analysis, and aggressive loop automatic optimizing and parallelizing transformations (e.g., loop tiling, loop skewing, loop interchange, etc. ).

The Apollo framework [START_REF] Martinez Caamano | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF][START_REF] Sukumaran-Rajam | The polyhedral model of nonlinear loops[END_REF] shows that even if loops that are outlined in a source code do not have the required shape, i.e., are not SCoPs, they may still be good candidates for polyhedral optimizations since their actual runtime behavior is similar to SCoPs' runtime behavior: actual loop bounds and memory references can be modeled as affine functions of surrounding loop indices. In this paper, we go further with this idea of discovering a "SCoP" runtime behavior for code extracts that are even not outlined as loops in the source code, but as recursive functions. Compared to what has been achieved with Apollo, the general goal is not only to build an affine modeling of the data flow, but also to build an affine (loop) modeling of the control flow.

Among the time-expensive structures, recursive functions also play an important role. They generally implement complex algorithms that may scan huge data structures such as graphs, trees and matrices. A recursive approach is usually adopted when solving any problem whose solution relies on smaller and simpler instances of its own. It also eases the expression of computations that are generic to parameters as search depths or problem dimensions.

It is well-known that recursive and iterative (loops) approaches can be used alternately to solve a problem. Replacing a recursive code by an equivalent iterative code and vice versa is usually possible [START_REF] Arsac | Some techniques for recursion removal from recursive functions[END_REF]. In particular, tail recursive functions are automatically handled by production compilers (CLANG, GCC) and transformed into equivalent loops. However, the so-generated loops are generally not suitable for advanced polyhedral optimizations. Moreover, unlike loops, recursive functions do not take advantage of very advanced automatic parallelization and optimization techniques [START_REF] Gupta | Automatic parallelization of recursive procedures[END_REF][START_REF] Mastoras | Ariadne -directive-based parallelism extraction from recursive functions[END_REF].

In this paper, we present our promising idea to make use of the polyhedral model's powerful capabilities to optimize and parallelize recursions, through a guided modeling of their runtime behaviors as polyhedral-compliant loops. Note that the proposed approach is different than static recursion elimination, since it mostly ignores the syntactic shape of the target code. It should rather be understood as a dynamic and automatic rewrite of the target recursive code, through the polyhedral modeling of a part of its data and control flow. Our approach mainly consists of a dynamic analysis technique to check for a linear looping memory and control behavior for a given recursion. Thus, before replacing the original recursive code, we make sure that, at runtime, it acts as an affine for-loop. Subsequently, we generate a semantically equivalent iterative code and, finally, apply further polyhedral optimizations.

The paper is organized as follows : Section 2 discusses some representative previous works on recursion optimizations and parallelizations. Section 3 describes the steps of our method and the tools utilized to implement our idea, and Section 4 supports our proposal by showing some preliminary results. Finally, Section 5 sums up the work and proposes some future works concerning our purpose.

Related Work

Various studies exist for the purpose of optimizing and parallelizing recursive functions. Nevertheless, the proposed parallelization techniques are mainly static and involve task parallelization where several recursive calls or invocations are run simultaneously.

Rugina and Rinard [START_REF] Rugina | Automatic parallelization of divide and conquer algorithms[END_REF] present a compiler to parallelize divide-and-conquer recursive functions using pointer analysis and symbolic analysis to detect independent recursive calls and generate code that executes these calls concurrently. Also, for divide-and-conquer implementations, Gupta et al. [START_REF] Gupta | Automatic parallelization of recursive procedures[END_REF] propose a compile-time framework that uses interprocedural symbolic array section analysis to detect the independence of multiple recursive calls of divide-and-conquer algorithms. It also proposes a speculative run-time parallelization technique when static analysis is not sufficient. Besides, a tool called Huckleberry [START_REF] Collins | Recursion-driven parallel code generation for multi-core platforms[END_REF] automatically parallelizes sequential recursive divide-and-conquer implementations for multi-core platforms.

Another approach [START_REF] Morihata | Automatic parallelization of recursive functions using quantifier elimination[END_REF] proposed by Morihata and Matsuzaki is an automatic parallelization method for recursive functions using quantifier elimination, where the input structure is partitioned into chunks that are run in parallel. In addition, Mizutani et al. [START_REF] Mizutani | Evaluation of a compiler with user-selectable execution strategies for parallel recursion[END_REF] propose a different strategy to parallelize recursive functions where programmers decide on using simple or dynamic load balancing depending on the workload of the recursive functions and on the threshold up to which each call is executed in parallel. Saougkos et al. [START_REF] Saougkos | Fine grained parallelism in recursive function calls[END_REF] propose an automatic fine-grained parallelism extraction method for recursive functions having integer variables updated in a systematic way.

DECAF [START_REF] Gupta | Optimizing recursive task parallel programs[END_REF] is a technique to optimize recursive task parallel programs by reducing the task creation and termination overheads. Tetzlaff and Glesner [START_REF] Tetzlaff | Static prediction of recursion frequency using machine learning to enable hot spot optimizations[END_REF] propose a machine learning based approach to statically predict the recursion frequency of functions used to guide various hot spot optimizations.

Sundararajah et al. [START_REF] Sundararajah | Locality transformations for nested recursive iteration spaces[END_REF] propose recursion twisting transformations for nested recursions that improve data locality. However, their technique is solely devoted to specific shapes of recursions where data is organized in two trees, called inner and outer trees, and where recursive calls are nested.

A different approach for optimizing recursive implementations is recursion elimination. In this context, the work presented by Bird [START_REF] Bird | Notes on recursion elimination[END_REF] proposes recursive elimination methods for some forms of recursions using stacks. Cohen [START_REF] Cohen | Characterization and elimination of redundancy in recursive programs[END_REF] presents four types of redundancy in recursive functions and suggests a solution for each of these types by eliminating recursion without using a stack. Another recursion elimination technique is tackled by Bird [START_REF] Bird | Tabulation techniques for recursive programs[END_REF], so-called tabulation, which addresses repetitive (redundant) evaluation of certain values by only computing values for once and, then, storing them for later uses. Also, Liu and Stoller [START_REF] Liu | Dynamic programming via static incrementalization[END_REF] describe a method for transforming some recursive programs into programs using the dynamic programming design technique. Stitt and Villarreal [START_REF] Stitt | Recursion flattening[END_REF] propose recursion flattening that can eliminate many instances of recursion by determining recursion depth, and then inlining recursive calls. However, it cannot eliminate all recursions.

Adriadne [START_REF] Mastoras | Ariadne -directive-based parallelism extraction from recursive functions[END_REF] is a compiler that extracts directive-based parallelism from recursive function calls. It extracts three forms of parallelisms and a transformation for each one: (1) recursion elimination: converting recursive function to iterative one, (2) parallel-reduction eliminating recursion and distributing workload into independent tasks, (3) thread-safe parallelizing recursive functions that contain independent recursive calls. So, Adriadne for one of its approaches performs recursion-to-iteration transformation, but it only supports recursive functions whose parameters remain unchanged among recursive calls, except for one integer parameter, the so-called index variable, participating in all the conditions and the termination of the recursive function.

Furthermore, some compilers today (e.g. clang/LLVM [START_REF] Lattner | LLVM: A Compilation Framework for Lifelong Program Analysis & Transformation[END_REF]) implement an optimization pass for transforming recursive functions into loops, however limited to tail recursions, i.e., to recursive functions where the recursive call is the last instruction of the function. Additionally, this transformation is not incremental and, thus, is not able to transform several nested recursive calls which would become tail after its application. Other well-known recursion elimination and recursion-to-loop transformation techniques are not yet automated in these compilers.

To sum up, recursive implementations may take advantage of parallelism and other optimizations as described in the studies above. However, the finest parallelism granularity is usually one recursive function invocation where data locality optimization cannot be handled very accurately. When transformations of recursive functions into loops are proposed, the resulting loops are generally too complicated for taking advantage of efficient compilers loop optimizations, thus not providing significant performance improvement regarding their recursive counterpart. This is mostly due to the employed static approaches that try to build a generic equivalent loop program, where either the call stack is simulated using a dedicated data structure, or the generated loops are while-loops with dynamic termination conditions, or the loop bodies embed many conditional branches to mimic the original cases of recursive calls.

In this paper, we show that dynamic analysis of recursive functions and guided modeling of the runtime actual behavior enable the generation of simpler equivalent loops that may be affine, i.e. polyhedral-compliant.

From Recursive Functions to Optimized Loops

In this section, we present a prototype of the implementation of our approach as a proof of concept for modeling recursive functions as affine loops and applying polyhedral optimizations.

Our framework involves the steps schematized in Figure 1 which can be summarized as three main phases: (1) Recursive control and memory behavior analysis, (2) Recursion to affine loop modeling and (3) Code generation and polyhedral optimizations. In the input recursive program, we distinguish two kinds of instructions:

1. Instructions whose role is exclusively dedicated to direct the control flow. In the case of recursive functions, such instructions are used to decide about function invocations and the setting of some parameters related to further potential (recursive) calls. 2. Instructions that participate in the computation and memory store of the final result of the program. Such instructions act directly or indirectly for the updating of the data structure that stores the final result. They constitute the minimal set of instructions of the original recursive program that, if conveniently scheduled and instantiated, have to be run to reach the (correct) Our framework is aimed at handling recursive programs whose control and memory runtime behavior can be characterized as follows:

• Basic blocks of the LLVM intermediate representation that contain impacting instructions are scheduled at runtime as if they would belong to nested loops whose respective unique index has affine loop bounds depending on the surrounding loop indices. • The target addresses of memory accesses regarding the main input and output data structures can be modeled as affine functions of surrounding loop indices; • The general affine loop structure modeling the control and memory behavior must be independent of the processed input. By "input-independent", we mean that the number and depths of the loop nests remains the same regardless of the input, while only some loop bounds and some coefficients of the memory reference functions may change accordingly.

We describe in details the steps of our framework in the following subsections. As illustrating example, we use a recursive implementation of the matrix product shown in Listing 1.

1 v o i d M a t r i x M u l t i p l i c a t i o n ( i n t A[N ] [ N ] , i n t B [N ] [ N ] ) 2 { 3
s t a t i c i n t row = 0 , column = 0 , i n d e x = 0 ;

Recursive Control and Memory Behavior

Analysis In this first phase, we execute the original recursive program given a relatively small input. If the recursion's behavior is affine, then the analysis result of the execution can be used to build a structure of equivalent loop nests. However, we may need to analyze several executions of the program to predict some loop bounds with respect to the input. The analysis results are then used for extrapolating the program's behavior for larger inputs. This phase involves several steps that are detailed in the following subsections.

Compile-Time Classical Optimizations

First of all, given an input recursive source code, we let the LLVM compiler apply some classical optimizations (e.g., promote memory to register) to prepare the intermediate representation (IR) of our input program for the later steps. We do not activate the tail call elimination pass which transforms tail recursive calls into loops for two reasons: (1) the way the target recursive function is transformed may not result in an affine loop and (2) if there are several recursive calls in the target code, only one tail call may be eliminated. a tail call. However, since the tail call elimination pass is not applied incrementally by the compiler, these remaining calls are not eliminated. Thus, we do not activate this pass in order to let our framework handle all the three recursive calls in the same way.

Compile-Time Recursion Detection

At this step, our framework checks, at compile time, if the program involves any recursions, and, if so, identifies these recursions and the functions participating in them. In order to detect recursions, it uses the call graph corresponding to the IR of the program. The call graph is a directed graph representing relationships between functions in a program, where a node represents a function and an edge from one node to another represents a call such that the source node constitutes the caller function and the destination is the callee. Figure 2 shows an example of a call graph of a program made up of seven functions: main, f , д, h, i, j, k such that:

• function main calls f • f invokes itself, h and д • д calls back f , k and i • i calls j
In this example, we notice that function f is in a direct recursion with itself and an indirect recursion through function д.

After obtaining the call graph, we search for strongly connected components (SCC) in it such that an SCC is a subgraph where every node is reachable from every other node. In this context, a cycle in SCC implies that a recursion exists among the functions represented by the nodes involved in this cycle. If the cycle is made up of only one node (loop), then it is a direct recursion. Else, it is indirect. In our example of Figure 2, we have an SCC involving the two gray nodes of functions f and д such that we have a loop over f signifying a direct recursion and a cycle from f to д and д back to f , implying an indirect recursion among f and д.

Recursion Reachability Recognition

Eventually, we need to trace the control and memory behavior of recursions and find out if they act as static control parts. Thus, we are interested in tracking basic blocks containing impacting instructions, whether they are executed directly or indirectly by the recursive functions. For this reason, in addition to the recursive functions themselves, our tool also determines their reachability in the program. By reachability, we mean all the functions that can be reached by a sequence of calls initiated by the recursive functions themselves. In Figure 2, the reachability of the recursive functions f and д includes:

• h (directly called by f )

• k and i (directly called by д)

• j (indirectly called by д through i) Example 3.2. In the matrix product example, function mul-tiplyMatrix is the only recursive function, and it does not call any other function than itself.

Impacting Instructions

We identify impacting instructions in the following way. In the IR of the program, our framework marks the writes to the data structure which is the final output of the recursion and their corresponding reachable functions. Then, for every write instruction of this kind, it marks all instructions leading and contributing to it, i.e., their backward static slices (BSS). A backward static slice, is a set of instructions existing in the code of a program that may affect a certain value, i.e., in our case, a value of the data structure which is the final output of the recursion. Slicing the important parts of a given recursive program (i.e., recursive functions and their reachability functions) is essential for the analysis phase of our implementation since it helps keeping track of the elemental parts of recursions that must be taken into consideration for further transformations and optimizations.

Example 3.3. In the LLVM IR of the recursive matrix product, the memory writes related to matrix C in function multi-plyMatrix, as well as their backward static slices are marked. These backward static slices include the memory reads related to matrices A, B and C, and the computations related to the updates of i, j and k, since they are involved in the computations of the memory addresses referenced by the memory instructions.

Profiling

Finally, printing statements are added to the IR to collect the IDs of executed basic blocks that include impacting instructions, and that would be enough to run to reach the final result if they are well initiated. Other blocks that are only dedicated to the recursive control are not instrumented. Additionally, writes to memory and reads from memory are to be instrumented to collect the corresponding target memory addresses, with the final purpose of building affine reference functions.

The program executes for a first time, and the output expected is the memory and control trace of recursions in the program. Subsequently, we investigate if this trace can be modeled as affine loop nests. For this purpose, we use the Nested Loop Recognition (NLR) algorithm originally presented in [START_REF] Ketterlin | Prediction and trace compression of data access addresses through nested loop recognition[END_REF].

Nested Loop Recognition (NLR) Algorithm

The NLR algorithm takes as input a trace of a program execution and constructs a sequence of loop nests that produce the same original trace when run. The applications of this algorithm include: (1) program behavior modeling for any measured quantity such as memory accesses, (2) execution trace compressing and (3) value prediction, i.e, extrapolating loops under construction (while reading input) to predict incoming values.

In our tool, not only do we use NLR to model memory accesses, which is one of its original goals, but also to model sequences of basic blocks IDs, which is more singular. Given our trace of a target recursive program, generated thanks to our instrumentation, if NLR builds affine loop nests including the interesting basic blocks IDs and memory addresses interpolated by the constructed loop indices, then the generation of equivalent affine loops may be performed.

An initial NLR model is built from running the instrumented code with small data input. Then, additional runs with different data inputs are performed and their associated NLR models are compared to the initial one, in order to find out if (1) the general affine loop structures are maintained across the models and (2) if the differences between the models can be identified as being solely loop bounds or coefficients of memory reference functions. Then, their respective values are interpolated relatively to some input parameters, as typically the problem size. Constant parts of the affine functions modeling memory references may obviously vary across the models, since they mostly represent base addresses of data structures.

Example 3.4. Listing 2 offers a shortened view of the result of NLR. The runtime behavior of function multiplyMatrix is modeled by NLR as two affine loop nests including memory accesses that are expressed as affine combinations of the loop indices. The first loop nest is obviously dedicated to the main computations of the matrix product, while the second loop nest models the returns of the recursive calls. Note that loop bounds and some memory reference coefficients are set to their generic value related to the problem size N , thanks to polynomial interpolation resulting from comparing the modeling of several runs with different problem sizes. Items MEM1, MEM2 and MEM3 represent the base addresses of matrices A, B and C respectively. v a l M a t r i x M u l t i p l i c a t i o n : : i f . end15 28 . . .

Listing 2. Recursive Matrix Multiplication NLR model

Recursion to Affine Loop Transformation

This phase requires both the NLR model and the original recursive program itself. Code generation and transformation are achieved as follows. According to the loop nests provided by NLR, our framework extracts their structures and generates corresponding loop nests in LLVM form using the suitable basic blocks, taking into consideration changing the accessed memory addresses to be in terms of the loop indices as it is in the NLR model. Then, they are inserted into the LLVM IR of the program replacing the targeted recursions.

Polyhedral Optimizations Application

The final resulting affine loop nests replacing the original recursive functions may already take advantage of standard compiler loop optimizations (loop invariant code motion, constant propagation, common subexpression elimination, etc.). However, they are also potential candidates for polyhedral optimizing and parallelizing transformations. For this purpose, our framework applies Polly [START_REF]Polly -llvm framework for high-level loop and data-locality optimizations[END_REF] on the LLVM IR of the affine loop nests.

Experiments

In this section, we show the results of our experiments by applying the described approach to the matrix product that was used as illustration all along the paper. The computer used is an Intel(R) Core(TM) i7-7Y75 CPU @ 1.30GHz running Ubuntu 16.04.5 LTS. We use Polly with llvm/clang version 8.0.0.

We compare the execution of the original recursive code against the iterative code, either handled solely by clang with flag -O3 or also handled by Polly. Our measures are represented by the curves in Figure 3 using a logarithmic scale for the execution times, and by varying the matrix sizes from 100 to 1,000. It clearly shows that the iterative code outperforms the recursive one, and particularly when polyhedral optimizations have been applied by Polly.

Conclusion and Perspectives

We have shown that recursive functions may act at runtime as affine loops, and hence be also good candidates for polyhedral optimizations. A dedicated automatic recursionto-affine-loop transformation made of static and dynamic analysis and transformation passes could be applied whenever possible, on any program embedding recursive functions. Such an approach brings handled recursive functions to a higher level of optimizations, compared to what has been achieved until now regarding recursive functions in compiler research. It also extends the polyhedral model applicability to non-loop control structures.

Yet, our work is still a proof of concept for this approach, and we aim at investigating more in this promising idea to implement an efficient specialized compiler for this sake covering many more and complicated recursive programs. Some of the future works, that should be done in order to achieve this purpose, include:

1. Performing dynamic analysis for recursive behavior at runtime on-the-fly, instead of analyzing an early program execution with a smaller input. This is because some recursive functions might not completely behave as affine loop nests. Instead they might involve non continuous large static control parts during their execution. Thus, at runtime, we keep checking for linear behavior and, accordingly, transform the code for executing these certain parts. 2. Inducing verification features to the framework so dynamic analysis can be possible and the iterative model constructed can be considered predictive. 3. Tackling input dependent recursive codes.

Figure 1 .

 1 Figure 1. General View of the Framework

Listing 1 .

 1 ] [ column ]+=A[ row ] [ i n d e x ] * B [ i n d e x ] [ column ] ; 11 i n d e x + + ; 12 M a t r i x M u l t i p l i c a t i o n ( A , B ) ; r i x M u l t i p l i c a t i o n ( A , B ) ; r i x M u l t i p l i c a t i o n ( A , B ) ; 21 } Recursive Matrix Multiplication C Function

Example 3 . 1 .Figure 2 .

 312 Figure 2. Example of a Call Graph of an Arbitrary Recursive Program

9 l o a d 10 v

 910 M a t r i x M u l t i p l i c a t i o n : : i f . t h e n 4 / / I R b a s i c b l o c MEM1 + 4 * N * i 0 + 4 * i 2 / / memory a d d r e s s i n t e r m s o f l o o p s i n d i c e s 8 . . . / / r e p e t i t i v e memory a c c e s s p a t t e r n s a l MEM2 + 4 * i 1 + 4 * N * i 2 / / 4 i s t h e s i z e o f an i n t e g e

17

 17

  v a l M a t r i x M u l t i p l i c a t i o n : : i f . end15 18 . . .19 v a l M a t r i x M u l t i p l i c a t i o n : : i f . end17 20 . . . 21 f o r i 0 = 0 t o N * N-1 22 f o r i 1 = 0 t o N-1 23 v a l M a t r i x M u l t i p l i c a t i o n : : i f . end17 . . . 25 v a l M a t r i x M u l t i p l i c a t i o n : : i f . end15 26 . . .

  27

Figure 3 .

 3 Figure 3. Performance results with the matrix product


V Bondarenko
email: valeria_bondarenko@yahoo.com

Member #94362473, IEEE S K Filiova-Petrakieva

D A Andreev

Keywords: stochastic model, fractional Brownian motion (fBm), parameters estimation, forecasting power consumption I

In the present paper will be discussed the problem related to the individual household electric power consumption of objects in different areas -industry, farmers, banks, hospitals, theaters, hostels, supermarkets, universities. The main goal of the directed research is to estimate the active P and full S power consumptions for all studied buildings. The defined goal is achieved by solving of the following three problems. The first problem studies which buildings increase their power consumption. The second one finds which objects have the greatest increase of power consumption. And the third problem regards if it is possible to make a short-term forecast, based on the solutions of previous two problems. The present research and solving of the aforementioned problems is conducted using fractional Brownian motion theory. The applicability of this approach is illustrated on the example with 20 real objects in different areas. The paper ends with conclusion notes about possibilities to make short-term forecasts about power consumption of the considered buildings.

approaches for forecasting: stochastic, intervals, based on artificial intelligent: cluster analysis, neural networks etc. To make a more adequate forecast, on one side it is necessary to have initial data as huge as possible. On the other side, it makes longer the processing time and it requires more powerful computers.

These forecasts could be very useful for different applications such as optimizing the active and the full power energy distribution as a function of the consumption.

There are two types of forecasts: short-time and long-time ones. It depends on the type of the considered object, respectively the associated forecasting parameter. In particular the reliability of the forecast depends on the type of the considered time series and on the amount of the primary data for the analyzed index.

The paper is organized as follows. In the next section are defined three basic problems for analysis of real data about the consumed active power P and the full power S for 20 different objects. The basis of fractional Brownian motion (fBm) theory is presented in Section 3. In Section 4 is discussed the algorithm for analyzing the considered time series. The real examples for the consumed active power P and the full power S are discussed in Section 5. The paper ends with concluding remarks about power consummation and its future forecasting for the considered objects.

II. PROBLEM STATEMENT

Let us consider the power consumption (active power P and full power S) in different administrative buildings where the consumption has different daily and nightly behavior. The analyzed objects correspond to different life's areas such as: 'Bank', 'Automobile Industry', 'BPO Industry', 'Cement Industry', 'Farmers 1', 'Farmers 2', 'Health Care Resources', 'Textile Industry', 'Poultry Industry', 'Residential (individual)', 'Residential (apartments)', 'Food Industry' ,'Chemical Industry', 'Handlooms', 'Fertilizer Industry', 'Hostel', 'Hospital', 'Supermarket', 'Theatre' and 'University'.

The main aim of the research is to estimate the actual power consumption (P and S) of each analyzed object and to make forecast about their future consumptions. This aim is realized by solving the following problems: Problem 1: Which buildings increase their consumption (P and S)?

Forecasting time series for power consumption data in different buildings using the fractional Brownian motion

Valeria Bondarenko, Simona Filipova-Petrakieva, Ina Taralova, Desislav Andreev T Problem 2: Which buildings have the greatest increase of their consumption (P and S)?

Problem 3: If it is possible to make a short-terms forecast, based on the results from the solutions of Problems

k = = 0 ( 1) , 0, 1 1
T t t k t k n - = + - = ( ( ))( ) x X t k k = Φ • is a model of the observed time series { } , , ..., 1 2 x x x n .
Let us call ξ as a basic process of the model. One of the most popular Markov models of time series is Gaussian random process, and fractional Brownian motion (fBm). The large use of this process is caused by its "convenient" properties, which are described below.

Fractional Brownian motion is defined as a Gaussian random process with characteristics:

( ) ( ) { } ( ) ( ) { } ( ) ( ) , 0, 0 0 1 2 2 2 . 2 B t
( ) ( ) ( ) ( ) 1 2 2 2 2 4 1 3 2 4 2 3 1 2 H H H H t t t t t t t t   - + - - - --     For discrete time: 1 k k B B k H H n n ξ -     = -        
, we obtain the correlation coefficient:

1 2 2 2 ( , ) ( 1 1 2 ), 2 H H H k j k j k j j k ρ ξ ξ = -+ + -- - - (2)
It means that increments form a stationary (in the narrow sense) sequence.

Let us mention some properties of fBm: 1. Changing time scale is equivalent to changing the "amplitude" of the process:

( ) ( ) ( ) ( ) H Law B at Law a B t H H = .
This equality denotes the coincidence of one-dimensional distributions of the processes:

( ) ( ) and H B at a B t H H .
This property is called self-similarity and it is useful for the analysis of time series.

Let us put in the formulae

(2) j = k + n . Then the correlation coefficient: ( ) ( ) 1 2 2 2 ( , ) 1 1 2 2 H H H n n n k k n ρ ξ ξ   = + + - - +     (3)
When n → ∞ it follows:

2 2 ( , ) ( 2
1) H H H n n k k n ρ ρ ξ ξ - = - + 
(4) So, the increments memory decrease exhibits a power law; the increments are independent for 0.5 H = . With 0.5 H < the increments form the sequence with short memory, and for 0.5 H > a sequence with a long memory. The sign of the correlation coefficient n ρ , defined by formula (4), depends on the value of H: 0, 0.5

H n ρ < < ; 0, 0.5 H n ρ > > .
For 0.5 H < the sequence n ξ of increments fBm is called "pink noise" and negativity of the variations indicates fast variability values. The process of fBm, 0.5 H < is known as "antipersistent". For 0.5 H > the sequence y n of increments fBm is called "black noise" and the fBm process is known as "persistent".

The original approach using stochastic and statistical analysis is based on real time series modeling assuming a hypothesis that their increments can be modeled by fBm. If the estimated Hurst exponent is bigger than 0.5 the process is called persistent and a forecast can be made.

IV. THE SELECTION OF THE CONVERSION OF TIME SERIES

( ) , 1,..., 1 1 y x x k n k k k = - = - +
and construct the statistics

(kurtosis) 2 1 . 1 2 1 2 . 1 2 n y k n k d n n y k n k       - =   = - = ∑ ∑ ,
z g y k k = .
Let us assume

( ) ( ) ( ) ( ) ( ) 2 1 1 . 1 lim 2 1 1 1 2 . 1 1 n g y k n k n n g y k n k π -       - =   = - → ∞ - = ∑ ∑
, where ( )

z g y k k =
is assumed as a Gaussian random value.

Let us demonstrate proposed algorithm with g as a power function. Assume ) ( sgn

1 1 k k k k y y y z - Φ = = λ , (5) sgn
( ), 0 y y z z k k k k λ λ = = Φ > , ( 5a
) then 2 1 1 . 1 1 1 1 2 . 1 1 n z k n k d n n z k n k λ λ -       - =   = - - = ∑ ∑ (6) lim d d n n = → ∞ is equal to the ratio of the corresponding mathematical expectations. If ( ) 2 2 1 2 0, . . 2 E α α α α ξ σ ξ σ π +   Ξ ⇒ = Γ      , then 1 2 1 2 . 1 2 d λ π λ +   Γ     =   Γ +    
, where the parameter λ is defined from the equations (5, 5a).

Thus, the proposed approximation leads to the following model of original time series:

sgn . 1 k x y z k j j j λ   =     = ∑ .
If we assume that the values of the sequence { } , , ..., 1 2 z z z n are increments of fBm, the Hurst exponent can be calculated by the following algorithm, as it can be seen from the proposed method:

1) Construct the statistic:

1 1 1 n R z z k n k = = = ∑ 2) Calculate the matrix 1 S H
-, where S H is a correlation matrix of increments fBm :

1 1 . ( , ) k k j j s B B B B jk H H H H j k n n n n ρ ξ ξ  -  -         = Ε - - =                         3) ( ) 1 , 0,8 1 1 S z z H Q R n - = - ,
The statistics Q is calculated for different Hurst exponent values with step (0.05 -0.1) and the corresponding Ĥ is calculated, such as: .

( ) 1 min, arg min ( ) 1. Q H H Q H -→ = - (6
k k H z c c n B B k k n n ξ  +      = = -             . Assume 1 1 k v z k j j - = =
∑ and construct the statistics [5, 6] ( ) ( ) ( )

1 3 . , 0; 0.5 1 1 2 3 . , 0; 0.5 1 1 1 3 . , 0.5; 1 2 1 k k k k n A v z H n k n k n B v z H n H n k n D v z H n k H n k = ∈ = = ∈ + = = ∈ = ∑ ∑ ∑ (7) If hypothesis T is true, there is convergence [1]: 3 3 2 2,5 2 2 ; 3 ; 2 2 A c B c D c B n n n η → - → → .
The decision about the veracity of the hypothesis T is accepted by comparing the real values of the statistics with their theoretical limit values. Therefore, let us determine the deviation from the limit values for the statistic A n

A A n A δ - = .
The limit distribution functions for the statistics ,

B D n n : { } 2.5 2.5 ( ) 3 ; 1 3 1 2 ( ) 2 1, 0, 2 3 x F x P c x c F x x x c η σ   - = < = Φ       = Φ - >      
, where Φ is the Laplace function, 0.5 ( 22) .

H σ - = +
The hypothesis T is accepted, if:

, , 0.5 0 1 or 0 , 0.5 2 B H n D H n δ β β β < < < < < > (

V. REAL EXAMPLES

Let us consider the consumed active power P and full power S for the following 20 objects. The data for individual household electric power consumption are saved each hour [START_REF]Individual household electric power consumption[END_REF]. The considered buildings are: 'Bank', 'Automobile Industry', 'BPO Industry', 'Cement Industry', 'Farmers 1', 'Farmers 2', 'Health Care Resources', 'Textile Industry', 'Poultry Industry', 'Residential (individual)', 'Residential (apartments)', 'Food Industry' ,'Chemical Industry', 'Handlooms', 'Fertilizer Industry', 'Hostel', 'Hospital', 'Supermarket', 'Theatre' and 'University'.

The research aim is to estimate the power consumption (P and S) of each analyzed object. This aim is realized by solving the following problems: Problem 1: Which buildings increase their consumption (P and S)? Problem 2: Which buildings have the greatest increase in their consumption (P and S)?

Problem 3: If it is possible to make a short forecasting, based on the results from the solutions of Problems 1 and 2?

Before starting the data processing all primary data have been normalized and the respective values of P and S take values in interval [0, 1]. Next, the primary data have been processed using linear approximation of the trend of the initial data and new sequences have been obtained in every example. This is realized by the following steps: 1. Calculation of the increments 5) and (5a).

( ) , 1,..., 1 1 y x x k n k k k = - = - + . 2. Construction of the new sequence { } , , ..., 1 2 1 z z z n- by (
3. Hurst exponent estimation by [START_REF] Nourdin | Central and non-central limit the orems for weighted power variations of fractional Brownian motion[END_REF]. [START_REF] Barndorff-Nielsen | I. A tutorial on Levy processes[END_REF]. Checking the quality of the approximation by A n and B n [START_REF]Individual household electric power consumption[END_REF].

The processing results are shown in Tables IA andIB.

As it can be seen from the TABLE IA, for all the data we have antipersistent parameter ˆ0.5 H < . For the examples (except Theater and Food Industry), the condition is clearly verified:

1 B n β < .
For all examples (except Theater and Textile Industry), the condition is clearly verified:

1 B n β <
and the process has short-memory.

According to the results gotten in the previous section it can make the following conclusions about three problems defined in the beginning of the paper: Problem 1: Which buildings increase their consumption (P and S)?

The approximation by fractional Brownian motion of the considered real data shows that Bank, Automobile Industry, BPO Industry, Textile Industry, Residential (individual), Chemical Industry, Hospital and University increase their consumption of the active power P, because for these buildings we've got the highest antipersistent Hurst parameter ˆ0.4 H = (which means that the processing data change more smoothly).

Problem 2: Which buildings have the greatest increase in their consumption (P and S)?

The approximation by fractional Brownian motion of the considered real data shows that Bank, BPO Industry and Textile Industry increase their full power S consumption because for these buildings we've got the highest antipersistent Hurst parameters: ˆ0.4 H = , ˆ0.46 H = and ˆ0.46 H = , respectively (which means that the processing data change more smoothly).

Problem 3: Is it possible to make a short forecasting, based on the results from the solutions of Problems 1 and 2?

All of the approximations have antipersistent character ( ˆ0.5 H < ) and it is adequate, if the conditions (5) are satisfied. Therefore, it is impossible to construct a short-term forecast. Then other methods will be used for approximation of the analyzed real data to make more reliable prognoses. But it will be the focus of the future investigations, where some newer techniques such as cluster analysis will be applied.

Numerical experiment has shown that for a lot of data it is possible to implement approximation by fractional Brownian motion, but without forecast and with short memory.

For the examples Food Industry, Textile Industry and Theater conditions [START_REF] Nourdin | Density formula and concentration inequalities with Malliavin calculus[END_REF] are not satisfied, so it means that they can not be assumed as Gaussian processes and it is possible to approximate them by other distributions.

VI. CONCLUSION

Fractional Brownian motion is used for approximation of real object data to estimate a capability of short and long-term prognosis about an active power P and full power S consumption. Different buildings consumptions have been analyzed, and the relevance of the chosen modeling by fBm increments has been evaluated using limit theorems. The used methodology is appropriated for short-term prognosis for all considered objects except of the Theater and Textile Industry.

ACKNOWLEDGMENT

The present survey has been carried out as a result of the collaboration between Bulgarian and French teams working together during the joint project "Cluster analysis of aperiodic time series signals" financed by the Bilateral program PHC "RILA" projects 2017-2018 -Bulgaria -France. Mr. Desislav Andreev has excellent communication and leading skills gained through his work in Visteon as a lecturer, team lead and a software engineer. His job requires multiple business travels to Santa Clara, Detroit and Karlsruhe, working together with clients and partners of Visteon. He is an advanced developer in C and C++ programing languages.

) 4 )

 4 Verify the following hypothesis T: The statistics of real data can be simulated by fractional Brownian motion increments. The algorithm with known H is the followingassume that hypothesis Т is verified 1 .

  methodology. She is currently a Researcher in the Laboratoire des Sciences du Numérique de Nantes (LS2N), which is a Joint Research Unit (U.M.R. 6004). Her scientific interests include control engineering, nonlinear system dynamics, pseudo-random chaotic generators design, secure information transmission. She is manager and scientific advisor of national and international research projects, and reviewer for International Journal ofBifurcation and Chaos, Control System Letters, Control Engineering Practice; IMA Journal of Mathematical Control and Information etc. Desislav Andreev was born in 1991 in Dobrich, Bulgaria. He received his B.S. in Computer Systems and Technologies at the Faculty of Computer Systems and Control of Technical University of Sofia. Next he received his second M.S. in Computer Systems and Technologies at the Faculty of German Engineering and Economics Education of Technical University of Sofia in 2016. He was top of the graduated. From 2012 to 2014 he works as a Software Developer in Atia Soft Ltd. and in the meantime starts as an Assistant Professor in the Faculty of Computer Systems and Control. I parallel he developed software for video editing tools, data monitoring and media database administration. From 2014 till now he work in Visteon Electronics as a Software Engineer. His main area of expertise is the system and infotainment states of automotive ECUs and he is one of the co-founders of the Visteon Engineering Academy. He engineered the architecture of the low-level software, communication layers, middle layer software and user interface protocols, thus assuming in 2016 the position of Senior Software Engineer. In 2016 he started as a Ph.D. student in Artificial Intelligence Systems at the Technical University of Sofia, preparing dissertation about Machine Learning Algorithms in Quantum Entanglement-based Environment. In 2017 he transferred his position in the ADAS department of Visteon, working on Autonomous Driving Technologies and becoming in 2018 lead of the Sofia branch of the department. He still works as a Assistant Professor. His main research interests are in Parallel Computing, Pattern Recognition, Sensor Fusion, Quantum Physics and Quantum Field Theory, Artificial Environments and World Simulations.

This work has been carried out as a result of the collaboration between Bulgarian and French teams working


Monique Aouad

Guillaume Borrel

Céline Brochier-Armanet
email: celine.brochier-armanet@univ-lyon1.fr

Simonetta Gribaldo
email: simonetta.gribaldo@pasteur.fr

Evolutionary placement of Methanonatronarchaeia 8

Methanonatronarchaeia, a newly discovered archaeal lineage of extremely halophilic

19

96

Correspondence and request for materials should be addressed to Simonetta Gribaldo 97 (simonetta.gribaldo@pasteur.fr) or Céline Brochier-Armanet (celine.brochier-armanet@univ-

  methanogens, were proposed to represent an evolutionary intermediate between 20 archaeal methanogens and the extremely halophilic Halobacteria. Here, we show that 21 the sistership between Methanonatronarchaeia and Halobacteria results from a tree 22 reconstruction artefact and that the divergence of Methanonatronarchaeia is in fact 23 much deeper. This sheds a new light on the adaptation to extreme halophilic lifestyle 24 in archaea and on the evolution of methanogenesis. 25 26 Sorokin and colleagues recently reported the identification of Methanonatronarchaeia, a 27 fascinating archaeal lineage of extremely halophilic, moderately thermophilic, methyl-28 reducing methanogens 1,2 . Similar to most recently discovered methanogens, 29 Methanonatronarchaeia perform methanogenesis based on H 2 and methyl compounds, a 30 metabolism not previously reported from hypersaline environments. Together with 31 Halobacteria and Nanohaloarchaea 3 , Methanonatronarchaeia represent the third discovered 32 lineage of extreme halophilic archaea and the most halophilic methanogens ever found. They33 have likely adapted to this lifestyle by employing a salt-in osmoprotection strategy 1 , unlike 34 previously known halophilic methanogens and similarly to the two other extreme halophilic 35 archaeal lineages 4 . Moreover, Methanonatronarchaeia rely on cytochromes for 36 methanogenesis 1 , a characteristic previously thought to be restricted to the 37 Methanosarcinales 5 . A maximum Likelihood (ML) phylogenetic analysis of a supermatrix 38 gathering ribosomal proteins indicated Methanonatronarchaeia as the closest relatives to 39 Halobacteria (Fig. 1A, red branches) 1 . They were therefore proposed to be evolutionary 40 intermediates on the path from methanogens to extreme halophiles 1 . However, multiple 41 substitutions occurring at the same site in sequences can mask the original phylogenetic 42 signal and provoke tree reconstruction artefacts 6 , a phenomenon particularly evident in 43 lineages that adapted to extreme salinity 7 . 44 supermatrix of ribosomal proteins used by Sorokin et al. 1 , through the progressive removal of 47 the fastest evolving sites, a method that is frequently used to reduce artefacts linked to 48 multiple substitutions 7 . This analysis, both by ML and Bayesian approaches including non-49 homogeneous evolutionary models, shows that the clustering of Halobacteria and 50 Methanonatronarchaeia (Fig. 1B-C, red line) was recovered only when the fastest evolving-51 sites are included in the analysis, while the progressive removal of these sites shifted the 52 position of Methanonatronarchaeia away from Halobacteria and to a deeper branching 53 position at the base of the superclass 'Methanotecta' 8 (Fig. 1B-C, green line). This placement 54 is also consistently and robustly recovered when Methanonatronarchaeia were included in 55 two recently published supermatrices comprising a larger number of markers 6 (over 250 56 conserved protein families) or a larger taxonomic sampling of the Methanotecta 9 (including 57 ANME1, Syntrophoarchaeales, Methanoliparia, and a third Methanonatronarchaeia 58 member). In contrast with the dataset of Sorokin et al. 1 , the new placement of 59 Methanonatronarchaeia was robust to the removal of the fastest-evolving sites for both these 60 supermatrices (Fig. 1D-G). 61 62 Our analyses indicate that the placement of Methanonatronarchaeia as the methanogenic 63 closest relatives of Halobacteria proposed in Sorokin et al. 1 is likely the consequence of a 64 tree reconstruction artefact induced by a multiple substitution-bias which is particulary strong 65 in their ribosomal protein dataset, but not in the other two datasets. The alternative position 66 of the Methanonatronarchaeia disclosed here provides a new perspective on the evolution of 67 this fascinating lineage. For example, it indicates that their adaptation to extreme halophily 68 would have occurred independently from the Halobacteria. Moreover, following the recent 69 proposal for the placement of Nanohaloarchaea as sister to the Methanocellales 6 , the salt-in 70 strategy used for thriving in hypersaline environments would have emerged three times 71 independently in the Archaea, a remarkable example of convergent evolution for adaptation 72 to similar environments. Finally, the new placement of Methanonatronarchaeia is highly 73 relevant for the evolution and diversity of methanogenesis, as their characteristics may 74 reflect those of the methanogenic ancestor of the whole 'Methanotecta' superclass. For 75 example, the fact that Methanonatronarchaeia rely on cytochromes for methanogenesis 3 76 raises the question of whether this feature may be ancestral to all Class II methanogens and 77 was retained only in Methanosarcinales while Methanomicrobiales, Methanocellales and 78 Methanoflorentaceae shifted secondarily to methanogenesis without cytochromes, or if 79 instead it emerged twice independently. 80 The current pace in the acquisition of genomic data and the discovery of new lineages 8,10 will 81 certainly allow to tackle these fundamental questions in the evolution and ecology of 82 methanogens and of Archaea in general. et al. (2017) Nat Microbiol 2:17081.86 2) Sorokin et al. (2018) Int J Syst Evol Microbiol 68(7):2199-2208.

87 3 )

 3 Narasingarao et al. (2012) ISME J 6, 81-93 88 4) Oren (2008) Saline Syst 4, 2.

89 5 )

 5 Thauer et al. (2008) Nat Rev Microbiol 6(8):579-91.

90 6 )

 6 Delsuc (2005) Nat Rev Genet 6(5):361-75.

91 7 )

 7 Aouad et al. (2018) Mol Phylogenet Evol 127:46-54 92 8) Adam et al. (2017) ISME J (11):2407-2425.

93 9 )

 9 Borrel et al. Nat Microbiol (accepted).

94 10 )

 10 Spang et al. (2017) Science 357(6351).

Acknowledgements 101 S.G., G.B., and C.B.A acknowledge funding from the French National Agency for Research,

102

Grant ArchEvol (ANR-16-CE02-0005-01). M.A. is funded by a doctoral fellowship from the 103 Région Rhône-Alpes-ARC 1 Santé. We thank the PRABI (Pôle Rhône-Alpes de 104 Bioinformatique) and the computational and storage services (TARS cluster) provided by the