Damien Chablat

Philippe Wenger

Stefan Staicu

DYNAMICS OF THE ORTHOGLIDE PARALLEL ROBOT

Keywords: Dynamics, Kinematics, Parallel robot, Virtual work

Articolul stabileşte relaţii matriceale pentru cinematica şi dinamica robotului paralel Orthoglide prevăzut cu trei acţionori prismatici concurenţi. Aceştia sunt aranjaţi în raport cu sistemul cartezian de coordonate astfel încât direcţiile lor să fie normale unele faţă de celelalte. Trei lanţuri cinematice identice, conectate la platforma mobilă, sunt localizate în trei plane perpendiculare unul pe celălalt. Cunoscând poziţia şi mişcarea de translaţie a platformei, se dezvoltă problema de cinematică inversă şi se determină poziţia, viteza şi acceleraţia fiecărui element al robotului. În continuare, principiul lucrului mecanic virtual este folosit în problema de dinamică inversă. Câteva ecuaţii matriceale oferă expresii recurente şi grafice pentru forţele active şi puterile mecanice ale celor trei acţionori.

Recursive matrix relations for kinematics and dynamics of the Orthoglide parallel robot having three concurrent prismatic actuators are established in this paper. These are arranged according to the Cartesian coordinate system with fixed orientation, which means that the actuating directions are normal to each other. Three identical legs connecting to the moving platform are located on three planes being perpendicular to each other too. Knowing the position and the translation motion of the platform, we develop the inverse kinematics problem and determine the position, velocity and acceleration of each element of the robot. Further, the principle of virtual work is used in the inverse dynamic problem. Some matrix equations offer iterative expressions and graphs for the input forces and the powers of the three actuators.

,  k k a

: orthogonal relative transformation matrix

Introduction

Generally, the mechanism of a parallel robot has two platforms: one of them is attached to the fixed reference frame and the other one can have arbitrary motions in its workspace. Some movable legs, made up as serial robots, connect the moving platform to the fixed platform. Typically, a parallel mechanism is said to be symmetrical if it satisfies the following conditions: the number of legs is equal to the number of degrees of freedom of the moving platform, one actuator, which can be mounted at or near the fixed base, controls every limb and the location and the number of actuated joints in all the limbs are the same (Tsai [START_REF] Tsai | Robot analysis: the mechanics of serial and parallel manipulator[END_REF]).

For two decades, parallel manipulators attracted to the attention of more and more researches that consider them as valuable alternative design for robotic mechanisms [START_REF] Asada | Robot Analysis and Control[END_REF], [START_REF] Fu | Robotics: Control, Sensing, Vision and Intelligence[END_REF], [START_REF] Craig | Introduction to Robotics: Mechanics and Control[END_REF]. As stated by a number of authors [START_REF] Tsai | Robot analysis: the mechanics of serial and parallel manipulator[END_REF], conventional serial kinematical machines have already reached their dynamic performance limits, which are bounded by high stiffness of the machine components required to support sequential joints, links and actuators.

The parallel robots are spatial mechanisms with supplementary characteristics, compared with the serial architecture manipulators such as: more rigid structure, important dynamic charge capacity, high orientation accuracy, stabile functioning as well as good control of velocity and acceleration limits. However, most existing parallel manipulators have limited and complicated workspace with singularities and highly non-isotropic input-output relations [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods and Algorithms[END_REF].

Research in the field of parallel manipulators began with the most known application in the flight simulator with six degrees of freedom, which is in fact the Stewart-Gough platform (Stewart [6]; Merlet [START_REF] Merlet | Parallel robots[END_REF]; Parenti-Castelli and Di Gregorio [START_REF] Parenti-Castelli | A new algorithm based on two extra-sensors for real-time computation of the actual configuration of generalized Stewart-Gough manipulator[END_REF]). The Star parallel manipulator (Hervé and Sparacino [START_REF] Hervé | Star. A New Concept in Robotics[END_REF]) and the Delta parallel robot (Clavel [10]; Tsai and Stamper [START_REF] Tsai | A parallel manipulator with only translational degrees of freedom[END_REF]; Staicu [START_REF] Staicu | Recursive modelling in dynamics of Delta parallel robot[END_REF]) equipped with three motors, which have a parallel setting, train on the effectors in a threedegrees-of-freedom general translation motion.

While the kinematics has been studied extensively during the last two decades, fewer papers can be focused on the dynamics of parallel robots. When good dynamic performance and precise positioning under high load are required, the dynamic model is important for their control. The analysis of parallel robots is usually implemented trough analytical methods in classical mechanics [START_REF] Li | Inverse dynamics and simulation of a 3-DOF spatial parallel manipulator[END_REF], in which projection and resolution of equations on the reference axes are written in a considerable number of cumbersome, scalar relations and the solutions are rendered by large scale computation together with time consuming computer codes. Geng [START_REF] Geng | On the dynamic model and kinematic analysis of a class of Stewart platforms[END_REF] developed Lagrange's equations of motion under some simplifying assumptions regarding the geometry and inertia distribution of the manipulator. Dasgupta and Mruthyunjaya [START_REF] Dasgupta | A Newton-Euler formulation for the inverse dynamics of the Stewart platform manipulator[END_REF] used the Newton-Euler approach to develop closed-form dynamic equations of Stewart platform, considering all dynamic and gravity effects as well as viscous friction at joints. However, to the best of our knowledge, these are no efficient dynamic modelling approach available for parallel manipulators. In recent years, several new kinematical structures have been proposed that possess higher isotropy [START_REF] Carricato | Singularity-Free Fully-Isotropic Translational Parallel Mechanisms[END_REF], [START_REF] Kong | A Class of 3-DOF Translational Parallel Manipulators with Linear I-O Equations[END_REF], [START_REF] Kim | Design Optimisation of a Cartesian Parallel Manipulator[END_REF], [START_REF] Zanganeh | Dynamics of a six-degree-of-freedom parallel manipulator with revolute legs[END_REF], [START_REF] Xi | Effect of dynamic balancing on four-bar linkage vibrations, Mechanism and Machine Theory[END_REF].

The objective of this paper is to analyse the kinematics and dynamics of the Orthoglide parallel robot, which is well adapted to the applications of precision assembly machines. In design, the three actuators are arranged according to the Cartesian coordinate space, which means that the actuating directions are normal to each other and the joints connecting to the moving platform are located on three planes being perpendicular to each other too. Proposed by Wenger and Chablat [START_REF] Pashkevich | Kinematics and workspace analysis of a three-axis parallel manipulator: the Orthoglide[END_REF], [START_REF] Wenger | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF], the prototype of the manipulator has good kinetostatic performance and some technological advantages such as: symmetrical design, regular workspace shape properties with a bounded velocity amplification factor and low inertia effects.

In the present paper we focus our attention on a recursive matrix method, which is adopted to derive the kinematics model and the inverse dynamics equations of the spatial Orthoglide parallel robot [START_REF] Pashkevich | Design Strategies for the Geometric Synthesis of Orthoglide-type Mechanisms, Mechanism and Machine Theory[END_REF], which has three translation degrees of freedom (Fig. 1).

Inverse kinematics

The mechanism input of the manipulator is made up of three actuated orthogonal prismatic joints. The output body is connected to the prismatic joints through a set of three identical kinematical chains.

The architecture of one of the three parallel closed chains of the Orthoglide manipulator is formally described as R PRP a mechanism, where R P, and a P denote the prismatic, revolute and parallelogram joints, respectively. So, the topological structure consists in an active prismatic system, a passive revolute joint, an intermediate mechanism with four revolute links that connect four bars, which are parallel two by two, ending with a passive revolute link connected to the moving platform. Inside each chain, the parallelogram mechanism is used and oriented in a manner that the end-effector is restricted to translation movement only. The arrangement of the joints in the chains has been defined to eliminate any constraint singularity in the Cartesian workspace [START_REF] Wenger | Architecture Optimization of a 3-DOF Parallel Mechanism for Machining Applications, the Orthoglide[END_REF], [START_REF] Pashkevich | Design Strategies for the Geometric Synthesis of Orthoglide-type Mechanisms, Mechanism and Machine Theory[END_REF], [START_REF] Liu | A Kind of Three Translational-DOF Parallel Cube-Manipulator[END_REF]. at the intersection point of three axes of actuated prismatic joints, about which the three-degrees-of-freedom manipulator moves. It has three legs of known dimensions and masses. To simplify the graphical image of the kinematical scheme of the mechanism, in the follows we will represent the intermediate reference systems by only two axes, so as is used in most of books [START_REF] Tsai | Robot analysis: the mechanics of serial and parallel manipulator[END_REF], [START_REF] Craig | Introduction to Robotics: Mechanics and Control[END_REF], [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods and Algorithms[END_REF], [START_REF] Merlet | Parallel robots[END_REF]. The k z axis is represented, of course, for each component element k T . We mention that the relative rotation or relative translation with

P x z y

1 ,  k k  angle or 1 ,  k k  displacement of k
T body most be always pointing about or along the direction k z . The first element of leg A is one of the three sliders of the robot. It is a homogenous rod of length  is denoted as 2 A . This link is connected to the frame

A A A z y x A 2 2 2 2 (called A T 2 )
and it has a relative rotation with the angle A 21  , so that

A A 21 21     and A A 21 21     
. It has the mass 2 m and the central tensor of inertia 2 Ĵ .

Further one, two identical and parallel bars Due to the special arrangement of the four-bar parallelograms and the three prismatic joints at points

1 1 1 , , C B A
, the mechanism has three translation degrees of freedom. This unique characteristic is useful in many applications, such as a z y x   positioning device. At the central configuration, we consider that all legs are initially extended at equal length and that the angles giving the orientation of the three sliders about their guide-ways are

       C B A .
In the followings, we apply the method of successive displacements to geometric analysis of closed-loop chains and we note that a joint variable is the displacement required to move a link from the initial location to the actual position. If every link is connected to least two other links, the chain forms one or more independent closed-loops. The variable angles where we denoted

                        0 0 1 0 1 0 1 0 0 , 0 0 1 0 1 0 1 0 0 2 1 a a ,              0 1 0 0 0 1 1 0 0 3 a              1 0 0 0 1 0 0 0 1 4 a , , 0 1 0 1 0 0 0 0 1 5            a             1 0 0 0 0 1 0 1 0 6 a                  1 0 0 0 cos sin 0 sin cos 1 , 1 , 1 , 1 , 1 , A k k A k k A k k A k k k k a      ,  . 5 ,..., 2 , 1 1 , 1 0        k a a k j j k j k k (2)
The translation conditions for the platform are given by the following identities

I c c b b a a T T T    50 50 50 50 50 50    , (3) where
, 0 0 1 0 1 0 1 0 0 , 1 0 0 0 0 1 0 1 0 50 50                               b a . 1 0 0 0 1 0 0 0 1 50               c (4)
From these relations, one obtains the following relations between angles

C C B B A A 21 54 21 54 21 54 , ,          . ( 5
)
The three concurrent displacements x 0 , G y 0 , G z 0 . Supposing, for example, that the known motion of the mass center G of the platform is expressed by the following relations

] [ 0 0 0 0 G G G G z y x r   ), 3 cos 1 ( ), 3 cos 1 ( ), 3 cos 1 ( * 0 0 * 0 0 * 0 0 t z z t y y t x x G G G G G G          (6) the inputs C B A 10 10 10 , ,   
of the manipulators and the variables

C C B B A A 32 21 32 21 32 21 , , , , ,      
will be given by the following geometrical conditions

G k GC T C k , k T k C k GB T B k k T k B k GA T A k , k T k A r r c r c r
                            , (7)
where, for example, one denoted

                                             0 0 0 0 0 1 0 1 0 , 1 0 0 , 0 1 0 , 0 0 1 3 3 2 1 u u u u    3 2 32 1 1 21 3 10 3 1 10 10 2 , ] cos sin 0 [ ) 2 cos ( u l r l l r u a l l l r A T A T A A                 (8) , 2 , 1 2 54 2 3 43 u l r u l r A A        T GA l l r ] 0 2 sin [ 1 5     .
Actually, these equations means that there is the inverse geometric solution for the manipulator, given through following analytical relations

C C G C l z       .
In that follows, we determine, the velocities and the accelerations of the robot, supposing that the translation motion of the moving platform is known.

The motions of the component elements of each leg (for example the leg A) are characterized by the following skew symmetric matrices [START_REF] Staicu | A novel dynamic modelling approach for parallel mechanisms analysis[END_REF] )

5 ,..., 2 ( , ~3 1 , 1 , 0 , 1 1 , 0        k u a a A k k T k k A k k k A k    , ( 10
)
which are associated to the absolute angular velocities given by the recurrence relations . ,

1 , 1 , 3 1 , 0 , 1 1 , 0 A k k A k k A k k A k k k A k u a                  (11)

Following relations give the velocities

A k v 0  of the joints k A   A - k k, A k A - k k k A k r v a v 1 0 , 1 1,0 1 , 0 ~        , 3 10 10 u v A A      . ( 12
)
If the other two kinematical chains of the manipulator are pursued, analogous relations can be easily obtained. Equations ( 3) and ( 7) can be differentiated with respect to time to obtain the following matrix conditions of connectivity [START_REF] Staicu | Power requirement comparison in the 3-RPR planar parallel robot dynamics[END_REF] , ) result as functions of the translation velocity of the platform. The relations [START_REF] Li | Inverse dynamics and simulation of a 3-DOF spatial parallel manipulator[END_REF] give the complete Jacobian matrix of the manipulator. This matrix is a fundamental element for the analysis of the robot workspace and the particular configurations of singularities where the manipulator becomes uncontrollable.

      i r u u u a u l u a u a u l u a u v u a u u a u G T i T T i A T T T i A T T i A T T i A T T i A                  (13)
Rearranging, above nine constraint equations ( 9) of the Orthoglide robot can immediately written as follows 2 3 J and 2 J are, respectively, the inverse and forward Jacobian of the manipulator and can be expressed as

2 10 3 0 2 0 2 0 ) ( l l x y z A G G G       2 3 2 10 3 0 2 0 2 0 ) ( l l y z x B G G G       ( 14
) 2 3 2 10 3 0 2 0 2 0 ) ( l l z x y C G G G       ,
} { 3 2 1 1    diag J  ,            3 0 0 0 2 0 0 0 1 2    G G G G G G y x z x z y J , ( 16
)
with

A G l x 10 3 0 1      ; B G l y 10 3 0 2      ; C G l z 10 3 0 3      . ( 17
)
The three kinds of singularities of the three closed-loop kinematical chains can be determined through the analysis of two Jacobian matrices 1 J and 2 J . Let us assume that the robot has a first virtual motion determined by the linear velocities

0 , 1 10 10   Bv a Av a v v , 0 10  Cv a v
. The characteristic virtual velocities are expressed as functions of the position of the mechanism by the general kinematical constraints equations [START_REF] Li | Inverse dynamics and simulation of a 3-DOF spatial parallel manipulator[END_REF]. Other two sets of relations of connectivity can be obtained if one considers successively:

As for the relative accelerations

A A A 32 21 10 , ,    and A A 21 54   
of the manipulator, the derivatives with respect to time of the relations [START_REF] Li | Inverse dynamics and simulation of a 3-DOF spatial parallel manipulator[END_REF] give other following conditions of connectivity [START_REF] Staicu | Inverse dynamics of the 3-PRR planar parallel robot[END_REF] ).

           i u u a u a u l u u u a u l u a u u a u l r u u u a u l u a u a u l u a u u a u u a u T T T i A A T T i A A T T T i A A G T i T T i A T T T i A T T i A T T i A T T i A                                (18)
The angular accelerations A k 0   and the accelerations A k 0   of joints are given by some relations, obtained by deriving the relations [START_REF] Clavel | Delta: a fast robot with parallel geometry[END_REF], [START_REF] Tsai | A parallel manipulator with only translational degrees of freedom[END_REF] and ( 12):

3 1 , 0 , 1 1 , 1 , 3 1 , 0 , 1 1 , 0 ~u a a u a T k k A k k k A k k A k k A k k k A k                      3 1 , 0 , 1 1 , 1 , 3 1 , 3 3 1 , 1 , 1 , 0 , 1 0
, 1 0 , 1 1 , 0 0 0 2 ũ a a u u u a a T k k A k k k A k k A k k A k k A k k T k k A k A k A k k k A k A k A k                                   A k k A k A k A k A k k k A k r a 1 , 0 , 1 0 , 1 0 , 1 0 , 1 1 , 0 ~                 , 3 10 10 u A A        (19)
The relations ( 13), [START_REF] Kim | Design Optimisation of a Cartesian Parallel Manipulator[END_REF] represent the inverse kinematics model of the Orthoglide parallel robot. As application let us consider a manipulator, which has the following characteristics

s t l l m l m l m l m l m z m y m x G G G 2 , 4 , , 85 . 0 , 08 . 0 , 15 . 0 , 20 . 0 20 . 0 , 10 . 0 , 05 . 0 2 4 3 2 1 * 0 * 0 * 0               kg m kg m 2 . 0 , 35 . 0 2 1   , . , 15 , , 5 . 2

Inverse dynamics model

In the context of the real-time control, neglecting the frictions forces and considering the gravitational effects, the relevant objective of the dynamics is to determine the input forces, which must be exerted by the actuators in order to produce a given trajectory of the effector.

There are three methods, which can provide the same results concerning these actuating forces. The first one is using the Newton-Euler classic procedure [START_REF] Li | Inverse dynamics and simulation of a 3-DOF spatial parallel manipulator[END_REF], [START_REF] Dasgupta | A Newton-Euler formulation for the inverse dynamics of the Stewart platform manipulator[END_REF], [START_REF] Zanganeh | Dynamics of a six-degree-of-freedom parallel manipulator with revolute legs[END_REF], [START_REF] Guegan | Modélisation dynamique d'un robot parallèle à 3-DDL: l'Orthoglide[END_REF], the second one applies the Lagrange's equations and multipliers formalism [START_REF] Geng | On the dynamic model and kinematic analysis of a class of Stewart platforms[END_REF], [START_REF] Miller | The Lagrange-Based Model of Delta-4 Robot Dynamics[END_REF] and the third one is based on the principle of virtual work [START_REF] Tsai | Robot analysis: the mechanics of serial and parallel manipulator[END_REF], [START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods and Algorithms[END_REF], [START_REF] Staicu | Inverse dynamics of the HALF parallel manipulator with revolute actuators[END_REF], [START_REF] Staicu | A novel dynamic modelling approach for parallel mechanisms analysis[END_REF]. In the inverse dynamic problem, in the present paper one applies the principle of virtual work in order to establish some recursive matrix relations for the forces and the powers of the three active systems.

Three input spatial concurrent forces j f 10 and three powers required in a given motion of the moving platform will easily be computed using a recursive procedure. Some independent pneumatic or hydraulic systems that generate three input forces

    CA k A k A k A k A k A k inA k r m f    0 0 0 0 0         (20)
and the resulting moment of the forces of inertia

] [ 0 0 0 0 0 A k A k A k A k A k A k CA k A k inA k I I r m m             (
81 . 9 u a m f k A k A k    , 3 0 * 81 . 9 u a r m m k CA k A k A k    ) 6 , ... , 2 , 1 (  k . ( 22
, ~1 , 1 , 1 1 , 1 0 1 , 1 0 A k T k k A k k A k T k k A k A k A k T k k A k A k f a r m a m m f a f f                    (23)
where one denoted

A k inA k A k f f f        0 0 , A k inA k A k m m m        0 0 . ( 24
) Considering three independent virtual motions of the robot, all virtual displacements and virtual velocities should be compatible with the virtual motions imposed by all kinematical constraints and joints at a given instant in time. By intermediate of the complete Jacobian matrix expressed by the conditions of connectivity (13), the absolute virtual velocities

v k v k v 0 0 ,   
associated with all moving links are related to a set of independent relative virtual velocities

3 1 , 1 , u v k k v k k        .
Knowing the position and kinematics state of each link as well as the external forces acting on the robot, in that follow one apply the principle of virtual work for an inverse dynamic problem. The fundamental principle of the virtual work states that a mechanism is under dynamic equilibrium if and only if the total virtual work developed by all external, internal and inertia forces vanish during any general virtual displacement, which is compatible with the constraints imposed on the mechanism. Assuming that frictional forces at the joints are negligible, the virtual work produced by the forces of constraint at the joints is zero. The relations ( 23), [START_REF] Staicu | Inverse dynamics of the HALF parallel manipulator with revolute actuators[END_REF] represent the inverse dynamics model of the Orthoglide parallel robot.

Based on the algorithm derived from the above recursive relations, a computer program solve the inverse dynamics modelling of the robot, using the MATLAB software.

Assuming that the weights g m A k  of compounding rigid bodies constitute the external forces acting on the robot during its evolution, a numerical computation in the dynamics is developed, based on the determination of the three active powers

Conclusions

In the inverse kinematics analysis some exact relations that give in real-time the position, velocity and acceleration of each element of the parallel robot have been established in present paper. The dynamics model takes into consideration the masses and forces of inertia introduced by all component elements of the robot.

The new approach based on the principle of virtual work can eliminate all forces of internal joints and establishes a direct determination of the time-history evolution of forces and powers required by the actuators. The recursive matrix relations [START_REF] Staicu | Inverse dynamics of the HALF parallel manipulator with revolute actuators[END_REF] represent the explicit equations of the dynamics simulation and can easily be transformed in a model for the automatic command of the Orthoglide parallel robot. Also, the method described above is quit available in forward and inverse mechanics of all serial or parallel mechanisms, the platform of which behaves in translation, rotation evolution or general 6-DOF motion.

R E F E R E N C E S

  symmetric matrix of tensor of inertia of k T about the link-frame

Fig. 1

 1 Fig. 1 Orthoglide parallel robot Let us locate a fixed reference frame ) ( 0 0 0 0 T z y Oxat the intersection point of three axes of actuated prismatic joints, about which the three-degrees-of-freedom manipulator moves. It has three legs of known dimensions and masses. To simplify the graphical image of the kinematical scheme of the mechanism, in the follows we will represent the intermediate reference systems by only two axes, so as is used in most of books[START_REF] Tsai | Robot analysis: the mechanics of serial and parallel manipulator[END_REF],[START_REF] Craig | Introduction to Robotics: Mechanics and Control[END_REF],[START_REF] Angeles | Fundamentals of Robotic Mechanical Systems: Theory, Methods and Algorithms[END_REF],[START_REF] Merlet | Parallel robots[END_REF]. The k z axis is represented, of course, for each component element k T . We mention that the relative rotation or relative translation with

.Fig. 2

 2 Fig. 2 Kinematical scheme of first leg A of the mechanism

  the joint axes k z are the parameters needed to bring the next link from a reference configuration to the next configuration. We call the matrix In the study of the kinematics of robot manipulators, we are interested in deriving a matrix equation relating the location of an arbitrary k T body to the joint variables. When the change of coordinates is successively considered, the corresponding matrices are multiplied.So, starting from the reference origin O and pursuing the three legs B A, ,C we obtains the following transformation matrices[START_REF] Staicu | Inverse dynamics of the HALF parallel manipulator with revolute actuators[END_REF]

  that give the input vector 10   of the instantaneous position of the mechanism. But, the objective of the inverse geometric problem is to find the input vector 10   and the position of the robot with the given three absolute coordinates of the center G of the platform: G

.

  The derivative with respect to time of conditions (14) leads to the matrix equation

10 

 10 implements the suggested algorithm is developed in MATLAB to solve first the inverse kinematics of the Orthoglide parallel robot. For illustration, it is assumed that for a period of two second the platform starts at rest from a central configuration and is moving in a general translation. A numerical study of the robot kinematics is carried out by computation of the input displacements , , for example, of three prismatic actuators (Fig.3, 4, 5).

Fig. 3 10  10  10 

 3101010 Fig. 3 Input displacement A 10  of first actuator

1 ,

 1 control the motion of the three sliders of the robot. Now, the parallel mechanism can artificially be transformed in a set of three open four-bar parallelograms and taking their effects into account by introducing the corresponding constraint conditions. The first and more complicated open tree system includes the acting link and could comprise the moving platform.The force of inertia of an arbitrary rigid body A k T , for example

  21) are determined with respect to the centre of its fist joint k A . On the other hand, the wrench of two vectors A k f   and A k m   evaluates the influence of the action of the external and internal forces applied to the same element

Fig. 7

 7 Fig. 7 Input power B p 10 of second actuator

  time-history evolution of the input powers A p 10 (fig. 6), B p 10 (fig. 7), C p 10 (fig. 8) required by the actuators are plotted for a period of two second of platform's motion.

  ) Finally, two recursive relations generate the vectors

Prof., Institut de Recherche en Communications et Cybernétique de Nantes, France

Prof., Institut de Recherche en Communications et Cybernétique de Nantes, France

Prof., Département de Mécanique, Université Polytechnique de Bucarest, Roumanie, e-mail :stefanstaicu@yahoo.com


Putting eagle rays on the map by coupling aerial videosurveys and deep learning

Keywords: Automated species detection, Convolutional neural networks, Coral reefs, Elasmobranchs, New-Caledonia

Reliable and efficient techniques are urgently needed to monitor elasmobranch populations that face increasing threats worldwide. Aerial video-surveys provide precise and verifiable observations for the rapid assessment of species distribution and abundance in coral reefs, but the manual processing of videos is a major bottleneck for timely conservation applications. In this study, we applied deep learning for the automated detection and mapping of vulnerable eagle rays from aerial videos. A light aircraft dedicated to touristic flights allowed us to collect 42 h of aerial video footage over a shallow coral lagoon in New Caledonia (Southwest Pacific). We extracted the videos at a rate of one image per second before annotating them, yielding 314 images with eagle rays. We then trained a convolutional neural network with 80% of the eagle ray images and evaluated its accuracy on the remaining 20% (independent data sets). Our deep learning model detected 92% of the annotated eagle rays in a diversity of habitats and acquisition conditions across the studied coral lagoon. Our study offers a potential breakthrough for the monitoring of ray populations in coral reef ecosystems by providing a fast and accurate alternative to the manual processing of aerial videos. Our deep learning approach can be extended to the detection of other elasmobranchs and applied to systematic aerial surveys to not only detect individuals but also estimate species density in coral reef habitats.

Highlights

► Efficient techniques are needed to monitor vulnerable elasmobranchs in space and time. ► Deep learning applied to images is a powerful tool for automated wildlife monitoring. ► Our deep learning model successfully detected 92% of eagle rays on images. ► This study is a step forward for ray monitoring in coral reef ecosystems.

, Vigliola L.

Introduction

Elasmobranchs, a subclass of cartilaginous fishes composed of sharks, rays, skates and sawfish, are among the most endangered animal taxa in the oceans [START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF]). These species are intrinsically sensitive to human activities due to their slow growth rate and limited reproduction capacity, preventing them from quickly recovering from overexploitation [START_REF] Pacoureau | Half a century of global decline in oceanic sharks and rays[END_REF].

Elasmobranchs are primarily threatened by targeted fisheries and incidental catches, although habitat degradation is a growing threat for coastal species [START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF][START_REF] Yan | Overfishing and habitat loss drive range contraction of iconic marine fishes to near extinction[END_REF].

Within the 1,199 species of elasmobranchs assessed by the IUCN in 2021, 10.4% were listed as near-threatened, 15% as vulnerable, 10.1% as endangered (compared to 4.1% in 2010), 7.5% as critically endangered, and 12.9% as data deficient [START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF]. Rays are even more threatened than sharks with 36% of all species threatened compared to 31.2% [START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF].

Currently, the limited knowledge and monitoring of elasmobranch abundance and distribution is a major impediment to the implementation of targeted conservation measures [START_REF] Jabado | Troubled waters : Threats and extinction risk of the sharks, rays and chimaeras of the Arabian Sea and adjacent waters[END_REF].

To fill these knowledge gaps, new techniques are urgently needed to efficiently and rapidly monitor threatened elasmobranchs in space and time in order to identify their key habitats and provide abundance estimates at the basis of IUCN assessments.

Video-surveys from drones or light aircraft are increasingly used to assess the distribution, behavior and abundance of marine megafauna [START_REF] Hodgson | Unmanned Aerial Vehicles (UAVs) for Surveying Marine Fauna : A Dugong Case Study[END_REF][START_REF] Kelaher | Assessing variation in assemblages of large marine fauna off ocean beaches using drones[END_REF][START_REF] Schofield | Aerial and underwater surveys reveal temporal variation in cleaning-station use by sea turtles at a temperate breeding area[END_REF]. Such digital surveys are particularly suited to study sharks and rays in coral lagoons where clear and shallow waters facilitate their detection [START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF][START_REF] Rieucau | Using unmanned aerial vehicle (UAV) surveys and image analysis in the study of large surface-associated marine species : A case study on reef sharks Carcharhinus melanopterus shoaling behaviour[END_REF]. Video-surveys offer important advantages over traditional observer-based surveys by generating precise and verifiable observations that are free from observer fatigue and subjectivity [START_REF] Colefax | The potential for unmanned aerial vehicles (UAVs) to conduct marine fauna surveys in place of manned aircraft[END_REF][START_REF] Kelaher | Comparison of sampling precision for nearshore marine wildlife using unmanned and manned aerial surveys[END_REF]. However, manual video analysis is a major bottleneck for timely conservation applications, as visualizing hours of footage is both extremely time-consuming and error-prone [START_REF] Ditria | Automating the Analysis of Fish Abundance Using Object Detection : Optimizing Animal Ecology With Deep Learning[END_REF][START_REF] Norouzzadeh | Automatically identifying, counting, and describing wild animals in camera-trap images with deep learning[END_REF][START_REF] Villon | A Deep learning method for accurate and fast identification of coral reef fishes in underwater images[END_REF]. Deep learning algorithms offer great promises to overcome this limitation by allowing the automated identification and detection of species on images [START_REF] Christin | Applications for deep learning in ecology[END_REF][START_REF] Norouzzadeh | Automatically identifying, counting, and describing wild animals in camera-trap images with deep learning[END_REF][START_REF] Torney | A comparison of deep learning and citizen science techniques for counting wildlife in aerial survey images[END_REF][START_REF] Eikelboom | Improving the precision and accuracy of animal population estimates with aerial image object detection[END_REF]). Such models have been successfully applied for the detection of sea turtles [START_REF] Dujon | Machine learning to detect marine animals in UAV imagery : Effect of morphology, spacing, behaviour and habitat[END_REF][START_REF] Gray | A Convolutional Neural Network for Detecting Sea Turtles in Drone Imagery[END_REF], dugongs [START_REF] Mannocci | Leveraging social media and deep learning to detect rare megafauna in video surveys[END_REF], pinnipeds [START_REF] Dujon | Machine learning to detect marine animals in UAV imagery : Effect of morphology, spacing, behaviour and habitat[END_REF][START_REF] Padubidri | Counting sea lions and elephants from aerial photography using deep learning with density maps[END_REF] and whales [START_REF] Gray | Drones and convolutional neural networks facilitate automated and accurate cetacean species identification and photogrammetry[END_REF][START_REF] Guirado | Whale counting in satellite and aerial images with deep learning[END_REF]. Although there are a few applications for elasmobranchs, these are generally dedicated to monitoring shark risks [START_REF] Gorkin | Sharkeye : Real-Time Autonomous Personal Shark Alerting via Aerial Surveillance[END_REF] rather than conservation objectives requiring abundance and distribution estimates. Accurate deep learning models would drastically increase the efficiency of aerial monitoring for these threatened species.

In this study, we combined aerial video-surveys and deep learning to detect eagle rays and map their distribution throughout a lagoon in New Caledonia, Southwest Pacific. New Caledonia hosts exceptional coral reefs and lagoons, which form one of the three most extensive reef systems in the world [START_REF] Ceccarelli | The Coral Sea[END_REF]. Eagle rays are conspicuous rays of the Myliobatidae family that are easily spotted from the surface owing to their relatively large size and characteristic diamond shape [START_REF] Last | Descriptions of new sharks and rays from Borneo[END_REF], making them good candidates for automated detection on aerial images. Two species of Myliobatidae are present in New Caledonia, the spotted eagle ray, Aetobatus narinari which is common, and the rarer mottled eagle ray, Aetomylaeus maculatus [START_REF] Fricke | Checklist of the fishes of New Caledonia, and their distribution in the Southwest Pacific Ocean (Pisces)[END_REF]. These species have been classified as globally endangered by the IUCN in 2020 [START_REF] Dulvy | Aetobatus narinari : The IUCN Red List of Threatened Species 2021[END_REF][START_REF] Rigby | Aetomylaeus maculatus : The IUCN Red List of Threatened Species 2020[END_REF], stressing the urgent need to monitor the trends of their populations to feed global indicators like the Living Planet Index [START_REF] Pacoureau | Half a century of global decline in oceanic sharks and rays[END_REF]. We trained and evaluated a deep learning model to automatically detect eagle rays on aerial images collected from an ultra-light motor plane (ULM). We then mapped their distribution across the studied lagoon. Our study unravels the potential of deep learning applied to aerial surveys for monitoring the distribution of vulnerable elasmobranchs in coral reefs.

Material and methods

Video data collection

Video-surveys were conducted from an amphibious ULM (AirMax SeaMax) operating touristic flights over the Poé lagoon on the Western coast of New Caledonia (Supplementary Figure A). This lagoon is shallower than 5 m and characterized by shallow reef, seagrass and sandy habitats. A GoPro Hero Black 7 camera was mounted under the right wing of the ULM, pointing downward.

The camera was configured to record videos at a rate of 24 frames per second in linear field of view mode at a resolution of 2.7 K (2,704 x 1,520 pixels) with integrated image stabilization. The camera was manually triggered by the pilot before each flight. Telemetry data, including GPS coordinates and altitudes, were also recorded by the GoPro along each flight (at a rate of 8 to 12 positions per second). The mean altitude across all flights was 152 m (standard deviation SD= 52 m). At this altitude the image covered a mean surface area of 161 m × 287 m corresponding to a ground sampling distance of 11 cm per pixel. In total, over 42 hours of videos representing 36 fly days were collected from September 2019 to January 2020 in good weather conditions.

Image annotation

Image annotation is a crucial prerequisite before applying deep learning models [START_REF] Gray | A Convolutional Neural Network for Detecting Sea Turtles in Drone Imagery[END_REF][START_REF] Norouzzadeh | Automatically identifying, counting, and describing wild animals in camera-trap images with deep learning[END_REF][START_REF] Villon | A new method to control error rates in automated identification with deep learning algorithms[END_REF]. All videos were first visualized by a team of students who recorded the times at which they spotted eagle rays (and other megafauna species). Videos that contained at least one eagle ray (representing 114 videos from a total of 228 (Table 1)) were then imported into a custom online application (http://webfish.mbb.univ-montp2.fr/) (Supplementary Figure B). Next, images were extracted from all videos at a rate of one image per second, representing a compromise between image diversity and annotation time.

The annotation procedure consisted in manually drawing rectangle bounding boxes around identified eagle rays and associating labels to these individuals. Only individuals that could be identified without ambiguity as eagle rays, owing to their large size, diamond shape and dark colour, were annotated. Although Aetomylaeus maculatus is generally smaller than Aetobatus narinari, their color patterns are similar (light spots on a dark disc) so they are easily confused in situ. The presence of a long spine near the tail's base of A. narinari can help to differentiate it from A. maculatus which has a long but spineless tail (Froese and Pauly 2021). Since we could not distinguish one or the other species on aerial videos we built a generic eagle ray (i.e., Myliobatidae) detector, although most sightings were likely of A. narinari which is much more common in New Caledonia [START_REF] Fricke | Checklist of the fishes of New Caledonia, and their distribution in the Southwest Pacific Ocean (Pisces)[END_REF]. Each annotation yielded a text file containing the coordinates and label of the bounding box, along with the corresponding image in jpeg format (examples of images are provided in Supplementary Figure C).

Eagle ray detection model

We used a convolutional neural network (CNN), a class of deep learning models that is widely applied for image classification and object detection, i.e., the task of simultaneously localizing and classifying objects on images [START_REF] Lecun | Deep learning[END_REF]. CNNs represent by far the most commonly used category of deep learning models in ecology [START_REF] Christin | Applications for deep learning in ecology[END_REF].

They are formed by stacked groups of convolutional layers and pooling layers that are particularly suited to process image inputs. Convolutional layers extract local combinations of pixels known as 'features' from images. In the convolution operation, a filter defined by a set of weights computes the local weighted sum of pixels over the three colour channels of a given image [START_REF] Lecun | Deep learning[END_REF]. In practice, CNNs are fed with large amounts of images in which target objects have been manually annotated so they can be trained to associate labels to a given object. During this training phase, the weights are iteratively modified to obtain the desired answer by minimizing the error function between the output of the CNN and the correct answer through a process called backpropagation [START_REF] Lecun | Deep learning[END_REF]. The final output of the CNN is a confidence score for each of the learned objects.

We selected a Faster R-CNN network [START_REF] Ren | Faster R-CNN : Towards Real-Time Object Detection with Region Proposal Networks[END_REF] publicly available from the Tensorflow model zoo and tuned it for eagle ray detection on aerial images. The Faster R-CNN is a deep learning algorithm specialized for object detection that consists of two fully-convolutional networks:

(1) a region proposal network, which predicts object positions along with their 'objectness' scores and (2) a detection network, which extracts features from the proposed regions and provides class labels for the bounding boxes. We specifically used a Faster-RCNN with a ResNet-101 backbone, a deep architecture in which layers have been reformulated as residual functions of input layers, leading to better optimization and increased accuracy. Our eagle ray detection framework followed the three main steps detailed below: 1) Image pre-processing, 2) Model training and 3) Model accuracy assessment. The eagle ray detection framework is illustrated in Figure 1.

Image pre-processing

A total of 314 ULM images containing at least one eagle ray (representing 372 individual encounters) were extracted out of the 79,325 collected images (Table 1). Bounding boxes surrounding eagle rays spanned on average 25 x 25 pixels (pi) on the 2,704 x 1,520 pi images, corresponding to a ratio of 0.0002 between the bounding box area and the image area. To maximize the detection of small eagle rays on ULM images, we split each image into four images with half the original size (i.e., 1,352 x 760 pi). This yielded 308 images with eagle rays (353 individual encounters), as rays located across image boundaries were lost. Image splitting approaches are known to efficiently boost detection accuracy by increasing the relative pixel area of small objects with respect to the entire images, thereby limiting detail losses when images are processed throughout the network [START_REF] Unel | The Power of Tiling for Small Object Detection[END_REF] .

Next, images were randomly partitioned, using 80% of images for the training (and validation) subset (corresponding to approximately 250 images) and 20% of images (approx. 60 images) for the testing subset. Full independence between subsets was ensured by selecting images belonging to different videos between the subsets. The training subset was then artificially augmented by applying random transformations to images, including rotations (by -10 to +10 degrees), translations (by -10 to +10 %), scaling (by 80 to 120%), horizontal and vertical flipping, and contrast modification (i.e., multiplying all image pixels with a value ranging from 0.6 to 1.4).

Artificial data augmentation is a particularly efficient technique for improving the generalization performance and accuracy of object detection models [START_REF] Zoph | Learning Data Augmentation Strategies for Object Detection[END_REF].

Model training

We initialized our Faster R-CNN with pre-trained weights based on the COCO (Common Objects in Context) dataset [START_REF] Lin | Microsoft COCO : Common Objects in Context[END_REF] downloaded from the Tensorflow model zoo. This process of applying previously learned knowledge to solve a new problem, called transfer learning, improves model accuracy and generalization when a limited annotated dataset is available [START_REF] Chen | End-to-End Airplane Detection Using Transfer Learning in Remote Sensing Images[END_REF]. We then trained the Faster R-CNN using a stochastic gradient descent optimizer with a momentum of 0.9 for the loss function [START_REF] Qian | On the momentum term in gradient descent learning algorithms[END_REF]. We applied a learning rate of 10 -3 , a L2 regularization (with a lambda of 0.004), and a dropout of 50% to mitigate overfitting [START_REF] Srivastava | Dropout : A Simple Way to Prevent Neural Networks from Overfitting[END_REF]. The training was stopped after 50,000 iterations to prevent overfitting as indicated by an increasing loss function for the validation subset [START_REF] Sarle | Stopped training and other remedies for overfitting[END_REF].

Model accuracy assessment

The Faster R-CNN was then applied for eagle ray detection on the test subset and its accuracy was evaluated using a 5-fold cross-validation. K-fold cross-validation is a common procedure for evaluating machine learning models while preventing systematic biases due to the partitioning of data subsets [START_REF] Wong | Performance evaluation of classification algorithms by k-fold and leave-oneout cross validation[END_REF]. The initialized model was trained five times, each time with a different training subset and its accuracy was evaluated five times, each time on an independent test subset.

We applied lenient thresholds of 50% for both the confidence score of predictions and the overlap of predictions with observations, since minimizing false negatives is more crucial than avoiding false positives in the case of rare megafauna species [START_REF] Villon | A new method to control error rates in automated identification with deep learning algorithms[END_REF]. As such, a predicted bounding box that was associated with a confidence score of at least 50% and that overlapped at least 50% in surface with an annotated eagle ray was considered a true positive (TP). Predicted bounding boxes not corresponding to an annotated bounding box were false positives (FP), while annotated bounding boxes not corresponding to a predicted bounding box were false negatives (FN). For each cross-validation test subset, the number of TPs, FPs and FNs were computed and performance metrics were calculated as described below.

Precision is the percentage of TPs with respect to all predictions (Equation ( 1)). It represents the percentage of predictions that are correct (the closest to 1, the fewest false positives):

Precision= TP / (TP + FP)

(1)

Recall (or sensitivity) is the percentage of TPs with respect to all annotated objects (Equation ( 2)).

It represents the percentage of positives that are actually predicted (the closest to 1, the fewest false negatives):

Recall= TP / (TP + FN) (2)
Finally, the f1-score evaluates the balance between FPs and FNs. It is an overall measure of accuracy calculated as the harmonic mean of precision and recall (Equation ( 3)).

F1-score= 2 x Recall x Precision / (Recall + Precision) (3)

Finally, the mean and standard deviation of the performance metrics were computed across the 5fold cross-validations splits.

We used the open-source Tensorflow object detection API version 1 [START_REF] Abadi | TensorFlow : Large-Scale Machine Learning on Heterogeneous Distributed Systems[END_REF] in Python version 3 for the training and testing of our model. One training process lasted on average 3 hours on a NVIDIA Quadro P6000 GPU with 64 GB of RAM. The application of the model took on average 5 seconds per image.

Spatial distribution of eagle rays

Locations of eagle ray occurrences obtained from both manual annotation and the deep learning model were mapped in the study area by retrieving the GPS coordinates of their image identifiers.

Locations of all ULM tracks were also mapped by retrieving the GPS coordinates of all video images. To account for the heterogeneous sampling effort, the encounter rate (individuals/km) was mapped throughout the study area. To do so, we created a spatial grid of 0.005° longitude x 0.005° latitude and summed the number of eagle rays and the length of ULM tracks in each cell. The number of individuals was then divided by the length of ULM tracks per cell to obtain the encounter rate. All maps were produced in R (version 4.0.3) with the OpenStreetMap [START_REF] Fellows | OpenStreetMap : Access to Open Street Map Raster Images[END_REF]) and ggplot2 packages [START_REF] Wickham | ggplot2 : Create Elegant Data Visualisations Using the Grammar of Graphics[END_REF].

Results

Deep learning model accuracy

Our deep learning model trained with 255 images on average (range= 252 -259 between crossvalidations) accurately detected eagle rays on independent images from the same lagoon. The model reached a mean precision of 0.90 on test images (SD= 0.08), meaning that 90% of the model predictions corresponded to a manually annotated eagle ray (i.e., were TPs) (Figure 2).

False positives were primarily associated with coral patches. The mean recall was 0.92 (SD= 0.06), meaning that 92% of the annotated eagle rays were detected (Figure 2). The model successfully detected eagle rays in various contexts, as illustrated in Figure 2 and Supplementary

Spatial distribution of eagle rays

Eagle rays detected from the deep learning model were distributed throughout the study area, but appeared concentrated in a more intensively surveyed portion of the barrier reef near the easternmost channel (Figure 4-a). The few FPs and FNs were scattered across the lagoon and on the barrier reef (Supplementary Figure E-1). The encounter rate map, accounting for the heterogeneous sampling effort, confirmed the slightly higher occurrence of eagle rays on the barrier reef compared to the lagoon (Figure 4-b). The spatial distributions of detected eagle rays and their encounter rates were similar to that of all annotated eagle rays (Supplementary Figures E-2 and E-3).

Discussion

More than one third of all cartilaginous fishes are threatened with extinction, primarily due to overfishing [START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF]. Rays are no exception as they represent 56.3% of threatened chondrichthyan and 12.3% of ray species are still lacking sufficient data for assessment (Dulvy et al. 2021). As human activities continue to jeopardize these species [START_REF] Pacoureau | Half a century of global decline in oceanic sharks and rays[END_REF][START_REF] Yan | Overfishing and habitat loss drive range contraction of iconic marine fishes to near extinction[END_REF], there is an urgent need for reliable and efficient approaches for monitoring populations.

Our study revealed the potential of deep learning for the accurate detection of eagle rays in coral reef ecosystems. Our model trained with fewer than 260 aerial images was able to detect 92% of the eagle rays on independent images from the same lagoon. Our study paves the way towards automated ray population monitoring in coral reefs by providing a fast and accurate alternative to the manual processing of aerial images [START_REF] Kelaher | Assessing variation in assemblages of large marine fauna off ocean beaches using drones[END_REF][START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF]. While deep learning for elasmobranch aerial detection has been applied in the context of beach surveillance [START_REF] Gorkin | Sharkeye : Real-Time Autonomous Personal Shark Alerting via Aerial Surveillance[END_REF], we present its first implementation towards ecological and conservation applications, including species distribution mapping).

Eagle ray detection accuracy

Our model achieved a very good detection performance despite the modest size of the training dataset. Obtaining large amounts of images for training deep learning models is a major bottleneck for ecological and conservation applications [START_REF] Christin | Applications for deep learning in ecology[END_REF]. To overcome this limitation, we relied on transfer learning and artificial data augmentation, two efficient techniques that are widely used for training models in data-limited situations [START_REF] Schneider | Three critical factors affecting automated image species recognition performance for camera traps[END_REF]. The model was successful at avoiding missed occurrences (false negatives), which is most critical when the objective is to detect vulnerable species that occur in low numbers such as rays and sharks [START_REF] Villon | A new method to control error rates in automated identification with deep learning algorithms[END_REF]. Eagle rays were consistently detected across the diversity of habitats (e.g., soft bottom and barrier reef) and acquisition conditions (e.g., luminosity, altitude and camera angle) in our study area. The robustness of the model at detecting eagle rays in more contexts and its generalization to new data could be further increased by expanding the size of both the training and the test datasets and the contextual variety at new sites in New-Caledonia and beyond. Moreover, there is a need to test the model's generalizability to a larger dataset in the future, as the size of the test dataset is also limited. The model was equally successful at avoiding false positives, with few misdetections primarily associated with coral patches. To eliminate these false positives, coral patch annotations could be incorporated into the training dataset so that the model explicitly learns this class. As deep learning algorithms rapidly improve, we could further enhance our eagle ray detection method by using most up-to-date object detection CNNs such as the YOLOv3 that achieved a high performance on fish detection [START_REF] Jalal | Fish detection and species classification in underwater environments using deep learning with temporal information[END_REF].

Comparison with other monitoring methods

Effective conservation requires up-to-date and high quality data collected with limited monetary and human costs over repeated periods [START_REF] Fust | Development perspectives for the application of autonomous, unmanned aerial systems (UASs) in wildlife conservation[END_REF]. Previous studies on the distribution and movements of eagle rays have relied on acoustic [START_REF] Degroot | Fine-scale movement and habitat use of whitespotted eagle rays Aetobatus narinari in the Indian River Lagoon, Florida, USA[END_REF] and satellite telemetry [START_REF] Ajemian | Towed-float satellite telemetry tracks large-scale movement and habitat connectivity of myliobatid stingrays[END_REF]. Active acoustic telemetry implies following the individuals in order to determine their movements in the water column, but is generally restricted to few individuals and necessitates a large array of hydrophones [START_REF] Degroot | Fine-scale movement and habitat use of whitespotted eagle rays Aetobatus narinari in the Indian River Lagoon, Florida, USA[END_REF]. Satellite telemetry allows tracking rays over potentially large spatial scales, but is constrained by the frequency and precision of GPS data and associated costs [START_REF] Ajemian | Towed-float satellite telemetry tracks large-scale movement and habitat connectivity of myliobatid stingrays[END_REF]. Both methods are intrusive as they require catching and manipulating individuals to attach the tags properly. Surveys from scuba divers and baited remote underwater videos [START_REF] Rizzari | Diversity, abundance, and distribution of reef sharks on outer-shelf reefs of the Great Barrier Reef, Australia[END_REF][START_REF] Ward-Paige | Global evaluation of shark sanctuaries[END_REF] are also used for elasmobranch censuses, especially for species that live further from the surface . However, these underwater surveys are limited in their spatial extent and may fail to detect the most elusive species [START_REF] Juhel | Reef accessibility impairs the protection of sharks[END_REF]. Moreover, observations may not be precisely located and are not verifiable, unlike those derived from video footage.

In this study, aerial images collected from an off-the-shelf camera and processed with a deep learning algorithm allowed us to precisely locate eagle rays in a coral lagoon at low financial and operational costs. The opportunistic use of an aircraft dedicated to touristic flights led to an heterogeneous survey effort, preventing the estimation of abundance from the traditional strip transect methodology [START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF][START_REF] Sykora-Bodie | Quantifying Nearshore Sea Turtle Densities : Applications of Unmanned Aerial Systems for Population Assessments[END_REF]. Nevertheless, our accurate algorithm will be applicable to images collected along systematically-designed transects for abundance estimation in the future. Despite the heterogeneous survey effort, the current method suggests a widespread distribution of eagle rays across a variety of coral reef habitats, which is in accordance with previous study [START_REF] Ajemian | Towed-float satellite telemetry tracks large-scale movement and habitat connectivity of myliobatid stingrays[END_REF]. Future studies should seek to quantify habitat preferences of eagle rays by linking effort-corrected encounter rates to local habitat information [START_REF] Ajemian | Estimating the Potential Impacts of Large Mesopredators on Benthic Resources : Integrative Assessment of Spotted Eagle Ray Foraging Ecology in Bermuda[END_REF][START_REF] Degroot | Fine-scale movement and habitat use of whitespotted eagle rays Aetobatus narinari in the Indian River Lagoon, Florida, USA[END_REF].

Shark and ray monitoring requires detection and census methods that are adapted to the studied habitats. While aerial surveys are very efficient in coral reefs with clear and shallow waters, open water or turbid waters (e.g., estuaries, mangroves) require non-visual methods such as acoustic telemetry. Environmental DNA is also an innovative method at the species level that can be notably used to detect rare species including elasmobranchs [START_REF] Boussarie | Environmental DNA illuminates the dark diversity of sharks[END_REF]. Coral lagoons are major habitats for eagle rays [START_REF] Degroot | Fine-scale movement and habitat use of whitespotted eagle rays Aetobatus narinari in the Indian River Lagoon, Florida, USA[END_REF][START_REF] Ajemian | Estimating the Potential Impacts of Large Mesopredators on Benthic Resources : Integrative Assessment of Spotted Eagle Ray Foraging Ecology in Bermuda[END_REF] and our aerial approach proved efficient for monitoring populations in these habitats. Our approach can be complemented by other methods (e.g., eDNA, acoustic telemetry) in habitats where eagle rays can occur [START_REF] Ajemian | Towed-float satellite telemetry tracks large-scale movement and habitat connectivity of myliobatid stingrays[END_REF][START_REF] Sellas | Population Structure and Seasonal Migration of the Spotted Eagle Ray, Aetobatus narinari[END_REF] but waters are deep or turbid.

Implications for elasmobranchs monitoring in coral reefs

Data on population trends and distributions of rays and sharks are difficult to collect; yet, such information is critical to establish appropriate conservation and management actions [START_REF] Dwyer | Individual and Population Benefits of Marine Reserves for Reef Sharks[END_REF][START_REF] Macneil | Global status and conservation potential of reef sharks[END_REF][START_REF] Pacoureau | Half a century of global decline in oceanic sharks and rays[END_REF]. Our approach combining video-surveys and Finally, our approach could be extended to systematic video-surveys from manned aircraft or drones to not only detect individuals, but also count them to derive abundance estimates and species density maps in a study area. Drones make a viable alternative to manned aircraft for marine megafauna surveys [START_REF] Gray | A Convolutional Neural Network for Detecting Sea Turtles in Drone Imagery[END_REF][START_REF] Hodgson | Unmanned Aerial Vehicles (UAVs) for Surveying Marine Fauna : A Dugong Case Study[END_REF][START_REF] Kelaher | Assessing variation in assemblages of large marine fauna off ocean beaches using drones[END_REF][START_REF] Kiszka | Using unmanned aerial vehicles (UAVs) to investigate shark and ray densities in a shallow coral lagoon[END_REF], alleviating safety risks, monetary costs and carbon emissions [START_REF] Hodgson | Unmanned Aerial Vehicles (UAVs) for Surveying Marine Fauna : A Dugong Case Study[END_REF].

However, the use of drones is subject to strict airspace regulations, and legislation in many areas necessitates the pilot to maintain visual-line-of-sight with the drone [START_REF] Raoult | Operational Protocols for the Use of Drones in Marine Animal Research[END_REF]. The platform choice will ultimately depend on the study question and the required imagery characteristics. Using an aircraft dedicated to touristic flights allowed us to achieve greater spatial and temporal coverage than would have been possible with a single drone and with no need to acquire permits. This method could be implemented in other touristic locations (e.g., Australia, French Polynesia, the Caribbean) where local companies operate scenic, low altitude flights over coastal areas.

Overall, our cost-effective approach succeeded in collecting high-quality images for training a deep learning model able to detect 92% of eagle rays in coral reefs. This new eagle ray detector will be critical for deriving abundance estimates in order to closely monitor these vulnerable populations in the future.

Tables Table 1: Overview of the New Caledonia video database. Abbreviations: SD= standard deviation, pi= pixels.

Number of videos

Mean video duration

Total video duration

  The barrier reef includes three deep passes and channels reaching 30 m. Part of the Poé lagoon was declared as a natural reserve (IUCN category IV) in 2006 and is located within the broader South Province Park created in 2009 and the UNESCO World Heritage area established in 2008.

Figure D .

 D Figure D. The mean f1-score balancing FPs and FNs was 0.91 (SD= 0.06). Precision, recall and the f1-score showed little sensitivity to the prediction confidence score (Figure 3).

  deep learning offers a potential breakthrough for the automated monitoring of eagle rays in coral reef ecosystems. The ability of our model at detecting eagle rays in the variety of habitats and conditions encompassed by our data highlights its potential robustness in a broad range of contexts. Future work should assess the model transferability to other coral lagoons in New Caledonia and beyond. Robust detection models would be particularly beneficial for ray monitoring in the Indo-Pacific biodiversity triangle where conservation efforts are most urgent due to pronounced levels of human threats[START_REF] Dulvy | Overfishing drives over one-third of all sharks and rays toward a global extinction crisis[END_REF]. Our deep learning approach could be applied to other distinctive elasmobranchs, provided sufficient images of these species are available for training the model.

FiguresFigure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: Results of eagle ray detection on test images for a prediction confidence score of 50%. The left graph shows the mean percentage of true positives (TPs) and false positives (FPs) with respect to all predictions. The right graph shows the mean percentage of TPs and false negatives (FNs) in the observations. The error bars are the standard deviations from the means. Examples are provided below the graphs for a TP in green (prediction associated with an annotation shown in white), a FP in red (prediction not corresponding to an annotation; here a coral patch) and a FN (annotation not corresponding to a prediction). Further examples of detection results are provided in Appendix D.

Figure 3 :

 3 Figure 3: Mean precision, recall and f1-score on the test images for varying prediction confidence scores. The standard deviation is represented by the shaded area.

Figure 4 :

 4 Figure 4: Spatial distribution of (a) eagle ray detections (dots) from the deep learning model mapped by retrieving the GPS coordinates of their image identifiers and the corresponding ULM flight tracks (black lines) and (b) the encounter rate (individuals/km) of detected eagle rays calculated on a spatial grid of 0.005° longitude x 0.005° latitude (the calculation is detailed in section 2.7).

Acknowledgements

We are grateful to Jugurtha Ifticen, Lola Romant, Nacim Guellati, Gwendal Quimbre and Matis Toitot for their help with visualisation of ULM videos and annotation of images. We thank Air Paradise for their collaboration in collecting ULM video sequences in New Caledonia. This project received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 845178 ('MEGAFAUNA').

Data statement

The eagle ray image database will be made available on Zenodo.

Author contributions

LM, DM, LV and MC conceived the ideas and designed the methodology; LM, LV and DM collected the data; LM and LD analysed the data and developed scripts; LD and LM led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.


Laurent Bulteau
email: laurent.bulteau@ens.fr

Guillaume Fertin

Irena Rusu
email: irena.rusu@univ-nantes.fr

Maximal Strip Recovery Problem with Gaps: Hardness and Approximation Algorithms

Keywords: algorithmic complexity, approximation algorithms, comparative maps, genome comparison, synteny blocks

Given two comparative maps, that is two sequences of markers each representing a genome, the Maximal Strip Recovery problem (MSR) asks to extract a largest sequence of markers from each map such that the two extracted sequences are decomposable into non-overlapping strips (or synteny blocks). This aims at dening a robust set of synteny blocks between dierent species, which is a key to understand the evolution process since their last common ancestor. In this paper, we add a fundamental constraint to the initial problem, which expresses the biologically sustained need to bound the number of intermediate (non-selected) markers between two consecutive markers in a strip. We therefore introduce the problem δ-gap-MSR, where δ is a (usually small) non-negative integer that upper bounds the number of non-selected markers between two consecutive markers in a strip. Depending on the nature of the comparative maps (i.e., with or without duplicates), we show that δ-gap-MSR is NP-complete for any δ ≥ 1, and even APX-hard for any δ ≥ 2. We also provide two approximation algorithms, with ratio 1.8 for δ = 1, and ratio 4 for δ ≥ 2.

Introduction

In comparative genomics, nding synteny blocks (that is, regions with similar content and gene order) of two genomes is a crucial task, as the decomposition of genomes into synteny blocks allows to estimate the nature of genome rearrangement events that hold during the evolution process since the last common ancestor of the genomes.

In addition to the diculty to dene a synteny block precisely, another difculty is introduced by the quality of genome annotation. Zheng et al. [START_REF] Zheng | Removing noise and ambiguities from comparative maps in rearrangement analysis[END_REF] make a list of possible errors and ambiguities introduced by the mapping technology, which is used to obtain a representation of a genome as a sequence of markers, called a genomic map. Each marker represents a small, specic element which has been identied on the genome, at a specic position which is the marker's position. Comparing two genomes is then possible using their genomic maps, assuming that the pairs of identical markers on the two genomes are known (the maps are then called comparative maps). Comparative maps are less precise than genome sequences (either as DNA sequences or as sequences of genes), but still allow the identication of synteny blocks.

The problem that needs to be solved when no error occurs is the following: Given two comparative maps, decompose them into non-intersecting synteny blocks. In case of errors or ambiguities, Zheng et al. [START_REF] Zheng | Removing noise and ambiguities from comparative maps in rearrangement analysis[END_REF] propose to switch to the following problem: Given two comparative maps, nd a longest (possibly non-contiguous) subsequence of markers in each comparative map, such that the subsequences are decomposable into non-intersecting synteny blocks. The idea behind this maximization problem is that true synteny is possibly interrupted by erroneous or ambiguous markers, which should be discarded before searching for synteny blocks.

The problem, called Maximal Strip Recovery (MSR), is obtained from this maximization problem using comparative maps with signed, but not duplicated, markers, and a specic denition of synteny blocks. Synteny blocks are dened as strips, which are contiguous sequences of at least two markers that occur on each genome either in the same order, or in reverse order and with a reversed sign. Zheng et al [START_REF] Zheng | Removing noise and ambiguities from comparative maps in rearrangement analysis[END_REF] and Choi et al. [START_REF] Choi | Algorithms for the extraction of synteny blocks from comparative maps[END_REF] propose two heuristics to solve the MSR problem. Chen et al. [START_REF] Chen | On recovering syntenic blocks from comparative maps[END_REF] devise a 4-approximation algorithm for it, propose its extension, called MSR-d, to an arbitrary number d ≥ 2 of genomes and show that MSR-3 is NP-complete. The NP-completeness of MSR (or equivalently MSR-2) is a result obtained by Wang et al. [START_REF] Wang | On the tractability of maximal strip recovery[END_REF], who also propose FPT algorithms for MSR-d (with arbitrary d) and MSR-DU, the variant of MSR where duplicated markers are allowed in the maps and in dierent synteny blocks.

The MSR problem takes into account the need to keep as much of the data as possible from the initial comparative maps and the need to have conict-free synteny blocks. However, it is too permissive as it allows two consecutive elements from one strip to be separated by an arbitrary long gap (in terms of intermediate markers) on the initial comparative maps, and possibly to be very close on one map and very far from each other on the other. As the discarded elements are supposed to be errors and ambiguities (which are rather the exception than the rule), and the elements kept in the subsequences are supposed to be the safe information (which is the major part of the comparative information), it follows that a safe synteny block should not allow arbitrarily long gaps.

We therefore introduce and study in this paper the δ-gap-MSR problem, a restriction of the MSR problem where the allowed gaps along the comparative maps between two consecutive elements in a strip are upper bounded by parameter δ, where δ is a given (usually small) non-negative integer. We investigate the algorithmic complexity of δ-gap-MSR depending on the allowed multiplicity for a marker and prove the results given in Table 1. For the NP-complete or APX-hard cases, we provide two approximation algorithms, whose approximation ratios are given in Table 2. The organization of the paper is as follows. In Section 2, we introduce some notations, and we dene formally MSR, MSR-DU, δ-gap-MSR and δ-gap-MSR-DU. We prove in Section 3 the hardness results (NP-completeness for δ = 1 in Section 3.1, APX-completeness for δ ≥ 2 in Section 3.2). We then give approximation algorithms in Section 4: a 1.8-approximation for 1-gap-MSR in Section 4.1, and a general 4-approximation in Section 4.2. Due to space constraints, most of the proofs are omitted from this paper.

Notations and Denitions

A comparative map M is a sequence of signed integers, where the absolute value of each integer represents a specic marker, and the sign represents the orientation of the marker on the chromosome. A marker may appear several times in a comparative map, possibly with dierent orientations: in this case, we say that the comparative map M has duplicates (the presence of duplicates is useful if we do not want to distinguish paralogs in the comparative map). A sequence M is denoted M= m 1 , m 2 , . . . , m l , and its i th element m i is (also) denoted M [i].

A subsequence σ of M is a sequence σ 1 , . . . , σ h of markers from M with h ≥ 2 and positions i 1 < i 2 < . . . < i h respectively on M. The vector

(i 1 , . . . , i h ) is denoted idx(σ, M). The gap of σ in M is max{i k+1 -i k -1 : 1 ≤ k < h}, its length |σ| is h.
Two subsequences σ and τ are non-overlapping in M if one appears strictly before the other (i.e., if the last element of idx(σ, M) is strictly smaller than the rst element of idx(τ, M) or vice-versa). The reversed opposite

of σ 1 , . . . , σ h is -σ h , -σ h-1 , . . . , -σ 1 .
Given two comparative maps M 1 and M 2 , a prestrip is a subsequence σ of M 1 such that either σ or its reversed opposite is a subsequence of M 2 , and such that the markers in σ are pairwise dierent. The gap of a prestrip is the maximum of the gaps of the two corresponding subsequences in M 1 and M 2 . Two prestrips are non-overlapping if the corresponding subsequences are non-overlapping, both in M 1 and M 2 . A strip is a prestrip with gap 0. Strips represent synteny blocks between two comparative maps. A prestrip can also be seen as a synteny bock, but only if we consider that there is noise in the comparative maps (false markers appear between two consecutive markers of the true synteny block). A set of prestrips S is said to be feasible if it contains pairwise non-overlapping prestrips, and we write ||S|| for its total size: ||S|| = σ∈S |σ|.

We nally dene some notions of graph theory: a graph G = (V, E) is cubic if every vertex u ∈ V has degree exactly 3. A set X ⊂ V is said to be independent if for every edge (u, v) ∈ E, u / ∈ X or v / ∈ X. The cardinality of a maximum independent set of G is written α(G).

The problems MSR (for Maximal Strip Recovery, see [START_REF] Zheng | Removing noise and ambiguities from comparative maps in rearrangement analysis[END_REF]) and MSR-DU [START_REF] Chen | On recovering syntenic blocks from comparative maps[END_REF] The idea behind both those problems is that, if we nd a set of compatible prestrips with maximum total size, the elements appearing in no prestrip are considered as noise: we can remove them to clean the data. Indeed, once those elements are removed, the comparative maps can be partitioned into common strips, i.e. we have decomposed both genomes into synteny blocks with the same set of blocks in both genomes. Heuristics for the rst problem have been given in [START_REF] Zheng | Removing noise and ambiguities from comparative maps in rearrangement analysis[END_REF][START_REF] Choi | Algorithms for the extraction of synteny blocks from comparative maps[END_REF]. They have been improved in [START_REF] Chen | On recovering syntenic blocks from comparative maps[END_REF] into a 4-approximation algorithm. Finally, those problems have been proved NP-complete in [START_REF] Wang | On the tractability of maximal strip recovery[END_REF], where an FPT algorithm is also provided.

The variant we introduce, δ-gap-MSR, takes into account the fact that it is unlikely that long sequences of markers can appear only from noise and errors. If a large number of elements are inserted between two consecutive elements of a prestrip (thus, if it has a large gap), then they are not errors, and the prestrip cannot be considered a synteny block of the original genomes. Thus we dene the following two problems: Problem: δ-gap-MSR Input: Two comparative maps M 1 and M 2 without duplicates, k ∈ N. Question: Is there a feasible set S of prestrips of M 1 and M 2 , such that every σ ∈ S has gap at most δ, and ||S|| ≥ k ? Problem: δ-gap-MSR-DU Input: Two comparative maps M 1 and M 2 (possibly with duplicates), k ∈ N.

Table 1 .

 1 Complexity of variants of MSR.

	Problem	Without duplicates	With duplicates
			(-DU variant)
	0-gap-MSR	P (Section 4.2)	?
	1-gap-MSR	NP-complete (Section 3.1) NP-complete (Section 3.1)
	δ-gap-MSR (δ ≥ 2) APX-complete (Section 3.2) APX-complete (Section 3.2)
	MSR	NP-complete [8]	NP-complete [8]

Table 2 .

 2 Best approximation ratios of variants of MSR.

	Problem	Without duplicates With duplicates
			(-DU variant)
	0-gap-MSR	-	4 (Section 4.2)
	1-gap-MSR	1.8 (Section 4.1) 4 (Section 4.2)
	δ-gap-MSR (δ ≥ 2) 4 (Section 4.2)	4 (Section 4.2)
	MSR	4 [3]	4 [3]

  are dened, in their decision formulation, as follows: Problem: MSR Input: Two comparative maps M 1 and M 2 without duplicates, k ∈ N. Question: Is there a feasible set S of prestrips of M 1 and M 2 , s.t. ||S|| ≥ k ? Problem: MSR-DU Input: Two comparative maps M 1 and M 2 (possibly with duplicates), k ∈ N. Question: Is there a feasible set S of prestrips of M 1 and M 2 , s.t. ||S|| ≥ k ?

Question: Is there a feasible set S of prestrips of M 1 and M 2 , such that every σ ∈ S has gap at most δ, and ||S|| ≥ k ?

With the gap constraint we introduce, we keep only prestrips which are nearly contiguous, while tolerating some noise in the input data. Note that those problems are dened for uni-chromosomal genomes. However, algorithms can easily be adapted to handle multi-chromosomal instances.

3 Hardness Results

NP-hardness of 1-gap-MSR

In this section, we prove the following theorem.

Theorem 1 1-gap-MSR and 1-gap-MSR-DU are NP-hard.

Note that we need to consider only 1-gap-MSR (without duplicates) since NPhardness of 1-gap-MSR-DU directly follows from NP-hardness of 1-gap-MSR.

The proof uses a reduction from a variant of Maximum Independent Set, 3-colored-MIS, which is dened below. A 3-edge-coloring (also known as Tait Coloring) of a cubic graph G = (V, E) is a partition of its edges in three classes E = E A ∪ E B ∪ E C such that if two edges e 1 , e 2 ∈ E are incident to a common vertex, they belong to dierent classes.

Problem: 3-colored-MIS

Starting from any instance of 3-colored-MIS, we construct two comparative maps as follows. First, we assign a list of 4 positive integers (or 4 markers) to each vertex u ∈ V : they are denoted y A1 u , y A2 u , y B1 u and y B2 u . We also assign a list of 10 integers x 1 uv , . . . , x 10 uv to each edge (u, v) ∈ E C , in such a way that no integer appears in two dierent lists. We will also use peg markers: written with the symbol ×, they are integers appearing only once, either in M 1 or in M 2 (and thus cannot belong to any prestrip).

We construct the comparative maps with the following iterative procedure. Suppose we have arbitrarily ordered the vertices in V . In that case:

where Γ 1 and Γ 2 are dened as:

Property 3 Let G = (V, E) be an n-vertex cubic graph with a 3-edge-coloring, and let M 1 and M 2 be the two comparative maps obtained by the construction dened above. Then the optimal value of 1-gap-MSR over (M 1 , M 2 ) equals 4n + 2α(G).

Proof (of Theorem 1). The above property directly implies that our construction (which can clearly be done in polynomial time) leads to a reduction from 3colored-MIS to 1-gap-MSR, which proves Theorem 1.

δ-gap-MSR and δ-gap-MSR-DU are APX-hard

In this section, we prove the following theorem.

Theorem 4 δ-gap-MSR and δ-gap-MSR-DU are APX-hard for any δ ≥ 2.

As in the previous section, we note that we need to consider only δ-gap-MSR (without duplicates) since APX-hardness of δ-gap-MSR-DU directly follows from APX-hardness of δ-gap-MSR. For this, we use an L-reduction [START_REF] Papadimitriou | Optimization, approximation, and complexity classes[END_REF] from the variant of Maximum Independent Set restricted to cubic graphs, that we call 3-MIS here. Note that the L-reduction refers to the optimization versions of problems δ-gap-MSR and δ-gap-MSR-DU, which are easy to deduce from the decision versions presented here.

Problem: 3-MIS

It is proved in [START_REF] Alimonti | Hardness of approximating problems on cubic graphs[END_REF] that 3-MIS is APX-hard. Given a cubic graph G = (V, E), our reduction consists in constructing two comparative maps M 1 and M 2 , having properties P1, P2 and P3 described below, where Ω denotes the set of all prestrips of M 1 and M 2 having gap at most δ:

P1. There exists a bijection Φ between V and Ω P2. Every prestrip in Ω has length 2 P3. Two prestrips σ 1 and σ 2 of Ω are overlapping i Φ -1 (σ 1 ), Φ -1 (σ 2 ) ∈ E Let P k denote the path graph with k vertices.

Lemma 5 Given a cubic graph G = (V, E), one can compute in polynomial time a partition of E into two classes E B and E W (for Black and White edges), such that [START_REF] Alimonti | Hardness of approximating problems on cubic graphs[END_REF]

The rst step of the reduction is to compute a partition of E into two classes E B and E W according to Lemma 5. We then construct two comparative maps M 1 and M 2 , satisfying properties P1, P2 and P3. Moreover, incompatibilities in M 1 (resp. M 2 ) will correspond to black (resp. white) edges. We begin by assigning a dierent pair of integers (x a , x a ) to every vertex a ∈ V (G); we write Φ(a) = x a , x a .

Then, for every black component B i of order k, let V (B i ) = {a h : 1 ≤ h ≤ k} and let (a h , a h+1 ) ∈ E B for 1 ≤ h < k; we construct the following sequence:

where × l represents l consecutive peg markers. The full comparative map M 1 is given by

For M 2 , we use a similar construction, but we need to take the reversed opposite of some subsequences to avoid creating undesired prestrips. For a white component W j having 4 vertices, say a, b, c and d with (a, b), (b, c), (c, d) ∈ E W , we create the following sequence:

If W j is of order three (resp. two), we remove the extra elements from J j , i.e.

Lemma 6

The set Ω of the prestrips of M 1 and M 2 with gap less than or

The consequence of this lemma is that M 1 and M 2 satisfy the three properties P1, P2 and P3 dened above. The reduction we have described is an L-reduction from 3-MIS to δ-gap-MSR: indeed, Φ transforms an independent set of size l into a feasible set of prestrips with gap δ of total size 2l, and Φ -1 does the reverse operation. So δ-gap-MSR, like 3-MIS, is APX-hard for δ ≥ 2.

Approximation Algorithms 4.1 1.8-approximation for 1-gap-MSR

In this section, we present an approximation algorithm for 1-gap-MSR. Our result is the following.

Theorem 7 There exists a factor-1.8 approximation algorithm for 1-gap-MSR. Proof. Our algorithm makes uses of an exact algorithm to solve Maximum Weight Independent Set (MWIS) on claw-free graphs. A claw is the 4-vertex graph (V, E) with V = {a, b, c, d} and E = {(a, b), (a, c), (a, d)}. A graph is said to be claw-free if none of its induced subgraphs is isomorphic to a claw. The variant of MWIS on claw-free graphs, Claw-Free-MWIS (which is known to be in P, [START_REF] Minty | On maximal independent sets of vertices in claw-free graphs[END_REF]), is stated as follows:

Problem: Claw-Free-MWIS Input: A claw-free graph G = (V, E), a weight function w : V → R + , k ∈ R + Question: Is there an independent set X of G such that x∈X w(x) ≥ k ? Our 1.8-approximation algorithm (given in Algorithm 1) works as follows. Given two comparative maps M 1 and M 2 , compute the set Ω of all prestrips with length 2 or 3 (and gap at most 1). Longer prestrips are ignored, since they can be split into smaller ones appearing in Ω. Select a subset V λ ⊆ Ω (according to some parameter λ: see the selection process described below), and create E λ , the set of all overlapping pairs of prestrips of V λ . The pair (V λ , E λ ) forms a graph which is claw-free (see Lemma 8). An independent set for this graph (computable in polynomial time) yields a feasible set of prestrips V λ Ind . The selection of V λ amongst Ω is done as follows: given a prestrip σ of M 1 and M 2 , take the values of idx(σ, M 2 ) -λ modulo 9. This is done by the arithmetic function π 9 , which takes the values of a list modulo 9: for example, if σ has indices (30, 32, 33) in M 2 , and λ = 5, then idx(σ, M 2 ) -λ = (25, 27, 28), and π 9 (idx(σ, M 2 ) -λ) = (7, 9, 1). If the result of π 9 (idx(σ, M 2 ) -λ) belongs to some list (the list T in Algorithm 1), add σ to V λ . We only need to test the 9 dierent values of λ to obtain 9 dierent feasible sets of prestrips.

Finally, Lemma 9 proves that there exists some λ for which the total size of the corresponding V λ Ind is at least 5/9 th of a maximum feasible set of prestrips of M 1 and M 2 . Thus, Algorithm 1 is a polynomial-time algorithm giving a 1.8approximation to 1-gap-MSR, and Theorem 7 is proved.

Algorithm 1 A factor-1.8 approximation algorithm for 1-gap-MSR Input: Two comparative maps M1, M2 without duplicates.

}; Ω ← set of all prestrips of M 1 and M 2 of length 2 or 3, with gap at most 1;

Ind ←Maximum Weight Independent Set of (V λ , E λ ) with weight w;

end for return max{||V λ Ind || : 1 ≤ λ ≤ 9};

Lemma 8 For each λ, the graph

provides a solution of total size at least 5||O||/9.

Reduction to Maximum Weight Independent Set

In this section we consider the variants of Maximum Weight Independent Set on two classes of graphs: interval graphs and 2-interval graphs.

An interval graph is a graph G = (V, E), where every vertex in V is seen as an interval I of R, and such that (I, J) ∈ E i (1) I and J are distinct intervals from V , and (2)

, where every vertex in V is seen as a pair of disjoint intervals (I 1 , I 2 ) of R (also called a 2-interval), and such that ((I 1 , I 2 ), (J 1 , J 2 )) ∈ E i (1) (I 1 , I 2 ) and (J 1 , J 2 ) are distinct 2-intervals from V , and (2)

The problem Interval-MWIS is known to be polynomial [START_REF] Golumbic | Algorithmic Graph Theory and Perfect Graphs[END_REF]. On the other hand, 2-Interval-MWIS is APX-hard, and we know a 4-approximation for it [START_REF] Bar-Yehuda | Scheduling split intervals[END_REF].

Theorem 10 There exists a factor-4 approximation algorithm for δ-gap-MSR for all δ ≥ 2, and for δ-gap-MSR-DU for all δ ≥ 0.

Proof. In this proof, we describe a reduction from δ-gap-MSR to 2-Interval-MWIS. Given a pair of comparative maps and a maximal gap δ, we construct a set of 2-intervals in the following way. First, compute the set Ω of all prestrips of M 1 and M 2 having gap at most δ. Then, to each prestrip σ ∈ Ω, assign the following 2-interval (where l is |M 1 | + 1):

with the weight:

We denote G δ (M 1 , M 2 ) the weighted 2-interval graph with vertex set {I σ : σ ∈ Ω} and weight w. It has the following property: Property 11 The set {I σ : σ ∈ S} is an independent set of G δ (M 1 , M 2 ) with weight W i S is a feasible subset of Ω with total size W .

The 4-approximation algorithm for δ-gap-MSR and δ-gap-MSR-DU is the following (adapted from the 4-approximation algorithm for MSR and MSR-DU [START_REF] Chen | On recovering syntenic blocks from comparative maps[END_REF]):

Property 11 tells us that the total size of S is the weight of X, and that δ-gap-MSR-DU(M 1 , M 2 ) and 2-Interval-MWIS(G δ (M 1 , M 2 )) have the same optimal values: so S is indeed a 4-approximation of the optimal solution of δ-gap-MSR-DU(M 1 , M 2 ). We have proved Theorem 10.

Theorem 12 There exists an exact polynomial-time algorithm for 0-gap-MSR.

Proof. We consider the case where M 1 has no duplicates and the maximum gap is 0 (we only consider strips instead of prestrips): this is the case for instances of 0-gap-MSR.

We use the same reduction as for Theorem 10, with the dierence that now, G 0 (M 1 , M 2 ) is in fact an interval graph. It can be seen by considering intervals I σ = [min(idx(σ, M 1 )), max(idx(σ, M 1 ))].

We no longer need to consider the interval coming from M 2 for the following reason. If two strips overlap in M 2 , since they have gap zero, they must have a common marker m appearing in M 2 . But since m can appear only once in M 1 , they also overlap in M 1 . Thus I σ and I τ intersect i I σ and I τ intersect: G 0 (M 1 , M 2 ) can thus be seen as an interval graph. Hence, we can adapt the previous algorithm to obtain an optimal solution, and complete the proof of Theorem 12:

1. Compute the weighted interval graph G δ (M 1 , M 2 ). 2. Compute X, an optimal solution to Interval-MWIS(G δ (M 1 , M 2 )). 3. Deduce a maximal feasible set of prestrips S = {σ : I σ ∈ X}.