Jan Michalczyk

Maurice Brémond

Vincent Acary

Roger Pissard-Gibollet

Exomars Rover Mechanical Modeling with Siconos

Keywords: robotics, dynamical systems, nonsmooth dynamics, multibody systems, simulation Modélisation robotique, systèmes dynamiques, dynamique nonrégulière, systèmes multicorps, simulation

This document contains specification and documentation of the tests performed on the exomars planetary rover model using Siconos software. Siconos is an open source scientific software targeted at modeling and simulating nonsmooth dynamical systems. The advantage of using Siconos for mechanical simulations is that it offers efficient friction-contact nonsmooth models. Performing mechanical simulations of the exomars rover appears to be of great importance with respect to verifying its design and estimating energy consumption. Seven different scenarios have been defined in order to simulate the mechanical behaviour of the exomars rover. Dynamics and kinematics models of the exomars rover have been created using HuMAnS software and implemented inside Siconos using its C code generator as plugin functions. The intention of this report is to create a common reference and unification of mechanical tests of the exomars rover across different simulation platforms.

Introduction

This document contains the specification and documentation of tests performed on the ESA (European Space Agency) Exomars rover model using Siconos software. It can also be considered an expansion and complement of work outlined in [START_REF] Acary | Mechanical simulation of the exomars rover using siconos in 3drov[END_REF]. Siconos is an open source software code capable of dealing with idealized contact interaction based on unilateral (rigid) contact, 3D Coulomb friction and impacts. Using such software in planetary rover's simulations may help in designing control laws for the locomotion and simulating the global behavior of the system. Model of the system has been created using HuMAnS software and implemented using its C code generator. The intention thereof is a common reference and unification of mechanical tests of the Exomars rover across different platforms. The document is structured in such way that firstly some basics on the Nonsmooth Mechanics and Siconos are given, then motivations for using the Nonsmooth approach are laid down and finally the first chapter ends with two subsections on how the simulations have been set up in Siconos and how the model of the system has been obtained. Second chapter contains detailed test cases along with comments on each of them.

Nonsmooth Mechanics at a glance

In this section, we recall the main ingredients of the nonsmooth mechanics approach for the modeling and the simulation of mechanical systems with contact and friction. For more details, we refer to [START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF]. Let us start with the following Lagrangian formulation of the equations of motion

                 q = v, (1a)
M (q) v + N (q, v) + F (t, q, v, u) = G ⊤ (q)λ, (1b)
u = d(t, q, v, u, λ), (1c)
g k (q) = 0, k ∈ E (1d)

g k (q) ≥ 0, λ k ≥ 0, λ k g k (q) = 0 k ∈ I, (1e)
where q(t) ∈ IR n is a set of coordinates that gives in a unique way the configuration of the system. The vector v(t) ∈ IR n describes the velocity. The matrix M ∈ IR n×n is a symmetric, positive and (semi-)definite matrix representing the inertia. The vector of gyroscopic effects is denoted as N ∈ IR n . The vector of forces F ∈ IR n applied to the system includes the internal and external forces and the effect of the vector u ∈ IR nu which might be considered as a control input vector. The dynamics of this control can be described by another first order dynamical system (1c). The sets E ⊂ IN and I ⊂ IN are the index sets of bilateral constraints (E stand for equalities) and of unilateral constraints (I stand for inequalities), respectively. The matrix G(q) = ∇ ⊤ q g(q) is the Jacobian of the constraints. The condition (1e) is called a complementarity condition, or equivalently the Signorini condition. It can be equivalently written as 0 ≤ g k (q) ⊥ λ k ≥ 0 for all k ∈ I.

(

) 2
This condition models the basic unilateral contact condition. The complementarity comes from the fact that the reaction force at contact represented by λ can only act when the contact is closed, that is when g k (q) = 0.

In order to complete the model at contact, let us introduce the local contact variables. Usual kinematic relations between the local velocity U k ∈ IR 3 at one contact k and the generalized variables, and by duality, the generalized forces r due to contact forces R k ∈ IR 3 give

U k (t) = G k (q) q = G k (q)v, r = G ⊤ k (q)R k . (3)
Let us drop the subscript k to lighten the notation. By introducing a local frame at contact (N, T, S) where N is a unit normal vector, we can decompose the local variables as

U = U N N + U T , U T ∈ IR 2 , R = R N N + R T , R T ∈ IR 2 .
The Lagrangian dynamics in (1) is usually not so smooth. If the normal relative velocity at contact U N (t ⋆ ) is positive when the bodies hit at time t ⋆ , in other words, when an impact occurs, a velocity jump must occur to satisfy the constraints after t ⋆ and an impact law has to be added. Let us simply choose the Newton impact law

U + N (t ⋆ ) = -eU - N (t ⋆ ), for all t such that g(q(t ⋆ )) = 0. (4)
By formulating the unilateral contact in terms of local variables at contact we get the impact law at the velocity level

if g N ∆ = g(q) ≤ 0, then 0 ≤ U + N + eU - N ⊥ R N ≥ 0. (5)
Coulomb's friction can be expressed in a disjunctive form as

if g N ≤ 0,      if U T = 0 then R ∈ C if U T = 0 then R ∈ ∂C, and it exists a ≥ 0 such that R T = -aU T ( 6
)
where C is the Coulomb friction cone, C = {R, | R T ≤ µR N } with µ the coefficient of friction.

Coulomb's friction can also be formulated compactly as a complementarity condition on second-order cones as [START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF][START_REF] Acary | A formulation of the linear discrete Coulomb friction problem via convex optimization[END_REF][START_REF] Acary | Recent Advances in Contact Mechanics[END_REF]]

- Û ∆ = - U N + µ U T U T and C * ∈ Û ⊥ R ∈ C (7)
where C * is the dual cone of C [START_REF] Rockafellar | Convex Analysis[END_REF].

Siconos

SICONOS is an Open Source scientific software primarily targeted at modeling and simulating nonsmooth dynamical systems in the most generic sense. In particular, it can deal with the following systems: The developers team focused on the development of the simulation and modeling tools which are as reusable as possible for different scientific contexts. Any specific pre/post-processing tools has been developed. Dynamical systems and interactions (unilateral constraints) classes provided by Siconos/Kernel may be configured with the help of plugins or by class derivation. For high-level applications SICONOS relies on external software and provides with an efficient computational engine. Concerning the mechanical modeling, it is possible to simulate linear or non-linear rigid or deformable bodies.

The numerical simulation strategies available in Siconos/Kernel are quasi-statics, smooth dynamics, nonsmooth dynamics, each one being either linear or non-linear and simulated with event-capturing or event-detecting schemes. Available time integration schemes are Moreau-Jean scheme, Schatzman-Paoli scheme, nonsmooth Newmark scheme [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF][START_REF] Chen | A nonsmooth generalized-α scheme for flexible multibody systems with unilateral constraints[END_REF], Odepack suite and DAE solvers developed by E. Hairer [START_REF] Hairer | Solving Ordinary Differential Equations I[END_REF].

Inside SICONOS/KERNEL, the description of the user model relies internally on a primary graph with dynamical systems as nodes and interactions as edges. During time evolution, the determination of active Inria constraints corresponds to the building of sub-graphs of the primary graph. In the case of problems formulated in local coordinates an adjoin graph transformation is dynamically applied to sub graphs in order to drive the computation of the components of the optimization problem given to Numerics solver.

The primary graph structure may be entirely specified once by the user, or it may also be rebuilt during the simulation, as, for example, the result of a contact detection broad-phase. Links with pre/post processing software for collision detection (Bullet contact detection [START_REF]Bullet Physics Library[END_REF]) and CAD libraries (OpenCascasde [START_REF]CASCADE Technology[END_REF] and PythonOCC [START_REF]PythonOCC. 3D CAD/CAE/PLM development framework for the Python programming language[END_REF]) has been set-up. The latter link enables the contact detection between bodies defined by B-Rep (splines, nurbs, ...) and then to compute gaps and local frame at contact. Such a geometrical representation of the data allows to use realistic modeling of surfaces with large sliding. The possible impact laws for mechanics are complementary condition, impact with or without friction.

SICONOS is written in C++ and has an object oriented architecture, there is no graphical user interface (GUI). A simulation is run through a C++ or Python script in which are gathered the objects creations and the method calls on these object that make up the computation. Furthermore, not only the coarse level of command driving is available, but the time loop can be exploded in every step to allow for an interactive simulation. In the same way the whole database is accessible through copy or reference, allowing for a efficient access and aimed at designing any post-processing functions.

More details can be found at http://siconos.gforge.inria.fr.

Motivation for using the Nonsmooth Mechanics approach

Here are listed some of the main interests of using an approach based on nonsmooth Mechanics when we want to deal with multi-body systems with unilateral contacts, friction and impacts:

1. The use of complementarity formulation allows us to deal with hard contacts in a very efficient way avoiding the issues related to the time-integration of stiff problems. When we want to simulate contacts on hard surfaces, avoiding deep inter-penetrations obliges to use large stiffness parameters in a standard approach. This is overcame in the presented approach. 2. Nonsmooth contact laws are based on few parameters (coefficient of restitution and friction). Contact parameters are often very difficult to measure and have a great variability. Proposing contact laws with few parameters that represents the main effects at contact (inter-penetrability avoidance, friction with real sticking phase and plastic dissipation) provides the user with robust simulation tools in view of Control applications. Furthermore, nonsmooth contact laws can be easily extended to take into account more enhanced effects such as stiffness, damping or cohesion [START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF]. For soft-soil simulation, piecewise continuous Bekker model could be implemented in this context.

3. The associated simulation methods are proven to be efficient on a large number of applications from medium to large scale multi-body systems [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF]. Granular materials with large number of contacts points (between 10 4 and 10 6 ) can be simulated with such algorithms. 4. Such an approach has been proved to be robust and efficient on several important industrial applications. Siconos is for instance used for more than 5 years in a strong collaboration between Schneider Electric and INRIA for the virtual prototyping of the mechanisms of circuit breakers [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and coulomb's friction[END_REF]. The robustness with respect to contact parameters and the efficiency of the simulation that avoids problems due to stiffness provides Schneider Electric with a robust predictive tool for new designs. This numerical methods are also the object of a closed collaboration with ANSYS for the rigid body dynamics toolbox [START_REF] Haddouni | Comparison of index-3, index-2 and index-1 dae solvers for nonsmooth multibody systems with unilateral or bilateral constraints[END_REF].

Rover simulation in Siconos

The Rover simulation in the Siconos platform allows us to simulate the robot motion taking into account rigorously dynamics and non-regular models of contact friction. Roughly, coding a Siconos application consists in two main steps : the first is to describe the dynamic Rover model and the contact model interaction between wheels and the environment; the second is to define appropriate numerical parameters for the model. More precisly, any Siconos code (C++ or Python) is composed of the following sub-steps:

1. define some DynamicalSystem instances RT n°448

2. define some Interaction instances, and for each Interaction: define Relation define NonSmoothLaw 3. build a NonSmoothDynamicalSystem : all DynamicalSystem + all Interaction objects arranged into a directed Graph 4. build a Model that holds the NonSmoothDynamicalSystem. 5. define a Simulation: the way the behavior of the NonSmoothDynamicalSystem will be computed: define a TimeDiscretisation choose a strategy: TimeStepping or EventDriven define some integrators for the DynamicalSystems (OneStepIntegrator) 6. define a way to formalize and solve the possible non smooth problems (formulation + solver, the OneStepNSProblem)

7. associate the Simulation to the Model 8. launch the Simulation (time loop)

Modeling of the system

The Rover dynamical model is a LagrangianDS siconos object (lagrangian formulation). For 3D Rover modeling, the inertia matrix and contact model interaction relations are difficult to set-up by hands. We use a Maple™ code to automatically generate these relations, and to export into C files. In Siconos, we import these relations as a plug-in function that are dynamically loaded. The HuMAns toolbox( [START_REF] Wieber | The humans toolbox, a homogenous framework for motion capture, analysis and simulation[END_REF]) using Maple™ software tool enable us to easily do complex modeling easily with high accuracy.

Originally developed for humanoid robotics research, the HuMAnS toolbox (for Humanoid Motion Analysis and Simulation) proposes a wealth of state -of -the -art algorithms from the field of robotics research for the modeling, the analysis and the simulation. In our case, we used it to generate the multibody dynamic model for the Rover. For this, the user need to describe the Rover parameters through several Maple™ input files. One file must describe kinematic parameters of the Rover, it indicates the relative location of key-points of the model and the degree of freedom of every connecting point. An other file describes the dynamical parameters, it describes the inertials parameters: the gravity vector, the mass of the segment, the position of the segments centers of mass and the inertia matrix of these segments relative to the centers of their attached frames.

Then, the Maple™ generator included in HuMAns produces C code for M (q) the system inertia, G(q) the forces which derived from a potential energy (for example the gravity), N (q, q) a non-linear effects matrix.

For the simulation we used a simplified model inspired from the ExoMars rover chassis with six wheels on three simple passive bogies located forward at each side of the rover and one on the rear (see [START_REF] Michaud | Development of the exomars chassis and locomotion subsystem[END_REF]). The payload is simply modelised with a mass uniform parallelepiped. The visual representation and joints notation are given on the figure 2. The rover has 21 Degree Of Freeedom (DOF) : 6 for the location of the rover, 6 for the rotation of each wheel, 6 for each wheel steering and 3 for each passive joint of the bogies.

To describe the contact of the Rover with the solid ground or granular soil in Siconos, we need two plugin functions: h to compute the distance between contact points and G which is the transformation matrix of velocity of contact points from the global frame to the local frame.

To simplify the model, we only consider point contact. In this contact model, we simplify the axis of the wheel into a point, and create contact model between this point and the sphere, a plane or a mesh.

For the simulation results given below, the ground is modeled by meshes. In order to be able to perform a simulation of the rover roaming on a triangular mesh in Siconos one needs to have the appropriate interaction and relation classes defined. Relation class is an interface providing distances between the surfaces being checked for collision as well as jacobians, which in turn provide mappings from the local to global coordinate frames. While distances are needed for collision detection, jacobians are crucial for resolving impacts and nonsmooth friction problems. Distances calculation is a broad topic in itself, briefly speaking these can be delivered by a third-party collision detection library or calculated according to an explicit algorithm built-in in the respective relation class. In the standalone simulation the latter approach has been chosen. Namely, distance function has been defined explicitly in order to calculate the gap between a contact point on each wheel of the robot and a triangle arbitrarily oriented in space. Everytime this function returns zero (with certain tolerance) a collision is declared. Jacobian is defined within the relation class as an operator resulting from the multiplication of the local triangle's rotation matrix with an inverse of robot's geometric jacobian. Relation object together with a non-smooth law define the interaction between systems. In the standalone simulation an interaction is instantiated between each wheel and each triangle. Thus, there is six times as many interactions as wheels. This approach can become inefficient as the number of triangles (and complexity of the mesh) grow. To compensate for this problem a space filter should be used in the future which in essence filters out and discards interactions between objects which are far enough.

The C code generated, for the rover model and the contact model interaction, is then imported in Siconos code plugin function. It enables us to easily do complex modeling with high accuracy. This work of modeling is detailed in the report [START_REF] Acary | Simulation of a mars rover on granular soils[END_REF].

2 Mechanical tests of the rover

Parameters used in simulations

For below presented simulations following parameters have been manipulated:

• Vector of initial positions and velocities of the joints and the base q 0 and v 0 • Friction and normal restitution coefficient µ and eps n

• Simulation time and timestep T and h

• Torques applied to the joints τ For all simulations gravity coefficient on Mars have been used (3.7 m s 2 ).

Tests scenarios specification

Seven tests have been defined and performed:

1. Rover stabilization on a horizontal plane after a free-fall phase.

In the first test setting the rover falls freely from the height of 2 meters onto a horizontal plane. No torques are applied in any of the joints. The only external forces acting on the rover are the gravity and the ground reaction forces. Initial position of the rover has been set to (x, y, z) = (5, 5, 2) [m].

Friction coefficient has been set to 0.7. This case has been divided into two sub-cases: first one with the normal restituion coefficient set to zero and the second one with the normal restitution coefficient set to 0.2. Tangential restitution coefficients have been set to zero in both sub-cases.

2. Rover stabilization on a horizontal plane with linearly varying velocity of one steering axis.

In the second test setting rover is dropped onto a horizontal plane and stands idle on it during the first 50 seconds of the simulation. After 50 seconds from the beginning of the simulation a constant torque τ = 0.00002N m is applied in one of the steering axes (steering axis FL) causing its rotational motion with linear velocity. Wheel makes two full rotations around its vertical steering axis (FL).

Other external forces acting on the rover are the gravity and ground reaction. Initial position of the center of mass of the robot has been set to (x, y, z) = (5, 5, 2) [m]. Friction coefficient has been set to 0.4. Restitution coefficients (tangential and normal) have been set to zero.

3. Rover stabilization on a horizontal plane with acceleration and deceleration phases.

In the third test setting rover's parcour is divided into six time intervals with different values of torques applied to the wheels:

• 0s < t s < 20s, τ = 0N m RT n°448 • 20s < t s < 40s, τ = 0.0007N m • 40s < t s < 80s, τ = 0N m • 80s < t s < 105s, τ = -0.0007N m • 105s < t s < 150s, τ = 0N m • 150s < t s < 200s, τ = -0.0005N m
Other external forces acting on the rover are the gravity and ground reactions. Initial position of the center of mass of the robot has been set to (x, y, z) = (5, 8, 0.57) [m]. In the 5 th phase rover effectively stops its motion until negative torques are applied. Friction coefficient has been set to 0.7.

Restitution coefficients (tangential and normal) have been set to zero.

4. Rover stabilization on an inclined plane with a phase of upward motion.

In the fourth test setting rover is set to stand steadily on an inclined plane. Inclination angle of the slope has been set to 10 • . Torques applied to the wheels counterbalance torques caused by the gravity. After initial period higher torques are applied to the wheels which cause the robot to drive upwards. Blocking torques are equal to τ b = -0.87072N m. During the second phase of motion linearly varying torques are applied to all wheels. Rover's parcour has been devided into two phases:

• 0s < t s < 30s, τ bconstant blocking torques

• 30s < t s < 40s, τ mlinearly increasing torques

Friction coefficient has been set to 0.8. Restitution coefficients (tangential and normal) have been set to zero.

Rover stabilization on an inclined plane with acceleration and deceleration phases.

In the fifth test setting rover is dropped onto the inclined plane. Inclination angle of the slope has been set to 10 • . After the 10 th second of motion a simple PD controller is set using position and velocity of the center of mass in order to stop the rover on the slope. Rover stops effectively at t = 25s. Friction coefficient has been set to 0.8. Restitution coefficients (tangential and normal) have been set to zero.

6. Spherical obstacle crossing on a horizontal plane.

In the sixth test setting rover is dropped onto the horizontal plane from the height of 2 m. After 10 s linearly varying torques are applied to all wheels. A spherical obstacle has been set in front of the rover. The obstacle is in the form of a sphere which protudes outside the plane. Center of the sphere has been set so that the protrusion is equal to 20 cm. Coefficients of friction and restitution (normal) are in this case respectively: 0.3 and 0.

7. Vertical step obstacle crossing.

In the seveth test setting rover is dropped onto a horizontal plane. An obstacle in the form of a horizontal step of height 0.1m (roughly equal to the wheel radius) has been set in front of the rover. At certain point of time (t = 10s) rover starts moving towards the step and crosses it mounting on the higher plane. Friction coefficient has been set to 0.6. Restitution coefficient (normal) has been set to 0.1.

Results of each of the above tests will be detailed in the following sections.

In mechanical simulations with nonsmooth approach one is mostly interested in the following quantities:

• x COM -mass center coordinates The following essential quantities have been plotted in this case:

Comments

As seen in the figure 4, coordinate z of the center of mass of the rover decreases parabolically until the system is stable on the plane. This represents the free-fall phase under gravity. After reaching the plane there is no rebound as the coefficient of restitution in the nonsmooth contact law has been set to zero. Impact corresponds to a sharp peak in the force as seen in the third sub-figure .  Also, the following complimentary quantities have been plotted:

Comments

Plots in the figure 5 correspond to forces (normal and tangential impulsions) computed in the local contact frame. Contact forces are always computed in the local frame and then expressed in the coordinates of the system using jacobian matrix of the constraint function. Constraint function (gap function) used to detect the contact and the relative velocity are plotted in the figures 6. Gap function is a relative distance between two contacting bodies and it imposes geometric constraints on the system. Once slightly below zero it will trigger a nonsmooth interaction between those bodies.

Normal restitution coefficient set to eps

n = 0.2

Comments

As the restitution coefficient has been set to a non-zero value in this case one observes multiple impacts. Impacts corresponds to peaks of forces in the third sub-figure and nonsmooth changes of velocity in the second sub-figure. As the energy is dissipated the magnitudes of impacts reduce in time.

Also, the following complimentary quantities have been plotted: In this case, following essential quantities have been plotted:

Comments

In the figure 11 one can see the variation of the steering axis displacement and velocity. After a constant torque is applied the axis is rotating with a linear velocity which corresponds to the quadratic shape of the displacement curve.

Also, the following complimentary quantities have been plotted:

Comments

In the figure 12 one can observe one impulsion (force) peak corresponding to the impact when the rover falls onto the ground. In the figure 13 the gap function and relative velocity are displayed.

Test 3 -Rover stabilization on a horizontal plane with acceleration and deceleration phases

In the third test setting rover's parcour is divided into six time intervals with different values of torques applied to the wheels:

1. 0s < t s < 20s, τ = 0N m 2. 20s < t s < 40s, τ = 0.0007N m 3. 40s < t s < 80s, τ = 0N m 4. 80s < t s < 105s, τ = -0.0007N m 5. 105s < t s < 150s, τ = 0N m 6. 150s < t s < 200s, τ = -0.0005N m
Other external forces acting on the rover are the gravity and ground reactions. Initial position of the center of mass of the robot has been set to (x, y, z) = (5, 8, 0.57) [m]. In the 5 th phase rover effectively stops its motion until negative torques are applied. Friction coefficient has been set to 0.7. Restitution coefficients (tangential and normal) have been set to zero.

Comments

In the figure 16 one can see the step pattern of the reaction torques resulting from the fact that the rover is moving along the x axis. Jumps in value represent changes in the friction force along x axis as torques are applied to the wheels. In the figure 17 all components of the local contact forces are displayed.

Step pattern along x axis can be observed.

6 Test 4 -Rover stabilization on an inclined plane with a phase of upward motion

In the fourth test setting rover is set to stand steadily on an inclined plane. Inclination angle of the slope has been set to 10 • . Torques applied to the wheels counterbalance torques caused by the gravity. After initial period higher torques are applied to the wheels which cause the robot to drive upwards. Blocking torques are equal to τ b = -0.87072N m. During the second phase of motion linearly varying torques are applied to all wheels. In this case, following essential quantities have been plotted:

Figure 19: x COM -mass center coordinates

Comments

Two phases are distinguished in this case. In the figure 19 one can see that in the first phase rover is motionless because of the constant blocking torques applied to all wheels. These torques through friction nullify the effect of gravity on the system. Coordinates of the center of mass are constant until the second phase. This means that the friction-contact model is efficient in the sense that it matches closely physical behavior of the phenomenon. One can see this more clearly in the velocity curve where the velocity of the center of mass is zero in the first phase of motion. Nonsmooth friction-contact model allows for this effect as it does not rely on regularization of friction laws. Last curve in the figure 19 represents the ground reaction forces applied to the mass center of the system. This force increases very slightly as the rover starts moving upwards because of the friction force.

Following additional quantities have been plotted:

Comments

In the figure 20 one can see normal and tangential components of the local contact forces. Immediate observation is that the projections of forces oscillate during the first phase of motion. This is due to the fact that contact solver finds a different solution to the friction-contact problem at each time instant.

7 Test 5 -Rover stabilization on an inclined plane with acceleration and deceleration phases

In the fifth test setting the rover is posed on and inclined plane. Inclination angle of the slope has been set to 10 • . After the 10 th second of simulation a simple PD controller is set using position and velocity of the center of mass, in order to stop the rover on the slope. Rover stops effectively at t = 25s. Friction coefficient has been set to 0.8. Restitution coefficients (tangential and normal) have been set to zero. In this setting we are mostly interested in two phases:

1. 0s < t s < 10s 2. 10s < t s < 40s

In this case, following essential quantities have been plotted:

Comments

One can clearly see two phases of motion of the rover. In the end of the second phase rover stops its motion on the inclined plane which once again demonstrates the efficiency of the nonsmooth friction-contact law. This can be seen in the figure 23 where in the second phase velocity components converge towards zero. This effect corresponding to a true nature of friction and can be captured in numerical experiments using nonsmooth laws. Using regularized laws to describe friction would not allow the system to retain non-zero torques with zero velocity.

Comments

As in the previous test setting, an immediate observation comes into mind as to the oscillations in the figure 24. Indeed, it is the same reason as previously that they appear in the plot.

Test 6 -Spherical obstacle crossing on a horizontal plane

In the sixth test setting rover is dropped onto the horizontal plane from the height of 2 m. After 10 s linearly varying torques are applied to all wheels. A spherical obstacle has been set in front of the rover. The obstacle is in the form of a sphere which protudes outside the plane. Center of the sphere has been set so that the protrusion is equal to 20 cm. Coefficient of friction has been set to 0.3 whereas the coefficient of restitution to 0.

Comments

In this setting attention can be drawn to the spherical obstacle crossing. In the figure 27 in the first sub-plot on can see z coordinate evolution of the mass center of the robot. It is clearly seen how around 11 th minute rover mounts the sphere. In the sub-plot below one can see the evolution of velocities of the mass center where several small jumps occur in the y and z dimension as a result of impacts. Velocity in the x direction is clearly larger as the rover moves along the x axis. Impacts are more clearly seen in the third sub-plot with each peak corresponding to an impact.

Following additional quantities have also been plotted:

•

  x wheels -wheels angular displacement • v COM -mass center velocity Inria • v wheels -wheels angular velocity • R -reaction forces (impulsions) in lagrangian (global) coordinates • λ N (λ n) -normal component of the contact force (impulsion) in local coordinates • λ Tx (λt) -tangential component of the contact force (impulsion) in local coordinates in the x direction • λ Tz (λ s) -tangential component of the contact force (impulsion) in local coordinates in the z direction • y N (y n) -gap function (distance between contact point and the constraint function) • ẏN ( ẏn ) -normal component of the local contact velocity • ẏTx ( ẏt) -tangential component x of the local contact velocity • ẏTz ( ẏs ) -tangential component z of the local contact velocityFor each scenario a subset of the above quantities has been plotted.

Figure 1 :Figure 2 :

 12 Figure 1: Assumed model of the robot

3. 1

 1 Normal restitution coefficient set to eps n = 0.0 In the first test setting rover falls freely from the height of 2 meters onto a horizontal plane. No torques are applied in any of the joints. The only external forces acting on the rover are the gravity and the ground reaction forces. Initial position of the rover has been set to (x, y, z) = (5, 5, 2) [m]. Friction coefficient has been set to 0.7. In this sub-case restitution coefficients (tangential and normal) have been set to zero.

Figure 3 :

 3 Figure 3: First test setting (arrows depict contact forces)

Figure 4 :

 4 Figure 4: Position, velocity and reaction forces of the center of mass

Figure 5 :

 5 Figure 5: Normal and tangential components of the local contact force (impulsion λ) for each wheel during impact

Figure 7 :

 7 Figure 7: Position, velocity and reaction forces of the center of mass

Figure 8 :

 8 Figure 8: Normal and tangential components of the local contact force (impulsion λ) for each wheel during impact

Figure 10 :

 10 Figure 10: Second test scenario (arrows depict contact forces)

Figure 11 :

 11 Figure 11: x F L , v F L -angular displacement and velocity of the FL axis

Figure 12 :

 12 Figure 12: λ N -normal component of the contact force (impulsion) for each wheel

Figure 14 :Figure 15 :

 1415 Figure 14: Third test scenario (arrows depict contact forces)

Figure 18 :

 18 Figure 18: Fourth test scenario

Figure 20 :

 20 Figure 20: λ N , λ Tx , λ Tz -normal and tangential components of the contact force (impulsion) for each wheel

Figure 22 :

 22 Figure 22: Fifth test scenario

Figure 23 :

 23 Figure 23: position, velocity and reaction force of the mass center

Figure 26 :

 26 Figure 26: Sixth test scenario

Figure 28 :

 28 Figure 28: λ N , λ Tx , λ Ty -normal component of the contact force (impulsion) for each wheel (interaction with the plane)

Figure 30 :

 30 Figure 30: Seventh test scenario (white arrows depict the reaction contact forces)

1 . 3 1 . 4 1 . 5 6 2 2 . 1 2 . 2 7 3 14 4

 13141562122714 Figure 32: Rover climbing the step

Michalczyk & Brémond & Acary & Pissard-GibolletSimplified kinematic structure of the robot is as follows:

RT n°448

Michalczyk & Brémond & Acary & Pissard-Gibollet

Inria

Comments

In the figure 9 one can see the gap function and relative velocity plots. With respect to the previous sub-case rebounds off the ground are visible.

Following additional quantities have also been plotted:


Julio Toss
email: jtoss@inf.ufrgs.br

João Comba
email: comba@inf.ufrgs.br

Bruno Raffin
email: bruno.raffin@inria.fr

Parallel Shortest Path Algorithm for Voronoi Diagrams with Generalized Distance Functions

Keywords: Voronoi Diagram, GPGPU, Parallel Programming, Physics Based Simulation

Stiffness gradient from soft (blue) to stiff (red) (b) Compliance distance: far (blue) to close (red) (c) Voronoi diagram Fig. 1. Computing the Voronoi diagram with generalized distance functions on a voxel grid with dimensions 100x40x20 for 20 seeds: The distance function is defined on a non-Euclidean space which depends on specific material properties, in this example a gradient that encodes material stiffness (a). The resulting distance map (b) is obtained after computing the shortest path on the voxel grid from each Voronoi seed (c). The GPU grid-based parallel shortest path algorithm we propose was 27 times faster than the reference sequential CPU implementation.

Abstract-Voronoi diagrams are fundamental data structures in computational geometry with applications on different areas. Recent soft object simulation algorithms for real time physics engines require the computation of Voronoi diagrams over 3D images with non-Euclidean distances. In this case, the computation must be performed over a graph, where the edges encode the required distance information. But excessive computation time of Voronoi diagrams prevent more sophisticated deformations that require interactive topological changes, such as cutting or stitching used in virtual surgery simulations. The major bottleneck in the Voronoi computation in this case is a shortestpath algorithm that must be computed multiple times during the deformation.

In this paper, we tackle this problem by proposing a GPU algorithm of the shortest-path algorithm from multiple sources using generalized distance functions. Our algorithm was designed to leverage the grid-based nature of the underlying graph used in the simulation. Experimental results report speed-ups up to 65x over a current reference sequential method.

I. INTRODUCTION

The Voronoi diagram is a classical subdivision of space that is suitable for answering proximity problems, such as finding the nearest site, facility location, motion planning, coverage in sensor networks, etc. There are many variations of Voronoi diagrams, which are often associated to the underlying distance measure used. In this work we are particularly interested in the computation of Voronoi diagrams for non-Euclidean spaces, where a generalized distance function is used. This motivation came from the area of physics-based simulation, where Voronoi diagrams were used by Faure et at. [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF] to simulate meshless deformable objects with heterogeneous material properties and complex geometries. Their proposal relies on a novel method which uses material-aware shape functions to describe the composition of simulated bodies, which can be composed of both soft and stiff materials.

The computation of Voronoi diagrams plays a central role in their proposal. To have accurate and realistic deformations, the underlying deformation space is discretized into a grid, and the material stiffness is defined for each vertex of the grid (material map). The deformation algorithm considers grid vertices as simulation nodes, with an associated Voronoi kernel function that limits the region of influence of the node. Since the distance function is not computed in a standard Euclidean plane, it must be scaled according to compliance values in the material map. Therefore, points connected by similar materials (i.e. inside the same Voronoi cell) will deform in a similar way.

The computation of the Voronoi kernel function and hence the Voronoi diagram is done during the setup phase of the simulation and remains unchanged during the whole simulation, as long as the topology of the object and material properties do not change. This allows the subsequent simulation phase to be performed in real-time, a necessary requirement for interactive applications. However, to enable interactive changes in the topology of the object (e.g. cutting or stitching), the Voronoi diagram must be recomputed during the simulation. This is not possible in their solution, since they report initialization times ranging from less that 1 second for a grid of 100x40 voxels to 10 minutes for a 500x200 grid. However, their implementation is strictly sequential leaving plenty of room for optimization and parallelization.

Contributions: In this work we explore certain properties of the deformation problem to speed up the computation of shortest paths in graphs, which are used in the Voronoi diagram computation. Our proposal describes a shortest path algorithm using a parallel implementation that leverages the processing power of Graphics Processing Units (GPUs). Current GPU proposals for shortest path algorithms consider the Single-Source-Shortest-Path (SSSP) problem [START_REF] Harish | Large Graph Algorithms for Massively Multithreaded Architectures[END_REF], using classic algorithms such as Dijkstra [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF] or Bellmann-Ford. Unlike these algorithms, our algorithm considers the Multiple-Source-Shortest-Path (MSSP), since multiple shortest path computations are triggered at each simulation node. Moreover, we leverage the fact that our shortest-path algorithm can be computed over a grid instead of general graphs (used in current algorithms), which allows the algorithm to be more efficient.

We implemented our algorithm in CUDA and tested different design decisions with a collection of simulation examples. We prepared an experimental evaluation comparing our GPU implementation against the sequential CPU reference method, and obtained speed-ups ranging form 3x for small inputs up to 65x for larger ones, both with synthetic and real datasets.

II. RELATED WORK

The Voronoi diagram is a data structure extensively studied in the context of computational geometry for many different applications. Originally it defines a region of proximity for a set of k sites (seeds) in a plane where the distance between points is defined by their Euclidean distance. Most works use methods to efficiently compute them on a contiguous Euclidean space [START_REF] Rong | GPU-assisted computation of centroidal Voronoi tessellation[END_REF], or to compute their discrete approximation [START_REF] Hoff | Fast computation of generalized Voronoi diagrams using graphics hardware[END_REF], [START_REF] Rong | Jump flooding in GPU with applications to Voronoi diagram and distance transform[END_REF].

Although most works in the literature compute Voronoi diagrams on Euclidean spaces, there are generalizations in the context of graphs [START_REF] Erwig | The graph Voronoi diagram with applications[END_REF], [START_REF] Hurtado | The weighted farthest color Voronoi diagram on trees and graphs[END_REF], sometimes called the Graph Voronoi Diagram [START_REF] Erwig | The graph Voronoi diagram with applications[END_REF]. In this context, the distance metric considered corresponds to the shortest path between nodes. This formulation of Voronoi diagrams often arises in the field of Facility Location, where clients and suppliers lie in an interconnection network. Computing these Voronoi diagrams basically consists in concurrently computing shortest paths from multiple sources on a weighted graph.

A. Parallel Voronoi diagrams computation

Over the last decade, the advent of parallel processors motivated the scientific community to put its effort in creating parallel algorithms for classical problems. At the same time, GPUs evolved to become general purpose massive parallel processors attracting attention from other fields of computing. Early study made by [START_REF] Hoff | Fast computation of generalized Voronoi diagrams using graphics hardware[END_REF] used graphics hardware (pre-CUDA) to compute an approximation of the Voronoi diagram. Most works found in the literature rely on a discrete approximation, usually either on 2D pixel-maps or on surfaces.

The many parallel approaches vary in the way the information of proximity from the Voronoi centers is propagated to each pixel. Jump Flooding Algorithm (JFA) is proposed by [START_REF] Rong | GPU-assisted computation of centroidal Voronoi tessellation[END_REF], [START_REF] Rong | Jump flooding in GPU with applications to Voronoi diagram and distance transform[END_REF] as an algorithmic paradigm for GPGPU with application on Voronoi diagram computation. In JFA, the seeds start propagating their coordinates to neighbor pixels according to a pattern that halves the offset at each step. Each pixel compares the new information received with the current one and keeps the coordinates of the closest seed. In this case, distances can easily be computed on the Euclidean plane.

Weber et al. [START_REF] Weber | Parallel algorithms for approximation of distance maps on parametric surfaces[END_REF] introduces an interesting parallel algorithm called parallel marching method (PMM) to compute distances on surfaces with application on Voronoi diagrams. This method is indeed an extension of the fast marching method which is based on a priority-queue. However, this kind of data structure is difficult to efficiently parallelize. Instead, their method uses a specific traversing order of the grid, called raster scan, which shows an efficient parallelization algorithm.

More recently, [START_REF] Reem | On the possibility of simple parallel computing of Voronoi diagrams and Delaunay graphs[END_REF] proposed a substantially different method from the previous ones. It uses a combinatorial approach to compute, in parallel, the exact polygons that form each cell of the Voronoi diagram. This work however is mainly theoretical, showing formal proofs without bringing experimental results.

All of these works consider distance computation on the Euclidean space, none of them deals with Voronoi diagrams on the graph space. This means that the shortest distances are always straight lines, hence these methods cannot be directly applied on graph problems. To the best of our knowledge there are no algorithm addressing the parallelization of the Graph Voronoi Diagram.

B. Parallel graph algorithms

Distance computation for the construction of the Voronoi diagram defined over graphs require finding shortest paths (with multiple sources plus concurrent search) [START_REF] Erwig | The graph Voronoi diagram with applications[END_REF]. The original Dijkstra's sequential algorithm for solving the Single-Source Shortest-Path (SSSP) [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF] has O(V 2 ) complexity on the number of vertices, while the min-priority-queue based version has complexity O(E + V logV ) , where E is the number of edges.

Several parallel approaches to solve the SSSP problem have been proposed on the literature. Crauser et Al. [START_REF] Crauser | A parallelization of Dijkstra's shortest path algorithm[END_REF] proposed a parallel PRAM algorithm of Dijkstra's which needs O(n 1/3 log n) time. However, Dijkstra is an inherently sequential algorithm, with lots of synchronizations with no efficient PRAM implementation [START_REF] Harish | Large Graph Algorithms for Massively Multithreaded Architectures[END_REF]. Alternatively, other works parallelize the SSSP based on Bellman-Ford's algorithm which is less efficient than Dijkstra on sequential implementations, but has a higher degree of parallelism [START_REF] Kumar | A modified parallel approach to Single Source Shortest Path Problem for massively dense graphs using CUDA[END_REF], [START_REF] Nasre | Atomic-free irregular computations on GPUs[END_REF].

Most of the existing parallel SSSP algorithms have to deal with a trade-off between the amount of parallelism exposed and the extra work generated. The parallel deltastepping method, proposed on [START_REF] Madduri | Parallel shortest path algorithms for solving large-scale instances[END_REF], has a good compromise between these two factors. They report an implementation exhibiting 30x speed-up on a CRAY MTA-2 shared memory architecture with 40 processors.

1: procedure RELAX(u, v, w) 2: if v.d > u.d + w(u, v) then 3: v.d ← u.d + w(u, v) 4:
end if 5: end procedure

C. GPU implementation of SSSP algorithms

Several GPU implementations have been proposed over the last years for different graph algorithms [START_REF] Harish | Large Graph Algorithms for Massively Multithreaded Architectures[END_REF]. For the shortest path problem specifically, Dijkstra-based parallelizations are more frequently used [START_REF] Harish | Large Graph Algorithms for Massively Multithreaded Architectures[END_REF], [START_REF] Madduri | Parallel shortest path algorithms for solving large-scale instances[END_REF], [START_REF] Ortega-Arranz | A new GPU-based approach to the shortest path problem[END_REF]. Although, other approaches exist (e.g [START_REF] Kumar | A modified parallel approach to Single Source Shortest Path Problem for massively dense graphs using CUDA[END_REF]), which proposes a modified Bellman-Ford algorithm on GPU for dense graphs.

In general, Dijkstra based algorithms use a technique known as edge relaxation. In this technique each vertex maintains an estimate shortest-path with distance v.d. The process of relaxation consist of trying to improve this estimate by going from vertex u to v through an edge of weight w(u, v) (Fig. 2).

When done in parallel, each vertex u is assigned to a thread which may update v.d concurrently, thus creating a critical section. Consequently, lines 2 -4, of algorithm in Fig. 2, have to be protected in an atomic region. In modern CUDA devices this atomic region can be efficiently implemented by the single atomic instruction atomic min(addr, val) 1 .

III. PARALLEL GRAPH VORONOI ALGORITHM

As mentioned on the previous section, the Graph Voronoi can be seen as an extension of the shortest path problem. However its parallelization poses additional problems of concurrent access on shared variables. In the Voronoi diagram problem, each voxel has to keep the distance estimate value to the seed and an extra variable for its Voronoi cell index. These variables would then be updated serially in the relaxation procedure, which, if executed by two threads in parallel, could lead to any combination of results in these variables. On concurrent programming this is a classical case of race condition. The straightforward solution for this problem would be to enclose the whole critical section (Fig. 2) within mutex locks. However, mutexes are expensive structures to implement on GPUs. To deal with this problem, we choose to encode both variables, Voronoi index and distance estimate, in a single 32-bit word (Fig. 3) which can then be atomically updated in a single atomic min() instruction. Our encoding can be adjusted to balance distance precision and maximum number of Voronoi cells. In our implementation, we reserved 24 bits for the distance and 8 bits for the Voronoi region index. 1 atomic min(addr, val) reads word old located at the address addr, computes the minimum of old and val, and stores the result back to memory at the same address. These three operations are performed in one atomic transaction [START_REF]Cuda c programing guide[END_REF].

A. Data structure

Our data representation in memory substantially differs from the classical graph data structures. Instead of using adjacency matrices or lists for the shortest path computation, like in [START_REF] Harish | Large Graph Algorithms for Massively Multithreaded Architectures[END_REF], we are dealing directly with images which are 3D node matrices. Each node keeps its compliance value, Voronoi cell index, and distance to a Voronoi source. The connectivity between nodes is given by their natural neighbors in the 3D volume (26 neighbors). The weight of each edge is given by some generalized distance function, dist(C v , C u ), relating adjacent neighboring nodes. In this work specifically, we employ the compliance scaled distance function used by Faure et Al. [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF]. In this case, the distance between two adjacent nodes is a function of the measure of compliance of the material at each node.

B. Base algorithm

Our algorithm uses four internal arrays, C 0 , C 1 , Vor and Mask , stored on the GPU global memory and with same size of the input volume (Fig. 5). The cost arrays C 0 and C 1 are used to keep the shortest-path estimates of each voxel. They are initialized with 0 at the seeds and ∞ (maximum unsigned integer value) everywhere else. The Voronoi diagram, stored on array Vor , is initially empty on every voxel, except for those corresponding to the seed's coordinate which are initialized with a unique Voronoi cell index. Finally, the boolean array Mask is used as activity mask to mark which voxels have an updated cost estimate indicating that it will be relaxed on the next step.

We assign one thread to every voxel . The execution then follows a scatter approach (Fig. 4) where each active thread, marked on Mask , will relax the cost estimate of its neighbors and set their correct Voronoi index.

The algorithm is divided in two parallel phases: relaxation for all v ∈ Mat do 3:

C 0 [v] ← ∞; C 1 [v] ← ∞ 4:
end for

5:

for all s ∈ Seeds do 6:

C 0 [s] ← 0; C 1 [s] ← 0 7: Mask [s] ← true 8:
Vor [s] ← idx + + until TERM 15: end procedure tid ← getT hreadIndex()

3: if Mask [tid] then 4:
for all neighbors nid of tid do 5: and update. The host code (Fig. 5) initializes the data structures and then iteratively calls the GPU kernels RELAXKERNEL (Fig. 6), add UPDATEKERNEL (Fig. 7), until the termination condition is satisfied. The distance function, at line 5 in RELAXKERNEL, computes the local distance between two neighbor voxels based on their compliance values in the material map [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF]. At each iteration, C 1 maintains the intermediate values computed during the relaxation. In the UPDATEKER-NEL procedure, the values from C 1 are copied back to C 0 and the activity mask is updated. The duplication of these cost matrices is needed to avoid read-after-write hazards when writing to global memory. The algorithm finishes when the diagram reaches a fixed point, where no more voxels are updated.

d new ← C 0 [tid] + localDist(tid, nid, M at

C. Algorithm with stream compression

As an enhancement to our base algorithm, we tried to reduce the number of idle threads by applying stream compression [START_REF] Harish | Large Graph Algorithms for Massively Multithreaded Architectures[END_REF], [START_REF] Hoberock | Stream compaction for deferred shading[END_REF]. Fig. 8 shows how the computation propagates to neighbors in a form similar to a wave. The black front indicates, at each step, which threads have True in the activity mask at the beginning of the relaxation kernel (at line 3 of algorithm in Fig. 6). Over the execution of the algorithm, we note that the number of active threads is much lower than the 1: procedure UPDATEKERNEL 2:

tid ← getT hreadIndex() grid size and also varies considerably along time (Fig. 12). This causes our thread blocks to be very inefficient as most of the threads will actually evaluate the conditional to False, without computing anything (line 3, Fig. 6) .

3: if C 0 [tid] > C 1 [tid] then 4: C 0 [tid] ← C 1 [tid]
Compact blocks solve this problem by grouping all the active thread in fewer blocks, thus reducing branch divergence, as well as the runtime overhead of scheduling idle threads.

Implementation: Stream compression is performed by a Scan operation over the activity mask followed by a Scatter. These operations can be easily implemented in CUDA using the Thrust template library [START_REF] Bell | Thrust: A productivity-oriented library for cuda[END_REF]. The result of the compression is an array mapping thread indexes to pixel coordinates, that are used in Kernel1 to retrieve the correct data.

This process adds a non-negligible overhead which sometimes can actually supersede the gains of performance. To be able to balance the trade-off between performance gains of compression and time spent by the scan+scatter process, we implemented a variable grain compression.

Variable grain compression: This method defines a coarser subdivision over the activity mask as show in Fig. 10. The coarser mask is parametrized by its grain size, which is set by their x, y and z dimensions. The algorithm then scans the coarser mask, identifying which grains contain active threads and launches only this amount of threads. We will use the notation dim x ×dim y ×dim z to refer to different grains used.

IV. EXPERIMENTS

Several benchmarks were performed to evaluate the performance of our algorithm. In the following sections we describe our test environment and input instances used for the experiments.

A. Testing environment

The platform used for the CPU benchmarks was an Intel CoreTM i7 CPU model 930 with 4 cores running at 2.89Ghz and 12 GB memory. Despite the multi-core architecture, the CPU implementation is strictly sequential. The results of our GPU algorithm were obtained on an NVIDIA GPU GTX480 with 1.5 GBytes of global memory and 15 Multiprocessors with 32 cores each, totaling 480 CUDA cores. The CPU codes were compiled with GCC 4.8 using -O2 optimization flags. The CUDA driver is version 6 while the run-time is version 5.5.

B. Input instances

The input data set used differs on 3 different parameters: volume size, material map topology and number of Voronoi seeds. For the material map topologies we considered both synthetic and real-application data. The synthetic topologies represents a cube volume with (a) an uniform constant stiffness and (b) a gradient stiffness varying uniformly from left to right (called Gradient-Fig. 1a). In these topologies, we variate the volume from 32 3 to 256 3 voxels, which are the common discretization sizes used for physics simulation in [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF]. The seeds are randomly distributed on each map. We note that, due to the compliance-scaled distance function employed, the same set of seeds actually generate very different Voronoi diagrams, depending on the topology of the material map used (see Fig. 11).

The real-application data-set is the discretized material map of the T-bone steak (Fig. 9) from paper [START_REF] Faure | Sparse meshless models of complex deformable solids[END_REF]. The map of the steak has a volume of size 64x64x15 voxels and exhibits nonuniform stiffness distribution. The data of the steak is freely available for download with the SOFA framework [START_REF] Sofa | Simulation open framework architecture[END_REF], [START_REF] Faure | SOFA: A Multi-Model Framework for Interactive Physical Simulation[END_REF].

V. RESULTS AND DISCUSSION

This section presents the results obtained on several test cases. For a proper analysis, we divide our experimental results in three parts presented bellow.

A. Base algorithm speed-up

We start by comparing our parallel algorithm with its sequential reference implementation on CPU. The results are shown in the form of parallel speed-ups on Fig. 13. Each bar represents a different input instance where labels cube32 3 10s, cube64 3 10s, cube128 3 10s and cube256 3 10s denote a cube with gradient topology with dimensions 32, 64, 128 and 256 respectively, each with 10 seeds. The plate 100x40x10 20s input is a plate of stiffness gradient with 20 seeds (shown in Fig. 1). Both steak instances have a bounding volume of 64x64x15 voxels. In this benchmark the speedup achieved varies from 3.8x for small volumes (Steak) up to almost 40x for bigger ones. These results show that our algorithm benefits from bigger input sizes, because they expose more parallelism. This is also confirmed in Fig. 12, where we note that the maximum amount of active parallel threads is higher for bigger volumes. We summarize the execution times obtained with this benchmark in Table I .

B. Voronoi seeds

On a second scenario, we investigated the impact of the number of seeds in the performance. We used a volume size of 128 3 and same gradient topology. For each quantity of seeds, we randomly generated 10 different seed sets. Each result shown on Fig. 14 is the average speed-up obtained with these 10 instances. These results suggest that for larger amounts of seeds the speedup increases. Indeed, having more Voronoi seeds has the effect of allowing more active threads at the first iteration. Moreover, the number of iterations of the algorithm tends to reduce as more Voronoi cells expand concurrently.

C. Stream compression

Our last set of experiments evaluates the stream compression optimization described in section III-C. This optimization can be parametrized by setting the grain dimensions used for the subdivision of the coarser mask like shown on Fig. 10. We used the CUDA profiling tools to analyze the trade-off between overhead of stream compression and gained performance at several grains. We summarized the most representative results in the stacked histogram of Fig. 15. The figure presents results for a volume size of 128 3 with gradient topology and a set of 10 fixed seeds. The bars are sorted by total execution time and each grain size is indicated on the x axis, where "static" refers to the base algorithm. Finner grains generate larger masks, therefore add more overhead for generating the map of active threads. A finner grain, however, provides better compression, which reduces the amount of idle threads and of useless thread-blocks. This positively impacts the time spent on the RELAXKERNEL procedure. See column 1 × 1 × 1 in Fig. 15.

On the other-hand, more compact thread-blocks will also increase the number of threads accessing non-coalesced memory locations. In our application the volume is stored linearly in memory, which means that neighbour voxels on the x dimension are stored contiguously in memory. Favoring a larger x grain-dimension increases the number of threads accessing the same memory segment, thus achieving a better memory throughput. This fact can be observed comparing grains of same sizes, but different shapes like 4 × 1 × 1 and 1 × 1 × 4.

Experiments with smaller volume sizes, like 32 3 and 64 3 , showed worse total execution time than the base algorithm. Nevertheless, for the 128 3 volume, the technique of stream compression led to a 23.29 % performance gain over the base implementation.

VI. CONCLUSION AND FUTURE WORK

In this work we presented a GPU algorithm for computing the Voronoi diagram with generalized distance functions. Our method adapts a graph algorithm, for the SSSP problem, to compute the Voronoi diagram on a 3D grid of voxels. We have shown through experimental evaluation that our base parallel implementation significantly speeds-up Voronoi computation. Additionally, we applied an optimization strategy called Stream compression that allows to increase utilization of the GPU on large volumes.

Regarding our implementation, there is still room for optimization on the data representation in memory. Its current linear representation cannot benefit from the locality present in the neighborhood computation. A better memory layout, like Z-curves, could further enhance the performance.

Our algorithm has a direct application on physics-based simulation algorithms. As a next step towards this direction, we plan to address the dynamic scenario simulating cuts in deformable objects. In this case, local updates of the Voronoi diagram would be needed to handle dynamic changes in topology.

Fig. 2 .

 2 Fig. 2. Relaxation Algorithm

Fig. 3 .Fig. 4 .

 34 Fig. 3. Encoding information of distance and Voronoi index in a single word. Values d and k can be changed to adjust precision.

Fig. 5

 5 Fig. 5. Host Code 1: procedure RELAXKERNEL 2:

  )

6 :

 6 AtomicMin((C 1 [nid]|Vor [nid]), 7: (d new |Vor [tid]))

Fig. 6 .

 6 Fig. 6. Relaxation Kernel: updates the current shortest path estimates and the closest Voronoi seed.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Update kernel: verifies the termination condition and updates the activity mask.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Use case example: The T-bone steak (a) contains a mixture of flexible meat, softer grease and a rigid bone. As input we take the voxelized material map of stiffness values (b) and the coordinates of the simulation nodes. Our method computes the Voronoi diagram (d) rooted at each node using a compliance-scaled distance metric (c).

Fig. 11 .

 11 Fig. 11. Comparison of Voronoi diagrams generated with the same set of seeds on two different material maps. Left: with an uniform stiffness. Right: with a stiffness gradient.

Fig. 12 .

 12 Fig. 12. Thousands of active threads at each iteration of RELAXKERNEL.

Fig. 13 .Fig. 14 .Fig. 15 .

 131415 Fig.[START_REF] Nasre | Atomic-free irregular computations on GPUs[END_REF]. Speed-up for different input sizes. Gradient and constant topologies are presented for synthetic benchmarks only. Steak's topology corresponds to the real data-set of Fig.9

TABLE I BENCHMARK

 I RESULTS WITH EXECUTION TIMES AND SPEED-UP OBTAINED WITH THE BASE GPU ALGORITHM.

				CPU			GPU
	Topology	Volume	#Seeds	Iterations	Time (ms)	Iterations	Time(ms)	Speed-up
		32 3	10	32768	32.24	30	1.96	16.46
		64 3	10	262144	296.77	70	12.13	24.46
	Gradient	128 3	10	2097152	3863.54	126	105.91	36.48
		256 3	10	16777216	38756.80	195	985.80	39.31
		100 × 40 × 20	20	80000	77.04	21	2.51	30.68
		32 3	10	32768	21.37	20	1.27	16.80
		64 3	10	262144	190.81	52	9.20	20.73
	Constant	128 3	10	2097152	1957.65	79	61.85	31.65
		256 3	10	16777216	22103.64	159	812.03	27.22
		100 × 40 × 20	20	80000	52.31	22	2.47	21.15
	Steak	64 × 64 × 15 64 × 64 × 15	3 10	12276 12276	10.80 10.39	38 45	2.16 2.72	5.01 3.82

ACKNOWLEDGMENT

We would like to thank Franc ¸cois Faure from INRIA-Grenoble for its insightful help and discussion about the method of material-aware distance function. We also thanks CAPES and CNPq (process 476685/2012-5 and 309483/2011-5) for the financial support.


Michail Misyrlis

Anna B Konova

Matthew B Blaschko

Jean Honorio

Nelly Alia-Klein

Rita Z Goldstein

Dimitris Samaras

Predicting cross-task behavioral variables from fMRI data using the k-support norm

Keywords: Functional magnetic resonance Imaging (fMRI), Regularization, Sparse representations

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Functional magnetic resonance imaging (fMRI) is a widely used modality within the field of neuroimaging, that measures brain activity by detecting associated changes in blood oxygenation. One of the goals of fMRI data analysis is to detect correlations between brain activation and a task the subject performs during the scan.

The main challenges in statistical fMRI data analysis [START_REF] Honorio | Can a single brain region predict a disorder?[END_REF][START_REF] Song | Comparative study of svm methods combined with voxel selection for object category classification on fmri data[END_REF][START_REF] Bartels | Brain dynamics during natural viewing conditionsa new guide for mapping connectivity in vivo[END_REF][START_REF] D R Hardoon | Unsupervised analysis of fmri data using kernel canonical correlation[END_REF] are (i) the curse of dimensionality (ii) a small number of samples, due to the high cost of fMRI acquisition, and (iii) high levels of noise, such as system noise and random neural activity.

Sparsity regularizers are key statistical methods for improving predictive performance in the event that the number of observations is substantially smaller than the dimensionality of the data, as is the case in fMRI analysis. In this paper we compare the most frequently applied sparsity regularizer developed in the statistics literature, LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], with the k-support norm [START_REF] Argyriou | Sparse prediction with the k-support norm[END_REF], a recently introduced method which is less biased towards sparse solutions.

The k-support norm can be viewed as a generalization of LASSO when k = 1 and ridge regression when k = d, where d is the dimensionality of the data. The k-support norm has previously been used in [START_REF] Argyriou | Sparse prediction with the k-support norm[END_REF] for classification. It was first used for fMRI data modelling in [START_REF] Gkirtzou | fmri analysis of cocaine addiction using k-support sparsity[END_REF] with a specific choice k parameter of the norm, that is, the method was not tested in a model selection setting. In both cases the k-support norm was used with the squared loss function.

We focus on comparing LASSO with the k-support norm in order to establish the latter regularizer's superiority in analyzing fMRI data. We use two loss functions, namely the squared error and the absolute error. The advantage of the absolute error loss is that it is more robust, in that it penalizes outliers less than squared loss, while still retaining convexity which guarantees finding the global optimum. We compare the methods not only in their predictive accuracy but also in the interpretability and stability of their results.

Our contribution in this paper is threefold. First, we introduce a novel method, the k-support norm with absolute error. Second, this is the first attempt to compare LASSO with the k-support norm in a principled model selection setting. Finally, to the best of our knowledge this is the first application of the k-support norm to a real valued response variable in a challenging clinical setting where the fMRI signal collected during one task is used to predict behavioral responses collected at a different time during a second task.

The neuroscientific motivation for our experiments is the exploration of human drug addiction. Basic studies have led to a theoretical model of human drug addiction, characterized by Impaired Response Inhibition (RI) and Salience Attribution (SA) (hence, I-RISA) [START_REF] Rz Goldstein | Drug addiction and its underlying neurobiological basis: neuroimaging evidence for the involvement of the frontal cortex[END_REF]. According to the model, the skew in SA is predictive of impaired RI, together contributing to excessive drug use and relapse, core clinical symptoms of cocaine addiction. We use the fMRI data from a SA task (drug Stroop) in order to predict behavioral data in a RI task (colorword Stroop) collected at a different time, hence providing further evidence to support the I-RISA model.

Methods

We denote by X ∈ R n×d the design matrix of n samples each with d dimensions; we denote by y ∈ R n the vector of targets.

A basis of statistical inference is the application of regularized risk, in which a loss function is evaluated over a sample of data and is linearly combined with a regularizer that penalizes some norm of the prediction function as in (Eq. ( 1)), where the first term is the loss function and the second is the penalty term:

min β f (β, X, y) + λJ(β). (1)
The scalar parameter λ > 0 controls the degree of regularization and J is a scalar valued function monotonic in a norm of β ∈ R n . Sparsity regularization is a key family of priors over linear functions that prevents overfitting and aids interpretability of the resulting models [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Argyriou | Sparse prediction with the k-support norm[END_REF]. Key to the mathematical understanding of sparsity regularizers is their interpretation as convex relaxations to quantities involving the ℓ 0 norm, which simply counts the number of non-zero elements of a vector. One of the most important sparsity regularizers is the LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], where λJ(β) = λ β 1 . In many learning problems of interest, LASSO has been observed to shrink too many of the β variables to zero. In the presence of a group of highly correlated variables, LASSO may prefer a sparse solution. However including all correlated variables in the model could potentially lead to higher predictive accuracy [START_REF] Argyriou | Sparse prediction with the k-support norm[END_REF] and the k-support norm provides a way of calibrating the cardinality of the regression vector β so as to include more variables.

The k-support norm can be computed as

β sp k =   k-r-1 i=1 (|β| ↓ i ) 2 + 1 r + 1 d i=k-r |β| ↓ i 2   1 2 (2)
where |β| ↓ i is the ith largest element of the vector and r is the unique integer in {0, . . . , k -1} satisfying

|β| ↓ k-r-1 > 1 r + 1 d i=k-r |β| ↓ i ≥ |β| ↓ k-r . (3)
In this paper, we consider LASSO and the k-support norm with two loss functions: the squared error f (β, X, y) = y -Xβ 2 2 and the absolute error

f (β, X, y) = y -Xβ 1 .
In practice, we approximate the absolute error with a Huber type smoothing around zero to ensure differentiability.

Experimental Set-up

In this section we present our experiments and the data sets used in them. Our experiments aim at providing empirical evidence for the support of the I-RISA model.

We use the fMRI drug-word task described in [START_REF] Rz Goldstein | Role of the anterior cingulate and medial orbitofrontal cortex in processing drug cues in cocaine addiction[END_REF][START_REF] Ab Konova | Structural and behavioral correlates of abnormal encoding of money value in the sensorimotor striatum in cocaine addiction[END_REF]. The neuropsychological experiment for cocaine addiction data set has a block design, which includes eight sessions, with each of them having different conditions. The two varying conditions are the monetary reward (50➣, 25➣, 1➣ and 0➣) and the cue shown (drug words, neutral words). The session consists of an initial screen displaying the monetary reward and then presenting a sequence of forty words in four different colors (yellow, blue, red or green). The subject was instructed to press one of four buttons matching the color of the word they had just read. The subjects were rewarded for correct performance depending on the monetary condition.

We use the behavioral responses of the same subjects in a color-word task [START_REF] Sj Moeller | Dopaminergic involvement during mental fatigue in health and cocaine addiction[END_REF], a classic task of inhibitory control. In this task the subjects pressed for ink color of color words printed in either their congruent or incongruent colors. Four colors and words (red, blue, yellow and green) were used in all possible combinations. Both congruent and incongruent stimuli were presented randomly. The subjects performed four consecutive runs of this task. As there were 12 incongruent events in each run of 200 events, each subject's data contained up to 48 incongruent events.

For 38 control subjects and 74 cocaine abusers, we use the fMRI data from the drug-word task, to predict color-word behavioral variables.

Our experimental setting consists of 500 trials with an 85% / 15% random split between training and test sets. We perform model selection on the training set. That is, for each combination of parameters (λ

∈ {10 i : i = -2, • • • , 8} for LASSO, λ ∈ {10 i : i = -2, • • • , 8}, k ∈ {1, 2, 3, 6, 12
, 100, 200, 300, 600} for k-support norm), we do a leave-one-subject-out cross validation on the samples that constitute the training set. We measure the correlation between the predicted and the true response variables on the training set. The parameter setting that leads to the highest correlation is used on the whole training set in order to learn a set of weights for each method, which are then applied on the test set. Finally, we measure the correlation between the predicted and the true response variables on the test set. We report the mean correlation on the holdout test samples and its standard error across the 500 random permutations in Sec. [START_REF] D R Hardoon | Unsupervised analysis of fmri data using kernel canonical correlation[END_REF]. We note that the same sample randomization is used for both LASSO and ksupport norm.

In experiment 1 we use the fMRI contrast drug > neutral words, averaged over monetary reward condition, to predict the conflict effect in the subjects' reaction time on the color-word task, defined as the difference in time between correctly performing the task for congruent and incongruent events. We use the Insula, Hippocampus Complex, Amygdala and ACC, part of the brain's limbic (emotion) circuit, as regions of interest (ROIs) for this experiment. These regions are chosen on the basis of previous studies on independent datasets that showed limbic system modulation by drug-related cues, eg. drug words [START_REF] Hw Chase | The neural basis of drug stimulus processing and craving: an activation likelihood estimation meta-analysis[END_REF].

In experiment 2 we use the fMRI contrast 50➣ > 0➣, averaged over word type condition, in order to predict the subjects' responses on the color-word task, defined as the difference in percent accuracy between performing the task for congruent and incongruent events. We use the Basal Ganglia and Thalamus, part of the brain's reward circuit, as ROIs for this experiment. We chose these ROIs on the basis of previous studies on independent datasets that showed reward system modulation by primary and secondary reinforcers, including money [START_REF] Liu | Common and distinct networks underlying reward valence and processing stages: a meta-analysis of functional neuroimaging studies[END_REF].

Results

We compare the performance of the two methods in Table 1 for the first experiment and Table 2 1. Mean (SE) correlation over 500 random permutations of the samples between the predicted and the actual conflict effect on the reaction times for drug > neutral using the limbic ROI, for all combinations of regularizers and loss functions. The pvalues were computed with a Wilcoxon signed rank test between the 500 correlation values for the two combinations of regularizer and loss function in the preceding rows or columns. Based on the p-values, there is a statistically significant difference between absolute loss predictions and squared loss predictions and between LASSO and ksupport norm with the squared loss function in both cocaine and control subjects.

With the squared loss function, the k-support norm outperforms LASSO for almost all cases, while when combined with the absolute loss function, the regularizers do not significantly differ in their predictive performance. The absolute loss function, for both regularizers, leads to correlations that are significantly higher than those with the squared loss function in almost all cases.

We report the fraction of non-zero weights that were selected by each method for over 50% of the 500 trials in Tables 3 and4 for the first and the second experiment respectively.

We average the weights assigned to the voxels over the 500 permutations and then compute the cumulative distribution function (CDF) for those weights. We threshold the CDF at 0.9 and visualize the weights of the voxels up to that threshold 1 in Fig. 1. The overly sparse solutions of the LASSO (Fig. 1(b),1(d)) lead to models that cannot be interpreted as easily as the solutions of the ksupport norm method (Fig. 1(a),1(c)).

In the presence of correlated features, the degree of sparsity of the solution can be tuned with the k-support norm in order to include several highly correlated features. In contrast, LASSO tends to pick one representative feature with 2. Mean (SE) correlation over 500 random permutations of the samples between the predicted and the actual response variables for 50➣> 0➣ using the Basal Ganglia, Thalamus ROI, for all combinations of regularizers and loss functions. The p-values were computed with a Wilcoxon signed rank test between the 500 correlation values for the two combinations of regularizer and loss function in the preceding rows or columns. Based on the p-values there is a statistically significant difference between absolute loss predictions and squared loss predictions and between k-support and LASSO with the squared loss in cocaine subjects only. 3. Voxel Selection stability over 500 random permutations of the samples for drug > neutral using the limbic ROI, for all combinations of regularizers and loss functions. The fraction of voxels which are selected for more than 50% of the 500 trials are presented. The higher values reported for k-support norm indicate that it makes more stable voxel selection than LASSO over different training sets.

no guarantee of consistency in feature selection across different splits of the data samples into training and test sets. In all cases the fraction of non-zero weights selected by the k-support norm is higher than that of LASSO, indicating that the k-support norm method leads to more stable solutions as compared to those obtained with LASSO. 4. Voxel Selection stability over 500 random permutations of the samples for 50➣> 0➣ using the Basal Ganglia, Thalamus ROI, for all combinations of regularizers and loss functions. The fraction of voxels which are selected for more than 50% of the 500 trials are presented. The higher values reported for k-support norm indicate that it makes more stable voxel selection than LASSO over different training sets.

Discussion

In our experiments, in almost all cases, the k-support norm outperforms LASSO in predicting the behavioral measures given fMRI data when combined with squared loss, while when combined with the absolute loss, the predictive accuracy of the two regularizers does not differ significantly. The absolute loss led to higher predictions than squared loss for both regularizers for almost all cases. The LASSO leads to sparse solutions, since it tends to pick one feature per group of correlated features. On the other hand, the k-support norm allows calibrating the cardinality of the solutions and thus can select more interpretable groupings of correlated features and also leads to more stable results across different training sets. Thus, our results support the further exploration of the k-support norm for fMRI analysis.

We also provide further evidence to support the I-RISA model of drug addiction, whereby the skew in SA in cocaine abusers, as indexed by fMRI response to drug words and monetary rewards, two motivationally salient stimuli, is predictive of RI, as indexed by response slowing and accuracy on a task requiring inhibitory control (the color-word Stroop). Specifically, we show that in cocaine users, response to drug words in voxels located in limbic brain regions, such as the anterior insula and ACC implicated in emotion processing and emotion regulation, was predictive of slower responses on the RI task (Exp. 1), while response to money in voxels located in reward-related brain regions, such as the putamen implicated in habits, was predictive of lower accuracy on the RI task (Exp. 2).

  (a) Most predictive voxels in Exp. 1 using the k-support norm with the Absolute Loss (b) Most predictive voxels in Exp. 1 using the LASSO with the Absolute Loss (c) Most predictive voxels in Exp. 2 using the k-support norm with the Squared Loss (d) Most predictive voxels in Exp. 2 using the LASSO with the Squared Loss

Fig. 1 .

 1 Fig. 1. Visualization of the most predictive voxels in Exp. 1 (top row) and Exp. 2 (bottom row) over the 500 permutations. The degree of sparsity of the solution can be tuned with the k-support norm, thus leading to models ((a), (c)) that are easier to interpret than those of LASSO ((b), (d)). (Best viewed in color)

  for the second experiment.

	Mean Correlation, D>N, Conflict effect on Reaction Time
		Control Subjects
	Norm / Loss Squared Absolute p
	LASSO	0.16 (0.02) 0.27 (0.02) <0.01
	k-support	0.22 (0.02) 0.24 (0.02) <0.05
	p	<0.001	0.21
		Cocaine Subjects
	Norm / Loss Squared Absolute p
	LASSO	0.27 (0.01) 0.37 (0.01) <0.001
	k-support	0.33 (0.01) 0.36 (0.02) <0.001
	p	<0.001	0.96
	Table

1

  Due to space constraints we include one representative example out of two for each experiment. The omitted results are qualitatively similar.

	Mean Correlation, 50➣>0➣, Conflict effect on Accuracy
		Control Subjects
	Norm / Loss Squared Absolute p
	LASSO	0.25 (0.02) 0.09 (0.02) <0.001
	k-support	0.26 (0.02) 0.09 (0.02) <0.001
	p	0.42	0.78
		Cocaine Subjects
	Norm / Loss Squared Absolute p
	LASSO	0.22 (0.02) 0.42 (0.02) <0.001
	k-support	0.27 (0.01) 0.41 (0.02) <0.001
	p	<0.001	0.78
	Table