Nivine Abbas
email: nivine.abbas@telecom-paristech.fr

Thomas Bonald
email: thomas.bonald@telecom-paristech.fr

Berna Sayrac
email: berna.sayrac@orange.com

How Mobility Impacts the Performance of Inter-cell Coordination in Cellular Data Networks

Keywords: Cellular data networks, mobility, sector coordination, flow-level modeling, queuing theory

In this paper, we assess the performance of intercell coordination in the presence of mobility. This performance depends primarily on the resource allocation scheme. Indeed, a scheduling strategy which may seem efficient when users are static can lead to bad performance when users are mobile. Several scheduling policies are investigated. Their performance critically depends on their ability to predict users' mobility. The results are based on the analysis of flow-level traffic models.

I. INTRODUCTION

Inter-cell interference is one of the key challenges faced in mobile communication systems. It restricts the re-usability of the radio resource and limits spectral efficiency. Since the days of GSM, various techniques are used to cope with interference. Coordinated MultiPoint (CoMP) [START_REF]Requirements for Further Advancements for E-UTRA (LTE-Advanced)[END_REF], [START_REF] Irmer | Coordinated multipoint: Concepts, performance, and field trial results[END_REF], a main feature on the LTE-A roadmap, is mentioned as a promising approach to mitigate its effects through the coordination of multiple cells.

In this paper, we focus on joint processing (JP) schemes, which can be achieved through two main techniques: joint transmission (JT) and dynamic point selection (DPS). JT consists in transmitting data simultaneously from multiple cells to a single user; DPS consists in dynamically selecting the best transmitting cell in the coordination cluster and in muting the other cells. Under both JT and DPS, all cells of the coordination cluster are involved in the transmission and thus not available for other users. Thus there is a tradeoff between the performance of cell-edge users and the ability of the network to process all traffic.

It is well known that CoMP schemes are only advantageous for cell-edge users [START_REF]Further Discussion on Inter-Cell Interference Mitigation Through Limited Coordination[END_REF]. It has recently been shown that they can also be detrimental for cell-edge users at high load, due to the inefficient utilization of radio resources [START_REF] Khlass | Flow-level performance of intra-site coordination in cellular networks[END_REF]. Their efficiency in fact depends critically on the scheduling strategy. It is proposed in [START_REF]Discussion on CoMP-SU-MIMO[END_REF], [START_REF]Further discussion of frequency plan scheme on CoMP-SU-MIMO[END_REF] to allocate a dedicated frequency band to cell-edge users in order to perform CoMP operations. In [START_REF] Liu | A novel transmission scheme and scheduling algorithm for CoMP-SU-MIMO in LTE-A system[END_REF], the authors propose a joint proportional fairness scheduling algorithm that treats cell-center and cell-edge users equally without any frequency band partitioning. The scheduling algorithm may also depend on the clustering technique, from static clusters with centralized scheduling [START_REF] Brueck | Centralized Scheduling for Joint Transmission Coordinated Multi-Point in LTE-Advanced[END_REF] to dynamic clusters with distributed scheduling [START_REF] Zhou | Joint scheduling algorithms for LTE-A comp system[END_REF].

It is common to assume static users in the performance evaluation of CoMP schemes. We show in this paper that mobility has a strong impact on the performance of CoMP schemes. A scheduling strategy which may seem efficient when users are static can lead to bad performance when users are mobile. This is a rather surprising result since mobility is generally thought as improving throughput performance, see for instance [START_REF] Abbas | Opportunistic gains of mobility in cellular data networks[END_REF], [START_REF] Bonald | How mobility impacts the flow-level performance of wireless data systems[END_REF], [START_REF] Bonald | Flow-level performance and capacity of wireless networks with user mobility[END_REF], [START_REF] Borst | Capacity of wireless data networks with intra-and inter-cell mobility[END_REF], [START_REF] Borst | Mobility-driven scheduling in wireless networks[END_REF], [START_REF] Karray | User's mobility effect on the performance of wireless cellular networks serving elastic traffic[END_REF], [START_REF] Kim | Stability of flow-level scheduling with markovian time-varying channels[END_REF]. We study the impact of mobility on the performance of CoMP under several scheduling strategies. The results are based on the analysis of flow-level traffic models.

The rest of the paper is organized as follows. In the next section, we present the reference model under various scheduling schemes in the absence of mobility. We then analyze the impact of mobility and conclude the paper.

II. REFERENCE MODEL

We first present the reference model in the absence of mobility under various scheduling schemes: a static allocation, the proportional fair allocation and two prioritization strategies.

A. Cellular network

We consider a cluster composed of K coordinated cells, where each cell is modeled by two main zones: a noncoordination zone and a coordination zone, which is part of the global coordination zone of the cluster. The non-coordination zones k = 1, ..., K, as illustrated by Figure 1, are those where users are typically close to the cell center, experience good radio conditions and thus served only by their own serving base station. The coordination zone (k = 0) is that where users are close to the cell edge, suffer from bad radio conditions and receive signals from the K coordinated cells. This corresponds to the JT scheme. In fact, the model applies to the DPS scheme as well since the K coordinated cells are involved in the transmission in both cases; only the transmission rates achieved by the coordination depend on the JP scheme, either JT or DPS.

We model each zone k by a set of N k regions. In each region, radio conditions are supposed to be homogeneous and thus users are served at the same physical data rate on the downlink. We model each region by a queue with a specific service rate corresponding to the physical data rate in this region. Consequently, the considered cluster can be viewed as a set of K k=0 N k queues with K coupled processors. The simple case of two coordinated cells and one region in each zone is illustrated by Figure 1. Figure 2 shows an example of radio resource sharing between CoMP and non-CoMP users.

B. Traffic model

We consider elastic traffic only, corresponding to data transfers. Specifically, we assume that new data flows are generated in region i of zone k at the random times of a Poisson process of intensity λ ik . We denote by λ k = i λ ik the total flow arrival rate in zone k and by λ = k λ k the total flow arrival rate in the coordination cluster. The probability that a new data flow is generated by a user in zone k is p k = λ k /λ while the probability that a new flow in zone k is generated by a user in region i is p ik = λ ik /λ k . Physically, p k corresponds to the relative area of zone k compared to the area of the cluster, while p ik represents the relative area of region i of zone k compared to the area of zone k, assuming uniform traffic distribution in the cell.

We use a flow-level model where each data flow is viewed as a fluid of random volume to be transmitted. The volumes have an exponential distribution with mean σ (in bits). When region i of zone k is served, flows are completed at rate µ ik in the absence of fast fading, corresponding to the physical rate µ ik σ (in bit/s).

We denote by X ik (t) the number of active flows in region i of zone k at time t and by X k (t) = i X ik (t) the total number of active flows in zone k. The vector X(t) = (X ik (t)) ik defines a Markov process of dimension k N k . The load of region i in zone k is ρ ik = λ ik /µ ik . The total load of zone k is given by ρ k = i ρ ik = λ k /µ k , where µ k is the weighted harmonic mean service rate:

µ k = 1 i p ik /µ ik .
The resources of each cell are shared between CoMP users (in coordination zone k = 0) and non-CoMP users (in noncoordination zones k = 1, . . . , K), as illustrated by Figure 2. The actual service rate in zone k is modulated by φ k , the fraction of time spent by the scheduler on users in zone k. This depends on the system state x and on the scheduling policy.

For work-conserving policies, we have ∀k :

x k > 0, φ k (x) + φ 0 (x) = 1.
(1)

We focus on how to allocate resources between CoMP and non-CoMP users' data flows. In each zone, users are assumed to share the allocated radio resources equally, independently of their radio conditions. However, when the network operates under an opportunistic scheduler which exploits the fast fading, a rate gain [START_REF] Abbas | Opportunistic gains of mobility in cellular data networks[END_REF] should be taken into account in each region of each zone. Note that this gain depends on the number of flows in the considered region.

A necessary condition for the ergodicity of the Markov process X(t) is ρ < 1 in the absence of mobility, with

ρ = ρ 0 + max k=1,...,K ρ k , that is ρ < λ/µ with µ = p 0 µ 0 + max k=1,...,K p k µ k -1
.

The sufficient stability condition depends on the scheduling policy. We shall see in the following sections that mobility may increase or decrease the stability region, depending on the scheduling strategy.

C. Throughput metrics

We measure performance in terms of mean throughputs in the different zones. In any state x such that x ik > 0, each user in region i of zone k has throughput µ ik σφ k (x)/x k . Now the distribution seen by users in region i of zone k is the sizebiased distribution [START_REF] Bonald | Network Performance Analysis[END_REF]:

π ik (x) ∝ x ik π(x).
We denote by E ik the corresponding expectation. The mean throughput of users in region i of zone k is then given by

γ ik = E ik µ ik σφ k (X) X k = E(µ ik σφ k (X) X ik X k ) E(X ik ) . (2)
By the traffic conservation equation

λ ik = E µ ik φ k (X) X ik X k , (3)
we deduce

γ ik = λ ik σ E(X ik ) . (4)
This is the ratio of the traffic intensity in region i of zone k to the mean number of data flows in this region. Observe that, by Little's law, this is also the ratio of mean flow size σ to mean flow duration in region i of zone k.

By a similar argument, we obtain the mean throughput in zone k,

γ k = λ k σ E(X k ) , (5)
and the mean throughput in the cluster,

γ = λσ E( k X k ) . (6)
This is the harmonic mean of the mean throughputs γ k weighted by p k , for k = 0, 1, . . . , K.

D. Static allocation

A simple strategy consists in allocating a fixed fraction of resources to CoMP and non-CoMP users' data flows, as proposed by 3GPP in [START_REF]Discussion on CoMP-SU-MIMO[END_REF], [START_REF]Further discussion of frequency plan scheme on CoMP-SU-MIMO[END_REF]: a CoMP frequency sub-band is introduced to perform CoMP operation, so that all users selected to be served in CoMP mode use this pre-configured frequency sub-band. CoMP users' data flows are allocated some fixed fraction φ 0 of the resources of each coordinated cell while non-CoMP users' data flows in cell k are allocated the fraction φ k = 1-φ 0 of the considered cell radio resources.

The corresponding Markov process X(t) is reversible, with stationary distribution

π(x) = K k=0 (1 -ρ k /φ k ) φ x k k x k x 1 k , . . . , x N k k ρ x 1 k 1 k . . . ρ x N k k N k k , provided ρ 0 < φ 0 and ρ k < 1 -φ 0 for all k = 1, . . . , K.
Observe that this condition may be more restrictive than the natural stability condition ρ < 1.

From ( 5), the mean throughput in zone k is

γ k = φ k µ k σ(1 -ρ k /φ k ).
Observe that the mean throughput in zone k decreases linearly from the physical rate φ k µ k σ to 0 when the load of the corresponding zone increases from 0 to φ k . Indeed, the stability condition of zone k is ρ k < φ k .

E. Proportional fairness

Under proportional fairness [START_REF] Kelly | Rate control for communication networks: shadow prices, proportional fairness and stability[END_REF], resources are shared so as to maximize

K k=0 x k log φ k .
We get φ 0 (x) = x 0 K k=0 x k and φ k (x) = 1 -φ 0 (x) for all k and all x such that x k > 0.

We generalize the model presented in [START_REF] Khlass | Flow-level performance of intra-site coordination in cellular networks[END_REF] for K = 2 coordinated cells to any K ≥ 2. The stationary distribution of the Markov process X(t) is:

π(x) = π(0) x 0 + x 1 . . . + x K x 0 × K k=0 x k x 1 k , . . . , x N k k ρ x 1 k 1 k . . . ρ x N k k N k k ,
under the stability condition ρ < 1, with

π(0) = (1 -ρ 0 -ρ 1 )(1 -ρ 0 -ρ 2 ) . . . (1 -ρ 0 -ρ K ) (1 -ρ 0 ) K-1 .
By ( 5)-( 6), we deduce the mean throughput in the noncoordination zone k,

γ k = µ k σ(1 -ρ 0 -ρ k ),
and the mean throughput in the coordination zone,

γ 0 = µ 0 σ K k=1 1 1 -ρ 0 -ρ k - K -1 1 -ρ 0 -1
.

F. Priority to non-CoMP users

Under this policy, non-CoMP users in each cell are scheduled first and are allocated all the radio resources whenever active; CoMP users wait untill resources in all coordinated cells become available. The stability condition is:

ρ 0 < K k=1 (1 -ρ k ), (7)
which is more restrictive than the natural stability condition ρ < 1. This is due to the fact that some cells may be idle when serving non-CoMP users.

The mean throughput in non-coordination zone k is

γ k = µ k σ(1 -ρ k ).
There is no explicit expression for the mean throughput in the coordination zone.

G. Priority to CoMP users

Under this policy, CoMP users are scheduled first and are allocated all radio resources whenever active. Non-CoMP users are served only when there are no active CoMP users. The stability condition is the natural condition ρ < 1. Since CoMP users are not affected by non-CoMP users, we have

E(X 0 ) = ρ 0 1 -ρ 0 .
By ( 5), the mean throughput in the coordination zone is:

γ 0 = µ 0 σ(1 -ρ 0 ).
Applying known results of queuing theory [START_REF] Kleinrock | Queueing Systems[END_REF], the mean number of flows in non-coordination zone k is

E(X k ) = ρ k 1 -ρ 0 1 + ρ k + ρ 0 µ k /µ 0 1 -ρ k -ρ 0
We deduce the mean throughput in non-coordination zone k,

γ k = µ k σ (1 -ρ k -ρ 0 )(1 -ρ 0 ) (1 -ρ k -ρ 0 ) + ρ k + ρ 0 µ k /µ 0 .
Observe that the mean throughput of CoMP users is positive whenever ρ 0 < 1 while the mean throughput of non-CoMP users in zone k is positive whenever ρ 0 + ρ k < 1.

H. Numerical results

Figure 3 compares the throughput performance of the four policies for K = 2 coordinated cells. We assume that there is only one region in each zone with p 0 = p 1 = p 2 = 1/3, µ 1 = 2, µ 2 = 3/2 and σ = 1. This is the case of δP = 12dB: the threshold for the difference in the received power between the cooperating cells, and a 50% coordination gain. We take φ 0 = 1/2 for the static allocation (SA). Observe that proportional fairness (PF) outperforms the three other strategies. The priority to CoMP (Pri-C) strategy largely outperforms the priority to non-CoMP (Pri-NC) strategy at high load. In terms of stability, the Pri-NC strategy is the worst. This is due to the fact that prioritizing Non-CoMP users adversely affects the stability condition of CoMP users. The SA policy benefits from low complexity at the expense of degraded performance. The maximum load is approximately equal to 0.77 under the Pri-NC scheme and to 0.88 under the SA scheme, in accordance with the presented results. We observe that inter-zone mobility improves the mean throughput in the cluster for all strategies except for the Pri-C strategy. Under this strategy, mobility leads to a throughput degradation both for non-CoMP users and for CoMP users. This is due to the fact that non-CoMP users may then suffer from bad radio conditions (served with CoMP mode) leading to extra waste of resources. Based on this observation, the following conclusion may be drawn: it is not worth losing a resource to serve a CoMP user if that user is moving to the center and could benefit from better channel conditions. The strategy of prioritizing non-CoMP users proves to be the best strategy when users move. This can be explained by the fact that prioritizing Non-CoMP users gives the chance to celledge users to move and to be served in better radio conditions where coordination is not required. However, according to the numerical and analytical results the PF strategy remains a good compromise in a network where the users' mobility cannot be predicted. In accordance with the presented results, the maximum loads are approximately equal to 1.07, 1.23, 1.55 and 0.6 under the SA, PF, Pri-NC and Pri-C strategies, respectively. This shows that in terms of stability Pri-NC is the best strategy while Pri-C is the worst.

IV. CONCLUSION

This paper studies the impact of mobility on the performance of cell coordination, showing the critical impact of the scheduling strategy. By prioritizing CoMP users, mobility surprisingly incurs a throughput degradation. The proportional fairness strategy is a good compromise when the mobility of users is not known a priori. However, when mobility can be predicted accurately, it is better to prioritize non-CoMP users. Following this conclusion, we plan to consider more advanced scheduling strategies in the future. We shall also consider other CoMP techniques such as coordinated beamforming [START_REF]Coordinated Multi-point Operation for LTE Physical Layer Aspects[END_REF].

Fig. 1 .

 1 Fig. 1. A simple model with two coordinated cells.

Fig. 2 .

 2 Fig. 2. Resource allocation between CoMP and non-CoMP users.

Fig. 3 .

 3 Fig. 3. Throughput performance for two coordinated cells without mobility.

Fig. 5 .

 5 Fig.5. Throughput performance for two coordinated cells with mobility (solid line) and without mobility (dashed line).

III. IMPACT OF MOBILITY

In this section, we add mobility to the previous model. We use the same mobility model as in [START_REF] Abbas | Opportunistic gains of mobility in cellular data networks[END_REF]. Users move from a non-coordination zone of a cell to a non-coordination zone of another cell through the coordination zone, as shown in Figure 4 for two coordinated cells.

A. Mobility model 1) Inter-zone mobility: We assume that each user in region N k of non-coordination zone k moves to region N 0 of the coordination zone and vice versa after exponential durations, at respective rates ν k,0 and ν 0,k . In state x, the total mobility rate from non-coordination zone k to the coordination zone 0 and from coordination zone 0 to non-coordination zone k are equal to x N k k ν k,0 and x N00 ν 0,k , respectively.

2) Intra-zone mobility: We assume that in each zone k, each user in region i moves to region i -1 (for i > 1) and to region i + 1 (for i < N k ) after exponential durations, at respective rates ν (k) i,i+1 and ν (k) i,i-1 . In state x, the total mobility rate from region i to region i + 1 and from region i + 1 to region i

i+1,i , respectively. The probability that a user is in region i of zone k (for k > 0) then satisfies

, while the probability that a user is in region i of zone k = 0 satisfies

.

Note that

The probability that a user is in zone k is then given by:

It is worth noting that, for mobile users, the traffic conservation equation (3) no longer applies: the traffic arriving in region i of zone k is not equal in general to the traffic served in region i of zone k. Thus the mean throughput in region i of zone k is given by (2) but not by [START_REF] Khlass | Flow-level performance of intra-site coordination in cellular networks[END_REF]. Similarly, the mean throughput in zone k is given by

but not by [START_REF]Discussion on CoMP-SU-MIMO[END_REF]. Now the overall traffic conservation equation in the cluster still applies so that the mean throughput in the cluster is still given by [START_REF]Further discussion of frequency plan scheme on CoMP-SU-MIMO[END_REF]. It is still the arithmetic mean of the mean throughputs γ k weighted by the probability p k that an active user is in zone k, given by

for all k = 0, 1, . . . , K.

B. Stability condition

The stability condition follows from the limiting regime of infinite mobility where ν

i+1,i → ∞ for all i < N k and k ≥ 0, and ν k,0 , ν 0,k → ∞ for all k > 0 [START_REF] Bonald | Flow-level performance and capacity of wireless networks with user mobility[END_REF]. In this regime, the mean service rate in zone k becomes

The overall mean service rate is

under the PF scheme, Note that in the absence of inter-zone mobility, the stability condition is:

under the PF and Pri-C schemes, and

under the Pri-NC scheme, with

.

C. Throughput in light traffic

The performance in light traffic (that is, when ρ → 0) is the same for all policies except the static allocation: a user when alone in the system is always allocated all radio resources. For instance, we consider that there is only one region in each zone. The probability that a user in non-coordination zone k moves to the coordination zone before leaving the system is

Similarly, the probability that a user in the coordination zone moves to non-coordination zone k before leaving the system is α 0,k = ν 0,k µ 0 + ν 0 , where

The mean duration of a flow initiated in a non-coordination zone is

while the mean duration of a flow initiated in the coordination zone is given by:

We deduce the mean duration of a flow initiated in the cluster,

and the average throughput in light traffic γ = σ/T. Under the SA policy, we get the mean throughput similarly by replacing µ k by φ k µ k (∀k ≤ K).

D. Numerical results

Consider the same scenario of two coordinated cells as in the previous section, with ν 01 = ν 10 = ν 02 = ν 20 = 1. The results are obtained by the numerical evaluation of the stationary distribution of the Markov process X(t) and shown in Figure 5.


The Mayan Long Count Calendar Thomas Chanier * 1 1 Université de Bretagne Occidentale, 6 avenue Victor le Gorgeu, F-29285 Brest Cedex, France

The Mayan Codices, bark-paper books from the Late Postclassic period (1300 to 1521 CE) contain many astronomical tables correlated to ritual cycles, evidence of the achievement of Mayan nakedeye astronomy and mathematics in connection to religion. In this study, a calendar supernumber is calculated by computing the least common multiple of 8 canonical astronomical periods. The three major calendar cycles, the Calendar Round, the Kawil and the Long Count Calendar are shown to derive from this supernumber. The 360-day Tun, the 365-day civil year Haab' and the 3276-day Kawil-direction-color cycle are determined from the prime factorization of the 8 canonical astronomical input parameters. The 260-day religious year Tzolk'in and the Long Count Periods (the 360-day Tun, the 7200-day Katun and the 144000-day Baktun) result from arithmetical calculations on the calendar supernumber. My calculations explain certain aspect of the Mayan Calendar notably the existence of the Maya Epoch, a cycle corresponding to 5 Maya Eras of 13 Baktun. Modular arithmetic considerations on the calendar supernumber give the position of the Calendar Rounds at the Mayan origin of time, the Long Count Calendar date 0.0.0.0.0 4 Ahau 8 Cumku. Various long count numbers identified on Mayan Codices and monuments are explained by my approach. In particular, the results provide the meaning of the Xultun numbers, four enigmatic long count numbers deciphered in 2012 by Saturno et al. on the inner walls of a small masonry-vaulted structure in the extensive Mayan ruins of Xultun, Guatemala. 1 The results show a connection between the religious sites of Xultun and Chichén Itzá, Mexico. This is the first study linking unambiguously Mayan astronomy, mathematics and religion, providing a unified description of the Mayan Calendar.

Mayan priests-astronomers were known for their astronomical and mathematical proficiency, as exemplified in the Dresden codex, a bark-paper book of the XI or XII century CE containing many astronomical tables correlated to ritual cycles. However, due to the zealous role of the Inquisition during the XVI century CE Spanish conquest of Mexico, number of these Codices were burnt, leaving us with few information on Pre-Columbian Mayan culture. Thanks to the work of Mayan archeologists and epigraphists since the early XX century, the few codices left, along with numerous inscriptions on monuments, were deciphered, underlying the importance of the concept of time in Maya civilisation. This is reflected by the three major Mayan Calendars, reminiscent of the Mayan cyclical conception of time: the Calendar Round, the Long Count Calendar and the Kawil-direction-color cycle. The Calendar Round (CR) represents a day in a nonrepeating 18980-day cycle, a period of roughly 52 years, the combination of the 365-day civil year Haab' and the 260-day religious year Tzolk'in. The Tzolk'in comprises 13 months (numerated from 1 to 13) containing 20 named days (Imix, Ik, Akbal, Kan, Chicchan, Cimi, Manik, Lamat, Muluc, Oc, Chuen, Eb, Ben, Ix, Men, Cib, Caban, Etznab, Cauac, and Ahau). This forms a list of 260 ordered Tzolk'in dates from 1 Imix, 2 Ik, ... to 13 Ahau. 2 The Haab' comprises 18 named months (Pop, Uo, Zip, Zotz, Tzec, Xul, Yaxkin, Mol, Chen, Yax, Zac, Ceh, Mac, Kankin, Muan, Pax, Kayab, and Cumku) with 20 days each (Winal) plus 1 extra month (Uayeb) with 5 nameless days. This forms a list of 365 ordered Haab' dates from 0 Pop, 1 Pop, ... to 4 Uayeb. 3 The Tzolk'in and the Haab' coincide every 73 Tzolk'in or 52 Haab' or a Calendar Round such as the least common multiple (LCM) of 260 and 365: 73 × 260 = 52 × 365 = 18980 days. In the Calendar Round, a date is represented by αXβY with the religious month 1 ≤ α ≤ 13, X one of the 20 religious days, the civil day 0 ≤ β ≤ 19, and Y one of the 18 civil months, 0 ≤ β ≤ 4 for the Uayeb. To reckon time in a period longer than 52 years, the Maya used the Long Count Calendar (LCC), describing a date D in a 1872000-day Maya Era of 13 Baktun, a period of roughly 5125 years, counting the number of day elapsed from the Mayan origin of time, the mythical date of creation 0.0.0.0.0 4 Ahau 8 Cumku, corresponding to the Gregorian Calendar date 11 August 3114 BC according to the Goodman-Martinez-Thompson correlation. The XXI century saw the passage of the new Era on 21 December 2012, a date related to several apocaliptic or world renewal New Age theories, relayed by mass-media. Whereas Mayan mathematics are based on a vigesimal basis, the LCC is a 18-20-mixed radix basis representation of time: a date D is defined by a set of 5 numbers (C 4 .C 3 .C 2 .C 1 .C 0 ) such as D ≡ mod(D, 13 × 144000) = C 0 +C 1 ×20+C 2 ×360+C 3 ×7200+C 4 ×144000 where C 4 is the number of Baktun (144000 days) in the 13 Baktun Era (0 ≤ C 4 < 13), C 3 the number of Katun (7200 days) in the current Baktun (0 ≤ C 3 < 20), C 2 the number of Tun (360 days) in the current Katun (0 ≤ C 2 < 18), C 1 the number of Winal (20 days) in the current Tun (0 ≤ C 1 < 20) and C 0 the number of Kin (days) in the current Winal (0 ≤ C 0 < 20). The Kawil-direction-color cycle or 4-Kawil is a 3276-day cycle, the combination of the 4 directions-colors and the 819-day Kawil. 4 Table I gives the different calendar cycles with their prime factorizations.

Cycle P [day] Prime factorization Haab' 365 5 × 73 Tzolk'in 260 2 2 × 5 × 13 Winal 20 2 2 × 5 Tun 360 2 3 × 3 2 × 5 Katun 7200 2 5 × 3 2 × 5 Baktun 144000 2 7 × 3 2 × 5 Kawil 819 3 2 × 7 × 13 4-Kawil 3276 2 2 × 3 2 × 7 × 13 TABLE I.
Calendar cycles and their prime factorizations.

The origin of the Long Count Periods is unknown. A common assumption is the desire of the calendar keeper to maintain the Tun in close agreement with the tropical/solar year of approximately 365.24 days. 5 There is no consensus concerning the origin of the Tzolk'in, which has been associated with various astronomical cycles. 3 Tzolk'in correspond to Mars synodic period, 16 Tzolk'in equal 11 of Saturn synodic periods (+2 days), and 23 Tzolk'in are equivalent to 15 Jupiter synodic periods (-5 days). 6 It has been tentatively connected to the eclipse half-year (173.31 days) because 2 Tzolk'in are very close to 3 eclipse half-years. [START_REF] Aveni | Skywatchers: A Revised and Updated Version of Skywatchers of Ancient Mexico[END_REF] Finally, it has been noted that the Tzolk'in approximates the length of time Venus is visible as a morning or evening star. [START_REF] Bricker | Astronomy in the Maya Codices[END_REF] The Kawil cycle has been attributed to the observation of Jupiter and Saturn 9,10 because 19 (6) Kawil correspond to 39 (13) Jupiter (Saturn) synodic period. However, these interpretations fail to link the Tzolk'in and the Kawil to the Long Count Periods. The reason why the initial state of the Calendar Round at the LCC origin of time 0.0.0.0.0 is 4 Ahau 8 Cumku remains unexplained up to now.

In 2012, four LCC numbers, the Xultun numbers (Table II), have been discovered on the walls of a small painted room in the Mayan ruins of Xultun, dating from the early IX century CE. 1 These numbers have a potential astronomical meaning. Indeed, X 0 = LCM(260,360,365) is a whole multiple of the Tzolk'in, Haab', Tun, Venus and Mars synodic periods: 341640 = 1314 × 260 = 936 × 365 = 949 × 360 = 585 × 584 = 438 × 780, X 1 = 365 × 3276 is the commensuration of the Haab' and the 4-Kawil cycle. The greatest common divisor of the X i 's is 56940 = 3 CR corresponding to the commensuration of the Haab' and Mars synodic period of 780 days. However, the meaning of X 2 and X 3 , related to X 0 by the equation X 3 = X 2 + 2X 0 is unknown. Mayan astronomers-priests, known for their astronomical proficiency, may have observed with a naked eye the periodic movements of the five planets visible in the night sky: the moon, Mercury, Venus, Earth (solar year), Mars, Jupiter, and Saturn. Their respective canonical synodic periods are given in Table III. Evidences have been found in different Mayan Codices for Mercury, Venus, and Mars, but it is unclear whether they tracked the movements of Jupiter and Saturn. [START_REF] Bricker | Astronomy in the Maya Codices[END_REF] In particular, on page 24 of the Dresden codex is written the so-called Long Round number noted 9.9.16.0.0 or 1366560 days, a whole multiple of the Tzolk'in, the Haab', the Tun, Venus and Mars synodic periods, the Calendar Round and the Xultun number X 0 : LR = 1366560 = 5256 × 260 = 3744 × 365 = 3796 × 360 = 2340 × 584 = 1752 × 780 = 72 × 18980 = 4 × 341640. The relevant periods for the prediction of solar/lunar eclipses are the lunar semesters of 177 or 178 days (6 Moon synodic periods), which are the time intervals between subsequent eclipse warning stations present in the Eclipse Table in the Dresden Codex and the lunar tables inscribed on the Xultun walls. 1, [START_REF] Bricker | [END_REF] From their prime factorizations (Table III), we calculate the calendar supernumber N defined as the least common multiple of the P i 's:

N = 20757814426440 (1) = 2 2 × 3 2 × 7 × 13 × 19 × 29 × 59 × 73 × 89 = 365 × 3276 × 2 × 3 × 19 × 29 × 59 × 89 = LCM(360, 365, 3276) × 3 × 19 × 29 × 59 × 89
Equation 1 gives the calendar supernumber and its prime factorization. It is expressed as a function of the Haab' and the 4-Kawil. The Haab', canonical solar year, is such as the Haab' and the P i 's are relatively primes (exept Venus and Mars): the {LCM(P i ,365)/365, i = 1..8} = {116, 8, 1, 156, 399, 378, 177, 178} (Table III). The 4-Kawil and the Haab' are relatively primes: the LCM(365,3276) = 365 × 3276 and their largest common divisor is 1. The 4-Kawil has the following properties: the {LCM(P i ,3276)/3276, i = 1..8} = {29, 146, 365, 5, 19, 3, 59, 89}. That defines the 4-Kawil. The commensuration of the 4-Kawil and the Haab' X 1 = 365 × 3276 gives: {LCM(P i ,X 1 )/X 1 , i = 1..8} = {29, 2, 1, 1, 19, 3, 59, 89}. 360 is the integer closest to 365 such that the LCM(360,3276) = 32760 and the {LCM(P i ,32760)/32760, i = 1..8} = {29, 73, 73, 1, 19, 3, 59, 89}. The number 32760 or 4.11.0.0 has been derived from inscriptions on the Temple of the Cross in Palenque, Chiapas, Mexico. [START_REF] Lounsbury | A Rationale for the Initial Date of the Temple of the Cross at Palenque. Second Palenque Roundtable[END_REF] The Tun has the following properties: Y = LCM(360,365,3276) = 2391480 such as

{LCM(P i ,Y)/Y, i = 1..8} = {29, 1, 1, 1, 19, 3, 59, 89}.
The commensuration of the Haab', the 4-Kawil and the Baktun (400 Tun = 144000 days) gives rise to a calendar Grand Cycle GC = LCM(365,3276,144000) = 400 × LCM(360,365,3276) = 7 × 400 × X 0 = 956592000. The Euclidean of the calendar supernumber N by GC gives:

N = GC × Q + R (2) Q = 21699 R = 724618440 = 101 × 126 × 56940 = 126 × 3 i=0 X i .
If we pose A = 13 × 73 × 144000 = 400 × X 0 = 100 × LR such as GC = 7 × A, the Euclidean division of N by A gives:

N = A × Q + R (3) Q = 151898 R = 41338440 = 121 × 341640 = 121 × X 0
The only P i 's commensurate to A = 13 × 73 × 144000 are the Haab', Venus and Mars canonical periods (Table III) (9,13,20). This 2340-day cycle is present in the Dresden Codex on pages D30c-D33c and has been attributed to a Venus-Mercury almanac because 2340 = 20 × 117 = 5 × 585 is an integer multiple of Mercury and Venus mean synodic periods (+1 day). [START_REF] Bricker | Astronomy in the Maya Codices[END_REF] Another explanation may be of divination origin because 117 = 9 × 13. In Mesoamerican mythology, there are a set of 9 Gods called the Lords of the Night [START_REF] Aveni | Skywatchers: A Revised and Updated Version of Skywatchers of Ancient Mexico[END_REF][START_REF] Boone | Cycles of Time and Meaning in the Mexican Books of Fate[END_REF][START_REF] Lounsbury | Dictionary of Scientific Biography[END_REF][START_REF] Joyce | Mexican Archaeology: An Introduction to the Archaeology of the Mexican and Mayan Civilizations of Pre-Spanish America[END_REF] and a set of 13 Gods called the Lords of the Day. [START_REF] Joyce | Mexican Archaeology: An Introduction to the Archaeology of the Mexican and Mayan Civilizations of Pre-Spanish America[END_REF] Each day is linked with 1 of the 13 Lords of the Day and 1 of the 9 Lords of the Night in a repeating 117-day cycle. We can rewrite Equ. 2 and 3 as:

N -121 × X 0 = 151898 × A (4) N -126 × 3 i=0 X i = 151893 × A
The Long Count Periods appear in Equ. 4: 151898 = 338 + 360 + 7200 + 144000 and 151893 = 333 + 360 + 7200 + 144000. Adding the two equations in 4, we obtain:

5 × A = 5 × X 0 + 95 × 126 × 56940 (5) 5 × A = 5 × X 0 + LCM(X 1 + X 2 + X 3 , X 1 + 2X 2 + X 3 )
Since X 1 is known (the commensuration of the Haab' and the 4-Kawil), that defines X 2 and X 3 . The relationship X 3 = X 2 + 2X 0 may be a mnemotechnic tool to calculate Y = X 3 + 2X 2 = 105 × 56940. The four Xultun numbers provides a proof that Mayan astronomerspriests determined the canonical synodic periods of the five planets visible with a naked eye: Mercury, Venus, Mars, Jupiter and Saturn. The Mayan Calendar is then constituted by a calendar Grand Cycle GC = 7 × A such as 5 × A = 5 × 13 × 73 × 144000 = 12000 × 56940 = 73 × E with E = 5 × 13 × 144000 is the Maya Epoch corresponding to 5 Maya Eras of 13 Baktun. This corresponds to the interpretation of the Aztec Calendar stone (Fig. 1): the four squares surrounding the central deity represent the four previous suns or Eras and the center deity represents the Fifth Sun or the common Era. The Tzolk'in is defined by the LCM (13,20). It is such that the LCM(260,365) is 73 × 260 = 52 × 365 = 18980 days or a Calendar Round.

A question arises at this point to know how the Maya, as early as the IX century CE, were able to comput tedious arithmetical operations on such large numbers with up to 14 digits in decimal basis. Here is a possible method. They determined the prime factorizations of the canonical astronomical periods P i (Table III) and listed each primes p i with their maximal order of multiplicity α i . They calculated the calendar supernumber N (the LCM of the P i 's) by multiplying each p i 's α i time. The Euclidean division of N by GC = 7 × A = 7 × 400 × X 0 (Equ. 2) is equivalent to a simplification of N by 7 × 341640 and the Euclidian division of the product of the 5 left primes (3 × 19 × 29 × 59 × 89 = 8679903) by 400. The Euclidean division of N by 13 × 73 × 144000 = 400 × X 0 (Equ. 3) is equivalent to a simplification of N by 341640 and the Euclidian division of the product of the 6 left primes (3 × 7 × 19 × 29 × 59 × 89 = 60759321) by 400. It is to be noted that the prime factorization of the calendar supernumber only includes prime numbers < 100 which facilitates the operation (there are only 25 prime lower than 100).

Modular arithmetic considerations on the calendar supernumber allows to calculate the position of the Calendar Round at the Mayan mythical date of creation. For that purpose, we first create ordered lists of the Haab' and the Tzolk'in, assigning a number on the Haab' month and the Tzolk'in day. 2,3 For the Haab', the first day is 0 Pop (numbered 0) and the last day 4 Uayeb (numbered 364). For the Tzolk'in, the first day is 1 Imix (numbered 0) and the last day 13 Ahau (numbered 259). In this notation, the date of creation 4 Ahau 8 Cumku is equivalent to {160;349} and a date D in the Calendar Round can be written as D ≡ {mod(D + 160,260);mod(D + 349,365)}. The Calendrical Supernumber is such that: mod(N /13/73,260) = 160, mod(N /13/73,13) = 4, mod(N /13/73,20) = 0 and mod(N /13/73,73) = 49. The choice of 73 instead of 365 may be because the largest common divisor of 260 (2 2 × 5 × 13) and 365 (5 × 73) is 5. The date {160;49} corresponds to 4 Ahau 8 Zip, the day 0 (mod 18980), the beginning/completion of a Calendar Round. We now consider the Xultun number X 0 = LCM(260,365,360) = 18 × 18980 = 341640 which corresponds to the time interval between two days of the same Haab', Tzolk'in, and Tun date, for example between the date origin 0.0.0.0.0 4 Ahau 8 Cumku and the date 2.7.9.0.0 4 Ahau 8 Cumku corresponding also to a completion of a 13 Tun cycle, a period of 4680 days Finally, we discuss an important religious site in Mesoamerica, the pyramid of Kukulkan built sometime in the X century CE at Chichén Itzá (Figure 2) where various numbers in the architecture seem to be related to calendar considerations. The pyramid shape may be linked to the Long Count Calendar and the planet canonical cycles which draws a pyramid-like structure (Table IV). The pyramid is constituted of 9 platforms with 4 stairways of 91 steps each leading to the platform temple corresponding to the 3276-day cycle: 3276 = LCM(4,9,91) = 4 × 819, the coincidence of the 4 directions-colors with the Kawil. The Haab' is represented by the platform temple making the 365 th step with the 4 × 91 = 364 steps of the 4 stairways. Each side of the pyramid contains 52 panels corresponding to the Calendar Round: 52 × 365 = 73 × 260 = 18980. The dimensions of the pyramid may be of significance: the width of the top platform is 19.52 m (13 zapal), the height up to the top of the platform temple is 30 m (20 zapal) and the width of the pyramid base is 55.30 m (37 zapal), taking into account the Mayan zapal length measurement such that 1 zapal ≈ 1.5 m. [START_REF] O'brien | [END_REF] The pyramid height and the width of the top platform represents the Tzolk'in (13 × 20). The stairways divide the 9 platforms of each side of the pyramid into 18 segments which, combined with the pyramid height, represents the 18 Winal of a Tun (18 × 20). The width of the base 37 and the 4-Kawil 3276 are such that 37 × 3276 = LCM(148,3276) = 121212 which represents the coincidence of the 4-Kawil and the pentalunex used for solar/lunar eclipse prediction in the Dresden Codex. [START_REF] Bricker | [END_REF] If we include the pentalunex to calculate the calendar supernumber M = LCM(148, N ) = 37 × N = 36 × i p i , with {p i , i = 1..11} = {2, 3, 5, 7, 13, 19, 29, 37, 59, 73, 89}, the result of the Euclidean division of M by B = 7 × 37 × A = 37 × GC = 7 × 13 × 37 × 73 × 144000 can be expressed as a function of the Xultun numbers:

M = B × Q + R (6) Q = 21699 R = 26810882280 = 2 × 3 2 × 7 × 37 × 101 × 56940 = 7 × 18 × 37 × 3 i=0 X i .
During an equinox, the Sun casts a shadow (7 triangles of light and shadow) on the northern stairway representing a serpent snaking down the pyramid (Figure 2). The 7 triangles, the 18 segments, the width of the pyramid base (37 zapal) and the Xultun numbers (Table II) may be interpreted as a representation of Equ. 6. That gives a connection between the religious sites of Xultun and Chichén Itzá.

In conclusion, this study presents a unified description of the Mayan Calendar based on naked-eye astronomy. A calendar supernumber N is calculated by taking the least common multiple of 8 naked-eye astronomy canonical periods describing the planet synodic movements and the apparition of solar/lunar eclipse. This calendar supernumber defines the three major Mayan Calendar cycles: the 4-Kawil cycle, combination of the 4 directions-colors and the 819-day Kawil, the Calendar Round and the Long Count Calendar. The 360-day Tun, the 365-day Haab' and the 3276-day 4-Kawil are issued from the prime factorization of the 8 canonical astronomical input parameters. The 260-day Tzolk'in, the 360-day Tun, the 7200-day Katun and the 144000-day Baktun are obtained from arithmetical calculations on N . The correlation of the three major calendar cycles represents a calendar grand cycle GC. The two Euclidean divisions of N by GC and GC/7 = 13 × 73 × 144000 show the existence of a Maya Epoch constituted of 5 Maya Eras of 13 Baktun and lead to the Xultun numbers, dating from the early IX century CE and deciphered by Saturno et al. in 2012 inside a small room of the extensive Mayan ruins of Xultun, Guatemala.1 Modular arithmetic considerations on the calendar supernumber determine the position of the Calendar Round at the Mayan mythical date of creation 0.0.0.0.0 4 Ahau 8 Cumku, reflecting the Mayan cyclical concept of time. The results show a connection between the religious sites of Xultun, Guatemala and Chichén Itzá, Mexico. This study constitutes a breakthrough towards the understanding of Mayan ethnomathematics of divination, used to correlate ritual cycles with astronomical events in order to rythm political life and religious practices, embedding Maya civilization in the endless course of time.

FIG. 1 .

 1 FIG. 1. Mayan/Aztec Calendar stone representing the Five Suns, discovered in 1790 at El Zócalo, Mexico City, Mexico.

  , such as the LCM(584,365) = 37960 and LCM(780,365) = 56940 the length of the Venus and Mars Table in the Dresden codex. We have A = LCM(260,365,144000) = 100 × LR = 7200 × 18980 = 3600 × 37960 = 2400 × 56940 = 1000 × 234 × 584 = 60 × 2920 × 780. The commensuration of the Winal and 234 is the LCM(20,234) = 2340 = 9 × 260 = 20 × 117 = LCM

  . The completion of a Calendar Round corresponds to 18980 days elapsed such as mod(18980,4680) = 260. Starting the CR count at 4 Ahau 8 Zip, the next date in the ordered CR list such as mod(D,4680) = 0 is the date 4 Ahau 8 Cumku, 4680 days later. A date D is then expressed as {mod(D + 4680 + 160,260);mod(D + 4680 + 49,365)} = {mod(D + 160,260);mod(D + 349,365)}. The Calendar Round started therefore on a 4 Ahau 8 Zip, 4680 days earlier than the Long Count Calendar such that the starting date of the LCC is 0.0.0.0.0 4 Ahau 8 Cumku {160;349}.

FIG. 2 .

 2 FIG. 2. Pyramid of Kukulkan during an equinox. The pyramid is situated at Chichén Itzá, Yucatán, Mexico.

TABLE III .

 III Planet canonical cycles 1,12,13 and their prime factorizations.

	Xi LCC	D [day] Xi/56940
	X0 2.7.9.0.0 341640	6
	X1 8.6.1.9.0 1195740	21
	X2 12.5.3.3.0 1765140	31
	X3 17.0.1.3.0 2448420	43
	TABLE II. Xultun numbers Xi. 1 56940 = LCM(365,780) =
	2 2 × 3 × 5 × 13 × 73 is their largest common divisor and
	X3 = X2 + 2X0.
	Planet i Pi [day] Prime factorization
	Mercury 1 116 2 2 × 29
	Venus	2 584 2 3 × 73
	Earth	3 365 5 × 73
	Mars	4 780 2 2 × 3 × 5 × 13
	Jupiter 5 399 3 × 7 × 19
	Saturn 6 378 2 × 3 3 × 7
	Lunar	7 177 3 × 59
	senesters 8 178 2 × 89

TABLE IV .

 IV Divisibility of the calendar supernumber N by a polynomial expression of the type Pi = 13 × 73 × ( 6-i n=0 18 × 20 n ). Qi is the quotient of the Euclidean division of N by Pi.

	Name	i Ci [day]	N /Pi	Qi
	-Tun Katun Baktun Pictun Calabtun 5 57600000 N /13/73/ 1 0 18 N /13/73/ 6 0 Ci 18 1 360 N /13/73/ 5 0 Ci 360 2 7200 N /13/73/ 4 0 Ci 7215 3 144000 N /13/73/ 3 0 Ci 144304 4 2880000 N /13/73/ 2 0 Ci 2886428 0 Ci 57866020
	Kinchiltun 6 1152000000 N /13/73/C0	1215186420

W. A. Saturno, D. Stuart, A. F. Aveni and F. Rossi, Science 336, 714 (2012).

2 A. F. Aveni, in Skywatchers: A Revised and Updated Version of Skywatchers of Ancient Mexico (University of Texas Press, 2001), p.

143. 3 A. F. Aveni, in Skywatchers: A Revised and Updated Version of Skywatchers of Ancient Mexico (University of Texas Press, 2001),

p. 147. 4 H. Berlin and D. H. Kelley, The 819-day Count and direction-color Symbolism among the Classic Maya. Middle American Research Institute

Publication 26 (1961). 5 A. F. Aveni, in Skywatchers: A Revised and Updated Version of Skywatchers of Ancient Mexico (University of

Texas Press, 2001), p. 136. 6 J. S. Justeson, The ancient Maya Ethnoastronomy: An Overview of Hieroglyphic Sources, in World Archaeoas-


L Albera

H Becker

A Karfoul

R Gribonval

A Kachenoura

S Bensaid

L Senhadji

A Hernandez

I Merlet

Localization of spatially distributed brain sources after a tensor-based preprocessing of interictal epileptic EEG data

This paper addresses the localization of spatially distributed sources from interictal epileptic electroencephalographic data after a tensor-based preprocessing. Justifying the Canonical Polyadic (CP) model of the space-time-frequency and space-time-wave-vector tensors is not an easy task when two or more extended sources have to be localized. On the other hand, the occurrence of several amplitude modulated spikes originating from the same epileptic region can be used to build a space-time-spike tensor from the EEG data. While the CP model of this tensor appears more justified, the exact computation of its loading matrices can be limited by the presence of highly correlated sources or/and a strong background noise. An efficient extended source localization scheme after the tensorbased preprocessing has then to be set up. Different strategies are thus investigated and compared on realistic simulated data: the "disk algorithm" using a precomputed dictionary of circular patches, a standardized Tikhonov regularization and a fused LASSO scheme.

I. INTRODUCTION

ElectroEncephaloGraphy (EEG) is the noninvasive technique most indicated to capture transient events directly related to epilepsies (like interictal spikes, in particular). In the context of epilepsy surgery, the use of EEG requires to solve an electromagnetic inverse problem, that is to say to estimate the current generators underlying noisy EEG data. Theoretically, a specific electromagnetic field pattern may be generated by an infinite number of spatially distributed (extended) sources, leading to an ill-posed Brain Source Imaging (BSI) problem. In the last decades, many different regularization schemes [START_REF] Becker | A performance study of various brain source imaging approaches[END_REF] have been proposed such as the 4-ExSo-MUSIC technique [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the exso-music approach[END_REF]. However, their performance can be limited in the presence of several extended sources and background activity [START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF].

To overcome these problems, a tensor-based preprocessing can be applied. It consists in constructing a threeway array from the data and decomposing it using the Canonical Polyadic (CP) model [START_REF] Hitchcock | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF]. Such a strategy aims at separating the different sources of interest, leading to an easier estimation of their location and spatial extent. In the past, several authors have proposed to decompose a Space-Time-Frequency (STF) tensor, which can be obtained by computing a wavelet transform [START_REF] Miwakeichi | Decomposing EEG data into space-time-frequency components using parallel factor analysis[END_REF]- [START_REF] Deburchgraeve | Neonatal seizure localization using ParaFac decomposition[END_REF] over the time L. Albera, A. Karfoul, A. Kachenoura, L. Senhadji, A. Hernandez and I. Merlet are with the INSERM U1099, Rennes F-35000, France, and the Université de Rennes 1, LTSI, Rennes F-35000, France. L. Albera and R. Gribonval are with Inria, Centre Inria Rennes -Bretagne Atlantique, Rennes, F-35000, France. H. dimension of the measurements. In [START_REF] Becker | Multi-way space-time-wave-vector analysis for eeg source separation[END_REF], authors have considered a Space-Time-Wave-Vector (STWV) tensor, which is derived from a local spatial Fourier transform of the EEG measurements. An advantage of this approach compared to the STF analysis consists in its better robustness to correlated source activities, which is of particular interest in the context of interictal spikes when spreading of epileptic spikes is suspected between two regions.

The topic of extended source localization after the tensorbased preprocessing has been originally addressed in [START_REF] Becker | Multi-way space-time-wave-vector analysis for eeg source separation[END_REF] through the Binary Coefficient Matching Pursuit (BCMP) algorithm in the context of a spherical head model and of simplistic simulated EEG signals. More recently, a more efficient method, named Disk Algorithm (DA), has been designed [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF] to localize extended sources from the tensor decomposition. A detailed study of the tensor-based preprocessing methods has also been conducted [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF], including an analysis of their theoretical foundation, their computational complexity, and their performance for realistic simulated data in comparison to conventional source localization algorithms such as sLORETA [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF] and 4-ExSo-MUSIC.

However, justifying the CP model of the STF and STWV tensors is not an easy task when two or more extended sources have to be localized as deeply analyzed in [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF]. On the other hand, the occurrence of several amplitude modulated spikes originating from the same epileptic region can be used to build a space-time-spike tensor from the EEG data [START_REF] Deburchgraeve | Neonatal seizure localization using ParaFac decomposition[END_REF]. While the CP model of this tensor appears more justified, the exact computation of its loading matrices can be limited by the presence of highly correlated sources or/and a strong background noise. An efficient extended source localization scheme after the tensor-based preprocessing has then to be set up. Different strategies are thus investigated and compared on realistic simulated data in terms of estimation accuracy and CPU time: DA, sLORETA and a fused LASSO scheme named SVB-SCCD [START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF].

II. BSI PROBLEM FORMULATION

The electrical activity observed at the surface of the scalp originates mostly from pyramidal cells that are arranged in the gray matter perpendicularly to the cortex. These neuronal populations are generally modeled by a grid of dipole sources located at the surface of the cortex [9, Dale and Sereno's reference]. The electric potential data, which is recorded at the M electrodes of an EEG sensor array for T time samples, then constitutes the superposition of all dipole signals that are transmitted to the surface of the scalp and some instrumentation noise N (i) :

X = G (e) S (e) +G (b) S (b) +N (i) = G (e) S (e) +N (i,b) (1)
where the matrices G (e) ∈ R M ×P (e)

and G (b) ∈ R M ×P (b) model the propagation of the epileptic and background dipole activities, respectively, to the scalp electrodes. Note that the global propagation matrix G = [G (e) , G (b) ] ∈ R M ×P with P = P (e) + P (b) is called the lead field matrix. It depends on spatial parameters of the head, such as the geometry of the brain, skull, and scalp as well as their conductivities. For a given head model and source space, the lead field matrix can be computed by solving numerically Poisson's equation. The matrices S (e) ∈ R P (e) ×T and S (b) ∈ R P (b) ×T contain the epileptic and background dipole signals, respectively. The matrix N (i,b) ∈ R M ×T is the superposition of the background activity transmitted to the surface of the scalp and the instrumentation noise.

In the context of epilepsy, the brain sources of interest are spatially distributed and their spatial extent Ω (e) can be described as the union of contiguous areas Ω (e) r (r ∈ {1, . . . , R}) of the cortex (so-called patches) with highly correlated source activities [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the exso-music approach[END_REF]. As the activity is highly synchronized within one extended source, we can assume that all signals emitted by dipoles of the r-th (r ∈ {1, . . . , R}) extended source are approximately equal to a signal denoted by S r,: in the following. Consequently, we can approximate the data matrix by:

X ≈ R r=1 A :,r • S r,: + N (i,b) = AS + N (i,b) (2)
where • denotes the outer product operator [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. The spatial mixing vector A :,r of the epileptic patch Ω (e) r corresponds to the sum of the lead field vectors of all dipoles of Ω (e) r :

A :,r = pr∈Ω (e) r G :,pr = GZ :,r (3)
where the P -dimensional vector Z :,r of ones and zeros is

P (e)
r -sparse with P (e) r the number of dipoles of the epileptic patch Ω (e) r . The BSI problem consists then in recovering the coefficient matrix Z ∈ R P ×R (3) from the EEG data X (2) using the known lead field matrix G.

III. TENSOR-BASED BSI SOLUTIONS

The analysis of space-time EEG signals using the CP decomposition of a specific third order tensor has been originally proposed in [START_REF] Mocks | Topographic components model for event-related potentials and some biophysical considerations[END_REF]. More precisely, the author canonically decomposed a space-time-subject tensor built from the data by considering that there were observations of the same phenomenon on several subjects. If such a consideration is valid for the analysis of event-related potentials, it does not hold in the processing of epileptic EEG signals since the latter are patient specific. Nevertheless, as described in [START_REF] Deburchgraeve | Neonatal seizure localization using ParaFac decomposition[END_REF] and briefly recalled hereafter some information redundancies can be used in interictal epileptic EEG data allowing us to build a Space-Time-Spike (STS) tensor. Three different approaches, namely DA, sLORETA and SVB-SCCD, are then described hereafter in order to derive the location and the spatial extent of the epileptic regions from the CP decomposition of the considered tensor.

A. Tensor-based preprocessing

Spikes related to the same brain region appear several times in the interictal epileptic EEG data. The latter can then be segmented in time slots corresponding to different spikes. The J different space-time events of size (M × K) can thus be mutually synchronized with respect to the top of each spike and concatenated in a STS tensor X ∈ R M ×K×J enjoying the following low-rank CP model:

X = R r=1 A :,r • B :,r • C :,r + N (4)
where each A :,r • B :,r • C :,r is a rank-1 three-way array and where the third order array N represents the model residual. More particularly, A :,r and B :,r denote the spatial mixing vector and the spike-like pattern, respectively, of the r-th epileptic extended source. The weights of the length-K spike-like pattern B :,r characterizing the spikes of the r-th source are given by the J elements of C :,r . Sometimes we use the notation X (A, B, C) to refer to the third order tensor X with its loading matrices A =

[A :,1 , • • • , A :,R ] ∈ R M ×R , B = [B :,1 , • • • , B :,R ] ∈ R K×R and C = [C :,1 , • • • , C :,R ] ∈ R J×R .
The CP model ( 4) defines a rank-R approximation of tensor X . If R = rank(X ) then the model residual is null and the CP decomposition is exact. In the latter case, there is almost surely a finite number of decompositions of the tensor X into the three loading matrices A, B and C, up to scale and permutation indeterminacies, if R < (M KJ)/(M + K + J -2) [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]; tighter bounds also give sufficient conditions for uniqueness such as the condition based on Kruskal rank [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. Note that Alternating Least Squares (ALS) algorithm [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF] was used to compute the CP decomposition of X .

B. Extended source localization

The objective of extended source localization consists in recovering the coefficient matrix Z (3). This can be achieved from an estimate A = A :,1 , . . . , A :,R of the spatial mixing matrix A computed by means of the CP decomposition of X (4), which permits to identify the grid dipoles for each patch separately. In principle, any localization algorithm such as the fused LASSO approach SVB-SCCD [START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF] and sLORETA [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF] that act on a vector of spatial measurements can be employed to this end. Here, we briefly describe DA [START_REF] Becker | EEG extended source localization: tensor-based vs. conventional methods[END_REF], a surface grid based version of sLORETA and SVB-SCCD.

1) DA: The concept underlying DA consists in recovering the extended source from a number of small circular-shaped patches of grid dipoles, the disks. Each extended source is then obtained by the union of several disks. For each grid dipole, several disks composed of the 0 to D max -1 nearest dipoles and the current grid dipole as central point are determined. Assuming equal coefficients for all grid dipoles belonging to the extended source, the spatial mixing vector of each disk is computed by adding the lead field vectors of all grid dipoles belonging to the disk. The reconstructed spatial mixing vector of the j-th disk (j ∈ {1, . . . , D max P }) is given by G u j where the elements of the coefficient vector u j that are associated with dipoles belonging to the disk are equal to 1 whereas all other elements are 0. To determine which disks of the parameter space best describe the r-th patch, the spatial mixing vectors u j of all disks are then compared to the estimated spatial mixing vector A :,r of the r-th source by means of the squared normalized inner product:

ϕ 1 (u) = A T :,r Gu 2 A T :,r A :,r u T G T Gu -1 (5)
The r-th patch is then estimated as the union of a certain number of disks which are associated with the highest values of the metric ϕ 1 [START_REF] Miwakeichi | Decomposing EEG data into space-time-frequency components using parallel factor analysis[END_REF].

2) sLORETA: sLORETA as all minimum norm solutions is based on the assumption that among all possible solutions that explain the data, the one with minimal energy is the most likely. In this context, the idea of the conventional sLORETA algorithm consists in standardizing the minimum norm solution Z :,r by weighting its components by their estimated variances. Thanks to an inversion lemma, the minimum argument of the following objective function:

ϕ 2 (u) = || A :,r -Gu|| 2 + λ||u|| 2 (6)
where λ is a regularization parameter, can be written as

u = Q A :,r with Q = G T (GG T + λI M ) -1 ∈ R P ×M
the Tikhonov-regularized inverse matrix of G and where I M denotes the identity matrix of R M ×M . The idea of sLORETA consists in weighting solution u by the inverse of the diagonal matrix (diag{diag{F }})1/2 where the matrix F of the form:

F = G T (GG T + λI M ) -1 G (7)
is the covariance matrix of u. The p-th coefficient (p ∈ {1, . . . , P }) of Z :,r is then computed as:

Z p,r = (F p,p ) -1 2 Q p,: A :,r (8)
The resulting coefficient vector Z :,r contains components with continuous values. The r-th patch is estimated by identifying the components of Z :,r which exceed a certain threshold.

3) SVB-SCCD: The SVB-SCCD algorithm [START_REF] Becker | Fast, variation-based methods for the analysis of extended brain sources[END_REF] computes a sparse piece-wise constant spatial solution, by imposing sparsity on both the coefficient vector Z :,r and its variational map. The latter can be computed by applying a linear transform, characterized by the matrix V , to the coefficient vector, which is equivalent to computing the total variation on the discretized cortical surface. The components V ,p of V ∈ R L×P where L is the number of edges of the triangular grid, are given by:

V ,p =    1 if p = p ,1 -1 if p = p ,2 0 otherwise (9)
where p ,1 and p ,2 are the indices of the dipoles sharing the -th edge. This leads to the following cost function:

ϕ 3 (u) = 1 2 || A :,r -Gu|| 2 F + λ(||V u|| 1 + α||u|| 1 ) (10)
whose minimization giving the coefficient vector Z :,r is equivalent to the sparse TV (or fused LASSO) approach. The latter permits us to adjust the size of the reconstructed source regions by varying the new regularization parameter α. In addition, it prevents the estimated signal vector from being amplitude-biased, which is a problem that frequently arises using only the TV norm. Setting α = 1 leads to very focal source estimates, whereas small α only avoid the amplitude bias, but do not influence the size of the reconstructed source regions. The regularization parameter λ regulates the impact of the source priors. Note that the Alternating Direction Method of Multipliers (ADMM) [3, Boyd et al.'s reference] was used to minimize ϕ 3 [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF].

IV. SIMULATIONS

In this section, we compare the performance of STS-DA, STS-SCCD and STS-sLORETA (i.e. the three tensor-based BSI methods described in section III) on realistic simulated interictal epileptic EEG data.

A. General settings and performance criteria

For STS-SCCD, the regularization parameter λ was adjusted such that the reconstruction error lies within a confidence interval of 95 to 99% of the noise power and we used a fixed parameter α = 0.67 because we found that this leads to reasonable results for the considered scenario.

The EEG signals were generated for N = 91 electrodes using a realistic head model with three compartments that represent the brain, the skull, and the scalp. The source space consists of D = 19626 dipoles corresponding to the triangles of the cortical surface mesh with orientations perpendicular to the cortical surface. A BEM method 1 was used to compute the lead field matrix.

We considered three patches, located in the inferior frontal gyrus (patch InfFr), the occipital temporal junction (patch OccTe), and the mid temporal gyrus (patch MidTe). Each one is composed of 100 adjacent dipoles. The patches InfFr and OccTe were attributed different epileptic spike signals comprising 200 time samples (at 256 Hz sampling frequency) that were segmented from stereotactic EEG (SEEG) recordings of a patient suffering from epilepsy. Then we generated 40 different realizations for each of these two signals, one for each patch dipole, by introducing small variations in amplitude and delay. Assuming that the patch MidTe was activated due to a propagation of the epileptic activity of the patch OccTe, we used the same signals for the dipoles of these two patches, but introduced a delay of 8 ms due to the distance between the patches.

All source dipoles that do not belong to a patch were attributed Gaussian background activity with an amplitude that was adjusted to the amplitude of the SEEG signals

B. Computer results

The performance achieved with the three BSI tensor-based algorithms in terms of DLE and CPU time averaged over 40 realizations is summarized in table I show that the STS-SCCD algorithm clearly gives the best compromise between estimation accuracy and processing time for the considered multi-patch scenario. Recall that the latter involves two highly correlated patches and one decorrelated patch. Consequently, two columns of the loading matrix B are quasi-collinear and the ALS-like rank-R CP decomposition of X does not succeed in perfectly estimating the spatial mixing matrix A. In fact, one column of A corresponding to the spatial mixing vector of patch InfFr appears to be well estimated contrarily to the two other columns of A. In other words, the tensor-based preprocessing allows us to extract only one patch among the three patches. Note that it is a good issue provided that the two nonseparated sources can be also localized from the estimated mixing matrix A . The good behavior of STS-SCCD can also be seen in figure 1, where we illustrate an example of the BSI results obtained with the different methods. Contrarily to STS-SCCD, STS-sLORETA gives a blurred solution while STS-DA finds only two patches.

V. CONCLUSION

This paper addresses the localization of spatially distributed sources from interictal epileptic electroencephalographic data after the canonical polyadic decomposition of a space-time-spike tensor built from the data. Different strategies of extended source localization are investigated, namely the "disk algorithm" using a precomputed dictionary of circular patches, a standardized Tikhonov regularization and a fused LASSO scheme. The best compromise between estimation accuracy and processing time is offered by the fused LASSO approach. Future work will include i) the use of other tensor decompositions such as the block term decomposition [START_REF] Hunyadi | Block term decomposition for modelling epileptic seizures[END_REF], to deal with the case of extended sources and ii) the design of constrained tensor decomposition algorithms (i.e. the structure of the spatial mixing matrix is imposed during the tensor decomposition step).

Fig. 1 . 2 F ≈ 1 .

 121 Fig. 1. BSI results obtained with the three tested algorithms.

ASA, ANT, Enschede, Netherlands