Jan Sanromà Sánchez

Andreas Gehrmann

Dr Paul Thevenon

Philippe Brocard

Ben Amani

Olivier Afia

Julien

Dr Olivier Julien

Use of a FishEye Camera for GNSS NLOS Exclusion and Characterization in Urban Environments

His research focused is on navigation in urban environments and the use of fisheye camera to increase positioning accuracy.

INTRODUCTION

Urban environments present a significant challenge for modern Global Navigation Satellite Systems (GNSS)based systems. These systems were originally designed to work in open outdoor areas, where direct line-of-sight view of the satellite is always possible. In cities however, surrounding buildings block the view of the satellite and their coming signals. The signals from these satellites can however arrive at the receiver trough different reflections, or multipath, still with strength. These reflections cause an increase in the travel distance of the signal, and thus these NLOS signals create large errors in the pseudorange measurements.

A usual solution to this problem is to mix GNSS information together with INS. INS estimate the user's dynamics and provide a position solution without using any external agent, and thus they do not suffer from problems when being used in urban environments. By mixing them with GNSS it is possible to produce a higher accuracy solution, reducing the importance of NLOS and achieving a better navigation solution. The problem of multipath and NLOS is however still present. What is proposed in this paper is to use a wide field of view (FOV) camera mounted on the roof of a car and looking upward. This camera will take a snapshot synchronized with GNSS measurements. Then, by projecting the satellites on that image, it is intended to determine which satellites are in the section of visible sky, and which are hidden behind buildings and sending NLOS signals. All the sensors used will represent cheap, mass market devices, in order to show the performance of a generic configuration. The main objective is thus to develop an algorithm that is capable to analyze each picture and detect which satellites are NLOS. Those satellites can then be removed from the computation, ideally eliminating the effect of multipath. By combining INS and camera exclusion with the GNSS measurements, a higher accuracy and more robust solution should be achieved.

GNSS/INS INTEGRATION

In this section the integrated INS/GNSS architecture is presented. Sensor fusion will allow combining the benefits of both GNSS and INS, and achieve a high level of accuracy and availability in urban environments. The integration algorithm uses a tightly coupled integration scheme based on a Kalman filter as presented in [START_REF] Groves | Principles of GNSS, inertial, and multisensor integrated navigation systems[END_REF]. The general strategy of a tightly-coupled architecture is to use GNSS measurements (pseudorange and Doppler) in order to produce an estimate of the inertial navigation system errors. These errors are mainly due to the inertial measurement unit (IMU) biases. Once estimated, these biases can be feedback to the IMU in a closed-loop correction. Figure 1 shows the general scheme of the proposed integration algorithm. In this architecture, the inertial navigation system estimates the vehicle position, velocity and attitude using the IMU measurements. This navigation solution is input to the Kalman filter which estimates, in addition to the IMU biases, the errors of this inertial solution using GNSS pseudoranges and Doppler measurements and outputs the corrected navigation solution.

GNSS measurement models

Pseudoranges and Doppler measurements were used in the Kalman filter. Both types must thus be correctly modeled in order to properly tune the Kalman filter. The pseudorange measurement is an estimate of the distance between the GNSS satellite and the receiver antenna, but also includes the mis-synchronization of satellite and receiver local oscillators. This measurement is computed from the propagation time between the transmission of the signal by the satellite and its reception by the receiver, and the multiplication of this time by the speed of light c. The fact that the clocks of the satellites and the receiver are not perfectly synchronized to GPS time, introduces a clock bias that drifts with time. Therefore, the pseudorange (PR) measurement can be modeled as:

𝜌 ̃𝐺 𝑖 = 𝑟 𝑖 + 𝑐𝛿𝑡 + 𝜂 𝜌𝐺 𝑖 (2) where  𝑑 ̃𝐺 𝑖 is the Doppler measurement of satellite i,  𝑟̇𝑖 is the true pseudorange rate between the satellite and the receiver,  𝛿 ̇𝑡 is the clock drift and  𝜂 𝜌̇𝐺 𝑖 is the range rate error modeled as a zero-mean Gaussian white noise with a variance denoted by 𝜎 𝑃𝑅𝑅 2 .

In cities, the biggest source of measurement errors will be the local effects affecting the tracking loops, on top of which will be multipath-and NLOS-induced errors. In order to try to correctly model this type of errors for both the pseudorange and pseudorange rate, the models proposed in [START_REF] Carcanague | Low-cost GPS/GLONASS Precise Positioning Algorithm in Constrainted Environments[END_REF] and [START_REF] Kuusniemi | User level reliability and quality monitoring in satellite based personal[END_REF], based on data collection in urban environments and commercial mass market receivers, were used. These models were specially proposed for use in urban environments, under conditions of heavy multipath.

For the pseudorange measurements, the covariance of the measurements' errors is defined as a function of the 𝐶/𝑁 0 as:

𝜎 𝑃𝑅 2 = 𝑎 + 𝑏10 -𝐶/𝑁 0 10 (3)
Parameters 𝑎 and 𝑏 of (3) are defined as in Table 1. For the Doppler measurements, this variance is defined as a function of the user velocity and received 𝐶/𝑁 0 according to "the Doppler measurement standard deviation as a function of velocity and 𝐶/𝑁 0 " table defined in [START_REF] Carcanague | Low-cost GPS/GLONASS Precise Positioning Algorithm in Constrainted Environments[END_REF].

Inertial Measurements Models

When fusing the data from GNSS and an INS, the position, velocity and attitude are predicted with the INS. Therefore, the Kalman filter process model shall characterize the evolution of the error on the states as propagated with the INS. These error models can be derived from the inertial sensors error models.

GNSS/INS Integrated Kalman Filter

Based on the integration scheme and the sensors' measurement models previously described, the Kalman state vector can be written in the Earth-Centered-Earth-Fixed (ECEF) frame as:

𝑿 𝑒 = [𝛿𝝍 𝑏𝑒 𝛿𝒗 𝑒 𝛿𝒑 𝑒 𝒃 𝑎 𝒃 𝑔 𝑐𝛿𝑡 𝑐𝛿 ̇𝑡] 𝑇 (6)
where  𝛿𝝍 𝑏𝑒 is the inertial attitude (orientation of the vehicle body frame (b) with respect to the ECEF (e) frame) error  𝛿𝒗 𝑒 is the vehicle velocity error expressed in the (e) frame  𝛿𝒑 𝑒 is the vehicle position error in the (e) frame  𝒃 𝑎 and 𝒃 𝑔 are respectively the IMU accelerometer and gyroscope biases  𝑐𝛿𝑡 and 𝑐𝛿 ̇𝑡 are respectively the GNSS clock bias and clock drift. If multiple GNSSs are used, then the time difference between the two systems should be also estimated in the state-vector. However, the clock drift between the two systems could be neglected, since GNSS time is very stable.

The Kalman filter propagation step is based on the computation of the navigation solution using the inertial system. The state transition model of the Kalman filter is derived by computing the state-vector rate of change and expressing it as a function of the state-vector and the system noise generated by the IMU and GNSS clock noises. This state transition equation is written as: 𝑿 ̇𝑒 = 𝑭 𝑒 𝑿 𝑒 + 𝑮 𝑒 𝒖 (7) where  𝑭 𝒆 is the state transition matrix,  𝑮 𝑒 is the system noise distribution matrix and  𝒖 is the covariance of the system noise vector The expressions of all these matrices and vectors are given in [START_REF] Groves | Principles of GNSS, inertial, and multisensor integrated navigation systems[END_REF].

The Kalman filter update uses the GNSS pseudoranges and pseudorange rate measurements in order to correct the inertial system outputs as well as the IMU biases. The observation model relates the GNSS measurement innovation 𝛿𝒛 𝐺 to the state-vector 𝑿 𝑒 as follows

𝛿𝒛 𝐺 = [ 𝜌 ̃𝐺 1 -𝜌 ̂𝐺 1 ⋮ 𝜌 ̃𝐺 𝑁 -𝜌 ̂𝐺 𝑁 𝜌̇𝐺 1 -𝜌̇𝐺 1 ⋮ 𝜌̇𝐺 𝑁 -𝜌̇𝐺 𝑁 ] = 𝑯 𝑒 𝑿 𝑒 + [ 𝜂 𝜌𝐺 1 ⋮ 𝜂 𝜌𝐺 𝑁 𝜂 𝜌̇𝐺 1 ⋮ 𝜂 𝜌̇𝐺 𝑁 ] (8)
where  𝜌 ̃𝐺 𝑖 is the i th satellite measured pseudorange  𝜌 ̂𝐺 𝑖 is the i th satellite pseudorange predicted from the inertial navigation solution  𝜌̇𝐺 𝑖 is the i th satellite measured pseudorange rate  𝜌̇𝐺 𝑖 is the i th satellite pseudorange rate predicted from the inertial navigation solution  𝑯 𝑒 is the measurement matrix. Its expression can be found in [START_REF] Groves | Principles of GNSS, inertial, and multisensor integrated navigation systems[END_REF] NLOS DETECTION ALGORITHM By joining INS and GNSS information in a tight architecture, the inertial sensors would guide the solution on those situations where very few GNSS satellites are in view, while GNSS satellites correct the drift present in the INS solution.

While navigating through a city, however, GNSS measurements will be heavily affected by multipath and NLOS satellites. This multipath can cause large errors on the GNSS pseudoranges and Doppler measurements. If these faulty GNSS measurements are not good enough to be used to correct for the INS drift, over time the position accuracy worsens considerably, as errors in the IMU start accumulating in the solution.

As explained before, the worse impact of multipath happens when direct view of the transmitting satellite is blocked by a solid object, like a building. In this situation, the receiver could start tracking a reflected signal coming from the same satellite as the original one. It is with the objective to block multipath, that this paper proposes the use of a high field-of-view (185 °) fish-eye camera. This camera takes pictures of the sky above (Figure 2), synchronized with the GNSS measurements. It is then possible to compute the satellite elevation and azimuth (which comes from the GNSS ephemeris and approximate receiver location) to project those satellites on the image, as long as the attitude of the camera, provided by the INS, is well estimated. The satellites that are placed in the sky part of the image (in green) can thus be treated as LOS with confidence. Those satellites placed behind a building or any other obstructing object will be considered as NLOS, or at least as less trustworthy signals. Once the distinction between LOS and NLOS (or suspicious) signals is made, it is possible to exclude or down-weight the corresponding pseudorange and Doppler measurements in the Kalman filter. This would ideally allow obtaining a better position and attitude solution.

Image Processing Algorithm Description

The objective of the image processing is to tell apart the sky from the surrounding masking buildings in each image.

The following operations are applied to an image to perform the identification of the sky vs blocked areas. Mistakes in the image processing can cause NLOS satellites to be used, or LOS satellites to be discarded. Different environments and weather conditions will present different challenges for the image processing and the camera, as they have a great effect on the way the image is captured and its quality. Some of them are:  Sun glares, clouds and trees. These three factors are very limiting in image segmentation as they appear as a solid object in an image. As they have edges in the picture, they limit sky detection, even though they have no effect on pseudorange measurements.  Rain: Although not tested, as the camera is not insulated against water, rain would cause a huge problem as water droplets would block the image. Some hydrophobic spray could be used, and then then the images would only show a cloudy day.  Night and tunnels:

Image processing is impossible during night time, as the image cannot properly detect the edge between a building and the sky. Long bridges and tunnels would present a similar problem, with the added problem that fast changes in brightness (like lamps) would blind the camera for a few seconds.

Satellite Selection based on 𝑪/𝑵 𝟎

In order to compensate for these problems, the proposed NLOS detection algorithm also uses 𝐶/𝑁 0 measurements. It is known that high elevation satellites should have a high 𝐶/𝑁 0 , while NLOS satellites should have a lower received 𝐶/𝑁 0 , although NLOS with high 𝐶/𝑁 0 can occurr. This information can then be used to make a first estimation of where the visible sky area is in the picture. The overview of the algorithm is as follows: them on the image based on the estimated attitude of the camera given by the IMU. 4) Choose the satellites with a 𝐶/𝑁 0 above 45 𝑑𝐵𝐻𝑧 (marked in green). Do a flood fill starting from these satellites. This operation spreads color until a border is found. Sky is then detected as the black region, and any other satellite present in it will be detected as LOS.

Once the sky area has been detected, it is possible to decide whether the rest of the satellites are LOS or NLOS based on the picture information: if they fall in the dark area (representing the sky) they are considered as LOS, if not they are discarded. In order to avoid the biggest errors, and based on experience, satellites with 𝐶/𝑁 0 above 40 𝑑𝐵𝐻𝑧 will be automatically considered as LOS, while those with power levels below 30 𝑑𝐵𝐻𝑧 will be considered as NLOS.

The different thresholds used can be seen in Figure 4. This is similar to what is proposed in [START_REF] Suzuki | N-LOS GNSS Signal Detection Using Fish-eye camera for vehicle navigation in urban environments[END_REF].

Figure 4: Thresholds for the camera algorithm.

The result of this algorithm is that the camera can be used in the variety of environments described before.  Clouds or sun glares do not affect the output even though they are detected as edges in the image.  Trees have a higher effect, as they reduce the 𝐶/𝑁 0 measurements a bit, but still they would not fully block the detection.  In places where no satellites are above the sky detection or LOS thresholds, all satellites are excluded, trusting only the inertial output. This allows the method to be used in tunnels, as any GNSS value received (if any) will be erroneous.  During night time the 𝐶/𝑁 0 will be the same, but the camera information will be useless (another type of camera can be used).

As per processing time, this algorithm runs close to realtime, even though it has only been tested in post processing. This allows fast testing for each of the thresholds and the algorithm accuracy.

TESTS Test campaign and conditions

Two test campaigns were performed on the same day to show the performance of the NLOS exclusion. One campaign consisted of a mix of suburban and urban environments and lasted for around 40 minutes. The second test included a highway section and lasted 30 minutes. The routes are shown in Figure 5. Toulouse city center consists of very narrow streets with very short houses, just three or four story high. It is worth mentioning also that the day was very sunny and the test was done soon after midday in the summer. It was seen that the sun would blind the camera, making segmentation much more difficult. A significant sun glare can also be seen in the center of most of the images.

Equipment Set Up

The equipment used consists of cheap to mid-cost sensors with normal mass-market performance:  U-Blox M8T: GPS/GLONASS receiver. Also used for timing the IMU and camera measurements.  Xsens MTi IMU  IDS uEye black and white camera with a fish-eye lens.  Novatel Span: reference equipment used to evaluate accuracy. [START_REF]Novatel[END_REF] Figure 6: Equipment mounted inside (Left) and on top (right) of a van.

All this equipment was mounted on top of a van for the data collection. The u-Blox receiver being a dual constellation mono-frequency GPS/GLONASS receiver, the increase in satellite availability caused by the dual constellation showed that at least one or two satellites were always in direct LOS view, even in the narrowest streets.

The receiver was also used as a timing reference for the IMU and the camera, which do not include accurate clocks of their own. Pictures and GNSS measurements were taken at 1 𝐻𝑧, while the IMU has an output rate of 100 𝐻𝑧.

The reference equipment produced centimeter level positioning and it was used to estimate the positioning accuracy of our algorithm. The heading estimation used to project the satellites on top of the camera image was initially based on the reference equipment, in order to first only test the selection algorithm based on the correct projection of the satellites on the image.

Initial Considerations on Code and Doppler Measurements' Filtering

The Kalman Filter innovation gives an indication on the discrepancy between reconstructed GNSS measurements based on the propagated solution and the GNSS measurements. If a large (abnormal) innovation is seen, it can be assumed that it comes from a faulty GNSS measurement, caused by a NLOS satellite. Innovation filtering then consists in setting a threshold, over which one of the measurements is considered as faulty. Table 3 shows the accuracy obtained when adding Doppler innovation filtering with 2 thresholds 4𝑚/𝑠 and 1𝑚/𝑠 (user velocity).

In this initial test, the camera was not used. It can be seen that Doppler innovation filtering brings a significant improvement in the position accuracy. These results indicate us that Doppler measurements can be extremely untrustworthy in urban environments, although when correctly filtered, they improve the accuracy. As a consequence, in the forthcoming results, a specific Doppler innovation filtering test is used with threshold of 4𝑚/𝑠. Besides, initial test using the camera exclusion detection were performed and a parametric analysis was made. In particular, a trade-off between a full exclusion and a simple down-weighting of the NLOS detected measurements was performed. It was shown that an optimal solution (with respect to the test campaign) was to fully exclude all the Doppler measurements from NLOS detected measurements and to only down-weight, by a factor of 1.5 on the variance values, the NLOS-detected pseudorange measurements. This is thus used in the following.

Urban Test

During the urban test, the influence of the camera can be clearly seen. The first 10 minutes of data are located in a suburban area between ENAC and Toulouse center, while the rest was located in the narrow streets of the city center. Table 4 shows the accuracy results for three cases:

 The first one is the case when only the INS/GNSS algorithm is used (no camera exclusion).  The second column represents the accuracy results when we turn on the camera for NLOS exclusion, and we use the extremely accurate reference heading to project the satellites on the images. This is obviously outside the scope of the system, as it relies on the very expensive reference receiver. Nonetheless it offers a good approximation on how good the NLOS exclusion algorithm is, as there are no errors when projecting the satellites on the image.  The final column represents what happens when, instead of using the reference heading for projecting the satellites, we use our own estimated heading using our IMU. This would be the performance if we used our GNSS, IMU and camera in an integrated, stand-alone device. It can be seen that the NLOS exclusion when using the reference heading does reduce the horizontal error. Both the mean error, and the 95 th percentile see a reduction of around one meter when the camera is activated. Figure 7 and Figure 8 show this tendency over the duration of the test. It can be seen that use of the camera exclusion reduces much of the errors present in the INS/GNSS integration. The heading error, which suffers a lot from the lack of satellite visibility in cities, also decreases to a very usable level.

However, it is worth noting that the maximum horizontal error does increase when using the camera. This is caused by NLOS detection errors in very low visibility areas. These errors cause a decrease in horizontal accuracy for an epoch or two, but then the solution improves again. In terms of the number of satellites excluded, Figure 9 shows a huge drop in the number of satellites used when the NLOS detection algorithm is used. The mean number of satellites seen by both configurations shows a difference of more than 7 satellites in Toulouse city center (from 15 to less than 8), and about 4 to 5 in sub-urban areas, showing how strong this exclusion is in an urban environment. In the narrowest streets, it can be seen that most of the time there are less than 4 satellites available when using the camera, when more than 10 satellites are generally available.

Highway test

During highway test, satellite visibility was much higher, as there were fewer obstacles along the way. Figure 10 shows the decrease in the number of available satellites when using the camera. The graph shows that most of the time only two or three satellites are excluded. However, there are many local drops in the number of satellites. These drops have many causes: bridges and traffic signs overhead would block the view of the camera, and so very few satellites would be used. Trucks driving next to the car would also block the visibility on one side. The camera then does not add much to the solution accuracy, as seen in Table 5. However some special environments could also be observed in the test. The performance of the algorithm was tested when driving under a short bridge. Figure 11 shows that once inside a tunnel, the algorithm decides (based on C/N0 values) that the van is no longer under open sky (no black areas in the image), and so it discard all satellites (seen in red). It then goes to fully inertial navigation for the duration of the tunnel. Figure 12 also shows the performance of the sky detection under trees. Trees cause troubles for the camera exclusion, as they are detected as solid objects in the picture: they have many edges and different colors than the sky, and thus they limit the image segmentation. Since the 𝐶/𝑁 0 information is used for the image floodfill, the algorithm knows that those satellites behind the tree are in fact LOS, as their 𝐶/𝑁 0 is high. However, many satellites are still detected as NLOS, as the dense foliage does have an effect over the received power levels. One solution could be to interpret such a fragmented image in a specific way associated to the presence of a tree. Figure 13 shows how satellites above 40 dB-Hz are always taken as LOS, while satellites with received power below 30 𝑑𝐵𝐻𝑧 are taken as NLOS. A gray area exists then between 30 -40 𝑑𝐵𝐻𝑧 in which the camera performs the main part of its detection works.

CONCLUSIONS AND FUTURE WORK

This article has shown an algorithm in which INS, GNSS and fisheye camera data are integrated in order to produce accurate positioning in urban environments. The tests performed with the algorithm show that it is a viable option, and that it does increase the accuracy of the solution. This increase in accuracy happens with respect to a GNSS/INS solution which is already quite precise, as the use of double constellation already allows for very good The most important factor in these results is that the integration was proved to work as a stand-alone device.

The u-Blox receiver can be used together with the xSense IMU and the IDS fish-eye camera in order to produce a more accurate solution for urban environments.

It is also important to mention the fact that the test route was in no way thought or planned to prove any performance. The route included the normal characteristics of urban environments and the weather conditions were not especially favorable for the camera detection, as the camera would be blinded by the sun.

Future work on this subject is still necessary, as the increase of accuracy does not necessarily justify the cost of adding a camera into the integration. Three main areas of quick improvement are seen for the current status of the integration algorithm:  fine tuning of the thresholds used for the NLOS exclusion, as seen in Figure 4. The current threshold was not based on hard data of the receiver's performance. The thresholds could probably be raised in order to obtain a more robust and accurate solution.  use of a color camera instead of a black and white one. The increased quality and range of colors will make segmentation between sky and buildings much easier, reducing the errors.

Finally, it is worth mentioning that such an equipment could be used to evaluate the NLOS measurements characteristics. Indeed, most of the safety-critical applications that are using GNSS-based location devices working in urban or sub-urban environments required a knowledge on the NLOS occurrence rate and NLOS measurement error distribution. The proposed scheme is fully able to do so since it can be integrated with a reference trajectory system. Therefore, measurement errors can be accurately estimated and can be associated to LOS and/or NLOS situations. It can also be related to the satellite elevation and/or azimuth. Initiating such measurement campaigns is thus also part of the future work.

Figure 1 : 2 The

 12 Figure 1: Integration scheme between GNSS, INS and camera exclusion

Figure 2 :

 2 Figure 2: Principle of NLOS exclusion and sky detection (green area).

1 ) 2 ) 4 ) 5 )

 1245 Morphological closing operation: This operation clears the image of small details like cables or lamp posts that would affect the image segmentation. Canny segmentation: This function [4] detects the edges present in the image. Edges are detected in highgradient regions of the image. The output is a binary image, in which the white zones represent the edges. 3) Hough line detection: This function looks for line patterns in the image, based on a test checking if different edge pixels follow a geometrical constraint described by the equation of a line [5]. It connects a set of edge points if they are aligned, allowing to close gaps left by the Canny segmentation. Morphological opening operation: This operation further homogenizes the segmented image, by removing small holes in the image introduced by image noise. Floodfill operation: The floodfill operation changes the color of a certain region to black until an edge is found. As described later, the starting point of the floodfill operation is a satellite perceived as a LOS signal. The filled area will correspond to open sky (LOS), while the unfilled area will correspond to blocked (NLOS) areas.

Figure 3 :Figure 3 : 1 )

 331 Figure 3: Overview of the exclusion algorithm. According to

Figure 5 :

 5 Figure 5: Route to Toulouse centre-ville (left). Highway route (right)

Figure 7 :

 7 Figure 7: Comparison of horizontal accuracies when using the estimated heading in urban environment.

Figure 8 :

 8 Figure 8: Comparison between the cumulative distribution function of the horizontal error.

Figure 9 :

 9 Figure 9: Number of LOS satellites as detected by the camera

Figure 10 :

 10 Figure 10: Number of satellites in the highway environment.

Figure 11 :

 11 Figure 11: Inside a tunnel. Left: original Image. Right: Processed image.

Figure 12 :

 12 Figure 12: Performance under trees. Left: original Image. Right: Processed image. NLOS satellites are in red.

Figure 13 :

 13 Figure 13: 𝐶/𝑁 0 depending on visibility. Green represents LOS and red NLOS.

Table 1 :

 1 Parameters for Pseudorange Error Variance

	Measurement	𝑎	𝑏
	Pseudorange	-1.5	731 2

Table 3 :

 3 Test with Doppler measurement selection in urban environment

	Parameters	Integration with Doppler	𝟒𝒎/𝒔 max innovation	𝟏𝒎/𝒔 max innovation
	Mean
	Horizontal	8.43	5.92	3.33
	Error [m]
	95%
	Horizontal	20.90	14.43	7.84
	Error [m]
	Max.
	Horizontal	33.29	23.74	23.40
	Error [m]
	Mean Yaw Error [deg]	7.79	5.26	2.79
	95% Yaw Error [deg]	24.78	12.72	6.77
	Max Yaw Error [deg]	40.10	19.67	14.26

Table 4 :

 4 accuracy results in the urban environment.

	Configuration	INS/GNSS No camera	Camera with Reference heading	Camera with Estimated heading
	Mean
	Horizontal	4.24	3.74	3.81
	Error [m]
	95%
	Horizontal	10.60	8.96	9.63
	Error [m]
	Max.
	Horizontal	18.67	21.34	22.17
	Error [m]
	Mean heading Error [°]	3.95	2.55	3.71
	95% heading Error [°]	10.11	7.45	8.40
	Max heading Error [°]	18.671	11.52	15.13

Table 5 :

 5 Accuracy results in the highway environment.

		INS/GNS	Camera	Camera
	Configuratio	S	with	with
	n	No	Referenc	Estimate
		camera	e heading	d heading
	Mean
	Horizontal	1.17	1.16	1.18
	Error [m]
	95%
	Horizontal	2.33	2.21	2.25
	Error [m]
	Max.
	Horizontal	5.09	7.03	4.55
	Error [m]
	Mean
	heading	3.78	3.21	3.45
	Error [°]
	95% heading Error [°]	17.69	17.08	17.08
	Max heading Error [°]	39.30	26.90	29.09

funded by GSA (European GNSS Agency), ENAC and Egis Avia.

Amani Ben Afia graduated in 2012 as an Electronics engineer from ENAC (French National Civil Aviation School),


Theresia Gschwandtner
email: gschwandtner@cvast.ac.at

Visual Analytics Meets Process Mining: Challenges and Opportunities

Keywords: Visual Process Mining, Visual Analytics, Challenges

Event data or traces of activities often exhibit unexpected behavior and complex relations. Thus, before and during the application of automated analysis methods, such as process mining algorithms, the analyst needs to investigate and understand the data at hand in order to decide which analysis methods might be appropriate. Visual analytics integrates the outstanding capabilities of humans in terms of visual information exploration with the enormous processing power of computers to form a powerful knowledge discovery environment. The combination of visual data exploration with process mining algorithms makes complex information structures more comprehensible and facilitates new insights. In this position paper I portray various concepts of interactive visual support for process mining, focusing on the challenges, but also the great opportunities for analyzing process data with visual analytics methods.

Introduction

Today we are confronted with an overload of data and information. The amount of information produced on a daily basis is just too much for anyone to process in order to answer specific questions. Examples include data from electronic health records, real time sensors, communication logs, and financial transactions. However, somewhere within this huge amounts of data, there is valuable information to solve specific tasks and answer specific questions. Yet, in order to make sense of the data, we need to find appropriate ways to process it. There are two opposed approaches to tackle this problem: (1) machine learning and data mining, and (2) visualization of the data.

The first approach, i.e. using machine learning and data mining, taps the computational power of the computer for statistical analysis of the data. As a result, a report is generated, giving -at best -the answer to the user's question. However, the second approach, i.e. visualizing the data, utilizes the power of human perception to simultaneously process large amounts of data. When visualizations are carefully designed (e.g., using a suitable combination of preattentive visual attributes to encode important information), the human perception is capable of processing large amounts of data in parallel and of quickly identifying patterns and outliers [START_REF] Mazza | Perception. In Introduction to Information Visualization[END_REF]. Thus, visualization helps to provide an overview of large and complex data. Furthermore, interactive controls allow the user to explore different aspects of the data in more detail.

Visual analytics is defined as 'the science of analytical reasoning facilitated by visual interactive interfaces' [45, p. 4]. It exploits both, the computational power of the computer and the human's perception system to facilitate insights and enable knowledge discovery in large and complex bodies of data. Thus, visual analytics combines both of the aforementioned approaches: machine learning and data mining, and visualization of the data.

Benefits of Visualization

One example that motivates the use of visualizations nicely, is the Anscombe's Quartet: In 1973 the statistician Francis Anscombe constructed four data sets [START_REF] Anscombe | Graphs in statistical analysis[END_REF] to demonstrate the necessity to analyze data with visual means as well as the impact of outliers on statistical properties. The summary statistic properties of these four data sets are nearly identical which could lead to the conclusion that these data sets are nearly identical. However, visualizing these data sets with scatter-plots reveals very different behavior of all four sets.

Another example is given by Spence and Garrison [START_REF] Spence | A remarkable scatterplot[END_REF]: The Hertzsprung Russel Diagram shows the temperature of stars plotted against their magnitude. When looking at this diagram humans can easily distinguish different clusters of types of stars, while automatic means fail due to the noise and artifacts within the data. Card et al. summarize the ways is which information visualization amplifies human cognition like this [START_REF]Readings in Information Visualization: Using Vision to Think[END_REF]:

-Increasing working resources such as by using a visual resource to offload work from the cognitive to the perceptional system -Reducing search such as by representing large amounts of data and grouping information -Enhancing the recognition of patterns such as by visually organizing data by structural relationships -Supporting perceptual inference of relationships by suiting the human's perceptual abilities -Perceptual monitoring of a large number of potential events such as by making changes stand out visually -Providing a manipulable medium that, unlike static diagrams, allows for interactively exploring the data

In contrast to pure question answering tools, such as IBM Watson [START_REF] Ferrucci | Introduction to "This is Watson[END_REF] (a computing system that uses machine learning technologies to answer specific questions), visual analytics combines these technologies with interactive presentations of data and results, and thus, it enables users to derive additional insights from interacting with the data. Hence, users are able to browse the data and they may find answers to questions they did not even know they were looking for.

Visual Approaches for Process Mining

In their process mining manifesto [START_REF] Van Der Aalst | Process mining manifesto[END_REF] van der Aalst et al. emphasize the potential of visual analytics, i.e. the combination of visualizations, interactions, and mining techniques, to enhance process mining. We conducted a literature research starting with the process mining manifesto and references to relevant systems we already knew. Moreover, we used libraries and search engines such as Google Scholar [16] and IEEE Xplore [24] to search for various related keywords and we went through proceedings of important conferences such as IEEE VIS [23], EuroVis [START_REF]EuroVis (EG/VGTC Conference on Visualization[END_REF] conference, and the BPM conference [START_REF]Conferences on Business Process Management[END_REF], as well as through the list of references of relevant papers.

Most process mining approaches employ some kind of visual representation due to the wealth and complexity of the data. Yet, they usually do not explicitly consider visualization as a key factor and rather put emphasis on the mining methods [START_REF]Process Mining: Discovery, Conformance and Enhancement of Business Processes[END_REF][START_REF] Van Der Aalst | Process mining manifesto[END_REF]. Visualizations are used, for instance, for process discovery: to explore the event log data and to understand, reason about, and fine-tune the derived models (e.g. [START_REF] Jagadeesh | Process diagnostics using trace alignment: Opportunities, issues, and challenges[END_REF][START_REF] Günther | Fuzzy mining -adaptive process simplification based on multi-perspective metrics[END_REF][START_REF] Wongsuphasawat | Exploring flow, factors, and outcomes of temporal event sequences with the outflow visualization[END_REF][START_REF] Maguire | Visual compression of workflow visualizations with automated detection of macro motifs[END_REF][START_REF] Hipp | Navigating in process model collections: A new approach inspired by google earth[END_REF][START_REF] Vrotsou | Activitree: Interactive visual exploration of sequences in event-based data using graph similarity[END_REF][START_REF] Günther | Disco: Discover your processes[END_REF]). To this end, flow charts [36] and directed graphs are used most frequently. Other examples use visualizations to identify interesting behavior and patterns in the event log data (e.g., [START_REF] Song | Supporting process mining by showing events at a glance[END_REF][START_REF] Bodesinsky | Exploration and assessment of event data[END_REF]), or to investigate process conformance (e.g. [49,[START_REF] Adriansyah | Conformance checking using cost-based fitness analysis[END_REF][START_REF] Bodesinsky | Visual analysis of compliance with clinical guidelines[END_REF]). In this position paper I give by no means a comprehensive list of visualization approaches in the context of process mining. Yet, I introduce a selection of our own approaches, highlighting some alternative visual representations and interaction techniques in order to spark ideas how to combine process mining with visual analytics.

Visual Analysis of Guideline Conformance

A related topic to the analysis of process conformance is the analysis of conformance of clinical actions with the recommended care process. Clinical practice guidelines give specific recommendations on how to treat a patient in a given clinical situation. They are aimed at assisting the care process and ensuring its quality [START_REF] Field | Clinical Practice Guidelines: Directions for a New Program[END_REF]. Hence, they can be seen as process models. Bodesinsky et al. [START_REF] Bodesinsky | Visual analysis of compliance with clinical guidelines[END_REF] presented a visual interactive approach to analyze the compliance of clinical care with these guidelines. The approach is aimed at supporting (1) physicians at the point of care to display clinical guidelines and point them to omitted clinical actions, and (2) medical experts in retrospective analysis of the quality of the treatment and also of the quality of the clinical practice guideline. It visualizes the recommended sequence of actions (see Figure 1 (a)) together with the actual clinical actions that were applied (see Figure 1 (c)). Applied clinical actions are represented by diamonds on a time axis (a pop-up window gives the user details about these actions when hovering them). Moreover, the patient's parameters (see Figure 1 (b)) and the compliance of actions with the clinical guidelines (see Figure 1 (c)) are represented on a time axis. To aid the analysis of conformance, they visually distinguish conform and non-conform actions, as well as missing actions. Conform actions are indicated by diamonds and non-conform actions are indicated by diamonds marked with an X. Any delay of application is indicated by two parallel colored bars, marking the time span between the intended time of application and the actual time of application (see Figure 1 (c)).

Plan Strips

Another example from the medical domain is the Plan Strips visualization [START_REF] Seyfang | Visualizing complex process hierarchies during the modeling process[END_REF]. This interactive approach visualizes the complex hierarchies of processes -again using the example of clinical practice guidelines -by a set of nested rectangles. These guidelines recommend treatment plans (groups of clinical actions) which are composed of sub-plans. The execution order of a plan determines in which of the following ways its sub-plans are to be carried out: in sequence, in parallel (all sub-plans start at the same time), unordered (in sequence, in parallel, or a mixture of both), or in any-order (in sequence but the order is not defined) [START_REF] Shahar | The Asgaard Project: A task-specific framework for the application and critiquing of time-oriented clinical guidelines[END_REF]. In addition, cyclical plans represent plans which are carried out in loops. Plan Strips encode these different kinds of sub-plan synchronization by color (see Figure 2): Gray rectangles indicate recommended clinical actions or activities while colored containers represent the corresponding parent-plans as well as the execution order. Interactive selection and highlighting techniques are used to ease navigation. One limitation of this approach is that the visualization demands for some learning effort, however, Plan Strips demonstrate an alternative way of representing highly specified and complex conditions of a process model in a space efficient way.

EventExplorer

EventExplorer [START_REF] Bodesinsky | Exploration and assessment of event data[END_REF] was designed to aid the analyst in exploring event log data, in order to gain some understanding of the data and the dependencies between events before deciding for concrete analysis methods or process mining algorithms. To this end, it combines pattern mining with interactive visualizations. Moreover, it supports the exploration of events on a temporal scale to aid the identification of important temporal dependencies of events (e.g., two events always occur together). Figure 3 shows a screenshot of EventExplorer. Each case (i.e., sequence of events) of the dataset is represented by a horizontal line composed of colored rectangles which represent the different events: The color indicates the type of event, while the order indicates the sequence of events.

The user can switch the representation to lay out events along a time axis to investigate their precise temporal behavior. EventExplorer uses pattern mining techniques to identify patterns within sequences of events. When a pattern is selected, all instances of this pattern within one case are highlighted and connected by gray arcs which emphasize the recurrence of this pattern within the case. EventExplorer also supports fuzzy pattern mining by allowing the use of wildcards when entering a pattern of interest. The approach presents interactive visualizations to (1) get familiar with the data, (2) to cluster and sort cases, and (3) to identify temporal patterns. The arcs visually emphasize recurring patterns as well as temporal distances (when using the time axis view). Fig. 3. EventExplorer [START_REF] Bodesinsky | Exploration and assessment of event data[END_REF]. In (a) EventExplorer represents all cases of the dataset as horizontal lines with color coded events (scrollable). (b) shows the different patterns (together with their number of occurrences) which match the settings in (c), i.e., the length of the pattern and the minimum number of occurrences. Alternatively, the user can enter a specific pattern of interest into the text field in (b). Selected patterns (in this case 'V-R-V') are highlighted by gray arcs in (a). The user can switch the x-axis in (a) from a sequential axis to a time axis (d) in order to explore the precise temporal characteristics of event sequences.

The three approaches outlined in detail were designed for very diverse purposes and give only examples of applications of visual analytics in process mining. However, there are still many open challenges but also opportunities when looking at the combination of these two research fields. We derived these challenges from the literature review, our previous work, and our long lasting experience in the research field of visual analytics.

Challenges and Opportunities

Researchers have identified challenges of today's process mining research in various works (e.g., [START_REF] Mans | Process mining in healthcare: Data challenges when answering frequently posed questions[END_REF][START_REF] Van Der Aalst | Process mining manifesto[END_REF]). In [START_REF] Van Der Aalst | Process mining manifesto[END_REF] van der Aalst et al. even mention the potential of visual analytics to support process mining. Others again have summarized open challenges in the field of visual analytics (e.g., [START_REF] Keim | Visual Analytics: Definition, process, and challenges[END_REF][START_REF] Thomas | Illuminating the Path: The Research and Development Agenda for Visual Analytics[END_REF][START_REF] Wong | The top 10 challenges in extreme-scale visual analytics[END_REF]). When focusing on the combination of visual analytics and process mining, we face related challenges.

C1: Intertwining Process Mining with Visual Analytics

Van der Aalst et al. identify the combination of process mining with other types of analysis as one of the challenges in process mining [START_REF] Van Der Aalst | Process mining manifesto[END_REF]. Moreover, they explicitly mention the potential of process mining combined with visual analytics for leading to further insights. Yet, how to intertwine process mining with visual analytics techniques is still an open problem. There are a couple of approaches which already tackle this problem. Among these visual approaches to aid process mining, the use of flow charts [36] or directed graphs to represent the process model are most common (e.g., [START_REF] Günther | Fuzzy mining -adaptive process simplification based on multi-perspective metrics[END_REF][START_REF] Maguire | Visual compression of workflow visualizations with automated detection of macro motifs[END_REF][START_REF] Hipp | Navigating in process model collections: A new approach inspired by google earth[END_REF][START_REF] Günther | Disco: Discover your processes[END_REF][START_REF] Seyfang | Visualizing complex process hierarchies during the modeling process[END_REF]). Visualizing the process model helps to understand and fine-tune the model, but also to check conformance. Other visual approaches show an overview of individual event sequences (e.g., [START_REF] Song | Supporting process mining by showing events at a glance[END_REF][START_REF] Jagadeesh | Process diagnostics using trace alignment: Opportunities, issues, and challenges[END_REF][START_REF] Wongsuphasawat | Exploring flow, factors, and outcomes of temporal event sequences with the outflow visualization[END_REF][START_REF] Vrotsou | Activitree: Interactive visual exploration of sequences in event-based data using graph similarity[END_REF][START_REF] Bodesinsky | Exploration and assessment of event data[END_REF]) which can be used to identify common paths and thus, to derive a process model. Visualizations that show the conformance of individual event sequences (e.g., [START_REF] Bodesinsky | Visual analysis of compliance with clinical guidelines[END_REF][START_REF] Adriansyah | Conformance checking using cost-based fitness analysis[END_REF]), however, help to identify shortcomings of the existing model, but also variations of individual cases. Yet, there is still much unexploited potential of visual analytics to better support pattern discovery and the derivation of process models, fine-tuning and enhancing these models, as well as identifying conformance problems and finding alternative solutions. To this end, visualizations need to be tightly coupled to analytical techniques and mining algorithms in an iterative loop, giving immediate visual feedback to changes and adjustments made to the algorithms. Appropriate views and interactions need to be tailored to specific user groups with respect to their specific tasks and data [START_REF] Miksch | A matter of time: Applying a data-users-tasks design triangle to visual analytics of time-oriented data[END_REF].

C2: Scalability and Aggregation

The scalability of visualizations is a well-known challenge in visual analytics [START_REF] Keim | Visual Analytics: Definition, process, and challenges[END_REF][START_REF] Thomas | Illuminating the Path: The Research and Development Agenda for Visual Analytics[END_REF][START_REF] Wong | The top 10 challenges in extreme-scale visual analytics[END_REF]. This is also true in the context of process mining, which usually deals with huge amounts of event sequences that need to be analyzed. On the one hand, process mining algorithms tackle this problem by reducing these huge event logs to a manageable amount of patterns and process models. On the other hand, the analysis of the raw data and single event sequences may yield some important insights, and thus, visualizations need to support the simultaneous analysis of huge amounts of individual event sequences. There are some approaches that tackle this problem to some extent (e.g., [START_REF] Jagadeesh | Process diagnostics using trace alignment: Opportunities, issues, and challenges[END_REF][START_REF] Wongsuphasawat | Exploring flow, factors, and outcomes of temporal event sequences with the outflow visualization[END_REF][START_REF] Bodesinsky | Exploration and assessment of event data[END_REF]) by simultaneously showing a number of individual event sequences at minimum display space. Still the number of individual event sequences that can be displayed simultaneously is limited and it is an open challenge how to design interactive visualizations to support the whole analysis cycle from single event sequences to multiple event logs. There is a need for scalable visualizations able to represent a big number of cases while allowing zoom-in and drill-down to investigate details about single event sequences, in analogy to the famous Information Seeking Mantra by Shneiderman: 'overview first, zoom and filter, then details-on-demand' [41, p.337].

C3: Interaction to Support Process Discovery and Enhancement

This challenge is tightly coupled with the last challenge of finding appropriate visualization techniques to support the different tasks. Data exploration is an iterative process with trial and error loops. Suitable interaction techniques need to support drilling down from an overview visualization to the investigation of single events. Commonly used techniques include zoom and filter functionalities, however, other interaction techniques may be needed depending on the data, user, and task [START_REF] Miksch | A matter of time: Applying a data-users-tasks design triangle to visual analytics of time-oriented data[END_REF]. Further interaction techniques include select, explore, reconfigure, encode, abstract/elaborate, and connect [START_REF] Yi | Toward a deeper understanding of the role of interaction in information visualization[END_REF]. Allowing the analyst to interactively investigate the raw event data as well as the effects of changes made to models and mining algorithms may foster new insights and improve process mining results.

EventExplorer [START_REF] Bodesinsky | Exploration and assessment of event data[END_REF] allows for sorting cases by pattern count or by time. Another interesting way of sorting would be by similarity of event sequences in order to aid the identification of different groups of cases. Moreover, it supports the investigation of events in sequence or on a temporal scale. An interesting addition would be allowing the analyst to apply different ways of aligning cases vertically, for instance, aligning cases with the first appearance of a specific event.

C4: Data Quality and Uncertainty

Data quality is a general challenge in visual analytics [START_REF] Kandel | Research directions in data wrangling: Visuatizations and transformations for usable and credible data[END_REF], as well as in process mining [START_REF] Van Der Aalst | Process mining manifesto[END_REF], and thus, data quality is just as important when it comes to the combination of both fields. While process mining algorithms require good data quality and well-structured event logs, in practice these can be erroneous and badly structured. Instances of the same event type may have different names, event granularities may vary, and quite often the event log contains missing, incorrect, imprecise, uncertain, or irrelevant data. Data quality control can be divided into (1) data profiling (i.e., identifying and communicating quality problems), (2) data cleansing (i.e., correcting erroneous data), and (3) data transformation (i.e., transforming the data into appropriate formats for automatic processing) [START_REF] Gschwandtner | A taxonomy of dirty timeoriented data[END_REF][START_REF] Kandel | Research directions in data wrangling: Visuatizations and transformations for usable and credible data[END_REF]. There are visual approaches that tackle the problem of data quality (e.g., [START_REF] Gschwandtner | TimeCleanser: A visual analytics approach for data cleansing of timeoriented data[END_REF][START_REF] Bors | Qualitytrails: Data quality provenance as a basis for sensemaking[END_REF][START_REF] Kandel | Profiler: Integrated statistical analysis and visualization for data quality assessment[END_REF][START_REF] Kandel | Wrangler: Interactive visual specification of data transformation scripts[END_REF]), however, the specific needs in process mining are hardly considered, such as case heterogeneity, event granularity, or concept drift.

On the other hand, data often contains some amount of uncertainty. Event log data may, for instance, contain uncertainties about which event type corresponds to which log entry or about the exact time of an event. Visually communicating these uncertainties to the analyst (e.g., [START_REF] Gschwandtner | Visual encodings of temporal uncertainty: A comparative user study[END_REF]) is important for a better informed reasoning.

C5: Complexity of Time-Oriented Data

Usually event log data has an inherent temporal structure. For any such event log at least the sequence of events is known and used to analyze paths and processes. Yet, there are also event logs with precise timestamps for each event, allowing for a detailed analysis regarding the temporal dimension of events and event sequences. Time, however, is a complex data domain with very special features that require special consideration [START_REF] Aigner | Visualization of Time-Oriented Data[END_REF][START_REF] Andrienko | Exploratory Analysis of Spatial and Temporal Data: A Systematic Approach[END_REF][START_REF] Shneiderman | The eyes have it: A task by data type taxonomy for information visualizations[END_REF]. Some examples of these special characteristics are: time can be given as time points or intervals, conventionally time is aggregated irregularly (i.e., 60 minutes per hour, 24 hours per day, 28-31 days per month, 356/366 days per year, ...), leap seconds and days, different time zones, recurrences (e.g., seasonal cycles), and it has a strong social component (e.g., weekends, holidays) [START_REF] Aigner | Visualization of Time-Oriented Data[END_REF]. Ailenei et al. [START_REF] Ailenei | Definition and validation of process mining use cases[END_REF] outline use cases with special consideration of the temporal behavior of processes that need to be considered in process mining. When it comes to the visualization of time-oriented data, Aigner et al. [START_REF] Aigner | Visualization of Time-Oriented Data[END_REF] give a comprehensive overview. However, considering these peculiar characteristics of events (e.g., type of day, some cyclic behavior synchronized with the calendar, or events that co-occur with a certain delay) and adequately visualizing them for the task at hand is still an open challenge.

C6: Evaluation

Since visual analytics aims at supporting human reasoning and gaining new insights, the quality of visual analytics approaches is hard to quantify. Plaisant [START_REF] Plaisant | The challenge of information visualization evaluation[END_REF] outlines different challenges of evaluating visualizations and emphasizes that evaluation strategies need to take the exploratory nature of tasks and the added value of visualization such as overall awareness and potential discoveries into account. Moreover, praxis related aspects such as a successful adoption of these visualizations need to be considered. There are different approaches for evaluating visualizations [START_REF] Shneiderman | Strategies for evaluating information visualization tools: Multi-dimensional in-depth long-term case studies[END_REF][START_REF] Isenberg | A systematic review on the practice of evaluating visualization[END_REF][START_REF] Forsell | Handbook of Human Centric Visualization, chapter An Introduction and Guide to Evaluation of Visualization Techniques Through User Studies[END_REF] which usually involve target user participation and are highly task dependent (see [START_REF] Andrienko | Exploratory Analysis of Spatial and Temporal Data: A Systematic Approach[END_REF] for a task framework and [START_REF] Lammarsch | Developing an extended task framework for exploratory data analysis along the structure of time[END_REF] for a timespecific task framework). In any case it is recommended to pursue an iterative design approach which involves feedback of the target user as early as possible.

On the other hand, evaluation in process mining is an open challenge too. There are different approaches on how to evaluate derived process models (e.g., [START_REF] Rozinat | The need for a process mining evaluation framework in research and practice[END_REF][START_REF] Ardagna | New Quality Metrics for Evaluating Process Models[END_REF]), and Ailenei et al. [START_REF] Ailenei | Definition and validation of process mining use cases[END_REF] present a set of task specific use cases for the evaluation of process mining tools. These use cases may as well help to evaluate visual analytics solutions for process mining from a technical point of view. However, the evaluation of gained awareness and insights is still an open challenge.

Conclusion

Handling huge amounts of process data can be tackled in different ways: automatic mining techniques, visualizations, or by visual analytics which combines the benefits of both fields. Especially the latter has potential to improve the results of process analysis in many ways. In this paper I highlighted some examples of how this was approached in recent years. Moreover, I outlined six open challenges which reveal opportunities for further research. I believe, visual analytics is beneficial to many tasks of process mining and -if done right -will have substantial impact on the quality of process analysis.

Fig. 1 .

 1 Fig. 1. Visual analysis of guideline conformance [8]. The recommended care process is visualized in a flow chart like graph in (a). In this screenshot the graph is zoomed in to show two repeated actions that should be carried out in parallel. (b) shows the patient's parameters that are affected by these actions. (c) represents actual clinical actions on a time axis. Actions that conform to the guideline are indicated by diamonds, nonconform actions are diamonds marked with an X, and the delay of actions is indicated by colored horizontal bars. A pop-up window gives the user details about these actions when hovering them.

Fig. 2 .

 2 Fig. 2. Plan Strips [39] represent the hierarchical structure (a) of clinical plans and subplans as nested rectangles (b). User performed actions are indicated by gray rectangles (light gray when selected), while the color of parent containers indicate the execution order.

ACKNOWLEDGMENTS

The research leading to these results has received funding from the Centre for Visual Analytics Science and Technology CVAST, funded by the Austrian Federal Ministry of Science, Research, and Economy in the exceptional Laura Bassi Centres of Excellence initiative (#822746).


Alain Jean-Marie
email: alain.jean-marie@inria.fr

marmoteCore: a Markov Modeling Platform

Keywords: • Mathematics of computing → Markov processes, Solvers, • Networks → Network performance modeling, Network performance analysis, Markov modeling, software environment

come

and general modeling constructions and solution methods, all accessible using an uniform programming interface.

marmoteCore is the prototype of such a platform. It consists in a programming library 1 with an object-oriented point of view, that allows programmers to create Markov models and "solve" them by using many available algorithms. These solvers are often embedded algorithms or applications written by third parties. Another feature of the library is the availability of ready-made Markov models of several levels of generality, with adapted solution methods.

What is not marmoteCore. The intent of the project is not to program anew all methods that already exist for creating or analyzing Markov chains. The idea is not either to develop an integrated modeling and analysis tool: this sort of functionality is offered by Workflow Management Systems (WMS, see Section 5.4) or integrated platforms such as DTK. 2The purpose of this document is to demonstrate the potentialities of marmoteCore's architecture. This is not a programming guide, yet several code excerpts will give a flavor of the way programs are written with this library. Indeed, while still being in development, marmoteCore is mature enough so as to be usable in practice.

The paper is organized as follows. Section 2 presents the general purpose and architecture of the software. Section 3 presents the core of the programmer's interface. In Section 4 we describe the principal methods of analysis currently available. In Section 5, we discuss the way marmoteCore interacts with existing scientific software. We conclude with our plans for the future in Section 6.

Acknowledgment. The existence of marmoteCore owes a lot to the contribution of Issam Rabhi, Emmanuel Hyon and Hlib Mykhailenko. The author wishes to acknowledge their contribution to the architecture and the code of the software.

MAIN CHARACTERISTICS 2.1 Background

The idea of providing users with software environments to help them realize their modeling experiments is, of course, not original. The most popular mathematical modeling environments, such as Matlab, Scilab, Mathematica, Maple, R, Sage, all provide packages with functions specialized to families of models (e.g. Matlab's SimuLink for system modeling), or analysis techniques (e.g. statistical packages).

To this day however, there does not exist standard packages specifically devoted to modeling with discrete space Markov chains, considered as entities richer than simple matrices. In the particular context of performance modeling, the list of software maintained at [START_REF] Hlynka | List of queueing theory software[END_REF] is typical of the fact that there is a large variety of solutions techniques adapted to specific models or families of models, but also that there is no uniform presentation or unique access point to these tools, that would help the newcomer to apply one to his/her needs. On the other hand, analysis tools are sometimes presented as "packaged" behind a graphical user interface, as were pioneers like GreatSPN, Tangram-II, and are modern tooks like Java Modeling Tools, CosyVerif, OpenAlea, to name just a few. In some cases, these tools and their advanced algorithms can hardly be used outside of their application context.

We believe there is a need for libraries giving access to efficient numerical algorithms at a lower programming level.

The general purpose of marmoteCore

The development of the marmoteCore software pursues several objectives simultaneously. On the one hand, it aims at providing to the general scientific user a "modeling environment" for Markovian systems. This environment should provide at the same time an access to basic models of the literature (an thus serve as a repository of models, together with their "classical" results), and an access to state-of-the-art solution algorithms. Through a programming interface, this Markov user should be able to use marmoteCore's library as "routine" in a scientific project.

In the practice of Markov modeling for the performance evaluation of concrete situations, it may happen that the modeler defines a family of parametric models and develops a specific experimental setup centered around this specific family. In that case, modeling with the aid of a graphical interface, or developing a set of ad-hoc scripts is convenient. In other situations, and especially in other fields of science, "solving" a Markov chain is just a subroutine in a much larger process involving statistical estimation, learning, decisions, etc. In those situations, providing solution methods through an API is a requirement of the end user.

At a second level, the environment should also be useful to Markov developers who are working on advanced solution methods, be they generic or adapted to specific classes of Markov models. marmoteCore will then be used as a library of algorithms as well as benchmark cases, and will serve as an experimental platform for comparing the performance of methods. It would then be seen as some middleware, providing a communication between some generic scientific software or some experiment management environment (e.g. a WMS) and solution algorithms/programs.

In coping with these simultaneous objectives, the main challenge is to produce a framework as open and flexible as possible, in order to welcome software contributed by external parties, while at the same time allowing the general user to use their favorite modeling language/software.

Architectural choices

2.3.1 Object-oriented organization. The objectives of the project make the choice of object-oriented programming almost obvious. Object-oriented languages give the possibility of designing high-level abstractions for mathematical concepts, representing objects with common properties. Yet, through the mechanism of inheritance, the user has the flexibility of controlling the implementation of specific instances of the model. The organization of marmoteCore relies a lot on hierarchies of models, and polymorphism. We illustrate the power of this idea in Section 3.

We have chosen the C++ language, partly to ease up the integration of legacy code available in C or C++. The first phase of the development, presented here, consisted in the development of a solid code base demonstrating the potentialities of the architecture. This choice is consistent with the longerterm objectives of the project: in a second phase, we will make the library available as packages, plugins or libraries to other languages/ platforms. The most popular ones (R, Python, Sage, ...) indeed all provide a way to incorporate functions written in C++.

Separation between the core and packages through wrapping.

The idea that marmoteCore should provide a core of objects and functionalities, and give access to more functionalities through wrapping of packages, comes from the following factors. On the one hand, there exists multiple incompatible data formats for models (typically, the matrix representation of Markov generator) and results of analyses (distributions, trajectories, etc.). On the other hand, there exists a wealth of software programmed in a variety of languages, and it is out of question to re-program all these applications.

A solution is to provide a small number of high-level abstractions that are common to any Markov model, and group them into a "core" which can be seen as an exchange hub. The core should be able to handle Markov models in various formats and call solution methods from different sources. This is done using three principal paradigms: (1) native programming within the core, using marmoteCore objects; (2) reusing methods already programmed, typically C or C++;

(3) wrapping of external calls to separate applications. The latter usually consists in a) generating a file suitable for some external application; b) execute this application through a "system call", directing the output to some other file; c) parsing of results into marmoteCore. We explain in Section 5 how these techniques are used concretely.

THE MARMOTECORE INTERFACE

The interface of marmoteCore is based on only four principal abstractions: MarkovChain, TransitionStructure, Marmote-Set and Distribution. These high-level classes are intended to be specialized for specific purposes. We present now their principal features: attributes and functionality, and give examples of this specialization.

Class MarkovChain

Obviously, the class MarkovChain is the central one for a Markov modeling software. Technically however, as it appears from Figure 1, it is just a container for the descriptive elements: a state space, a transition structure and an initial distribution. Those are described in the following sections. The class has a large number of methods, which are reviewed in Section 4.

Implementations: a hierarchy of Markov Chains.

The principal novelty and full potential of the software resides in its organization in a hierarchy of concepts. Markov Chain and Transition Structure objects in particular can be organized from the most specific to the most generic.

A target hierarchy for Markov chains in the marmoteCore environment is displayed in Figure 2 (for continuous-time chains: a similar hierarchy exists for discrete-time ones). It illustrates the project to populate the environment with This proposal is by no means exhaustive, and can be adapted to host more "high level" modeling paradigms (Stochastic/Timed Petri Nets, Stochastic Automata Networks, Stochastic Process Algebras) as well as "low level" models with few parameters and strong structure stemming from various fields: population models, epidemics, etc.

Figure 3 illustrates the use of polymorphism, that enables to use, for some model in the hierarchy, methods that have been attached to models higher in the hierarchy. In this toy example, several ways of computing or approximating the stationary distribution are used. This example also illustrates that the same functionality can (and should) be offered with several implementations.

Class MarmoteSet

The MarmoteSet class implements discrete sets of states, with elementary objects such as integer intervals, and elementary constructions such as Cartesian products and unions.

3.2.1 Functionalities. The principal methods of Marmote-Set objects are represented in Figure 4. 3Methods Index() and DecodeState() provide the twoway maps from possibly complex state spaces, represented as vectors of numbers, to a linear numbering of states. They are necessary since numerical methods usually involve the manipulation of arrays with a standard numbering of entries. In some sense, they are also sufficient to implement functionalities using state spaces.

Methods FirstState() and NextState() respectively set some state vector to the conventional "initial" state of the set, and transform this state vector into the following state in the conventional order. Together with IsFirst() which tests equality of a state with the initial state, they are enough to enumerate all states of some set, as shown in Figure 7.

Implementation.

Currently implemented sets include MarmoteInterval, the standard integer interval, Marmote-Box, its multidimensional extension, BinarySequences, the representation of {0, 1} N suitable for many models in computer science and statistical physics, and more complex geometries such as Simplex and BinarySimplex.

Class TransitionStructure

The TransitionStructure class is an abstraction for the labeled state-to-state transitions that are common to discretetime and continuous-time Markov chains. The mathematical representation of this concept is the matrix, with the usual convention that origin states are rows and destination states are columns.

Functionalities.

The principal methods/functionalities for this object are listed in Figure 5. There are basic methods for constructing the structure: setting and accessing the entries. A somewhat advanced access to the structure is through method TransDistrib() which returns a probability distribution over the destination states. It is particularly useful for Monte-Carlo simulation.

Methods Uniformize() and Embed() provide the standard ways for passing from continuous to discrete time.

Methods EvaluateMeasure() and EvaluateValue() represent respectively the action of the operator on row vectors (usually, probability measures) and column vectors (usually, rewards or weights of some sort, attached to states). These correspond to left-and right-vector/matrix multiplications in the matrix algebraic representation. if ( i == j ) return r ; else if ( i == j + 1 ) return p ; else if ( i == j -1 ) return q ; else return 0.0; }; spaces. Some analysis methods, e.g. simulation, can work with these objects. The "block transition" class, suitable for modeling QBDs, currently being implemented, has also this feature.

The class EventMixture is adapted to "event modeling", where elementary events, occuring at constant rates/probabilities provoke elementary transitions (births, deaths, movements, etc.). These are the basic elements for building many complex models, in the view of Poisson Systems. They can also be used to construct Markov decision models or Markov games.

Class Distribution

The Distribution class implements classical attributes for probability distributions, including the calculation of moments, Laplace transforms, and sampling of pseudo-random values. The specific objects implemented in marmoteCore include the ubiquitous Dirac, Bernoulli, Geometric, Poisson and Exponential distributions. The class DiscreteDistribution is adequate for generic distributions on state spaces.

ANALYSIS METHODS

We survey in this section some analysis and "solution" methods presently available in marmoteCore. Those are originally present in existing software packages: the presentation also serves for illustrating the way existing Markov software can be integrated into the platform.

Structural Analysis

Structural analysis methods are related to the communication structure in the transition graph. They include Absor-bingStates(), Recurrent/CommunicatingClasses(), which compute state space decompositions into corresponding classes. Their results are used by IsIrreducible() and IsAccessible(i,j), as well as Period(). The method SubChains() extracts the set of independent MarkovChain objects corresponding to recurrent classes of some Markov chain.

One implementation of these methods uses algorithms of [START_REF] Jarvis | Graph-theoretic analysis of finite Markov chains[END_REF] (see also [START_REF] Stewart | Introduction to the numerical Solution of Markov Chains[END_REF]) suitably modified to be parallelized on multi-core machines when available. Another implementation uses the R markovchain package [START_REF] Spedicato | The markovchain Package: A Package for Easily Handling Discrete Markov Chains in R[END_REF]. See also Section 5.1. This illustrates again the possibility of having the same functionality implemented with different algorithms.

Monte-Carlo simulation

The standard forward Monte-Carlo simulation is available with the SimulateChain() method. The method includes several controls for the printing of the trajectory and the collection of empirical occupancy measures. Perfect sampling of Markov chains using the method of Propp and Wilson [START_REF] Propp | Exact Sampling with Coupled Markov Chains and Applications to Statistical Mechanics[END_REF] is also available. This method is imported from the PSI package [START_REF] Vincent | PSI, Perfect Simulator. Online[END_REF], see also Section 5.2.

Steady-state distributions

One of the aims of the marmoteCore project is to give access to several methods for solving the same problem, using a relatively common interface. This is the case for the computation of (approximate) steady-state distributions for Markov chains, the basic problem of the field.

In the current version, five methods are available for the top-level class MarkovChain. One is the standard power method. Two provide an interface to the linear-algebraic methods GTH and SOR (see Stewart [START_REF] Stewart | Introduction to the numerical Solution of Markov Chains[END_REF]): these are borrowed from the XBorne package [START_REF] Fourneau | XBorne 2016: A Brief Introduction[END_REF], see Section 5.3. The last one is an interface to the functionality offered by the R package markovchain already mentioned [START_REF] Spedicato | The markovchain Package: A Package for Easily Handling Discrete Markov Chains in R[END_REF]. In addition, the hierarchical organization of Markov models allows to derive specific methods for models with wellknown structure. This is the case for classes TwoState, Homo-geneous1DBirthDeath (aka M/M/1/K), Homogeneous1DRan-domWalk, Felsenstein81, representing (a partial list of) models for which the steady-state distribution can be found in closed form. Abstraction allows to handle some infinite-state models: for the infinite-state random walk and birth-death, an object of type GeometricDistribution is returned by StationaryDistribution().

Transient distributions

Computing transient distributions for general Markov models is a highly technical task and marmoteCore applies its "reuse" philosophy to delegate it to specialized tools. Currently, the TransientDistribution() method of the toplevel MarkovChain uses the package described in [START_REF] Cerdà-Alabern | Closed form transient solution of continuous time Markov chains through uniformization[END_REF].

Similarly to stationary distributions, closed form solutions that exist can be implemented at lower levels of the hierarchy. This is the case for Homogeneous1DRandomWalk, Homogeneous-1DBirthDeath, Poisson, TwoState, Felsenstein81.

Hitting times

A method is provided to obtain (approximate) average hitting times of an arbitrary subset of the state space, starting from every initial state. Hitting time distributions can be sampled using Monte-Carlo simulations. The exact value for this distribution is returned for some elementary models, e.g. currently: TwoState, Poisson, Felsenstein81.

INTEGRATION WITH EXISTING SCIENTIFIC SOFTWARE

In this section, we explain how marmoteCore can (and does) interface with existing software systems. Indeed, it is essential to the project that existing methods of analysis (solution algorithms, data assimilation and model construction, e.g. [START_REF] Horváth | PhFit: A General Phase-Type Fitting Tool[END_REF], statistical analysis, ...) be accessible easily through the programming interface. We explain through several examples that this idea works well. The functionalities of mar-moteCore are limited only by the time needed to realize the interfaces, and the accessibility of the software to be integrated.

R

The case of R is particularly interesting because there exist two ways of interacting with it. One first way is through wrapping (see Section 2.3), as for any external program. A second, more intimate way, is to create an instance of R's execution engine inside the program, using the Rinside library. Instructions are passed to this engine, and results recovered as C++ objects. We have used this technique to wrap some functionalities of the markovchain package [START_REF] Spedicato | The markovchain Package: A Package for Easily Handling Discrete Markov Chains in R[END_REF] and the transient probabilities solver [START_REF] Cerdà-Alabern | Closed form transient solution of continuous time Markov chains through uniformization[END_REF]. The method Sample() of the class PoissonDistribution also uses R: a validation that the non-trivial statistical methods need not be reimplemented.

PSI

PSI, the Perfect SImulator, is a unix tool that provides a simulation kernel for continuous-or discrete-time Markov chains, based on backwards coupling [START_REF] Vincent | PSI, Perfect Simulator. Online[END_REF]. For marmoteCore, this represents the example of wrapping the call of external applications. Markov matrices are written in the MARCA format [START_REF] Stewart | Introduction to the numerical Solution of Markov Chains[END_REF] suitable for executing the psi sample application that generates exact samples of the stationary distribution.

XBorne

XBorne is a set of independent programs, written in C, for building and transforming stochastic matrices, and computing bounds or exact results for their steady-state and transient distributions [START_REF] Fourneau | XBorne 2016: A Brief Introduction[END_REF]. For marmoteCore, this is an example of direct reuse of C code. Courtesy of the package's maintainers, applications from XBorne have been slightly modified so as to provide their functionality through function calls. 4 Communication of the Markov model is done through files, although in principle the translation of data structures could be performed. The methods SOR and GTH for computing stationary distributions [START_REF] Stewart | Introduction to the numerical Solution of Markov Chains[END_REF] are included this way in marmo-teCore.

Workflow Management Systems

marmoteCore will ultimately offer a large variety of analysis tools. The effective use of those features within a complex experiment will probably be easier by interfacing marmote-Core with a Workflow Management System (WMS).

We have experimented this possibility with a well-known WMS: Kepler. The Kepler Project 5 is an open source workflow application, designed to help scientists, for designing, executing, reusing, evolving, archiving, and sharing scientific 4 Technically: using extern "C" declarations. 5 https://kepler-project.org/.

workflows. Kepler is a java-based application: it was necessary to wrap our C++ code into it. FC2K (Fortran/C/C++ to Kepler) is a tool which offer this functionality by generating a Kepler actor that executes a given C/C++ function or Fortran routine via JNI (java native interface).

CONCLUSION AND FUTURE WORK

We have presented the general orientations of the marmote-Core software project. These features are implemented in the current release of the core, and are fully operational: marmo-teCore has been used to realize numerical and simulation experiments in [START_REF] Dimitriou | A Markovian queueing system for modeling a smart green base station[END_REF][START_REF] Vallet | Semiconductor laser Markov models in the micro-canonical, canonical and grand-canonical ensembles[END_REF]. The current development plans include the development of QBD-like formalisms, interfacing available QBD software. The extension of modeling capabilities towards Markov Reward /Accumulation processes (MMPP, MAP, MArP, aka Markov Arrival Processes, of Figure 2), controlled Markov chains and Markov games is also foreseen. Finally, we will start the second phase of the project with the development of a Python interface.

Figure 1 :

 1 Figure 1: Attributes of the MarkovChain class

Figure 2 :

 2 Figure 2: Hierarchy of continuous-time Markov chains; arrows represent the "is a" relationship

Figure 3 :

 3 Figure 3: Illustration of polymorphism for the call to solution methods

Figure 4 :

 4 Figure 4: Minimal methods of the MarmoteSet class 3.3.2 Implementations. This class has currently several implementations. The one that is dedicated to numerical computation, uses a sparse matrix structure and consistently is named SparseMatrix.Other implementations are not based on the comprehensive storage, for instance MultiDimHomTransition that represents the multi-dimensional birth/death transitions with constant rates or probabilities. In the code excerpt in Figure6(for the 1-D case), no bounds on state numbers appear: it is possible to represent this way processes on infinite state

Figure 5 :

 5 Figure 5: Methods of a TransitionStructure

Figure 6 :

 6 Figure 6: Implementation of an unbounded birthdeath transition structure

Figure 7 :

 7 Figure 7: Generic way of constructing an infinitesimal generator

Available at http://marmotecore.gforge.inria.fr.

http://dtk.inria.fr/

In this figure and following ones, methods written in italics are virtual and must be implemented in derived classes.

This research was funded by the French National Research Agency under grant ANR-12-MONU-00019, project "MARMOTE" (MARkovian MOdeling Tools and Environments), https://wiki.inria.fr/ MARMOTE/Welcome. The author is also affiliated to LIRMM, Université de Montpellier/CNRS.