N Boukhelifa

W Cancino

A Bezerianos

E Lutton

Evolutionary Visual Exploration: Evaluation With Expert Users

We present an Evolutionary Visual Exploration (EVE) system that combines visual analytics with stochastic optimisation to aid the exploration of multidimensional datasets characterised by a large number of possible views or projections. Starting from dimensions whose values are automatically calculated by a PCA, an interactive evolutionary algorithm progressively builds (or evolves) non-trivial viewpoints in the form of linear and non-linear dimension combinations, to help users discover new interesting views and relationships in their data. The criteria for evolving new dimensions is not known a priori and are partially specified by the user via an interactive interface: (i) The user selects views with meaningful or interesting visual patterns and provides a satisfaction score. (ii)

The system calibrates a fitness function (optimised by the evolutionary algorithm) to take into account the user input, and then calculates new views. Our method leverages automatic tools to detect interesting visual features and human interpretation to derive meaning, validate the findings and guide the exploration without having to grasp advanced statistical concepts. To validate our method, we built a prototype tool (EvoGraphDice) as an extension of an existing scatterplot matrix inspection tool, and conducted an observational study with five domain experts. Our results show that EvoGraphDice can help users quantify qualitative hypotheses and try out different scenarios to dynamically transform their data. Importantly, it allowed our experts to think laterally, better formulate their research questions and build new hypotheses for further investigation.

Introduction

The purpose of visual exploration is to find meaningful patterns in the data which can lead to insight. In a highdimensionality context, this task becomes rather challenging as viewers may be faced with a large space of alternative views on the data. One way to help navigate such a space is the "grand tour" method [START_REF] Asimov | The grand tour: a tool for viewing multidimensional data[END_REF] which offers a complete view of the search space in a smooth sequence of projections showing various viewpoints of the data. However, the time required to inspect all these views may be prohibitive [Hub85]. A related approach that improves on this is "projection pursuit" [Fri87] where the aim is to visit only the most interesting views; interesting referring to projections that deviate more from a normal distribution. The criteria for deciding whether a projection is interesting have mostly been defined prior to user exploration, using objective measures such as the quality metrics surveyed in [BTK11].

We present a novel visual analysis tool to explore multidimensional datasets where the system proposes interesting views based on both objective measures, such as different vi-sual patterns in the two-dimensional projections of the data, and subjective measures corresponding to user satisfaction with the presented view. These subjective measures are not known prior to user exploration. To demonstrate our ideas, we built a prototype (EvoGraphDice) as an extension of an existing scatterplot matrix inspection tool. We use low dimension projection to handle data multi-dimensionality, and linear and non-linear combinations of dimensions for an axis of the projection plane to propose alternative views. User exploration is guided by an Interactive Evolutionary Algorithm (IEA) which can both generate new views and adapt to user interest. Below, we provide background for the topic of evolutionary computation before listing our contributions.

Evolutionary Algorithms (EAs) are stochastic optimisation heuristics that copy, in a very abstract manner, the principles of natural evolution that let a population of individuals be adapted to its environment [Gol89]. They have the major advantage over other optimisation techniques of making only few assumptions on the function to be optimised. In short, an EA considers populations of potential so- lutions exactly like a natural population of individuals that live, fight, and reproduce, but the natural environment pressure is replaced by an "optimisation" pressure. In this way, individuals that reproduce are the best ones with respect to the problem to be solved. Reproduction consists of generating new solutions via variation schemes (the genetic operators), that, by analogy with nature, are called mutation if they involve one individual, or crossover if they involve two parent solutions. A fitness function, computed for each individual, is optimised by the EA. Evolutionary optimisation techniques are particularly efficient to address complex problems (irregular, discontinuous) where classical deterministic methods fail [Ban97, PLM08], but they can also deal with varying environments [JB05], or non computable quantities [Tak08]. More specifically, Interactive Evolutionary Algorithms (IEAs) are focussed on the optimisation of subjective quantities captured via a user interface.

Evolutionary Visual Exploration (EVE): we feel that

Interactive Evolutionary Algorithms (IEA) are convenient for guiding the user in exploring complex datasets. This opinion is founded by the following characteristics of EAs: (i) focus: an IEA performs an optimisation, i.e. it drives the exploration towards "interesting" areas of the search space (areas of high fitness function and good user satisfaction), (ii) diversity: by nature, an IEA has a stochastic behaviour, and its population-based scheme allows to display a variety of solutions to the user at any time, (iii) adaptation: EAs are able to deal with time varying environments and are able to follow changes of user interest and focus [Lut06].

The contributions of this paper are: (1) a framework for Evolutionary Visual Exploration (EVE) that marries techniques from visual analysis and evolutionary computation to guide user exploration towards interesting views on the data; (2) a prototype tool (EvoGraphDice) [Evo] to demonstrate our framework; and (3) an observational study with five domain expert users to evaluate EvoGraphDice.

Related Work

Related work is organised as follows; (1) a brief overview of quality metrics used to describe specific properties of data projections; (2) description of quality metrics we use in this work as part of the automatic evaluation of scatterplots; and

(3) a summary of related work to IEA.

Quality Metrics: faced with the overwhelming possibilities of exploration paths in multidimensional visualization, researchers in the field designed quality metrics that evaluate the various projections of the data, in the hope of focusing user search on the most promising views. In a recent survey, Bertini et al. [BTK11] used the data flow model to classify quality metrics into three types: metrics that draw information from the data space, from the image space or from both.

Amongst metrics calculated at the data space are clustering and outliers. The rank-by-feature framework [SS05], for instance, visualises an optimal set of features according to a user selected quality metric such as correlation or uniformity. They use axis-parallel projections to produce 1D or 2D views and color brightness to denote ranking scores. Amongst image based metrics are scagnostics [WW08] which describe measures of interest for pairs of dimensions based on their geometrical appearance on a scatterplot. The mixed metrics combine information from the data and image spaces at the same time. Peng et al. [START_REF] Peng | Clutter reduction in multi-dimensional data visualization using dimension reordering[END_REF], for example, combine data features such as correlation information with view features such as axes adjacency to measure clutter as a result of reordering visualization axes [BTK11].

When interaction with quality metrics is available, it is either to select a metric amongst others, or to set threshold values [BTK11]. Having to specify the type of ranking criteria requires users to be familiar with advanced statistical concepts. In our case, the quality metric is pre-defined as a vector of nine image-based measures (described in the next section) and the threshold values are adapted according to user feedback. Thus, our method leverages automatic tools to detect interesting features and human interpretation to derive meaning, validate the findings and guide the exploration without having to grasp advanced statistical concepts.

Scagnostics † are based on geometric graphs which are calculated from areas, perimeters and lengths of these graphs. They include nine measures to characterise scatterplots (Fig. 2) and are useful for quickly discovering regularities and anomalies in scatterplot matrices. The underlying algorithm detects different types of point distributions including multivariate normal, log normal, multinomial, sparse, dense, convex and clusters. It does so by binning, detecting outliers and computing measures based on the following three statistical properties: shape for convex, skinny and stringy distributions; trend for monotonic distributions; and † Available as a free downloadable package in R from http:// www.rforge.net/scagnostics/ Visualization and IEA: visualization tools have been used in IEA both as representation and exploration tools to help users better evaluate the output of interactive evolutionary algorithms [HT00, LSA * 06]. Despite efforts to design good user interfaces for IEA, human interaction with these systems usually raises several problems, mainly linked to the "user bottleneck" [PC97], human fatigue and slowness. Various solutions have been considered [PC97, Tak98, Ban97] such as reducing the population size (micro-EAs), constraining the search space to focus on a priori "interesting" areas, and deploying approximated user models (also called surrogate functions) to filter obvious bad solutions [LPLV05]. In the visualization community, work on parameter space exploration and optimisation relates to ours. Matkovic et al. [START_REF] Matkovic | Interactive visual analysis supporting design, tuning, and optimization of diesel engine injection[END_REF], for instance, tried to interactively find an optimal combination of input parameters for a complex diesel engine injection system using visual analysis techniques. However, to our knowledge, we are the first to propose using IEA as optimisation tools to help navigate large search spaces.

EvoGraphDice

Since our main contribution in this work does not lie in a novel visualization system, but in enabling an IEA to guide user exploration, we used an existing visualization tool (GraphDice [EDF08, BCD * 10]) to manage the various projections of the data. Views are organised in a scatterplot matrix (SPLOM) of 2D projections, Fig. 1(a). Users can do brushing and linking using a lasso tool. EvoGraphDice displays the dimensions proposed by the IEA as additional rows (and columns) in the SPLOM. The system initially displays dimensions returned by a PCA, after which the user can evolve new dimensions by pressing the "evolve" button, Fig. 1(d). The proposed views are displayed in yellow background; the darker the color the more interesting the view. The system provides an initial score (1 to 5) for each new view but the user can adapt this score using the slider in Fig. 1(d). User evaluated cells are flagged (small black square) to distinguish them from system evaluated N. Boukhelifa et al. / Evolutionary Visual Exploration cells. EvoGraphDice can be initialised at any time using the "restart" button which resets parameters of the IEA. Users can save views (Fig. 1(f)) and bring them back into the SPLOM if they have been replaced during the exploration.

The current population is also displayed as a table (Fig. 1(h)) where each row corresponds to a combined dimension described by a mathematical expression and various components of the fitness function such as the scagnostics measures. The user can edit an individual using the "dimension editor" in Fig. 1(j), and limit the dimension search space Fig. 1(i), which results in a system reset similar to precessing the "restart" button. Note that many EA parameters can be tuned, such as the fitness threshold and crossover/mutation/replacement rates (see [CBL12]).

Our prototype has been developed from a first version [CBL12] based on an IEA that only manipulated linear combinations of dimensions. Our new extensions are: (i) a Genetic Programming (GP) algorithm allowing the manipulation of non-linear combinations of dimensions as variable size mathematical formula, (ii) user assessment of proposed views is explicitly captured via a slider, (iii) a surrogate function based on scagnostics measurements is used to predict and simplify the interactions of the user with the IEA, (iv) color highlighting of cells is used to draw user attention to the most interesting views.

Search Space: The space searched by the evolutionary process is the set of all dimensions that can be built by combining the initial dimensions with operators and constants, encoded as trees according to the Genetic Programming (GP) framework [Koz92]. These combinations can be complex mathematical expressions containing quadratic, exponential or logarithmic terms (evolved expressions can be any combination using +, -, * , /, (.) (.) , exp and log operators).

Genetic Engine:

We have chosen to evolve a small set of combined dimensions, in order to let the user see all individuals of the population at a glance: if n is the number of initial dimensions, a population of another n combined dimensions is evolved. At each iteration, that is each time the user clicks on the "evolve" button, a new generation is produced by application of selection/crossover/mutation operators and then presented to the user whose judgment (evaluation) is explicitly collected via a slider.

Initialisation: A set of a priori interesting dimensions has been chosen as a starting point. A PCA analysis is performed [Smi02] on the original data and the corresponding n linear combinations form the initial population.

The fitness function, that is optimised by the genetic engine, is a sum of three terms:

1. A surrogate function fsc, that plays the role of a predictor, and helps the system to better adapt to user needs. It is based on scagnostics measurements computed for every cell of each dimension y i , the corresponding fitness term is a linear combination of the highest values of the scagnostics (SC k (y i , x j )) of each scatterplot cell (y i , x j ):

fsc(y i ) = ∑ k=1..9 w k (max j SC k (y i , x j )) (1)
The weights w k that govern the relative importance of each scagnostic measurement are initialised to a uniform weight (1/9). Then, as soon as enough interactions are recorded (n, the number of variables), w k are updated via a simple multilinear regression on the m past interactions (m ≥ n corresponds to the length of the "memory" of the system). 2. A Complexity term that favours dimensions made of a small number of variables and simple mathematical expressions :

fc(y i ) = 1 - nvars(y i ) n × 1 depth(y i ) , (2)
nvars(y i ) is the number of original variables involved in the mathematical expression of y i , and depth(y i ) is the depth of the GP tree representing y i . 3. A user evaluation term, fu(y i ), that is an average of the user evaluation for each cell corresponding to y i (range of 1 to 5 from "bad" to "excellent").

Diversity management :

The evolutionary mechanisms naturally tend to concentrate the population around good solutions. So for small populations sizes, there is a risk of premature convergence if no diversity preservation mechanism exists. In EvoGraphDice, each time a new dimension y i is generated, its Euclidean distance to the current population is computed. If y i is too close to one of the individuals of the current population, it is replaced by a random individual.

Case Studies with Expert Users

We conducted an observational study with five domain experts. During the study sessions, we encouraged participants to think-aloud and share their findings with the study facilitator. We wrote observations, conducted semi-structured interviews and questionnaires, video-recorded the sessions and logged user interactions. The following sections describe the study setup, observations and findings for each expert.

Method

Due to the open-ended style of exploration using Evo-GraphDice, and the subjective nature of our fitness function, we chose a qualitative observational study methodology [Car08, MSM12, SMM12] that better suits our evaluation needs. We wanted to evaluate the usability and utility of our tool. In particular, we attempted to answer the following three questions: (i) is our tool understandable and can it be learnt; (ii) are experts able to confirm known insight in their data; and (iii) are experts able to evolve views that contain new insight or allow them to generate a new hypothesis, and if so how easy or difficult is it to reach those findings.

Participants and Apparatus

We evaluated our prototype with 5 domain expert users (2 female), ages 27 -42 (mean 34.2). Experts were academics and practitioners who had multidimensional datasets related to their domain of expertise (scientific simulation, medicine and geography) and were interested in further exploration. They consisted of one graduate student, four senior researchers and one medical surgeon. Participants had previously explored their datasets using graphical tools (e.g. Excel and JMP) or used statistical methods (PCA and regression analysis) but felt there is more to discover in their data than their current tools allowed them to. Experience with advanced multidimensional visualization tools varied from none, to experts who already used GraphDice or other SPLOM-based tools (two experts). None of our participants previously used dimension combination to analyse their data but three performed PCA-type analysis. The first three case studies ran at our research lab on an HP Z800 workstation PC with a 19 dual monitor (1280 x 1024 screen resolution). The last two case studies ran on a similar setup at the experts' institutions. Each session lasted on average 2.5 hours.

Tasks and Procedure

Participants were asked to carry out two main tasks: (T1) show in the tool what they already know about their data, hypothesis and questions they wanted answered; and (T2) explore their data in light of these hypotheses and research questions. The first task (and a training game) was designed to test if the tool is understandable, easy to learn, and can help experts rediscover known findings. The second openended one explores how domain experts use our tool to answer questions about their data and gain new insights.

Prior to the actual study, participants filled by e-mail a prequestionnaire to elicit their background, knowledge about the dataset they want to explore, and experience with multidimensional data visualizations. In particular, they were asked to describe the dimensions of the data sets they provided, known relationships between variables, and hypotheses they wanted to investigate. The main study ran in two parts; first training then open exploration as follows:

Training: participants played a game designed to teach them how to operate the tool. A 5D dataset was synthesised with two enclosed curvilinear dependencies between two variables (x0 and x1) and random data for the rest of the dimensions. Participants were asked to evolve a scatterplot where it is possible to separate the two curves in Fig. 3 (left) with a straight line and were given around 20 minutes to complete the task (this task is equivalent to separating the two convex hulls in Fig. 3). Two participants successfully separated the two curves, while the remaining experts evolved views very close to a correct solution within the allocated time.

Open Exploration: the second part of the study ended after about one hour of exploration (a maximum limit of two hours was set based on a pilot to avoid user fatigue), and participants were encouraged to take breaks. A facilitator was present to answer experts' questions and discuss their findings. Throughout the study, a second screen with an open text editor and pen and paper were provided to the experts as means of writing down their exploration findings. At the end, participants filled in a short questionnaire rating aspects of the tool (5-point Likert scales), such as the ease of performing the two main tasks, and open ended questions regarding their exploration strategy and helpful features of the tool.

Data Collection and Analysis

Participants were video-taped and log data of user interactions was gathered for further analysis (table1). Live and video observations, the results of the questionnaire, and the log analysis are described separately for each case study.

Expert 1: Electrical Consumption Profiles

Dataset: (9D) describes the electrical consumption of 900 anonymised businesses during non-peak (npk) and non-plan (npn) hours ('plan' refers to an agreed unit rate for a defined period of time) for winter (W) and summer seasons (S), their geographical altitude and the total consumption cost.

Goal: the expert wanted to investigate electricity consumption patterns of these businesses and their impact on the total cost of consumption. The expert noted in the prequestionnaire that he would like to sum-up some dimensions in twos in order to focus on one aspect of consumption (e.g. non-plan and non-peak for summer), and therefore had a clear motivation for combining dimensions. This, he argued, may allow him to see interesting consumption profiles.

Observations: the expert hypothesised that altitude has an influence on electricity consumption during both summer and winter seasons. He also had some prior knowledge about existing outliers in the data. During the study he was able to quickly verify both of these hypothesis.

The most important new finding made by the expert, which was not part of the original search space, is a view showing a linear combination of the four parameters of interest to the expert (npn and npk for summer and winter) which brought to evidence in a quantitative manner that npnW consumption is the more correlated to the total consumption. 1: Log data showing: (G) the generation when a solution for the game was found, (T1&T2) experts' scores for ease of completing tasks T1&T2 on a 5-point Likert scale, 5 signifies "very easy", (Q) score for user agreement with EvoGraphDice cell evaluations on 5-point Likert scale, 5 indicates strong agreement, (Data&Size) type and size of dataset, (D) duration (hh:mm) of T2, (LimitSearch) breath of exploration indicated by the number of times the expert limited the search space, (Evolve) depth of exploration indicated by the maximum reached generation, (Eval) how many new cells were evaluated by the user, (OVisits&NVisits) number of times the expert visited the original cells and the new cells respectively, and (Insight) number of times the expert limited the search space and the generation (between parenthesis) where the insight was found.

In the user's own words: "we always talk about this qualitatively. This is the first time I see concrete weights ... To understand what is a better fare, it is necessary to find a good approximation of the consumption profile", like the one found in Fig. 4. Thus this insight can lead to electricity plans that are better suited to clients' consumption needs. According to this participant, his exploration strategy was to look at propositions in detail along a row, e.g., to examine proposed dimensions plotted against the total consumption. Overall, the expert did not evolve many generations (depth of exploration was three generations at most), but used "limit the search space" facility three times, indicating that he was trying to formulate an interesting hypothesis more than he investigated one in depth. The solution he found was after limiting the search space for the second time.

The expert liked the ability to limit the search space and to enter formulae for the combined dimensions using the dimension editor, e.g. to invert a weight.

Expert 2: Biscuit Baking Process

Dataset: (7D) describes 78 data points recorded from several industrial biscuit training processes taken by experts in the industry. There are two input parameters relating to temperature settings and three output parameters relating to biscuits (weight loss, height and colour) and a timestep.

Goal: the expert wanted to visualise dependencies in the data between input and output parameters (intuitively such correlation should exist but its exact nature is not clear).

Observations: the expert was able to quickly verify known profiles in the data, for instance the influence of temperature on height and color of the biscuit.

The more general profile of the relationships between input and output parameters was not evident from the original dimensions, thus the expert looked at a wider space using combined dimensions. He observed that there might be some exponential factors that link outputs and inputs, in particular an exponential dimension of one of the input parameters (proposed by the GP) was linearly correlated between all output variables: "we would probably not have considered looking at exponential relationships" indicating a surprising finding and, thus, the ability of the tool to encourage lateral thinking. Further investigation showed that the exponential of temperatures has a specific meaning in thermally activated processes (explained by the Arrhenius law [Wik]).

To look for new relationships, the expert's strategy was to evolve a few generations and choose a visualization that showed linear or quadratic relationships. Like the first expert, he edited promising views using the "dimension editor" to see if this had better or worse impact on the relationship between variables. Most importantly he tried to reduce formula complexity to make better sense of the relationship.

The expert found the evolution functionality and the preview matrix useful to reach interesting cells. However, he wanted to examine scores per dimension as well as per views, e.g. to see if a dimension is always highly ranked; and for the tool to highlight new dimensions that have not been visited before.

Expert 3: Anatomical Planning for Surgery

Dataset: (12D) of which half describe anatomical and geometrical values related to a 3D planning of a surgical operation (total hip arthroplasty) for 67 patients, and the other half represent values for the same parameters after surgery.

Goal: the expert wanted to investigate whether there is a correlation between the planned values and the final values for each of the investigated parameters, and, if it exists, how strong is this correlation. Since there are many parameters to examine with potentially many interactions between them, the expert focussed on examining the offsets for the cup anteversion parameter (AntvCupSupine) which corresponds to the orientation angle of the cup of the hip prosthesis.

Observations: the expert already knew that there is a relationship between the planned and real values for the AntvCupSupine parameter. This was easily verified in the original dimension space. To explore this further, the expert examined the view showing the before and after values, then made three lasso selections corresponding to over-fit, bestfit and under-fit values (Fig. 5 left), then examined brushed cells in the original search space.

In terms of new insight, the expert found a new cell where the two problematic groups (in red and blue) were separated from the well-restored group (in green) with the exception of one data point. The proposed dimension had a simple formula that involved two original dimensions. Views showing such separation may correspond to special geometrical settings or anatomical features for the observed patients. The expert noted that he needs to examine these patients more carefully with special attention to the selected parameters which can then lead to better pre and post-surgery results.

The expert followed the training game example as his exploration strategy, which may explain the big depth of exploration (21 generations): he made lasso selections of groups of data points, and evolved views that he scored highly depending on whether the overlap between the clusters is minimised. He examined the proposed dimensions in relation to the AntvCupSupine parameter and made use of the "favorite" facility to compare interesting cells that were replaced in the next generations. Notably, he evaluated more than 26% of new visited cells. This seemed to be an important part of his exploration strategy. The expert commented that he liked the direct visual interaction with the data but he did not like the uncertainty in whether a solution existed and whether the tool will find it. For example, he was interested to see if the degree of separation between his data groups became smaller between evolutions. He suggested adding more adapted tools for selecting data clusters and including statistical information.

Expert 4: Pareto Front Exploration

Dataset: (10D) describes the outpt of a genetic algorithm that was used to calibrate a city growth and emergence model. The data represents a set of parameter values (7 dimensions) and their objective fitness scores (3 dimensions). The explored dataset only includes the first best 200 parameter values that the algorithm found according to the three objectives of the calibration model (i.e. the Pareto front of the global parameter space exploration).

Goal: the expert wanted to explore the dataset from the two different perspectives (parameter and objective space) as well as the interaction between the two spaces, e.g. does a special profile in the parameter space correspond to a special profile in the objective space?

Observations: prior to the study, the expert had an idea about some characteristics of the data, e.g. there are two large clusters that can be differentiated by the value of one parameter (pAdoption). This type of calibration was also known to produce a characterisable response in the objective space. This hypothesis was easily verified using Evo-GraphDice via brushing and linking between cells in the parameter space and the objective space. In terms of new insight, the expert was able to find an interesting combined dimension that gave a good correlation for two parameters of the objective space (Fig. 6). The expert commented that this combination may be an important finding because it involves parameters that affect only one part of the simulation model. This indicates that those parameters, at least for these two output indicators, work together; and that this linear combination could be one way to reduce the complexity of the model.

As for strategies, the expert mentioned primarily limiting the search space (7 times), evolving (13 generations) and examining cells that had monotonic or striated distributions. She also made good use of visual queries and cell evaluation; 27% of new visited cells were evaluated by the participant.

The expert liked the ability to limit the search space and the freedom the tool offers to explore the data by "evolution". She added that the tool was helpful in reaching insight because of its ability to visualise and suggest combinations of dimensions that actually had a visual pattern: "in Excel it is difficult to find a formula that would give a nice pattern".

Expert 5: Urban Organisation and Perception

Dataset: (11D) describes geo-spatial information about 653 inhabitants and their urban environments including their profiles (e.g. gender, age), perception of their neighbourhood and objective variables describing their street such as the distance to the nearest metro station and type of district.

Goal: as with other experts, this participant was interested in finding relationships between variables; specifically she wanted to identify groups of inhabitants having similar profiles and to find the most discriminating variables for these individuals in order to make sense of the formed groups.

Observations: the expert had already found interesting correlations between the different original variables using her own statistical tools. She was able to confirm these findings. For instance, that an individual's perception of the size of their neighbourhood was dependent on the distance to the next metro station, but this was only true if context (type of district) was taken into account. However, the expert was aware of correlations requiring an interaction between two variables against a third but found it difficult to see them using EvoGraphDice. The expert noted two major difficulties that may have hindered the exploration and thus lack of early insight or hypothesis generation: (a) difficulty in determining the criteria for scoring patterns without knowing what a good pattern is in advance; and (b) the nature of data about human behaviour and perception has high variability, thus examining averages, for instance, is more appropriate.

Despite the aforementioned difficulties, the expert found two interesting views where clusters and outlier groups seem to correspond to a known profile (Fig. 7). However, she was not able to fully interpret the proposed combined dimensions as the choice of variables made sense, but the overall interpretation of the pattern was not clear to the participant. This expert's exploration strategy was to limit the search space to 3 -4 variables and examine their interaction with one original dimension (e.g. perception of space). She also made selections and examined the brushed views in the original space. Since the expert did not evolve many generation (5 max) and only evaluated a few cells (20 overall)-due to the aforementioned difficulties-the system did not learn well the type of distributions the expert was looking for.

The expert tended to agree with the system's proposed scores (Table 1 Q), which she found interesting because of the choice of variables and the simplicity of the proposed formulae. As interpretation of results was difficult using the current point-based presentation, the expert noted that showing aggregated values and variance would help her better understand the views.

Summary of Results

Almost all participants were able to easily confirm prior knowledge about their datasets (2 x 'very easy', 2 x 'easy', 1 x 'neutral'). One expert found this task challenging because of the lack of data aggregation that her type of analysis requires. Overall, participants confirmed known correlation, clusters or outliers in their data. In the remainder of this section, we summarise our study findings concerning new found insight, successful tasks and exploration strategies.

Insight Generation and Tasks

If we include hypothesis formation as part of insight generation, similar to work by Saraiya et al. [SND05], Evo-GraphDice helped our participants generate new insight in the form of distinct observations about the data (4 experts), new hypothesis (1 experts) and better formulation of research questions (4 experts). Distinct observations found by the experts were either clustering, linear or non-linear relationships, and similarly to generated hypotheses, they always linked a dimension in the original data set and a new proposed dimension. The subjective evaluation of ease of task T2 (table 1) shows most experts found it easy to find new insight: 1 x 'very easy', 3 x 'easy' and 1 x 'not easy'. Not surprisingly, those who reached a concrete new finding scored the tool highly in comparison to those who did not.

The found solutions were regarded by the experts as interesting because they had one or more of the following properties: (i) a visual pattern such as those modeled by the scagnostics measures; (ii) a simple formula involving few dimensions; (iii) a selective choice of dimensions (corresponding to an unformulated hypothesis or an inherent aspect of their data model); and (iv) a domain value. Regarding the latter point, not all participants were able to state the immediate domain value but in general, our participants stated that EvoGraphDice helped them:

• interact visually with data (experts 3) • try out alternative scenarios by editing dimensions (experts 1,2) • think laterally (expert 2) • quantify a qualitative hypothesis (expert 1) • formulate a new hypothesis or refine an existing one (1-4)

Exploration Strategies

Overall, participants followed the same exploration pattern consisting of first examining the original dimensions then inspecting and evaluating the first generation of the proposed dimensions (returned by the PCA) followed by one or more iterations of the following steps: (i) limit the search space; (ii) select and rank cells; (iii) evolve; and (iv) interpret and verify. However, the frequency of using some tools (e.g. "evolve" vs. "limit the search space") varied depending on whether the expert had an a priori focused hypothesis (i.e. a research question involving typically 3 -4 dimensions). We observed that the looser the initial hypothesis, the more often they tried to change the search space; and the more focused the hypothesis the more generations they inspected. Indeed, these two strategies of exploration and exploitation are supported by EAs [Ban97] where on the one hand the user wants to visit new regions of the search space and on the other hand they want to explore solutions (combined dimensions) close to one region of the search space.

Discussion

Most of our experts were able to formulate interesting hypothesis or reach new insight requiring looking at data in terms of a combination of dimensions. Our approach consists of proposing new views based on automatically calculated metrics and user feedback. On the one hand, our method is complimentary to PCA, clustering and regression analysis that automatically find data patterns, and optimise a fit. On the other hand, we allow users to interactively select examples of visual pattern types they are interested in, and that may not be easy to express mathematically. Users can then verify the new relationships they find in EvoGraphDice by using the dedicated automatic data analysis tools.

In comparison to automatic analysis such as in statistics and data mining, our approach offers: (i) Intuitiveness: a visual approach to interact with data requiring no prior statistical knowledge; (ii) Interactivity: rather than fitting the data to pre-defined shapes in a static manner, using an IEA the user can dynamically steer the exploration process towards a pattern of interest. These patterns can involve dimension concatenations that are not obvious at the outset of the exploration; (iii) Flexibility: ability to edit and try out alternative dimension combination scenarios, or limit the search space. (iv) and Adaptability: the system can adjust to user change of interest over time.

There are limitations to using our tool, such as the types of datasets to explore and issues related to the interpretation of combined dimensions which we discuss below. Moreover, the issue of the convergence of the genetic algorithm is an interesting one given that the IEA deals with optimisation. However, this is not easy to study in our case as there is no unique solution to converge to, rather the optimisation is dynamically adapted to follow user interest over time.

First, we are constrained by the SPLOM representation of EvoGraphDice which does not provide a natural way to interact with some dataset types such as timeseries. Data with high variability provides additional challenges that we do not currently address, such as detecting and evolving aggregated patterns. In addition, we tested our prototype with user-provided datasets that are small to medium sized, having dimensions between 7 -12. Although our algorithm can deal with a large number of data points, it may not handle well larger number of dimensions as complex combined dimensions may be difficult to avoid. In this case, a dimension reduction technique can be applied to the dataset before feeding the results to EvoGraphDice.

Second, not all variables can be combined, therefore the user should as soon as possible limit the search space to "combinable" dimensions. This in a sense requires the user to have some domain knowledge and to make an initial hypothesis about the data. The proposed dimensions can involve complex or unforeseen combinations yielding a visual pattern but one that can be difficult to interpret. To help address this issue, we used "complexity" of a combined dimension as a component of the IEA fitness function. Nonetheless, our method can still yield complex dimensions that are difficult to interpret. We noticed that our participants only looked at the combined dimensions in relation to an original dimension, most likely to ground the observation. This problem of interpretation, however, is common to all tools that offer dimension combination.

Conclusion and Future Work

We presented a prototype tool (EvoGraphDice) for supporting Evolutionary Visual Exploration (EVE) that combines visual analysis with interactive evolutionary computation to help steer the exploration towards interesting views on the data. Our method complements PCA, clustering and regression types of analysis, offering additional features such as interactivity and adaptability. We conducted an observational study with domain experts and found that our tool allowed users to evolve characteristics that are not visible in the original dimensions space. Our experts were able to try out different scenarios, think laterally, quantify qualitative hypotheses and formulate new ones.

Future work for our tool includes longitudinal user studies to explore in detail the long term evolution of user focus, as well as addressing issues such as improving the IEA to detect more complex visual patterns (beyond those currently detected by Scagnostics), handling data with highdimensionality and bridging EvoGraphDice and existing statistical packages to combine powerful statistical analysis with flexible and intuitive visual exploration.

Our work demonstrated that tightly combining visualization and optimisation techniques can yield exciting results in data analysis and opens new venues for research, but also highlights challenges such as monitoring algorithm convergence, history visualization of diverging exploration paths, and appropriate methodologies for evaluation.

Figure 1 :

 1 Figure 1: EvoGraphDice prototype showing an exploration session of a synthetic dataset. New extensions to the GraphDice system are indicated by coloured label arrows. Widgets: (a) an overview scatterplot matrix showing the original data set of 5 dimensions (x0..x4) and the new dimensions (1..5) as suggested by the evolutionary algorithm. (b) main plot view. (c) tool bar for main plot view. (d) a tool bar with (top to bottom)"favorite" toggle button, "evolve" button , a slider to evaluate cells and a restart (PCA) button. (e) the selection history tool. (f) the favorite cells window. (g) the selection query window. (h) IEA main control window. (i) window to limit the search space. (j) dimension editor.

Figure 2 :

 2 Figure 2: Nine scagnostics measures from [WW08].

Figure 3 :

 3 Figure 3: Two different solutions (screenshots of plots) for the training game problem (left) that involve a simple dimension combination (middle) and a complex formula (right).

Figure 4 :

 4 Figure 4: Confirmed findings (left and centre) and new insight found by the expert (right): a linear combination of four parameters that approximates customer consumption.

Figure 5 :

 5 Figure 5: Selections of over-fit (red), best-fit (green) and under-fit (blue) parameter values (left), and (right) a finding by the expert showing a separation between the two groups of interest in relation to a new parameter.

Figure 6 :

 6 Figure 6: An interesting combined dimension from the parameter space and its impact on two objective dimensions.

Figure 7 :

 7 Figure 7: Two interesting combined dimensions (centre and right) found by the system and their impact on one objective dimension (aireha). Brushing and linking to an original view (left) shows interesting profiles.

Table

						N. Boukhelifa et al. / Evolutionary Visual Exploration
	Expert	G	T1	T2	Q	Data	Size	D	LimitSearch	Evolve	Eval	OVisits	NVisits	Insight
	1	-	4	4	3	business	9x900	1:10	3	3	16	40	105	2(1)
	2	-	4	4	3	timeseries	7x78	1:33	4	3	8	114	115	4(3)
	3	9	5	5	3	geometrical	12x67	0:49	4	21	90	99	344	2(1)
	4	7	5	4	3	statistical	10x200	2:23	7	13	83	110	309	6(1)
	5	-	3	2	4	geospatial	11x653	1:27	5	5	20	64	229	-

c 2013 The Author(s) c 2013 The Eurographics Association and Blackwell Publishing Ltd.

Acknowledgements:

We thank all our expert users for their participation in this study, and J-D Fekete for his useful feedback and comments.


Simon Labarthe

Edward Vigmond

Yves Coudière

Jacques Henry

Hubert Cochet

Pierre Jaïs

A Computational Bilayer Surface Model of Human Atria

A bilayer surface model of human atria is presented. A rule-based bi-layer physiological fibre arrangement is defined on an imaged in-vivo atrial geometry. This fibre architecture includes the main structures of fibres in the right and left atria, including the main transmural heterogeneities and transseptal connexions.

The precision of the corresponding bilayer mathematical model is assessed by comparing the case of two sheets of orthogonal fibres to a three dimensional slab of tissue. A comparison of bilayer and monolayer simulations of sinus propagation is finally proposed.

This model allows inclusion of transmural heterogeneities while maintaining small computational costs associated with surface models. It then proposes a good trade-off between model precision and computing efforts for performing complex atrial simulations for a clinical use.

Introduction

Atrial computational models have proven their utility in fundamental or clinical research by giving insight into atrial fibrillation (AF) initiation and 1 perpetuation, radio-frequency ablation, and the development of signal analysis techniques [START_REF] Dössel | Computational modeling of the human atrial anatomy and electrophysiology[END_REF].

As the electrical activity spreads across the atria, propagation patterns and P-wave duration are strongly influenced by atrial fibre structure [START_REF] Krueger | Modeling atrial fiber orientation in patient-specific geometries: A semi-automatic rule-based approach[END_REF][START_REF] Zhao | An image-based model of atrial muscular architecture / clinical perspective[END_REF]. Modeling the anisotropy of the atria is, therefore, a major step in the construction of atrial models. This construction can result from an imagebased process on ex-vivo hearts [START_REF] Zhao | An image-based model of atrial muscular architecture / clinical perspective[END_REF] or from a semi-automatic rule-based approach for acquired geometries of in-vivo atria [START_REF] Krueger | Modeling atrial fiber orientation in patient-specific geometries: A semi-automatic rule-based approach[END_REF][START_REF] Labarthe | A Semi-Automatic Method To Construct Atrial Fibre Structures: a Tool for Atrial Simulations[END_REF].

Atrial three dimensional models are usually divided into two main categories: surface and volumetric models. Whereas the former propose a homogeneous transmural fibre distribution and computationally light simulations, the latter capture important transmural heterogeneities or layer dissociations as observed in animal models [START_REF] Derakhchan | Method for simultaneous epicardial and endocardial mapping of in vivo canine heart: Application to atrial conduction properties and arrhythmia mechanisms[END_REF][START_REF] Eckstein | Time course and mechanisms of endo-epicardial electrical dissociation during atrial fibrillation in the goat[END_REF][START_REF] Michowitz | Electrophysiological differences between the epicardium and the endocardium of the left atrium[END_REF] in exchange for increased computational costs [START_REF] Dössel | Computational modeling of the human atrial anatomy and electrophysiology[END_REF].

As a compromise between the surface two dimensional and volumetric three dimensional models, we propose a two layer computational model of the atria that can reproduce complex three dimensional patterns while keeping surface computational load. Such a two layer approach has also been used in [START_REF] Gharaviri | A computer model of endo-epicardial electrical dissociation and transmural conduction during atrial fibrillation[END_REF] on idealized geometries, but has not been demonstrated to reproduce proper behavior.

We first describe the construction of the model. We then give some mathematical arguments to the validity of the bilayer approach. We finally compare simulations of a normal action potential (AP) propagation in a bilayer and a monolayer model.

Method

Imaging the Atria

Bi-atrial geometry was acquired in vivo using ECG-gated contrast-enhanced multi-detector computed tomography (MDCT) in a patient with no history of atrial disorder. Imaging was performed using a 64-slice CT scanner (SO-MATOM Definition, Siemens Medical Solutions, Forchheim, Germany) during intravenous injection of a 120 mL bolus of iodinated contrast agent at the rate of 4 mL s -1 . ECG-gating was set for the acquisition window to occur at end-systole, when atrial motion is minimal. Acquisition parameters were: slice thickness 0.5 mm, tube Voltage 120 kV, maximum tube current 850 mAs, and gantry rotation time 330 ms. An imaging volume comprising the whole atria was reconstructed with a voxel size of 0.5x0.4x0.4 mm. Transaxial images were imported in DICOM format in a local database of the soft-ware OsiriX 3.6.1 (OsiriX foundation, Geneva, Switzerland). The endocardial surface was segmented automatically using region growth segmentation. Volume rendering reconstructions were used to remove non atrial structures from the segmentation. Pulmonary veins were cut several centimetres from their ostia because myocardium is known to extend to the first centimeter beyond the pulmonary veins ostia. The resulting segmentation was used to compute a 3D mesh in vtk format using the software CardioviZ3D (IN-RIA Asclepios, Sophia Antipolis, France). The mesh was smoothed with weighted Laplacian smoothing and optimized with the software freeYams to obtain a mesh suitable for calculation that contained 361483 nodes and 715034 elements with a mean diameter of 0.492 mm.

Construction of the Geometry of the Bilayer Model

Algorithms.

Two methods were used for fibres definition and mesh completion.

Fibre definition. The fibre structure was defined with a rule-based semiautomatic method [START_REF] Labarthe | A Semi-Automatic Method To Construct Atrial Fibre Structures: a Tool for Atrial Simulations[END_REF] that was founded on two atrial characteristics : (1) atrial geometry can be divided in basic cylinder-like structures with circumferential fibre orientation. (2) Fibre arrangement is continuous. The atrial mesh was first manually decomposed in cylinder-like elements. The local circumferential structure was semi-automatically obtained by determining the cylindrical axis as an eigen vector of an inertia tensor. Finally, an inpainting algorithm extended continuously the fibre direction on the whole geometry.

Manual mesh completion. Small structures that were partially acquired, such as transseptal connexions, were manually constructed with the following workflow : the points of insertion of the new structure in the mesh were identified. Plane surfaces connecting those points were defined and meshed with the software Gmsh. The structure was then smoothed by a Laplacian method with the software MeshLab (3D-CoForm Project, Italy).

Main modeling options.

In order to describe the atria as a bilayer structure, we need to define some physiological a-priori modeling options, specially focusing on transmural heterogeneities.

Right Atrium (RA). According to histological descriptions of the atria [START_REF] Siew | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF], the RA can be seen as a monolayer structure, except in the appendage (RAA) and in the roof where the crista terminalis (CT) and the pectinate muscles (PM) bring thickness and fibres inhomogeneities. A first layer modeled the epicardial part of the RA including its main fibre structures : the circular fibres that surrounded the superior veina cava (SVC) crossed the CT and covered the RAA to blend with the circular bundle of the vestibule. A second layer was applied on the location of PM and the CT to include the longitudinal fibre direction along those structures in the endocardial part of the RA. A poor conduction zone was introduced along the left side of the sino-atrial node (SAN) as assessed in [START_REF] Vadim | Structural and functional evidence for discrete exit pathways that connect the canine sinoatrial node and atria[END_REF].

Left Atrium (LA). The left atrium (LA) fibres structure can be described as a two layer geometry. The main transmural heterogeneities are concentrated in the anterior wall and on the pulmonary veins (PV). The circumferential and oblique fibres of the anterior epicardium blend with the Bachmann's bundle whereas the endocardial fibres of the septoatrial bundle run from the vestibule to the dome of the LA. The fibre structure of the PV can show strong heterogeneities with non uniform arrangement and transmural discontinuities of fibres directions [START_REF] Saito | Becker Left atrial myocardial extension onto pulmonary veins in humans: anatomic observations relevant for atrial arrhythmias[END_REF].

Transseptal connexions. The chambers of this atrial model were electrically insulated and connected by the three main interatrial connexions : the BB, the fossa ovalis (FO) and the coronary sinus (CS). The BB was segmented manually as a three-dimensional structure and then reduced to two surfaces : the first one corresponded to the posterior face of the BB and inserted into the septal wall of LA while the second one was constructed as a mean surface of the volumetric BB and blended into the septopulmonary bundle of the LA. The FO and the CS were identifiable after the segmentation step and connected manually.

Modeling the AP Propagation

The monodomain equation was used to describe the propagation of electrical conductivity

div(σ i ∇V m ) = β C m ∂ t V m + I ion (V m , t, ζ) (1)
where V m is the transmembrane voltage, σ i is the conducitivity tensor, β is the ratio of membrane surface to volume in which it is contained, C m is the membrane capacitance (1F cm -2 ), and I ion is the transmembrane current which is a function of voltage, time, and ζ, a set of state variables.

Mathematical assessment of the bilayer model validity.

A bilayer expression of this mathematical model has been asymptotically derived in the limit of vanishing thickness in [START_REF] Bendahmane | An asymptotic two-layers monodomain model of cardiac electrophysiology in the atria[END_REF] from a bilayer volumetric monodomain model. That reads:

β(C m ∂ t V (k) m + I ion (V (k) m , ζ (k) )) = div(σ ′(k) ∇V (k) m ) + (-1) k σ c (V (1)
m -V (2) m ) (2)
in the layer k = 1, 2, where σ c = 3

h 2 σ (1) 3 σ (2) 3 σ (1)
3 +σ

(2) 3

. σ ′(k) and σ

(k)
3 are respectively the diffusion coefficients in the layer tangent plane and in the transmural direction. h is the thickness of the layers.

Three dimensional, bilayer and monolayer finite element P 1 -Lagrange simulations have been performed on, respectively, a 3D slab and its 2D equivalent for various h, with the Beeler-Reuter model and a semi-implicit Euler scheme. The activation maps of the mean potential in the thickness of each simulation have been compared. Corrections of the coupling coefficient σ c were uniformly applied to improve the robustness of the model under physiological thickness, far from the asymptotic regime.

Bilayer atrial model. The ionic model was described by the Courtemanche human atrial action potential model [START_REF] Courtemanche | Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model[END_REF]. The sodium conductance was increased by a factor of four from the published value to achieve realistic propagation speeds. σ i and β were set to (0.3, 0.04) S m -1 and 0.14 m -2 . The structural block zone close to the SAN [START_REF] Vadim | Structural and functional evidence for discrete exit pathways that connect the canine sinoatrial node and atria[END_REF] was implemented with a passive ionic model with a resting level -80 mV and an intracellular conduc- tivity reduced by a factor of 100. The SAN is located between the block zone and the crista terminalis. Accordingly, a small region was stimulated there to simulate the site of first excitation.

The equations were solved using the finite element method. Each surface was composed of linear triangular elements while connections between surfaces were implemented as one dimensional finite elements. 134510 line elements were added to couple the surfaces together. The system was solved using the CARP cardiac simulator [START_REF] Edward | Computational tools for modeling electrical activity in cardiac tissue[END_REF] with a fixed time step of 25 µs using a Crank-Nicolson approach to solve eq. 1. On a desktop workstation with Intel(R) Xeon X5650 CPUs running at 2.67GHz, it took seven seconds to simulate one ms on ten cores, while on a desktop with an NVidia Quadro 4000 (256 cores) it took 18 seconds per millisecond running it on the GPU, achieving a better computing time/cores ratio than CPU calculation.

Results

Fibres architecture

The main fibre structures as described in [START_REF] Siew | Atrial structure and fibres: morphologic bases of atrial conduction[END_REF] are correctly reproduced by our model. In fig. 2, we can see that the differences between the endo and epicardial fibres orientations are mainly localized in the CT and PM in the RA, and in the anterior and posterior wall, together with the PVs in the LA, where the fibres in both layers can be orthogonal.

Simulations

Validity of the bilayer approach. The bilayer mathematical model has been shown to have second-order convergence towards the volumetric model when h -→ 0. For physiological h, a correction can be brought to the coupling value σ c so that the accuracy of the bilayer surface model is preserved: complex depolarisation patterns can be obtained with a relative error of 3% respectively to 3D simulations (figure 3). We refer to [START_REF] Bendahmane | An asymptotic two-layers monodomain model of cardiac electrophysiology in the atria[END_REF] for a deeper study of this model.

Bi-atrial simulations. The propagation of electrical activity showed a normal pathway: it first propagated from SAN to the BB along the CT. LA was first activated by the BB, and then by the FO and the CS. Both atria were totaly activated after 196 ms -fig. 4. The sinus wave ended between the left PVs and in the LAA. The effect of transmural heterogeneities are visible along the PM on the activation map. The comparison of the monolayer and bilayer activation maps showed that the main differences of propagation are in accordance with the map of fibres differences -fig. 2. The bilayer wave is faster in the PM direction, around the mitral valve and on the left PVs, and slower on the right PVs and on the anterior and posterior wall of the LA.

Discussion and conclusions

Clinical uses of atrial modeling or acurate atrial models calibrated for fundamental research need personalized fibres orientations, and then efficient methods of construction. The work-flow presented here is still time-consuming and has to be improved in view of such a systematic use. Other methods that are being developped for 3D models [START_REF] Krueger | Modeling atrial fiber orientation in patient-specific geometries: A semi-automatic rule-based approach[END_REF] could be adapted to bilayer surface models. Statistical methods of registration to transport this model to personalized geometries could be another way to shorten this process [START_REF] Toussaint | In vivo human 3d cardiac fibre architecture: Reconstruction using curvilinear interpolation of diffusion tensor images[END_REF]. Modeling assumptions were set prior to model construction. The bilayer representation of some zones (LA,CT,PM) and the monolayer representation of others (RA) were guided by histological descriptions. If more complex transmural heterogeneities are observed, the model can be locally improved by adding layers or by simulating small entities with a volumetric model.

A normal propagation was obtained. The main qualitative characteristics of an atrial normal AP propagation were correctly reproduced: the transseptal pathway and the endpoint of the sinus wave are in a physiological range. The total activation time is larger than those reported by the literature [START_REF] Harrild | A computer model of normal conduction in the human atria[END_REF] but this can be tuned by further adjusting a combination of sodium channel conductances and tissue conductivities.

The bilayer model demonstrated its capacity to capture some perturbations of action potential propagation induced by transmural heterogeneities; that represents an important improvement over monolayer surface models. We could show that those perturbations were reliable on idealized geometries, indicating that this model could be efficiently used on physiological geometries. The demonstration of the potential of the bilayer approach could be completed with a comparison between a volumetric model and its equivalent bilayer surface model. This bi-atrial bilayer model proved to be a good trade-off between surface and volumetric models; it reproduced transmural heterogeneities that were specific to volumetric models with a small computational load. It is an

Figure 1 :

 1 Figure 1: Left:Inferior view of the atria. Right: Superior view. All distinct surfaces are colored uniquely. SVC-superior vena cava. IVC-inferior vena cava. BB-bachmanns bundle. CS-coronary sinus. BZ-block zone. SAN-sinoatrial node. RSPV-right superior pulmonary vein. RIPV-right inferior pulmonary vein.LSPV-left superior pulmonary vein. LIPV-left inferior pulmonary vein. RAA-right atrial appendage. LAA-left atrial appendage. CT-crist terminalis. PM-pectinate muscles.

Figure 2 :

 2 Figure 2: Top: Anterior view of the epicardial fibres arrangement of both atria and zoom on the left inferior part of the septopulmonary bundle. Bottom:Map of the angle in radians formed by the endo and epicardial fibre direction. Dark blue (negative values): only the epicardial layer is modeled. Light blue: 0 rad. Same fibre direction in the endo and epicardial layer. Red: π/2 rad. Orthogonal directions.

Figure 3 :

 3 Figure 3: Left: Snapshots of transmurally averaged transmembrane potential maps at t = 16ms for a 1mm thick slab of tissue with two layers of orthogonal fibres. Top: Volumetric model. Bottom: bilayer surface model (corrected coefficient). Right: Relative L 1 error on activation map versus the thickness of the tissue in log-log scale for the monolayer (blue, stars), the bilayer (green, x) and the corrected bilayer model (red, +)

Figure 4 :

 4 Figure 4: Left: posterior view of the activation map with 5 ms isolines. Dark blue: 0 ms. Red: 196 ms. The deformation of the front is clearly visible in the PM. Right: posterior view of the difference between activation map of the bilayer (epicardial layer, CT, PM, endocardial layer of LA) and monolayer (epicardial layer) surface models. Colors are ranging from -10 ms (dark blue) to +10ms (red). 2 ms isodifference lines are shown.

Acknowledgements: This work was partially supported by an ANR grant part of "Investissements d'Avenir" program reference ANR-10-IAHU-04. The research leading to these results has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under Grant Agreement HEALTH-F2-2010-261057.


Radu Mateescu
email: radu.mateescu@inria.fr

Gwen Salaün
email: gwen.salaun@inria.fr

PIC2LNT: Model Transformation for Model Checking an Applied Pi-Calculus

Introduction

The π-calculus [START_REF] Milner | A calculus of mobile processes[END_REF] was proposed by Milner, Parrow, and Walker about twenty years ago for describing concurrent systems with mobile communication. The π-calculus is equipped with operational semantics defined in terms of Ltss (Labelled Transition Systems). Although a lot of theoretical results have been achieved on this language (see, e.g., [1, chapter 8] for a survey), only a few verification tools have been designed for analysing π-calculus specifications automatically. The two most famous examples are the Mobility Workbench (Mwb) [START_REF] Victor | The mobility workbench -a tool for the pi-calculus[END_REF] and Jack [START_REF] Ferrari | An automated based verification environment for mobile processes[END_REF], which were developed in the 90s.

Our objective is to provide analysis features for π-calculus specifications by reusing the verification technology already available for value-passing process algebras without mobility. Contrary to existing verification tools for the π-calculus, which rely on specific algorithms and intermediate models, such as Hd-automata [START_REF] Ferrari | An automated based verification environment for mobile processes[END_REF], our approach is based on a novel translation [START_REF] Mateescu | Translating pi-calculus into lotos nt[END_REF] from the finite control fragment of π-calculus to a standard process algebra called Lotos Nt (Lnt for short) [START_REF] Champelovier | Reference manual of the LOTOS NT to LOTOS translator (version 5.1)[END_REF]. Lnt is a value-passing process algebra with imperative programming flavour accepted as input by the Cadp verification toolbox [START_REF] Garavel | CADP 2010: A Toolbox for the Construction and Analysis of Distributed Processes[END_REF]. It supports the specification of data structures (constructed types, pattern-matching, recursive functions) and concurrent processes. Lnt has a user-friendly syntax and a formal operational semantics defined in terms of Ltss. To the best of our knowledge, this is the first π-calculus translation having a standard process algebra as target language.

In this work, we go a step further by extending the original polyadic πcalculus with data-handling features. This results in a general-purpose applied π-calculus, which offers a good level of expressiveness for specifying mobile concurrent systems, and therefore for widening its possible application domains. As language for describing data types and functions, a natural choice was Lnt itself: in this way, the data types and functions used in the π-calculus specification can be directly imported into the Lnt code produced by translation. We generalized our previous translation [START_REF] Mateescu | Translating pi-calculus into lotos nt[END_REF] to handle applied π-calculus specifications, and we automated it in the tool Pic2Lnt 2.0. This enables the analysis of applied π-calculus specifications using all verification tools of Cadp, in particular the Evaluator 4.0 on-the-fly model checker [START_REF] Mateescu | A model checking language for concurrent valuepassing systems[END_REF], which evaluates temporal properties involving channel names and data values.

Applied Pi-Calculus

We designed our applied π-calculus by extending the original polyadic πcalculus [START_REF] Milner | A calculus of mobile processes[END_REF] equipped with the early operational semantics defined in [START_REF] Montanari | Checking bisimilarity for finitary pi-calculus[END_REF]. We consider π-calculus specifications satisfying the finite control property [START_REF] Dam | Model checking mobile processes[END_REF], which amounts to forbid recursive agent calls through parallel composition operators. When the set of channels is bounded, this results in finite-state systems that can be analyzed using existing model checking techniques. We extended the original π-calculus with constructs for manipulating data variables and expressions. Agents can be parameterized by data variables in addition to channel names, and the polyadic communication was extended to handle emission/reception of data values. The guard operator was generalized to handle arbitrary Boolean expressions (in addition to the comparison of channel names), and a new operator was added for declaring and initializing data variables. The replication operator was restricted to a bounded version (in order to satisfy the finite control property), which instantiates n parallel copies of an agent, and therefore enables to describe mobile systems containing a finite amount of dynamic control. Data types and functions are specified in Lnt [START_REF] Champelovier | Reference manual of the LOTOS NT to LOTOS translator (version 5.1)[END_REF] as external modules, which are imported in the applied π-calculus specification. The concrete syntax (which is compatible with Mwb for dataless π-calculus specifications) and semantics of the applied π-calculus are described in [START_REF] Mateescu | Pic2Lnt: A translator from the applied pi-calculus to lotos nt[END_REF].

We present below a code sample to illustrate our applied π-calculus on a load balancing system, which is a networking method to distribute workloads across multiple servers. The specification (Main agent) given below consists of five agents: a client, the load balancer, and three servers. The client corresponds to a possible environment and is used to simulate various scenarios. The load balancer receives new tasks (task) with a private name (com), and then interacts with the three servers to know their current load. To do so, a public channel (e.g., al for the first server) is used for sending the request and receiving the result. The load balancer compares the different loads and forwards the private name originally submitted by the client to the server with the minimum load. A server has three possible behaviours: it can be asked by the load balancer to return its current load; it can receive a request for a new task (reception of a private name from the load balancer and interaction with the client on this private channel to receive the new load); or it can execute part of its work if its total load is greater than zero. We can see with this simple example how data expressions (natural numbers, comparison, addition, etc.) appear as parameters of agents and channels to specify loads and their manipulation.

Client(task ) = (ν com1) task com1 .com1 2 of Nat .(ν com2) task com2 .com2 1 of Nat . (ν com3) task com3 .com3 1 of Nat .(ν com4) task com4 .com4

Translation to LNT

Most of the π-calculus constructs are translated quite straightforwardly into Lnt because of its high level of expressiveness. Nevertheless, we faced some subtle difficulties in obtaining a translation as succinct as possible while still preserving the Lts semantics, i.e., mapping each transition of a π-calculus agent to a transition of the resulting Lnt process. One of the main problems was to encode the binary, unidirectional, and mobile communication of π-calculus into a specification language enabling multi-way and bidirectional communication on static channels. Since mobile communication cannot be described directly using Lnt static channels, we overcome this issue by exploiting the data types and synchronization features of Lnt. We represent π-calculus channel names as values of a Lnt datatype Chan, which defines all the public and private names appearing in the specification. Then, we model channel mobility between π-calculus agents by communicating Chan values along Lnt channels. Binary unidirectional communications and two-among-n synchronizations, which cannot be directly described in Lnt, are encoded by means of dedicated Lnt channels (one for each π-calculus parallel composition operator), on which the sender and receiver are indicated explicitly using process identifiers and placeholders. Communication on a πcalculus channel is translated in Lnt as a choice on all Lnt channels connecting the current agent to its environment. The translation of the original π-calculus to Lnt is detailed in [START_REF] Mateescu | Translating pi-calculus into lotos nt[END_REF].

Tool Support and Verification with CADP

The translation from our applied π-calculus to Lnt has been automated by the translator Pic2Lnt 2.0, implemented using the Syntax+Traian compiler construction technology [START_REF] Garavel | Compiler construction using lotos nt[END_REF]. The tool consists of about 2, 300 lines of Syntax code, 4, 800 lines of Lotos Nt code, and 700 lines of C code 3 .

Figure 1 gives an overview of the complete tool chain. Given a specification in applied π-calculus, possibly containing data types and functions described in Lnt libraries, Pic2Lnt translates it into an equivalent Lnt specification, which is accepted as input by the Cadp tools. The resulting Lnt specification is connected by the Lnt.Open tool (via an intermediate translation to Lotos) to the Open/Caesar environment [START_REF] Garavel | Open/caesar: An open software architecture for verification, simulation, and testing[END_REF], which gives access to all the on-the-fly verification tools of Cadp. The pic2lnt dyn.tnt file (static code) contains external C functions for generating fresh channel names and process identifiers. [START_REF] Mateescu | A model checking language for concurrent valuepassing systems[END_REF], an extension of alternation-free µ-calculus with regular expressions, data-based constructs, and fairness operators. Mcl is suitable for analyzing applied π-calculus specifications, because the properties can involve both the channel names and the data values transmitted. The Lts actions, which carry additional information introduced during the translation to Lnt, are renamed on-the-fly to retrieve the original π-calculus format using a predefined label renaming file.

Going back to the load balancing system specified in Section 2, the Lts of the resulting Lnt specification contains 2, 007 states and 5, 450 transitions. As an example, we can check that this Lts satisfies the Mcl data-based response property below, which states that every time a server has begun an execution, it will eventually exhaust its workload by executing it one unit at a time: The action predicates enclosed between curly braces enable to capture the information present on Lts actions, i.e., the channel names (interpreted as character strings) and the data values transmitted. The box modality matches all sequences that end up, after zero or more steps, with an execute action carrying a channel name ld and a workload load. These values are captured and used later in the parameterized minimal fixed point operator µX, which expresses the inevitable reachability of consecutive execute actions that carry decreasing workloads.

The Pic2Lnt 2.0 translator is currently provided as a Cadp plug-in. The manual page and the executable files for several architectures (Mac computers, Pcs running Linux or Windows, Solaris workstations) are available on-line [START_REF] Mateescu | Pic2Lnt: A translator from the applied pi-calculus to lotos nt[END_REF].

Experimental Evaluation

We applied Pic2Lnt 2.0 on a benchmark of π-calculus specifications, which includes most of the examples provided with Mwb, as well as applied π-calculus examples that we specified ourselves. Our benchmark currently contains 284 files, totalizing about 5, 200 lines of π-calculus, which were translated in about 50, 000 lines of Lnt. This expansion in size, which is negligible given the speed of the Lnt compiler, is caused partly by the complexity of the translation (one new Lnt channel per parallel composition operator) and partly by the verbosity of Lnt w.r.t. the compact algebraic notation of the π-calculus.

The table below shows a few examples from the Mwb distribution. For each example, the table gives the number of agents, the size of the specification before and after translation, and some quantitative information (size, time) about the Lts generated using Pic2Lnt 2.0 and the Cadp exploration tools. Our tool support for the applied π-calculus is already used for teaching purposes at the University of Saarbrücken (Germany). It is also currently used for specifying and verifying self-deployment and other self-management protocols designed in the context of the OpenCloudware4 project, which aims at building an open software engineering platform for the collaborative development of distributed applications to be deployed on multiple cloud infrastructures. Since the applied π-calculus is convenient for specifying many kinds of mobile systems (e.g., Web services, autonomic applications, cloud computing protocols, software architectures, biological systems, cryptographic protocols, etc.), we believe that our tool support can provide a useful service in a wide range of application areas.

Concluding Remarks

We introduced in this paper an applied π-calculus equipped with data-handling features, and proposed a translation of this language into the Lnt value-passing process algebra. This translation, automated by the Pic2Lnt 2.0 tool, enables the analysis of applied π-calculus specifications using all verification tools of Cadp. As far as we are aware, this results in one of the few operational frameworks for verifying an applied π-calculus. ProVerif [START_REF] Blanchet | Proverif: Automatic cryptographic protocol verifier, user manual and tutorial[END_REF] is an alternative approach focused on the verification of cryptographic protocols and security properties (secrecy, authentication, etc.). In contrast, our solution is independent of any application domain and provides a larger panel of verification techniques.

  Main = (ν task , al, ar, bl, br, cl, cr) ( Client(task ) | LoadBalancer (task , al, ar, bl, br, cl, cr) | Server (al, ar, 0 of Nat) | Server (bl, br, 0 of Nat) | Server (cl, cr, 0 of Nat) )

Fig. 1 .

 1 Fig. 1. Overview of the tool chain

[

  true * .{execute ?ld :String ?load :Nat ... where load > 1}] µX(crt load :Nat := load -1).( true true ∧ [¬{execute !ld !crt load ...}] X(crt load ) ∧ [{execute !ld !crt load ... where crt load > 1}] X(crt load -1) )

  2 of Nat .0 LoadBalancer (task , al, ar, bl, br, cl, cr) = task (com). al.al (v1 : Nat ). bl .bl (v2 : Nat). cl .cl (v3 : Nat ). ( [(v1 ≤ v2) and (v1 ≤ v3)] ar com .LoadBalancer (task, al, ar, bl, br, cl, cr) + [(v2 ≤ v1) and (v2 ≤ v3)] br com .LoadBalancer (task, al, ar, bl, br, cl, cr) + [(v3 ≤ v1) and (v3 ≤ v2)] cr com .LoadBalancer (task, al, ar, bl, br, cl, cr) ) Server (ld, rq , totalload : Nat) = ld.ld totalload .Server (ld, rq , totalload ) + rq(req ).req (newload : Nat ).Server (ld, rq, totalload + newload ) + [totalload > 0] execute ld, totalload .Server (ld , rq, totalload -1)

The version 1.0 of Pic2Lnt, which handled the original π-calculus (without data manipulation), consisted of about 3, 700 lines of code.

http://opencloudware.org

Acknowledgments. We are grateful to Hubert Garavel for his valuable feedback about the applied π-calculus and the connection of Pic2Lnt 2.0 to Cadp.