Erwan Bocher
email: erwan.bocher@univ-ubs.fr

Gwendall Petit

Nicolas Fortin

Judicaël Picaut

Gwenaël Guillaume

Sylvain Palominos

OnoM@p : a Spatial Data Infrastructure dedicated to noise monitoring based on volunteers measurements

Keywords: Noise, GIS, SDI, OGC, spatial analysis, crowdsourcing, VGI, smartphones

The present paper proposes an ideal Spatial Data Infrastructure (SDI) dedicated to noise monitoring based on volunteers measurements. Called OnoM@P, it takes advantage of the geospatial standards and open source tools to build an integrated platform to manage the whole knowledge about a territory and to observe its dynamics. It intends also to diffuse good practices to organize, collect, represent and process geospatial data in field of acoustic researches. OnoM@p falls within the framework of the Environmental Noise Directive (END) 2002/49/CE. The system relies on the NoiseCapture Android application developed for allowing each citizen to estimate its own noise exposure with its smartphone and to contribute to the production of community noisemaps.

INTRODUCTION

In 2015, more than 7 billion people (i.e. 96,4% of the world's population 1 ) owns a mobile phone subscription 2 . While in that year only 31% of sold phones were smartphones, the smartphone sales position in the phone market is widening with 446 million units sold only during the second quarter of 2015 3 . This growth affects all countries. For example, Deloitte company considers that "The smartphone market will grow by over 40% in 2015 with 70 million copies will be sold in Africa". These communication devices benefit from the technology advances in consumer electronics, and integrate an increasing number of sensors, including for the latest versions position sensors (GPS), motion sensors (accelerometer, gyroscope) and environmental sensors (microphone, camera, temperature sensor, photometer, barometer). In addition, the significantly growing computational power of smartphones and their Internet connectivity increase their attractiveness at once for the public, the public authorities and business.

In this context, the idea to rely on smartphones for initiating a new relation between the citizens and the environmental stakes is naturally emerging, for which the citizen would play a dual role of consumer and producer of environmental information. The participatory sensing concept rests on user-centric monitoring and environmental sensing by means of smartphones [START_REF] Burke | Participatory sensing[END_REF]. This notion recently arisen as a low-cost alternative to large-scale and costly infrastructures sensing based on sensor networks [START_REF] Cuff | Urban sensing: Out of the woods[END_REF]. Thus, numerous approaches were developed based on citizen-centric surveys (Andrew T. [START_REF] Campbell | People-centric urban sensing[END_REF], [START_REF] Adeel | Self-optimizing citizen-centric mobile urban sensing systems[END_REF] and opportunities to volunteer to take part into scientific research projects in environmental monitoring [START_REF] Conrad | A review of citizen science and community-based environmental monitoring: issues and opportunities[END_REF]. The OpenStreetMap project (OSM4 ) is the leading example of the effectiveness of citizen participatory with more than 2.8 million registered members by july 2016 5 . With the emerging demands of open geodata, OSM made one of the most used non-proprietary online maps and an input for numerous third-parties [START_REF] Neis | Analyzing the contributor activity of a volunteered geographic information project -the case of openstreetmap[END_REF]. The citizen can thus become a main actor by monitoring his environment and sharing his measurements or observations concerning his territory of life.

Several applications were recently developed for the purpose of noise levels measurements [START_REF] Leao | Monitoring exposure to traffic noise with mobile phones in china: A review of context[END_REF], as for example NoiseTube [START_REF] Maisonneuve | Participatory noise pollution monitoring using mobile phones[END_REF], [START_REF] D'hondt | Participatory noise mapping works! an evaluation of participatory sensing as an alternative to standard techniques for environmental monitoring[END_REF], [START_REF] Calabrese | Privacypreserving computation of participatory noise maps in the cloud[END_REF], WideNoise [START_REF] Becker | Awareness and learning in participatory noise sensing[END_REF], NoiseSpy [START_REF] Kanjo | Mobsens: Making smart phones smarter[END_REF], [START_REF] Kanjo | Noisespy: A real-time mobile phone platform for urban noise monitoring and mapping[END_REF], NoizCrowd [START_REF] Wisniewski | NoizCrowd: A Crowd-Based Data Gathering and Management System for Noise Level Data[END_REF], Noise-Watch [START_REF] Nugent | Noisewatch[END_REF], SoundOfTheCity [START_REF] Ruge | Soundofthecity -continuous noise monitoring for a healthy city[END_REF], etc. However, [START_REF] Guillaume | Noise mapping based on participative measurements[END_REF] point out some scientific and technological bottlenecks :

• mobile phones agents would mainly measure their daily sound exposure while the Directive 2002/49/EC fixed at least two indicators that are much less covered,

• lack of proper validation both for sound measures and GPS location,

• components of these applications does not rely sufficiently on standards, particularly those of the GIS community.

In this context, a Spatial Data Infrastructure (SDI) integrating a Geographic Information System (GIS) enables the use of spatial and environmental data in an efficient and flexible way and remains an ideal framework to fulfill this last point [START_REF] Mohammadi | Spatial Data Integrability and Interoperability in the Context of SDI[END_REF]. Indeed, since 2000s, the GIS community made progress and efforts to define open standards to enhance geospatial data exchange and sharing. These standards, mainly defined by the Open Geospatial Consortium are the main fuel to set up integrated platform, to facilitate the interconnection of systems and the carrying out of systemic approaches [START_REF] Steiniger | Free and Open Source GIS Software for Building a Spatial Data Infrastructure[END_REF]. They offers protocols, conceptual models and format to:

• store and query data : Simple feature access -Part 1 and 2, GeoSPARQL...

• style data : Symbology Encoding (SE) 6

• deliver data : Web Mapping Service (WMS) 7 , Web Feature Service (WFS) 8 ...

• search data : Catalogue Service (CSW) 9 , Gazetteer Service (WFS-G)...

• deliver and manage processing : Web Processing Service (WPS 10 ) and Web Coverage Processing Service (WCPS).

The present paper proposes an ideal Spatial Data Infrastructure (SDI) dedicated to noise monitoring based on volunteers measurements. Called OnoM@P, it takes advantage of the geospatial standards and open source tools to build an integrated platform to manage the whole knowledge about a territory and to observe its dynamics. It intends also to diffuse good practices to organize, collect, represent and process geospatial data in field of acoustic researches. Indeed, through this platform and its technical issues, several scientific locks are addressed :

• How to store and describe multiscale and temporal data collected from smartphones?

• What kind of language to query the data that makes sense for geographer, acousticians, developers... ?

• How to enable the integration between noise measures and GIS layers as OSM?

• How to aggregate data to produce human understandable and scientific results?

• How to map raw data or results that are scale and time dependant?

THE ONOM@P SDI

OnoM@p is an open source Spatial Data Infrastructure dedicated to noise monitoring based on volunteers measurement. OnoM@p is developed within the framework of the ENERGIC-OD 11 project that aims at deploying a set of Virtual Hubs (VH) based on a broker approach to offer to both end-users (through geoportals) and machines (web services, applications), unique and mutually consistent access points to heterogeneous data sources, including INSPIRE-compliant systems and Copernicus 12 /GMES 13 services.

System components and standards

The OnoM@p system is divided in two components (Figure 1) : The "NoiseCapture" application and the server side services (Virtual hub).

Figure 1. The OnoM@p system components.

The NoiseCapture application

NoiseCapture is an Android application to collect noise measure at given locations. Noise indicators are computed on the smartphone device during a period of time specified by the user. The main functionalities are the following (Figure 2) :

• Measurement (Figure 2(a)): Once the sound level calibration is done, the user start the measurement in order to record each second the LAeq (a weighted equivalent sound level), an average of the sound energy over a period of 1s. The noise spectrum repartition is also analysed using the Fourrier transform. The device locations are recorded while measuring the sound level. The user has also the ability to provide his own feedback about the feeling of the noise environment.

• Extended report: Advanced statistics are computed locally on the phone and shown to the user (Figure 2(b)). For each user's measurement the locations of the noise levels are displayed in a map (Figure 2(c)).

• Share results with the community: Anonymized results are transferred to Virtual Hubs (web server) and post-processed in order to build a noise map that merge all community results.

The server side services

The server side services are hosted on the French part of the Energic-OD Virtual Hub. They are organized around a Geoservices component, a relational database management system and an experimental processing module.

The Geoservices

The Geoservices component uses the GeoServer application and a plugin to support Web Processing Services (WPS). This WPS plugin is used to push from the smartphone to a relational database a The zipped file contains 3 files:

• a geojson file, that stores all track coordinates and the noise indicators computed along the track,

• a metadata file that concatenates some metadata concerning the smartphone and the measurement such as the calibration method,

• a readme file that contains some details, definitions about the noise indicators.

Once the data uploaded, a 4 sub-steps process is run, thanks to the Groovy language:

• Extract: First of all, the zipped file is analysed to avoid security problems, then the file is unzipped and 3 files are stored in a temporary folder,

• Transform: The geojson and the metadata files are parsed and transformed into SQL scripts,

• Export: The SQL scripts are uploaded on a Noise database,

• Archive: The input zipped file is archived into a specific folder.

The rendering module of the Geoserver application is connected with the Noise database to offer on demand a set of maps using the WMS standard and the raw data using the WFS / GeoJson standard.

The relational database management system Two spatial databases are built to store the data collected and needed by the OnoM@p system. These databases use the RDMS PostgreSQL14 with the Postgis15 extension.

The first database stores all data collected from the NoiseCapture application and some experimental indicators (i.e. noise population exposure). The second one is shared by the French Energic-OD hub and delivers common data set required by the applications. For example, the French statistical geographic zones (IRIS) or the French census statistics (INSEE 16 ).

The experimental processing module

This module is implemented on top of the H2GIS database [START_REF] Bocher | H2GIS a spatial database to feed urban climate issues[END_REF] and the OrbisGIS WPS server [START_REF] Bocher | ORBISGIS: Geographical Information System Designed by and for Research[END_REF]. It contains an experimental chain of spatial and statistical analysis methods to compute new indicators from research work. This module is compliant with the last version of the WPS standard (WPS 2.0).

Operating mode

The OnoM@p system is an integrated platform which covers all workflow and lifecycle stages to collect, compute, share and display noise data and at the end inform the user about its noise exposure. In the OnoM@p ecosystem, we consider four levels of stakeholders who will play a role (Figure 3).

VGI (Volunteered geographic information)

The volunteer collects noise data from its smartphone and publishes it on the OnoM@p database. The volunteer must respect some guidelines to obtain a good measure. This guidelines are not constrained by the application. It's up to the volunteer to follow best practices such as:

• the smartphone has to be outside the clothe's pockets,

• the smartphone has to be in the hand (e.g. not on a table),

• the microphone has to be free (no mask),

• the volunteer should walk (instead of running or biking),

• the measure must be done outside a building,

• the measure must be done when the GPS signal is correct (e.g. not in a tunnel, ...).

Nevertheless, the WPS operation in charge of indicators computation integrates a post-process analysis to control the quality of the measure. The uncertainty of the microphones and the precision of the GPS locations are estimated, taking account of the time period of each records. Finally bad located measures are excluded from the chain of indicators and noise maps.

The Decision Maker

The Decision Maker is a stakeholder who can be in a public authority or in a company. He is using the visualization services (web map viewer) to see and understand the noise map exposure on its territory. He can also use this kind of information to take or support a decision.

The Public

The Public is made of citizens or of civil society organisations. They mainly use the visualization services to be aware of noise issues.

The Expert

The expert is someone (geographer, acoustician . . . ) who is able to understand and manage the raw data, extracted from the database, and to use it in its domain of application to produce added value data (e.g. identify places where the noise exposure is too high, or analyse the evolution of the noise exposure on several years...).

Mapping sound indicators

Both acoustical and statistical indicators are computed for the whole duration of the measurement. E.g. L A10 that represents the peak noise, L A50 the mean value of the noise levels and L A90 the mean value of the background noise (Figure 2 The statistical indicators can also be evaluated by connecting the noise data with geographical and population data. For example, population exposure can be reckoned by confronting noise data with land-use rates data. Thus, an estimation of the inhabitants exposed to ranged noise levels (i.e. the distribution of noise exposure in the population) is allowed by combining the sound levels (e.g. the Lden values) comprised within a range of values over the studied territory and its demographic distribution (Figure 5). Compute indicators from collected smartphone data is one of the big challenge of the OnoM@p platform and need new methods and techniques:

• to express evolution and changes,

• to build consistent analysis of spatial and temporal variation,

• to find the best geographic unit to maximize cross-temporal comparability.

CONCLUSION

The OnoM@p Spatial Data Infrastructure is a very promising tool for the environmental noise monitoring. It offers solutions both for the experts and the community, in order to evaluate the noise annoyance, as well as the noise exposure. Comparing with classical noise evaluation methods, based on numerical simulations (with a limited number of sound sources, in a simplified urban area, with approximated noise propagation models), the proposed methodology enabled to present the 'real' state of the noise exposure, based on real measurements. In addition, since measurements can be carried out everywhere, it allows to produce a noise evaluation that is not restricted to urban areas.

In an acoustics point of view, one challenge (i.e. one criticism too) is still the quality of the noise measurement, since smartphone capabilities can not be compared to professional sonometers.

Although recent technical developments seems to show that, in a few years, this will no longer be a problem, one can proposed original solutions to post-process measurements, using for example a cross-calibration procedure of all mobile measurements, using an in-situ sensors network [START_REF] Can | Cross-calibration of participatory sensor networks for environmental noise mapping[END_REF]. Beyond the acoustic measurement, the last challenge could be to offer to the community and to the experts, simple and relevant representations of sound environment, which are not limited to noise indicators (i.e. physical measurement), but to the perceived environment. The choice of appropriate and simple spatial representations, associated with the ability to cross acoustic indicators with other geo-referenced data, clearly shows that geomatics has become today an essential component in the assessment, thus the reduction, of noise in the environment.

FUNDING STATEMENT

The ENERGIC-OD Project (European Network for Redistributing Geospatial Information to user Communities -Open Data) is partially funded under the ICT Policy Support Programme (ICT PSP) as part of the Competitiveness and Innovation Framework Programme by the European Community. The present work is also supported by the French geographic portal G ÉOPAL of the Pays de la Loire region.

Figure 2 .

 2 Figure 2. Screenshots of the NoiseCapture application: (a) measurement screen, (b) results screen and (c) noise level visualization on a map.

Figure 3 .

 3 Figure 3. The four levels of stakeholders.

  (b)). A first native map is rendered to represent a classification of five sound levels as follows (Figure 4): [<45 dB(A)], [45-55 dB(A)], [55-65 dB(A)], [65-75 dB(A)], [>75 dB(A)]

Figure 4 .

 4 Figure 4. Aggregated WMS layers : OSM and Noise measurements in the italian city of Lecco.

Figure 5 .

 5 Figure 5. Aggregated WMS layers: OSM and noise population exposure, at IRIS (aggregated units for statistical information) scale in the french city of Nantes.

https://www.openstreetmap.org/ (accessed in july 2016)

http://www.openstreetmap.org/stats/data stats.html (accessed in july 2016)

https://www.postgresql.org/ (accessed in july 2016)

http://postgis.net/ (accessed in july 2016)

http://www.insee.fr/en/default.asp (accessed in july 2016)


Keywords: propagation de croyances, approximation de Bethe, reconstruction de trafic, prédiction, systèmes de transport intelligent, véhicules traceurs

Propagation de croyances et approximation de Bethe pour la prédiction de trafic

Résumé : On définit et étudie un algorithme de reconstruction utilisant l'algorithme Belief Propagation (propagation de croyances, BP) et l'approximation de Bethe. L'idée est d'encoder dans un graphe des données a priori composées de corrélations ou de lois marginales et d'utiliser une procédure de passage de messages pour estimer l'état réel à partir d'informations temps-réel. Cette méthode, développée pour des besoins de prédiction de trafic, est particulièrement adaptée au cas où la seule information disponible provient de véhicules sonde (Floating Car Data). Nous proposons une discrétisation binaire du trafic s'appuyant sur le modèle d'Ising de physique statistique, permettant de reconstruire et de prédire le trafic en temps réel. Des propriétés générales de l'algorithme BP sont discutées dans ce contexte. En particulier une étude détaillée des propriétés de stabilité fonction des données a priori et de la topologie du graphe est fournie. Une étude numérique sur un modèle de trafic simplifié permet d'illustrer le fonctionnement de l'algorithme. La façon de généraliser cette approche pour encoder une superposition de plusieurs états de trafic est discutée.

Introduction

With an estimated 1% GDP cost in the European Union (i.e. more than hundred billions euros), congestion is not only a time waste for drivers and an environmental challenge, but also an economic issue. This is why the European commission financed the REACT project, where new traffic prediction models have been developed. These predictions are to be used to inform the public and possibly to regulate the traffic.

Today, some urban and inter-urban areas have traffic management and advice systems that collect data from stationary sensors, analyze them, and post notices about road conditions ahead and recommended speed limits on display signs located at various points along specific routes. However, these systems are not available everywhere and they are virtually non-existent on rural areas. With rural road crashes accounting for more than 60% of all road fatalities in OECD (Organization for Economic Cooperation and Development) countries, the need for a system that can cover these roads is compelling if a significant reduction in traffic deaths is to be achievable.

The REACT project combines a traditional traffic prediction approach on equipped motorways with an innovative approach on non-equipped roads. The idea is to obtain floating car data from a fleet of probe vehicles and reconstruct the traffic conditions from this partial information. To understand why it is not possible to fuse these two parts, we have to go a bit more into prediction algorithms details.

Two types of approaches are usually distinguished, namely data driven (application of statistical models to a large amount of data, for example regression analysis) and model based (simulation or mathematical models explaining the traffic patterns). As we stated before, the choice is largely led by the availability of data. In our case, since little data is available on non-equipped roads (only the equipped vehicles driving along the observed roads), the model driven approach is the only feasible one. For more information about traffic prediction methods, we refer the reader to [START_REF] Benz | Information supply for intelligent routing services -the INVENT traffic network equalizer approach[END_REF][START_REF]Technology assessment and expected targets[END_REF][START_REF] Versteegt | PredicTime -state of the art and functional architecture[END_REF].

Most current traffic models are deterministic, described either at a macroscopic level by a set of differential equations linking variables such as flow and density, or by Newton's law at a microscopic level where each individual car is considered. Intermediate descriptions are essentially kinetic models, like for example cellular automata [START_REF] Nagel | A cellular automaton model for freeway traffic[END_REF], which are very well adapted to freeway traffic modeling and adapted to some extent to urban traffic modeling [START_REF] Essen | Microscopic simulation of urban traffic based on cellular automata[END_REF]. Traffic flow models are quite adapted and efficient on motorways where fluid approximation of the traffic is reasonable; they tend to fail for cities or rural roads. The reason is that the velocity flow field is subject to much greater fluctuations induced by the nature of the network (presence of intersections and short distance between two intersections) than by the traffic itself. These fluctuations are both spatial and temporal (a red or green traffic light at a cross-road, a road-work, etc). There is no local stationary regime for the velocity, the dynamics are dominated by the fluctuations.

We propose in this paper a hybrid approach in the continuation of [START_REF] Furtlehner | Belief-propagation algorithm for a traffic prediction system based on probe vehicles[END_REF], by taking full advantage of the statistical nature of the information, in combination with a stochastic modeling of traffic patterns. In order to reconstruct the traffic and make predictions, we propose a model-the Bethe approximation (BA)-to encode the statistical fluctuations and stochastic evolution of the traffic and an algorithm-the belief propagation (BP) algorithm-to decode the information. Those concepts are familiar to the computer science and statistical physics communities since it was shown [START_REF]Constructing free-energy approximations and generalized belief propagation algorithms[END_REF] that the output of BP is in general the Bethe approximation [START_REF] Bethe | Statistical theory of superlattices[END_REF].

The paper is organized as follows: Section 2 describes the model and its relationship to the Ising model and the Bethe approximation. The inference problem and our strategy to tackle it using the Belief Propagation algorithm are stated in Section 3. The implementation of these ideas requires some new results about the BP algorithm, which are the subject of Section 4; this concerns in particular the effect of the normalization of the messages, the parameterization of the model and the stability of the fixed points. Section 5 is devoted to implementation details of the decoding algorithm and to some numerical results illustrating the method. Finally, some new research directions are proposed in Section 6.

Traffic description and statistical physics

The graph onto which we apply the belief propagation procedure is made of spacetime vertices that encode both a location (road link) and a time (discretized on a few minutes scale). More precisely, the set of vertices is V = L ⊗ Z + , where L corresponds to the links of the network and Z + to the time discretization. To each point α = (ℓ, t) ∈ V, we attach an information τ α ∈ {0, 1} indicating the state of the traffic (1 if congested, 0 otherwise). Each cell is correlated to its neighbors (in time and space) and the evaluation of this local correlation determines the model. In other words, we assume that the joint probability distribution of τ The homogeneous Ising model (uniform coupling constants) is a well-studied model of ferro (positive coupling) or anti-ferro (negative coupling) material in statistical physics. It displays a phase transition phenomenon with respect to the value of the coupling. At weak coupling, only one disordered state occurs, where spins are randomly distributed around a mean-zero value. Conversely, when the coupling is strong, there are two equally probable states that correspond to the onset of a macroscopic magnetization either in the up or down direction: each spin has a larger probability to be oriented in the privileged direction than in the opposite one.

V def = {τ α , α ∈ V} ∈ {0, 1} V is of the form p({τ α , α ∈ V}) = α∈V φ α (τ α ) (α,β)∈E ψ αβ (τ α , τ β ) (2.1) INRIA (a) (b)
From the point of view of a traffic network, this means that such a model is able to describe three possible traffic regimes: fluid (most of the spins up), congested (most of the spins down) and dense (roughly half of the links are congested). For real situations, we expect other types of congestion patterns, and we seek to associate them either to the p-state Potts model if we extend the binary to p-ary description, or to the possible states of an inhomogeneous Ising model with frustration (i.e. with possibly negative coupling parameters), referred as spin glasses in statistical physics [START_REF] Mézard | Spin glass theory and beyond[END_REF]. When such a system is frustrated because some negative couplings, leading to a certain number of contradictions, a proliferation of meta-stable states occurs, which eventually scales exponentially with the size of the system.

On a simply connected graph, the knowledge of the one-vertex and two-vertices marginal probabilities is sufficient [START_REF] Pearl | Probabilistic reasoning in intelligent systems: Network of plausible inference[END_REF] to fully determine the measure (2.1).

p(τ V ) = (α,β)∈E pαβ (τ α , τ β ) α∈V p(τ α ) qα-1 = α∈V pα (τ α ) (α,β)∈E pαβ (τ α , τ β ) pα (τ α )p β (τ β ) , (2.2)
where q α denotes the number of neighbors of α. Since our space time graph G is multi-connected, this relationship between local marginals and the full joint probability measure can only be an approximation, which in the context of statistical physics is referred to as the Bethe approximation. This approximation is provided by the minimum of the so-called Bethe free energy, which, based on the form (2.2), is an approximate form of the Kullback-Leibler distance,

D(b(τ V ) p(τ V )) def = τ V b(τ V ) ln b(τ V ) p(τ V ) ,
and which rewrites in terms of a free energy as

D(b(τ V ) p(τ V )) = F(b(τ V )) -F(p(τ V )),
where

F(b(τ V )) def = U(b(τ V )) -S(b(τ V )), (2.3)
with the respective definitions of the energy U and of the entropy

S U(b(τ V )) def = - (α,β)∈E b αβ (τ α , τ β ) log ψ αβ (τ α , τ β ) - α∈V b α (τ α ) log φ α (τ α ), S(b(τ V )) def = - (α,β)∈E b αβ (τ α , τ β ) log b αβ (τ α , τ β ) - α∈V b α (τ α ) log b α (τ α ).
In practice, what we retain from an Ising description is the possibility to encode a certain number of traffic patterns in a statistical physics model. This property is actually shared also by the Bethe Approximation (BA) and this is the reason for us to directly encode the traffic patterns in a BA rather than the inhomogeneous Ising model itself, based on historical data, and to avoid therefore an intermediate approximation step. BA simply provides us with a set of marginals probabilities that we INRIA try to match with the historical data. But this set, which is the result of an iterative procedure, is not necessarily unique (see for example [START_REF] Mooij | On the properties of the Bethe approximation and loopy belief propagation on binary network[END_REF]) and the proliferation of possible solutions depends on the frustration induced by the historical correlations used to define the ψ's of (2.1). The setting of our model consists therefore into an optimization procedure of the matching between the set of historical values obtained from probe vehicles and the set of marginal probabilities of the BA.

The data collected from the probe vehicles is used in two different ways. The most evident one is that the data of the current day directly influences the prediction. In parallel, this data is collected over long periods (weeks or months) in order to estimate the model (2.1). Typical historical data that is accumulated is • pα (τ α ): the probability that vertex α is congested (τ α = 1) or not (τ α = 0);

• pαβ (τ α , τ β ): the probability that a probe vehicle going from α to β ∈ V(α) finds α with state τ α and β with state τ β .

The edges (α, β) of the space time graph G are constructed based on the presence of a measured mutual information between α and β, which is the case when pαβ (τ α , τ β ) = pα (τ α )p β (τ β ).

3 The reconstruction and prediction algorithm

Statement of the inference problem

We turn now to our present work concerning an inference problem, which we set in general terms as follows: a set of observables τ V = {τ α , α ∈ V}, which are stochastic variables are attached to the set V of vertices of a graph. For each edge (α, β) ∈ E of the graph, an accumulation of repetitive observations allows to build the empirical marginal probabilities {p αβ }. The question is then: given the values of a subset τ V * = {τ α * , α * ∈ V * }, what prediction can be made concerning V * , the complementary set of V * in V? There are two main issues:

• how to encode the historical observations (inverse problem) in an Ising model, such that its marginal probabilities on the edges coincide with the pαβ ?

• how to decode in the most efficient manner-typically in real time-this information, in terms of conditional probabilities P (τ α |τ V * )?

The answer to the second question will somehow give a hint to the first one.

The belief propagation algorithm

BP is a message passing procedure, which output is a set of estimated marginal probabilities, the beliefs b αβ [START_REF] Pearl | Probabilistic reasoning in intelligent systems: Network of plausible inference[END_REF]. The idea of the BP algorithm is to factor the marginal probability at a given site in a product of contributions coming from neighboring sites, which are the messages. The messages sent by a vertex α to β ∈ V(α) depends on the messages it received previously from other vertices:

m α→β (τ β ) ← τα∈{0,1} n α→β (τ α )φ α (τ α )ψ αβ (τ α , τ β ), (3.1)
where

n α→β (τ α ) def = γ∈V(α)\{β} m γ→α (τ α ). (3.2)
In practice, the messages will be normalized so that

τ β ∈{0,1} m α→β (τ β ) = 1. (3.3)
We will come back to the effects of this in Section 4.2.

The output of the algorithm is a set of beliefs, which are an approximation of the one-vertex and two-vertices marginals of p(τ V ). The beliefs b α are reconstructed according to

b α (τ α ) ∝ φ α (τ α ) m β→α (τ α ), (3.4)
and, similarly, the belief b αβ of the joint probability of (τ α , τ β ) is given by

b αβ (τ α , τ β ) ∝ n α→β (τ α )n β→α (τ β ) × φ α (τ α )φ β (τ β )ψ αβ (τ α , τ β ). (3.5)
In the formulas above and in the remainder of this paper, the proportionality symbol ∝ indicates that one must normalize the beliefs so that they sum to 1.

A simple computation shows that equations (3.4) and (3.5) are compatible, since

(3.1)-(3.2) imply that τα∈{0,1} b αβ (τ α , τ β ) = b β (τ β ).
It has been realized a few years ago [START_REF] Yedidia | Generalized belief propagation[END_REF] that the fixed points of the BP algorithm coincide with local minima of the Bethe free energy (2.3). This justifies that we can use this algorithm to approximate our Ising model.

We propose to use the BP algorithm for two purposes: estimation of the model parameters (the functions φ and ψ) from historical data and reconstruction of traffic from current data.

INRIA

Setting the model with Belief Propagation

The fixed points of the BP algorithm (and therefore the Bethe approximation) allow to approximate the joint marginal probability p αβ when the functions ψ αβ and φ α are known. Conversely, it can provide good candidates for ψ αβ and φ α from the historical values pαβ and pα .

To set up our model, we are looking for a fixed point of the BP algorithm satisfying

(3.1)-(3.2) and such that b αβ (τ α , τ β ) = pαβ (τ α , τ β ) and therefore b α (τ α ) = pα (τ α ).
It is easy to check that the following choice of φ and ψ,

ψ αβ (τ α , τ β ) = pαβ (τ α , τ β ) pα (τ α )p β (τ β ) , (3.6)
φ α (τ α ) = pα (τ α ), (3.7)
leads (2.1) to coincide with (2.2). They correspond to a normalized BP fixed point for which all messages are equal to 1/2. There is however no guarantee that this fixed point is a stable fixed point; actually, for an Ising-type system below the critical temperature, we often observe that this point is unstable (see Section 5 for the simulation results). It will be shown however in Section 4.1 that this form of φ and ψ is in some sense canonical.

Traffic reconstruction and prediction

Let V * be the set of vertices that have been visited by probe vehicles. Reconstructing traffic from the data gathered by those vehicles is equivalent to evaluating the conditional probability

p α (τ α |τ V * ) = p α,V * (τ α , τ V * ) p V * (τ V * ) ,
where τ V * is a shorthand notation for the set {τ α * } α * ∈V * . The BP algorithm applies to this case if a specific rule is defined for vertices α * ∈ V * : since the value of τ α * is known, there is no need to sum over possible values and (3.1) becomes

m α * →β (τ β ) ← n α * →β (τ α * )φ α * (τ α * )ψ α * β (τ α * , τ β ).
The resulting algorithm is supposed to be run in real time, over a graph which corresponds to a time window (typically a few hours) centered around present time, with probe vehicle data added as it is available. In this perspective, the reconstruction and prediction operations are done simultaneously on an equal footing, the distinction being simply the time-stamp (past for reconstruction or future for prediction) of a given computed belief. The output of the previous run can be used as initial messages for a new run, in order to speedup convergence. Full re-initialization (typically a random set of initial messages) has to be performed within a time interval of the order but smaller than the time-scale of typical traffic fluctuations.

Some general properties of the Belief Propagation algorithm

This section contains several theoretical results on the BP algorithm. Although they are stated in the context of Section 2, most of these results can be trivially extended to a general factor graph and variables taking more than two values (transforming the Ising model into a Potts model), except possibly for Section 4.3.

Building the model from its fixed points

The particular use that we make of the Bethe approximation, as outlined in Section 3.3, means that the output of the algorithm takes precedence over the underlying Ising model, which is an unusual situation. The following propositions shows how to estimate φ α and ψ αβ from the historical values pα and pαβ .

Let us start with a direct consequence of the BP fixed point equations. The following straightforward proposition extends (2.2) to the case of a non-tree structure.

p(τ V ) = α,β b αβ (τ α , τ β ) α b qα-1 α (τ α ) = α∈V b α (τ α ) (α,β)∈E b αβ (τ α , τ β ) b α (τ α )b β (τ β ) .
Proof. This is a simple consequence of (3.4)- (3.5).

What this proposition means is that different BP fixed points correspond to different factorizations of the joint measure (2.1). The knowledge of a set of beliefs is thus sufficient to determine the underlying Ising model and consequently the other fixed points of the algorithm.

The following proposition gives more insight on how the different components of (2.1) can be written in terms of the BP fixed points. Then the following equalities hold

INRIA

ψ αβ (τ α , τ β ) = pαβ (τ α , τ β ) pα (τ α )p β (τ β ) m α→β (τ β )m β→α (τ α ), (4.2)
φ α (τ α ) = pα (τ α ) β∈V(α) m β→α (τ α ) . (4.3)
Conversely, assume that there exist boolean functions µ αβ (τ β ) such that

ψ αβ (τ α , τ β ) = pαβ (τ α , τ β ) pα (τ α )p β (τ β ) µ αβ (τ β )µ βα (τ α ), (4.4)
φ α (τ α ) = pα (τ α ) β∈V(α) µ βα (τ α ) . ( 4

.5)

Then m α→β = µ αβ is a fixed point of the BP algorithm and (4.1) holds.

Proof. Relation (4.2) is obtained by rewriting (3.4) and (3.5) as

ψ αβ (τ α , τ β ) = b αβ (τ α , τ β ) φ α (τ α )n α→β (τ α )n β→α (τ β )φ β (τ β ) = b αβ (τ α , τ β ) b α (τ α )b β (τ β ) m α→β (τ β )m β→α (τ α ). (4.6)
To prove the second assertion, the first step is to show that µ αβ is a BP fixed point:

τα∈{0,1} γ∈V(α)\{β} µ γα (τ α ) φ α (τ α )ψ αβ (τ α , τ β ) = τα∈{0,1} pα (τ α ) µ βα (τ α ) pαβ (τ α , τ β ) pα (τ α )p β (τ β ) µ αβ (τ β )µ βα (τ α ) = τα∈{0,1} pαβ (τ α , τ β ) pβ (τ β ) µ αβ (τ β ) = µ αβ (τ β ).
For this fixed point, (3.5) reduces to (4.1), which concludes the proof of the proposition.

While Proposition 4.2 seems to indicate that there is some leeway in choosing ψ αβ , a proper change of variables shows that all the choices are equivalent. Let us define the following new set of messages

x α→β (τ β ) def = m α→β (τ β ) µ αβ (τ β ) .
Equation (3.1) then becomes

x α→β (τ β )µ αβ (τ β ) = τα∈{0,1} γ∈V(α)\{β} x γ→α (τ α )µ γα (τ α ) φ α (τ α )ψ αβ (τ α , τ β ) = τα∈{0,1} γ∈V(α)\{β} x γ→α (τ α ) pα (τ α ) µ βα (τ α ) pαβ (τ α , τ β ) pα (τ α )p β (τ β ) µ αβ (τ β )µ βα (τ α ) = τα∈{0,1} γ∈V(α)\{β} x γ→α (τ α ) pαβ (τ α , τ β ) pβ (τ β ) µ αβ (τ β ),
and therefore

x α→β (τ β ) = τα∈{0,1} γ∈V(α)\{β} x γ→α (τ α ) pαβ (τ α , τ β ) pβ (τ β ) .
This version of the BP algorithm is thus equivalent to the heuristic choice (3.6)-(3.7), which corresponds to the trivial fixed point x α→β (τ β ) ≡ 1.

Since it is equivalent in terms of convergence to the original choice of ψ αβ and φ α , this can be seen as the canonical choice of functions to define our Ising model.

The freedom we have in the definition of φ and ψ yields the following possibility:

Proposition 4.3. Assume that the schema (3.1)-(3.2) admits a set {m i }, i ∈ I, of fixed points with corresponding beliefs {b i }. For any i 0 ∈ I, choosing i 0 as a reference state by changing φ and ψ according to

ψ i 0 αβ (τ α , τ β ) = b i 0 αβ (τ α , τ β ) b i 0 α (τ α )b i 0 β (τ β ) , φ i 0 α (τ α ) = b i 0 α (τ α ),

INRIA

yields a new BP scheme, with unchanged beliefs {b i }, but with a new set of fixed points

m (i/i 0 ) α→β (τ β ) = m (i) α→β (τ β ) m (i 0 ) α→β (τ β ) .
In particular, the new reference fixed point {m (i 0 /i 0 ) } has all its components identically equal to 1.

Normalization and fixed points

We discuss here a feature of the algorithm which did not get that much attention in the literature, which is the possibility of normalizing the messages and its consequences on the results. In most studies, it is assumed that the messages are normalized so that (3.3) holds. The update rule (3.1) indeed indicates that there is an important risk to see the messages converge to 0 or diverge to infinity. It is however not immediate to check that the normalized version of the algorithm has the same fixed points as the original one (and therefore the Bethe approximation).

In order to make the definition of normalization clear, define the mapping

Θ αβ (m)(τ β ) def = τα∈{0,1} γ∈V(α)\{β} m γ→α (τ α ) φ α (τ α )ψ αβ (τ α , τ β ),
Then the normalized version of BP is defined by the following update rule

mα→β (τ β ) ← Θ αβ ( m)(τ β ) Θ αβ ( m)(0) + Θ αβ ( m)(1) . (4.7)
The relation between the fixed points of BP and normalized BP can be described as follows.

Proposition 4.4. Any normalized fixed point (except 0) of the BP algorithm is a fixed point of the version of BP algorithm with normalized messages.

Conversely, a fixed point of the BP algorithm with normalized messages corresponds (through multiplication by a proper constant) to an unique fixed point of the basic BP algorithm, except possibly when the graph G has exactly one cycle.

Proof. Let m be a non-null fixed point of the BP algorithm, that is

m α→β (τ β ) = Θ αβ (m)(τ β ), ∀(α, β) ∈ E and let mα→β (τ β ) = m α→β (τ β ) m α→β (0) + m α→β (1) = k αβ m α→β (τ β ).
From its definition, Θ αβ is multilinear and

Θ αβ ( m)(τ β ) = γ∈V(α)\{β} k γα Θ αβ (m)(τ β ),
and therefore m is a fixed point of the schema (4.7).

Conversely, let m be a fixed point of (4.7). Then there exists a set of constants

K αβ satisfying Θ αβ ( m)(τ β ) = K αβ mα→β (τ β ).
Let us find a set of constants c αβ such that

m α→β (τ β ) = c αβ mα→β (τ β ),
be a non-zero fixed point of (3.1). We have

Θ αβ (m)(τ β ) = γ∈V(α)\{β} c γα Θ αβ ( m)(τ β ) = γ∈V(α)\{β} c γα K αβ mα→β (τ β ) = 1 c αβ γ∈V(α)\{β} c γα K αβ m α→β (τ β ),
and therefore log c αβ -

γ∈V(α)\{β} log c γα = log K αβ . (4.8)
Solving this equation amounts to invert a matrix I -A where A is an incidence matrix on the dual factor graph (the graph which connects oriented pairs (α, β) ∈ E, see Figure 5.1). Let v αβ = log c αβ . The homogeneous equation rewrites

v αβ + v βα = γ∈V(α) v γα . (4.9)
INRIA When a non-zero solution exists, then a simple symmetry argument shows that the right-hand side does not depend on either α or β and therefore can be set to 1 without loss of generality. Therefore, summing over all oriented edges, 2|E| = 2

(α,β)∈E (v αβ + v βα ) = α∈V β∈V(α) v αβ + β∈V α∈V(β) v βα = 2|V|,
with |E| and |V| respectively the number of edges and vertices. Since G has only one component, by the well-known formula [START_REF] Berge | Théorie des graphes et ses applications[END_REF], the number of cycles in the graph is |E| -|V| + 1, only graphs with one cycle give possibly rise to a non-zero solution to (4.9). Conversely, when a graph has one unique cycle, it is possible to provide an partial ordering of vertices such that each vertex has exactly one neighbor greater than itself, and v αβ = 1 1 {α>β} is a solution to (4.9). This proposition does not describe what happens when G has exactly one cycle. The existence of a solution to (4.8) actually depends on the value of log K, which itself depends on the fixed point m. However, since BP is known to converge in a finite number of steps for graphs with at most 1 cycle, normalization is not useful in this situation.

From now on reference to the BP algorithm is to be understood as its normalized version.

Stability of BP fixed points

The next issue to tackle regarding the fixed points of BP is their stability. The following definition of conditional belief will be useful

b αβ (τ α |τ β ) def = b αβ (τ α , τ β ) b β (τ β ) .
For the general case we have the following Proposition 4.5. The stability of any fixed point of the BP algorithm is determined by the set of beliefs {b} of that fixed point: the fixed point is stable if, and only if, the matrix defined, for any pair of oriented edges (α, β) ∈ E, (α ′ , β ′ ) ∈ E, by the elements

J α ′ β ′ αβ = b αβ (1|1) -b αβ (1|0) 1 1 {α ′ ∈V(α)\{β}, β ′ =α} = 1 -b αβ (0|1) -b αβ (1|0) 1 1 {α ′ ∈V(α)\{β}, β ′ =α} , (4.10)
has a spectral radius smaller than 1.

A sufficient condition for this stability is therefore

b αβ (1|1) -b αβ (1|0) < 1 q α -1
, for all α ∈ V, β ∈ V(α).

Proof. Since we are dealing with binary variables, messages are vectors with two components, and it is easier to set

η α→β def = m α→β (1) m α→β (0) .
This normalization is equivalent to the one proposed in Section 4. after performing the change of referential of Proposition 4.3 with reference point {b}. We look for small perturbations around the fixed point η α→β ≡ 1 for all (α, β). The Jacobian at the point η = 1 reads:

∂η αβ ∂η α ′ β ′ η=1 = b αβ (1|1) -b αβ (1|0) 1 1 {α ′ ∈V(α)\{β}, β ′ =α} ,
which proves (4.10). The rest of the proposition is a consequence of basic inequalities on the spectral radius of a matrix.

INRIA 0 0.5 1 1.5 2 0 0.5 1 1.5 2 f(η) η a) (q -1)κ < -1 0 0.5 1 1.5 2 0 0.5 1 1.5 2 f(η) η b) (q -1)κ ∈ (-1, 0
) and (q -1)κ ∈ (0, 1) Remark For a totally symmetric graph with connectivity q, (4.11) reduces to If one considers the dual graph formed by function nodes, where links relate pairs of function nodes having a variable in common, then on this graph the Jacobian matrix has the structure of the incidence matrix A already encountered in the preceding section. This matrix is not symmetric, but eigenvalues greater than 1 in modulus indicate anyway an instability. These are obtained by forming the new matrix J(λ) def = J -λI with I the identity matrix and finding roots of det J(λ) = 0.

0 0.5 1 1.5 2 0 0.5 1 1.5 2 f(η) η c) (q -1)κ > 1
η ← f (η) def = b(0|1) + η q-1 b(1|1) b(0|0) + η q-1 b(1|0) ,
The expansion of det J(λ) involves permutations which are compatible with circuits of the dual graph, where each vertex is visited once. Each permutation is uniquely represented by a product of permutation cycles (orbits) with disjoint support and is attached to a sub-graph of the dual graph. Let us call maximal permutation, a permutation such that the complementary graph of its associated sub-graph is cycle free. Adapting results of [START_REF] Harary | The determinant of the adjacency matrix of a graph[END_REF], det J(λ) may be expanded according to the following, Proposition 4.6.

det

J(λ) = σ∈C ω i ⊂σ (det ω i + (-λ) |ω i | ), (4.12)
where the sum runs over all possible maximal permutations σ, each one being expressed as a product of n ≥ 1 circular permutations (cycles) ω i , i = 1 . . . n, of size |ω i |, with determinant given by

det ω = -(-1) |ω| (α,β)∈ω b αβ (1|1) -b αβ (1|0) .
On a tree, as expected, zero is the only eigenvalue, in fact J is a nilpotent matrix of index the size of the longest directed path in the graph. If there is only one cycle ω, (4.12) reduces to det J(λ) = λ N -|ω| (det ω + (-λ) |ω| ), which yields the eigenvalues

λ k = (α,β)∈ω b αβ (1|0) -b αβ (1|1)
1 ω e (2k+1)iπ/|ω| , with modulus obviously smaller than one. As a consequence, the following proposition holds. Proposition 4.7. BP fixed points for a graph containing at most one oriented loop are stable.

This has been remarked by different means in [START_REF] Heskes | On the uniqueness of loopy belief propagation fixed points[END_REF]. Unstable modes correspond to eigenvalues larger than 1, and might reveal vertices or cycles mostly responsible for the instabilities. An interesting case occurs when cycles of the dual graph have disjoint supports, because then only one maximal permutation σ exists and expansion (4.12) reduces to one term,

det J(λ) = ω i ⊂σ (det ω i + (-λ) |ω i | ).
As a result, since the modulus of the Jacobian coefficients are always smaller than 1, to each cycle is associated an eigenvalue smaller than 1 and the state is stable.

On a graph which is locally a tree (Bethe lattice), a mean-field equation can be used to evaluate the stability of a given fixed point. The idea is to consider the iterated Jacobian matrix in a statistical manner, by looking at the distribution P (n) (v) of components v of an iterated vector starting from a non-degenerate initial condition,

V (n) = J n V (0) .

INRIA

The mean-field stability equation then simply reads (after assuming the usual independence property of parent messages)

P n+1 (v) = c>1 Q(c) {v i },{κ i } c-1 i=1 P n (v i )R(κ i )δ(v - c-1 i=1 v i κ i ), (4.13)
with Q the connectivity distribution in the dual graph and R the Jacobian coefficient distribution (see Figure 5.3 for example). The instability is therefore fully characterized by the statistical properties of the considered BP fixed point and by the statistical properties of the graph (connectivity), which sometimes can be an adjustable parameter.

5 Toy Model simulations

From theory to practice

We illustrate these ideas on a simulated traffic system which has the advantage to yield exact empirical data correlations. For real data, problems may arise because of noise in the historical information used to build the model. This additional difficulty will be treated in a separate work.

The model consists of a queueing network system. Each queue represents a link of the traffic network (a single-way lane) and has a finite capacity; to each link we attach a variable ρ ∈ [0, 1], the car density, which is represented by a color code in the picture (Figure 5.2 on page 23).

As already stated in Section 2, the physical traffic network is replicated, to form a space time graph, in which each vertex α = (ℓ, t) corresponds to a link ℓ at a given time t of the traffic graph. To any space-time vertex α, we associate a binary congestion variable τ α ∈ {0, 1}.

The statistical physics description amounts to relating the probability of saturation P (τ α = 1) to the density ρ α . For the sake of simplicity, we consider a linear relation and build our historical p according to the rules

α β γ 1 γ 2 γ 3 γ 4 γ3 → α γ2 → α γ1 → α α → β γ4 → α
pα (1) = µ ∞ (ρ ℓ ), pα (0) = µ ∞ (1 -ρ ℓ ), pαβ (1, 1) = µ ∞ (ρ ℓ ρ ℓ ′ ), pαβ (1, 0) = µ ∞ (ρ ℓ (1 -ρ ℓ ′ )), . . .
where µ ∞ is simply a frequency estimator. Note that, to follow some realistic statistical constraints, we use here only data aggregated in time. More realistic data collection and modeling would work the same way. The structure of the factor-graph on which we propagate the information is depicted in Figure 5.1. Some fine tuning is required to let the algorithm work correctly. First, from Proposition 4.5, we know that the stability of the reference point p encoded in (3.6)-(3.7) is not guaranteed; this may be evaluated on the basis of distributions depicted in Figure 5.3, using equation (4.13). In absence of negative correlations, it is likely that our system is either a paramagnetic-like (in the Ising-model terminology) system, with small fluctuations around average values, or a ferromagnetic-like system, in the sense that positive correlation drive the system to a state where links are in a similar state, i.e. mostly fluid (low state) or congested (high state). This which means that eigenvalues are rescaled by a factor ǫ. For our purpose, this provides us with an adjustable mean-field parameter, to correct some artificial amplification of correlations caused by closed loops in the graph. We expect that there exists a critical value of ǫ c corresponding to the ferromagnetic phase transition point (high temperature means here small ǫ). In addition, since for small ǫ we recover in one sweep the bare mean results, this parameter can be used for a simulated annealing procedure, by letting it converge from zero to the desired value during the BP iterations.

The second adjustment concerns the encoding of real-time information. The probe vehicle is assumed to send an information for some space-time vertex α * , typically in the form of an instantaneous velocity, from which is estimated the probability p α * of saturation. Instead of projecting this information on one of the two states (τ α * = 0 or τ α * = 1), which turns out in practice to be too coarse, we use a procedure which amounts to bias the messages sent by α * in proportion to the observed belief p α * . In the statistical physics language, this amounts to impose an external local field on the observed variable.

The last issue concerns the situation where the system is below the transition point, in which case we have two separate states, and it is always possible that BP converges towards the wrong fixed-point. In this simple ferromagnetic situation, it is in fact easy to enforce the convergence of the algorithm to a specified fixed point by applying a slowly decaying external field, enforcing either the fluid or the congested state. As a result of this procedure, we obtain two sets of beliefs {b 0 } and {b 1 }, with corresponding free energies F 0 and F 1 , from which we build the superposition belief,

b α = e -F 0 b 0 α + e -F 1 b 1 α e -F 0 + e -F 1 .
RR n 6144 which in practice, for sufficiently large systems, because F is extensive, turns out to be the set of beliefs corresponding to the lowest free energy.

In the following we refer the sets of belief {b 0 α }, {b 1 α } and {b α } respectively to the low, the high and the combined inference state. Accordingly the set {p α } is referred to as the historical state. In addition, the combination of observations with historical data (by replacing the historical value with the last observation in the window time) yields the actual state. To estimate the quality of the traffic restoration we use the following estimator:

reconstruction rate def = 1 |V| α∈V 1 1 {|bα-ρα|<0.2} ,
which computes the fraction of space-time nodes α for which the belief b α does not differ by more than an arbitrary threshold of 0.2 from ρ α .

Numerical results

We have tested the algorithm on the toy traffic network shown on the program's screen-shot of Figure 5.2. The characteristics of this network are summarized in Table 5.1. Two types of traffic conditions have been used, that both correspond to periodic oscillation superimposed with noise (see blue curve of Figure 5.4); they simply differ by the level of the noise.

The two values of ǫ c in Table 5.1 that have been computed for the two different traffic regimes using (4.13) are close to the observed values, which indicates that the space-time graph on which BP is run is close to the conditions of a dilute graph (Bethe lattice).

The simulation run of Figure 5.4 compares the policies of using only the low state (green), the high state (red) or the combined state w.r.t. the free energy in the low-noise case. Abrupt changes of the combined state prediction correspond to the crossing of the Bethe free energies. In the transition regimes, which correspond to out-of-equilibrium situations, the free energy criteria sometimes select the wrong state. The reason for this is that the present design of our algorithm encodes only statistical information at equilibrium. Time correlations should be incorporated in some way, to encode transition rates between the macro-states (here the low and high traffic density).

Distributions of performance errors shown in Figure 5.5 are based on a simulation run of 10000 traffic time units where a belief propagation is run every 3 units of time for both low and high inference states, to reconstruct the traffic. Varying the

Conclusion and perspectives

We have presented a novel methodology for reconstruction and prediction of traffic using the Belief Propagation algorithm on Floating Car Data. We have shown how the underlying Ising model can be determined in a straightforward manner and that it is unique up to some change of variables. In addition, the effect of message normalization and the stability properties can be asserted from the original data. The unfortunate fact that the BP fixed point corresponding to the historical data may be unstable can be circumvented by rescaling of the correlations. The algorithm has been implemented and illustrated using a toy traffic model. Several generalizations are considered for future work:

• firstly, the binary description corresponding to the underlying Ising model is arbitrary. Traffic patterns could be represented in terms of p different states. A Potts model with p-states variables would leave the belief propagation algorithm and its stability properties structurally unchanged. Actually this number p should be subject to an optimization procedure.

• secondly, our way of encoding traffic network information might need to be augmented to cope with real world situations. This would simply amount to redefine the factor-graph used to propagate this information. In particular it is likely that a great deal of information is contained in the correlations of local congestion with aggregate traffic indexes, corresponding to sub-regions of the traffic network. Taking these correlations into account would result in the introduction of specific variables and function nodes associated to these aggregate traffic indexes. These aggregate variables would naturally lead to a hierarchical representation of the factor graph, which is necessary for inferring the traffic on large scale network. Additionally, time dependent correlations which are needed for the description of traffic, which by essence is an out of equilibrium phenomenon, could be conveniently encoded in these traffic index variables.

Ultimately, for the elaboration of a powerful prediction system, the structure of the information content of a traffic-road network has to be elucidated through a RR n 6144

Fig. 2 . 1 :

 21 Fig. 2.1: Traffic network (a) and Ising model (b) on a random graph

Proposition 4 . 1 .

 41 A set of beliefs {b α , b αβ } corresponding to a BP fixed point of (3.1)-(3.5) always satisfies

Proposition 4 . 2 .

 42 Assume that there exists a fixed point of the BP algorithm satisfying (3.1)-(3.5) and such that b αβ (τ α , τ β ) = pαβ (τ α , τ β ), and therefore b α (τ α ) = pα (τ α ). (4.1)

  α→β , and the scaled BP algorithm update rule (4.7) can be rewritten as η α→β ← b αβ (0|1) + γ∈V(α)\{β} η γ→α b αβ (1|1) b αβ (0|0) + γ∈V(α)\{β} η γ→α b αβ (1|0) , (4.11)

Fig. 4 . 1 :

 41 Fig. 4.1: Different possible graphs of f (η) versus η depending on the value of κ = b(1|1) -b(1|0): (a) one unstable fixed point, (b) one stable fixed point and (c) one unstable and two stable fixed points.

  and the classification in terms of b(1|1) -b(1|0) is pictured in Figure 4.1. Note that b(1|1)-b(1|0) > 0 (resp. < 0) corresponds a ferromagnetic (resp. anti-ferromagnetic) system.

Fig. 5 . 1 :

 51 Fig. 5.1: Structure of the factor graph (left), 3 time-layers portions are represented, black circles correspond to crossroads and blue squares to factor-vertices. Corresponding graph for the Jacobian matrix (right).

INRIA

  scenario corresponds to the regimes pictured on Figure4.1 where case (c) is the usual ferromagnetic phase transition in the Ising model. It is also a well-known fact that this transition is driven by the temperature. To introduce the equivalent of a temperature in our equations, since its effect is essentially to reduce correlations, let us consider modified pairwise marginal lawspαβ (ǫ) = ǫp αβ + (1 -ǫ)p α pβ .The high temperature regime corresponds here to ǫ → 0 and the vanishing of the correlations. The Jacobian coefficients κ αβ def = b αβ (1|1)-b αβ (1|0) are modified according to καβ (ǫ) = ǫκ αβ ,

Fig. 5 . 2 :

 52 Fig. 5.2: Traffic network as produced by the simulator. The continuous color code represents the traffic index from 0 (full green) to 1 (full red).

  Fig. 5.3: Connectivity (black), correlation coefficients (red) and Jacobian coefficients (blue) histograms for low (top) and high (bottom) noise level.

Fig. 5 . 6 :Fig. 5 . 7 :

 5657 Fig. 5.6: (a) Reconstruction rates obtained with 0, 1, 5 and 10 probe vehicles and various possible inference states; (b) variance of the reconstruction rate obtained with 10 vehicles also for the various possible states; (c) reconstruction rates obtained for the noisy network, again with 10 probes.

  Reconstruction rates for the various possible inference states as a function of time with corresponding free energies, with 10 probe vehicles and ǫ = 0.75.

				1
		1					combined state 1,5,10 probes low state 1,5,10 probes
				0,9			high state 1,5,10 probes
		0,8					historical state
				0,8
	prediction rate/ traffic index/ energy	-0,4 -0,2 0 0,2 0,4 0,6	prediction rate	0,4 0,5 0,6 0,7 0,3	low Bethe free energy combined state prediction low state prediction high state prediction traffic level 0,6	epsilon	0,8	1	(a)
					high Bethe free energy
					3000		4000		5000
						time
	3: Connectivity (black), correlation coefficients (red) and Jacobian coefficients (blue) histograms for low (top) and high (bottom) noise level. 35 122 43 5246 0.67 Fig. 5.5: Distribution of reconstruction errors for the various possible inference epsilon 1.29 reconstruction error 0,6 0,8 1 1,2 0,4 nodes links time steps graph size ǫ c 1 (oscillating) ǫ c 2 (noisy) 0,2 0,4 0,6 0,8 1 0 0,15 0,2 0,6 0,8 1 epsilon 0 0,1 0,2 prediction rate variance combined state, 10 probes low state, 10 probes high state, 10 probes historical state (b) 0,5 0,8 Fig. 5.4: 0 0,05 0,1 actual state, 10 probes 0,6 prediction rate historical state distribution combined state, 10 probes combined state, 10 probes low state, 10 probes low state, 10 probes high state, 10 probes high state, 10 probes 0,7 historical state
	Tab. 5.1: Toy model characteristics states for ǫ = 0.75.
									INRIA

RR n 6144

β∈V(α)

specific statistical analysis. The use of probe vehicles, based on modern communications devices, combined with a belief propagation approach, is in this respect a very promising approach.