Aude Bolopion
email: aude.bolopion@isir.upmc.fr

Barthélemy Cagneau

Stephane Redon

Stéphane Régnier

Comparing position and force control for interactive molecular simulators with haptic feedback

Keywords: molecular simulation, haptic coupling, force feedback, transparency, stability

This paper presents a novel tool for the analysis of new molecular structures which enables a wide variety of manipulations. It is composed of a molecular simulator and a haptic device. The simulation software deals with systems of hundreds or thousands of degrees of freedom and computes the reconfiguration of the molecules in a few tenths of a second. For the ease of manipulation and to help the operator understand nanoscale phenomena, a haptic device is connected to the simulator. To handle a wide variety of applications, both position and force control are implemented. To our knowledge, this is the first time the applications of force control are detailed for molecular simulation. These two control modes are compared in terms of adequacy with molecular dynamics, transparency and stability sensitivity with respect to environmental conditions. Based on their specificity the operations they can realize are detailed. Experiments highlight the usability of our tool for the different steps of the analysis of molecular structures. It includes the global reconfiguration of a molecular system, the measurement of molecular properties and the comprehension of nanoscale interactions. Compared to most existing systems, the one developed in this paper offers a wide range of possible experiments. The detailed analysis of the properties of the control modes can be easily used to implement haptic feeback on other molecular simulators.

Introduction

The conception of new medicines, the synthesis of alloys, the prototyping of bio-nanorobots [START_REF] Hamdi | Prototyping bionanorobots using molecular dynamics simulation and virtual reality[END_REF] as well as the transmission of information between two nanomachines [START_REF] Moore | Interfacing with nanomachines through molecular communication[END_REF], [START_REF] Atakan | An information theoretical approach for molecular communication, Bio-Inspired Models of Network[END_REF], and the development of bio-microelectromechanical system devices (bio-MEMS) [START_REF] Rossi | Molecular and supramolecular systems for sensing and actuation[END_REF], [START_REF] Grove | Packaging actomyosin-based biomolecular motor-driven devices for nanoactuator applications[END_REF] are made possible thanks to molecular simulators. Since ease of manipulation is one of the key points of such tools, haptic devices are widely used. They provide a natural manipulation mode of the molecules and force information is easily interpreted by operators. It has been demonstrated that it helps operators to understand nanoscale phenomena. Its use is recommended for educational purposes [START_REF] Persson | Designing and evaluating a haptic system for biomolecular education[END_REF], [START_REF] Murayama | An educational environment for chemical contents with haptic interaction[END_REF], [START_REF] Sankaranarayanan | Role of haptics in teaching structural molecular biology[END_REF], but it is not limited to academic courses. It can also help chemists and biologists to find specific locations in complex molecular systems (for example docking sites) [START_REF] Brooks | Project GROPE -Haptic displays for scientific visualization[END_REF], or to characterize the molecular properties [START_REF] Hamdi | Molecular mechanics simulation of bio-nanorobotic components using force feedback[END_REF].

The benefits of haptic feedback depends on the coupling used. In particular, how to reach the entire virtual environment using a haptic device with a limited workspace is a key issue to get an interactive system. Several techniques are proposed, from the concept of clutching (freezing the displacement of the virtual object while enabling the user to modify the position of the haptic handle), to the Bubble technique (position control is used for fine positioning, and rate control is used for large displacements) [START_REF] Dominjon | The "Bubble" technique: interacting with large virtual environments using haptic devices with limited workspace[END_REF]. As an alternative, the concept of Active Haptic Workspace is considered [START_REF] Subasi | A new haptic interaction and visualization approach for rigid molecular docking in virtual environments[END_REF]. As in the Bubble technique position control is used inside the AHW, and objects are rotated and translated when the cursor is outside. In all cases, when the user is close to the point of interest, a common position coupling scheme is used to control the object and send the interaction forces to the user. Thus, even if complex techniques are used, the choice of the haptic coupling used remains unchanged when dealing with a precise area of the workspace. The haptic coupling structure depends on the available inputs and outputs. In almost all the works dealing with molecular systems, position control is used (the user sets the position of the molecule, and feels the interaction forces through the haptic interface). However, stability is difficult to ensure due to long computation times, scaling factors used to link the macro and nano worlds, and the high variation of the forces. In most of the works, the system stability is guaranteed at the expense of the fidelity of the force feedback. Either the accuracy of the molecular interactions computed is decreased by using simpler models [START_REF] Lee | Smoothing haptic interaction using molecular force calculations[END_REF], or the damping added to the coupling deteriorates the transparency [START_REF] Daunay | Stable six degrees of freedom haptic feedback for flexible ligand-protein docking[END_REF]. Scaling factors are also a major issue for stability [START_REF] Wollacott | Haptic applications for molecular structure manipulation[END_REF]. Position control may also lead to instabilities as the position of the molecule is first set according to the user input, before being corrected during the update step of the simulation using the equation of motion to take into account physical constraints, as it is explained in this paper. In addition to the haptic coupling specificity described above, the system must fulfill several requirements to be usable for the analysis of molecular structures. It should enable the manipulation of the molecules, the measurement of their inherent properties and should provide a deep understanding of the nanoscale interactions. According to potential users (including bio-physicist researchers and researchers working in pharmacological firms), the system should present:

• R 1 : a fast and accurate manipulation and reconfiguration computation,

• R 2 : different manipulation choices (both the whole molecule and a single atom), different force feedback (both the internal and external interaction forces),

• R 3 : transparent and stable haptic coupling schemes despite the computation time.

The system we propose is composed of a molecular simulator, SAMSON, and two haptic coupling schemes. SAMSON enables a fast computation of the interaction forces and of the reconfiguration of the system thanks to a tree representation. The molecules can be manipulated by setting both a position and a force thanks to the methodology used to compute the solution of the equation of motion. Regarding the haptic coupling, a force control scheme is used in addition to position control to compensate for the lack of stability of the position mode while dealing with large displacements. It will also be seen that this control mode is more suited to molecular simulators. Moreover, new applications are possible thanks to this coupling. The control schemes presented are widely used for macroscale teleoperations. Even if force control has already been mentioned for molecular simulation in particular in [START_REF] Stone | A system for interactive molecular dynamics simulation[END_REF] and [START_REF] Delalande | Complex molecular assemblies at hand via interactive simulations[END_REF], to our knowledge, this is the first time it is detailed and analyzed for molecular simulation. It can be implemented in our system thanks to the molecular dynamics equation solved.

We did a preliminary work where some of the potential applications of our system (conception of new molecules, analysis of molecular properties, ...) are described [START_REF] Bolopion | Haptic feedback for molecular simulation[END_REF]. We now present a comparison of position and force control. A detailed description of the molecular simulator and the control schemes, as well as how they are connected, is made to compare the adequacy of the two control modes to molecular dynamics simulators. The stability sensitivity of the coupling schemes is considered for the specific application of molecular interactions. The sensitivity is defined as the variation of the magnitude of the roots of the control schemes' characteristic equation with respect to given parameters. Relevant experiments of all the stages of the analysis of a molecule are performed and highlight the specificity of each control mode. Based on the experimental results, it seems that the system we propose could be used for the analysis of new molecules, thanks to the two complementary control modes, and to its adequacy with the requirements (R 1 -R 3 ).

The rest of the paper is organized as follows. The molecular simulator is presented in Section 2. Based on the available inputs and outputs, the two control modes (position and force control) are detailed and compared in Section 3 in terms of adequacy with the simulator. The transparency and the sensitivity of the coupling schemes' stability is analyzed in Section 4. In Section 5, examples of manipulations that can be handled by our tool are presented. The advantages and drawbacks of each control scheme discussed in previous sections are compared with these experimental data. Section 6 concludes the paper. A table summarizing the notations is provided in Section 7.

Molecular simulator

As mentioned in Section 1, the molecular simulator must deal with reconfigurable objects of hundreds or thousands of degrees of freedom. Several approaches have been used, such as simplifying the force fields, or considering only rigid molecules. The resulting loss of accuracy limits the possible applications. The simulator we detail hereafter uses a special representation of the molecules, and a specific algorithm based on a quasistatics method to simulate the motion of the molecular system. This leads to a fast computation of both the internal and external forces and a fast reconfiguration of the molecule [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF], [START_REF] Grudinin | Practical modeling of molecular systems with symmetries[END_REF]. Since the simulator has an open architecture, and most functions are implemented as plug-ins, the haptic coupling can be directly integrated into the equations solved by the simulator. This increases the variety of control modes that can be proposed compared to most of the commercially available software.

Representation of a molecule

The method relies on a recursive, divide-and-conquer representation of the molecular system, where any non-rigid group of atoms is considered to be the union of two groups of atoms. Thus, any molecular system is associated to a binary assembly tree, in which the leaf nodes represent the user-defined rigid bodies (atoms, or group of atoms). The internal nodes represent the sub-assemblies, and the root node represents the complete molecular system (Figure 1). Using the spatial notation, the equation of motion of each rigid body is [START_REF] Featherstone | A divide-and-conquer articulated-body algorithm for parallel O(log(n)) calculation of rigid-body dynamics. Part 1: basic algorithm[END_REF]:

I • a = F ext + F int -v × I • v (1) 2
where a (resp. v) is the spatial acceleration (resp. velocity) of the molecule. I is its inertia. F ext is the sum of the forces applied by other molecules and F int is the sum of the internal forces applied to the considered rigid body. These forces depend on the force field being used in the simulation to model interatomic interactions. In our current implementation, the force field is derived from a well known molecular mechanics force field, CHARMM [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF], which models interactions through van der Waals, electrostatic and dihedral contributions.

Computation of the solution of the equation of motion

To solve the equation of motion for each rigid body, the simulator follows traditional approaches: at each time step, the interatomic forces are first computed, then the acceleration of the system is determined and used to update its state and the internal configuration of the molecule. This operation is made easier thanks to the tree representation, which enables to recursively compute the forces and accelerations from the root node to the leaves after having recursively computed the coefficients of the equation in a bottom-up pass. Thanks to the chosen representation of the molecule, and the way of solving the equation of motion, systems of hundreds of degrees of freedom can be simulated in tenths of a second.

Manipulation modes

To enable the user to manipulate the molecule (or one of the rigid bodies), its action must be taken into account while solving the equation of motion. Thanks to the open architecture of the software, two solutions are possible. They are detailed hereafter:

Position mode. One way to control the simulation through the haptic device is to control the position of the manipulated rigid body or molecule. In the algorithm, the instantaneous position is set by the user, instead of being determined by the computed acceleration. The internal degrees of freedom are simulated using Equation (1). This manipulation mode is intuitive for the user, since we are used to set the position of macroscopic objects. However, between two simulation steps, the position set by the user may be physically unacceptable (e.g. atom clashes may happen).

Force mode. Another way to control the simulation is to apply a force on the manipulated objects. This force F i is added to the interatomic forces that are applied to the object, and the total force is used to determine the accelerations in the molecular system and update the positions of the rigid bodies. The equation of motion is then:

I • a = F i + F ext + F int -v × I • v (2)
Based on these manipulation modes, two corresponding control schemes may be used for position and force control.

Force sent to the user

The simulator determines the force F m applied to the manipulated molecule. It is composed of van der Waals and electrostatic contributions from other molecules F ext , as well as those resulting from other rigid bodies in the selected molecule F int :

F m = F ext + F int (3)
The user will thus feel both internal and external forces. Note that, in this sum, the contributions from rigid bodies within the selected molecule cancel out if the user manipulates a whole molecule (F int = F j int = 0). When only one rigid body from a molecule is selected, the molecular force F j m applied to the j th rigid body by the environment is computed as before, but the kinematic constraint forces, i.e. the forces which maintain the kinematic constraints between the selected rigid body and its topological neighbors [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF] do not cancel out. In particular, if the rigid body is part of a molecule on which no external force is applied, the user will only feel the kinematic constraints. Since no modification on our system are necessary the user can easily switch from one mode to the other.

Note that, although we have used SAMSON simulator in our implementation, the control schemes we will present may be connected with any simulator which provides access to the required variables (positions, forces, etc.).

Coupling between the simulation software and the haptic interface

Before considering the two control schemes some notations must be introduced for both the simulator and the haptic coupling. These notations are summarized in Section 7.

Notations 3.1.1. Simulator

The n rigid bodies (atoms or groups of atoms) which are inside the manipulated molecule will be noted R j m , j = 1, ..., n and the q ones which are part of the environment will be noted R l e , l = 1, ..., q (see Figure 2). The force applied to each rigid body R j m by all the other rigid bodies of the simulation (including internal forces F j int and external forces F j ext ) is:

F j m = F j int + F j ext . The total force applied to the manipulated molecule F m is: F m = j F j m .
The position of R j m after the reconfiguration step is x j m , and x m is the position of the center of mass of the manipulated molecule. Similarly, F l e = F l int + F l ext is the force applied to each rigid body R l e by all the other rigid bodies of the simulation, x l e is its position after the reconfiguration step. m p is the mass of the p th rigid body in the simulation. m is the total mass of the manipulated molecule. The molecular dynamics do not allow to use the same time basis as the one we use [START_REF] Stone | A system for interactive molecular dynamics simulation[END_REF]. A single simulation step takes around 100 milliseconds (moving the molecule, reconfiguring the system, sending the interaction forces to the user), and this corresponds to a period of around T m = 2 femtoseconds for the evolution of the molecular system. The molecular integration time T m is thus introduced.

Haptic coupling

When connecting a haptic (or master) device to a molecular simulator, the delays due to long computations of complex interactions must be taken into account while analyzing the stability of the control scheme. The transfer functions are thus represented using the discrete time variable z. Since the variables coming from the simulator are delayed by one sampling period T s (defined as the largest simulation period for a given experiment), an explicit one step delay is introduced between the simulation and the coupling. Note that we consider here the worst case by taking the largest delay. This is modeled by the delay operator z -1 (see Figure 3 and Figure 4). The discrete time transfer function H(z) of the haptic interface (which input, resp. output, is the force applied to the haptic handle, resp. its position) is computed from the continuous time domain using the Z-transform function [START_REF] Ogata | Discrete-time control systems[END_REF], and taking into account the effects of the sampler and zero order hold:

H(z) = 1 B h T s (z -δ) -(1 -δ)(z -1) M h B h (z -1)(z -δ) (4)
where δ = e -B h Ts M h . M h and B h are respectively the master's inertia and viscosity, as given in [START_REF] Diolaiti | Stability of haptic rendering: Discretization, quantization, time delay, and coulomb effects[END_REF]. Due to the scale difference between the macro and the nanoworlds, two scaling factors must be introduced to connect the haptic interface and the molecule. The force (resp. position) scaling factor is denoted by A f (resp. A d ). Depending on the manipulation mode, the user may either set the desired instantaneous position of the molecule x i , or apply an additional force F i on it.

Based on the two manipulation modes described in Section 2 and the notations given above, position and force control are compared.

Position control

The first control scheme, the classical impedance display or Direct Force Feedback (DFF) coupling, is based on position control.

The Direct Force Feedback control scheme

The DFF coupling, represented in Figure 3, is the simplest structure to connect the haptic device to the simulation using the scaled position of the haptic handle x i as the input of the molecular dynamics simulator. It is also the most intuitive, since the user directly controls the position of the molecule and feels the interaction forces F h between the molecule and its environment. Only two scaling factors are necessary to control the molecule.

F m F h x i F op H(z) A -1 d A f z -1 + -
x h

Details of the position control algorithm

As mentioned in Section 2, the calculation of positions and forces performed by the simulation software are made on each rigid body. However, for ease of manipulation, the user manipulates the molecule as a whole, controlling the position of the center of mass x i . As explained in Section 2.2.1, the instantaneous position x i of the rigid body is set by the user. Internal and external forces are computed and used to update the positions of the molecules in the environment and their internal degrees of freedom (to simulate their flexibility). This reconfiguration step involves a modification of the rigid bodies positions, and thus the actual position of the center of mass of the manipulated molecule x m might be different from the one set by the user using the haptic interface x i . The user feels the force applied on the manipulated molecule F m composed of internal and external forces (Equation ( 3)). This methodology is summarized in Algorithm 1 1 .

The DFF coupling has one main drawback, due to the manipulation mode used. The user sets an arbitrary position using the haptic interface. The molecular system first takes the corresponding configuration before being corrected during the update step by the simulator. This can lead to instabilities, and the forces sent to the user may misrepresent the correct molecular interactions.

Force control

As an alternative to position control, we can consider a second manipulation mode, force control.

Algorithm 1 Position control algorithm (see Section 7 for the notations)

1. Compute the desired position of the center of mass x i (k+1) set by the operator using the control scheme (position of the haptic handle scaled by A d ). 2. Set the position of the rigid bodies to x j i (k + 1), computed based on the relative positions with respect to x i . 3. Update all interatomic forces F 1 m (k +1), ..., F n m (k +1) (resp. F 1 e (k + 1), ..., F q e (k + 1)) applied to the rigid bodies which are inside (resp. outside) the manipulated molecule using the molecular simulator.

Compute the sum

F m (k + 1) = n j=1 F j m (k + 1)
of the forces applied to the molecule, and send it back to the user. 5. Update the positions x j m (resp. x l e ) of the rigid bodies which are inside (resp. outside) the manipulated molecule using the quasi-statics simulator. Note that in the following equations x j i (k + 1) (resp. F j m (k + 1) and F l e (k + 1)) are known since they have previously been computed Step 2 (resp. 3) of the algorithm:

x j m (k + 1) = x j i (k + 1) + T 2 m m j F j m (k + 1)
x l e (k

+ 1) = x l e (k) + T 2 m m l F l e (k + 1)

The Force-Force control scheme

The Force-Force (FF) control scheme is depicted in Figure 4. The input of the simulation is the force applied by the user on the molecule. In addition to the scaling factors, a proportional gain k i is introduced. This gain is used to adjust the magnitude of the force applied to the molecule, which is based on the difference between the molecule and the haptic handle positions.

F m x m F h x h F i F op H(z) z -1 z -1 A -1 d k i A f + + - -

Details of the force control algorithm

Contrary to the DFF algorithm, the simulator's input is the force F i that the user wants to apply to the entire molecule. This force is set by controlling the distance between the haptic handle and the molecule positions. The corresponding forces applied on each rigid body F j i of the manipulated molecule as well as the environment forces are used to update the position of all the rigid bodies (inside and outside the manipulated molecule) using Equation [START_REF] Moore | Interfacing with nanomachines through molecular communication[END_REF]. As for position control, the user feels the force applied on the manipulated molecule F m composed of internal and external forces (Equation ( 3)). This is summarized in Algorithm 2.

Algorithm 2 Force control algorithm (see Section 7 for the notations)

1. Compute the force applied by the operator using the control scheme:

F i (k + 1) = k i x h (k+1) A d -x m (k)
2. Compute the force F j i (k + 1) = F i (k + 1)/n, applied to all the j th rigid bodies of the manipulated molecule. 3. Update the positions x j m (resp. x l e ) of the rigid bodies which are inside (resp. outside) the manipulated molecule using the quasi-statics simulator. Note that in the following equations, F j i (k + 1) is known since it has previously been computed (Step 2 of this algorithm):

x j m (k + 1) = x j m (k) + T 2 m m j (F j m (k) + F j i (k + 1)) x l e (k + 1) = x l e (k) + T 2 m m l F l e (k) 4. Update all interatomic forces F 1 m (k +1), ..., F n m (k +1) (resp. F 1 e (k + 1), ..., F q e (k + 1
)) applied to the rigid bodies which are inside (resp. outside) the manipulated molecule using the molecular simulator. The user input might be unrealistic in both position and force control modes. However, given the user input, its integration into the equation of motion differs. For position mode, the system configuration is first set according to the user input, and in a second step, it is corrected based on the equation of motion. On the contrary, using force control enables to integrate the user's input while solving the equation of motion (Equation ( 2)). The system behaves as if these forces were actually applied. Thus, it is more adapted to molecular dynamics than position control. As force mode presents a better adequacy with the solved equations, instabilities due to molecule manipulations are reduced.

Comparison of transparency and stability properties

In this section, the transparency and the sensitivity of the control schemes stability are compared. The analysis is made on one axis: scalar variables are used. The action of the operator is not taken into account, and a simplified model of the environment is considered (an equivalent spring constant k e ). This is obviously not enough to model the complex interactions between molecules, but comparisons of the control schemes can be made. Transparency analysis will be carried out in continuous time, without considering the delays. These assumptions are made to derive comprehensive relationships to compare the two control modes [START_REF] Bolopion | Analysis of stability and transparency for nanoscale force feedback in bilateral coupling[END_REF]. On the contrary, stability will be studied in discrete time as we will be interested on the sensitivity of the control schemes with respect to the sampling period.

Transparency analysis

A haptic coupling is transparent if it perfectly renders molecular forces and molecular dynamics. It is measured by comparing the impedance of the environment Z e = F m /V m and the one felt by the operator Z op = F op /V h [START_REF] Lawrence | Stability and transparency in bilateral teleoperation[END_REF]. Ideal transparency is achieved when:

Z op (s) = Z e (s) (5)
where s is the parameter of the Laplace transform. For submicron scales, the scaling factors A f and A d must be considered so that the impedances can be compared. In our context, perfect transparency will be achieved if:

F op (s) V h (s) = A f A d F m (s) V m (s) ⇔ Z op (s) = A f A d Z e (s) (6)
Two distinct modes will be considered: 1. non-contact: there is no interaction between the manipulated molecule and its environment (Z e = 0) 2. contact: there are repulsive forces, modeled as a first approximation as a spring (Z e (s) = k e s ) Note that, in theory, non-contact configurations do not exist in molecular systems since interactions can be felt at infinite distance (e.g. electrostatic interactions). However, molecular dynamics simulators frequently use cutoffs, i.e. distances over which interactions are considered to be zero. This makes relevant the analysis on non-contact configurations. An approximated impedance will be derived in the frequency domain (s = jω) considering low (ω << 1) or high (ω >> 1) frequencies.

The Direct Force Feedback control scheme

Using the control scheme depicted in Figure 3, the impedance felt by the operator is computed:

Z DFF op (s) = F op (s) V h (s) = A f A d Z e (s) + 1 H(s) (7)
Non-contact. When no force is applied to the molecule, Equation ( 7) is simplified to:

Z DFF op ( jω) = M h jω + B h (8)
At low frequencies, the operator will mainly feel the viscosity of the haptic device. At high frequencies, the inertia of the master arm will be predominant.

Contact. When repulsive forces are applied to the manipulated molecule, the impedance felt by the operator is given by Equation [START_REF] Murayama | An educational environment for chemical contents with haptic interaction[END_REF]. This corresponds to the impedance he or she should ideally feel (

A f A d Z e
) modulated by the haptic device characteristics. For low frequencies, the impedance Z DFF op,LF is approximated by:

Z DFF op,LF ( jω) ≈ ω<<1 A f A d Z e ( jω) + B h (9)
The viscosity of the haptic interface deteriorates the transparency. For high frequencies, the impedance Z DFF op,HF is:

Z DFF op,HF ( jω) ≈ ω>>1 M h jω ( 10
)
The transparency of the coupling is only affected by the inertia of the haptic device.

The Force-Force control scheme

For the FF coupling, the impedance felt by the operator is:

Z FF op (s) = 1 H(s) + A f k i A d (ms 2 + k i + k e ) Z e (s) (11)
Non-contact. As for the DFF coupling, the operator's impedance is given by Equation ( 8), and he or she feels the viscosity of the haptic device for low frequencies, and its inertia for high ones.

Contact. When the contact is established, Z FF op,LF can be computed for low frequencies:

Z FF op,LF ( jω) ≈ ω<<1 k i k i + k e A f A d Z e ( jω) + B h ( 12
)
Equation ( 12) differs from the ideal impedance

A f
A d Z e ( jω) due to the static gain k i k i +k e , and the viscosity of the haptic interface. For high frequencies, the impedance felt by the user is the inertia of the haptic device.

Summary

Table 1 summarizes the approximated impedances for the two control schemes. For the DFF control scheme, the quality of the force rendering only depends on the haptic device characteristics. The transparency of this coupling is close to the ideal transparency. The FF transparency differs from the DFF one for the contact mode, for low frequencies. The ideal impedance is scaled by k i k i +k e for the FF coupling. To tend to perfect transparency, the value of k i must be greater than k e . This means that the force applied by the user is greater than the ones coming from the molecular system. For an infinite value of k i the molecule position will only depend on the force applied by the user, and the manipulation mode will tend to the behavior of the position manipulation mode. On the contrary, if k i is much lower than k e , the perfect impedance will be scaled by k i k e . We will see in Section 5.1 that it corresponds to divide the perfect impedance by a factor less than two as long as k e ≤ 100N.m -1 . The transparency of this coupling is thus only slightly deteriorated.

Stability analysis

Since passivity is conservative (it ensures stability for any operator and any environment as long as they are themselves passive) [START_REF] Khalil | Nonlinear systems[END_REF], it is not used in this work. We chose instead to compare the sensitivity of the control schemes' stability with respect to the sampling period. Thus discrete time must be used. The computational time is taken into account by introducing a one period delay on the variables computed by the simulator (see Figure 3 and Figure 4). The sensitivity analyzed is the variation of the magnitude of the roots of the characteristic equations with respect to given parameters, which is here the sampling period.

Characteristic equations

The discrete time transfer function of the DFF control scheme is derived from the coupling presented in Figure 3:

Z x h (s) F op (s) = X h (z) F op (z) (13) = T s z-δ z-1 -(1 -δ) M h B h B h (z -δ) + T s z-δ z-1 -(1 -δ) M h B h A f A d k e z -1 ( 14
)
where Z is the Z-transform function, and δ is the parameter defined in Equation ( 4). The corresponding characteristic equation is:

a 3 z 3 + a 2 z 2 + a 1 z + a 0 = 0 ( 15
)
where:

a 3 = B h , a 2 = -B h (1 + δ) a 1 = δB h + T s -(1 -δ) M h B h A f A d k e a 0 = -δT s + (1 -δ) M h B h A f A d k e (16)
The transfer function of the FF coupling is computed according to Figure 4:

X h (z) F op (z) = H(z)

+ H(z)G(z)

A f A d [START_REF] Delalande | Complex molecular assemblies at hand via interactive simulations[END_REF] where:

G(z) = F m (z) X h (z) = k i k e T 2 m mz(z-1)+k i T 2 m +k e T 2
and m is the mass of the whole molecule. The characteristic equation is:

b 4 z 4 + b 3 z 3 + b 2 z 2 + b 1 z + b 0 = 0 ( 18
)
where:

b 4 = mB h , b 3 = -(2 + δ)mB h b 2 = (k i + k e )T 2 m B h + mB h (1 + 2δ) b 1 = -(1 + δ)(k i + k e )T 2 m B h -δmB h + T s -(1 -δ) M h B h k i k e T 2 m A f A d b 0 = B h δ(k i + k e )T 2 m + -T s δ + (1 -δ) M h B h k i k e T 2 m A f A d ( 19
)
Analytical stability criteria can be obtained for both of these control schemes using the Jury criterion [START_REF] Ogata | Discrete-time control systems[END_REF]. However, the relations derived are complex and the influence of the control scheme parameters on stability is not highlighted. In the next section, a numerical analysis is performed.

Numerical comparison of the control schemes' stability

A numerical analysis of the influence of the environment parameters k e and T s on each control scheme is performed. The variations of the highest magnitude of the roots of the characteristic equations ( 14) and ( 17) ( r max ) are considered in Figure 5. The numerical values correspond to the experiments carried out in Section 5.1 (see Table 3). The plots of Figure 5 represent the value of r max for both the control schemes, depending on k e and T s . The variations of the magnitude of the roots are of particular interest. For the two control schemes, long computations (large values of T s ) and/or stiff environments may lead to an unstable system. Moreover the magnitude of the roots of the DFF control scheme quickly increases when T s or k e are increased. On the contrary, the roots obtained with the FF control scheme are almost constant. Thus, the DFF control scheme stability is more sensitive to the variations of T s than the FF one. These results will be experimentally confirmed in the next section.

Comparison of the control schemes' specificity

As seen from Section 2 to Section 4, the position and force control modes present different characteristics in terms of adequacy with the algorithm, transparency and stability, which are summarized in Table 2. Position control is more transparent. However, force control is better suited for the integration into the simulator since the user's input is directly taken into account while solving the equation of motion, and its stability is less sensitive to environmental conditions.

Next sections detail the possible applications of our system, and highlight the advantages and drawbacks of both position and control mode depending on the given task.

Molecular manipulations and haptic feedback of interaction phenomena

As seen in Section 1, a wide variety of applications must be considered. To illustrate the possibilities of our system, the following experiments are performed in the next sections: • manipulation of whole molecules to modify the global structure of the molecular system,

• deformation of molecular structures to measure their inherent properties,

• determination of equilibrium positions which is a common example of the specific phenomena of nanoscales that users would like to understand.

Experiments with varying molecules and parameters have been carried out. For the sake of clarity, only some of them which illustrate the three types of operation listed above are presented. Similar conclusions were derived for the other cases we considered.

The experiments are performed on a 3GHz Intel PC processor, using a Phantom Omni haptic interface from Sensable2 . This master arm is a 6 degrees of freedom (dof ) movement and a 3 dof force-feedback device, with an inertia M h = 0.072kg and a viscosity B h = 0.005N.s.m -1 [START_REF] Diolaiti | Stability of haptic rendering: Discretization, quantization, time delay, and coulomb effects[END_REF]. The maximum force that it can provide is 3.3N. In the following experiments, the maximum force sent to the haptic interface is limited to 3N. The OpenHaptics development toolkit provided by Sensable is used to communicate with the haptic interface. The setup is presented in Figure 6.

Manipulation of whole molecules, the HIV example

The haptic device must enable the user to manipulate the whole molecule into the environment, so that it would be equivalent for him or her to manipulate their center of mass in the case of rigid molecules. This mode has been chosen since it enables an untrained user to precisely and intuitively move the molecules. Two solutions are possible. The first one is to set the position of its center of mass, and the second one is to apply a force on each of the rigid bodies. To compare these two manipulation modes and the associated control schemes, the HIV protease and its inhibitors (Protein Data Bank code: 2AZ8) is considered. This system is composed of 848 degrees of freedom. The parameters used for the simulator and the haptic coupling schemes are summarized in Table 3. These values are the ones used for the numerical comparison of stability in Section The experiment is repeated twice for both the DFF and FF coupling schemes, with varying number of activated dof (172 for the first experiment, 848 for the second one), and thus varying simulation periods (T s = 92 and 117 ms). This is possible thanks to the adaptive simulation algorithms used by our simulator which enable to decrease the number of active degrees of freedom with a minimal loss of accuracy [START_REF] Rossi | Adaptive torsion-angle quasi-statics: a general simulation method with applications to protein structure analysis and design[END_REF]. In order to compare the results, the movement of the user has been limited to one axis (with 3D force feedback). The forces returned to the operator on the y axis of the haptic interface with respect to the position of the molecule's center of mass are depicted in Figure 7 and8.

Three phases can clearly be distinguished. The first one is the extraction of one of the inhibitors, when it is situated at less than 16 Å of its initial position. The second phase consists in manipulating the inhibitor away from the rest of the structure, when the inhibitor is further than 16 Å. No force is applied on it. The reinsertion into the system occurs when the position is again less than 16 Å. The forces felt using the DFF and FF control schemes are sim- ilar. For the two coupling schemes, the magnitude of the haptic forces are in the order of the newton so that the interaction forces can clearly be felt by the user. The differences on the force profile are due to the flexibility of the molecules. Depending on the manipulation conditions (velocity of the molecule, force exerted, etc.) the deformations are not the same, leading to slightly different configurations and interaction forces. As predicted in Section 4.2, the DFF coupling is very sensitive to the environment characteristics. Only a few oscillations can be seen in Figure 7 for very stiff interactions. When the sampling period increases, oscillations become important and more frequent. On the opposite, the FF coupling remains stable even for stiff environments, and with long simulation periods. Adding damping to the DFF coupling would improves the stability, but transparency would be deteriorated. Since molecular interactions involve complex variation of the forces and efforts of small magnitude (in particular for attractive interactions), transparency is a major concern to provide to the user a high fidelity force feedback.

Deformation of molecular structures

The measurement of the molecule inherent properties (stiffness, inertia, ...) is of primary importance [START_REF] Hamdi | Prototyping bionanorobots using molecular dynamics simulation and virtual reality[END_REF]. To deform such structures, it is necessary to apply an external force on the molecule. This is done using force control.

The experiment presented in this paragraph consists in manipulating a single rigid body of a flexible molecule using the FF coupling, and feeling the internal forces and the inertia. The same methodology as the one described in Algorithm 2 is used, except that the whole molecule is replaced by a single rigid body R j m and the kinematic constraint forces are added as mentioned in Section 2. Since no other molecules are present, the external force applied to the rigid body are null (F ext = 0), and the operator will only feel the internal forces F m = F int . Using the haptic device the operator applies a force on one of the rigid bodies of the HIV inhibitor to unfold it. The other extremity is fixed. When the molecule is stretched, the user moves back the haptic device to the initial position. Doing so, the molecule is refolded. Figure 9 represents the norm of the forces felt by the user against the norm of the position of the center of mass. The forces increase since the molecule is stretched (part A). When the molecule is deployed, the user applies an additional force and the internal forces increase (part B). After that, the user moves back the rigid body and he or she releases the constraint applied on it. The forces quickly decrease because of the reduction of internal forces (part C). Finally, the molecule is moved toward its initial position. The interaction forces due to the other rigid bodies are felt (part D). This shape is consistent with the one described in [START_REF] Hamdi | Prototyping bionanorobots using molecular dynamics simulation and virtual reality[END_REF]. Our system is able to render small variations of the forces to the user (in the range of the nanonewton). The FF coupling is well suited to deform molecules by applying an external force, and can be used to measure their inherent properties.

Determination of equilibrium positions

Using force control, the user applies an effort on the manipulated molecule, and the equilibrium state of the system takes into account this external force. However, to better understand molecular interactions, it is necessary to know the equilibrium state around a given position without considering the action of the user. Position control with force feedback enables the user to get such a knowledge.

An experiment involving two water molecules is performed using the DFF coupling. Due to van der Waals forces, there is an equilibrium position when the distance of the two water molecules is around 2.6Å. Attractive and repulsive forces tend to keep the molecules in this position. In the following experiment, one of the molecules is fixed. The user manipulates the second one around the equilibrium position, while feeling the interaction forces which tend to keep the system into the equilibrium state. Experimental results are given in Figure 10. Position and force are plotted with respect to the time to evaluate the time taken by the system to reach the equilibrium position. The user can feel both attractive and repulsive forces and can clearly determine the equilibrium position since it corresponds to a null force feedback. This helps him or her understand one of the most specific molecular interactions, the long range attractive and repulsive phenomena.

Conclusion

In this work, we have presented a new tool for the analysis of molecular structures. The simulator computes the reconfiguration of molecular systems composed of hundreds or thousands of degrees of freedom in a few tenths of a second. Thanks to its opened architecture, the haptic coupling can directly be integrated into the equation of motion, leading to two possible control modes. The associated DFF and FF haptic coupling schemes are complementary since high transparency is achieved for position control. Force control presents a better adequacy to molecular dynamics and its stability is not sensitive to the environmental conditions. Based on the control mode specificity, different operations can be carried out. Position control enables to get a better understanding of nanoscale interactions. Force control is used to deform molecular structures to measure their inherent properties. The performed experiments are very promising since they highlight the diversity of applications that can be handled by our system. The detailed analysis of the properties of the control modes can also easily be used to implement haptic feeback on other molecular simulators.

The stability conditions are not easily verified. Based on the experimental results, it can be seen however that force control remains stable for given values of the control scheme parameters. This is possible since its stability is not sensitive to environmental conditions. Future work include determining precisely how to choose the parameters to get a stable system. As position control stability is sensitive, non-constant gain might be a solution to get a stable system for a wide variety of environmental conditions. A user study will then be conducted to determine the benefits of the different haptic feedbacks from the two control modes. To enable the design of new molecular structures, our tool should offer the possibility to add and/or remove atoms. Haptic feedback would in this case assist the user while breaking the energy barriers. This functionality will be investigated.

Notations

Acknowledgment

This work was supported by the French National Agency of Research, through the PACMAN project.

Figure 1 :

 1 Figure 1: Representation of the tetra-alanine molecule. (a): The leaf nodes represent the user defined rigid bodies, (b): The binary assembly tree where the internal nodes represent the sub-assemblies, and the root node represents the complete molecular system.

Figure 2 :

 2 Figure 2: Notations related to the molecular simulator

Figure 3 :

 3 Figure 3: DFF control scheme

Figure 4 :

 4 Figure 4: FF control scheme

5 .

 5 Compute the sum F m (k + 1) + 1) of the forces applied to the molecule, and send it to the operator.

Figure 5 :

 5 Figure 5: Highest magnitude of the roots of the characteristic equation r max with respect to the stiffness of the interaction k e and to the computation time T s

Figure 6 :

 6 Figure 6: The HIV molecule and its inhibitors, 848 degrees of freedom. A haptic device is used to control the inhibitors and to feel the interactions with the HIV molecule.

4. 2 .

 2 The experiment consists in extracting an inhibitor of the molecule, and then reinserting it. All the molecules are flexible. Table 3: Parameters of the simulator and haptic coupling schemes for the HIV experimentSimulator -Mass of the manipulated molecule: m = 7.5 • 10 -25 kg -Cut off distance: d=8Å -Integration time for the equation of motion: T m = 2 • 10 -15 s. Haptic coupling -Force scaling coefficient: A f = 0.3 • 10 9 -Displacement scaling coefficient: A d = 0.25 • 10 9 -Stiffness of the FF coupling: k i = 100N.m -1

Figure 7 :3Figure 8 :

 78 Figure 7: Force felt by the user while manipulating an HIV's inhibitor (interaction forces between molecules). 172 dof activated -T s = 92ms

Figure 9 :

 9 Figure 9: Forces felt by the user while unfolding an inhibitor of the HIV virus (A), keeping the molecule stretched (B), and refolding it (C and D). The user feels the internal forces. 17 degrees of freedom, T s = 14ms, A d = 5 • 10 9 , A f = 0.1 • 10 9 , k i = 20N.m -1 .

Figure 10 :

 10 Figure 10: Manipulation of two water molecules. The black cross indicates the position of the equilibrium position. T s = 18ms, A d = 0.25 • 10 9 , A f = 1.5 • 10 9 .

F

  op /F h User/haptic force x h /v h Haptic position/velocity H(s)/H(z) Haptic device continuous/discrete transfer function M h /B h Haptic device inertia/viscosity Coupling. A f /A d Force/position scaling factor k i Stiffness of the FF coupling T s Sampling period T m Simulation integration period Molecule. x m /x j m Molecule/ j th rigid body actual position x i /x j i Molecule/ j th rigid body input position F m /F j m Force applied to the molecule/ j th rigid body F i /F j i User force applied to the molecule/ j th rigid body x l e /F l e Position/force of rigid bodies in the environment m/m p Molecule/p th rigid body inertia

Table 1 :

 1 Approximated values of Z op for the two control schemes

			Low frequencies	High frequencies
	DFF FF	Non-contact B h	Contact A d Z e ( jω) + B h A f k i k i +k e A f A d Z	Non-contact	Contact M h jω

e ( jω) + B h

Table 2 :

 2 Comparison of the DFF and FF coupling schemes

		Advantages	Drawbacks
			-Molecule's control not
	DFF Good transparency	optimal for the simulation -Stability sensitive to en-
		-Molecule's control	vironmental conditions
	FF	suited for the simulation -Good transparency during non-contact	Ideal transparency scaled by k i ke+k i during contact
		-Stability not sensitive to
		environmental conditions

For the sake of clarity only the translations are considered in this algorithm. Rotations are implemented using the same methodology. This remark also applies for

Algorithm 2.

http://www.sensable.com


Sarah Cohen-Boulakia

Christine Froidevaux

Jiuqiang Chen

Réécriture de workflows scientifiques et provenance

Keywords: workflows scientifiques, provenance, graphes série-parallèle

Les systèmes de workflow sont nombreux et disposent de modules de gestion de provenance qui collectent les informations relatives aux exécutions (données consommées et produites) permettant d'assurer la reproductibilité d'une expérience. Un grand nombre d'approches s'est développé pour aider à la gestion de ces masses de données de provenance. Un certain nombre de ces approches ont une bonne complexité parce qu'elles sont dédiées à des structures de workflows série-parallèles. Réécrire un workflow en un workflow série-parallèle permettrait donc de mieux exploiter l'ensemble des outils de provenance existants. Nos contributions sont : (i) introduction de la notion de réécriture de workflow provenance-equivalence, (ii) revue de transformations de graphes, (iii) conception de l'algorithme de réécriture SPFlow préservant la provenance (iv) évaluation de notre approche sur un millier de workflows.

Introduction

Les systèmes de gestion des workflows scientifiques (e.g. [1], [2], [START_REF] Foster | Chimera: A virtual data system for representing, querying, and automating data derivation[END_REF], [START_REF] Goecks | Galaxy: a comprehensive approach for supporting accessible, reproducible, and transparent computational research in the life sciences[END_REF], [START_REF] Gil | Wings for pegasus: Creating large-scale scientific applications using semantic representations of computational workflows[END_REF], [START_REF] Cohen-Boulakia | Selecting biomedical data sources according to user preferences[END_REF]) sont de plus en plus utilisés pour spécifier et gérer les expériences bioinformatiques. Dans ces systèmes, une expérience est représentée par un workflow dans lequel un grand nombre de tâches bioinformatiques (outils bioinformatiques, accès aux données) s'enchaînent les unes aux autres. Un workflow se représente alors naturellement sous la forme d'un graphe orienté où les sommets sont les tâches bioinformatiques à effectuer et les arcs représentent l'ordre dans lequel les tâches seront exécutées. La plupart des systèmes de workflows considèrent des graphes acycliques orientés : c'est ce que nous choisirons aussi dans cet article. Chaque workflow peut être exécuté plusieurs fois en faisant varier les données prises en entrée ou le paramétrage des outils (sans perte de généralité, les paramètres d'outils sont considérés comme des données). Nous considérerons ici le cas où une exécution se représente à son tour sous forme de graphe qui a la même structure que la spécification, où chaque 1 sommet représente l'exécution d'un outil et où les arcs sont étiquetés par les données consommées et produites à chaque étape. Aussi, nous considérons le cas de tâches reproductibles : deux exécutions d'une même tâche avec les mêmes données d'entrée fournira toujours le même résultat (en particulier, le comportement d'une tâche ne doit pas dépendre de variables aléatoires).

La Figure 1 propose (a) un exemple de workflow réel issu du système Taverna, (b) sa représentation sous forme de graphe et (c) un exemple de graphe d'exécution de ce workflow. Tout graphe (représentant un workflow ou une exécution) a une source unique s et un puits unique t, tels que tout sommet soit sur un chemin de s à t (accessibilité). Face à la complexité croissante des exécutions de workflows, au besoin de reproductibilité des résultats et au besoin d'associer une qualité aux données obtenues, la gestion de la provenance dans les workflows scientifiques est devenue un domaine de recherche en pleine expansion, en particulier dans la communauté Bases de Données. Un nombre important d'outils d'aide à la gestion des masses de données de provenance ont été conçus : certains proposent d'aider au stockage des données de provenance (réduction au maximum du volume de données nécessaire à la reproduction d'une expérience, indexation), d'autres permettent d'interroger ces données (pour mieux comprendre une différence entre deux exécutions, pour rechercher des motifs dans un ensemble (d'exécutions de) workflows disponibles), d'autres encore permettent de visualiser la provenance ou encore de (re)planifier les exécutions, un domaine particulièrement d'actualité avec l'avènement du cloud computing. Ces outils font tous intrinsèquement des opérations complexes sur des structures de graphes (recherche de sous-graphe dans un graphe, comparaison de graphes, ...), qui, si elles sont effectuées sur des graphes acycliques (DAGs) sans autre restriction de structure conduisent à des problèmes NP-difficiles. Au contraire, ces problèmes peuvent être résolus en temps polynomial lorsque certaines restrictions sont imposées aux graphes, comme par exemple qu'ils aient une structure série-parallèle (SP) [START_REF] Bein | Optimal reductions of two-terminal directed acyclic graphs[END_REF]. Certaines approches de gestion de provenance [START_REF] Bao | Differencing provenance in scientific workflows[END_REF][START_REF] Biton | Querying and managing provenance through user views in scientific workflows[END_REF][START_REF] Callahan | Vistrails: visualization meets data management[END_REF] ont donc choisi de se restreindre à des graphes de workflows ayant des structures SP. Néanmoins, en général, les workflows générés par les systèmes de workflows sont des DAGs de structure quelconque.

Proposer une procédure de réécriture de tout graphe de workflow en un graphe série-parallèle tout en assurant que la provenance du graphe transformé soit la même que dans le graphe original permettrait de mieux exploiter l'ensemble des outils de gestion de provenance existants. C'est l'objectif de cet article.

Nos contributions sont les suivantes : (i) nous introduisons la notion de réécriture de workflow préservant la provenance (section 3), (ii) nous déterminons parmi les stratégies de transformations de graphes existantes celles qui préservent la provenance et proposons de nouvelles stratégies (section 4), (iii) nous introduisons l'algorithme de réécriture de graphes SPFlow (section 5), (iv) nous fournissons des chiffres clés sur les résultats obtenus par notre approche sur un ensemble de workflows scientifiques réels et nous démontrons la faisabilité de notre approche et en particulier montrons que le coût de la réécriture est raisonnable (section 6). Nous discutons nos résultats et concluons en section 7.

Structures SP et workflows

Dans cette section, nous introduisons une formalisation des workflows et des exécutions sous forme de graphes et donnons les fondements sur les structures sérieparallèle. Nous présentons ensuite trois scénarios motivant l'intérêt de l'utilisation de ces structures pour les workflows scientifiques.

Workflows série-parallèle

Modèle de Workflow

Intuitivement, une spécification de workflow est un cadre pour des exécutions de workflows, qui spécifie l'ensemble des tâches qui sont exécutées et l'ordre qui doit être respecté entre les différentes exécutions des tâches. Formellement, nous modélisons une spécification de workflow comme un multigraphe orienté étiqueté, dont les sommets représentent les tâches du workflow et les arcs représentent le flot des données entre les tâches.

Considérons un multi-graphe acyclique G = (V, E), où V est un ensemble de sommets et E un ensemble d'arcs avec E ⊆ V × V .

On utilise le terme st-dag [START_REF] Bein | Optimal reductions of two-terminal directed acyclic graphs[END_REF] pour désigner un multi-graphe acyclique ayant deux sommets particuliers, s la source, et t le puits. On note s(G) (resp. t(G)) la source (resp. le puits) du st-dag G. Définition 2.1 [START_REF] Bao | Differencing provenance in scientific workflows[END_REF] Une spécification de workflow est un st-dag G spec = (V spec , E spec ) dont les sommets sont étiquetés par une fonction L vs : V spec → L V S , où L V S est un ensemble de valeurs de sommets.

Une spécification de workflow peut être exécutée plusieurs fois. Ainsi, selon les données en entrée du workflow et les assignations de valeurs aux paramètres des tâches, on obtiendra des exécutions de workflows différentes. Chaque exécution de workflow exécute les tâches et fait transiter entre elles des données.

Nous donnons ci-dessous la définition formelle d'une exécution de wokflow, qui est adaptée de celle de [START_REF] Bao | Differencing provenance in scientific workflows[END_REF] : Définition 2.2 Une exécution de workflow G run = (V run , E run ) est un st-dag dont les sommets et les arcs sont étiquetés, respectivement par une fonction L vr : V run → L V R , où L V R est un ensemble de valeurs de sommets, et une fonction L er : E run → L ER , où L ER est un ensemble de valeurs d'arcs. On notera x l'étiquette du sommet x, i.e. L vr (x) = x et d i l'étiquette de l'arc e i , i.e. L er (d i ) = e i . Une exécution de workflow G run = (V run , E run ) est valide par rapport à la spécification de workflow G spec = (V spec , E spec ) si :

(i) V run = V spec , (ii) E run = E spec et (iii) L vr = L vs avec L V R = L V S .
On introduit dans la sous-section suivante la notion de graphe série-parallèle.

Graphes série-parallèle

Les graphes série-parallèle ont été introduits par Duffin [START_REF] Duffin | Topology of series-parallel networks[END_REF] comme un modèle mathématique pour les réseaux électriques. Duffin a prouvé qu'un st-DAG est sérieparallèle si et seulement si il ne contenait pas un sous-graphe homomorphique au graphe "interdit" représenté en Figure 5(a). Des algorithmes pour détecter qu'un graphe est série-parallèle ont été proposés dont l'algorithme de Valdes et al. [START_REF] Valdes | The recognition of series parallel digraphs[END_REF] qui est en temps linéaire.

Un graphe série-parallèle (graphe SP) est un st-dag G que l'on obtient itérativement de la façon suivante : Un st-dag G, qui n'est pas un graphe SP est appelé un graphe non-SP. Déterminer qu'un graphe a ou non une structure SP se fonde sur un ensemble d'opérations de réduction qui ont été introduites dans [START_REF] Bein | Optimal reductions of two-terminal directed acyclic graphs[END_REF], reprises ci-après, et illustrées en Figure 2. Nous complétons ces définitions en spécifiant les étiquettes des sommets et des arcs. Les arcs introduits par les opérations de réduction sont étiquetés de sorte à garder la trace des sous-graphes qu'ils remplacent.

Définition 2.3 Soit G 1 = (V 1 , E 1 )
un st-dag dont les sommets et les arcs sont étiquetés, respectivement par une fonction

L 1vr : V 1 → L V R , et une fonction L 1er : E 1 → L ER . L'opération élémentaire op transforme G 1 en un st-dag op(G 1 ) = G 2 = (V 2 , E 2 )
, dont les sommets et les arcs sont étiquetés, respectivement par une fonction

L 2vr : V 2 → L V R , et une fonction L 2er : E 2 → (L V R ∪ L ER ) * , où (L V R ∪ L ER ) * désigne l'ensemble des expressions régulières construites sur L V R ∪ L ER .
1. Réduction Série Soient u, v, w trois sommets de V 1 tels que e = (u, v) est l'unique arc entrant de v et f = (v, w) est l'unique arc sortant de v. L'opération op de réduction en série en v remplace e et f par g = (u, w). Le graphe G 2 = (V 2 , E 2 ) obtenu à partir de G 1 est tel que :

V 2 ⊂ V 1 ,
V 2 = V 1 , L 2vr = L 1vr , L 2er identique à L 1er sur E 1 ∩ E 2 , et L 2er (g) = L 1vr (e 1 ) + ... + L 1vr (e k ).
3. Réduction de sommet Soit v un sommet de V 1 ayant un unique arc entrant e = (u, v) et ayant k arcs sortants f 1 = (v, w 1 ), ..., f k = (v, w k ). L'opération op de réduction de sommet en v remplace les k arcs {e, f 1 , ..., f k } par les k nouveaux arcs {g 1 , ..., g k } où g i = (u, w i ), 1. On appelle MaxRed(G) le graphe obtenu par application successive de toutes les opérations de réductions série et parallèle possibles sur G.

pour i ∈ [1, k]. Le graphe G 2 = (V 2 , E 2 ) obtenu à partir de G 1 est tel que : V 2 ⊂ V 1 , L 2vr est la restriction de L 1vr sur V 2 , L 2er est identique à L 1er sur E 1 ∩ E 2 , et L 2er (g i ) = L 1er (f i ) • L 1vr (v) • L 1er (e).

Un st-dag

G est dit maximalement réduit si G = M axRed(G). 3. G est SP si et seulement si M axRed(G) = BSP .
La Figure 3 décrit les opérations de réductions effectuées depuis le graphe de la Figure 1(b) jusqu'au graphe BSP. Aucune réduction de sommet n'a été utilisée : le graphe initial G 0 est SP. L'opération de réduction de sommet permet de donner une idée de la distance qui sépare un graphe non SP d'une structure SP. Intuitivement plus cette opération est utilisée (i.e. plus il y a de sommets de réduction) plus le graphe est loin d'une structure SP.

Dans la suite de cette section, nous présentons trois scénarios pour mettre en évidence l'intérêt de considérer des structures SP dans les workflows scientifiques.

Scénarios

Dans cette section, nous proposons trois scénarios pour illustrer les avantages à travailler sur des structures SP plutôt que sur des graphes acycliques quelconques.

Conception de Workflows

Alors que les workflows scientifiques ont pour objectif d'aider au partage et à la réutilisation de procédures d'analyses, une étude récente [START_REF] Starlinger | use in public scientific workflow repositories[END_REF] a montré que les auteurs de workflows réutilisaient facilement leurs propres workflow mais très rarement les workflows d'une tierce personne. Une des explications est que la structure de graphe d'un workflow scientifique peut être particulièrement complexe, rendant les principales étapes de l'analyse difficiles à appréhender. Guider les concepteurs de workflows à construire des workflows simples à comprendre est fondamentalement important pour améliorer le partage et la réutilisation de workflows.

Nous pensons que les structures SP peuvent être d'une aide précieuse dans ce contexte. D'abord intuitivement, et d'un point de vue purement visuel, les structures SP sont simples; les graphes SP sont en étage, leurs arcs ne se croisent pas, rendant les phases principales du workflow plus faciles à distinguer. Par exemple, considérons le workflow de la Figure 4

Interrogation de workflows

Une autre façon de concevoir un workflow est de se baser sur des workflows existants. L'utilisateur peut donc être amené à interroger les entrepôts de workflows pour rechercher un workflow ayant une structure particulière ou contenant un motif donné.

Le besoin de pouvoir faire ce type de recherche dans les entrepôts de workflow est exprimé depuis plusieurs années [START_REF] Goderis | Benchmarking workflow discovery: a case study from bioinformatics[END_REF][START_REF] Cohen-Boulakia | Search, adapt, and reuse: the future of scientific workflows[END_REF][START_REF] Gater | A graph-based approach for semantic process model discovery[END_REF] mais n'est aujourd'hui toujours pas pris en compte dans les entrepôts de workflows où les workflows sont interrogeables uniquement par mots clés présents dans leur description textuelle. La raison en est simple : ce type de recherche est directement associé à des problèmes connus pour être NP-difficiles (isomorphisme de sous-graphes) sur des DAGs classiques. Les structures SP présentent un avantage clair : le problème de recherche d'un sousgraphe isomorphe à un graphe se traite en temps polynomial sur de telles structures.

Un autre exemple de type de requêtes faisant appel à des opérations sur les graphes est la recherche de différences entre structures de workflows ayant des points communs. Là encore le problème du calcul de la différence de deux sous-graphes d'un même graphe est NP-difficile dans le cas général et polynomial pour des graphes SP.

Planification de workflows

De façon orthogonale, les structures SP peuvent être aussi particulièrement intéressantes dans le contexte de la planification d'exécutions de workflows. Dans le domaine plus général de la planification de tâches dans les programmes, les structures SP sont exploitées depuis des décennies [START_REF] González-Escribano | Performance implications of synchronization structure in parallel programming[END_REF], notamment parce qu'elles ont montré leurs avantages pour l'analyse de programme [START_REF] Lodaya | Series-parallel posets: Algebra, automata and languages[END_REF], l'estimation de coût [START_REF] Van Gemund | The importance of synchronization structure in parallel program optimization[END_REF], et l'efficacité de la planification [START_REF] Finta | Scheduling uet-uct series-parallel graphs on two processors[END_REF].

Avec le développement des grilles de calcul et plus récemment du cloud computing, la planification de tâches de workflow sur des ressources multiples et distribuées est d'importance croissante. Dans l'exemple simple de la Figure 4 les tâches composites A, B, C, D peuvent être ordonnancées dans cet ordre dans la structure modifiée SP, alors qu'en considérant le graphe d'origine non SP il y a une dépendance en plus à considérer entre les tâches A et D, rendant plus complexe le découpage du workflow pour la planification.

Provenance

Modèle de provenance

La notion de provenance d'une donnée produite par une tâche fixée a été introduite dans [START_REF] Cohen | Towards a model of provenance and user views in scientific workflows[END_REF] pour une exécution de workflow. Nous reprenons cette notion et la formalisons ici. Intuitivement, la provenance d'une donnée est la suite ordonnée des tâches exécutées pour produire cette donnée (c-à-d, les tâches dont elle dépend), et les données d'entrée de ces tâches. En ce sens, notre modèle est naturellement compatible avec des approches suivant OPM (Open Provenance Model1 ). Nous distinguons la provenance immédiate, qui décrit la dernière étape de production d'une donnée et la provenance profonde qui donne la succession des étapes ayant permis de produire la donnée. Rappelons que dans la représentation sous forme de multi-graphe étiqueté d'une exécution de workflow, une donnée produite par une tâche est l'étiquette d'un arc sortant du sommet représentant cette tâche. Formellement, on a la définition suivante : Définition 3.1 Provenance d'une donnée Etant donné une exécution de workflow G run = (V run , E run ) soit :

• u ∈ V run qui n'est pas la source de (G run ), avec L vr (u) = ũ ;

• f ∈ E run un arc sortant de u avec L er (f ) = d; • e i ∈ E run , 1 ≤ i ≤ p les arcs entrants de u, avec L er (e i ) = d i .
La Provenance immédiate de d dans G run (où d est l'étiquette de l'arc f ) est définie par la fonction imP rov :

E run → (L V R ∪ L ER ) * , avec : imP rov(f ) Grun = ũ • (d 1 + . . . + d p ) (1)
La Provenance Profonde de d dans G run (où d est l'étiquette de l'arc f ) est définie comme la suite ordonnée des tâches et des données d'entrée des tâches qui ont été successivement utilisées pour la produire. Elle est formalisée par la fonction DP rov :

E run → (L V R ∪ L ER ) * , qui est définie récursivement comme suit : DP rov(f ) Grun = ũ • (d 1 • DP rov(e 1 ) Grun + . . . + d p • DP rov(e p ) Grun ) (2)
Le cas de base est donné lorsque u égal à s, la source de G run (st-dag), avec f ∈ E run un arc sortant de s, s = L vr (s) et L er (f ) = d. La provenance de d est alors donnée par :

DP rov(f

) Grun = imP rov(f ) Grun = s (3)
Les expressions régulières permettent d'exprimer de façon concise le fait qu'une donnée est produite par une tâche à partir de l'utilisation conjointe de l'ensemble de données d'entrée de cette tâche (opération "+") et le fait que les tâches en aval consomment les données produites par les taches en amont (opération "•").

Exemple 3.1 Considérons le graphe G r de la Figure 5 (a). La provenance immédiate de la donnée d 4 transitant sur l'arc e 4 est donnée par la tâche qui l'a produite, v et les données directement utilisées en entrée de cette tâche pour la produire, soit d 3 et d 1 . Cette information s'écrit de la façon suivante : imP rov(e 4 ) Gr = ṽ

• (d 1 + d 3 ).
On a aussi : DP rov(e 4 ) Gr = ṽ

• (d 1 • DP rov(e 1 ) Gr + d 3 • DP rov(e 3 ) Gr ) = ṽ • (d 1 • s + d 3 • [ũ • d 2 • s]) = ṽ • [d 1 + (d 3 • ũ • d 2 )] • s.
•ṽ•[d 1 +(d 3 •ũ•d 2 )]•s peut aussi s'écrire d 4 •ṽ•d 1 •s+d 4 •ṽ•d 3 •ũ•d 2 •s.
Dans cette deuxième formulation, on écrit jusqu'à la source du graphe tous les chemins qui ont produit les données d'entrées à la dernière tâche ayant produit d 4 plutôt que d'avoir une formulation factorisée montrant l'aspect récursif du processus. D'autre part, étant donnée une exécution de workflow G run et un sommet u, la provenance de tous les arcs sortants de u est la même. En d'autres termes, toutes les sorties d'une même tâche ont la même provenance, ce qui exprime bien qu'elles ont été (récursivement) obtenues de la même façon.

Dans certains cas d'étude de la provenance d'une donnée, il peut être intéressant de connaître simplement les données impliquées dans la production d'une donnée sans se soucier de l'ordre dans lequel elles ont été consommées ni des tâches qui ont été exécutées. Pour cela, nous introduisons la définition suivante de U sef ul qui fournit toutes les données utilisées pour produire une donnée fixée.

Définition 3.2 Données utilisées

Étant donnée une exécution de workflow G run = (V run , E run ) soit d l'étiquette de l'arc f dans G run . Les données utilisées pour produire d sont fournies par U sef ul(f ) qui est l'ensemble de toutes les étiquettes d i apparaissant dans DP rov(f ) Grun . On dira que d dépend de d i .

Il est important de noter que U sef ul ne suffit pas à montrer que deux provenances sont équivalentes (puisqu'on perd -en plus des tâches -à la fois l'ordre d'utilisation des données et le nombre de fois qu'elles peuvent être utilisées). Néanmoins, si les données dépendantes d'une donnée d sont différentes alors les provenances de d ne sont pas équivalentes. Nous introduisons maintenant la notion d'histoire d'une tâche dans une exécution de workflow qui va nous permettre de définir ensuite la notion fondamentale de provenance d'une sortie d'exécution. L'histoire d'une tâche est liée à la notion de provenance, dans la mesure où la provenance d'une donnée produite par une tâche est obtenue par la concaténation de l'étiquette de la tâche et de l'histoire de la tâche. Formellement : Définition 3.3 Histoire Soit G run = (V run , E run ) une exécution de workflow. L'histoire de u dans G run est donnée par une fonction Hist : V run → (L V R ∪ L ER ) * , qui est définie comme suit :

• si u=s(G run ), Hist(s) Grun = ε (mot vide) • si u=t(G run ), soient e i ∈ E run , les arcs entrants de t (1 ≤ i ≤ p), avec L er (e i ) =
G r ) = [d 4 • ṽ • [d 1 + (d 3 • ũ • d 2 )] • s] + [d 5 • ũ • d 2
• s] qui exprime que la provenance de sortie du graphe est l'histoire de t c'est-à-dire dans cet exemple, la provenance des données d 4 et d 5 .

Nous disposons maintenant de l'ensemble des notions nécessaires pour définir la provenance-équivalence de deux exécutions qui fait l'objet de la sous section suivante.

Conservation de la Provenance

Nous cherchons à transformer une exécution de workflow non SP en une autre qui soit SP et qui soit le plus proche possible du graphe original. Pour cela nous voulons conserver la provenance de sortie du graphe. Nous introduisons donc ici la notion d'équivalence des provenances de sortie de deux graphes.

OutP rov(G ′ r ) = [d 4 • ṽ • [d 1 + (d 3 • ũ • d 2 )] • s] + [d 5 • ũ • d 2 • s] qui n'est rien d'autre que OutP rov(G r ). D'où : G r prov ⇔ G ′ r .
Propriété 3.1 Conservation de la provenance des opérations de réduction : Les opérations de réduction conservent la provenance. Plus précisément, en reprenant les notations de la définition 2.3 et en considérant la Figure 2 on a :

1. Réduction série :

Hist(w) G 1 = Hist(w) G 2 ;
2. Réduction parallèle : Hist(w) G 1 = Hist(w) G 2 ;

3. Réduction de sommet :

Hist(w i ) G 1 = Hist(w i ) G 2 pour tout i ∈ [1, k].
Cette propriété résulte du fait que nous stockons dans l'étiquette des arcs restants,dans l'ordre inverse du flot de données, les étiquettes des arcs et des sommets qui ont été supprimés.

4 Réécriture de graphes

Dans cette section nous passons d'abord en revue les principales stratégies qui existent dans la littérature pour réécrire un graphe non-SP en un graphe SP. Pour chacune d'elles nous déterminons si elles préservent la provenance ou pas. Dans un second temps, nous introduisons notre stratégie de réécriture de graphe, fondée sur la duplication de sommet en précisant dans quel contexte cette stratégie préserve la provenance. 4.2 Notre approche : Duplication de sommets L'approche que nous proposons se fonde sur une opération de duplication de sommet. Nous résolvons le problème des sommets de réduction en proposant de les dupliquer ce qui permet de garder le même ensemble de tâches. Le workflow de la Figure 4 de la section 2 a été transformé en suivant notre procédure. Deux opérations de duplication de sommets sont possibles et représentées en Figure 7, chacune fondée sur l'une des deux opérations de réduction de sommet. Nous définissons formellement l'opération de duplication de type (a) (fondée sur les entrées de la tâche à dupliquer).

(Re)synchronisation

Définition 4.1 Soit G 1 = (V 1 , E 1 ) un st-dag dont les sommets et les arcs sont éti- quetés par L 1vr et L 1er , respectivement. L'opération op transforme G 1 en un st-dag op(G 1 ) = G 2 = (V 2 , E 2 ), dont les sommets et les arcs sont étiquetés, respectivement par une fonction L 2vr : V 2 → L V R , et une fonction L 2er : E 1 → (L V R ∪ L ER ) * .
Duplication de sommet: Soit v un sommet de V 1 ayant un unique arc entrant e = (u, v) et ayant k arcs sortants f 1 = (v, w 1 ), ..., f k = (v, w k ). L'opération op de duplication du sommet v ajoute de nouveaux sommets v 1 ,..., v k-1 , qui sont des copies du sommet v, ajoute les arcs {e 1 , ..., e k-1 }, avec

e i = (u, v i ) pour i ∈ [1, k -1], et remplace les arcs {f 2 , ..., f k } par de nouveaux arcs {g 2 , ..., g k } avec g i = (v i-1 , w i ) pour i ∈ [2, k]. Le graphe G 2 = (V 2 , E 2 ) = op(G 1 ) est tel que : V 2 ⊃ V 1 , L 2vr est une extension de L 1vr sur V 2 , qui coïncide avec L 1vr sur V 1 ∩ V 2 avec L 2vr (v i ) = L 1vr (v) pour tout i ∈ [1, k -1]. De plus, L 2er est identique à L 1er sur E 1 ∩ E 2 , L 2er (e i-1 ) = L 1er (e), et L 2er (g i ) = L 1er (f i ) pour i ∈ [2, k].
Propriété 4.1 Conservation de la provenance par l'opération de duplication de sommet. L'opération de duplication de sommet conserve la provenance. Plus précisément, en reprenant les notations de la définition ci-dessus on a :

Hist(w i ) G 1 = Hist(w i ) G 2 pour tout i ∈ [1, k].
5 Algorithme de transformation

L'algorithme SPFlow

Dans cette section nous introduisons l'algorithme SPFlow qui réécrit un graphe G non-SP en un nouveau graphe SP, appelé SP G, obtenu à partir de G en dupliquant certains sommets de G, tout en s'assurant que G et SP G sont provenanceéquivalents. Nous allons appliquer successivement des réductions série et parallèle, et effectuer des réductions de sommet sur G construisant ainsi un graphe intermédiaire G red , et ce, jusqu'à ce que G red soit le graphe basique BSP . Les sommets à dupliquer sont déterminés par les sommets à réduire dans G red . Tout d'abord, nous rappelons la notion de sous-graphe autonome, présentée dans [START_REF] Bein | Optimal reductions of two-terminal directed acyclic graphs[END_REF] et utilisée dans notre algorithme. Intuitivement, les sous-graphes autonomes permettent de se limiter à des composantes plus petites du graphe initial pour effectuer les duplications de sommets, sans interaction avec le reste du graphe, puisqu'aucun arc n'entre dans le sous-graphe autonome ou n'en sort. Définition 5.1 Sous-graphe autonome Soit G un st-dag. Nous notons G[(v, w) le sous-graphe de G de source v et de puits w. On dit que G(v, w) est un sous-graphe autonome de G s'il satisfait la propriété suivante : Pour tout chemin P de s vers t dans G, l'ensemble des arcs de P ∩G[v, w] est vide ou forme un chemin de v vers w. Notons que v peut être s, et que w peut être t ou les deux, mais que G(v, w) ne peut pas être un arc unique ni tout le graphe G. Si G(v, w) est un sous-graphe autonome de G, nous appelons (v,w) une paire de coupure de G.

Une décomposition de G en sous-graphes autonomes peut être obtenue en temps linéaire [START_REF] Bein | Optimal reductions of two-terminal directed acyclic graphs[END_REF]. Dans la suite, nous notons G[v, w] le sous graphe de G qui contient tous les noeuds et arcs de tous les chemins de v ) w dans G. successives de G (y compris des réductions de sommets). Dans SP G certains sommets de G sont dupliqués; ces sommets sont déterminés à partir des réductions de sommet effectuées dans G red .

Initialisation : (i) SP G ← G; s ← s(G); t ← t(G) (ii) G red ← M axRed(G) (iii)

Expérimentations

Les données utilisées dans cette section ont été collectées dans myExperiment2 , l'entrepôt de workflows scientifique publique le plus important. Nous avons focalisé notre étude sur les workflows conçus dans le cadre du système Taverna, qui est l'un des systèmes de workflows les plus utilisés aujourd'hui et qui est le principal système dont les workflows de myExperiment sont issus. En Mai 2012, 1421 workflows de Taverna étaient présents dans la base. Nous avons retiré de cet ensemble les workflows mal formés ainsi que les doublons; le nombre total de workflows ainsi sélectionnés est alors de 1014 wokflows. La présente section a deux objectifs : (1) présenter l'étude que nous avons faite sur la structure de workflows réels, (2) fournir des chiffres clés relatifs à l'application de notre approche sur ces workflows. [START_REF] Wang | Study of the managment of the provenance information in taverna[END_REF]. Nous avons implémenté l'algorithme de détection de structures SP dont nous avons donné un exemple de déroulement dans la section 2, Figure 4. Nous indiquons ci-après des informations sur (i) la proportion de structures non-SP vs SP dans les workflows, (ii) l'évolution de cette proportion et du nombre de workflows dans les dernières années, (iii) la distance qui sépare les structures non-SP des structures SP et (iv) des caractéristiques des workflows non-SP.

Proportion de workflows SP et non-SP

Notre premier résultat donné en Table 1 montre qu'il y avait une majorité de structures SP dans l'ensemble des workflows disponibles lorsque myExperiment a débuté et que la proportion de graphes SP est stable. Plus précisément nous avons étudié la répartition des workflows SP et non-SP sur quatre familles de workflows définies en fonction de leur nombre de sommets : les workflows simples ont de 1 à 3 sommets, les intermédiaires de 4 à 10, les complexes de 11 à 20, et les très complexes strictement plus de 20 sommets.

Sans surprise, on trouve moins de structures SP parmi les workflows les plus complexes. Les chiffres obtenus sont particulièrement tranchés : alors que les workflows simples ou intermédiaires sont quasiment tous SP, la proportion de structures de workflows SP chute pour les workflows de plus de 10 sommets : 55% de structures non-SP pour les workflows complexes et 93% pour les workflows ayant plus de 20 sommets.

Dans la suite, nous étudions la distance qui sépare les structures non-SP des structures SP en nous basant sur le nombre d'opérations de réduction de sommet qu'il faut appliquer au graphe de départ pour obtenir le graphe basique BSP [START_REF] Bein | Optimal reductions of two-terminal directed acyclic graphs[END_REF].

Caractéristiques des workflows non-SP

Dans ce second ensemble d'expériences, notre objectif est d'évaluer à quel point les structures non SP sont éloignées des structures SP. La Figure 9 présente le pourcentage de workflows ayant un nombre donné de sommets de réduction dans les structures non-SP. Par exemple, on peut lire que 31% des workflows non SP ont un sommet de réduction alors que 27% en ont deux.

Le principal résultat donné par la Figure 9 est que 67% des workflows non SP ont seulement de 1 à 3 sommets de réduction.

Evaluation de SPFlow

Dans cette sous-section, nous évaluons le comportement de l'algorithme SPFlow sur les données réelles. Nous travaillons ici sur un ensemble de 390 workflows, puisque nous ne considérons que les graphes non-SP que nous réécrivons en graphes SP.

La première information est celle du rapport entre la taille du graphe initial (en nombre de sommets) et la taille du graphe réécrit. La Figure 10 fournit cette information. Bien que le nombre de sommets de réduction ne soit pas directement associé à la taille (en nombre de sommets) d'un graphe, on peut constater que les workflows de moins de 30 sommets ont au pire une taille qui triple lors de l'opération de réécriture. Même si trois graphes ont un nombre de sommets dupliqués qui peut être très important, d'autres très gros graphes (plus de 300 sommets) peuvent avoir un ratio faible.

Enfin, on peut noter que le temps de réécriture est négligeable pour les structures actuelles de workflows (qui ont relativement peu de sommets de réduction). Sur une machine dual core @2.2GHz et 2GB de RAM, le temps minimal d'exécution est de 1 ms, le temps maximum de 434 ms. Même dans le cas de quelques workflows très complexes, ce temps est très raisonnable.

Figure 1 :

 1 Figure 1: Exemple de workflow simple issu du système Taverna; (b) représentation sous forme de graphe; (c) exemple de graphe d'exécution

  Le sommet v est appelé sommet de réduction. Le cas où v a un degré sortant égal à 1 et un degré entrant strictement supérieur à 1 est défini de façon analogue (cf. 2(d)).

Figure 2 :

 2 Figure 2: Opérations de réductions : (a) réduction série; (b) réduction parallèle, (c) réduction de sommet fondée sur les sorties; (d) réduction de sommet fondée sur les entrées

Figure 3 :

 3 Figure 3: Exemple d'opérations de réduction appliquées successivement à G 0 .

  (a) dont la structure représentée en Figure4(b) n'est pas SP. Il est relativement complexe de distinguer visuellement les grandes étapes de ce workflow ou de construire des sous-workflows indépendants. La Figure4(c) représente le même workflow dans lequel le sommet 2 a été dupliqué en 2' (ce qui revient à proposer d'interroger deux fois le catalogue de données). Le workflow restructuré est SP et déjà plus facile à appréhender alors que sa structure originale n'est pas très éloignée d'un graphe SP puisque seul le sommet 2 est un sommet de réduction.

Figure 4 :

 4 Figure 4: (a) Exemple de workflow simple issu du système Taverna, (b) structure graphe du workflow (non SP); (c) transformation de la structure en graphe SP; (d) proposition de sommets composites; (e) graphe de haut niveau obtenu

Figure 5 :

 5 Figure 5: Deux graphes pour illustrer les notions relatives à la provenance

Exemple 3 . 2 Définition 3 . 4 Exemple 3 . 3

 323433 d i : Hist(t) Grun = d 1 • DP rov(e 1 ) Grun + . . . + d p • DP rov(e p ) Grun • pour tout sommet u ∈ V run tel que u = s(G run ) et u = t(G run ) avec f ∈ E run un arc sortant de u, Hist(u) Grun est tel que DP rov(f ) Grun = ũ • Hist(u) Grun . Considérons à nouveau le graphe G r de la Figure 5 (a), l'histoire de u est la provenance des données prises en entrée par u, soit uniquement d2 : Hist(u) Gr = d 2 • s. Provenance de sortie d'un workflow Etant donnée une exécution de workflow G run , sa provenance de sortie, notée OutP rov(G run ), est définie comme l'histoire du sommet puits du graphe. Formellement : OutP rov(G run ) = Hist(t(G run )) Grun . Toujours sur la Figure 5 (a), on a OutP rov(

Définition 3 . 5 Exemple 3 . 4 5 .

 35345 Provenance-équivalence Soit G r1 , G r2 deux exécutions de workflows. Nous disons que G r1 et G r2 sont provenance-équivalents si et seulement s'ils ont même provenance de sortie, i.e. OutP rov(G r1 ) = OutP rov(G r2 ). Nous utiliserons la notation G r1 prov ⇔ G r2 dans ce qui suit. Reprenons la Figure 5 et considérons les graphes G r (a) et G ′ r (b). Le graphe G ′ r a été obtenu à partir de G r en dupliquant le sommet u en un sommet u ′ : ils ont donc la même étiquette. De même, les arcs e 2 et e ′ 2 ont même étiquette, ainsi que les arcs e 5 et e ′ OutP rov(G ′ r ) = [L er (e 4 ) • L vr (v) • [L er (e 1 ) + (L er (e 3 ) • L vr (u) • L er (e 2 ))] • L vr (s)] + [L er (e ′ 5 ) • L vr (u ′ ) • L er (e 2 ) • L vr (s)] Or L vr (u ′ ) = L vr (u), L er (e 5 ) = L er (e ′ 5 ) et L er (e 2 ) L er (e 2 ). On a donc :

  Une stratégie pour rendre SP des structures non SP a été proposée par Escribano[START_REF] González-Escribano | Synchronization architecture in parallel programming models[END_REF][START_REF] Escribano | Conversion from nsp to sp graphs[END_REF] en se fondant sur la notion de (re)synchronisation. De façon informelle, l'idée est "d'étager" le graphe en ajoutant des sommets (et des arcs) artificiels, qui jouent le rôle de tâches de synchronisation. Ces algorithmes combinent trois stratégies principales de synchronisation illustrées en Figure6. À partir du graphe représenté en (a), la stratégie haute (up-synchronisation) fournit le graphe (b), la stratégie basse (down-synchronisation) fournit le graphe (c), et la stratégie croisée (across-synchronisation) le graphe (d). Dans la suite nous donnons une intuition de la raison pour laquelle chacune de ces stratégies transforme un sous-graphe non SP en un sous-graphe SP et préserve ou non la provenance. Ces transformations rendent-elles le graphe final SP ? Oui. Considérons la synchronisation haute (graphe 2.(b)). Une réduction parallèle supprime les doubles arcs entre u et v d'une part et entre v et t d'autre part. Une réduction série supprime le sommet u qui n'a plus qu'un arc entrant et un arc sortant suite à la réduction précédente; une réduction parallèle et une réduction série permettent enfin d'obtenir le graphe BSP . Le même type de raisonnement peut se faire sur les cas (c) et (d). Ces transformations conservent-elles la provenance ? Non. Si la provenance était conservée, aucune dépendance de données ne serait ajoutée ou perdue. Or, dans le graphe (b), la donnée d5 dépend de la donnée d 4 (d 4 ∈ U sef ul(d 5 )), ce qui n'est pas le cas dans le graphe (a) (d 4 / ∈ U sef ul(d 5 )). Dans le graphe (c), c'est la donnée d 3 qui dépend de d 2 ce qui n'est pas le cas dans le graphe (a). Enfin, dans le cas (d), les données en sortie du nouveau sommet de synchronisation x dépendent de toutes les données entrantes de ce sommet, en particulier d 5 dépend de d 1 , ce qui n'est pas le cas dans le graphe (a).

Figure 6 :

 6 Figure 6: Resynchronization: À partir du graphe interdit en (a), les opérations de up-, down-et across-synchronisation fournissent respectivement les graphes (b), (c) et (d).

Figure 7 :

 7 Figure 7: Duplication de sommets (a) à partir des entrées, (b) à partir des sorties.

Exemple 5 . 1

 51 Le graphe G de la Figure 8 a plusieurs sous-graphes autonomes. Par exemple, G[s, x] et G[x, t] dont les paires de coupure sont (s, x) et (x, t).

  Décomposer G red en sous-graphes autonomes. (iv) Appeler T ransf oRec(G, G red , SP G, s, t). Procédure TransfoRec (IN: G; IN/OUT: G red ; IN/OUT: SP G; IN: u, p) Tant que G red = BSP faire : Etape 1 : Choisir dans G red un sommet v successeur de u dans G red sur lequel on puisse appliquer une réduction de sommet selon Figure 2(c) (cf. Section 2). Etape 2: (i) v est la source d'un sous-graphe autonome de puits w ∈ G red . Rappeler TransfoRec (G, G red , SP G, v,w), en considérant G red [v, w] au lieu de G red , c-a-d, en considérant v comme la nouvelle source du graphe réduit. (ii) v n'est pas la source d'un sous-graphe autonome. 1. Dupliquer dans SP G le sommet v k -1 fois, si v a k successeurs dans G red . Lors de cette duplication, au lieu de simplement dupliquer l'arc e de G red , dupliquer le sous-graphe SP G[u, v] en SP G[u, v 1 ], ..., SP G[u, v k-1 ]. De même, au lieu de simplement considérer les arcs f i de G red , considérer les sous-graphes SP G[v, w i ], qui deviennent SP G[v 1 , w 2 ], ..., SP G[v k-1 , w k ].2. Appliquer dans G red l'opération de réduction de sommet à v.

3 .

 3 G red ← M axRed(G red )Fin Tant que Fin TransfoRec

Figure 8 :

 8 Figure 8: Exemple d'une étape d'exécution de TransfoRec

Figure 9 :

 9 Figure 9: Pourcentage des workflows ayant un nombre fixé de sommets de réductions

Figure 10 :

 10 Figure 10: Nombre de sommets dupliqués dans le graphe réécrit par rapport au nombre de sommets dans le graphe initial

  Composition Série : étant donnés G 1 (avec une source s 1 et un puits t 1 ) et G 2 (avec une source s 2 et un puits t 2 ) deux graphes SP, G est obtenu en identifiant s = s 1 , s 2 = t 1 et t = t 2 ; 3. Composition Parallèle : étant donnés G 1 (avec une source s 1 et un puits t 1 ) et G 2 (avec une source s 2 et un puits t 2 ) deux graphes SP, G est obtenu en identifiant s = s 1 = s 2 et t = t 1 = t 2 .

	1. Graphe SP Basique : G, le graphe qui contient un unique arc entre s et t, est
	un graphe SP (appelé "BSP" pour "basic SP-graph");
	2.

  Réduction Parallèle Soient v et w deux sommets de V 1 reliés par k arcs e 1 ... e k . L'opération op de réduction en parallèle en v et w remplace les k arcs par un unique arc g = (v, w) et laisse tous les autres arcs inchangés. Le graphe G 2 = (V 2 , E 2 ) obtenu à partir de G 1 est tel que :

	2.
	L 2vr est
	la restriction de L 1vr sur V 2 , L 2er est identique à L 1er sur les arcs communs
	aux deux graphes, et L 2er (g) = L 1er (f ) • L 1vr (v) • L 1er (e).

  Cette formule exprime que la donnée transitant sur l'arc e 4 a été obtenue par l'utilisation de la tâche v qui a pris pour entrées d 1 et d 3 dont on donne la provenance : d 1 ; qui est obtenue directement de la source s alors que la provenance de d 3 nécessite une étape de récursion supplémentaire : elle est obtenue en sortie de la tâche u qui a utilisé d 2 , avec d 2 obtenue par la source s.Remarques: Compte tenu des propriétés de commutativité, associativité, et distributivité de • et +, la provenance d'une donnée est définie de façon unique pour une exécution de workflow fixée. Intuitivement, la notation sous forme d'expression régulière permet de représenter l'exécution de façon concise. Que l'on cherche à factoriser ou développer l'expression, on exprime la même provenance. Par exemple, l'expression d 4

Table 1 :

 1 Évolution de la proportion des structures SP parmi les workflows scientifiques Chaque workflow a d'abord été représenté par un st-DAG, i.e. un DAG dans lequel nous avons ajouté une unique source pointant vers tous les sommets n'ayant pas de prédécesseurs et un unique puits sur lequel pointent tous les sommets n'ayant pas de successeurs

	de myExperiment
	Date Nombre de graphes SP graphes SP
	workflows (proportion) (proportion)
	2010	681	429 (63%)	252 (37%)
	2011	879	554 (63%)	325 (37%)
	2012	1014	624 (61,5%) 390 (38,5%)
	Table 2: Proportion des structures Non-SP vs SP par famille
	Famille	#workflows % de structures non-SP
	Simple		92	100 %
	Intermédiaire	480	87 %
	Complexe		214	45 %
	Très complexe	228	6,7 %
	6.1 Structure des workflows

http://openprovenance.org/

Principe de l'algorithme: L'algorithme SP F low prend en entrée le graphe G et donne en sortie le graphe G red et le graphe SP G. G red est obtenu par réductions

http://www.myexperiment.org