Applying advanced machine learning models to classify electro-physiological activity of human brain for use in biometric identification

Iaroslav Omelianenko

NewGround LLC Kiev, Ukraine Email: yaric@newground.com.ua Abstract-In this article we present the results of our research related to the study of correlations between specific visual stimulation and the elicited brain's electro-physiological response collected by EEG sensors from a group of participants. We will look at how the various characteristics of visual stimulation affect the measured electro-physiological response of the brain and describe the optimal parameters found that elicit a steadystate visually evoked potential (SSVEP) in certain parts of the cerebral cortex where it can be reliably perceived by the electrode of the EEG device. After that, we continue with a description of the advanced machine learning pipeline model that can perform confident classification of the collected EEG data in order to (a) reliably distinguish signal from noise (about 85% validation score) and (b) reliably distinguish between EEG records collected from different human participants (about 80% validation score). Finally, we demonstrate that the proposed method works reliably even with an inexpensive (less than $100) consumer-grade EEG sensing device and with participants who do not have previous experience with EEG technology (EEG illiterate). All this in combination opens up broad prospects for the development of new types of consumer devices, [e.g.] based on virtual reality helmets or augmented reality glasses where EEG sensor can be easily integrated.

The proposed method can be used to improve an online user experience by providing [e.g.] password-less user identification for VR / AR applications. It can also find a more advanced application in intensive care units where collected EEG data can be used to classify the level of conscious awareness of patients during anesthesia or to automatically detect hardware failures by classifying the input signal as noise.

Keywords-biometric identification; electro-physiological brain activity; machine learning; auto-encoder; visually evoked potential; SSVEP; EEG

I. INTRODUCTION

A steady increase in the computing power of mobile devices allows them to be fully integrated with a variety of portable or wearable consumer devices such as virtual reality (VR) helmets, augmented reality (AR) glasses, monitors of the health and fitness, and so on. At the same time, they have enough computational resources to execute advanced machine learning models for processing of biometric data received from various sensors in the real time. One of these non-intrusive sensors that can be used to collect biometric information is the sensor to monitor the electro-physiological activity of the brain -electroencephalography (EEG). To collect EEG data it's enough to place electrodes along the scalp of the user. Moreover, the EEG sensors can be easily integrated into VR helmets or AR glasses, also providing easy means for specific visual stimulation. As has been demonstrated in recent studies [START_REF] Ding | Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency[END_REF], [START_REF] Regan | Steady-state evoked potentials[END_REF], exposing a person to external stimuli, such as aural tone or a specific visualization, can cause a specific response in the electro-physiological activity of the brain. And this stationary cerebral cortex response, called steady-state visually evoked potentials (SSVEP) [START_REF] Regan | Steady-state evoked potentials[END_REF], demonstrates a strong correlation with the parameters of visual stimuli, such as flicker frequency and visual appearance (configuration).

In [START_REF] Ding | Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency[END_REF], it was found that a reliable SSVEP response can be detected in the Delta and Alpha (Low and High) frequency bands of the EEG signal collected at stimulation with a frequency of visual stimuli in the range: 4 -10Hz. That is, visual stimuli with a certain frequency and configuration can elicit a reliable SSVEP response in certain parts of the cerebral cortex, which leads to a resonance effect that can be easily detected by EEG sensors placed in specific positions on the scalp of the participant's head. (This can be partially explained if we consider neural networks of the cerebral cortex as a structural dynamic network of oscillators [START_REF] Tang | Control of Dynamics in Brain Networks[END_REF] which can be entrained by specific stimulus; where different regions of cerebral cortex having specific resonance frequencies amplifying its response on stimulation.) Later, the amplified response detected by EEG sensors can be processed using advanced machine learning classification methods to extract useful patterns.

In this work, we use a consumer-grade EEG sensing device with a single electrode located above the FP1 position, which collects the electrical activity of the frontal and pre-frontal cerebral cortex. Because of this limitation, we are most interested in finding specific visual stimuli capable of eliciting a reliable SVEP response in these areas. From the work of other researchers [START_REF] Ding | Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency[END_REF], we know that the frontal cortex has a strong response in the High Alpha band, so we focused our efforts on creating a main visual stimulus with a frequency in the range of 8-10Hz and auxiliary visual stimulus with the frequency of one harmonic down. In [START_REF] Ding | Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency[END_REF] it was found that with simultaneous use of several visual stimuli with frequencies related on the harmonic scale, the SSVEP response of certain parts of neural networks of the cerebral cortex can be amplified.

The purpose of this study is to create an advanced machine learning classification model suitable for processing of collected EEG records and capable of revealing the differences:

-between signal and noise, -between EEG records collected from different human participants.

This paper is organized as follows: In Section II, we describe the raw input signals received from the EEG device, visual stimulation details, and basic information about the involved participants. It is followed in Section III by details about collected data records and how it was pre-processed before the analysis. In Section IV, we provide description of machine learning classification model applied to the data corpus. Then, the Section V follows, describing the methodology used to find the optimal configuration of visual stimulation and to find optimal hyper-parameters of the machine learning classification model. After that in Section VI we consider achieved results (validation scores) for both signal / noise and participants identification tasks. Finally, in Section VII and VIII we describe our plans for future experiments and review the results obtained during the research.

II. EXPERIMENT SETUP DESCRIPTION

A. EEG monitoring device and raw input signals

In our experiments, we use NeuroSky MindWave device that monitors the brain's electro-physiological activity using one electrode at the FP1 position and other at the ear lobe. It has a sampling rate of 512 Hz and decomposes the raw input signal on several frequency bands at the hardware level using a fast Fourier transform (FFT) [START_REF] Heideman | Gauss and the history of the fast Fourier transform[END_REF]. After FFT decomposition, the MindWave device broadcasts the decomposed EEG signal once per second (i.e., the refresh rate is 1Hz). See Table I We use the mentioned complex synthetic signals as control signals that trigger the recording / processing of EEG data when a certain threshold value is reached.

B. Visual stimulation and data collection

Following research [START_REF] Ding | Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency[END_REF], we created an advanced software library to produce a certain visual stimulation. The visualization library provides tools for configuring various visualization options and applying various visual modeling schemes: WISP and WAVE. This allows a series of experiments with different parameters to be performed to find the optimal parameters for eliciting proper SSVEP response.

The visual stimuli created by the visualization library have a strong correlation with the collected EEG data in real time, providing positive feedback to the user on the activity of the cerebral cortex. In the WISP rendering mode, this correlation is represented as the relative position of points in polar coordinates and as the color of points and wisp lines (see Figure 1). The position of the points relative to the center, determined by the strength of the EEG signal for a particular frequency band. Each frequency band has a certain color assigned to it. To elicit appropriate SSVEP response the visualization has two flickering areas:

1) wisp of lines with dots -flickers with selected primary frequency 2) background -flickers with secondary frequency (one harmonic down)

The background color saturation is determined by the current value of selected control signal -Attention or Meditation, where higher values result in a more saturated color. Thus, the saturation of the background color provides positive feedback for the participant about current value of the control signal.

The central idea of the described scheme of visual stimulation is to achieve aesthetically pleasing and attractive visualization, which can be used in consumer devices: helmets VR, AR glasses, etc. But at the same time, it should provide appropriate visual stimulation to elicit proper SSVEP response from cerebral cortex.

Another studied mode of visual stimulation, we called WAVE. Its main difference from WISP is that it displays only dots, skipping lines that connect them to the center point (see Figure 2). We found that different visual modes have a variable effect on elicited SSVEP response in different frequency ranges of stimuli. In the WISP rendering mode the most effective primary stimulation frequency is 10Hz (High Alpha) and secondary -5Hz. And for the WAVE this is: primary -8Hz (Low Alpha) and secondary -4Hz. The graphic interface of the created visualization library also allows you to set a certain control signal threshold that starts / stops recording the received EEG signals. We learned that with higher threshold values, it is much more difficult for a participant to maintain that level of control signal (Attention / Meditation), but the collected data is of better quality. Thus, it is important to find the right balance between the quality of the collected data and ease of use for the participants. In our experiments, we applied four threshold values from 80 down to 50 in steps of ten. The optimal threshold value for the WISP rendering mode is 70 providing a good balance between the quality of the collected data and ease of use for the participant.

We measure quality of the collected data and the applicability of a specific configuration of visual stimuli, by checking Pearson product-moment correlation [START_REF] Gain | The frequency distribution of the product moment correlation coefficient in random samples of any size draw from nonnormal universes[END_REF] between the processed data from the recorded sessions for each participant per each configuration. (For more information on calculating correlation coefficients, see Section: " Finding Optimal Configuration for Visual Stimulation".) And for further analysis, the configuration of visual stimuli creating the most correlated results is considered. Our hypothesis for this is that the maximum correlation between recording sessions is the result of properly elicit SSVEP response from a specific parts of the cerebral cortex that has a steadily recognizable pattern among all recorded sessions.

Applying the above-mentioned method to estimate quality of the collected EEG data, following visualization parameters were chosen for all participants as optimal:

-WISP rendering mode -10Hz primary and 5Hz secondary flicker frequencies -Attention complex synthetic signal as a control signal -70 as the control signal threshold value In a series of experiments, we found that reliable data processing requires at least forty collected data samples (EEG events) per session. The MindWave device broadcasts data event signals (with 8 frequency bands values) once per second. This means that one recording session should take at least 40 seconds, if participant can maintain a continuous high-level value of the control signal. But this can take much longer if the value of the control signal drops below the threshold and the recording pauses until participant can again reach the threshold value.

C. Experiment Participants

In our experiments we collected EEG data records from three participants: 1) male, 43 years, -1 myopia, astigmatism 2) female 28 years, normal sight 3) male 35 years, -7 myopia For each participant, we collected from six to eighteen sessions of EEG records with various schemes of visual stimulation and at different times. It is very important that EEG sessions be recorded at different times in order to observe the electro-physiological response of the participants' brain in different conditions.

III. DATA CORPUS DESCRIPTION

The data collected from the EEG device during a particular recording session is stored as comma-separated records with 40 rows, 11 columns wide. Each line contains three metadata variables (time stamp of record, frequency of primary and secondary stimuli) and eight values for each frequency band.

The frequency bands data collected from EEG device represent a power spectrum and its values vary exponentially in the range [0, 32767]. But performance of a majority of machine learning methods is best with small floating point values centered around zero. Thus, in order to improve the quality of the collected data and to balance the values of the data points for different frequency bands, we scaled it to fit in range [0, 1] during the preprocessing. And the preprocessed data is later used as input for the machine learning classification model.

IV. MACHINE LEARNING CLASSIFICATION PIPELINE MODEL DESCRIPTION

Essentially, the collected data corpus of EEG records is a time series. Therefore, it may seem natural to analyze it using machine learning methods, which are usually considered to have the best architecture for time series data analysis, e.g. Recurrent Neural Networks (RNN) and its modern reincarnation -Long Short Time Memory Neural Network (LSTM) [START_REF] Hochreiter | Long Short-Term Memory[END_REF]. But for RNN model to learn useful patterns from input data at least two conditions should be met: (a) the data samples should be collected at equal time intervals between events; (b) it is important to collect a significant number of data samples for each recording session. In our experiment, no above condition can be fulfilled. First of all, it is unpractical to collect a lot of data samples for consumer device setup, because people easily get bored during long recording sessions. Furthermore, we collect EEG data events only when the level of the control signal exceeds a certain threshold, which makes our time series rather inconsistent with different intervals between time stamps of collected events.

Taking into account the specific conditions of our experiment mentioned above, we decided to explore a novel approach Fig. 3: The machine learning processing pipeline scheme from left to right: the EEG input data source, the auto-encoder, and the classifier.

to the analysis of time series, which allows us to extract useful patterns (essence) from the analyzed EEG signal records and classify it later. Our idea is to create advanced machine learning classification pipeline, where the input EEG data is processed in two stages (see Figure 3):

-at the first stage, raw input data are analyzed using unsupervised machine learning method based on the auto-encoder neural network with single hidden layer [START_REF] Bengio | Learning Deep Architectures for AI[END_REF] -at the second stage, the learned weights of the autoencoder hidden layer (data encoding) are used as input for the classifier of choice to train the supervised classification model

In the following we describe considered methods of machine learning used at both stages, with more details.

A. The First stage -Auto-Encoder Selection

In our experimental setup -per each individual sessionwe have the data corpus consisting of eight frequency bands values collected for forty data events. That gives us the features space dimension width of 320 for each EEG recording session. At the same time, we have very limited number of recording sessions (6 to 18) for every participant, which is much less than the dimensionality of the features space. Such characteristics of the input data make it unacceptable for analysis by almost any machine learning classification model, since ideally the number of data samples must significantly exceed the number of features and not vice versa. This forced us to look for a specific method of machine learning that can extract its reliable representation (encoding) from the raw EEG input signal and accordingly reduce the dimensionality of the features space.

Our attention was captured by family of unsupervised machine learning methods called auto-encoders [START_REF] Bengio | Learning Deep Architectures for AI[END_REF]. An autoencoder takes an input x ∈ [0, 1] d and first maps it (with an encoder) to a hidden representation h ∈ [0, 1] d through a deterministic mapping, e.g.:

h = σ(Wx + b) (1)
where σ is a non-linearity such as the sigmoid.

The latent representation h, or encoding is then mapped back with a decoder into a reconstruction z of the same shape as x. The mapping happens through a similar transformation, e.g.:

z = σ(W h + b ) (2)
z should be seen as a prediction of x, given the code h.

Optionally, the weight matrix W of the reverse mapping may be constrained to be the transpose of the forward mapping: W = W T . This is referred to as tied weights. The parameters of this model (W, b, b and, if one doesn't use tied weights, also W ) are optimized such that the average reconstruction error is minimized.

The reconstruction error can be measured in many ways, depending on the appropriate distributional assumptions on the input given the code. The traditional squared error L(xz) = ||x -z|| 2 , can be used as well. If the input is interpreted as either bit vectors or vectors of bit probabilities, the crossentropy of the reconstruction can be used:

L H (x, z) = - d k=1 [x k log z k + (1 -x k ) log(1 -z k )] (3)
The h can be viewed as lossy compression of x and optimization makes it a good compression for training examples and hopefully for all other examples as well. But this can not be guaranteed for arbitrary inputs. The simple auto-encoder gives low reconstruction error on test samples from the same distribution as train samples, but usually gives higher error values for inputs randomly chosen from input space.

To prevent auto-encoder from learning identity function of input signal and to force it to learn something useful about input signal in its hidden units, several methodologies can be applied: (a) addition of sparsity to inputs, (b) addition of randomness to transformation (encoding) or (c) introducing specific explicit regularizer to the objective function.

In this research we study two types of auto-encoders with various learning enhancements applied:

-Denoising auto-encoder (addition of sparsity to inputs) -Contractive auto-encoder (the explicit regularizer to the objective function)

The denoising auto-encoder [START_REF] Vincent | Extracting and Composing Robust Features with Denoising Autoencoders[END_REF] is a stochastic version of simple auto-encoder where stochastic corruption process randomly sets many of inputs (maximum half of them) to zero. Then we train it to reconstruct correct input from a stochastically corrupted version of it. As result, the higher level representations of input data (encoding) are relatively stable and robust to the corruption of input data.

The contractive auto-encoder [START_REF] Rifai | Contractive Auto-Encoders: Explicit Invariance During Feature Extraction[END_REF] adds explicit regularizer to the objective function that forces model to obtain a robust representation of the input space not sensitive to slight variations of input data. This regularizer corresponds to the Frobenius norm (L2) of the Jacobian matrix of the hidden representation with respect to the input. The final objective function has the following form:

L = - d k=1 [x k log z k + (1 -x k ) log(1 -z k )] + λ d i=1 n j=1 J 2 ij (4) where z = σ(W h + b ) is the reconstruction of input vector, J i = h i (1 -h i ) • W i is
Jacobian of h with respect to x, and h = σ(Wx + b) is the projection of the input into the latent space h.

Both mentioned auto-encoder variations are robust to the slight corruption of input data. This is especially useful for our experiment, since it is well known that the electrophysiological activity of the brain is non-stationary and varies greatly, depending on the circumstances. I.e. the collected EEG data for each subsequent recording session under the same type of visual stimuli can give the same essence (compressed representation), but can introduce various non-stationary deviations from previous sessions. Using an auto-encoder to analyze the input data at the first stage, we hope to reduce non-stationary deviations and extract a useful compressed representation. Then, at the second stage, the extracted compressed representation with significantly reduced dimensionality of features space will be fed into the classifier of choice.

The above-mentioned auto-encoders was implemented in Python programming language using Theano framework [START_REF] Al-Rfou | A Python framework for fast computation of mathematical expressions[END_REF].

B. The Second Stage -Classifiers Evaluation

At the end of the first stage of the processing pipeline, the dimensionality of data corpus features space is considerably reduced. After that, at the second stage, the data is fed to a set of classifiers from simple to advanced in order to find the one that has the most powerful prediction model.

To automate the evaluation of classifiers against different sets of hyper-parameters, we use exhaustive grid search with cross-validation. This method exhaustively generates candidates from the grid of parameter values specified for each classifier, and then performs a 3-fold cross-validation (CV) of the predictive performance of each classifier with respect to the input data [START_REF] Kohavi | A study of cross-validation and bootstrap for accuracy estimation and model selection[END_REF]. In the k-fold cross-validation method, the training set is split into smaller sets and the following procedure is performed for each of the folds: 1) a model is trained using k -1 of the folds as training data 2) the resulting model is validated on the remaining part of the data (i.e., it is used as a test set to compute a performance measure such as accuracy)

The resulting classifier performance measure is done by averaging scores computed for each fold. Thus, we can avoid the effect of model overfitting, when it shows perfect score for train data but fails to predict anything useful for unseen data. Refer to Table II for a list of all tested classifiers and all tested values of corresponding hyper-parameters per each estimator.

Hereafter we provide short description of each tested classifier. [START_REF] Breiman | Random Forests[END_REF] is a meta estimator that fits a number of decision tree classifiers on various sub-samples of the dataset and use averaging to improve the predictive accuracy and control overfitting. [START_REF] Zhu | Multi-class AdaBoost[END_REF] is a meta-estimator that begins by fitting a classifier on the original dataset and then fits additional copies of the classifier on the same dataset but where the weights of incorrectly classified instances are adjusted such that subsequent classifiers focus more on difficult cases. [START_REF] Breiman | Classification and Regression Trees[END_REF] is a flow-chart-like structure, where each internal (non-leaf) node denotes a test on an Gaussian Process classifier [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] based on Laplace approximation, which is used for approximating the non-Gaussian posterior by a Gaussian. [START_REF] Hinton | Connectionist learning procedures[END_REF] is an artificial neural network based classifier optimizing log-loss function using LBFGS or stochastic gradient descent. [START_REF] Altman | An introduction to kernel and nearest-neighbor nonparametric regression[END_REF] use a type of instance-based learning or non-generalizing learning where classification is computed from a simple majority vote of the nearest neighbors of each point: a query point is assigned the data class which has the most representatives within the nearest neighbors of the point. [17] is a model assuming that the continuous values associated with each class are distributed according to a Gaussian distribution. [START_REF] Hastie | The Elements of Statistical Learning[END_REF] is a classifier with a quadratic decision boundary, generated by fitting class conditional densities to the data and using Bayes rule. [START_REF] Cortes | Support-vector networks[END_REF] builds a representation of the examples as points in space, mapped so that the examples of the separate categories are divided by a clear gap that is as wide as possible. New examples are then mapped into that same space and predicted to belong to a category based on on which side of the gap they fall.

Random Forest classifier

Ada Boost classifier

Decision Tree classifier

Multi-layer Perceptron (MLP) classifier

K-neighbors classifier

Gaussian Naive Bayes

Quadratic Discriminant Analysis

Support Vector Machine classifier

We implemented exhaustive grid search with 3-fold crossvalidation for all above mentioned classifiers in Python programming language using Scikit-learn framework [START_REF] Pedregosa | Scikitlearn: Machine Learning in Python[END_REF]. It can be seen that majority of sessions highly correlate with each other and only session #6 is a bit outlier. Such high correlation values is a good sign that applied configuration of the visual stimuli is optimal to elicit robust SSVEP response which can be compressed to extract useful pattern (encoding).

V. FINDING OPTIMAL EXPERIMENT CONFIGURATION

The details of the experiment described above suggest that in order to obtain reliable classification, it is important to find optimal configuration parameters for both parts of the experiment: visual stimulation and machine learning classification pipeline.

Hereafter we consider in more details how optimal configuration was found.

A. Finding Optimal Configuration for Visual Stimulation

To achieve the objectives of the experiment, first, we need to find the optimal configuration of visual stimuli that elicit a reliable SSVEP response in the cerebral cortex of participants. The optimal configuration should give comparable results for different EEG sessions recorded on different times. As mentioned in Subsection: "Visual stimulation and data collection" to do this we conducted series of experiments (recording sessions) and evaluate collected data by calculating Pearson product-moment correlation [START_REF] Gain | The frequency distribution of the product moment correlation coefficient in random samples of any size draw from nonnormal universes[END_REF] between outputs from the first stage of machine learning pipeline (auto-encoder). Our hypothesis is that the visual stimulation scheme leading to outputs with the maximal correlation can be considered optimal.

The correlation coefficients can be calculated as following:

1) The collected EEG raw input data preprocessed by scaling it down to fit range [0, 1] (see Section: "Data Corpus Description"). 2) After that, the scaled data corpus is fed into autoencoder of choice. The input data is filtered to include only Delta and High Alpha bands. These bands were chosen based on the results of [START_REF] Ding | Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency[END_REF], which state that the elicited SSVEP response should have bands frequencies close to the frequencies of the visual stimuli (10/5Hz). It was found that for best results the autoencoder should have only two units in the hidden layer.

3) The learned weight matrix of the auto-encoder's hidden layer (compressed input data representation) per each recording session then collected as representation data matrix with columns holding compressed data for each session. 4) Finally, the Pearson product-moment correlation coefficients calculated from representation matrix and heat map plot with correlations are rendered (see Figure 4).

B. Finding Optimal Configuration for Machine Learning Classification Pipeline Model

The optimal configuration of classification pipeline was found in a series of experiments by applying various hyperparameters to auto-encoders and classifiers. After that, the classification model with the maximum validation score was chosen as optimal. To automate the process of selecting the optimal model, a software framework was created that performs the following procedures:

1) starts preprocessor to process the raw input data corpus consisting of collected EEG recording sessions, 2) uses auto-encoder against preprocessed input data to extract useful patterns, 3) perform classifiers evaluation against extracted patterns as described in Section: "The Second Stage -Classifiers Evaluation".

The experimental framework allows to separately define various configuration parameters for each processing stage. Thus, the configuration of each experiment can be encapsulated in a separate Python script, which makes it easy to reproduce a specific experiment in the future.

We tested a variety of hyper-parameters combinations per each auto-encoder with different sets of filtered input data samples (only values for specific frequency bands are included). As result, we found that learning rate, contraction / corruption level and the number of learning epochs can be fixed by the values presented in Table III. Only the batch size, the number of auto-encoder's hidden units and frequency bands included into input data corpus, can be changed when searching for the best validation score.

VI. RESULTS

In this paper, we considered the possibility of creating reliable machine learning classification pipeline model able to:

1) confidently distinguish between signal and noise in collected EEG records 2) confidently classify EEG records of different participants, that is, provide means to determine which EEG record belongs to which participant.

Hereafter we consider achieved results in more details.

A. Signal / Noise Classification

For this task our main goal is to build classifier capable of distinguishing between signal and noise records collected under different conditions and at different times. The signal records include SSVEP affected records plus EEG records obtained from participants performing other cognitive tasks (sky watching, reading, etc). The basic requirement for the proper signal record -is to be collected when the value of the control signal (Attention or Meditation) is greater than or equal to the specified threshold value. The noise records are collected when the electrodes of the EEG device are improperly installed on the participants (for example, the base electrode is not connected) or when the device is turned on but not connected to the participant at all. At total, we studied data samples collected over 26 signal and 26 noise sessions. We began experiments with finding a best validation score for a specific set of auto-encoder's hyper-parameters when only two frequency bands were included in the input data corpus and continued until all bands were involved. We found that in all experiments the classification pipeline based on contractive auto-encoder was superior to that based on the denoising autoencoder. Therefore, we will provide only the best results for the contractive auto-encoder together with the classifier name and the included frequency bands (see Table IV).

The best signal / noise classification score (0.849) obtained for the machine learning classification pipeline, including the contractive auto-encoder, followed by the classifier based on multi-layer perceptron architecture. See Table III, IV for used auto-encoder hyper-parameters. Optimum hyper-parameters of the most effective classifier (MLP) -momentum: 0.999, activation: 'relu', solver: 'sgd', alpha: 0.0001, learning rate init: 0.001, hidden layer sizes: 8

As a result of the experiments, we found that for the signal / noise classification the SSVEP response elicited by visual stimulation is not important. Only the level of the control signal (Attention or Meditation) during the recording of EEG events is of primary importance. The best result was achieved when almost all frequency bands is included in the analyzed data corpus, which confirms that the elicited by visual simulation SSVEP response is not important for this task. Thus, we can conclude that the proposed method of signal / noise classification can be used without any visual stimulation. The scope of this method application can be related both to the monitoring of EEG hardware failures and to control the user's level of consciousness in intensive care units.

B. The Classification of Participants

The second goal of our research is to build a model of machine learning classification pipeline that allows correctly labeling EEG records for specific participants, that is, allowing to classify a participant on the basis of her unique fingerprint in the collected EEG data.

In our experiments we decided to apply visual stimulation with main flicker frequency at 10Hz and secondary -at 5Hz, which corresponds to Theta (Θ) and High Alpha (Hα) bands. In accordance with [START_REF] Ding | Attentional modulation of SSVEP power depends on the network tagged by the flicker frequency[END_REF], we assumed that elicited SSVEP response should be detected at these or adjacent frequency bands and that the best validation score will be obtained if only these frequency bands are included in the input data corpus.

A series of experiments was performed on the collected signal-only EEG records with various frequency bands included and the best result (0.800) was achieved with data corpus consisting of signal-only records with Delta (∆) and High Alpha (Hα) bands included. It almost perfectly fits into the expected range -the Delta (∆) band is adjacent to the Theta (Θ) band. This congruence is a good indication that SSVEP stimulation is important when collecting EEG records to be classified by participants. See Table V for details about the best results found for various data corpus and hyper-parameters configurations.

VII. FUTURE WORKS

Our current research based on very limited data corpus collected from only three participants. Because of this limi- tation, we may have missed some important patterns in EEG records, which can further improve the validation scores of studied machine learning classification pipeline methods. In the future, we plan to collect and process EEG records from a much wider audience of participants.

VIII. CONCLUSION

In this paper, we demonstrated how to build machine learning classification pipeline capable of processing and reliably classifying EEG records collected from different participants. We described a visual stimulation scheme suitable to elicit specific SSVEP response from visual cortex and hypothesized how to measure its quality using the collected EEG recordings. In addition, special attention was paid to creating an aesthetically attractive scheme of visual stimulation in the assumption of its use in consumer devices.

We managed to get a fairly good estimate of confidence score with signal / noise classification task (best validation score: 0.849, see Table IV). This confirms that the proposed method is suitable for automatic processing of EEG records, when it is important to distinguish the signal from noise.

Fig. 1 :

 1 Fig. 1: The screen shoot of visual stimulator window in WISP rendering mode.

Fig. 2 :

 2 Fig. 2: The screen shoot of visual stimulator window in WAVE rendering mode.

Fig. 4 :

 4 Fig.4: The heat map with correlation coefficients between eight recording sessions for one of the experiment participants. It can be seen that majority of sessions highly correlate with each other and only session #6 is a bit outlier. Such high correlation values is a good sign that applied configuration of the visual stimuli is optimal to elicit robust SSVEP response which can be compressed to extract useful pattern (encoding).

TABLE I

 I

	: The ranges of frequency bands for FFT decomposed
	EEG signal
	Band name	Frequency range, Hz
	Delta, ∆	1-3
	Theta, Θ	4-7
	Low Alpha, Lα	8-9
	High Alpha, Hα	10-12
	Low Beta, Lβ	13-17
	High Beta, Hβ	18-30
	Low Gamma, Lγ	31-40
	High Gamma, Hγ	41-50
	Additionally, it provides two complex synthetic signals
	calculated by combining data from several frequency bands,
	namely:
	-Attention to indicate users attention level in range
	[1, 100]
	-Meditation to indicate users contemplation level in
	range [1, 100]

TABLE II :

 II The classifiers used in the study with range of hyper-parameter values applied for exhaustive grid search.

	Classifier	Parameters space
	Random Forest	max depth: [3, 5, 8, None], max features: [1, 3, 4],
		min samples split: [2, 3, 4], min samples leaf: [1,
		3, 10], bootstrap: [True, False], criterion: ["gini",
		entropy], n estimators: [10, 20, 50, 100, 200]
	Ada Boost	learning rate: [0.01, 0.1, 1], n estimators: [10, 20,
		50, 100]
	Decision Tree	max depth: [3, 5, 8, None], min samples split: [2,
		3, 4], min samples leaf: [1, 3, 10, 20, 30]
	Gaussian Process	warm start: [True, False]
	Multi Layer Perceptron	alpha: [0.0001, 0.001, 0.01], learning rate init:
		[0.001, 0.01, 0.1, 0.5], momentum: [0.9, 0.99,
		0.999], solver: ["lbfgs", "sgd", adam"], activation:
		["logistic", "tanh", relu"], hidden layer sizes: [4,
		8, 10]
	K-neighbors	n neighbors: [2, 3, 5], algorithm: ["ball tree",
		"kd tree", "brute"]
	Gaussian Nave Bayes	priors: [None]
	Quadratic Discriminant	priors: [None], reg param: [0.0, 0.01, 0.1, 0.9]
	Analysis
	Support Vector Machine	C: [0.1, 0.5, 1.0], gamma: [0.1, 0.5, 1.0, 2.0, 3.0,
	with RBF kernel	'auto']
	Support Vector Machine	C: [0.001, 0.005, 0.01, 0.025, 0.05, 0.1, 0.25, 0.5,
	with Linear kernel	1.0]

attribute, each branch represents the outcome of a test, and each leaf (or terminal) node holds a class label. This structure allows to go from observations (branches) to the conclusions about the item's target value (leaves).

TABLE III :

 III The optimal hyper-parameters for auto-encoders for signal / noise classification task.

	Hyper-parameter	Tested values range	Contractive auto-encoder	Denoising auto-encoder
	Batch size	[1, 10]	1 or 10	5
	Learning rate	[0.001, 0.1]	0.1	0.1
	Number of units in hid-den layer	[2, 10]	5	5
	Number of learning epochs	[10 4 , 5 • 10 4 ]	5 • 10 4	5 • 10 4
	Contraction level (con-tractive auto-encoder)	[0.1, 0.3]	0.1	-
	Corruption level (denois-ing auto-encoder)	[0.1, 0.3]	-	0.1
	Included frequency bands	∆, Θ, Lα, Hα, Lβ, Hβ, Lγ, Hγ	∆, Θ, Lα, Hα, Lβ, Hβ	∆, Θ, Lα, Hα, Lβ, Hβ

TABLE IV :

 IV The comparison of the best validation scores for signal / noise classification per particular set of included frequency bands and specific values of contractive auto-encoder hyper-parameters. The score value is given as the average of 3fold cross validation with standard deviation among individual values for each fold.

	Frequency bands included	Batch size	# of hid-den units	Best classifier	Validation score
	∆, Hα	5	2	Multi Perceptron Layer	0.736 (std: 0.184)
	∆, Lα, Hα	5	4	Multi Perceptron Layer	0.736 (std: 0.050)
	∆, Θ, Lα, Hα	5	3	SVM with lin-ear kernel	0.755 (std: 0.034)
	∆, Θ, Lα, Hα	1	2	Gaussian Nave Bayes	0.833 (std: 0.068)
	∆, Θ, Lα, Hα, Lβ	5	4	Multi Perceptron Layer	0.774 (std: 0.098)
	∆, Θ, Lα, Hα, Lβ, Hβ	10	5	Multi Perceptron Layer	0.849 (std: 0.058)
	∆, Θ, Lα,
	Hα, Lβ, Hβ,	5	9	Random Forest	0.792 (std: 0.075)
	Lγ

TABLE V :

 V Comparison of the best validation scores for signalonly classification per particular set of included frequency bands and specific values of hyper parameters of contractive auto-encoder. The score value is given as an average of 3-fold cross validation with standard deviation between individual values for each fold.

	Frequency bands included	Batch size	# of hid-den units	Best classifier	Validation score
	Θ, Hα	5	3	Multi Perceptron Layer	0.750 (std: 0.165)
	∆, Hα	5 or 1	2	Random Forest	0.800 (std: 0.041)
	∆, Lα, Hα	5	2	Multi Perceptron Layer	0.750 (std: 0.065)
	∆, Θ, Lα	5	2	Multi Perceptron Layer	0.600 (std: 0.267)
	∆, Θ, Lα, Hα	5	5	Random Forest	0.650 (std: 0.138)
	∆, Θ, Lα, Hα, Lβ	5	4	Random Forest	0.700 (std: 0.135)
	∆, Θ, Lα, Hα, Lβ, Hβ	1	5	Random Forest	0.750 (std: 0.102)

We also demonstrate that, using the proposed method, it is possible with sufficient confidence to classify EEG records for each participant (best validation score: 0.800, see Table V). We believe that further improvement of the method can provide sufficiently reliable means for use even for user identification based on the collected EEG data.


Nazim Fatès
email: nazim.fates@loria.fr

Biswanath Sethi

Sukanta Das

On the reversibility of ECAs with fully asynchronous updating: the recurrence point of view (draft paper)

The reversibility of classical cellular automata is now a well-studied topic but what is reversibility when the evolution of the system is stochastic? In this context, we study a particular form of reversibility: the possibility of returning infinitely often to the initial condition after a random number of time steps. This corresponds to the recurrence property of the system.

We analyse this property for the 256 elementary cellular automata with a finite size and a fully asynchronous updating, that is, we update only one cell, randomly chosen, at each time step. We show that there are 46 recurrent rules which almost surely come back to their initial condition. We analyse the structure of the communication graph of the system and find that the number of the communication classes may have different scaling laws, depending on the active transitions of the rules (those for which the state of the cell is modified when an update occurs).

Introduction

The question of reversibility is an important topic in the study of discrete dynamical systems because it is intimately linked to the question of time. In the context of cellular automata, a great amount of efforts have been deployed to understand how the local rule which defines the interaction between cells may generate a global behaviour which is reversible [START_REF] Morita | Reversible computing and cellular automata -a survey[END_REF]. However, the majority of the results concern deterministic cellular automata where reversibility has been interpreted as "backward determinism" [AP72, [START_REF] Richardson | Tessellations with local transformations[END_REF][START_REF] Toffoli | Computation and construction universality of reversible cellular automata[END_REF], that is, the fact the each global state of the system, the configuration, has only a unique predecessor. In the specific case of finite cellular automata, this implies that the system will necessarily return to its initial state after some finite number of iterations.

The question can then be raised to know what is the equivalent of reversibility when the evolution of the cellular automaton is submitted to a form of noise. The non-determinism of the evolution implies that each configuration may have several successors. Consequently, the one-to-one correspondence between a configuration and its image can no longer be taken as a criterion of reversibility. Instead, the criterion which consists in looking if the system comes back to its initial condition may lead to interesting classifications of rules, depending on the conditions under which this "eternal return" phenomenon happens.

In this contribution, we consider the case of elementary cellular automata with fully asynchronous update. In other words, we study finite-size systems, with a local rule that involves nearest-neighbour interaction and binary states, and at each time step, we update one cell chosen randomly and uniformly in the set of cells.

So far, only a few results on the reversibility of asynchronous cellular automata have been obtained. In the case where one is allowed to choose freely the order of update of the cells, it was shown that there are some rules which allow one to return to the initial condition and some that do not always allow this [START_REF] Das | Synthesis of reversible asynchronous cellular automata for pattern generation with specific hamming distance[END_REF][START_REF] Sarkar | Reversibility in asynchronous cellular automata[END_REF]. With a different perspective, Worsch and Wacker examined how to construct an "inverse" rule, in the sense that its transition graph would be the "inverse" of the transition graph of the original inverse [START_REF] Wacker | On completeness and decidability of phase space invertible asynchronous cellular automata[END_REF]. The design of asynchronous circuits with reversible gates is also a current important topic of research (see e.g. Ref. [PLC + 13, LPC + 16] and references therein).

In a previous work, we analysed the set of elementary cellular automata with the aim of knowing in which case a cellular automaton may not return to its initial condition [START_REF] Sethi | Reversibility of elementary cellular automata under fully asynchronous update[END_REF]. This resulted in the identification of the irreversible rules in the sense that such rules would not always return to their initial condition: out of the 256 rules, there are 132 strongly irreversible rules, which possess at least configuration that is not reachable from itself, and 78 rules for which the system would return to its initial configuration a finite number of times and then stay away from it for ever. It was conjectured that the rest of the rules would be recurrent in the sense that that would always return to their initial condition an infinite number of times. Our goal is to formally prove this recurrence property for the set of the remaining 46 candidate rules.

Notations

basis

We consider one-dimensional binary cellular automata with periodic boundary conditions. The cells are arranged in a ring and the set of indices that represent each cell is denoted by L = Z/nZ, where n > 1 is the number of cells. We exclude the case n = 1 to avoid having a cell connected only to itself. At each time step t ∈ N, a cell is assigned a state in {0, 1}. The collection of all states at given time is called a configuration and the set of configurations is denoted by E n = {0, 1} L . We will often assimilate a configuration to a circular word and we use |x| w to denote the number of occurrences of the sub-word w ∈ Q * in the circular word formed by the states of x. A 1-region (resp. 0-region) corresponds to a maximal set of contiguous cells with the state 1 (resp. 0). Note that for a configuration different from 0 and 1, the number of 0-regions and 1-regions are equal.

An elementary cellular automaton (ECA) is a function f : {0, 1} 3 → {0, 1}. A rule f describes the local transformations (x, y, z) → f (x, y, z) with (x, y, z) ∈ {0, 1} 3 ; we call each quadruplet x, y, z, f (x, y, z) , a transition. We say that a transition is active if it changes the state of a cell (f (x, y, z) = y); it is passive otherwise. The decimal code associated to f is the number:

W (f ) = f (0, 0, 0) • 2 0 + f (0, 0, 1) • 2 1 + • • • + f (0, 0, 1) • 2 7 .
Since its introduction by Wolfram, this code has become a standard way to identify the ECAs. We will also use another code, called the transition code, to identify an ECA. This code consists of a concatenation of letters in {A, . . . , H}, where each letter denotes an active transition of the rule. The association of a letter to a transition follows the mapping of Table 1.

In the following, a rule will be denoted both by its decimal code and its transition code. For instance, majority rule is 232:DE (only 101 and 010 are active; this corresponds to RMT 5 and 2 in the notation used by Das et al. [START_REF] Sarkar | Reversibility in asynchronous cellular automata[END_REF]). The transition code of special rule that has no active transition is 204:I, where I stands for identity.

Given a local rule f and a ring size n, we define the global transition rule with fully asynchronous updating as F : E n × L → E n which maps a configuration x ∈ E n and a cell to update c ∈ L to the configuration y such that:

∀i ∈ L, y i = f (x i-1 , x i , x i+1 ) if i = c, x i otherwise.
Starting from an initial condition x, the evolution of the system can thus be described by the Markov chain (X t ) t∈N , where the random variables X t are obtained recursively with X 0 = x (with prob. 1) and X t+1 = F (X t , U t ) where (U t ) t∈N is a series of uniformly and identically distributed random variables on L.

Defining recurrence

Recall that our aim is to know, given a cellular automaton f and a ring size n, if this system always returns to its initial condition. Since the updating is stochastic (one cell is chosen randomly and uniformly at each time step), one may believe at first sight that advanced tools from probability theory are needed to answer this question. In fact, because we are dealing with finite and memoryless systems, we only need some simple notions form the Markov chain theory. We now introduce these tools.

Definition 1. For two arbitrary configurations x, y ∈ E n , we say that y is a successor of x if there exists an update u ∈ L such that F (x, u) = y. For x ∈ E n , we denote by succ(x) = {F (x, u), u ∈ L} the set of successors of x and succ(X) = ∪ x∈X succ(x) the set of successors of a set X ⊆ E n .

By iterating the successor relation, we obtain a new relation, "being reachable". Definition 2. For two arbitrary configurations x, y ∈ E n , we say that y is reachable from x if y can be obtained from x by following a sequence of successor relations, that is, ∃(x i ), i ∈ {0, . . . , k} such that x 0 = x, x k = y and x i+1 ∈ succ(x i ) for i ∈ {0, . . . , k -1}.

We denote by x y that y is reachable from x and the sequence of successors that leads from x to y is denoted by x = x 0 , x 1 , . . . , x k = y . We call this sequence a chain of successors, or simply a chain.

We can now define the recurrence property: Definition 3. A configuration x ∈ E n is recurrent if for every configuration y that is reachable from x, x is also reachable from y. A configuration that is not recurrent is transient.

Intuitively, a transient configuration is such that there exists a particular sequence of updates which brings the system to a configuration from which it will never be possible to return back to the initial configuration.

More generally, if y is reachable from x and x is reachable from y, we say that x and y communicate. By convention, all states communicate with themselves. Clearly, the relationship "communicate" is an equivalence relation; this relation partitions the set of configurations into communication classes. Consequently, two major behaviours exist: for the transient configurations, the system remains a finite time in its communication class, then "escapes" this class and never returns back to it. In contrast, when the system starts from a recurrent configuration, it remains in the communication class of this configuration for ever.

Definition 4. A rule is recurrent for size n if each configuration of E n is recurrent; otherwise, it is irreversible (for this size). A rule is recurrent if it is recurrent for all the sizes n > 1.

The set of the 256 ECA can be partitioned into recurrent and irreversible rules according to the following theorem:

Theorem 1. The rules with an asterisk are directly recurrent (see below).

Direct recurrence

Definition 5. We say that a rule has the property of direct recurrence if for any two configurations x, y ∈ E n , if y is a successor of x, then x is also a successor of y. Formally: ∀x ∈ E n , x ∈ succ(succ(x)).

Proposition 1. A rule that is directly recurrent is recurrent.

Proof. Consider a directly recurrent rule and assume that x and y are two configurations such that x y with the chain of successors x 0 = x, . . . , x k = y . The direct recurrence property gives x i ∈ succ(x i+1 ) for all i ∈ {0, . . . , k -1}, which implies that y x with the chain x k = y, . . . , x 0 = x . Proposition 2. A rule is directly recurrent if and only if it contains 0 or 2 active transitions in each of the four sets of transitions: {A, E}, {B, F}, {C, G}, {D, H}.

Proof. The proof is clear by looking at Table 1 (p. 3).

If the rule has no active transition, that is, if we have the identity rule, then the rule is recurrent because every point is a fixed point. Let us now assume for example that we go from x to y by applying transition A (000 → 1). Then, on the configuration y, on the cell where the transition was applied, it is now possible to apply transition E (010 → 0) to return to the previous configuration x. Symmetrically, if E was applied first, then it is possible to apply A to reverse the effect of the transition. The same arguments apply for the other couples of transitions {B, F}, {C, G}, and {D, H}. Reciprocally, if a rule is directly recurrent, one of these cases apply.

We thus determine the minimal rules which are directly recurrent as: 51:ABCDEFGH, 54:BCDFGH, 57:ACDEGH, 60:CDGH, 105:ADEH, 108:DH, 150:BCFG, 156:CG, 204:I.

A lemma for determining the communication classes

Lemma 1. For a given rule f and a given ring size n, let (E i ) i∈{1,...,k} be a partition of the configuration space E n and let (r i ) be a sequence of configurations such that ∀i ∈ {1, . . . , k}, r i ∈ E i . If for each i ∈ {1, . . . , k}:

(a) Any configuration of E i communicates with r i , (b) E i is closed under the application of F (that is, succ(E i ) ⊆ E i ),
then f is recurrent and its communication classes are formed by the sets (E i ).

To verify the lemma, one simply needs to observe that for each i ∈ {1, . . . , k}, the set E i is recurrent and maximal. Consider x ∈ E i , if we take a configuration y ∈ E i , then x and y communicate through r i and if we take a configuration y out of E i than y is not reachable from x, otherwise this would contradict the second hypothesis which stipulates that the set E i are closed.

Since communication classes are essentially equivalence classes, we call the configurations (r i ) the representatives of the partition (E i ). In our proofs, we will often show that two different representatives do not communicate, which is equivalent to the second property of the lemma.

Analysis of the recurrent rules

Let us now analyse the recurrent ECAs and their communication classes. We start with the rules that have a large number of communication classes and finish with the rules where the communication class is the whole set of configurations E n .

Rules with an exponential number of communication classes

Clearly, the identity rule 204:I is recurrent and because each configuration is a fixed point and thus constitutes its own communication class.

Proposition 3. Rules 156:CG is recurrent. The number of its communication classes C CG scales exponentially with the ring size n. More precisely, we have

C CG (n) ∼ φ n , where φ = (1 + √ 5)/2 is the golden ratio.
Proof. This rule is directly recurrent.

To count the number of communication classes, let us define the set of configurations which do not contain the pattern 11 as R = {x ∈ E n ; ∃i ∈ L; (x i , x i+1 ) = (1, 1)}. We want to show that the elements of R and the fixed point 1 constitute a set of representatives. Let us take a configuration x / ∈ R ∪ {1}. This configuration necessarily contains the pattern 11 and, since it is different from 1, it contains the pattern 110. The transition G (110 → 0) is active; x thus communicates with a configuration x that has a 1-region which has been decreased by 1 (e.g. x = 0011101 and x = 0011001). By transitivity ECA 108:DH ECA 156:CG Figure 1: Space-time diagrams of two ECA rules with an exponential number of communication classes. Blue and white squares represents cells in state 0 and 1, respectively. Time goes from bottom to top. Each row represents the state of the system after n random updates of cells. We have here n = 40 and the evolution is shown during 30 • n time steps. This convention is kept.

x communicates with x ∈ R which is such that all the 1-regions of x have been reduced to a size 1 (e.g. x = 0010001).

Moreover, if we take two different configurations in R, they do not communicate. Indeed, notice that the transitions C and G do not modify the positions of the 01-frontiers (see Fig. 1). The set R ∪ {1} is thus a set of representatives and we have C CG (n) = card R + 1. Now, it is well-known that the number f n of binary words of size n that do not contain the pattern 11 follows the Fibonacci sequence, defined by: f n = f n-1 + f n-2 with f 1 = 2 and f 2 = 3. This relationship can be proven with a recurrence on n by counting the number of words which end with a 0 or a 1, respectively. In order to obtain the size of R, we need to take into account the circularity of the configuration and thus to suppress the words of length n ≥ 3 which start with a 1 and end with a 1. The number δ n of such words also verifies δ n = δ n-1 + δ n-2 with δ 3 =1 (101) and δ 4 = 1 (1001). With the convention δ 2 = 0,

δ 3 = 1, for n ≥ 2, the number of classes of communication of rule CG reads C CG (n) = (f n -δ n ) + 1, which implies that it obeys the recursive equation C CG (n) = C CG (n -1) + C CG (n -2) -1 with C(1) = 2 and C(2) = 4.
As a result, we obtain that C CG (n) ∼ φ n , where φ = (1 + √ 5)/2 is the golden ratio1 . Proposition 4. Rules 108:DH is recurrent. The number of its communication classes scales exponentially with the ring size n.

Proof. This rule is directly recurrent. Let us evaluate the number of its communication classes. Let R be the set of configurations which do not have an isolated 0: We want to show that R constitutes a set of representatives.

Since transition D (101 → 1) is active, it is clear that any configuration x which does not belong to R communicates with a configuration x with one less isolated 0, and, by transitivity, with a configuration which has no isolated 0, that is with an element of R. Moreover, notice that for any configuration, a 0-region whose length is strictly greater than 1 is stable. Consequently, two different configurations of R cannot communicate because they necessarily have one 0-region for which they differ (otherwise they would be equal). The set R is thus a set of representatives.

Let us now estimate the size of R as a function of the ring size n. As for rule CG (see above), we can first evaluate the number g n of (non-circular) words which do not contain the pattern 101. By looking at the recurrence relations that govern the number of (legal) words which end with the patterns 00,01,10, and 11, one can prove that (g n ) verifies the recurrence equation 2 : g n = g n-1 + g n-2 + g n-4 with g(0) = 1, g(1) = 2, g(2) = 4 and g(3) = 7. One can then count the number δ n of words which contain an isolated 0which is exclusively created by the circularity of the configuration. This corresponds to the words of the form 01w1 and 1w 10 where w and w are binary words of length n -3 that do not contain an isolated 0. By applying the same techniques as above, one can prove that δ n obeys the same recurrence equation as g n , but with different initial conditions. More precisely, we have:

δ 1 = 0, δ 2 = δ 3 = δ 4 = 2 and δ n = 2g n-4 for n > 4.
As we have

C DH (n) = card R = g n -δ n = g n -2.g n-4
, we obtain the recurrence equation: C(1) = 10. It should be noted that the case n = 1, which was excluded from our study, does not correspond to a case where the rule is recurrent but it is here given only to initiate the recurrence 3 . To solve this linear recurrence equation, we need to find the roots of the polynomial P (x) = x 4 -x 3 -x 2 -1 = (x + 1)(x 3 -2x 2 + x -1). Using a program, we found that the four roots of P are -1, λ ≈ 1.75488, and two complex conjugate solutions whose module is strictly smaller than 1. As a result, we obtain that C DH (n) ∼ C.λ n where C is a constant which is determined to be equal to 1 (see Appendix p. 20).

C DH (n) = C DH (n -1) + C DH (n -2) + C DH (n -4) with C(1) = 2, C(2) = 2, C(1) = 5,

Rules with a linear number of communication classes

Proposition 5. Rules 134:BFG, 142:BG, and 150:BCFG are recurrent. Their communication classes correspond to the fixed points and to the configurations with an equal number of 1-regions.

The number of their communication classes is:

CC(n) = 3 + n/2 if n ∈ 2Z, 2 + n/2 otherwise.
Proof. Let us first focus on rule BG. This rule has two fixed points: 0 and 1 and for n even, it has two additional fixed points: (01) n 2 and (10)

n 2 .
2 This corresponds to the sequence A005251 of the OEIS (with different initial conditions since their sequence (a(n)) n∈N is given by a(n) = g n-3 ).

3 This corresponds to the sequence A259967 of the OEIS.

8

ECA 134:BFG ECA 142:BG Let us now consider a configuration that is not a fixed point. For k ∈ {1, . . . , n/2 } and k = n/2, let E k be the set of configurations with k 1-regions, that is, E k = {x ∈ E, |x| 01 = k}. We show that s k = (10) k 0 n-k is a representative of E k . Let us take a configuration x ∈ E k and examine if it communicates with s k . To show that x s k , we define the following intermediary configurations (see Tab. 2):

• x is obtained from x by reducing the 1-regions to a size 1 from their right part.

• x is obtained from x by "stacking" the 1-regions.

We now show that: x x x s k . We advise the readers to follow the text below with an eye on Table 2.

Step (a): To see that x x , we only need to remark that the transition G (110 → 0) allows one to decrease the length of the 1-regions when this length is greater or equal than 2.

Step (b): To show that x x , let us first observe the following property. If a configuration x contains the pattern 001, then the 1 can "move one cell to the left". Formally, if we have z = w 1 001w 2 , with w 1 , w 2 two words of {0, 1} such that |w 1 | + |w 2 | = n -3, then z z = w 1 100w 2 . Indeed, as the transition B (001 → 1) and G are active, we have z w 1 011w 2 z .

By iterating this argument, we have that if x has the form x = 0 i (10) j 10 k w with i, j ≥ 0 and k > 1 then x 0 i (10) j 10 k w. In words, this means that there are chains of successors which allow one to move a 1 to its leftmost position. Again, by considering each 1-region (from left to right) and moving the 1's to the left, we can build a chain such that x x .

Step (c): By considering the possibility of moving the 1's to the left as described above, we can show that for any configuration x that is not a fixed point, x σ(x) where σ denotes the left shift. By repeating the shift operation as many times as needed, we can show that there exists a chain of successors such that x s k . In the same way, it can be shown that s k y, where y is any configuration that is not a fixed point. We use the same arguments as above but follow a different path. First, we have that s k y , where y is the configuration obtained from s k with a given number of shifts such that the rightmost 1 of y coincides with the rightmost 1 of s k (see Tab. 2). Let y be the configuration where all the ones are aligned with the rightmost 1's of y. Clearly, y y . Finally, the relationship y y can be obtained with the application of the transition B to make the 1-regions grow to the appropriate size.

As the application of the rule BG does not modify the number of regions, the sets E k are closed under the application of this rule, and, naturally, so are the fixed points. Consequently, by applying Lemma 1, the rule BG is recurrent.

The same conclusion holds for rules BFG and BCFG. Indeed, the presence of the active transitions C and F gives the chains of successors more possibilities of "movement" in the configuration space but does not change the communication classes. We now turn our attention to what is probably the most interesting rule of our study. The number of its communication classes depends on the divisibility of n by 4 and in the general case the sizes of the classes of communication strongly vary depending on how much "freedom of movement" the system has. Proposition 6. Rule 105:ADEH is recurrent. The number of its communication classes CC depends on the ring size n according to:

CC(n) =      2 if n ∈ 2Z + 1 n/2 + 3 if n ∈ 4Z n/2 otherwise.
Proof. To understand the structure of the communication graph of rule ADEH, let us introduce the transformation T of the configuration from the states 0 and 1 to the binary alphabet Q = {., X} where X and . respectively denote the presence of absence of a change of state between a cell land its right neighbour. Formally, we associate to a configuration x ∈ Q L a word x = T (x) ∈ Q L such that ∀i ∈ L, x (i) = . if x i = x i+1 and x (i) = X otherwise. Note that T is not injective since a configuration x and its invert x, that is, the configuration where the state of all cells has been inverted, have the same image. The function T is not either surjective because, since the number of 0-regions and 1-regions is equal on rings, only images with an even number of X's can be obtained (see Tab. 3). The two homogeneous fixed points 0 and 1 are mapped to the all-. configuration (. L ). Now, remark that any application of the rule ADEH amounts to transforming a pair .. into a pair XX or a pair XX into a pair .. in the image of T . Moreover, remark that a class of communication in the image process corresponds to one or two communication classes of the original rule ADEH. Indeed, given a communication class of the image, there are two possibilities: either any two pre-images of this set communicate or they do not communicate and form two distinct classes (see examples below). These properties of the images allow us to determine the classes of communication of the original rule more easily. An example of evolution of the rule and its image is shown on Fig. 3 and Tab. 3. Moreover, an example of the communication graph is shown on Fig. 4 for n = 6.

(a) Let us start with the case where n is odd. Let us consider a configuration x / ∈ {0, 1}, its image T (x), and let us examine the regions of this image. Since n is odd, it can be remarked that there exists at least one region of . or X in T (x) that has an even size, otherwise we would have an even number of regions of odd size and n would be even. This region can be made "to disappear" by changing the state of all its pairs. In other words, if for example x = w 1 .X 2k .w 2 and y = w 1 . 2k+2 w 2 , where k > 0 and w 1 and w 2 are arbitrary words, then x and y communicate. By iterating this process, we obtain a homogeneous configuration, that is, an all-. or all-X configuration, but this latter possibility is excluded because the images necessarily have an even number of X's. Consequently, x communicates with any of the two homogeneous configurations 0 or 1 (the pre-images of all-.). Now, consider the number of 00 pairs: by looking at the effect of the transitions A, D, E and H, it can be seen that its parity is a conserved quantity, that is, if y ∈ succ(x), then |x| 00 mod 2 = |y| 00 mod 2. This implies that the two configurations 0 and 1 do not communicate and the space of configurations can be partitioned into two communication classes, each one containing the configuration with an even or an odd number of 00 pairs, respectively.

(b) We can now consider the case where n is even and write n = 2k. First let us examine what happens when n is not a multiple of four, that is, when k = n/2 = 2j +1 is odd. We claim that there are k classes of communications, which representative images can be taken as the configurations (.) n , (.) n-4i (X.) 2i and (.) n-4i (.X) 2i , with i ∈ {1, . . . , j}. In words, we can always reach an image with an even number of X's which alternate with .'s.

Indeed, by reasoning on the images and by looking at the 1-regions, we can see that the transformation of pairs XX → .. makes an arbitrary configuration x communicate with a configuration x whose image T (x ) is either all-., or where all the X's are isolated. Then, if T (x ) contains two 1's which are separated by an even number of cells, then x communicates with x such that the two 1's of T (x ) are replaced by two .'s in T (x ). By iterating this argument, we find that x communicates with y, which is such that either all the X's of T (y) are on even cells or they are all on odd cells. Last, we can see that y communicates with the representative where the 1's have the same parity as in T (y).

We now examine whether two representatives communicate. Let us consider the following quantities: e and o, the number of even and odd cells with state X in the image. For a configuration x ∈ E n , we define: e = card {i ∈ L; i ∈ 2Z and T (x) i = X} and o = card {i ∈ L; i ∈ 2Z + 1 and T (x) i = X}. It is then clear that q = eo is a conserved quantity, that is, if y ∈ succ(x), then q(T (y)) = q(T (x)). Indeed, as the X's appear or disappear by pairs, the two functions e and o can only increase or decrease by one for each active update. Since all the representative have different values of q, they can not communicate (see Fig. 4).

The last detail we need to settle is whether each representative image corresponds to one or two communication classes. Indeed, remember that each image corresponds to two pre-images, which are conjugate. We claim that in the case where k = n/2 is odd, each representative corresponds to a single communication class as its two pre-images communicate. The representative all-. corresponds the two pre-images 0 and 1, which both communicate with (01) k , and thus communicate together. If we take s = (.) n-4i (X.) 2i with i ∈ {1, . . . , j} as a representative, its two pre-images are r = 0 n-4i (0011) i and its conjugate r = 1 n-4i (1100) i . To show that r and r communicate, one may for example observe that r communicates with r = 1 n-4i (0011) i (see Fig. 4). Then, by noting that r and r have the same number of 1's, one may construct a chain of communication from r to r by turning pairs of 0's into pairs of 1's (with two updates) to "shift" the pairs of 0's to their right position.

(c) We can now examine our last case: when n is a multiple of four. We write n = 2k = 4j. In fact, all the arguments stated above still apply, with an exception which concerns the two representative images (01) k and (10) k . These two images have four pre-images, the four fixed points of the rule: (0011) j , (0110) j , (1001) j and (1100) j . As a fixed point is its own class of communication, there are two more communication classes than in the previous case. There are thus 1 + 2(j -1) + 4 = n/2 + 3 communication classes, which respectively correspond to the following representatives: 0, 0 n-2i (01) 2i , 0 n-2i (10) 2i with i ∈ {1, . . . , j -1}, and the four fixed points.

Rules with two communication classes

In the following, in many cases we will prove that 1 is a representative configuration. It is thus useful to know in which case 1 will be reachable and which are the configurations that can be reached from 1. Definition 6. For a given rule, we say that there is a density-decreasing path from x to y if ∃(x 0 , . . . , x k ) such that x 0 = x, . . . , x k = y and:

|x i+1 | 1 < |x i | 1 for all i ∈ {0, . . . , k -1}.
We say that 1 is an Olympic point of E if there is a density-decreasing path from 1 to any configuration of E \ 1. The symbols X and . represent the image by a transformation which displays the differences between a cell and its right neighbour (see text for details). Selfloops are not shown, links are bidirectional.

We say that 0 is an Avernal point of E if there is a density-decreasing path from any configuration of E \ 0 to 0.

It should be noted that the two definitions above are not exactly symmetric: in the first case we start from the Olympic point to reach other configurations, in the second case we start from a given configuration and reach the Avernal Proof. The proof is straightforward. Let us assume that the conditions (a) and (b) are verified. Let us consider a configuration x / ∈ {0, 1}. Then, if x does not contain a pair 00, that is, if all the 0's are isolated, then, x has a predecessor y where a 0 has been replaced by a 1 and such that transition H can be applied in this cell. Formally: ∃i ∈ L, y ∈ succ(x); ∀j = i, x j = y j and x i = 0 and (y i-1 , y i , y i+1 ) = (1, 1, 1). If x contains a pair 00, as it is different from 0, it contains a triplet 001 and a triplet 100. In the same way as we did before, if transition F (resp. G) is active we can construct a configuration y where a 0 has been replaced by a 1 and such that transition F (resp. G) can be applied in this cell. Formally: ∃i ∈ L, y ∈ succ(x); x i = 0, (y i-1 , y i , y i+1 ) = (0, 1, 1) or (y i-1 , y i , y i+1 ) = (1, 1, 0), and ∀j = i, x j = y j .

As a consequence, every non-uniform configuration has a predecessor which has a higher number of 1's, which implies that this configuration can be reached from 1.

In case E is active, it can easily be seen that 0 is reachable from the configurations with only one 1, which are themselves reachable from 1. The uniform configuration 1 is thus an Olympic point for E n . It can also be verified that the number of 1's can always be decreased, which implies that 0 is an Avernal point. Proof. According to Prop. 7, as transitions G and H are active, we know that 1 is an Olympic point for E * n (any configuration of this set is reachable from 1, see Tab. 4). Now, let us take x ∈ E * n and show x 1. It can be observed that if x = 1, for all these rules, as either B and D are active or C and D are active, it is always possible to increase the number of 1's. Indeed, if x contains the pattern 00, then it also necessarily contains the pattern 001 and 100, and transition B or C can be applied. If x does not contains the pattern 00, this means that all the 0's are isolated, that is, in the form 101, and transition D can be applied. The configuration 1 is thus reachable from x (see Tab. 4). Consequently, 1 is a representative configuration of E * n and, given that 0 is a fixed point that can not be reached from E * n , we have succ(E * n ) ⊆ E * n , which allows us to use Lemma 1 to conclude that the four rules are recurrent.

Rules with a single communication class

Proposition 9. For a ring size n = 3, the rules 33:ADEFGH and 41:ADEGH are recurrent. Their communication class is the whole configuration space.

The case of these two rules is interesting: its proof of recurrence is less direct than the previous rules. In particular, the rules are not recurrent for n = 3. Indeed, if one starts from 011, the configuration 001 is reachable, but it is not possible to go back to the initial condition from this configuration as transitions B and C are inactive. The rule is however recurrent for n > 3. rule). It would be interesting to relate this complexity of analysis to some formal properties of the communication graph. The fact that its has fixed points of length 4 but no fixed points of smaller size may give some hints on why this rule has such a peculiar behaviour.

We can note that determining the recurrence properties of the system seems somehow an easier task than, for example, distinguishing classes of asynchronous rules as to whether they converge rapidly or not to a fixed point [START_REF] Fatès | A guided tour of asynchronous cellular automata[END_REF]. The situation is similar with synchronous systems: only recently was achieved the classification of the 256 elementary rules according to their dynamical properties [START_REF] Schüle | A full computation-relevant topological dynamics classification of elementary cellular automata[END_REF]. One may thus legitimately ask to which extent the results on the recurrence obtained above can extend to a large class of rules, radius-2 rules or two-dimensional rules for instance.

Another path or research is to estimate the average time for returning to the initial configuration. For the "flip" rule (ECA 51), one may verify easily that the symmetry properties of the rule imply that the stationary distribution, that is, the distribution one approaches as time goes to infinity, is uniform over the set of configurations. As a result, the average return time of a n-sized system is 2 n . This result can be generalised to the directly reversible rules: the average return time from an initial configuration is equal to the size of the communication class of this configuration. For the other rules, where a dissymmetry between the states 0 and 1 may for example create a non-uniformity between the configurations of a communication class, a more in-depth analysis would be required.

Finally, we ask what definitions of the recurrence property can be given for an infinite-size system. In this context, one may either use interacting particle systems with continuous time or α-asynchronous systems, where the cells have a probability to be updated at each time step (see Ref. [START_REF] Fatès | A guided tour of asynchronous cellular automata[END_REF] for more details). One may then think about defining recurrence as the ability to come back infinitely often arbitrarily close to the initial condition. However, if we take the asynchronous shift as an example, even though this rule has a tendency to make the configurations more and more uniform, it is always possible that the information that is far away influences the central zone to make it "close" to the initial configuration. In the infinite context, contrarily to what we have seen for the finite case studied here, it seems that the property to erase information locally would not be sufficient to make a rule non-recurrent. To sum up, the notion of reversibility in the context of stochastic cellular automata is a rich problem and its exploration is only at its beginning.

Figure 2 :

 2 Figure 2: Space-time diagrams of two ECA rules with a linear number of communication classes.

  Figure 3: Space-time diagrams of ECA 105:ADEH (linear number of communication classes) and ECA46:BDGH (two communication classes).

Figure 4 :

 4 Figure 4: Communication graph for rule 105:ADEH and a ring of size n = 6.The symbols X and . represent the image by a transformation which displays the differences between a cell and its right neighbour (see text for details). Selfloops are not shown, links are bidirectional.

point 4 .

 4 Proposition 7. For a given rule, if (a) transition H is active and (b) transition F or G is active, then 1 is an Olympic point of E n \ 0. If conditions (a) and (b) are verified and transition E is active, then 1 is an Olympic point of the whole configuration space E n and 0 is an Avernal point of E n .

Proposition 8 .

 8 Rules 38:BDFGH, 46:BDGH and 54:BCDFGH, 60:CDGH are recurrent. Their communication classes are formed by the fixed point 0 and by the rest of the configuration space E * n = E n \ {0}.

Figure 5 :

 5 Figure 5: Space-time diagrams of two ECA rules for which all the configurations communicate.

Table 1 :

 1 Table showing the transitions and their labels.

	A	B	C	D
	000 001 100 101
	010 011 110 111
	E	F	G	H

Table 2 :

 2

	x 1100111010111	s 3 1010100000000
	hgcbfhgdedchh	ededecaaaaaab
	x 0000100010100	y 0000001010100
	aaabecabedeca	aaaaabededeca
	x 0000101010000	y 0100000010100
	aaabededecaaa	becaaaabedeca
	s 3 1010100000000	y 1100111010111
	ededecaaaaaab	hgcbfhgdedchh
	to a representative point	from a representative point

Rule BG: sequence of configurations showing that the communication between an arbitrary configuration and a representative configuration. Below each configuration, we indicate which transition can be applied in each cell (a for transition A, b for transition B, etc.).

Table 3 :

 3 Rule ADEH: two sequences of configurations showing that the communication between an arbitrary configuration and a representative configuration. On the right of each configuration is shown its image according to the transformation T (see text for details). Times runs from top to bottom.

	0110110001 X.XX.X..XX	0110110001 X.XX.X..XX
	0111110001 X....X..XX	0111110001 X....X..XX
	0111110101 X....XXXXX	0001010001 X..XXX..XX
	0111111101 X......XXX	0101010001 XXXXXX..XX
	0111111111 X........X	(...)	(...)
	1111111111 ..........	0000000000 ..........

Table 4 :

 4 Rule 46:BDGH: sequence of configurations showing that the communication between an arbitrary configuration and the Olympic point 1.

		x 1100010100011
		hgcabedecabfh
		x 1111010111011
		hhhgdedfhgdfh
		1 1111111111111
		hhhhhhhhhhhhh
	ECA 33:ADEFGH	ECA 43:ABDEGH

The first elements of the sequence C CG (n) are

2, 4, 5, 8, 12, 19, 30 and the sequence has index A001612 in the On-Line Encyclopedia of Integer Sequences (OEIS).

This has some logic: as expected one may go from the Olympus to any point, but the reverse is not granted, and one may also go from any point to the Avernus (the entry of the underworld in Roman mythology) but the other way round should be a more adventurous...

Proof. The proof of recurrence is decomposed into three steps: (a) 0 is an Avernal point; (b) 1 is reachable from 0; (c) 1 is an Olympic point.

Steps (a) and (c) are directly obtained from Prop. 7. To prove step (b) for the two rules, let us first assume that the ring size is even and write n = 2k. The sequence of configurations (x i ) i∈{1,...,k} with x i = (10) i 0 2k-i forms a chain of successors since the active transition A can always be applied in the cells in state 0 surrounded by two 0's (odd positions). We obtain 0 (01) k . Similarly, we have a chain of successors (y i ) with y i = (11) i (10)k -i because the transition D is active. We obtain (01) k 1 and thus 0 1. Now let us consider an odd ring size n = 2k + 1. Using the same operations as above, we can show that 0 z = 1 2k-1 00 (e.g., for z = 11111100 for n = 7). Now, if we want to reach 1, it is necessary to "go backwards". With the active transition G, we have z z = 11111000 and then, by applying transition A on the penultimate 0, a successor z z = 11111000 with two isolated 0's (e.g., 11111010). With the transition D, we then have z 1 which gives us 0 1. The three properties (a), (b) and (c) show the recurrence property of the two rules 33:ADEFGH and 41:ADEGH. Indeed for any two configurations x, y ∈ E n , we have x 0 1 y.

Proposition 10. The four rules: 35:ABDEFGH, 43:ABDEGH, 51:ABCDEFGH, and 57:ACDEGH are recurrent. Their communication class is the whole configuration space.

Proof. For these three rules, by Prop. 7 we know 1 is an Olympic point. Moreover, as either the three transitions A, B, and D or the three transitions A, C, and D are active, the number of 1's can always be increased from any configuration, which implies that 1 can be reached from any configuration. It thus a representative point of the whole configuration space E n and the rule is recurrent by Lemma 1.

Discussion

We have analysed the complete set of elementary cellular automata under fully asynchronous updating. We have asked if the behaviour of the system was recurrent and used the Markov chain theory to partition the set of rules between 46 recurrent rules (or 18 equivalence classes), and the rest of the rules, for which the return to the initial condition either never happens, or happens a finite number of times. The main tool to make this partition was to analyse the reachability relationships in the communication graph of the system. Moreover, the analysis of the communication classes showed various scaling relation of the number of communication classes as a function of the number of cells: we found the presence of an exponential number, a linear number and a constant number (two or one) or classes. The rule that displayed the greatest complexity in the analysis of the communication graph was rule 105:ADEH, which is rather surprising given that this rule can be written f (x, y, z) = 1 ⊕ x ⊕ y ⊕ z (affine

Appendix

We calculate here an equivalent the number of communication classes of rule 105:ADEH for a rinf of size n. For simplicity, let us denote this number by

This quantity verifies the following linear recurrence equation:

The values u 0 and u 1 are defined only for convenience and do not correspond to C DH (0) and C DH (1), which are excluded form our study. To solve this equation, we need to find the roots of the polynomial

There are two real roots, λ and -1, and two complex roots, which are conjugate, and which we denote by h and h . The general solution can thus be written :

where A, B, C, D are four constants which belong to the set of complex numbers.

Looking at the first terms, we have:

Now let us combine these equations with the evaluation of q = u 3 -2u 2 +u 1 -u 0 . We find that q = Q(λ)A+Q(h

In a second step, we evaluate

The evaluation of q -q allows one to get rid of B and we have :

It is now useful to write:

from which we obtain by identification λ + h + h = 2 and λh h = 1. Now, by using the relationships h + h = 2 -λ and h h = 1/λ in Eq. 2, we obtain: (2λ 2 -2λ + 3)A = (2λ 2 -2λ + 3), which leads to A = 1.

The algebraic or numerical estimation of the roots of Q gives : λ ≈ 1.75488 and h = 0.12256 + 0.74486i, whose module is strictly smaller than 1.

We thus obtain u n ∼ λ n , which was the desired result.