Ioan D Landau
email: ioan-dore.landau@gipsa-lab.fr

Jaime Saavedra

Sophie Cerf

Bogdan Robu

Nicolas Marchand

Sara Bouchenak

Can adaptive feedforward control improve operation of cloud services?

Keywords: Control of computing systems, Adaptive control, Cloud computing, Feedback control, Feedforward control, MapReduce

. However since the disturbance is measurable an adaptive feedforward compensation can be added for improving the performance of a feedback controller. A general algorithm for adaptive feedforward control in the context of cloud computing is proposed and analyzed. Simplified versions are also presented and analyzed. Experimental results obtained on the Grid'5000 (French nation-wide cluster infrastructure) will be presented.

I. INTRODUCTION

Thanks to the improved capacity of the hardware storage and transport network, the cloud computing architecture has evolved. Cloud computing is particularly well suited for the processing of Big Data1 . Cloud computing is becoming a more and more attractive solution. One of the biggest appeals of cloud computing is the on-demand assigning of a large group of shared hardware resources to software applications, called elastic resource provisioning. Therefore, we need frameworks to efficiently deal with the dynamic aspect of cluster resources allocation [START_REF] Lorido-Botran | A review of auto-scaling techniques for elastic applications in cloud environments[END_REF].

One of issues with automatic resource provisioning approaches such as the ones currently deployed in public clouds is that they don't work well for real time applications. Currently, if an application running in the cloud has to meet runtime criteria, it is up to the human application manager to decide the amount of resources it needs. However, when the manager is notified that an application is running slowly, it requires high level of expertise to decide on how much to intervene. This decision is a highly difficult task due to workloads fluctuations over time [START_REF] Ari | Managing flash crowds on the internet[END_REF], [START_REF] Li | A realistic evaluation of memory hardware errors and software system susceptibility[END_REF], the use of shared hardware resources or interferency and concurrency issues among many others [START_REF] Koh | An analysis of performance interference effects in virtual environments[END_REF]. Therefore, there is a need for fully automatized cluster scaling algorithms. Furthermore, as cloud providers desire to maximize the resource utilization, they have mechanisms for the dynamic reallocation of unused resources in the cluster, which further adds to the variability of system performance. So even with the same workload and resource amount, an application performance depends also on the load induced by neighboring applications, making the ideal scaling even harder to achieve. Current approaches to ensure performance in cloud systems are basically static taking in account experience. Improvement of this approach is provided by a predictive estimation of the computer load. Another approach called reactive can be considered as using feedback but with no guarantee that the system will be stable Of course combination of these approaches (hybrid) is also considered [START_REF] Ali-Eldin | An adaptive hybrid elasticity controller for cloud infrastructures[END_REF].

The use of a control approach for driving the computing resources allocation has been also considered in the recent years. This new application of control brings novel challenges that enable to enrich both communities. A basic framework is considered in the book [START_REF] Hellerstein | Feedback control of computing systems[END_REF] however it is the paper [START_REF] Berekmeri | A Control Approach for Performance of Big Data Systems[END_REF] which represents a significant advance in applying control methodology to this problem since in addition of the theoretical analysis, experimental results obtained on the Grid'5000, a French nation wide infrastructure made up of a 5000 CPU's are presented. The concepts of control objectives (like service time), disturbances (user demand) and the control input (number of CPU) have been clarified. Models of the system have been identified. A PI control has been implemented for assuring steady state performances and a feedforward control has been considered since the disturbance is measurable. Unfortunately the characteristics of the model vary and the feedforward control is very sensitive to the knowledge of the gain of the system. Introduction of an adaptive feedforward control able to provide good results when the values of the parameters of the system are unknown has been considered in [START_REF] Cerf | Adaptive feedforward and feedback control for cloud services[END_REF]. The objective was to adapt only the static gain of the disturbance compensator in the absence of the feedback controller. Analysis of the algorithm has been provided as well as simulation results using the Berekmeri's identified models and a PI feedback controller.

In the present paper we present a general algorithm for adaptive feedforward control of cloud computing centers taking in account the full dynamic of the system and the presence of a feedback controller. Furthermore since the models of the system have delays, which can not be handled optimally by a PI controller, one considers feedback polynomial controllers with an integrator. This allows to assign the poles of the closed loop at desired values even in the presence of large delays. This approach enables more robustness compared with the state of the art solutions. The paper is organized as follows: Sections II, III and IV give a control perspective to the system considered. Sections V and VI give details on the control structure. Sections VII and VIII are dedicated to the development and analysis of adaptive feedforward control algorithms. Experimental results are presented in Section IX.

II. CLOUD COMPUTING OPERATION -A CONTROL

PERSPECTIVE

The use of the control approach for resource allocations in cloud computing can not ignore the computer structure and the programming techniques used. We have to be able do define the variables of interest, measure them, transmit them to the controller and implement the control. We are in the context of massive parallel processing over distributed platforms and one of the most popular programming approaches is MapReduce and its implementation via the Hadoop environment developed by Apache [START_REF] Dean | Mapreduce: simplified data processing on large clusters[END_REF].

From the cloud provider point of view, the service time (the time it takes for a user request to be treated) is thus a performance metric that has to be monitored (controlled variable i.e. the system output). One simple action that can be done to control the on-line service time is to modify the number of resources of the cloud allocated to the jobs (control input). For the MapReduce use case, adding resources, commonly known as nodes, to the cluster will increase the number of Map and Reduce functions processing the input data leading to a reduction of the service time. If the number of resources (nodes) diminishes, the results are converse. The cluster size is then our control signal. However, in the case of public clouds, multiple clients send requests at the same time (measurable disturbance) thus generating a varying input workload which influences the service performance. If multiple concurrent jobs are running, the amount of resources allocated for each job is reduced and thus the job service time increases. As the workload of the system is independent of the cloud provider we will consider it as a disturbance. Figure 1 gives an image of the system from a control perspective. This approach will be illustrated at the level of a cluster2 .

III. MAPREDUCE ENVIRONMENT

As an example of application of our methodology, we chose the programming model called MapReduce.

For a user to run a MapReduce job, at least three things need to be supplied to the framework: the input data to be treated, a Map function, and a Reduce function. From our point of view, the Map and Reduce functions can be only treated as black box models since they are entirely application-specific, and we have no a priori knowledge of their behavior. Without some profiling, no assumptions can be made regarding their runtime, resource usage or the amount of output data produced. On top of this, many factors (independent of the input data and of the Map and Reduce functions) are identified that influence the performance of MapReduce jobs: CPU, input/output skews, software failures [10], Hadoop's node homogeneity assumption not holding up [START_REF] Anjos | MRA++: Scheduling and data placement on mapreduce for heterogeneous environments[END_REF], and bursty workloads [START_REF] Ghit | Balanced resource allocations across multiple dynamic MapReduce clusters[END_REF] among others. All the sources of disturbances specific to clouds (concurrency, network skews [START_REF] Li | OFScheduler: A dynamic network optimizer for MapReduce in heterogeneous cluster[END_REF] or hardware failures among others) can be added to the variability of our system, especially when MapReduce is executed on a public cloud.

MapReduce resource provisioning for ensuring Service Level Agreement (SLA) 3 objectives is relatively a fresh area of research. Few approaches formulate the problem of finding the optimal resource configuration either as an offline optimization problem [START_REF] Zhang | Automated profiling and resource management of pig programs for meeting service level objectives[END_REF] or as an on-line one, based on past runtimes [START_REF] Papadopoulos | Peas: A performance evaluation framework for autoscaling strategies in cloud applications[END_REF], or by ensuring job level deadlines. Our work is related to these approaches but with several important differences, since we deal on-line with a workload of multiple data intensive jobs (which is the case in large production clusters), while taking advantage of the mathematical proofs of convergence and robustness of control theoretic techniques. The control is implemented in Matlab and all the measurements are made on-line in real time. The service time (the time it takes for a client request to be fulfilled) and the number of clients are measured from the cluster. The number of nodes in the cluster is a control signal that we use to ensure the service time deadlines, regardless the changes in the number of clients. All the actuators and sensors are implemented in Linux Bash scripts.

IV. MAPREDUCE MODELING

The experiments presented in this paper, were run using the MapReduce Benchmark Suite (MRBS) developed by [START_REF] Sangroya | MRBS: towards dependability benchmarking for Hadoop Mapreduce[END_REF], which is a performance and dependability benchmark suite for MapReduce systems. MRBS can emulate several types of workloads and inject different fault types into a MapReduce system. The workloads emulated by MRBS are selected to represent a range of loads, from the compute-intensive to the data-intensive (e.g. business intelligence -BI) workload.

One of the main objectives of MRBS is to emulate client interactions, which may consist of one or more MapReduce jobs run at the same time. A data intensive BI workload is selected as the workload. The BI benchmark consists of a decision support system for a wholesale supplier. Each client interaction emulates a typical business oriented query run over a large amount of data (10GB here). To generate the client interactions Apache Hive is deployed on top of Hadoop. This converts SQL like queries to a series of MapReduce jobs. All the nodes in the cluster are on the same switch to minimize network skews. All the experiments have been conducted on-line, on Grid'5000, on a single cluster of 60 nodes. A dynamic model that predicts MapReduce cluster performance (i.e. the average service time in the present paper) with respect to the number of nodes and clients was proposed and validated by experiments in [START_REF] Berekmeri | A Control Approach for Performance of Big Data Systems[END_REF]. Since the system is considered linear in the operating region, the principle of superposition can be applied to calculate the output. The model structure is recalled here in Figure 2, and the mathematical expression in the discrete time domain is:

+ + u(k) y (k) d(k
y(t) = D(q -1 )d(t) + G(q -1 )u(t) (1)
where y is the average service time of all jobs (the output), u is the changes in the number of nodes in the cluster (the control input) and d is the changes in clients (considered as a measurable disturbance). In the following we will use the error ε(t), which is the comparison of the measured service time with its reference value, as a performance indicator of our control. D(q -1 ) is the direct path transfer function i.e. the link between the disturbance and our performance metric, and G(q -1 ) represents the compensatory path, allowing to modify the cluster size for controlling the service time. Both models were identified as first-order transfer function with delays, as described in eqs. ( 2) and ( 3)

D(q -1 ) = b D q -1 (1 + β d q -1 ) 1 + a D q -1 q -r D = B d A d (2)
G(q -1 ) = b g q -1 (1 + β g q -1 ) 1 + a G q -1 q -r G = B g A g (3)
The terms (1 + β d q -1 ) and (1 + β g q -1 ) allow to take in account fractional delay (less than one sampling period).

V. FEEDBACK CONTROL

A feedback control is introduced to reduce the effect of the disturbances (number of clients) and to assure in steady state the desired value of the service time. As such, the feedback controller should incorporate an integrator. Since the system to be controlled has delays, it is reasonable to use a polynomial controller R/S in order to be able to assign conveniently the poles of the closed loop and to shape the sensitivity functions for robustness and performance. The structure of the controller is:

K(q -1 ) = B K (q -1 ) (1 -q -1 )A K (q -1 ) , (4)
where

B K (q -1 ) = b K 0 + b K 1 q -1 + ... + b K n B K q -n B K , (5)
A K (q -1 ) = 1 + a K 1 q -1 + ... + a K n A K q -n A K . (6)
For the models considered, the closed loop poles are the roots of the polynomial equation:

P(q -1 ) = A g (1 -q -1 )A K + B g B K (7)
This equation allows either to compute B K and A K given the desired closed loop poles defined by P or to compute the closed loop poles given B K and A K . A particular case of the polynomial R/S controller is the digital PI controller already used in [START_REF] Berekmeri | Feedback Autonomic Provisioning for Guaranteeing Performance in MapReduce Systems[END_REF]:

H PI = b K 0 + b K 1 q -1 (1 -q -1 ) (8)
The closed loop poles are in this case given by:

P PI (q -1 ) = (1 + a g q -1 )(1 -q -1 ) (9)
+b g (1 + β g q -1 )(b K 0 + b K 1 q -1 )q -(r g +1) VI. STRUCTURE OF THE FEEDFORWARD CONTROL Since the disturbance (number of clients) is measurable, a feedforward compensation can be added. In a linear context one has to assume that D and G in Figure 2 are perfectly known. This will allow to define the structure of the feedforward compensator and to compute its optimal parameters. Unfortunately in practice since the parameters of G and D are not perfectly known and in addition they can be time varying an adaptive approach has to be considered. For the linear case with known parameters one has the following result: Lemma 1: Assuming that the parameters of D and G are known it exists a feedforward compensator:

N(q -1 ) = q -r r 0 + r 1 q -1 + r 2 q -2 1 + s 1 q -1 + s 2 q -2 (10)
where: r = r dr g [START_REF] Anjos | MRA++: Scheduling and data placement on mapreduce for heterogeneous environments[END_REF] such that (perfect matching condition):

N.G = -D (12)
In this case perfect compensation of the disturbance is achieved.

Proof: Taking N as:

N = (1 + a g q -1 )(1 + β g q -1 ) (1 + a d q -1 )(1 + β d q -1 ) b d b g q -r ( 13
)
it can be straightforwardly verified that eq. 12 holds. Doing the corresponding multiplications in eq. 13 one gets the form of N given in eq. 10 Therefore N given in eq. 10 has the structure of an optimal feedforward compensator and the optimal values can be computed from eq. 13 . Since the parameters of the direct path as well as the parameter b g of the secondary path are not perfectly known and they also can vary, an adaptive feedforward compensator should be used. A global view of the control system is shown in Fig 3. The adaptive feedforward compensator will have the same structure as the linear compensator except that its parameters will be adapted. The corresponding structure is:

N(q -1 ) = q -d N (t, q -1 ) (14)
where:

N (t, q -1 ) = r0 (t) + r1 (t)q -1 + r2 (t)q -2 1 + ŝ1 (t)q -1 + ŝ2 (t)q -2 (15)
ri (t) and ŝi (t) will be estimated in real time. The input to the feedforward compensator is d(t), the output of the feedforward compensator is denoted u 1 (t) and the output of the block N is denoted u 1 (t).

The development of the adaptive control law will be done under the following hypotheses: (H5) There exists a filter N of the form given in eq. 14 such that eq. 12 holds. (H6) The system operates in the presence of a feedback polynomial controller of the form given in eq. 4

To proceed to the development of the adaptive feedforward control law, the key point is to establish a relation between the estimated parameters of the feedforward compensator and the measured error in terms of service time error with respect to the desired value. Using the results of [START_REF] Landau | Adaptive and Robust Active Vibration Control: Methodology and Tests[END_REF], Chapter 15, Section 15.17, eq. ( 15.86) for our particular scheme one gets:

ε(t + 1) = A k B g S(A g A k + BgB k ) [θ -θ ] T φ (t) (16)
where S is the denominator of the optimal linear feedforward compensator,

θ T = [s 1 , s 2 , r 0 , r 1 , r 2 ] (17)
is the vector of the parameters of the optimal filter N (unknown),

θ T = [ ŝ1 , ŝ2 , r0 , r1 , r2 ] (18)
is the vector of estimated parameters of the filter N ,

φ T (t) = [-u 1 (t), -u 1 (t -1), d(t + 1), d(t), d(t -1)] (19)
is the observation vector, where: d(t + 1) 4 , and u (t) are respectively the input and the output of the filter N .

It is useful for stability reasons (as it will be shown later) to filter the observation vector φ (t):

φ f (t) = L(q -1 )φ (t) (20)
Eq.16 will take the form:

ε(t + 1) = A k B g S(A g A k + BgB k )L [θ -θ ] T φ f (t) (21)
This equation when θ is time-varying (under the effect of adaptation) takes the form (neglecting the non-commutativity of the time operators when θ is time varying):

ε(t + 1) = A k B g S(A g A k + BgB k )L [θ -θ (t + 1)] T φ f (t) (22)
This equation has the standard form of an a posteriori adaptation error equation [START_REF] Landau | Adaptive control: algorithms, analysis and applications[END_REF] allowing straightforwardly to write the parameter adaptation algorithm and to get the stability condition of the full scheme. One has the following result: Lemma 2:(stability of the adaptive control law) For the system described by eqs. 2, 3, 4 and 14 using the parameter adaptation algorithm:

θ (t + 1) = θ (t) + FΦ(t)ε(t + 1) ; (23)
ε(t + 1) = ε • (t + 1) 1 + Φ T (t)FΦ(t) ; (24)
F = αI ; α > 0 ; (25)
Φ(t) = φ f (t) ; ( 26
)
where ε • (t + 1) is the measured a priori error on the service time and under the condition that:

H = A k B g S(A g A k + BgB k )L (27)
is a strictly positive real transfer function, one has:

lim t→∞ ε(t + 1) = lim t→∞ ε • (t + 1) = 0. ( 28
)
for any bounded initial conditions. Proof: The proof is a direct application of the results presented in [START_REF] Landau | Adaptive control: algorithms, analysis and applications[END_REF], chapter 3 and it is omitted. Remark: A time varying matrix adaptation gain can also be used (see [START_REF] Landau | Adaptive control: algorithms, analysis and applications[END_REF] for details) It follows that for the implementation of this adaptation algorithm the filter L should be designed in order to satisfy the strictly positive real condition on 27. Taking

L = A k Bg Ŝ( Âg A k + Bg B k ) (29)
The condition of eq. 27 becomes that the transfer function:

H L = B g Bg Ŝ S ( Âg A k + Bg B k ) (A g A k + BgB k ) (30)
should be strictly positive real. The satisfaction of this condition will depend upon the quality of the estimation of Ag, Bg, S. Conversely the allowed tolerance in the estimation of Ag, Bg, S results from the positive real condition on the transfer function of eq. VII.

VIII. AN INTERESTING PARTICULAR CASE

One can assume that the feedback controller has a very slow action and that its main task is to assure a zero steady state error. In this case one may neglect its dynamic effect and therefore one can consider A K = 1 and B K = 0 (absence of the controller). This lead to a filter L of the form:

L = Bg Ŝ Âg (31)
Furthermore we can make the assumptions :

r g = r d ; a g = a d ; β g = β d (32)
which means in practice that we look to the compensation of the disturbance in steady state only. In this case the feedforward commpensation law will take the form:

u(t) = gd(t) (33)
which corresponds to R = g = b d /b g and S = 1. The correponding adaptive feedforward law becomes:

u(t) = ĝ(t)d(t) ( 34
)
where ĝ is given by ĝ

(t + 1) = ĝ(t) + α f (t + 1)e(t + 1) = ĝ(t) + α f (t + 1) 1 + α f 2 (t + 1)
e 0 (t + 1) with α > 0 (35) where e 0 (t + 1) is the a priori measured service time error and f (t + 1) is given by:

f (t + 1) = L(q -1 )d(t) (36)
where:

L = (sign b g )q -(r G +1) (1 + βg q -1 ) 1 + âg q -1 , (37)
The stability condition becomes in this case:

H OL = 1 + âg z -1 1 + a g z -1 . 1 + βg z -1 1 + β g z -1 (38)
should be a strictly positive real transfer function. This implies first that βg < 1 for the stability of the transfer function. If âg is enough close to a g and βg is enough close to β g , this condition is always satisfied. 5 It is this algorithm which has been implemented in the presence of a digital PI controller.

IX. EXPERIMENTAL RESULTS

The experiments were conducted on Grid'5000, on a single cluster of 60 nodes. Each node from the cluster used for the test has a quad-core Intel CPU of 2.53GHz, an internal RAM memory of 15GB, 298GB disk space and infinite band network. The interconnection between the cloud and the controller implemented in Matlab has required an intensive effort for developing the appropriate interfaces and communication protocols [START_REF] Berekmeri | A Control Approach for Performance of Big Data Systems[END_REF]. The workload used for the experiments are some Business Intelligence task, taken from MRBS benchmark suite [START_REF] Sangroya | MRBS: towards dependability benchmarking for Hadoop Mapreduce[END_REF]. Only the number of tasks to be executed is varied here (#clients). Fig. 4 gives an block diagram of the full system. 6 shows the evolution of the service time and of the number of nodes when a digital PI controller is used. Fig. 7 shows the evolution of the service time and of the number of nodes when in addition to the digital PI controller an adaptive feedforward control law is used. One can see that adding adaptive feedforward compensation improves the performance obtained with just a PI controller. The mean absolute error between the achieved service time and its reference value is of 44.9sec with just a PI and of 8sec when adding a adaptive feedforward. Using adaptive feedforward control for cloud computing shows encouraging results. The methodology should be further tested for assessing the advantages of using more complex algorithms both for feedback and feedforward control.

Fig. 1 .

 1 Fig. 1. Cloud control system.

Fig. 2 .

 2 Fig. 2. MapReduce model.

Fig. 3 .

 3 Fig. 3. MapReduce model and control scheme

(

  H1) The disturbance d(t) (number of clients)is bounded (H2) The delays r d and r g are known (identified) (H3) r g ≤ r d (H4) The signs of b d and b g are known ( b d > 0; b g < 0)

Fig. 4 .

 4 Fig. 4. Computer system configuration Fig. 5 shows the evolution of the service time and of the number of nodes in open loop operation of the system. Fig.6shows the evolution of the service time and of the number of nodes when a digital PI controller is used. Fig.7shows the evolution of the service time and of the number of nodes when in addition to the digital PI controller an adaptive feedforward control law is used. One can see that adding adaptive feedforward compensation improves the performance obtained with just a PI controller. The mean absolute error between the achieved service time and its reference value is of 44.9sec with just a PI and of 8sec when adding a adaptive feedforward.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Performance in open loop operation

Fig. 7 .

 7 Fig. 7. Performance in the presence of a PI controller and an adaptive feedforward control

While there is no a consensus in defining "Big Data" one can say that amount of data that can not be processed reasonable fast in a single computing center constitutes "Big Data".

One can view the cloud as a hierarchical system where the upper level is the cloud, the intermediate level is the data center and the basic level is the cluster.

SLA is a contract between a cloud service provider and a service client, where the quality of service is formally defined.

adaptation for getting estimated parameters at t+1 starts once d(t + 1) is measured.

If β g ≥ 1 which correponds to a fractional delay larger than half of the sampling period, it is reasonable to approximate the model with an additional delay of one sampling.

ACKNOWLEDGMENTS

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).


Adrien Varet

Nicolas Prcovic

Cyril Terrioux
email: cyril.terrioux@univ-amu.fr

Denis Hagebaum-Reignier

Yannick Carissan
email: yannick.carissan@univ-amu.fr

BenzAI: A Program to Design Benzenoids With Defined Properties Using Constraint Programming

The BenzAI program can automatically generate benzenoids with defined structural constraints, like the number of hexagons, the number of carbon and/or hydrogen atoms, the existence of symmetries, the number of Kekulé structures, their diameter, or more elaborate criteria like their irregularity. BenzAI allows both automatic generation and manual building of benzenoids including or excluding defined patterns. For all benzenoids with less than 10 rings, it can extract IR spectra from a database that we generated. It computes local aromaticity for closed-shell and mono radical species using circuit-based algorithms. We present practical examples of what can be done with BenzAI for a diversity of cases.

Introduction

This article is devoted to the presentation of the BenzAI software that can generate any benzenoid, perform an analysis of its electronic properties in terms of Kekulé structures, and provide infrared spectra stored in an external database. A great feature of BenzAI is the fact that the generation of the structure can be done with various constraints especially with the constraint to contain or not one or multiple patterns. We believe BenzAI provides data of in-terest for a broad spectrum of fields in physics and chemistry, e.g., astrochemistry, material design, synthetic chemistry, among others.

Benzenoids are polycyclic aromatic hydrocarbons (PAHs) made of fused six-member rings only. They are well-studied in various domains because of their electronic and optical properties. [START_REF] Rakhi | Optoelectronic Properties of Polycyclic Benzenoid Hydrocarbons of Various Sizes and Shapes for Donor-π-Acceptor Systems: A DFT Study[END_REF] They are also studied in detail in interstellar chemistry because they are suspected to be present in interstellar clouds and are believed to act as catalysts for chemical reactions taking place in space. [START_REF] Bouwman | Effect of Molecular Structure on the Infrared Signatures of Astronomically Relevant PAHs[END_REF] The aromatic character of benzenoids can be studied using the Clar sextet rule for ground states. [START_REF] Wang | Quantitative Resonance Theory Based on the Clar Sextet Model[END_REF][START_REF] Anjalikrishna | Antiaromaticity-Aromaticity Interplay in Fused Benzenoid Systems Using Molecular Electrostatic Potential Topology[END_REF][START_REF] Jousselin-Oba | Aromaticity, and Application of peri-Pentacenopentacene: Localized Representation of Benzenoid Aromatic Compounds[END_REF] For large benzenoids, the shape of edges has a large influence on their properties. [START_REF] Lin | Theoretical study of nitrogendoped graphene nanoflakes: Stability and spectroscopy depending on dopant types and flake sizes[END_REF][START_REF] Pelloni | Polygonal current models for polycyclic aromatic hydrocarbons and graphene sheets of various shapes[END_REF] They are also a challenge for synthetic chemists, which lead to the synthesis of remarkable molecules on surfaces. [START_REF] Arikawa | Synthesis and Isolation of a Kinetically Stabilized Crystalline Triangulene[END_REF] All these points can be addressed with Ben-zAI. Aromaticity can be easily evaluated even for radical species, for which BenzAI provides Clar covering and ring bond order analysis. Infrared spectra of multiple species can be determined through the automated consultation of our database, and building of benzenoids with a defined irregularity parameter. As for properties defined by the shape of the edges or the presence of a defined pattern, we have integrated such features, for which practical examples are given in the manuscript.

In the literature, there some software and databases available to generate or study PAH structures. Among them, the CaGe software is based on graph theory and allows the construction of diverse 2D or 3D PAHs. [START_REF] Brinkmann | CaGe -a Virtual Environment for Studying Some Special Classes of Plane Graphs -an Update[END_REF] It allows to generate 3D structures with fused 5member rings and a variety of specific PAHs. Tománek and Frederick have setup a database, from which one can recover the geometries of all isomers of fullerenes with less than 720 carbon atoms. 10 This work intends to study the curvature of these objects. [START_REF] Guan | Local Curvature and Stability of Two-Dimensional Systems[END_REF] Another software to analyse geometry deformation in PAHs is available as a website and computes POAV (π Orbital Axis Vector) values and further developments to study geometrical deformation in PAHs. [12][START_REF] Haddon | Hybridization and the Orientation and Alignment of .Pi.-Orbitals in Nonplanar Conjugated Organic Molecules: .Pi.-Orbital Axis Vector Analysis (POAV2)[END_REF][START_REF] Sabalot-Cuzzubbo | Relating the Molecular Topology and Local Geometry: Haddon's Pyramidalization Angle and the Gaussian Curvature[END_REF] The NASA database is also available online and contains a large number of PAHs. [START_REF] Bauschlicher | The Infrared Spectra of Very Large Irregular Polycyclic Aromatic Hydrocarbons (PAHs): Observational Probes of Astronomical PAH Geometry, Size and Charge[END_REF] This work is intended for infrared spectra studies. Recently, a large set of PAH and their local aromaticity was computed. [START_REF] Alvarez-Ramirez | Database of Nuclear Independent Chemical Shifts (NICS) versus NICSZZ of Polycyclic Aromatic Hydrocarbons (PAHs)[END_REF] Some software are also available to study aromaticity in general such as the Aroma package, which scans NICS values on a grid above a molecular plane. [START_REF] Stanger | Nucleus-Independent Chemical Shifts (NICS): Distance Dependence and Revised Criteria for Aromaticity and Antiaromaticity[END_REF][START_REF] Stanger | Obtaining Relative Induced Ring Currents Quantitatively from NICS[END_REF][START_REF] Gershoni-Poranne | The NICS-XY-Scan: Identification of Local and Global Ring Currents in Multi-Ring Systems[END_REF] The newly developed ims3d.py suite does the same kind of work in a general fashion, which allows for the study of aromaticity in distorted molecules, for which the determination of a plane is ambiguous. [START_REF] Artigas | Visualizing Electron Delocalization in Contorted Polycyclic Aromatic Hydrocarbons[END_REF] In this context, Ben-zAI aims at being an easy to use software flexible and expandable. It handles fragment combination for the generation of benzenoids, which is not covered by the software previously cited and deals with the electronic structure analysis based on graph theory. Furthermore, it is connected to a rich database dedicated to aromaticity. The full feature of this database will be published soon. IUPAC defines aromaticity as a gain in energy of systems that provide for their enhanced thermodynamic stability (relative to acyclic structural analog). [START_REF] Mcnaught | IUPAC. Compendium of Chemical Terminology, 2nd Ed. (the "Gold Book[END_REF] Thus, strictly speaking, it is a global property. However, as chemists, we know that rings in molecules are not equivalent and some are more likely to react (less stabilized by aromaticity) than others. Thus, BenzAI can compute local indices of aromaticity. This is not new and the literature gives a plethora of such indices. The interested reader can look at the article by Solà, who deciphers the reasons for such a large number of criteria. 22 In the first part, we present the methods used to solve the benzenoid generation problem. The second part deals with the computation of local aromaticity and the third shows practical examples of what can be achieved with Ben-zAI. BenzAI is also connected to a database of benzenoids with less than ten cycles. For these molecules, the infrared spectra were calculated at the 6-31G/B3LYP level of theory. Whenever a molecule in this database is generated and if the database is reachable, the results are automatically retrieved and can be explored within BenzAI.

Underlying Concepts And Methods

Generation With Constraints

We briefly recall here the strategy used to generate all the benzenoids having a given number n of hexagons. For further details, please refer to our previously published work on the generation of benzenoids. [START_REF] Carissan | Using Constraint Programming to Generate Benzenoid Structures in Theoretical Chemistry[END_REF][START_REF] Carissan | Exhaustive Generation of Benzenoid Structures Sharing Common Patterns[END_REF] Constraint programming is a branch of artificial intelligence, which focuses on automated reasoning. A generic solver is used whatever the nature of the problem. The problem is expressed in terms of constraints among variables in a given formalism sent to the solver, which outputs the solutions. We use constraint programming, as it easily allows for the extension of the generation of such molecules under combined specific constraints as opposed to bespoke approaches, which could be faster but do not guarantee that additional constraints could be used without a complete rewriting of the code. [START_REF] Dias | Benzenoid Series Having a Constant Number of Isomers[END_REF][START_REF] Dias | The Most Stable Class of Benzenoid Hydrocarbons and Their Topological Characteristics -Total Resonant Sextet Benzenoids Revisited[END_REF] We define the graph of a benzenoid B as the undirected graph B h in which each hexagon is a vertex and an edge links two hexagons that are adjacent in the molecule. This is often referred to as the inner dual graph. Any benzenoid of n hexagons can be embedded in a coronenoid of size at most k(n) = n 2 + 1 . The model relies on the fact that the generation of the n hexagon benzenoids can be translated in graph theory as finding all the connected sub-graphs of the hexagon graph of a coronenoid of size k(n). Let B c(k(n)) h be the graph of the coronenoid of size k(n) and respectively n c and m c its number of vertices and edges. We use a Boolean variable x i per vertex and a Boolean variable y i,j per edge of

B c(k(n)) h . x i is true when the i th hexagon of B c(k(n)) h
is part of B h , while y i,j is true when the i th and j th hexagons of

B c(k(n)) h are fused in B h . The search of all valid B h contained inside B c(k(n)) h
is done under the following constraints:

• If two hexagons are fused in the benzenoid graph B h , they must be part of B h and conversely: y i,j =true ⇔ x i = x j =true;

• B h must have n hexagons:

x i = n;
• B h must be connected: any pair of hexagons can be linked by a path.

• If six hexagons form a ring, the ring in the middle must be part of the graph (coronene case).

Further constraints were added to avoid redundant solutions, or can be added to constrain the solutions to be of a certain symmetry. Many libraries are available to represent and solve efficiently such problems. In this work, we exploit the open-source Java library Choco Solver 27 that natively proposes graph variables and the more usual graph-related constraints (notably connected constraint). An example of such a generation is given as an example in Figure 1. As for the generation of benzenoids with a constraint on fragments, it is based on the concept of extended benzenoid graph. The extension consists in assigning the hexagons of the graph a value: plus (+ colored in green), minus (-colored in red), or neutral (0 colored in orange). When a hexagon is assigned a plus value, it must be in the graph of the solution. On the contrary, a minus value imposes the absence of the hexagon in the solution. Neutral hexagons are neither required nor excluded. This allows for the generation of benzenoid graphs, which represent fragments. The generation of the benzenoid is thus done with the constraint that the extended graph of one or more fragments must be in the solution or excluded from the solution. All details and models used for this are published elsewhere. [START_REF] Carissan | Exhaustive Generation of Benzenoid Structures Sharing Common Patterns[END_REF] As an example, we show in Figure 2 the generation of all benzenoids with seven hexagons made of the combination of two fragments. For fragments, green hexagons indicate that they must be present in the solution, whereas red hexagons must not be in the solution. More detailed information is available in the supplementary information (Figure SI 6).

Fragments to combine

Analysis of the Electronic Structure of a Benzenoid

The values obtained in the framework of this work are subject to limitations intrinsic to the method. The geometries used for the topological aromaticity calculations are not optimized (all hexagons are considered regular) and the aromaticity indices as well as the bond orders are quantities which differ from the ones extracted from quantum chemical calculations. However, the trends are usually good as shown in our recent work. [START_REF] Carissan | Computing the Local Aromaticity of Benzenoids Thanks to Constraint Programming[END_REF]

Local Aromaticity

In our program, local aromaticity is computed following the work of Randić 29 based on minimal conjugated circuits.

Definition 1 Let B be a benzenoid and K one of its Kekulé structure. C is a minimal conjugated circuit of the hexagon h of B (also called h-MCC) with respect to K if C is one of the covering circuits of h having the smallest size.

Such a circuit is a loop made of alternating single and double bonds in a Kekulé structure. In benzenoids, these circuits can only be of size 4n + 2 (with n a positive integer). The algorithm suggested by Randić is described in Algorithm 1. The difficulty lies in obtaining all possible Kekulé structures and their set of conjugated circuits. The detailed implementation of this model using the approaches by the Lin and Fan method [START_REF] Lin | Algorithms for the Count of Linearly Independent and Minimal Conjugated Circuits in Benzenoid Hydrocarbons[END_REF] or the improved method by Lin [START_REF] Lin | Efficient Method for Calculating the Resonance Energy Expression of Benzenoid Hydrocarbons Based on the Enumeration of Conjugated Circuits[END_REF] with constraint programming is thoroughly described in a previous work. [START_REF] Carissan | Computing the Local Aromaticity of Benzenoids Thanks to Constraint Programming[END_REF] In short, this approach consists in enumerating all possible h-MCCs in all Kekulé structures of a benzenoid and attributing to each of them a "resonance energy", R n , which value is tabulated according to their size (4n + 2). The sum of all these energy contributions normalized to the number of Kekulé structures gives the total resonance energy -or global aromaticity -of the benzenoid. The contribution of each ring of the benzenoid, called ring resonance energy (RRE), then corresponds to local aromaticity.

Algorithm 1: Algorithm used to determine local aromaticity of all benzenic cycles of a benzenoid Data: K the complete set Kekulé structures the benzenoid under study {h i } n h i=1 the benzenic cycles (h stands for hexagon) of the benzenoid under study {R i } the resonance energies associated to 6, 10, . . . 4i + 2 conjugated electrons Result: Local aromaticity of each cycle

1 foreach K ∈ K do 2 determine C, the complete set of conjugated circuit of K i 3 foreach C ∈ C do 4 for k ← 1 to n h do 5 if h k is involved in C 6 then 7 e ← number of electrons conjugated in C 8 E[k] ← E[k] + R e 9 end 10 end 11 end 12 for k ← 1 to n h do 13 E[k] ← E[k] |K| 14 end 15 end

Clar cover

From Solà 2013: [START_REF] Solà | Forty Years of Clar's Aromatic π-Sextet Rule[END_REF] Clar's rule states that the Kekulé resonance structure with the largest number of disjoint aromatic π-sextets, i.e., benzene-like moieties, is the most important for the characterization of properties of polycyclic aromatic hydrocarbons (PAHs). Aromatic πsextets are defined as six π-electrons localized in a single benzene-like ring separated from adjacent rings by formal CC single bonds. [START_REF] Hauptmann | The Aromatic Sextet[END_REF] For a given benzenoid, a Clar cover is obtained by maximizing the number of Clar sextets of its electronic structure. In our approach, we can generate radical Clar covers by imposing the constraint that the molecule must carry one, or more radicals.

Ring bond order

Clar formulas can be numerically represented by ring bond orders (RBO). [START_REF] Randić | Local Aromaticity and Aromatic Sextet Theory beyond Clar[END_REF] The bond order is a well-known concept among the chemist community and was first introduced by L. Pauling. [START_REF] Pauling | The Dependence of Interatomic Distance on Single Bond-Double Bond Resonance1[END_REF] In the framework of π conjugated systems, a "pure" π bond as in the ethylene molecule, has a bond order of 1. When one sums up all bond orders of the six C-C bonds forming each ring, one obtains an RBO. The archetypal RBO of a π system is the RBO of the benzene molecule and is equal to 3. Thus, the closest to 3 an RBO is, the more aromatic the ring is. In BenzAI, the bond order b of a benzenoid B is calculated as the ratio, in percent, of the number of Kekulé structures in which b is a double bond and the total number of Kekulé structures of B.

Interface

Molecules are grouped into collections. These collections can be filled by the user via drawing or importing a molecule or automatic generation using constraints. Within a collection, one can select a set of molecules to create a new collection or for particular calculations, for instance, calculation of local aromaticity. Molecules can be appended to an existing collection by unit or by batches of selected molecules. A filter can be applied to a collection to generate a new collection. Figure 3 shows screenshots of BenzAI for a typical input and the corresponding output.

Drawing and Import/Export

The drawing module allows for the easy drawing of a benzenoid on a hexagonal grid. The default grid is made of three concentric crowns of hexagons: the smallest crown is just a single ring, i.e. the benzene molecule, two crowns correspond to the coronene molecule, and so on. Such concentric crowns of hexagons are thus called coronenoids, as they contain a coronene substructure. A simple click on a hexagon of the grid adds it to the desired molecule. There can be holes in the benzenoid but all cycles must be fused. Once drawn, the molecule is added to the desired collection of molecules. Hence many molecules can be drawn by hand and analyzed.

Benzenoids can also be exported from the application using the export menu, either individually or by collection. Benzenoid structures can be exported from the application using the export menu, either individually or by collection, to the following formats: graph format, 36 CML (Chemical Markup Language), XYZ and PNG format (Portable Network Graphics). They can be imported later on in BenzAI only using the graph format. Importing a collection consists simply in selecting several files (these files must be located in the same directory).

Generation of Benzenoid Structures at Will

The BenzAI program can generate benzenoid structures with the following constraints, which can be combined: [START_REF] Carissan | Using Constraint Programming to Generate Benzenoid Structures in Theoretical Chemistry[END_REF] Based on the chemical composition:

• Total number of hexagons

• Total number of carbon atoms

• Total number of hydrogen atoms Based on the structure:

• Irregularity parameter ξ and/or number of solo, duo, trio, and quarto hydrogen atoms (see definition below)

• Diameter: it is defined as the maximal length of the shortest path connecting any hexagon to another. This length is equal to the number of edges that are crossed along the path. For instance, the molecule shown in Figure 4 has a diameter of 4, the maximal length of the shortest path is either 0-1-2-5-7 or 3-1-2-5-7.

• Coronenoid, rectangle, rhombus, catacondensed

• With/without holes (Coronoids) • With/without a combination of patterns.

Patterns can be joint, disjoint, or share hexagons, they can be saved and loaded in a user-friendly manner.

• Using symmetries: 13 types of symmetry are possible for benzenoids, using either reflection or rotation symmetry (or any combination of both of types) referred to with their Schoenflies symmetry point group symbol. [START_REF] Kovic | On Symmetries of Benzenoid Systems[END_REF] Based on generation itself:

• Number of solutions desired

• Total time for the generation Please note that the number of solutions can be extremely high. So we advise the user to make extensive use of the total time generation criterion as unwise requests can exhaust his patience or RAM, whichever comes first.

Filtering

Any collection can be filtered with the same criteria used for the generation, except for the number of solutions or total time. In addition to the generation criteria, the number of Kekulé structures can be used as a new filtering criterion. When set to zero, this allows a simple way to filter non-kekulean benzenoids from the solutions of a collection. The filtering option is thus a convenient way to generate numerous sub-collections according to some criteria from a common parent collection having given composition or structural characteristics. Hence, one can explore a collection through many aspects.

Sorting

Any generated collection can be sorted out in increasing or decreasing order with respect to structural parameters (number of carbon atoms, hydrogen atoms, rings, or irregularity) or with respect to the number of Kekulé structures. As an example, let us generate the set of all benzenoids with 32 carbon atoms. This leads to a collection of 330 molecules. Sorting by decreasing order of ξ shows that only on such molecule has a ξ value of 1. This molecule is shown in Figure 4.

Computations on a collection

Using the "Computations" tab of BenzAI, one can perform two types of analysis on either an entire collection or on manually selected solutions of a collection:

• an electronic structure analysis according to the aromaticity criteria previously mentioned (RRE, Clar cover, RBO).

• a histogram analysis with respect to the irregularity parameter These analyses will be illustrated below in the "Use Cases" section.

Database

If possible, BenzAI can retrieve data from the database we computed, which contains all possible benzenoids of less than 10 cycles and their vibrational spectra computed at the 6-31g/B3LYP level of theory. All geometries were optimized and do not present any negative frequency. These spectra can be seen individually or they are merged using Boltzmann statistics based on the DFT energies of each molecule. For example, let us generate the set of benzenoids with four hexagons and retrieve their vibrational spectra from the database, Figure 5. BenzAI allows for the visualization of 3 spectra, each of them being a Boltzmann's weighted spectrum of all molecules with the same number of carbon and hydrogen atoms.

Use Cases Aromaticity of a Set of Molecules

In order to show the capabilities of BenzAI for the computation of local aromaticity and re- The local aromaticity calculation and the ring bond order representation fully agree. Pentacene has five Clar structures, which differ by the cycle carrying the Clar sextet. Thus, the ring bond order representation gives an equal value for each ring. However, in acenes, the local aromaticity is larger for the central cycles. This is well captured by the local aromaticity calculation based on electronic circuits. The Benzo[a]tetraphene is similar to pentacene with two sliding sextets (Figure 7). The local aromaticity of the molecule is well described and fully in line with the Clar cover. The ring bond order description lacks accuracy for the rightmost ring (by construction) but the description of aromaticity with electronic circuits captures the fact that this ring is less aromatic than its neighbour. maticity, in agreement with the Clar cover, in which the central ring cannot carry a sextet. It is well-known that perylene can be seen as two naphthalene moieties linked by two bonds of essentially single character. However, the ring bond order of the central ring is not null, because it takes into account the Kekulé structures for which its bonds can have a double character. This is never the case for the two central bonds but the other ones are double in the Kekulé structures of naphthalene. Finally, the C 17 H 13 molecule is not kekuleable, because it has an odd number of π electrons. The local aromaticity cannot be evaluated by methods based on electronic circuits made of alternating single and double bonds because of the radical delocalization. The Clar cover picture shows a carbon, which carries the radical and all its second neighbors can have a radical character as well, Figure 8. Then the ring bond order picture allows for an evaluation of local aromaticity as well as the localization of the radical. In the case of a radical species, the ring bond order is obtained by summing on all possible electronic structures, which have one radical and maximize the number of Clar sextets.

Irregularity

The irregularity parameter ξ, is defined as follows:

2 ξ = N trio + N quarto N solo + N duo + N trio + N quarto (1)
It is used in astrochemistry to characterize the regularity of polycyclic aromatic hydrocarbons. BenzAI allows for the calculation of the distribution of ξ in a collection of molecule and the generation of molecules with a constraint ξ value or constraint values of N solo , N duo , N trio or N quarto . Recently, Bouwman suggested that the NASA PAH database lacks large molecules (with 34 carbon atoms or more) with high irregularity values. [START_REF] Bouwman | Effect of Molecular Structure on the Infrared Signatures of Astronomically Relevant PAHs[END_REF] In Figure 9, after generating the full set of benzenoids with nine hexagons or less and filtering them by the number of carbon atoms (less than or equal to 34 and more than 34 carbon atoms), we compare the distribution of ξ in both sets. As suggested, the NASA database over-represents molecules with a low value of ξ. This is of the highest importance because the kind of hydrogen of a PAH (solo, duo, etc.) was shown to have an influence on its IR absorption band. [START_REF] Bauschlicher | The Infrared Spectra of Very Large Irregular Polycyclic Aromatic Hydrocarbons (PAHs): Observational Probes of Astronomical PAH Geometry, Size and Charge[END_REF][START_REF] Hony | The CH Out-of-Plane Bending Modes of PAH Molecules in Astrophysical Environments[END_REF][START_REF] Ricca | Polycyclic Aromatic Hydrocarbons with Armchair Edges: Potential Emitters in Class B Sources[END_REF] Molecular Design

Generation of Large Symmetric Molecules

BenzAI allows for the generation of a determined number of benzenoids with a high number of benzene rings and specific properties. In the following, we shall generate five benzenoids with 24 cycles, a mirror symmetry and constrained ξ values (ξ < 0.2, ξ = 0.2, 0.4, 0.6, 0.8 and ξ > 0.8). Results of this generation are shown in Figure 10. It is interesting to see that higher ξ values lead to more branched structures (as opposed to more compact for low ξ values), which have non-aromatic cycles. Of course, such the generation of only five benzenoids per ξ value cannot be generalized, but if one combines these results with the distribution of ξ shown in Figure 9, it is of interest to extend the available databases by a tool such as BenzAI, to better represent the variety of benzenoids.

Prediction of Magnetic Properties

Some of us have recently studied chiral nonplanar PAHs based on triphenylene. These molecules were obtained by a trimerization of helicenes (formally a [2+2+2] reaction, Figure 11). [START_REF] Artigas | Visualizing Electron Delocalization in Contorted Polycyclic Aromatic Hydrocarbons[END_REF] The trimerization leads to the formation of a new C 6 cycle, which is referred to as the central cycle in the following discussion. Such molecules could exhibit magnetic properties if their π extensions would stabilize unpaired electrons. Thus, we search monomeric units with stable unpaired electrons and an even number of carbons, i.e. concealed non-kekulean benzenoids. [START_REF] Cyvin | Search for Concealed Non-Kekulean Benzenoids and Coronoids[END_REF] We know that the C 13 H 9 triangulene carries a radical localized on the second neighbor of the central carbon atom, Figure 12. [START_REF] Rico-Sotomayor | Triangulenes and theirs ions: reaching the limits of Clar's rule[END_REF][START_REF] Misra | Theory for Radical Benzenoids[END_REF] There must be two such units by monomer for the total molecule to have an even number of electrons and these two units must not share a ring. The combination of these constraints can be handled by BenzAI. In the generator, we define the C 13 H 9 as a pattern as shown in Figure 13 and ask for two such disjoint units. In the framework of this article, we decided to use only the two solutions shown in Figure 14: one solution A which should exhibit a diradical character (hence a trimer of septet ground state) and B, which should lead to a closed-shell wave-function. By construction, 3-A (Figure 15) should carry six unpaired electrons localized on the triangulene moieties whereas 3-B should be a closed-shell molecule.

In order to verify these hypotheses, we computed the energy of the spin states of both molecules, Table 1.

Hence, the use of BenzAI allowed us to find a

Conclusion

This article presents the BenzAI program, an open-source software available at the following URL: https://benzai-team.github.io/ BenzAI/. This program allows for the generation of benzenoids under structural constraints. Then, it allows for the determination of the local aromaticity of the hexagons and the calculation of IR spectra provided that the generated benzenoids have less than 10 hexagons. In this article, we show concrete applications of this software in diverse aspects of chemistry. The beauty of constraint programming lies in the ease to add new constraints in the future without having to modify a bespoke program. Thus, we are already working on new features like allowing for cycles with more or less than six carbon atoms, which could extend the possibilities of the program.

Supporting Information Available

Complementary to Figure 2, the generation of benzenoids with seven hexagons, made of all possible combinations of two fragments with the fragments colored for easier recognition.

Complementary to Figure 6, The full set of all benzenoids with five rings is available together with their local aromaticity computed with the Lin, Lin and Fan and ring bond order methods as well as their Clar cover.

Data and Software Availability

The BenzAI is available free of charge at https: //benzai-team.github.io/BenzAI/. The stared positions and equivalent carry a radical. The purple cycle will be the central cycle of the trimer, the blue and green triangulene units each carry a radical. These two electrons cannot be paired: it is a concealed non-kekulean benzenoid. Bonds surrounded in red are single in all triplet Kekulé structures.

Figure 1 :

 1 Figure 1: Generation of the only four benzenoids with six cycles and a mirror symmetry.

Figure 2 :

 2 Figure 2: Example of the generation of benzenoids with seven hexagons, made of all possible combinations of two fragments. For fragments, green hexagons indicate that they must be present in the solution, whereas red hexagons must not be in the solution. More detailed information is available in the supplementary information (Figure SI 6).

Figure 3 :

 3 Figure 3: Screenshots of the BenzAI graphical interface: generator interface (left) result of a generation (right).

Figure 4 :

 4 Figure 4: Only benzenoid with 32 carbon atoms with a ξ value of 1 obtained by filtering the collection of 330 benzenoids with 32 carbon atoms by decreasing ξ value.

Figure 5 :

 5 Figure 5: Vibrational spectra extracted from the database for the benzenoids with four hexagons. Three spectra are obtained as these benzenoids can have either 16, 17 or 18 carbon atoms. The spectrum obtained for 18 carbon atoms is a Boltzmann weighted spectrum of the spectra of the five benzenoids.

Figure 6 :

 6 Figure 6: From left to right: local aromaticity using Lin's method, Clar cover and ring bond order for a subset of 5-member ring benzenoids. Darker blue color indicates higher local aromaticity or higher ring bond order. For the Clar cover representation, only one among all possible Clar structures is given. Note that the circle around one atom indicates a radical.

Figure 7 :

 7 Figure 7: Clar structures with two sextets of benzo[a]tetraphene.

Figure 8 :

 8 Figure 8: Clar structures with three sextets of C 17 H 13 . The atoms surrounded by a circle carry the radical.

Figure 9 :

 9 Figure 9: Comparison of the distributions of the irregularity parameter (less or equal than 34 carbon atoms (a) and more than 34 carbon atoms (b)) as shown by Bouwman et al. 2 (blue) and this work (orange) for the full set of benzenoids with nine or less hexagons.

Figure 10 :

 10 Figure 10: Five first solutions obtained by six successive benzenoid generation done with BenzAI with the following constraints: 24 benzene rings, mirror symmetry, and an irregularity parameter ξ < 0.2, ξ = {0.2, 0.4, 0.6, 0.8} and ξ > 0.8. Blue color scale indicates Lin aromaticity. White color indicates no local aromaticity. All these solutions were obtained within a few seconds on an Intel core i7 laptop.

Figure 11 :

 11 Figure 11: Formal [2+2+2] trimerization of three monomeric units. The arrows indicate the possible π extensions of these monomers.

Figure 12 :

 12 Figure 12: Singly occupied molecular orbital at the Hückel level of theory (left) and spin density of the C 13 H 9 radical computer at the PBE0/def2-SV(P) level of theory (right).

Figure 13 :Figure 14 :

 1314 Figure 13: Pattern definition for the C 13 H 9 radical carrier fragment (left) and of the central cycle fragment (right). Green cycles must be in the solution whereas red cycles are excluded from the solution.

Figure 15 :

 15 Figure 15: Trimerization of compound A (left) and of compound B (right). The dashed lines indicate the bonds which will lead to the formation of the central cycle. The bold lines indicate that the part of the molecule comes towards the reader, whereas the hollow lines go away from the reader. Blue and green triangulene moieties carry one radical.

Table 1 :

 1 Energy of the spin states of 3-A and 3-B in kJ/mol obtained at the PBE0/def2-SV(P) level of theory.

		3-A 3-B
	singlet 241	0
	triplet	83	13
	quintet 80	34
	septet	0	53
	concealed non-kekulean PAH with a high spin
	ground state.

Acknowledgement The authors thank the ANR project DEMOGRAPH (ANR-16-CE40-0028) and the ministère de l'éducation nationale for funding the PhD of Adrien Varet.

Graphical TOC Entry

BenzAI, benzenoids made by A.I.


Feda Almuhisen

Nicolas Durand

Leonardo Brenner

Mohamed Quafafou

Sequential Patterns for Spatio-Temporal Traffic Prediction

Keywords: Closed sequential patterns, emerging patterns, Markov chains, moving object trajectory data, traffic prediction and monitoring

This paper proposes a method for predicting the traffic status of a city within time windows. The method takes advantage of spacepartitioning, closed sequential pattern extraction, emerging pattern detection, and Markov chain modeling. From trajectories, we identify active regions in which moving objects mostly visit. The traffic status of each region is detected based on continuous tracking of closed sequential patterns evolution over time. Based on the proposed Markov model, the near-future status of traffic is then predicted. The traffic status is reported on maps and can be used to enhance future city transportation. The experiments on realworld data sets show that the proposed method provides promising results.

INTRODUCTION

Smart city concept has become essential to manage the city and to provide innovative services to citizens. Traffic prediction is one of these services and aims to improve the mobility in the city [START_REF] Nagy | Survey on traffic prediction in smart cities[END_REF]. Indeed, the increasing of the urban population makes the creation of intelligent transportation systems strategic. IT platforms and sensor networks play a crucial role in the management of the spatiotemporal data, preprocessing, traffic prediction and services [6] [30]. Traffic prediction [START_REF] Yuan | A Survey of Traffic Prediction: from Spatio-Temporal Data to Intelligent Transportation[END_REF] refers to: (1) traffic status prediction: congestion prediction in future time, (2) traffic flow prediction: prediction of traffic flows in future time, and (3) travel demand prediction: prediction of people's travel demands in order to better dispatch taxi vehicles for different regions. The practical goal of this work is to deliver a global view of current or future traffic status in regions of interest (see Figure 6). Prior short term traffic predictions models can be used directly by experts to take relevant actions against prevailing or incipient traffic trends. Developing such prediction systems needs trajectory historical data analysis which can be also used in other applications such as improving route navigation and urban area planning [START_REF] Lin | Mining GPS data for mobility patterns: A survey[END_REF] [START_REF] Pan | Trace analysis and mining for smart cities: Issues, methods, and applications[END_REF]. Markov models have been widely used in the field of traffic prediction in order to improve traffic flow [START_REF] Zhao | Traffic flow prediction based on optimized hidden Markov model[END_REF], predict traffic conditions [START_REF] Qi | A Hidden Markov Model for short term prediction of traffic conditions on freeways[END_REF] [START_REF] Zaki | Traffic congestion prediction based on Hidden Markov Models and contrast measure[END_REF], forecast travel speed [START_REF] Zhang | A research on driving condition prediction for HEVs based on Markov chain[END_REF], and predict routes [START_REF] Lassoued | A Hidden Markov Model for Route and Destination Prediction[END_REF].

In this paper, we focus on traffic status prediction by analysing the evolution of traffic. In fact, we analyze the trajectories of vehicles in historical data to detect and visualize the evolution of traffic in the city. We propose a new method for predicting the traffic status (increasing, decreasing, etc.) in the next future. This method is called SP4TP (Sequential Patterns for Traffic Prediction) and its originality is to leverage space partitioning, closed sequential pattern mining [START_REF] Fournier-Viger | A survey of sequential pattern mining[END_REF], emerging pattern mining [START_REF] Dong | Efficient mining of emerging patterns: Discovering trends and differences[END_REF], and Markov chain modeling [START_REF] Stewart | Introduction to the numerical solution of Markov chains[END_REF]. The space area is divided into regions (uniformed squares) using a given spatial granularity value. Considering a selected region and a time windows, our method based on Markov models predicts the next traffic status. To train the proposed Markov model, we determine the traffic status by detecting the evolution types (Emerging, Decreasing, etc.) of sequential formal concepts (containing closed sequential patterns) extracted from trajectory data [START_REF] Almuhisen | Detecting behavior types of moving object trajectories[END_REF] [START_REF] Almuhisen | Sequential Formal Concepts over Time for Trajectory Analysis[END_REF]. The prediction results are visualized in geotagged maps. To the best of our knowledge, no prior methods based on sequential patterns, emerging patterns, and Markov models have been proposed to build traffic prediction models. Two real-world data sets from two cities, Beijing and San Francisco, are used in the experiments to evaluate our method and compare it to an itemset based approach and a baseline method which does not use patterns. We also assess the advantage of using sequential patterns in our context.

The rest of the paper is organised as follows. Section 2 presents the notions useful for understanding the paper. Related work is discussed in Section 3. The proposed method, SP4TP, is detailed in Section 4. Section 5 presents the experimental results. Finally, we conclude in Section 6.

PRELIMINARIES

In this section, we present the notions useful to understand the proposed method (in Section 4): frequent itemsets, sequential patterns, formal concepts, emerging patterns, and Markov chains. We illustrate some notions on an example: a set of moving object trajectories (see Figure 1 (a)).

Processing Trajectories to Sequences

A trajectory 𝑇𝑟 𝑗 𝑖 is an ordered set of points. A point is a timestamped GPS-location, i.e., (time, latitude, longitude) [START_REF] Zheng | Trajectory data mining: An overview[END_REF]. Figure 1 (b) presents 4 trajectories over two days. The space area is divided into uniform grid cells, e.g., 40 meters (corresponding to the spatial granularity value). Each trajectory is split according to the time granularity value, e.g., one day, and labelled by the corresponding time value. Let us note that a trajectory can be also split into subtrajectories reflecting several journeys in relation to the detected long stops, for instance. 𝑇𝑟 𝑗 1 has been slit into 5 sequences: 𝑠 1 , 𝑠 2 , . . . , 𝑠 5 . Trajectories are then mapped in the raster area to obtain the data set presented in Figure 1 (c). For instance, the sequence 𝑠 2 occurs in Day 1 and is composed of grid cells 20, 16 and 15.

Frequent Itemsets and Sequential Formal Concepts

Frequent itemsets have been introduced in [START_REF] Agrawal | Mining association rules between sets of items in large database[END_REF]. Their extraction corresponds to finding the sets of items (i.e., attribute values) that appear simultaneously in at least a certain given number of transactions (i.e., objects) recorded in a database. In the case of ordered items in the data set (like trajectory data), sequential patterns are extracted [START_REF] Fournier-Viger | A survey of sequential pattern mining[END_REF]. Let I = {𝑖 1 , ..., 𝑖 𝑘 } be a set of items. A subset 𝐼 ⊆ I is called an itemset. |𝐼 | denotes the number of items of 𝐼 . A sequence is an ordered list of itemsets. It can be represented as 𝑠 = <(𝐼 1 )(𝐼 2 )...(𝐼 𝑘 )> where each 𝐼 𝑖 is a subset of I and 𝐼 𝑖 comes before 𝐼 𝑗 if 𝑖 ≤ 𝑗. The aim is to find all the frequent subsequences in a sequential database. Let S = <𝑠 1 , ..., 𝑠 𝑛 > be a list of sequences in a database. A sequence 𝑠 1 = <𝑎 1 , . . . , 𝑎 𝑚 > is a subsequence of another sequence 𝑠 2 = <𝑏 1 , . . . , 𝑏 𝑛 >, denoted by (𝑠 1 ⊆ 𝑠 2 ), if there are integers 1 ≤ 𝑖 1 ≤ 𝑖 2 ≤ • • • ≤ 𝑖 𝑚 ≤ 𝑛, such as 𝑎 1 ⊆ 𝑏 𝑖 1 , 𝑎 2 ⊆ 𝑏 𝑖2 , . . . , and 𝑎 𝑚 ⊆ 𝑏 𝑖𝑚 . 𝑠 2 is called a super-sequence of 𝑠 1 (𝑠 2 contains 𝑠 1 ). The support of a sequence 𝑠 𝑛 , noted 𝑠𝑢𝑝𝑝𝑜𝑟𝑡 (𝑠), is the number of occurrences of 𝑠 𝑛 in S. A sequence 𝑠 is said to be a sequential pattern (also called a frequent sequence) if 𝑠𝑢𝑝𝑝𝑜𝑟𝑡 (𝑠) ≥ 𝑚𝑖𝑛𝑠𝑢𝑝, where 𝑚𝑖𝑛𝑠𝑢𝑝 is the minimum support threshold value set by the user. One of the limitations of using sequential pattern algorithms, is that it generates the full set of sequential patterns and scans the database many times when a longer pattern exists. This is why, closed sequential patterns are interesting. Closed sequential patterns are the set of sequential patterns that are not included in other sequential patterns having the same support because they represent the largest frequent subsequences common to sets of sequences. A sequential pattern 𝑠 1 is closed if 𝑠 2 with 𝑠 1 ⊂ 𝑠 2 and 𝑠𝑢𝑝𝑝𝑜𝑟𝑡 (𝑠 1 ) = 𝑠𝑢𝑝𝑝𝑜𝑟𝑡 (𝑠 2 ).

In Figure 1 (c), the sequence 𝑠 = <(19)(23) [START_REF] Stewart | Introduction to the numerical solution of Markov chains[END_REF]> is contained in 𝑠 3 end 𝑠 7 . If 𝑚𝑖𝑛𝑠𝑢𝑝 is set to 2, 𝑠 is frequent. Moreover, 𝑠 is the maximal subsequence shared by 𝑠 3 and 𝑠 7 , 𝑠 is thus a closed sequential pattern.

Given a sequential data set, there is a unique ordered set that describes the inherent lattice structure defining natural groupings and relationships among the objects and their related sequences [START_REF] Garriga | Lattice theory for sequences[END_REF]. This structure is called a sequence concept lattice. Each element of the lattice is a couple, called sequence formal concept, composed of a set of objects 𝑂 (the extent) and a set of sequences 𝑆 (the intent). We call sequential formal concepts, the sequence formal concepts that have at least 𝑚𝑖𝑛𝑠𝑢𝑝 objects in their extent, and sequential concept lattice, the lattice formed using the sequential formal concepts.

In the example (Figure 1 (c), 𝑚𝑖𝑛𝑠𝑢𝑝=2), the sequential formal concepts are: ({𝑠

Emerging Patterns

Given two data sets D 𝑖 , D 𝑗 , the quantitative evaluation of the contrast between data sets brought by a pattern is measured by its growth rate [START_REF] Dong | Efficient mining of emerging patterns: Discovering trends and differences[END_REF]. The growth rate of a pattern 𝑋 from D 𝑗 to D 𝑖 is defined by [START_REF] Agrawal | Mining association rules between sets of items in large database[END_REF].

𝐺𝑅 𝑗,𝑖 (𝑋 ) = |D 𝑗 | |D 𝑖 | × 𝑠𝑢𝑝𝑝𝑜𝑟𝑡 (𝑋, D 𝑖 ) 𝑠𝑢𝑝𝑝𝑜𝑟𝑡 (𝑋, D 𝑗 ) . (1)
The more 𝐺𝑅 𝑗,𝑖 (𝑋 ) is high, the more 𝑋 characterizes D 𝑖 compared to D 𝑗 . Given a threshold value 𝑚𝑖𝑛𝑔𝑟 > 1, a pattern 𝑋 is an emerging pattern from

D 𝑗 to D 𝑖 if 𝐺𝑅 𝑗,𝑖 (
= <(19)(23)(22)> is not present in Day2 (𝑠𝑢𝑝𝑝𝑜𝑟𝑡 (𝑠, 𝐷𝑎𝑦2) = 0) whereas its support in Day1 is equal to 2.
Thus, 𝑠 is a jumping emerging sequential pattern from Day2 to Day1.

Predictive Models and Markov Chains

Markov chains are a mathematical formalism from probability theory which is used to model the possible states of a system and the transitions among these states over time [START_REF] Stewart | Introduction to the numerical solution of Markov chains[END_REF]. Figure 2 The change from the state 𝑒 1 to state 𝑒 2 is activated by the occurrence of the transition with probability 𝜆 𝑒 1 𝑒 2 . Let us note that the transition rates are associated in case of continuous time, while we have probabilities in discrete time. The duration of residence in each state are random variables of an exponential distribution. At an instant of time, the change to a next state depends only on the current state and not on the time elapsed in that state. Markov chain models can also be represented by a transition matrix. For example, for the two states 𝑒 1 , 𝑒 2 , Four combinations are possible and can be represented in a transition matrix as follows:

1 𝑒 2 𝜆 𝑒 1 𝑒 1 𝜆 𝑒 2 𝑒 1 𝜆 𝑒 1 𝑒 2 𝜆 𝑒 2 𝑒 2
𝑒 1 𝑒 2 P = 𝑒 1 𝑒 2 𝜆 𝑒 1 𝑒 1 𝜆 𝑒 1 𝑒 2 𝜆 𝑒 2 𝑒 1 𝜆 𝑒 2 𝑒 2 (2)
The sum of all the probabilities of the transitions from a state must be 1.0. Two types of prediction can be obtained from such models: long term and short term predictions. A stationary analysis for Markov models allows to make a long-term prediction. This is the probability that remains unchanged as time progresses. To compute the stationary probabilities, a vector 𝜋 which satisfies the equation 𝜋.𝑃=0 is needed. This vector 𝜋 can be found by two different ways: iterative methods (Gauss-Seidel method, Power method, etc.) or direct methods [START_REF] Stewart | Introduction to the numerical solution of Markov chains[END_REF]. While the transient analysis allows to predict the evolution of the model from a given state and for a specified period, mainly it represent the distribution of the system at a given time (short-term prediction). In order to find the transient probabilities, the uniformization method [START_REF] Jensen | Markoff chains as an aid in the study of Markoff processes[END_REF] can be employed to compute how many 𝜋 = 𝑃(𝜋-1) are necessary to arrive to 1 time unit from the initial state. The transient distribution at time 𝑡 is obtained by computing an approximation to the infinite summation. Equation 3 presents the uniformization formula where 𝜋 (0) is an initial distribution, and (Γ𝑡 ) 𝑘 𝑘! 𝑒 -Γ𝑡 is the probability that a Poisson random variable with uniform rate parameter Γ takes the value 𝑘 (see Equation 4).

𝜋 (𝑡) = ∞ ∑ 𝑘=0 𝜋 (0)𝑃 𝑘 (Γ𝑡) 𝑘 𝑘! 𝑒 -Γ𝑡 (3)
P{𝑁 (𝑡) = 𝑘 } = (Γ𝑡) 𝑘 𝑘! 𝑒 -Γ𝑡 𝑘 ≥ 0, 𝑡 ≥ 0 (4)
The proposed method in this paper is concerned with Discrete Time Markov Chain (DTMC) where the state of a system is observed at discrete set of times, including periodicity and recurrence in a finite state-space case.

RELATED WORK

Markov chains are often used in a traffic prediction context contrary to sequential patterns. Sequential pattern mining is used in a more general context to analyze spatio-temporal data.

In [START_REF] Wang | Discovering closed frequent patterns in moving trajectory database[END_REF] the authors adapted the classical closed frequent sequence pattern mining algorithm CloSpan to moving trajectory frequent pattern mining. This algorithm, called MTCloSpan, is tested according to the 𝑚𝑖𝑛𝑠𝑢𝑝 value and the characteristics of trajectories (number and length). The authors of [START_REF] Maciag | Discovery of closed spatiotemporal sequential patterns from event data[END_REF] proposed an algorithm, called CST-SPMiner, to discover all participation index strong spatiotemporal sequential patterns from event data. Nevertheless, event data are quite different to moving object trajectory data and this work can not be used in our context. Sequential patterns are very usual. However, our method extracts more than sequential patterns, we compute sequential formal concepts which are more informative (transactions IDs, sequence sets) and allow easier exploitation.

Wang et al. proposed a periodic pattern mining of structural evolution in the case of complex trajectories [START_REF] Wang | Mining Evolution Patterns from Complex Trajectory Structures[END_REF]. Mining periodic patterns from spatio-temporal trajectories reveals information about people's regular and recurrent movements and behaviors. In [START_REF] Jindal | Spatiotemporal Periodical Pattern Mining in Traffic Data[END_REF], the authors proposed a method to mine spatio-temporal periodic patterns in the traffic data and use these periodic behaviors to summarize the huge road network by clustering. In our work, we do not interest in periodicity but in evolution of extracted patterns over time (between two studied time windows, are they emerging, jumping, lost, etc.?). The evolution types of our method characterize the traffic status trends. We want to predict the traffic status in the next future.

Many recent works has been devoted to develop methods based on Markov models in order to predict, for instance, traffic conditions, congestions, routes and destinations. In [START_REF] Qi | A Hidden Markov Model for short term prediction of traffic conditions on freeways[END_REF], the authors proposed a probabilistic approach using Hidden Markov Model (HMM) to model the stochastic variation of traffic conditions and predict the traffic conditions over short time periods. The model defines traffic states based on speed observations of vehicles. Zaki et al. presented in [START_REF] Zaki | Traffic congestion prediction based on Hidden Markov Models and contrast measure[END_REF] a model, based on HMM, to define the traffic states during peak hours in two dimensional space. The model uses mean speed and contrast to capture the variability in traffic patterns.

Most of works based on vehicles speed use sensors at specific road places. We can say that these methods are local while our approach is global. Our method operates at a more macroscopic level by using any mobility GPS data.

Bouyahia et al. [START_REF] Bouyahia | A two-stage road traffic congestion prediction and resource dispatching toward a self-organizing traffic control system[END_REF] presented a two stage method for predicting traffic congestions. The first stage uses a Markov random field to model and predict the propagation of traffic congestion. The road network is assimilated to a graph. The second stage allows to organize the deployment of the traffic regulation resources using a Markov decision process. In [START_REF] Rathore | A Scalable Framework for Trajectory Prediction[END_REF], Rathore et al. proposed a framework for trajectory prediction from GPS data. The framework is based on trajectory clustering to find frequent route patterns. For each cluster of trajectories, a Markov chain model is trained to make future trajectory prediction. As our method, a special step is done before training the Markov models. In our case, this is the computation of the sequential formal concepts and the detection of their evolution types.

The previously mentioned works show that Markov models are widely used in a traffic context. Nevertheless, to the best of our knowledge, there is no work combining sequential patterns, detection of patterns evolution, and Markov models to predict trends and traffic evolution.

SP4TP METHOD

In this section, we present the proposed method, called SP4TP (Sequential Patterns for Traffic Prediction), which considers extracting traffic evolution status in a study area to make predictions. Figure 3 shows a general overview of this method. The main steps are:

• Spatio-temporal preprocessing (see Section 4.1),

• Pattern mining and status detection (see Section 4.1),

• Region traffic model construction (see Section 4.2),

• Prediction (see Section 4.2).

For a given region and two time windows, the method computes the closed formal concepts (i.e., closed sequential patterns associated to sequence IDs) and detects their evolution types between the times windows. These evolution types (considered as traffic status) are used as input to train the proposed Markov model, which will use to predict the traffic status in the next future windows.

Pattern Mining and Status Detection

In order to extract specific region traffic status, we applied the method mentioned in Section 2.1. The studied area is divided into grid cells according to a spatial granularity value. Trajectories are then mapped in the discretized area. This discrete representation of trajectory is used to compute the sequential formal concepts (see Section 2.2) according to 𝑚𝑖𝑛𝑠𝑢𝑝 (the minimum threshold value to consider the traffic). Let (𝑂, 𝑆) be a sequential formal concept, 𝜃 be the minimum threshold value of emergence, and 𝜖 be the error tolerance. In this context, the pattern evolution of each concept (𝑂, 𝑆) between time 𝑖 and time 𝑖-1 is detected by computing 𝐾 𝑖 , an indicative value of 𝑆, between time 𝑖 and time 𝑖-1 (see Equation 5). 𝐾 𝑖 is related to the growth rate of a pattern (see Section 2.3).

𝐾 𝑖 (𝑂, 𝑆) = 𝑐𝑜𝑢𝑛𝑡 (𝑂, 𝑡 𝑖 ) 𝑐𝑜𝑢𝑛𝑡 (𝑂, 𝑡 𝑖-1 ) (5)
Where 𝑐𝑜𝑢𝑛𝑡 (𝑂, 𝑡 𝑗 ) is the number of trajectories of 𝑂 labeled by the corresponding 𝑡 𝑗 . We apply the following rules based on the 𝐾 𝑖 value. If 𝐾 𝑖 (𝑂, 𝑆) = 0 ± 𝜖 , then the type of 𝑆 is lost. If ((𝐾 𝑖 (𝑂, 𝑆) > 𝜃 ) ∧ (𝜃 = 1)) ∨ ((𝐾 𝑖 (𝑂, 𝑆) ≥ 𝜃 ) ∧ (𝜃 > 1))) then the type of 𝑆 is emerging. If 𝐾 𝑖 (𝑂, 𝑆) < 𝜃 , then the type of 𝑆 is decreasing. If 𝐾 𝑖 (𝑂, 𝑆) = 1 ± 𝜖, then the type of 𝑆 is latent. If 𝐾 𝑖 (𝑂, 𝑆) = +∞, then the type of 𝑆 is jumping. In a more intuitive way, Emerging means that the presence of the pattern increased in 𝑡 𝑖+1 compared to 𝑡 𝑖 ; Decreasing means that the presence of the pattern decreased in 𝑡 𝑖+1 compared to 𝑡 𝑖 ; Latent means that the presence of the pattern is quite similar in both time. Jumping means the pattern which was absent in 𝑡 𝑖 , appeared in 𝑡 𝑖+1 . Lost means that the pattern disappeared in 𝑡 𝑖+1 .

Table 1 presents some examples of evolution type detection from the sequential data in Figure 1 (c). The sequence < (14)(10) > has been detected as Latent, 𝐾 2 is equal to 1. This sequence appears in the same way at time 1 and time 2. The status of a grid cell is set by considering the evolution type of all the sequential formal concepts that contain this grid cell in its intent and by applying a majority vote. Let 𝑆𝑇 be the set of status labels: 𝑆𝑇 ={𝐿𝑎𝑡𝑒𝑛𝑡 (LA), 𝐸𝑚𝑒𝑟𝑔𝑖𝑛𝑔 (E), 𝐷𝑒𝑐𝑟𝑒𝑎𝑠𝑖𝑛𝑔 (D), 𝐽𝑢𝑚𝑝𝑖𝑛𝑔 (J), 𝐿𝑜𝑠𝑡 (LO), 𝑁𝑜𝑡ℎ𝑖𝑛𝑔 (∅)}. For each grid cell 𝑔 𝑖 , we determine the status vector 𝑉 𝑔 𝑖 =< 𝑆𝑇 𝑡 𝑚𝑖𝑛 , 𝑆𝑇 𝑡 2 , ...., 𝑆𝑇 𝑡 𝑚𝑎𝑥 > where 𝑆𝑇 𝑡 𝑗 is the detected status for 𝑔 𝑖 at time 𝑡 𝑗 , and [𝑡 𝑚𝑖𝑛 , 𝑡 𝑚𝑎𝑥 ] is the studied time interval. This vector expresses a sequence of changing states over a discrete time. Figure 4 presents an example to determine status vectors. The status vector of 𝑔 13 is 𝑉 𝑔 𝑖 =< 𝐽, 𝐸, 𝐸, 𝐸, ., ., . >. For 𝑡 3 , there are three detected status: 𝐸, 𝐸, and 𝐿𝐴. A majority voting process is needed and gives 𝐸 as result.

Status vectors are used in the next step of our method as input to train the proposed Markov model.

Modeling Traffic Status by Markov Chains

To predict traffic status, we propose a method based on Markov model representations. Our method uses a specific area and a time window as inputs for studying each grid cell pattern trends. The proposed algorithm that builds the model is composed of three main steps. The first step leads to a graph, where nodes are status of a grid cell and edges represent the transition between two status. The strength of a given transition between nodes A and B is represented by the parameter 𝜆 𝑎𝑏 . The role of the second step is to compute values of all parameters 𝜆. The last step computes the transits state probability vector for an initial state. More details about these three steps are given below.

Step 1: Considering the set of status 𝑆𝑇 (see Section 4.1), we construct the structure of the model presented as a graph in Figure 5, which represents the Discrete-Time Markov Chain model for each grid cell. The corresponding transition matrix (𝑛 × 𝑛 with 𝑛=6) is defined in Equation 6. Let us remark that the real-world modeling leads to ignore some transitions. For instance, there is no transition from LO to D because it is not possible to switch from Lost status (LO) to Decreasing status (D).

E D LA J LO ∅ P = E D LA J LO ∅             𝜆 𝑒𝑒 𝜆 𝑒𝑑 𝜆 𝑒𝑙𝑎 0 𝜆 𝑒𝑙𝑜 𝜆 𝑒 ∅ 𝜆 𝑑𝑒 𝜆 𝑑𝑑 𝜆 𝑑𝑙𝑎 0 𝜆 𝑑𝑙𝑜 𝜆 𝑑 ∅ 𝜆 𝑙𝑎𝑒 𝜆 𝑙𝑎𝑑 𝜆 𝑙𝑎𝑙𝑎 0 𝜆 𝑙𝑎𝑙𝑜 𝜆 𝑙𝑎 ∅ 𝜆 𝑗𝑒 𝜆 𝑗𝑑 𝜆 𝑗𝑙𝑎 0 𝜆 𝑗𝑙𝑜 𝜆 𝑗 ∅ 𝜆 𝑙𝑜𝑒 0 0 𝜆 𝑙𝑜 𝑗 0 𝜆 𝑙𝑜 ∅ 𝜆 ∅𝑒 0 𝜆 ∅𝑙𝑎 𝜆 ∅ 𝑗 0 𝜆 ∅∅             (6)

Traffic status regions vectors

Figure 4: Example of grid computed status vector.

We will use two indexes 𝑖 and 𝑗 for accessing to a transition value 𝜆 between the nodes 𝐴 and 𝐵 corresponding to 𝑖 and 𝑗, respectively. For example, 𝑝 12 corresponds to 𝜆 𝑒𝑑 (the transition from E to D).

Step 2: During this second step, we compute the transition probabilities 𝜆. A discrete-time Markov chain is a process whose state space is finite with time interval indexes as 𝑡 = (1, 2, ...). The generated discrete time vector of status for each grid cell (see Figure 4) has randomness and the probability to have a certain state that depends only on the current state and not the whole previous states in the generated vector. 𝑃 is the transition probability matrix where each 𝑝 𝑖 𝑗 representing the probability of moving from state 𝑖 to state 𝑗. The matrix elements are not negative, and 𝑛 𝑗=1 𝑝 𝑖 𝑗 = 1 ∀𝑖. Thus, the probability 𝑝 𝑖 𝑗 is given by Equation 7.

𝑝 𝑖 𝑗 = 𝑓 𝑖 𝑗 𝑛 𝑘=1 𝑓 𝑖𝑘 (7)
Where 𝑓 𝑖 𝑗 is the number of occurrences of the transitions from 𝑖 to 𝑗 observed in the data, and 𝑛 𝑘=1 𝑓 𝑖𝑘 is the sum of the frequencies of 𝑖 to all other states 𝑘.

In our context, a transition represents a change in traffic status from one time interval to another. For each transition, a probability is associated which is the ratio of the number of changes from state 𝑖 to state 𝑗 divided by the total number of changes from state 𝑖.

Step 3: In the last step, we calculate the transits state probability vector for a given initial state 𝑠 0 which is defined according to the requested time to be predicted. Our assumption is that the detected traffic status at some point in the future (t+1) can be determined as a function of current traffic status (t), more precisely: 𝑝 𝑖 𝑗 = P(𝑉 𝑛+1 = 𝑗 | 𝑉 𝑛 = 𝑖), which is one step transition probability.

Given a positive integer 𝑚, the m-step transition probability

𝑝 (𝑚)
𝑖 𝑗 is the probability to reach the state 𝑗 after 𝑚 steps starting from the state 𝑖, 𝑝

𝑖 𝑗 = P(𝑉 𝑛+𝑚 = 𝑗 | 𝑉 𝑛 = 𝑖). In this work, we are interested in computing the traffic status after 𝑚 steps. Thus, the uniformization method (see Equation 3 in Section 2.4) is used to solve the model and obtain the vector of transient probability.

Let us consider the example of the grid cell 𝑔 13 (see Figure 4). The vector 𝑉 𝑔 13 =< 𝐽, 𝐸, 𝐸, 𝐸, 𝐿𝐴, 𝐷, ∅ > contains the observed status for each time 𝑡. From this vector, the transition probabilities are computed as described in Step 2. Starting from 𝑠 0 =𝐽 in 𝑡 1 , the vector 𝜋 is computed by applying the power method and uniformization method (there are 16 iterations). The maximal value in this vector is obtained for the transition 𝑗 → 𝑒 (the value is 0.522977). Thus, the predicted status for the next time 𝑡 2 is 𝐸 (i.e., Emerging).

EXPERIMENTS

This section presents the experiments performed on two real-world data sets. After the description of the data sets, we describe the experimental protocol and discuss the results.

Data

In the experiments, we have used two real-world taxi mobility data sets from two cities: Beijing (China) and San Francisco (USA). Table 2 provides a summary of those data sets.

T-Drive data set contains a large amount of taxi GPS trajectories collected in Beijing by Microsoft Research Asia [START_REF] Yuan | Tdrive: Driving directions based on taxi trajectories[END_REF]. The original data set contains the trajectories of 10,357 taxis during one week from 2 February 2008 to 8 February 2008. The total number of points is about 15 million and the total distance of the trajectories reaches 9 million kilometers. The GPS coordinates were returned every 5 to The second data set is Taxi cabs San Francisco which is available from the CRAWDAD website [START_REF] Piorkowski | CRAWDAD dataset epfl/mobility (v. 2009-02-24[END_REF]. San Francisco data set includes 500 taxis trajectories over 24 days during the period from 17 May 2008 to 10 June 2008 in the San Francisco Bay Area. This data set contains approximately 11 million taxi GPS samples, with a median time gap between two consecutive GPS measures of 60 seconds. After a preprocessing step of the data, we have obtained 455 active taxis.

Protocol

The following experimental protocol has been applied to each data set. First, we started by preprocessing all the trajectory data. The trajectories were segmented according to the time granularity value set to 12 hours (AM, PM). The study areas were divided into regions by setting the spatial granularity value to 15, 30, 60, or 120 meters.

After this preprocessing step, we detected the traffic status of each region between each pair of contiguous time windows (see Section 4.1). For this, we set 𝑚𝑖𝑛𝑠𝑢𝑝 to 10 for T-Drive and 𝑚𝑖𝑛𝑠𝑢𝑝 to 45 for San Francisco data set, 𝜃 = 1, and 𝜖 = 0. Let us remark For the evaluation step, we have compared the predicted status to the real status. We have computed the local accuracy values, the accuracy values for each traffic status and the global average accuracy values. We have compared our proposition to two other approaches to detect traffic status: a baseline approach that does not consider patterns (only singletons), and an itemset approach that does not consider the order of items. In each case, we have used our Markov model and only the way to obtain the traffic status are different. So, we evaluated the gain of using our approach based on sequential patterns.

We have implemented our proposition mainly with JAVA. We have used the SPMF library [START_REF] Fournier-Viger | The SPMF Open-Source Data Mining Library Version 2[END_REF] (developed in JAVA) for computing patterns, and the PEPS software tool [START_REF] Brenner | PEPS 2007 -Stochastic Automata Networks Software Tool[END_REF] (developed in C++) for solving Markov models. All the experiments were performed on an Intel Xeon X5560 2.8GHz with 16GB of memory.

Results

Tables 3 and4 present the global average accuracy obtained by applying our proposition (noted SP4TP seq.) and the two other methods: Baseline and Itemsets (noted SP4TP it.) on T-Drive and San Francisco data sets. We can observe that if the spatial granularity value decreases, the accuracy tends to increase for SP4TP (it. and seq.). Nevertheless, the spatial granularity value can not be too low because there will be some difficulties to detect patterns (high number of regions, execution time, memory space). For the Baseline approach, if the spatial granularity value decreases, the accuracy tends to decrease. The absence of relations between grid cells (i.e., patterns) is the explanation. For T-Drive (Beijing), our proposition outperforms the frequent pattern and the baseline approach for all different spatial granularity values. For San Francisco data set, our proposition has obtained best results for low spatial granularity values but the gain is not very significant. We can explain that by the relative quality of the collected data of San Francisco, especially the high time between two GPS measures (this can be a problem to correctly detect patterns).

Tables 5 and6 present the average accuracy values per status obtained with T-Drive and San Francisco data set, respectively. As we can see, the Lost and Jumping status are particular. They are very few and well predicted (almost no errors). We obtained very high values. We can say that, here, these status are not very significant.

Let us focus on Emerging, Decreasing and Latent. SP4TP (especially with sequences) outperforms the Baseline method which has its worst results. We can remark that the Baseline method has obtained the best result for the status Nothin but this status is not very informative and useful for user applications. For example, Figure 6 presents a tagged maps of Beijing. Regions having the Nothing status are not colored because we cannot conclude on them. To analyse the traffic or to establish a travel, a user or an application will naturally lean on the other status (Emerging, Jumping, . . . ).

CONCLUSION

In this paper, we have proposed a new method for predicting the traffic status in the next future. This method, called SP4TP, leverages closed sequential patterns, emerging patterns and Markov chains. The experimental results on real-world data sets from Beijing and San Francisco have shown that our proposition based on sequential patterns (via sequential formal concepts) outperforms the alternative approaches using itemsets or not using patterns.

In future work, we will perform other experiments with more data sets and compare with other existing systems. Furthermore, we will improve and extend our model to encode the relations between regions (for instance, related to road networks). For this, we can use the SAN (Stochastic Automata Networks) formalism [START_REF] Plateau | Stochastic Automata Networks[END_REF] which are able to model the system by the composition of a set of sub-system and it is equivalent to Markov Chains.

  (a) Moving object trajectories. (b) Mapping on raster area.

Figure 1 :

 1 Figure 1: From trajectories to sequences.

  𝑋 ) ≥ 𝑚𝑖𝑛𝑔𝑟 . An important class of emerging patterns are jumping emerging patterns which correspond to a subset of patterns that is present in one class and absent from the other. If 𝐺𝑅 𝑗,𝑖 (𝑋 ) = +∞, 𝑋 is said to be a jumping emerging pattern from D 𝑗 to D 𝑖 . Considering sequential patterns and the example of Figure 1 (c), we observe that 𝑠

  presents a simple transition diagram of a Markov model. The nodes represent the states of the model, and the arrows represent the transitions.

𝑒

Figure 2 :

 2 Figure 2: Markov chain example.

Figure 3 :

 3 Figure 3: Overview of the proposed method.

Figure 5 :

 5 Figure 5: Proposed model for each grid cell.

Figure 6 :

 6 Figure 6: An example of tagged map of Beijing (Blue=Latent, Green=Emerging, Red=Decreasing).

  3 , 𝑠 7 } ; {<(19)(23)(22)>}), ({𝑠 3 , 𝑠 7 , 𝑠 10 } ; {<(19) (23)>}), ({𝑠 4 , 𝑠 10 , 𝑠 12 } ; {<(14)>}), ({𝑠 4 , 𝑠 12 } ; {<(14)(10)>}), etc.

Table 1 :

 1 Some sequential formal concepts and their evolution types obtained from Figure1(c) with 𝑚𝑖𝑛𝑠𝑢𝑝 = 2, 𝜃 = 1, and

	𝜖 = 0
	Extent (O)	Intent (S)	𝐾 2 Evolution type
	{𝑠 3 , 𝑠 7 }	{< (19)(23)(22) >} 0	Lost
	{𝑠 3 , 𝑠 7 , 𝑠 10 }	{< (19)(23) >}	0.5	Decreasing
	{𝑠 4 , 𝑠 10 , 𝑠 12 }	{< (14) >}	2	Emerging
	{𝑠 4 , 𝑠 12 }	{< (14)(10) >}	1	Latent

1 t 2 t 3 t 4 t 5 t 6 t 7 13

												Majority
												voting
												filter
	T 7 ……	21 17 13	21 17	22 18 14 21 22 23 19 15 22 18 17 18 21	24 20 16 23 23 19 19 22	24 20 24 20 23	T 4	T 3	24 T 2	T 1	g i t Time interval
		E 9	13 E,E , LA	10 13 E	11 17 14 14	12 18 15 15 E	16 16 19	20	LA
			9			10 13		11	14	12	15	16
		5		6	9		7 10	J	8 11	12
			5			6	9			7	10	8	11	12
		1		2	5		6	3		4 E 7	8
			1			2	5			3	6	4	7	8
					1		2			3	J	4
							1				2	3	4

J E E E LA D Ɵ Time interval ={t 1 ,t 2 ,t 3 …..t 7 } ST = {E, D, LA, J, LO, Ɵ }

Table 2 :

 2 Taxi mobility data sets

	Characteristics	Beijing (T-Drive) San Francisco
	Measurement	GPS	GPS
	Number of samples	15 million	11 million
	Duration	1 week	24 days
	Number of taxis	10,357	500
	10 seconds. In order to perform our experiments, we have chosen
	222 taxis randomly.

Table 3 :

 3 Global average accuracy for T-Drive (Beijing) Data We present here the best obtained results. Considering the previous detected traffic status as data, the Markov model has been built for each region (see Section 4.2). The proposed method have been applied by changing the start state 𝑠 0 during a time interval of 7 days for T-drive, and 24 days in San Francisco data. From 𝑡 𝑖 , the traffic state of 𝑡 𝑗 has been predicted, 𝑖 = 1 . . . 6 and 𝑗 = 2 . . . 7 for T-Drive. 𝑖 = 1 . . . 23 and 𝑗 = 2 . . . 24 for San Francisco data.

	Spatial gr. Baseline SP4TP (it.) SP4TP (seq.)
	15 m	0.747	0.921	0.962
	30 m	0.785	0.906	0.942
	60 m	0.807	0.902	0.929
	120 m	0.831	0.900	0.920
	that the parameter setting is proper to each data set and requires
	some tests.

Table 4 :

 4 Global average accuracy for San Francisco Data Spatial gr. Baseline SP4TP (it.) SP4TP (seq.)

	15 m	0.8984	0.922	0.924
	30 m	0.901	0.909	0.910
	60 m	0.899	0.897	0.899
	120 m	0.903	0.902	0.909

Table 5 :

 5 Average accuracy per status for T-Drive (Beijing) Data

	Status	Baseline SP4TP (it.) SP4TP (seq.)
	Emerging	0.653	0.775	0.849
	Decreasing	0.754	0.936	0.945
	Latent	0.784	0.969	0.987
	Lost	0.880	1	0.999
	Jumping	0.754	0.987	0.986
	Nothing	0.929	0.778	0.863

Table 6 :

 6 Average accuracy per status for San Francisco Data

	Status	Baseline SP4TP (it.) SP4TP (seq.)
	Emerging	0.732	0.827	0.836
	Decreasing	0.729	0.790	0.797
	Latent	0.944	0.990	0.988
	Lost	0.997	1	1
	Jumping	0.997	1	1
	Nothing	0.999	0.837	0.841