Amelia Carolina Sparavigna

The alignment of Casamari Abbey along the moonrise azimuth on Easter 1203

Keywords:

Casamari Abbey is a Cistercian abbey in Italy, about 10 kilometres from Veroli. The building of the church that we see today started in 1203. Here we show that a link can exist between the direction of the moonrise on Easter 1203 and that of the nave of the church. To determine the moonrise azimuth we use software CalSKY. To simulate it on satellite images, we use The Photographer's Ephemeris.

Casamari Abbey is a Cistercian abbey in Italy, about 10 kilometres from Veroli. The abbey was built at the site of Cereatae, the birthplace of Caius Marius. A chronicle of the abbey, written in the 13th century, tells that the Abbey had its origin from a Benedictine monastery of 9th century [1]. In the 12th century, the abbey suffered a long period of decline. In this period, St Bernard of Clairvaux, the promoter of the Cistercian reform of monasticism, incorporated Casamari in his new order. He officially enlisted it in the Cistercian directory as the 29th foundation of Citeaux [1]. Under the Cistercians, the abbey and its church were completely rebuilt, between 1203 and 1217, in accordance with their own standards [1]. Today, Casamari is a well preserved example of the Burgundian early-Gothic architecture. Of this style, another example in Italy is the Abbey of Fossanova

[1-3].
It is known that the building of the Cistercian Casamari abbey started in 1203 [1]. Filippo Rondinini in his Monasterii Sanctae Mariae Et Sanctorum Johannis Et Pauli de Casaemario Brevis Historia is telling that "da un antichissimo codice membranaceo conservato dai monaci consta nell'anno del Signore 1203, epacta sesta, indizione corrente, il 6 maggio, nell'anno sesto del pontificato di Innocenzo III, viene posta dal reverendo abate del monastero Gerardo la prima pietra per la nuova basilica, benedetta e rimessa dal detto sommo pontefice" [START_REF]Monasterii sanctae Mariae et sanctorum Johannis et Pauli de Casaemario brevis historia studio & opera Philippi Rondinini Faventini digesta & illustrata[END_REF][START_REF] Farina | La grande opera del grande abate Giraldo. 10 Settember[END_REF]. And also that "Per procura di fratel Rajnerio da Ponza, di santa memoria, il papa Innocenzo III getta le fondamenta per la nuova basilica".

Here we will show that a link can exist between the year of foundation and the orientation of the abbey. In particular we investigate the following possibility, that the nave of the church had been aligned close to the moonrise azimuth of the first full moon after the spring equinox of year 1203. That is, we consider the possibility that the abbey was aligned along the moonrise azimuth on Easter 1203. The proposed discussion is therefore an example of an astronomical analysis, applied to find a possible link of the foundation of a church to the moonrise on Easter.

First, let us note that an orientation of the nave of the Abbey, that is of its man axis, according to the sunrise on the day of the foundation (6 May 1203, Julian Calendar) is impossible. We can use software CalSKY (see Table 1). We have that the sunrise azimuth, astronomical horizon, on 6 May 1203 had an azimuth equal to 64.4°. Table 1 It is true that the natural horizon is different from the astronomical horizon. Therefore, let us consider the altitude, on the astronomical horizon, that the sun (or the moon) must have to be seen rising above the natural landscape. Using elevation profiles of Google Earth, in the direction of the nave of Casamari Abbey, we find that the altitude of sun or moon must be 2,3 degrees. To appreciate the difference, let us use the Photographers Ephemeris software. This software is using a proleptic Gregorian Calendar, then 6 May 1203, Julian Calendar, is corresponding to 13 May 1203, proleptic Gregorian Calendar. The result of the simulation -screenshot -is given in the following Figure 1. As we can see from Figure 2, we have a remarkable alignment of the nave along the moonrise direction on 30 March 1203. We decided to use this date because it was import for the Church: it was the Easter of 1203, as we can easily find by using the "computus". It is a calculation that determines the date of Easter in the Julian Calendar [6,7]. Therefore the case seen in the Figure 2 seems representing sunrise and moonrise on Easter 1203. However, let us further investigate the moonrise azimuth, using CalSKY too.

We know that Easter is the Sunday after the first full moon which happens after the spring equinox. According to CalSKY, the spring equinox was on 12 or 13 March 1203 (the sunrise azimuth was 90.3° on 12 March and 89.8° on 13 March -astronomical horizon). The full moon after the spring equinox is given on 29 March 1203 (Julian Calendar), the day before than that we obtained by means of the computus. CalSKY tells that 29 March was Saturday, so 30 March was Sunday and therefore 30 March was Easter 1203.

Software CalSKY is also giving, day by day, the moonrise azimuth, here in the Table 2.

Table 2 Let us consider the Figure 2 once more. From satellite images, the direction of the nave of the abbey is 108° (within a degree). The moonrise direction according to the natural horizon has the same azimuth of the church, 108°. The moonrise azimuth on the astronomical horizon is 105°( within a degree). So we are in the same astronomical case of 30 March 1203, given in the Table 2.

On that day, the moon was rising on the natural horizon in the direction of the nave of the abbey.

Usually, it is told that churches are oriented according to the sunrise on the day of their foundation (see for instance [START_REF] Sparavigna | The solar orientation of the Gothic cathedrals of France[END_REF]), but this is not the case of Casamari Abbey. We have the date of the ceremony of the foundation, 6 May 1203, during which the first stone was laid. The orientation of the church was not given by the sunrise on that day, without any doubt. Moreover, it is reasonable that planning and orientation of buildings were prepared and decided before the ceremony. In the case of Casamari, an astronomic link to the moon on Easter 1203 is possible.

Some examples of orientation according to the moonrise are existing (see [START_REF] Carolina | La Basilica di Sant'Andrea a Vercelli ed il Lunistizio Settentrionale del[END_REF] and references therein). So the case of the Casamari Abbey could be a remarkable example of a church oriented according to the moonrise on Easter, as shown by the simulations obtained by means of CalSKY and the Photographer's Ephemeris. However, it is necessary to note that the alignment of the nave could be simply an occurrence in the absence of any reference to the moon. For this reason, further studies are necessary to investigate other Cistercians abbeys and their possible astronomical orientation, with respect to sun, moon and stars.

Figure 1 -

 1 Figure 1 -Many thanks to The Photographer's Ephemeris for the simulation and Google for the satellite images. Let me stress that the image is used for scientific purposes. The sunrise azimuth on the astronomical horizon is represented by the thick yellow line. The sunrise azimuth on the natural horizon (sun altitude 2,3°) is represented by the thin yellow line. It is evident that the nave of the church was not aligned along the sunrise azimuth.

Figure 1 -

 1 Figure 1 -Many thanks to The Photographer's Ephemeris for the simulation and Google for the satellite images. Let me stress that the image is used for scientific purposes. Moonrise azimuth (thick pale blue line) on astronomical horizon is 105° (within a degree). Moonrise azimuth (thin pale blue line) on natural horizon is 108° (within a degree). The simulation is given on 6 April 1203 (proleptic Gregorian Calendar), corresponding to 30 March 1203 (Julian Calendar).

Zenodo, 13 September 2019 -DOI 10.5281/zenodo.3407199


Guangshuo Chen
email: guangshuo.chen@inria.fr

Aline Carneiro Viana

Marco Fiore

Carlos Sarraute

Complete Trajectory Reconstruction from Sparse Mobile Phone Data

Keywords: human mobility, mobile phone dataset, data enrichment, seamless trajectory reconstruction, location predictability, trajectory unicity

Mobile phone data are a popular source of positioning information in many recent studies that have largely improved our understanding of human mobility. These data consist of time-stamped and geo-referenced communication events recorded by network operators, on a per-subscriber basis. They allow for unprecedented tracking of populations of millions of individuals over long periods that span months. Nevertheless, due to the uneven processes that govern mobile communications, the sampling of user locations provided by mobile phone data tends to be sparse and irregular in time, leading to substantial gaps in the resulting trajectory information. In this paper, we illustrate the severity of the problem through an empirical study of a large-scale Call Detail Records (CDR) dataset. We then propose Context-enhanced Trajectory Reconstruction, a new technique that hinges on tensor factorization as a core method to complete individual CDR-based trajectories. The proposed solution infers missing locations with a median displacement within two network cells from the actual position of the user, on an hourly basis and even when as little as 1% of her original mobility is known. Our approach lets us revisit seminal works in the light of complete mobility data, unveiling potential biases that incomplete trajectories obtained from legacy CDR induce on key results about human mobility laws, trajectory uniqueness, and movement predictability.

Introduction

Over the past ten years, the proliferation of mobile devices indirectly fuelled many original studies on human mobility, as mobile phone data collected by network operators have become an invaluable source of information about individual movement patterns at scale. Data extracted from mobile communications have supported many breakthroughs [START_REF] Blondel | A survey of results on mobile phone datasets analysis[END_REF][START_REF] Naboulsi | Large-scale mobile traffic analysis: A survey[END_REF], including seminal works demonstrating that individual movements are regular [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF] and predictable [START_REF] Song | Limits of predictability in human mobility[END_REF], as well as unique [START_REF] De Montjoye | Unique in the crowd: The privacy bounds of human mobility[END_REF]. Mobile phone data have allowed unveiling major patterns in human movements, such as home-work commuting [START_REF] Ahas | Modelling home and work locations of populations using passive mobile positioning data[END_REF], or frequent and recurrent motifs [7]. The novel understanding of individual movements enabled by network operator data has facilitated innovation in, e.g., transportation systems [START_REF] Ferreira | Visual exploration of big spatio-temporal urban data: A study of new york city taxi trips[END_REF], urban planning [START_REF] Zhang | comobile: Real-time human mobility modeling at urban scale using multi-view learning[END_REF], communication infrastructure deployment [START_REF] Zang | Mining call and mobility data to improve paging efficiency in cellular networks[END_REF][START_REF] Oliveira | From routine to network deployment for data offloading in metropolitan areas[END_REF], and epidemics control [START_REF] Frias-Martinez | An agent-based model of epidemic spread using human mobility and social network information[END_REF].

Call Detail Records (CDR) are the current de-facto standard for mobile phone data used in human mobility studies. As illustrated in Table 1, CDR are logs of events generated during mobile communications, such as voice calls, text messages, and possibly mobile data traffic sessions, which are collected for billing purposes [START_REF] Naboulsi | Large-scale mobile traffic analysis: A survey[END_REF].

Those events are associated with individual mobile devices and are time-stamped and geo-referenced, which makes CDR an obvious proxy for the trajectories of the mobile network subscribers. Interestingly, CDR encompass huge populations (e.g., millions of users) over large (e.g., citywide or nationwide) geographic regions, and cover long periods (e.g., months or years): they thus enable human mobility investigations at unprecedented scales.

A drawback of CDR-based trajectories (or simply trajectories hereafter) is that they are often largely incomplete. As a matter of fact, telecommunication events are punctual and provide information at specific time instants, which are also sparse and irregularly distributed in time. Thus, CDR offer a very partial view of the overall mobility of each user, with substantial continuous periods where location information is entirely absent [START_REF] Chen | Enriching sparse mobility information in call detail records[END_REF]. From a research standpoint, the reduced completeness of CDR yields severe consequences: (i) as many studies require exhaustive knowledge of personal trajectories, they focus on few highly active individuals that generate a large number of CDR events, and discard the rest -in several important works, this amounts to ignoring over 98% of the subscriber population [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF][START_REF] Song | Limits of predictability in human mobility[END_REF]; (ii) when the whole user base is considered, the dependability of the results is questioned by the fact that they are derived from partial movement information, and potential biases due to missing locations are difficult to assess [START_REF] Ranjan | Are call detail records biased for sampling human mobility[END_REF].

A solution to the problem of mobile phone data sparsity is trajectory reconstruction, which aims at completing individual movement information by recovering the unspecified positions of each user. However, as detailed in Section 2, the literature on trajectory reconstruction for mobile phone data is relatively thin, and mostly focuses on a partial reconstruction that still leaves (possibly large) gaps in our knowledge of personal mobility. In this paper, we present a novel and more effective approach to trajectory reconstruction, named Context-enhanced Trajectory Reconstruction (CTR), which allows completing all missing trajectory data with reasonable accuracy. CTR enables human mobility analyses over much larger populations, and without biases due to incomplete information. We apply CTR to the problem of reconstructing trajectories from CDR, since they represent the primary type of movement data employed for large-scale human mobility analysis [START_REF] Blondel | A survey of results on mobile phone datasets analysis[END_REF][START_REF] Naboulsi | Large-scale mobile traffic analysis: A survey[END_REF]. However, our method is general and can be applied to other classes of mobile phone data; specifically, by running CTR on data originally collected at a frequency of a few hundreds samples per day, one could aim at completing individual trajectories at high temporal resolutions, e.g., of minutes. The design, evaluation, and application of our solution yield the following contributions.

• First, we provide evidence of the severe sparsity that affects mobile phone data, by analyzing the CDR of 1.8 million users collected during three consecutive months. We quantify the phenomenon by means of relevant metrics, including the duration, sampling frequency, and completeness of individual trajectories in the data. In particular, our results show that legacy preprocessing techniques that discard users based on arbitrary completeness thresholds, in fact, ignore a vast user population with substantial and potentially serviceable mobility information. Details are in Section 3. • Second, we introduce our novel approach, CTR, which leverages well-known features of human mobility to customize tensor factorization in a way that befits our problem. This original methodology lets CTR transform sparse mobile phone data into seamless individual trajectories that span the full dataset duration, for all users. Details are in Sections 4 and 5.

• Third, we validate the proposed strategy with ground truth data. Comprehensive performance evaluation shows that CTR achieves full reconstruction of individual trajectories on an hourly basis, with a median displacement between 1 and 2 network cells that depends on the sparsity of the original CDR data. This effectively means that, in the reconstructed data, a user is typically placed in the correct cell or in one that is very close to it. Such a level of accuracy is acceptable for metropolitan-scale analyses (where the urban surface is typically covered by hundreds of cells) and is excellent for national-scale studies (as inter-city mobility is perfectly captured). Details are in Section 6. • Fourth, we demonstrate the importance of trajectory reconstruction in CDRdriven human mobility analysis. Specifically, we revisit three seminal studies [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF][START_REF] Song | Limits of predictability in human mobility[END_REF][START_REF] De Montjoye | Unique in the crowd: The privacy bounds of human mobility[END_REF] by using the complete mobility of 1.7 million users, instead of the incomplete trajectories of a small fraction of especially active users as in the original works. Our results show that key results in these studies may change, even quite dramatically, in presence of complete mobility, proving that trajectory reconstruction is an indispensable first step for analyses of human mobility that rely on mobile phone data. Details are in Section 7.

We then draw conclusions and comment on the perspectives of our work in Section 8.

Related work

A number of strategies adopted in mobile phone data analysis aim at mitigating the sparsity and irregularity of the individual trajectories extracted from CDR.

We classify current solutions depending on whether they are targetting (i) time discretization, (ii) user filtering, or (iii) trajectory reconstruction, and discuss them below.

Time discretization

A very common practice in mobile phone data analyses, time discretization reduces the temporal resolution of the data so as to make trajectories appear complete. Indeed, discretizing time allows associating a single location to each time step, implying that the location is valid across the whole interval covered by the step. We remark that time discretization is basically unavoidable in CDR, due to the limited granularity of mobile communication events. Interestingly, however, it is performed at quite different temporal resolutions by works in the literature, spanning from 15 minutes [START_REF] Ranjan | Are call detail records biased for sampling human mobility[END_REF] to 1 hour [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF][START_REF] Song | Limits of predictability in human mobility[END_REF], or even to 2 hours [START_REF] De Montjoye | Unique in the crowd: The privacy bounds of human mobility[END_REF]. Also related to data alteration in time, some studies focus on specific periods (e.g., work hours, commuting hours, nighttime hours), which allows them to treat the (time-discretized) data as complete even if just such target periods are covered [START_REF] Chen | Enriching sparse mobility information in call detail records[END_REF]. Ultimately, time discretization can be regarded as a data preprocessing step, which does not explicitly address the problem of trajectory incompleteness in CDR: indeed, it does not offer any guarantee in terms of the fraction of time intervals during which location information is available for a given user.

User filtering

It is customary for works on mobile phone data to filter the available user base before further processing. User filtering can be adopted to eliminate outliers (e.g., call centers [START_REF] Sarraute | A study of age and gender seen through mobile phone usage patterns in mexico[END_REF]), but it is more often considered to retain trajectories that provide sufficient mobility information. Indeed, when extensive knowledge of the user mobility is required by the nature of the investigation, the straightforward way to deal with incomplete trajectory data is that of filtering users so as to exclude from the analysis those with insufficient location information. This is typically performed by imposing thresholds on the amount of mobility data available in CDR, though there is hardly a standard to define them. As a reference, several important works opt for a choice of users who have a minimum average sampling frequency of 0.5 events per hour, and who are observed at two unique locations at least [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF][START_REF] Song | Limits of predictability in human mobility[END_REF][START_REF] Ranjan | Are call detail records biased for sampling human mobility[END_REF].

User filtering comes with a cost in terms of below-threshold users that are excluded from the analysis. Unfortunately, in practical cases involving CDR datasets, even the loose requirements set by the thresholds mentioned above risk to lead to severe reductions in the examined user population. For instance, two very well-known studies analyze just 1.67% [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF] and 0.45% [START_REF] Song | Limits of predictability in human mobility[END_REF] of the total subscribers available in their datasets upon filtering, leaving substantial richness in the original data untapped.

Trajectory reconstruction

The objective of trajectory reconstruction is to infer the positions of the users at times during which the original mobile phone data does not provide such information. By recovering the missing location data, a successful reconstruction makes user movements more complete, hence effectively obviates the problem of discarding vast portions of users caused by filters on mobility data.

The literature on trajectory reconstruction from mobile phone data, including CDR, is relatively thin, and mostly focuses on partial reconstruction. A fairly standard solution is extending the validity of instantaneous locations over longer time intervals, by assuming that users do not move for some periods (e.g.,, in the order of hours) before and after every recorded event in the CDR dataset [START_REF] Jo | Spatiotemporal correlations of handset-based service usages[END_REF]. Enhancements to this approach have also been proposed, by dynamically adapting the duration of location validity periods via supervised learning [START_REF] Chen | Enriching sparse mobility information in call detail records[END_REF]17]. A different strategy hinges on probabilistic reconstruction, implemented by modeling positions between two CDR events as random variables [START_REF] Ficek | Inter-call mobility model: A spatio-temporal refinement of call data records using a gaussian mixture model[END_REF]. However, all these techniques incur into non-negligible spatial errors, and yet, as they target partial reconstruction, cannot return full trajectories that cover the entirety of time [START_REF] Chen | Enriching sparse mobility information in call detail records[END_REF].

Spatiotemporal interpolation [START_REF] Hoteit | Estimating human trajectories and hotspots through mobile phone data[END_REF] is the only method proposed to date for complete trajectory reconstruction, i.e., a reconstruction that returns seamless location information for each user at the target time discretization. Multiple interpolation methods can be adopted to this end, and the approach has been proven effective with dense CDR entries featuring more than 1000 events per day. However, users with such a high level of communication activity are infrequent, and indeed only 0.07% of the total subscribers in the original CDR dataset are retained for analysis in [START_REF] Hoteit | Estimating human trajectories and hotspots through mobile phone data[END_REF]. When applied to typical CDR data, interpolation methods are ineffective and provide no clear improvement over a naive error-prone method where users are associated with the last location they are seen at [17].

The solution we propose in this paper also falls in the trajectory reconstruction category, and, like interpolation, aims at inferring complete individual trajectories. However, we seek an approach that effectively also operates on very sparse positioning data, and thus can be applied to entire CDR datasets.

Sparsity in CDR-based trajectory data

As a preliminary step to our study, we carry out a data-driven investigation of sparsity in individual trajectories extracted from mobile phone data. Our goal is quantifying the phenomenon, so as to demonstrate the significance of the problem and to motivate the need for a dedicated solution. To this end, we employ a CDR dataset collected from 1.8 million prepaid mobile network subscribers of a major Latin America mobile network operator. The data consists of 778 million voice call entries generated by customers during a consecutive 90-day period. The customers appear within a large geographic area of several thousands of km 2 covered by 4, 200 cell towers. As illustrated in Table 1, each CDR entry contains the caller's and the callee's hashed identifiers, the call duration, the event time, and the identifier of the cell tower which the caller's device was connected to when the call originated. The mapping of each cell tower identifier and its location in latitude/longitude is also provided by the data source. The results presented below refer to the CDR-based trajectories of all users in the dataset.

Trajectory completeness

We first quantify the amount of missing information in the individual trajectories. To this end, we assume a total observation time T and discretize it into time intervals of duration τ , which we will refer to as the temporal resolution. We define the completeness of a specific trajectory as the fraction of time intervals for which we have at least one location sample. As an example, let us assume that T = 7 days and the temporal resolution is τ = 1 hour: then, a trajectory with locations in 80 different hours has a completeness of 80/(7 × 24) = 0.476.

In practice, T corresponds to the duration of the dataset from which the trajectories are extracted, and τ to the temporal resolution required for the mobility analysis. As both values may vary based on the nature of the mobility study, we investigate their impact on completeness through comprehensive parametric analysis. Instead, our notion of completeness is independent of how positions are determined in the presence of multiple samples in the same time interval. Hence we do not need to account for that aspect at this stage.

Figure 1 portrays the cumulative distribution function (CDF) of the completeness of 1.8 million trajectories from our reference dataset when different combinations of T and τ are considered. Each plot refers to a different temporal resolution τ , while curves identify different observation periods T . We observe the following.

• Seamless trajectories with completeness equal to 1 are extremely rare in mobile phone data, under any combination of T and τ . Even with short observation periods (T = 7 days) and low resolution (τ = 2 hours), not a single complete trajectory is extracted from our reference dataset, despite the large user population. In fact, less than 1% of users would satisfy a minimum completeness requirement of 0.5, which allows half of the locations to be unknown; and, more than 20% of trajectories do not even meet the degree of completeness 0.1, i.e., miss 9 locations out of 10. The fact that such poor figures refer to the best case, i.e., the less demanding combination of T and τ , illustrates well the severity of the sparsity problem in CDR datasets. • The duration of the observation period has a marginal impact on trajectory completeness. More precisely, and as one would expect, slightly higher completeness is recorded for trajectories that span a shorter observation period when T ∈ [START_REF] Paul | Understanding traffic dynamics in cellular data networks[END_REF]90] days; however, and quite interestingly, completeness hardly varies when T ≤ 30 days. We argue that the result can be linked to the weekly periodicity of human activities (which entail a reduced difference for T ∈ [7, 30] days), plus seasonal effects such as summer holidays that only intervene over long timescales (T ≥ 30 days in our dataset). • The temporal resolution has a stronger effect on completeness. For instance, when T = 30 days, only 8% of the trajectories have completeness above 0.1 when τ = 15 minutes, whereas the same percentage grows to 26%, 53% and 75% when τ = 30, 60 and 120 minutes, respectively. To better understand the scaling of completeness with the temporal resolution, we show in Figure 2 the mean and median values of completeness, computed across all users, versus τ . We conclude that completeness grows almost linearly (Pearson correlation coefficients are at 0.981 and 0.983, for the mean and median, respectively) with τ in the range 15 to 120 minutes, under all observation periods T . The takeaway message is that one cannot easily escape the tradeoff that higher temporal resolution incurs into curbed completeness in CDR data. Overall, our results highlight that the vast majority of trajectories from a typical CDR dataset only capture a rather small portion of the complete user mobility. Based on Figure 2, we quantify the fraction of known positions between 5% and 20% on average, depending on the temporal resolution needed for the analysis. These results are well aligned with the outcome of user filtering in previous studies that leverage mobile phone data, where only a small subset of users displays sufficient completeness and is retained for analysis [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF][START_REF] Song | Limits of predictability in human mobility[END_REF][START_REF] Hoteit | Estimating human trajectories and hotspots through mobile phone data[END_REF].

A statistical model of completeness

As a further contribution to the characterization of the sparsity problem, we derive a data-driven model of CDR-based trajectory completeness. The model can be leveraged to obtain a quick, approximate estimate of the completeness of trajectories extracted from a CDR dataset, by merely providing the timespan of the CDR data and the desired temporal resolution.

In order to derive the model, we fit theoretical distributions to the completeness probability density function (PDF) and link the distribution variables to the system parameters T and τ . Earlier studies have shown that mobile communication activity tends to follow long-tailed distributions across users [START_REF] Seshadri | Mobile call graphs: Beyond power-law and lognormal distributions[END_REF]. Since the usage patterns of mobile devices directly determine completeness, it makes sense to look at heavytailed theoretical distributions. We run a maximum likelihood estimation of the variables of six standard long-tailed distributions (Generalized Pareto, Lévy, Powerlaw, Lognormal, Gamma, and Weibull) on the empirical PDF of completeness. We then evaluate the quality of fitting via the coefficient of determination R 2 , and by running the Kolmogorov-Smirnov test to obtain the D KS statistic.

Table 2 summarizes the results. Goodness-of-fit tests consistently point at either of two distributions, i.e., Weibull and lognormal, depending on the combination of system parameters. An illustrative example, for T = 60 days and τ = 60 minutes, is provided in Figure 3(a). The Weibull distribution emerges as a clear winner, as it prevails in 32 out of 40 cases in Table 2.

The PDF of the Weibull distribution is expressed as 3 show that the estimated Weibull parameters k and λ scale linearly (Pearson correlation coefficients of 0.89 and 0.98, respectively) with the most important system parameter, τ . Also, the observation time T introduces an offset to the linear relationship. We thus model k = [τ × (4.569 -0.01037T ) + 1086.5-1.14515T ]×10 -3 and λ = [τ ×(1.814-0.00379T )+26.23-0.07659T ]×10 -3 , and the CDF of CDR-based trajectory completeness as f Weibull (x; k, λ).

f Weibull (x) = k λ ( x λ ) k-1 e -(x/λ) k . Plots (b) and (c) of Figure

Missing location inference

In this section, we formalize the problem of trajectory reconstruction from sparse and irregular sampling and discuss how popular features of human mobility can be leveraged to design a sensible solution to the problem.

Trajectory reconstruction problem

Let us consider mobile phone data spanning an observation period T . User trajectories are inferred by discretizing T into intervals of duration τ , hence T = {1, • • • , N } can be regarded as the set of all time intervals. The actual (i.e., ground-truth) trajectory of a generic individual in the dataset is denoted by L T = {l i | i ∈ T }, where l i is the representative location of the user during the i th time slot. When a single mobile communication event is associated to a time interval i, the representative location l i maps to the geographic location associated with the event. This is the geographical position of the antenna that is the most representative of the user location. Although our discussion is general and can accommodate any definition of antenna representativeness, we opt for a strategy of selecting the antenna that handles the most communications during the target time slot i. We remark that this approach allows mitigating the so-called "ping-pong" (also referred to as "flickering" or "oscillation") effect that affects mobile phone data. This phenomenon consists in the mobile device association bouncing across nearby antennas, so that the location proxy for a same user changes even in absence of physical movement. The effect occurs at fast timescales of seconds [START_REF] Iovan | Moving and calling: Mobile phone data quality measurements and spatiotemporal uncertainty in human mobility studies[END_REF][START_REF] Katsikouli | Characterizing and Removing Oscillations in Mobile Phone Location Data[END_REF], hence oscillations are highly likely to occur within a single time slot i, which spans tens of minutes or more in our case. In each slot, oscillations are then canceled by the antenna selection process above.

As discussed in Section 3, trajectories extracted from mobile phone data are usually incomplete, i.e., no event is recorded during many time intervals, and l i = ∅ for a large fraction of i's. Let Ω = {i ∈ T | l i = ∅} be the set of time instant for which we have location information for the considered user. We define the observed trajectory of non-null locations as L Ω = {l i | i ∈ Ω}; this is the user mobility information that can be derived from the mobile phone data. Let us also denote by Ω C = T -Ω the time slots set during which the user location is unknown.

The trajectory reconstruction problem is then that of inferring, based on the sole knowledge of L Ω , an estimated complete trajectory LT

= { li = ∅ | i ∈ Ω ∪ Ω C } such that arg min LT 1 |T | N i=1 l i -li , (1)
where • denotes the geographic distance between the two argument locations. The expression in [START_REF] Blondel | A survey of results on mobile phone datasets analysis[END_REF] implies that the desired solution minimizes the mean absolute error (MAE) between the positions of the user in the complete ground-truth and reconstructed trajectories, at the same time instant.

Design rationale

In order to solve the problem in (1), we have to rely on the known locations of a user, L Ω , to infer those missing. Several well-known properties of human mobility play in our favor and can be leveraged to solve the problem effectively, as follows.

• Prevalence of static phases. People tend to spend a substantial amount of their time at a few fixed locations (e.g., workplace, school, shopping mall, sports center), where they linger for long continuous periods in the order of hours. Transitions among these locations are instead fairly rapid: in fact, people typically perceive movement phases between points of actual interest as a waste of time and strive to reduce them to a minimum. This results in a pattern of long static phases with fast movements in between [START_REF] Ficek | Inter-call mobility model: A spatio-temporal refinement of call data records using a gaussian mixture model[END_REF]. The prevalence of static behaviors allows adopting temporal resolutions that are granular enough to capture the important stay points of each individual, yet are still tractable for reconstruction. For instance, in the example of Figure 4, the user tends to stay for long times at the same location during working hours, and considering τ = 1 hour as done in the plot does not lose substantial positioning information during such hours. • Overnight invariance. A user is typically at the same location, i.e., at home, during nighttime, as also demonstrated by the fact that even sparse CDR data can be effectively employed to identify dwelling units [START_REF] Ranjan | Are call detail records biased for sampling human mobility[END_REF]. The observation that most nighttime locations match is important for trajectory reconstruction since mobile phone events tend to be especially scarce overnight, yet knowledge gaps can be filled with a limited amount of observed data during night hours.

As an illustrative example, the user in Figure 4 is always found at the same place early in the morning, late in the evening, and once overnight: it is apparent that such a location can be sensibly extended to all night hours, which are otherwise very poorly sampled. • Regularity of movement. Human mobility is strongly regular, from multiple perspectives. People frequently return to a same, limited set of locations [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF][START_REF] Mucelli Rezende Oliveira | On the regularity of human mobility[END_REF], which results in repetitive sequences of visits to a few places [7,[START_REF] Mucelli Rezende Oliveira | On the regularity of human mobility[END_REF], and infrequent trips to other locations in a (typically large) geographical region [START_REF] Mucelli Rezende Oliveira | On the regularity of human mobility[END_REF]. Regularity also occurs in time, as the locations above are visited in highly periodic patterns [START_REF] Song | Limits of predictability in human mobility[END_REF]. When considered jointly, these results have critical implications for trajectory reconstruction: they suggest that sparse location data may still be highly informative of the mobility of one individual if such data are sufficiently distributed in time to capture diverse moments of the day and week. This is typically the case for mobile phone data, which feature a combination of temporal irregularity of the sampling and long observation periods T . As an example, in the case of Figure 4, the fact that the user generates CDR events at different moments of the morning, working hours, and evening allows depicting a fairly complete picture of her mobility pattern during a typical day, by observing mobile communication samples for a sufficiently long amount of time (15 days in the considered sample).

Context-enhanced Trajectory Reconstruction

Context-enhanced Trajectory Reconstruction (CTR) is a novel trajectory reconstruction approach that solves the problem in (1), by receiving the observed trajectory L Ω as the input, and generating an estimated complete trajectory LT . To this end, CTR builds upon the observations in Section 4.2, via the three steps detailed next. The notation employed throughout the discussion is summarized in Section 8.

Nighttime trajectory reconstruction

The first step in CTR aims at reconstructing missing portions of a trajectory during the night hours. The rationale is that, as seen in Section 4.2, these are easily reconstructed from the extractable knowledge about users' home locations, which is also straightforward to derive from mobile phone data. Specifically, we assume that nighttime spans between 10 PM and 9 AM of the subsequent day: this period is adapted to the typical schedule in the Latin America country where our data is collected, but are easily adjusted to any other settings with statistical data on the local population habits. For each user in the dataset, we proceed as follows.

We first identify the most frequent location visited by the user during nighttime, which we name lH . Then, we verify that lH accounts for a conservative 80% (or more) of the positions of the user during the night period. If this is not the case, we deem that the result is not consistent enough to establish with high confidence the actual home location of the user, and we skip the nighttime trajectory reconstruction step. Otherwise, we identify lH as the home location of the user. When applied to our reference dataset of 1.8 million subscribers, this strategy allowed assigning a home location to 0.72 million users, equivalent to 42% of the whole population.

If a home location lH is identified, each 10 PM -9 AM period is processed separately for every day. The user is considered to be at home only in periods where no location other than lH is observed in the data: in this case, all time intervals i in the period are added to a set Ω H . Otherwise, if at least one location different from lH is present during a 10 PM -9 AM period, we consider that the user may not have spent that specific night at home, and no modification is made to Ω H . Once all days have been processed, we assign the home location to all slots in Ω H , i.e., li = lH , ∀i ∈ Ω H . We then update the known trajectory as LΩ = {l i | i ∈ Ω} ∪ { li | i ∈ Ω H }, and the set of missing locations by Ω C = T -Ω -Ω H . In our reference dataset, an average of 82% of the nighttime slots of the aforementioned 0.72 million users could be reconstructed this way.

Seamless trajectory reconstruction with tensor factorization

The second step of CTR aims at reconstructing the trajectory during the remaining time intervals in Ω C , by tailoring state-of-the-art tensor factorization techniques to our problem. Tensor factorization exploits redundancy to recover missing data; in our context, redundancy is created by the regularity of human movement, which creates repeated patterns of visited locations over the many days and weeks covered by a CDR dataset, as discussed in Section 4.2.

We start by organizing the trajectory data of a user into a three-dimensional tensor format, as illustrated in Figure 5. The location tensor X has dimensions that reflect the weekly, daily, and instantaneous (with granularity τ ) mobility of the user. Therefore X ∈ R Nw×N d ×2Nτ , where N w is the number of weeks in the observation period, N d is the number of days in a week, and N τ is the number of time slots in a day.

A trajectory L T is transformed into X via the following two steps:

• First, the trajectory is split into multiple one-day sub-trajectories. Each oneday sub-trajectory is then converted into a one-dimensional vector: for instance, the sub-trajectory of the j th day of the i th week is denoted by x i,j and satisfies x i,j = (l

(x) ni,j,1 , l (y)
ni,j,1 , . . . , l

n i,j,Nτ , l (x)
n i,j,Nτ ) ∈ R 2Nτ , where l (x) k (y)
and l (y) k are the two coordinates that identify the location at time step k, for l k = (l

(x) k , l (y) k ).
• Second, we enter all of the one-day sub-trajectories of a user trajectory into the location tensor by organizing them into a matrix of one-day sub-trajectories for all N w weeks in the observation period

T , i.e., X = [x i,j ] Nw×N d = {X i,j,k } Nw×N d ×2Nτ .
The location tensor X is partially filled with the original incomplete trajectory from mobile phone data, as well as with the inferred overnight locations as per Section 5.1. At this stage, the location tensor includes information for time slots in Ω ∪ Ω H . However, the tensor is still relatively sparse, as the set Ω C of time intervals where locations are missing is still large. In our reference dataset, the location tensors have an average density of 0.63 (where 1 denotes a complete tensor) when first generated.

In order to complete the location tensor of an individual, we employ tensor factorization, which decomposes the tensor into hyper-parameters and then uses them to infer missing values. Factorization has proven very effective in recovering highly incomplete structural data in other contexts [START_REF] Kong | Data loss and reconstruction in wireless sensor networks[END_REF][START_REF] Karatzoglou | Multiverse recommendation: N-dimensional tensor factorization for context-aware collaborative filtering[END_REF]. Formally, we formulate an equivalent tensor factorization (TF) problem as follows. First, we perform a canonical polyadic decomposition (CPD) [START_REF] Kolda | Tensor decompositions and applications[END_REF] of the tensor X ∈ R Nw×N d ×2Nτ into three R-rank matrices A ∈ R Nw×R , B ∈ R N d ×R , and C ∈ R 2Nτ ×R , where each value X i,j,k in the tensor is approximated as X i,j,k = R δ=1 A i,δ B j,δ C k,δ . In our experiments, we set the rank R to the minimum cardinality of the tensor dimensions, i.e., R = min(N w , N d , 2N τ ). This value entails a small number of hyperparameters in tensor decomposition, and is well within the upper bound of the tensor rank min(N w N d , 2N w N τ , 2N d N τ ), which cannot be otherwise computed with a finite algorithm [START_REF] Kolda | Tensor decompositions and applications[END_REF]. For simplicity, we employ the concise notation A, B, C for the CPD using matrices A, B and C. The TF problem is then Â, B, Ĉ = arg min

A,B,C (i,j,k)∈Ω∪Ω H (X i,j,k -A, B, C i,j,k ) 2 +λ( A 2 F + B 2 F + C 2 F ), ( 2
)
where λ is a penalty parameter to avoid overfitting [START_REF] Kolda | Tensor decompositions and applications[END_REF] and • F is the Frobenius norm. The solution to the problem in ( 2) is an estimation of the complete tensor X = Â, B, Ĉ , which thus recovers the missing locations in Ω C . In other words, we are seeking a tensor that has R-rank CPD decomposition and that fits the observed locations as closely as possible.

In fact, ( 2) is a standard TF problem, which does not account for contextual specificities that can improve the solution. Specifically, the formulation in (2) treats each dimension and every single value of the tensor equally, while human mobility patterns exhibit non-uniform importance of dimensions and values. Indeed, (i) daily periodicities tend to be stronger than weekly ones [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF], and (ii) consecutive time slots, hence nearby values in a same one-day sub-trajectory vector, show a strong correlation in terms of locations of the user [START_REF] Ficek | Inter-call mobility model: A spatio-temporal refinement of call data records using a gaussian mixture model[END_REF]. In the light of these considerations, we customize the TF problem for location inference, by introducing two additional elements to the optimization problem in [START_REF] Naboulsi | Large-scale mobile traffic analysis: A survey[END_REF].

The first element emphasizes daily repetitive mobility patterns by means of Toeplitz matrices. Toeplitz matrices allow modelling relationships between specific elements in the location tensor. For instance, given a matrix P = [p i,j ] N ×N and a Toeplitz matrix Q = Toeplitz(1, -1, 0) N ×N , then the product PQ 2 F represents the sum of the differences between consecutive values in the original matrix P, i.e., expands to

N i=1 N j=0 (p i-1,j -p i,j ) 2
F . Therefore, we construct a matrix D = Toeplitz(1, -1, 0) N d ×N d as illustrated in Figure 6; then, the expression X × d D 2 F , where × d is the tensor-matrix product [START_REF] Kolda | Tensor decompositions and applications[END_REF] such that (X × d D) i,m,k = N d n=1 X i,n,k D m,n , represents the sum of squared differences of location coordinates at the same hour of consecutive days. In our case, X is decomposed via CPD into A, B, C , hence we have A similar approach is then adopted to ensure that the solution to the TF problem favors the similarity of locations in consecutive time intervals. Here, we construct a Toeplitz matrix T = Toeplitz(1, 0, -1, 0) 2Nτ ×2Nτ , and denote by × τ the tensormatrix product such that (X × τ D) i,m,k = 2Nτ n=1 X i,m,n T k,n . Recall that each daily sub-trajectory is a one-dimensional vector (l

( A, B, C × d D) i,m,k = N d n=1 R δ=1 A i,δ B n,δ C k,δ D m,n = R δ=1 A i,δ ( N d n=1 D m,n B n,δ )C k,δ = R δ=1 A i,δ (DB) m,δ C k,δ (3)
(x) 1 , l (y) 1 , l (x) 2 , l (y)
2 , . . . ) ∈ R 2Nτ , hence the longitude (respectively, latitude) of two consecutive locations are separated by a coordinate of the other type. Then, including in the minimization problem the term X × τ T 2 F = A, B, TC 2 F allows constraining the diversity of locations between consecutive hours in a same day.

Embedding both elements above in (2) leads to an enhanced TF problem arg min

A,B,C (i,j,k)∈Ω∪Ω H (X i,j,k -A, B, C i,j,k ) 2 + λ( A 2 F + B 2 F + C 2 F ) + λ d A, DB, C 2 F + λ τ A, B, TC 2 F (4)
where λ d and λ τ weight the importance of similarity across time slots that are 24-hour apart and adjacent, respectively. The problem in ( 4) is a combination of multiple least square problems, and is efficiently solved by means of an alternating least square (ALS) technique [START_REF] Kolda | Tensor decompositions and applications[END_REF].

The complete tensor estimation X contains the original locations of the user extracted from the mobile phone data during time instants in Ω, the locations inferred during night hours in Ω H as per Section 5.1, and those estimated by solving (4) for all time slots in Ω C . Therefore, X can be used to fill an updated trajectory

LT = {l i | i ∈ Ω} ∪ { li | i ∈ Ω H ∪ Ω C }
, which is seamless, i.e., covers the whole observation period T .

Homogeneous quantization of locations

The locations in the original mobile phone data, as well as those inferred for night hours, match the position of cellular network towers. Instead, the solution of the enhanced TF problem in (4) returns locations whose coordinates are real values, which may not match the position of actual cellular towers. Since it is desirable that the reconstructed trajectory has consistent properties (including a homogeneous quantization of the locations that is based on the cellular infrastructure deployment), the third step in CTR aims at associating each location estimated as per Section 5.2 to the position of a real-world cellular tower.

Let C = {c j | 1 ≤ j ≤ N c } be the set of locations of the N c cell towers deployed in the target geographic region. Our goal is that of replacing the estimated locations li such that i ∈ Ω C by one of the cell tower positions c j ∈ C. To this end, we first design a sensible cost function for determining the suitability of a match ( li , c j ) between a generic location and a cell tower position. The cost function, denoted by f ( li , c j ), accounts for (i) the geographic distance between li and c j , and (ii) the popularity of the tower c j , as:

f ( li , c j ) =    1/e j if c j , li ≤ r j +∞ otherwise . , (5)
In expression [START_REF] Ahas | Modelling home and work locations of populations using passive mobile positioning data[END_REF] above: e j is the number of communication events that are generated by all the users in the dataset and associated to the cellular tower c j in the mobile phone data; • denotes the geographic distance between the two location parameters; and, r j is the cellular radius of tower c j , i.e., the maximum distance between the actual tower position c j and the farthest corner of its Voronoi tessellation cell [START_REF] Portela | Cellular network as a multiplicatively weighted voronoi diagram[END_REF] as illustrated in Figure 7.

The cell tower replacing an estimated location li is that satisfying

arg min cj ∈C f ( li , c j ). ( 6
)
In other words, minimizing the cost function limits the selection of matching cell towers for an estimated location li to those within distance r j of li ; then, it picks the candidate having the highest probability of occurrence in the data, hence adopting a wisdom-of-the-crowd approach that has already proven effective in personal mobility analysis [START_REF] Jeong | Cluster-aided mobility predictions[END_REF].

Validation

In order to validate the proposed CTR solution, we employ a set of ground-truth trajectories with completeness equal to 1, presented in Section 6.1. We generate incomplete trajectories from the ground-truth data, following the same sampling process observed in our reference CDR dataset, as described in Section 6.2. Then, complete trajectories are reconstructed with CTR from the downsampled ones, and compared against the ground-truth, in Section 6.3: the level of agreement lets us comment on the quality of the reconstruction.

Fine-grained ground-truth trajectories

We have access to data about the mobile Internet sessions of a small subset of the 1.8 million users in the CDR dataset introduced in Section 3. As more events are generated from Internet access rather than voice calls, they entail a sensibly higher frequency in location sampling. When considered jointly with CDR, rich mobile Internet data allows obtaining the complete trajectories of 1, 450 users over an observation period T of two consecutive weeks, and with a temporal resolution of τ = 1 hour. If multiple samples are recorded in the same time slot, the unique location is determined as that of the cell tower which most events are associated with. In other words, we have 1, 450 complete trajectories having observation sets Ω = T . In each of them, all the locations are known, i.e., L Ω = {l i | i ∈ Ω} = L T . We remark that the available fine-grained ground-truth trajectories bound our validation to a two-week period. We leave to future analyses, based on mobility data covering longer periods, the assessment of the impact of the trajectory duration on the accuracy of the completion process.

Incomplete trajectories

Incomplete trajectories are generated by downsampling the ground-truth ones. For instance, a trajectory with completeness 0.3 is derived from a two-week ground-truth trajectory having 14 × 24 = 336 time slots by keeping the locations in 101 time slots and discarding the rest. Here, the sampling process must be chosen appropriately: for instance, a uniform random selection would result in a temporal distribution of retained locations that is inconsistent with that produced by mobile communication events in CDR. We ensure that the incomplete trajectories mimic the temporal sampling properties of those extracted from CDR as follows: (i) we construct an empirical distribution of the probability that events are recorded at each time slots, based on CDR data of all 1.8 million users in our dataset; (ii) we set the desired completeness value k that shall characterize the output incomplete trajectory; (iii) we retain a fraction k of the overall location samples in one complete trajectory according to a random selection of time slots that follows the empirical distribution above. Performing this approach once leads to an artificial missing set Ω C , and results in a down-sampled incomplete trajectory

L Ω = {l i | i ∈ Ω} (recall that Ω = T -Ω C ).
This approach lets us produce a large amount of downsampled trajectories out of the 1, 450 ground-truth ones, and gives us enough flexibility to carry out a comprehensive validation. Specifically, we generate 3 × 10 5 incomplete trajectories with completeness ranging from 0.01 to 0.5, which allows testing CTR trajectory reconstructions in a wide diversity of settings.

Quality of CTR trajectory reconstruction

We run CTR on the incomplete CDR-like data and update the incomplete CDR-like trajectories as LT = {l i | i ∈ Ω}∪{ li | i ∈ Ω C }. We then assess the quality of the reconstructed trajectories LT against the complete ground-truth known locations L T .

The performance metric we adopt is cell displacement, which expresses the MAE used in (1) in terms of mobile network cells, as

∆( LT , L T ) = 1 |L T | -|L Ω | i∈T \Ω l i , li r j . ( 7
)
We recall that • denotes the geographic distance between the two location parameters, and r j is the cellular radius of the tower j to which the user is associated at time step i, i.e., for which l i = c j . Expression (7) implies that the cell displacement is computed only based on locations that are missing in the original data (i.e., for time instants in T \ Ω), and known user positions do not affect it. Unlike MAE, cell displacement compensates for the very diverse density (hence radius) that cellular towers tend to have in different regions of a large geographic area and provides a relative measure that is comparable across all locations.

The validation results are summarized in Figure 8, where the cell displacement measured in the trajectories reconstructed by CTR is shown against the completeness of the input CDR-like data. The candlesticks report the mean and median cell displacement, as well as the 10 th , 25 th , 75 th , and 90 th percentiles. The horizontal red line highlights a cell displacement of 1: below this value, the estimation error is lower than spatial precision of the original mobile phone data, i.e., the geographical coverage radius of the cell tower the user is presently associated to. Therefore, the error of the completion process is smaller than that inherent to the original data.

The most important remark is that the median cell displacement is always between 1 and 2, which denotes a substantial accuracy in the reconstructed trajectories: even under completeness levels as low as 1%, estimated locations are typically placed in a cell that is immediately adjacent to the correct one, or one cell further. More complete original trajectories provide CTR with added information to fill gaps in the data, leading to an improved performance where the typical estimation error is pushed closer to the correct cell. Interestingly, the gain of accuracy is higher at low completeness levels, and a trajectory with just 5% completeness already reduces the median cell displacement to less than 1.5. Some variability is observed around the median, which is however natural, given the heterogeneity that characterizes the mobility patterns and communication activities of different users.

The level of reconstruction accuracy achieved by CTR is acceptable for metropolitanscale analyses: in urban areas, network cells usually span a few hundreds of square meters and cover, e.g., (portions of) individual neighborhoods, hence the displacement in Figure 8 would still allow locating users fairly precisely and investigating mobility flows at an inter-neighborhood level. The reconstructed trajectory precision is instead excellent for regional-or national -scale studies, since cell displacements of 1 or 2 -over large surfaces covered by tens of thousands of cell towers -allow to capture human mobility at, e.g., inter-city level, perfectly.

Revisiting key results in the literature

As mentioned in Section 1, mobile phone data are the cornerstone of many seminal works on large-scale human mobility analysis. All such works (i) employ raw CDR datasets that are highly incomplete, and (ii) possibly adopt substantial user filtering. However, if and how trajectory incompleteness and user filtering effectand possibly bias -the outcome of these studies is unclear. Our aim in this section is shedding light on this matter, which also lets us demonstrate the importance of trajectory reconstruction in mobile phone data analysis.

To this end, we run CTR on the reference dataset presented in Section 3 and reconstruct all trajectories which have completeness higher than 1% in the original data. Selecting such a lower threshold on completeness bounds the typical cell displacement to a small value according to the validation in Section 6, hence ensure high confidence in the quality of the estimated locations; moreover, a 1% completeness threshold allows retaining 95% of the user population. Overall, we reconstruct the complete trajectories of around 1.7 million individuals in a large geographic region, with a temporal resolution of τ = 1 hour during three months.

We then reproduce three key studies on human mobility, using both the original CDR dataset and the complete trajectory data returned by CTR, and discuss the eventual differences that we observe in the results. The three analyses are separately presented next.

Laws of individual mobility

In an important study appeared in Nature in 2008, Gonzalez et al. [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF] derive a seminal model of individual human mobility based on CDR data. The data-driven analysis reveals that the travel distances of single users between two consecutive locations follow a truncated power-law. Formally, for any given individual, her travel distance ∆r between two consecutive events follows a truncated power-law P (∆r|r g ) ∼ r -α g F (∆r/r g ), where r g is the radius of gyration [1] , and F (x) is a function such that F (x) ∼ x -α for x < 1 and F (x) rapidly decreases for x > 1. This expression indicates that human movements resemble a Lévy flight that follows a power-law P (∆r) ∼ ∆r -α within a geographical region bounded by r g ; long displacements beyond such region are instead increasingly rare.

These key insights are obtained from a random sample of 1.67% of the total 6 million users, and the trajectories of the selected users are largely incomplete, with less than one location per user and per day recorded on average. We run the same analysis on the seamless trajectories of 1.7 million users, and show in Figure 9 how [1] The radius of gyration is a scalar metric assessing the overall mobility of a user. For a trajectory L T , it is computed as rg

= 1 |T | i∈T li, l C 2 , where l C = 1 |T | i∈T li
is the center of mass of all the locations in L T . considering complete movement data of a full user population affects the results. Clearly, we look now at a complete set of events -one per time slot -and the travel distance ∆r is computed between the locations in any two subsequent time slots. This ultimately results in a more comprehensive view of individual travel distances. We make the following considerations.

• As a preliminary result in their analysis, Gonzalez et al. find that the travel distances and radiuses of gyration aggregated over the whole user population follow truncated power laws (x + x 0 ) -β e -x/k . We confirm that this is also the case under complete trajectory data, as shown in the top plots of Figure 9.

The exponent values are also consistent, as Gonzalez et al. find β ∆r = 1.77 and β rg = 1.65, whereas β ∆r = 1.78 and β rg = 1.68 from our complete trajectory data. However, we remark a sensible difference in the cutoff values. Travel distances can be sensibly higher in complete data; also, k rg is at 400 km in [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF], which is not far from the 340 km from our original (i.e., incomplete) CDR-based trajectories but is reduced in the complete data, where the cutoff occurs instead at 175 km. We ascribe the difference to the fact that completion often leads to add missing locations that are far from the linear interpolation between CDR positions, such as those generated by infrequent but recurrent long-haul trips. Our conclusion is that, CDR data sparsity risks to both underestimate long trips, and overestimate the region within which the Lévy flight behavior of human mobility occurs. • Concerning individual movements, we can reproduce the truncated power-law behavior found by Gonzalez et al. with both our sparse CDR-based trajectories and the complete ones reconstructed via CTR. This is illustrated in the bottom plots of Figure 9. In the original work, the authors find an exponent α ≈ 1.2, and we estimate the α parameter at 1.8 and 1.4 for original and complete trajectories. These exponents are qualitatively consistent since they all are in the (1, 3) range that characterizes Lévy walks [START_REF] Shlesinger | Lévy Flights and Related Topics in Physics[END_REF]. • We ascribe the quantitative differences with respect to [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF] to the inherent specificity of each dataset. Although we also use locations from voice calls as in [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF], the data refer to very different cultural and economic settings (i.e., a European country and a South American one) and geographical span (our country covers a territory that is four times wider than the largest country in Europe). Overall, our results corroborate the findings in [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF], and show that the resulting conclusions hold also with seamless trajectory data from larger populations.

Uniqueness of individual trajectories

The work by De Montjoye et al. [START_REF] De Montjoye | Unique in the crowd: The privacy bounds of human mobility[END_REF], published in Scientific Reports in 2013, is the first to raise significant concerns on the privacy of users whose communication activities are recorded in CDR datasets. The study unveils the very high uniqueness of individual trajectories: knowing a few random time-stamped locations of an individual allows pinpointing her trajectory within the mobile phone data of a large user population. Notably, the authors show that 4 random l i 's identify one specific user among 1.5 million individuals 95% of times. In practice, the result implies that an adversary would need minimal knowledge about the mobility of a victim in order to perform a re-identification attack on a CDR dataset. The tests leading to the conclusion above consider a large user population but rely on raw CDR data that is highly incomplete: on average, a user only has approximately 19 unique locations observed in a month from 4, 200 cell towers.

Figure 10 depicts how complete trajectories affect the results. The plots show how many trajectories in the dataset are uniquely identified by knowing a random number of time-stamped locations (in the abscissa), with incomplete CDR data, and with seamless reconstructed trajectories. The result with incomplete data is very well aligned with that of De Montjoye et al.. However, the difference with complete trajectories is striking: for the low number of locations considered in [START_REF] De Montjoye | Unique in the crowd: The privacy bounds of human mobility[END_REF], the uniqueness is reduced by around one order of magnitude when the location information is seamless. The contrast lets us argue that the extremely high uniqueness identified by De Montjoye et al. is mainly caused by the diverse temporal patterns of the mobile communications of each user, rather than by his/her distinctive mobility. In other words, many users tend to have similar complete trajectories, and sampling them randomly yields limited uniqueness; it is instead the sparse sampling of the trajectories inherent to CDR that introduces an artificial temporal diversity and leads to dramatically increased uniqueness.

We thus conclude that actual human mobility is sensibly less unique than what suggested by incomplete trajectories inferred from mobile phone data. Note, however, that this does not mean that user privacy is preserved in reconstructed trajectory data, nor it negates the privacy warning by De Montjoye et al.. As shown in Figure 10, re-identification is still very possible and yet becomes more difficult. For instance, pinpointing a user in 90% of cases requires knowledge of 80 random spatiotemporal samples, instead of the 5 suggested by the original study.

Next-location predictability

The high theoretical predictability of human mobility was first exposed by Song et al. [START_REF] Song | Limits of predictability in human mobility[END_REF], in a famous paper that appeared in Science in 2010. In that work, the authors model the movement of each user as a time series of locations extracted from mobile phone data. They then measure the entropy rate of the probability of finding a particular time-ordered subsequence in the trajectory of the user. Finally, they compute the maximum predictability of each outcome from the entropy rate, by using Fano's inequality. The analysis considers 45,000 users (out of a total population of 10 million) with trajectory completeness higher than 0.2, at least 0.5 calls/hour and more than 2 unique recorded locations during an observation period of 14 weeks. Based on these data, the theoretical predictability is found to have a very high upper bound at around 93%.

We repeat the experiment with our reference dataset and trajectory reconstruction method. We run the filtering adopted by Song et al. on our original CDR data and extract 8, 000 users: this maps to around 0.44% of the available population, which is very close to the 0.45% considered in [START_REF] Song | Limits of predictability in human mobility[END_REF]. The probability distribution of the maximum predictability across the retained users is labeled as original in Figure 11. The mean predictability is at 81%, which is high yet quite far from the 93% found by Song et al.: As users are selected in the exact same way in the two cases, we ascribe the difference to the diverse mobility habits of the user populations in the two datasets, which are collected in different countries, at different scales, and during different time periods.

More interesting to our study is thus the comparison of the predictability obtained with the incomplete trajectories of 0.44% of the users, and that resulting from the seamless trajectories of the full available population (i.e., 1.7 million users), reconstructed through CTR. The curve labeled as complete in Figure 11 is considerably shifted to the right, with a much-reduced variance around the peak, now at 94%. Our hypothesis for this outcome is that the filtering introduced in [START_REF] Song | Limits of predictability in human mobility[END_REF] dramatically reduces the set of users, favoring individuals who are very active from a mobile communication viewpoint. It has been repeatedly demonstrated that interactions with the cellular network are more frequent for users with higher mobility [START_REF] Paul | Understanding traffic dynamics in cellular data networks[END_REF][START_REF] Couronné | Chatty mobiles:individual mobility and communication patterns[END_REF][START_REF] Hess | Exploring communication and mobility behavior of 3g network users and its temporal consistency[END_REF]. As a result, the sparse nature of CDR data lets Song et al. introduce an unwanted bias in their study, which is ultimately focused on very mobile individuals whose displacements are more difficult to anticipate. Instead, considering a much larger population lets us account for the vast majority of fairly static users, and reveals that people's movements are on average even easier to predict than estimated in [START_REF] Song | Limits of predictability in human mobility[END_REF].

In order to verify our hypothesis, we show in Figure 11 the predictability distribution for the complete, reconstructed trajectories of the same 8, 000 users that are retained from the original data by the filtering proposed in [START_REF] Song | Limits of predictability in human mobility[END_REF]. The result, tagged as complete/filtered is very consistent with that derived with the sparse CDR-based trajectories of those subscribers, i.e., the original curve. The result supports our intuition that it is the population sampling and not the reconstruction process that determines the striking difference in the figure.

Conclusions

We have investigated the problem of sparsity in trajectories inferred from mobile phone data. We have quantified the issue in a representative real-world dataset, proposed a novel solution, named CTR, to reconstruct a seamless trajectory from sparse CDR data, and validated our methodology using ground-truth movement patterns. We have also demonstrated the importance of trajectory reconstruction by revisiting well-known results on human mobility based on raw CDR, and showing that complete trajectories can in some cases affect the outcome of those analyses substantially.

List of abbreviations

Quantity

Notation Observation period of a trajectory T Ground-truth trajectory

L T = {l i | i ∈ T } Observed period of a trajectory Ω = {i ∈ T | l i = ∅} Observed trajectory L Ω = {l i | i ∈ Ω} Unknown period of a trajectory Ω C = T -Ω Estimated trajectory LT = { li = ∅ | i ∈ Ω ∪ Ω C } Empty set ∅ Location of the i th time slot l i

Funding

This work is supported by the EU FP7 ERANET program under grant CHIST-ERA-2012 MACACO and is performed in the context of the EMBRACE Associated Team of Inria and STIC AmSud MOTIf. Figure 11 Revisited results from [START_REF] Song | Limits of predictability in human mobility[END_REF]. Distributions of the entropy and maximum theoretical predictability derived with (i) the incomplete original CDR-based trajectories of 8, 000 users filtered with the method in [START_REF] Song | Limits of predictability in human mobility[END_REF], (ii) the complete trajectories of all 1.7 million users, and (iii) the complete trajectories of the same 8, 000 users considered in the first case.

Tables Table 1 The sample of CDR dataset. In this case, each CDR entry describes a voice call. All user identifiers and cell tower ids are coded to hide personal information. The GPS location of the cell towers is provided in the data source.

UserID

  and equivalently, A, B, C × d D = A, DB, C . Minimizing A, DB, C 2 F implies reducing the diversity of locations at a same hour across individual days.

FiguresFigure 1

 1 Figures

Figure 2

 2 Figure 2 Completeness versus temporal resolution. The left and right plots show the mean and median of the completeness of trajectories, respectively, as the function of temporal resolution τ .

Figure 3

 3 Figure 3 Fitting completeness into long-tailed distributions. (a) Empirical PDF of the completeness of 1.8 million user trajectories, with Weibull and lognormal theoretical fittings, for T = 60 days and τ = 60 minutes. (b)(c) Weibull variables k and λ versus the temporal resolution τ and observation period T .

Figure 4

 4 Figure 4 Example of the observed trajectory of a user. The trajectory is recored during T = 15 days and has a degree of completeness in 0.27. Each grid element maps to the representative location associated to each time interval of τ = 1 hour. Colors map to different places: the most active daytime and nighttime locations are marked in pink and blue respectively; white blocks are missing locations. Figure best viewed in colors.

Figure 5

 5 Figure 5 Example of the location tensor of a user X . A trajectory L T is represented as daily and weekly combinations of one-day sub-trajectories.

Figure 6

 6 Figure 6 Example of matrix T oeplitz(1, -1, 0) ∈ R N d ×N d . The diagonal is filled by 1's, the first upper diagonal is filled by -1's, and all other values are 0's.

Figure 7

 7 Figure7Example of cellular radius r j for cellular tower c j . Solid black lines outline the Voronoi tessellation of space, while the solid blue arrow denotes the cell radius r j . The dashed blue circle outlines the region whose radius is the cell radius r j . Because the region covers the estimated locations l1 , • • • , l3 , the cell tower c j is a possible candidate for associating these estimated locations to actual cell towers.

Figure 8

 8 Figure 8 Performance of CTR in cell displacement. Cell displacement of locations inferred via CTR with respect to the actual ones, versus the completeness of the original trajectory. Candlesticks highlight the mean (light blue), median (dark blue), 25 th and 75 th percentiles (box), and 10 th and 90 th percentiles (errorbars). The horizontal line (red) highlights a one-cell displacement.

Figure 9

 9 Figure 9 Revisited results from[START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF]. Distributions of travel distance ∆r, radius of gyration rg, and conditional travel distance ∆r (with respect to ∆r/rg) for original and complete trajectories.

Figure 10

 10 Figure 10 Revisited results from [5]. Ratio of trajectories with uniqueness |U | = 1 for a number of spatiotemporal points in the abscissa, for original (left) and complete (right) trajectories.

Table 2

 2 Empirical fittings of the completeness of trajectories. Fitting quality of six standard long-tailed distributions to the empirical PDF of trajectory completeness, in terms of R 2 and D KS . Rows refer to different combinations of observation period T and time resolution τ . Best-fit R 2 and D KS are highlighted in bold.

				Event Time		Cell Tower	Caller	Callee	In/Out Duration (s)
	38DA6	2015-05-01 18:26:50	1921		38DA6	163B7	Out		52
	78EC3	2015-05-01 14:16:09	2189		53808	78EC3	In		600
	9FAFE	2015-05-01 23:20:09	2189		9FAFE	7BBF1	Out		41
	708A2	2015-05-01 08:21:10	1988		96EC4	708A2	In		37
	A27AD 2015-05-01 21:51:09	2189		EA33F	A27AD	In		108
	9F3C7	2015-05-01 13:21:25	20102		C5691	9F3C7	In		134
	D4578	2015-05-01 17:03:46	20103		D4578	5B9A3	Out		30
	F904A	2015-05-01 23:24:03	1998		F904A	F3C88	Out		10
	4CCEA	2015-05-01 20:11:38	21104		4CCEA	5EF18	Out		438
	A77B8	2015-05-01 09:40:26	21104		A77B8	BD3E5	Out		33
	D5761	2015-05-01 20:34:40	1999		DAA24	D5761	In		600
	Duration Resolution	Weibull	Lognormal	Gamma	Pareto		Levy		Power law
	T	τ		DKS	R 2	DKS	R 2	DKS	R 2	DKS	R 2	DKS	R 2	DKS	R 2
		15min	0.0334	0.9990	0.0318	0.9973	0.3710	0.4679 0.0495 0.9969 0.2509 0.8046 0.3538 0.5031
	7d	30min 60min	0.0302 0.0278	0.9993 0.9990	0.0345 0.0372	0.9967 0.9958	0.0548 0.0475	0.9962 0.1716 0.8973 0.2649 0.7848 0.3587 0.4894 0.9977 0.1198 0.9516 0.2888 0.7506 0.4162 0.2041
		120min	0.0375	0.9972	0.0443	0.9938	0.0629	0.9853 0.1043 0.9768 0.3238 0.6977 0.2456 0.7450
		15min	0.0264	0.9986	0.0233	0.9984	0.3594	0.5007 0.0627 0.9893 0.2620 0.7853 0.3949 0.4120
	15d	30min 60min	0.0208 0.0188	0.9995 0.9997	0.0264 0.0279	0.9980 0.9974	0.0271 0.0257	0.9991 0.0724 0.9851 0.2765 0.7647 0.3576 0.4975 0.9987 0.0935 0.9719 0.2997 0.7315 0.3176 0.6076
		120min	0.0254	0.9987	0.0332	0.9961	0.0913	0.9572 0.1880 0.8780 0.3349 0.6792 0.2721 0.7116
		15min	0.0239	0.9985	0.0207 0.9985	0.3514	0.5261 0.0700 0.9835 0.2619 0.7872 0.3829 0.4365
	30d	30min 60min	0.0205 0.0149	0.9992 0.9996	0.0216 0.0263	0.9983 0.9975	0.0203 0.9991 0.1528 0.8976 0.2763 0.7661 0.3912 0.4149 0.0289 0.9982 0.0895 0.9702 0.3003 0.7346 0.3131 0.6212
		120min	0.0239	0.9984	0.0315	0.9962	0.0995	0.9479 0.1069 0.9538 0.3337 0.6893 0.3154 0.6176
		15min	0.0266	0.9980	0.0239	0.9977	0.1990	0.8284 0.0458 0.9934 0.2527 0.8076 0.3850 0.4156
	60d	30min 60min	0.0233 0.0207	0.9985 0.9985	0.0234 0.0264	0.9976 0.9970	0.0286 0.0310	0.9974 0.1944 0.8669 0.2649 0.7903 0.3897 0.4212 0.9980 0.0772 0.9769 0.2883 0.7622 0.3451 0.5635
		120min	0.0245	0.9975	0.0316	0.9958	0.0754	0.9750 0.0946 0.9617 0.3209 0.7196 0.3491 0.5070
		15min	0.0298	0.9954	0.0329	0.9954	0.3619	0.5248 0.0322 0.9954 0.2371 0.8390 0.3866 0.4528
	90d	30min 60min	0.0336 0.0305	0.9958 0.9960	0.0335 0.0355	0.9950 0.9944	0.0583 0.0454	0.9864 0.0398 0.9942 0.2487 0.8264 0.3819 0.4774 0.9913 0.0611 0.9843 0.2705 0.8020 0.3231 0.6123
		120min	0.0369	0.9940	0.0379	0.9932	0.0486	0.9927 0.0760 0.9743 0.2998 0.7687 0.3359 0.5885

Availability of data and materials

Due to privacy consideration and non-disclosure agreements with the data owners, we cannot make our datasets publicly available as they contain sensitive personal info and traces. We understand and appreciate the need for transparency in research. We are ready to help researchers who meet the criteria for accessing to confidential data with a signed confidentiality agreement. Please send your queries to Carlos Sarraute, at charles@grandata.com

Ethics approval and consent to participate Data collection was approved by the data owners, and written informed consent has been obtained for all study participants.

Competing interests

The authors declare that they have no competing interests.

Author's contributions

Designed the study: GC, AV, and MF. Collected and analyzed the data: GC and CS. All authors wrote and approved the final manuscript.

Author details 1 INRIA, Université Paris Saclay, 91120, Palaiseau, France. 2 École Polytechnique, Université Paris Saclay, 91120, Palaiseau, France. 3 CNR-IEIIT, 10129, Torino, Italy. 4 Grandata Labs, 550 15th Street, San Francisco, 94103, California, USA.