Mikaël Martin
email: mikael.martin-7@etudiants.ensam.eu

François Malburet
email: françois.malburet@ensam.eu

Julien Gomand
email: julien.gomand@ensam.eu

Pierre-Jean Barre
email: pierre-jean.barre@ensam.eu

Multi-Physic System Simplification Method Applied to a Helicopter Flight Axis Active Control

Keywords: Energetic Macroscopic Representation, modelling, control, Internal Model Reference, Helicopter

A helicopter flight axis control, which is a complex multi-physic system, is modelled using an energetic based graphical tool: the Energetic Macroscopic Representation. Elements of the system are mainly composed of passive technologies and their number tends to increase year after year to improve the pilots comfort by adding new functions. A new methodology is proposed to transform the system into a new active one by replacing some hydro-mechanical elements by a new controllable active mechanical source. The challenge is to simplify the flight control architecture while preserving the global behaviour of the system.

I. INTRODUCTION

From the first flight in 1903 of Wright Flyer in Kitty Hawk (US) to the last new technology ones, aircrafts have been equipped with mechanical and hydro-mechanical actuators for flight control. Over time, designers have increased the number of elements of each flight axis control to improve pilots comfort. However, mechanical and hydro-mechanical systems are relatively heavy, difficult to adjust and without possibilities of evolution. They also have limited dynamic capabilities compared to electromechanical systems but are able to provide higher forces. Following that general and inevitable evolution, helicopters became complex systems as they are composed of numerous subsystems, involving many physical domains (mechanics, hydraulics, electronics, aerodynamics…). The matter of this paper is to present a new methodology which aims to help engineers designing new architectures. More specifically, it concerns complex systems which have to be improved or simplified by substituting passive elements for active technologies. The challenge is to preserve the global behaviour of the passive system. A helicopter flight axis control is described in this paper to illustrate the methodology. Specific graphical approaches can be used to help in modelling and controlling multi-physic systems [START_REF] Hautier | The causal ordering graph -A tool for system modelling and control[END_REF]. They offer specific powerful unified modelling formalisms to simplify complex multi-physics system analysis. The methodology presented in this paper is based on the Energetic Macroscopic Representation which brings a global overview of the system. The first section is devoted to the energetic description of the helicopter flight axis control. The second section presents the methodology to transform the system composed of passive elements into an active one. Simulation results using Matlab/Simulink finally validate the proposed methodology.

II. MULTI-PHYSICS SYSTEM MODELLING USING ENERGETIC GRAPHICAL DESCRIPTION

This section is dedicated to a helicopter axis control modelling. The main elements which contribute to the functioning of helicopter rotor blades pitch control are introduced. A lumped parameter model of an independent and simplified axis control is deduced. Using an energetic graphical description, a model organisation is proposed to highlight power transfers and energy storages.

A. Studied system description

Helicopter pilots initially had to control the rotor blades pitch through four direct mechanical control chains, without any assistance. Each axis corresponded to an independent flight axis: Pitch and Roll axes are dedicated to control the cyclic rotor swachplate angles (forward and lateral directions), while the collective axis is linked to the rotorcraft power (correlated to acceleration and speed), and the Yaw axis acts as an antitorque (that of the main rotor) control through the tail rotor. With the increase of the aircraft size, and thus their mass increase, irreversible hydraulic powered assistance actuators (see, Fig. 1) became necessary to provide enough force and control rotor blades pitch. This way, pilots workload is significantly reduced: they only provide a low level of force to control the hydraulic actuators control input. This leads to one drawback: the force feedback provided to the pilots is constant whatever the rotor blades pitch are. Spring based elements (see Fig. 1) mounted in parallel with the flight controls have been necessary to provide a feedback proportional to the control stick deviation with regard to its initial position. High dynamic control movements are forbidden to preserve the stability and the integrity of the helicopter. To prevent from pilots unintentional high dynamic stick movements and to damp the pilot's actions, mechanical dampers (see Fig. 1) are also placed in parallel with the control mechanical chains. The residual forces due to the presence of multiple rods, bell cranks and ball bearings generate undesired friction forces between stick controls and main hydraulic actuators. In order to improve pilots comfort, an additional irreversible hydraulic assistance (see, Fig. 1) has been introduced into each flight axis control. For the sake of security, a fast stabilisation actuator (see Fig. 1) with low stroke capabilities has been introduced in series with each flight control in order to correct small deviations of the helicopter from its set point. In the following, the considered system is part of one helicopter flight axis control, from the stick to the stabilisation actuator, described in previous part and highlighted in Fig. 1. Part of this control chain has already been represented in [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF], highlighting two energetic subsystems: the first one is composed of the stick, the damper and the spring based element; it is dedicated to create a stick force feedback and to control the power of the second energetic system through the control of the hydraulic booster valve. The second energetic subsystem is the remaining part of the flight axis control: it creates the adequate force to control the rotor blade pitch angles from a given set point imposed by the hydraulic booster valve. The valve, which is the link between the two energetic systems, represents an energetic break in the flight axis control. In [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF], the aim was to remove the spring based element and to replace it by an active force feedback electromechanical system. The force feedback reference which is function of the stick position was obtained by empirical method comparing both systems. In this paper, a method has been carried out to determine from the actual process the reference to be applied on the active system. This method is verified on the entire process: other elements such as the hydraulic booster than the spring based element are removed. As a consequence, the two energetic subsystems will not be independent anymore. The entire representation of the flight axis control is therefore required to consider the interaction between these two subsystems. The main purpose of this section is to obtain an energetic graphical description of the entire flight axis control. The first step consists in a structural analysis of the control chain, leading to a lumped parameter model.

B. Energetic description

Bond Graph (BG) [START_REF] Paytner | Analysis and Design of Engineering systems[END_REF]- [START_REF] Karnopp | System Dynamics -Modeling and Simulation of Mechatronic Systems[END_REF] and Energetic Macroscopic Representation (EMR) [START_REF] Barre | Inversion-Based Control of Electromechanical Systems Using Causal Graphical Descriptions[END_REF] are both graphical and energetic oriented tools which simplify the analysis of complex multiphysic systems. A complementary approach between BG and EMR is dealt in [START_REF] Chikhaoui | Complementary use of BG and EMR formalisms for multiphysics systems analysis and control[END_REF]. Both of them represent energy flows in the system highlighting the power exchanged between each component. Whereas BG defines a bidirectional connection (power bond) between two components, the more recent energetic tool EMR describes power exchanges through bidirectional connection bringing out the effort and flow state variables.

EMR is a control oriented tool: it describes physical process imposing natural integral causality contrary to BG which allows a mathematical derivative description of physical accumulative processes. Integral causality representation makes thereafter easier the design of control architectures using a systematic model based inversion methodology. A large part of our researches are focused on the control of an active element justifying why EMR has been adopted. Fundamental energy processes are defined by simple specific elements. Depending on their power function, four main types of elements can be distinguished: energy sources and accumulation elements, conversion elements (without energy accumulation), and coupling elements for energy distribution. Based on equations of below section II.C and from [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF], the EMR of the studied helicopter flight axis control is given in Fig. 3.

C. Flight axis control chain modelling 1) From stick to hydraulic assistance control valve

Main characteristics of this first part of the model (lower left part of Fig. 3) are mathematically defined and represented in [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF]. The stick is mainly characterised by the stiffness Kstick which is nearly half of the total flight axis stiffness. The stick force feedback is mainly provided by a physical rotational spring Kspring in the case of the spring based element configuration. The damper is composed of multiple gear boxes assimilated to a single equivalent inertia, Jdamper. The mechanical link between the spring based element and the input of the hydraulic actuator is assimilated to a spring/mass system (Kl_rod, Ml_axis), see mathematical equations in [START_REF] Hautier | The causal ordering graph -A tool for system modelling and control[END_REF]. The last element of this first control chain part is the hydraulic assistance actuator control valve which is directly linked to the stick displacement. This valve is assimilated to a mass/spring system (Mvalve, Kvalve).

                _ _ _ _ 0 _ _ _ _ 0 1 t k rod l rod l rod l axis t l axis k rod l axis l axis F t K v T v T dT v t F T F T dT M              (1)

2) Hydraulic assistance

The hydraulic actuator is represented as an energy transformer (from hydraulic energy into mechanical work. A model of a hydraulic actuator is described in [START_REF] Fitzsimons | Part I: Modeling of a onedegree-of-freedom active hydraulic mount[END_REF], including the effects of servo-valve dynamics and fluid compressibility. A simplex hydraulic system with two symmetrical chambers is Energetic Macroscopic Representation of a Helicopter flight axis control ("hands on" case) considered. Hydraulic pressure losses are supposed to be insignificant. The control valve routes the fluid from the hydraulic source at pressure Ps to one chamber and links the second one to the return tank (Pr) (Fig. 2). The corresponding equations matching the fluid flow (Q1, Q2) and the pressure (P1, P2) of each chamber depending on the position of the control valve ud are summarised in [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF]. kd [kg/m] is a positive constant which characterises an hydraulic control valve.

                            (2)
The dynamic behaviour of the piston depends on pressure evolution in both chambers. The pressure variation inside one chamber is proportional to the difference between input and output flows and to the variation of the chamber volume. If leakage flows are neglected and if the actuator is symmetrical (meaning the useful surface of both chambers is the same and equal to S), the compressibility (βe) dynamic equations of the pressures in both chamber are given in [START_REF] Fitzsimons | Part I: Modeling of a onedegree-of-freedom active hydraulic mount[END_REF]. V0 is the hydraulic cylinder dead volume; upiston and vpiston are respectively the position and the speed of the hydraulic cylinder. The variations of fluid flows in both chambers make the piston moving. The corresponding hydro-mechanical conversion is given by equations (4).

                    11 0 0 22 0 0 t e P t Q T S v T dT piston V S u T piston t e P t Q T S v T dT piston V S u T piston                      (3)

 

12

F S P P piston Q S v piston piston          (4)
The hydraulic assistance is associated with a mechanical position regulating loop (Fig. 2): The valve position is mechanically linked to both the control chain from the pilot side, with speed vlr_byp and from the hydraulic cylinder speed vl_bp so as to ensure a piston displacement proportional to the control input whatever the load on the piston is. The kinematic equation of the mechanism is given by ( 5), where a and b are rod lengths:

3) Stabilisation actuator

The stabilisation actuator (upper right part of Fig. 3) is introduced in series in flight control architecture. Stabilisation process consists in driving the position of the output shaft actuator with regard to the helicopter attitude, measured through inertial data centres. The position reference is provided by the Auto Pilot Module (APM) to the stabilisation actuator which is then self-monitored (current and speed loops). The actuator mainly consists of a motor inertia Jstab, an inductance Lstab and an irreversible gear box [START_REF] Barre | Inversion-Based Control of Electromechanical Systems Using Causal Graphical Descriptions[END_REF]. The ball-screw stiffness is included in the stiffness Ku_rod.

                1 0 1 0 stab stab t i t U T e T dT stab stab stab L stab t t kt i t Tu T dT screw axis J stab                 (6)

4) Upper flight axis properties

The output of the hydraulic assistance actuator is modelled with a rigid body of mass, M_piston, including the hydraulic cylinder and stabilisation actuator masses (except from the output ball-screw mass which is included in equivalent mass M_load). M_load mainly represents the load mass of the upper flight axis. Stabilisation ball-screw and upper flight axis stiffnesses are represented by K_uRod. The corresponding equations are given in [START_REF] Chikhaoui | Complementary use of BG and EMR formalisms for multiphysics systems analysis and control[END_REF].

                        1 _ 0 _ _ _ 0 1 _ 0 t v t F T F T dT piston piston u axis M piston t F t K v T v T dT u axis u axis u axis stab t v t F T F T dT load u axis load M load                    (7)

III. INTEGRATION OF AN ACTIVE CONTROL DEVICE

Over time, designers have increased the number of elements of each flight axis control to improve pilots comfort. They added multiple simple elements. Simple is in the sense that each improvement (i.e. new function) was realised by additional equipment. Moreover, these improvements were based on relatively heavy technologies (e.g. hydromechanical), difficult to adjust and without possibilities of evolution (passive devices).

Current researches are more oriented to the development of new active flight control architectures including more functions, more flexibility, more security and more comfort with fewer elements. The objective is to transform part of the hydro-mechanical system into a new active system. Some elements have to be removed but their functional properties should be preserved by an extra new multi-purpose active system. One of the main challenges is to preserve the stick force feedback in helicopter flight axis control. A methodology based on a variation of the internal model control theory is presented in this section.

A. Internal Model Reference generator

Internal Model Control philosophy [START_REF] Morari | Robust Process Control[END_REF] relies on the following principle: control can be achieved only if the control system encapsulates, either implicitly or explicitly, some representation of the process to be controlled. If the control scheme has been developed based on an exact model of the process, then perfect control is theoretically possible. A variant of this method is proposed to reach our goal.

The general idea is to use part of the current hydromechanical control chain model as a reference generator for the new active system. Functional contributions of the removed elements should thus be preserved. Contrary to Internal Model Control philosophy in which the model approximates the real process, Internal Model Reference (IRM) generator is the representation of a reference process whose behaviour has to be reproduced on the new active process.

1) Internal Reference Model structure

The first step of the transformation is to delimitate on the actual process a unique group containing elements to be removed. These devices are physically removed from the process and the corresponding group defines the IRM structure. Two cases are taken into account depending on the interaction of the IRM and the rest of the model: either the IRM is in parallel or in series with the process:

-If the IRM elements were connected with the remaining process through a single power bond (parallel configuration, ie. a coupling), the IRM output signal (flow or effort) becomes a reference for the active process, and its input has to be measured on the real system. -If the group of elements to be removed is connected with the remaining process through two power bonds (series configuration), the IRM design is more complicated. Indeed, it has to generate two references from two measured signals, one at each bond. The active element should thus control both the action and reaction state variables, which is physically impossible. A compromise is necessary, choosing one state variable as a reference, the other being preserved by the addition of an extra transformer element (e.g. a gear box) within the physical part of the active system. In the case of the flight axis control, a Brushless Direct Current (BLDC) motor is used to control the stick force feedback. The ambition is not only to create a force feeback to replace the current spring, but also to assist the pilot against friction forces and to damp the flight axis control (resp. functions of the hydraulic assistance and the damper). The initial spring-based force feedback element, hydraulic actuator, damper and associated flight axis properties (kinematics, stiffness) form the IRM structure. These elements are connected in series with the flight axis control, i.e. between the stick and the stabilisation actuator (Fig. 1). One important constraint is to preserve the kinematic ratio between the stick angular position and the rotor blades angle. Moreover, the stick force feedback has to be kept dependant from its angular position. The new architecture should ideally preserve both the stick force feedback correlated to its position (input of the system), and the rotor blades angle (output of the system). As it is impossible, the first priority is to preserve the stick force feedback since kinematic elements can be added ("kin" transformer, close to "Bond 2", in Fig. 4) to guarantee the right ratio between the stick position and rotor blades angle. The new active element is then inserted in parallel with the flight control, on the spring based element coupling ("Bond 1" on Fig. 4). The EMR of the active system showing the IRM structure is represented in Fig. 4.

2) Description of the active element

The active element is composed of a brushless DC motor described by equations ( 8) (where Lcoil is its inductance, ktactive, its electromechanical conversion coefficient and Jactive, its inertia) associated with a gear box.

              1 () 0 1 _ 0 t i t V T kt T dT active cmd active active L coil t t kt i T T T dT active active active t active J active                   (8)

B. Control architecture design

In the previous part, the Internal Reference Model dedicated to the control of an active element inserted in an axis flight control has been defined. This part is focused on the design of the control architecture of this active device. The output of the IRM structure becomes the reference for the active element control. Control architecture of multi-element systems are often separately elaborated and locally tuned. Their regulating loops are usually PID controllers where gain parameters are directly set up on dedicated test benches. Difficulties in tuning controllers then appear when the system contains several energetic interactions, e.g. in case of complex multi-elements systems.

Control architecture of a system can be designed determining the inverse of its model. A model inversion technique is exposed and illustrated in [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF]-[6]- [START_REF] Martin | Modelling and control of an effort feedback actuator in helicopter flight control using Energetic Macroscopic Representation[END_REF]. These inversion rules are applied to each element along the chosen tuning path linking the state variable to be controlled and the tuning variable. In the present application, the tuning path links the controllable voltage source to the output of the active device: the output shaft speed (See Fig. 4). The tuning path is going through two energy storage elements. The first is the BLDC motor coil, electrical kinetic storage element (Lcoil) Its indirect inversion [START_REF] Barre | Inversion-Based Control of Electromechanical Systems Using Causal Graphical Descriptions[END_REF] leads to a current primary control loop. The second is the mechanical kinetic energy storage element of inertia Jactive and is associated to a speed control loop. The IRM block role is equivalent to the force control loop exposed in [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF].

IV. MODEL COMPARISON SIMULATION

Both hydro-mechanical and active systems have been simulated with Matlab/Simulink software. The challenge is, at equal stick angular position inputs, to reach an active force feedback behaviour matching with the hydro-mechanical system behaviour. The active system control has been obtained following the methodology described in previous section. The Internal Reference Model (IRM) represents removed parts of the hydro-mechanical process and the deduced reference is the set point value for the active device controller. Parameters of the speed and current controllers have been calculated following the Naslin criterion as exposed in [START_REF] Martin | Methodology for Modelling and Control of a Complex Multi-physic System -Case of helicopter flight axis control[END_REF].

Because friction forces have a major impact on the correct system functioning, Karnopp friction model has been used [START_REF] Karnopp | Computer Simulation of Stick-Slip Friction in Mechanical Dynamic Systems[END_REF]. The Karnopp model has the advantage of defining a clear stop condition when both speed and torque of the moving element are low.

The system is supposed to interact with the pilot ("hands on"). Studies of the relation between muscular activation patterns and movements show that human beings dynamic capabilities are limited and can be modelled with minimum-Jerk trajectories [START_REF] Harris | Exploring smoothness and discontinuities in human motor behavior with Fourier analysis[END_REF]. Such a smooth motion is generated thanks to a trapezoidal velocity signal filtered by a sliding mean filter in order to limit the Jerk.

The stick force law represents the stick force as a function of the stick position. The breakout point is the force level required to start a stick movement. The force is then proportional with the stick position. When the stick is dropped at an approximately stick angle of 11 degrees (time = 3s), it comes back and stops to its initial position.

Simulations results (Fig. 6.a) show that during the first part of the simulation ("hands on" configuration, the stick position is imposed by the pilot), the active system perfectly matches the hydro-mechanical system. The stick force is therefore accurately reconstructed from the Internal Model Reference. After three seconds of simulation, time from which the stick starts returning to its initial position, the stick is less damped in the active case than in the hydro-mechanical case. But in this unusual case, the pilots do not handle the stick and the stick Fig. 4 EMR of the active system and the resulting IRM based control structure ("hands on" case) force feedback signal is not of interest. However, it remains important that the stick position fit with the one of the hydromechanical system, which is the case: there is no over-shoot on the stick position and the stick angle error remains insignificant (less than 0.08°, see Fig. The model of a helicopter flight axis control is represented using the Energetic Macroscopic Representation. The model is mainly composed of hydro-mechanical elements. An evolution of the system is proposed in order to simplify the helicopter flight control architecture with active technology. The solution leads to a reduction of the number of hydro-mechanical elements, which are generally heavy and without evolution possibility, while preserving the global flight axis control behaviour.

A methodology is proposed to help the designers transforming the hydro-mechanical system into a new active system. The methodology is based on the Internal Model Reference (IRM) principle. The general idea is to simulate the behaviour of the hydro-mechanical system through the control architecture of the active element. The computed reference is then used as a speed set-point to drive the motor of the active element. The EMR energetic method has been used to define a systematic model based control architecture respecting natural causality. The speed and current controllers deduced from the inversion have been designed using the Naslin polynomial methodology. Both hydro-mechanical and active systems have been simulated using Matlab/Simulink software. Simulation results show that the IRM contained in the new active helicopter flight axis control contributes to the conservation of the hydro-mechanical system behaviour. The stick force feedback is preserved, validating the proposed methodology.

Future works are oriented on the simplification of the IRM. The more elements are removed (from the process to be simplified), the higher the order of the IRM is. The order of the IRM shall be decreased to reduce the complexity of the control and the computation time. The methodology would be based on the comparison of energy activity of each energy storage element included in the IRM.

Fig. 1 .

 1 Fig. 1. Helicopter flight axis control, a multi-physical system

  Fig. 3Energetic Macroscopic Representation of a Helicopter flight axis control ("hands on" case)

Fig. 2 .

 2 Fig. 2. Hydraulic assistance mechanical position regulating loop principle

  Fig. 6aStick force feedback vs. stick angular potision0 0.511.5 2 2.5 3 3.5 4 4.5 5 0 0.02 0.04 0.06 0.08 Stick angle error -Hydro-Mechanical and active model comparison Time (s)

ACKNOWLEDGMENT

This work has been supported by Olivier Honnorat, innovation project manager at Eurocopter, Marignane (Fr).


Gregory Kucherov

Kamil Salikhov

Dekel Tsur

Approximate String Matching using a Bidirectional Index

. We introduce a formalism, called search schemes, to specify search strategies of this type, then develop a probabilistic measure for the efficiency of a search scheme, prove several combinatorial results on efficient search schemes, and finally, provide experimental computations supporting the superiority of our strategies.

Introduction

Approximate string matching has numerous practical applications and long been a subject of extensive studies by algorithmic researchers. If errors are allowed in a match between a pattern string and a text string, most of the fundamental ideas behind exact string search algorithms become inapplicable. The Approximate string matching problem comes in different variants. In this paper, we are concerned with the indexing variant, when a static text is available for preprocessing and storing in a data structure, before any matching query is made.

The quest for efficient approximate string matching algorithms has been boosted by a new generation of DNA sequencing technologies, able to produce huge quantities of short DNA reads. Those reads can then be mapped to a given genomic sequence, which requires very fast and accurate approximate string matching algorithms. Many algorithms and associated software programs have been designed for this task, we refer to [START_REF] Li | A survey of sequence alignment algorithms for next-generation sequencing[END_REF] for a survey, and many of them rely on full-text indexes.

The classical indexing paradigm consists in building a text index in order to quickly identify pattern occurrences, preferably within a worst-case time weakly dependent on the text length. In the context of approximate matching, even the case of one error turned out to be highly nontrivial and gave rise to a series of works (see [START_REF] Lam | Improved approximate string matching using compressed suffix data structures[END_REF] and references therein). In the case of k errors, existing solutions generally have time or space complexity that are exponential in k, see [START_REF] Sung | Indexed approximate string matching[END_REF] for a survey.

Some of the existing algorithms use standard text indexes, such as suffix tree or suffix arrays. However, for large datasets occurring in modern applications, these indexes are known to take too much memory. Suffix arrays and suffix trees typically require at least 4 or 10 bytes per character respectively. The last years saw the emergence of succinct or compressed full-text indexes that occupy virtually as much memory as the sequence itself and yet provide very powerful functionalities [START_REF] Navarro | Compressed full-text indexes[END_REF]. For example, the FM-index [START_REF] Ferragina | Opportunistic data structures with applications[END_REF], based on the Burrows-Wheeler Transform [START_REF] Burrow | A block-sorting lossless data compression algorithm[END_REF], may occupy 2-4 bits of memory per character for DNA texts. FMindex has now been used in many practical bioinformatics software programs, e.g. [START_REF] Langmead | Ultrafast and memoryefficient alignment of short DNA sequences to the human genome[END_REF][START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF][START_REF] Simpson | Efficient de novo assembly of large genomes using compressed data structures[END_REF]. Even if succinct indexes are primarily designed for exact string search, using them for approximate matching naturally became an attractive opportunity. This direction has been taken in several papers, see [START_REF] Russo | Approximate string matching with compressed indexes[END_REF], as well as in practical implementations [START_REF] Simpson | Efficient de novo assembly of large genomes using compressed data structures[END_REF].

Interestingly, succinct indexes can provide even more functionalities than classical ones. In particular, several succinct indexes can be made bidirectional, i.e. can perform pattern search in both directions [START_REF] Belazzougui | Versatile succinct representations of the bidirectional burrows-wheeler transform[END_REF][START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF][START_REF] Russo | Approximate string matching with compressed indexes[END_REF][START_REF] Schnattinger | Bidirectional search in a string with wavelet trees and bidirectional matching statistics[END_REF]. Lam et al. [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF] showed how a bidirectional FM-index can be used to efficiently search for strings up to a small number (one or two) errors. The idea is to partition the pattern into k+1 equal parts, where k is the number of errors, and then perform multiple searches on the FM-index, where each search assumes a different distribution of mismatches among the pattern parts. Lam et al. implemented the proposed algorithm and reported that it outperforms in speed the best existing read alignment software [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF]. Related algorithmic ideas appear also in [START_REF] Russo | Approximate string matching with compressed indexes[END_REF].

In this paper, we extend the search strategy of [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF] in two main directions. We consider the case of arbitrary k and propose to partition the pattern into more than k + 1 parts that can be of unequal size. To demonstrate the benefit of both ideas, we first introduce a general formal framework for this kind of algorithm, called search scheme, that allows us to easily specify them and to reason about them (Section 2). Then, in Section 3 we perform a probabilistic analysis that provides us with a quantitative measure of performance of a search scheme, and give an efficient algorithm for obtaining the optimal pattern partition for a given scheme. Furthermore, we prove several combinatorial results on the design of efficient search schemes (Section 4). Finally, Section 5 contains comparative analytical estimations, based on our probabilistic analysis, that demonstrate the superiority of our search strategies for many practical parameter ranges. We further report on large-scale experiments on genomic data supporting this analysis.

Bidirectional search

In the framework of text indexing, pattern search is usually done by scanning the pattern online and recomputing index points referring to the occurrences of the scanned part of the pattern. With classical text indexes, such as suffix trees or suffix arrays, the pattern is scanned left-to-right (forward search). However, some compact indexes such as FM-index provide a search algorithm that scans the pattern right-to-left (backward search).

Consider now approximate string matching, where k letter mismatches are allowed between a pattern P and a text T . In this paper, we present our ideas for the case of Hamming distance. However, they also apply to the edit distance, see Conclusions. Both forward and backward search can be extended to approximate search in a straightforward way, by exploring all possible mismatches along the search, as long as their number does not exceed k and the current pattern still occurs in the text. For the forward search, for example, the algorithm enumerates all substrings of T with Hamming distance at most k to a prefix of P . Starting with the empty string, the enumeration is done by extending the current string with the corresponding letter of P , and with all other letters provided that the number of accumulated mismatches has not yet reached k. For each extension, its positions in T are computed using the index. Note that the set of enumerated strings is closed under prefixes and therefore can be represented by the nodes of a trie. Similarly to forward search, backward search enumerates all substrings of T with Hamming distance at most k to a suffix of P .

Clearly, backward and forward search are symmetric and, once we have an implementation of one, the other can be implemented similarly by constructing an index for the reversed text. However, combining both forward and backward search within one algorithm results in a more efficient search. To illustrate this, consider the case k = 1. Partition P into two equal length parts P = P 1 P 2 . The idea is to perform two complementary searches: forward search for occurrences of P with a mismatch in P 2 and backward search for occurrences with a mismatch in P 1 . In both searches, branching is performed only after |P |/2 characters are matched. Then, the number of strings enumerated by the two searches is much less than the number of strings enumerated by a single standard forward search, even though two searches are performed instead of one.

A bidirectional index of a text allows one to extend the current string S both left and right, that is, compute the positions of both cS or Sc from the positions of S. Note that a bidirectional index allows forward and backward searches to alternate, which will be crucial for our purposes. Lam et al. [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF] showed how the FM-index can be made bidirectional. Other succinct bidirectional indexes were given in [START_REF] Belazzougui | Versatile succinct representations of the bidirectional burrows-wheeler transform[END_REF][START_REF] Russo | Approximate string matching with compressed indexes[END_REF][START_REF] Schnattinger | Bidirectional search in a string with wavelet trees and bidirectional matching statistics[END_REF]. Using a bidirectional index, such as FM-index, forward and backward searches can be performed in time linear in the number of enumerated strings. Therefore, our main goal is to organize the search so that the number of enumerated strings is minimized.

Lam et al. [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF] gave a new search algorithm, called bidirectional search, that utilizes the bidirectional property of the index. Consider the case k = 2, studied in [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF]. In this case, the pattern is partitioned into three equal length parts, P = P 1 P 2 P 3 . There are now 6 cases to consider according to the placement of mismatches within the parts: 011 (i.e. one mismatch in P 2 and one mismatch in P 3 ), 101, 110, 002, 020, and 200. The algorithm of Lam et al. [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF] performs three searches:

(1) A forward search that allows no mismatches when processing characters of P 1 , and 0 to 2 accumulated mismatches when processing characters of P 2 and P 3 . This search handles the cases 011, 002, and 020 above.

(2) A backward search that allows no mismatches when processing characters of P 3 , 0 to 1 accumulated mismatches when processing characters of P 2 , and 0 to 2 accumulated mismatches when processing characters of P 1 . This search handles the cases 110 and 200 above.

(3) The remaining case is 101. This case is handled using a bidirectional search. It starts with a forward search on string P ′ = P 2 P 3 that allows no mismatches when processing characters of P 2 , and 0 to 1 accumulated mismatches when processing the characters of P 3 . For each string S of length |P ′ | enumerated by the forward search whose Hamming distance from P ′ is exactly 1, a backward search for P 1 is performed by extending S to the left, allowing one additional mismatch. In other words, the search allows 1 to 2 accumulated mismatches when processing the characters of P 1 .

We now give a formal definition for the above. Suppose that pattern P is partitioned into p parts. A search is a triplet of strings S = (π, L, U ) where π is a permutation string of length p over {1, . . . , p}, and L, U are strings of length p over {0, . . . , k}. String π indicates the order in which the parts of P are processed, and thus it must satisfy the following property: For every i > 1, π(i) is either (min j<i π(j)) -1 or (max j<i π(j)) + 1. Strings U and L give upper and lower bounds on the number of mismatches: When the j-th part is processed, the number of accumulated mismatches between the active strings and the corresponding substring of P must be between L[j] and U [j]. Formally, for a string M over integers, the weight of

M is i M [i]. A search S = (π, L, U ) covers a string M if L[i+1] ≤ i j=1 M [j] ≤ U [i] for all i (assuming L[p+1] = 0).
A k-mismatch search scheme S is a collection of searches such that for every string M of weight k, there is a search in S that covers M . For example, the 2-mismatch scheme of Lam et al. consists of searches S f = (123, 000, 022), S b = (321, 000, 012), and S bd = (231, 001, 012). We denote this scheme by S LLTWWY .

In this work, we introduce two types of improvements over the search scheme of Lam et al.

Uneven partition. In S LLTWWY , search S f enumerates more strings than the other two searches, as it allows 2 mismatches on the second processed part of P , while the other two searches allow only one mismatch. If we increase the length of P 1 in the partition of P , the number of strings enumerated by S f will decrease, while the number of strings enumerated by the two other searches will increase. We show that for some typical parameters of the problem, the decrease in the former number is larger than the increase of the latter number, leading to a more efficient search.

More parts. Another improvement can be achieved using partitions with k +2 or more parts, rather than k + 1 parts.

Analysis of search schemes

In this section we show how to estimate the performance of a given search scheme S. Using this technique, we present a dynamic programming algorithm for designing an optimal partition of a pattern.

Estimating the efficiency of a search scheme

To measure the efficiency of a search scheme, we estimate the number of strings enumerated by all the searches of S. We assume that performing single steps of forward, backward, or bidirectional searches takes the same amount of time. It is fairly straightforward to extend the method of this section to the case when these times are not equal. Note that the bidirectional index of Lam et al. [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF] reportedly spends slightly more time (order of 10%) on forward search than on backward search.

For the analysis, we assume that characters of T and P are randomly chosen uniformly and independently from the alphabet. We note that it is possible to extend the method of this section to a non-uniform distribution. For more complex distributions, a Monte Carlo simulation can be applied which, however, requires much more time than the method of this section.

Let #str(S, X, σ, n) denote the expected number of strings enumerated when performing a search S = (π, L, U ) on a random text of length n and random pattern of length m, where X is a partition of the pattern and σ is the size of the alphabet (note that m is not a parameter for #str since the value of m is implied from X). For a search scheme S, #str(S, X, σ, n) = S∈S #str(S, X, σ, n).

Fix S, X, σ, and n. Let A l be the set of enumerated strings of length l when performing the search S on a random pattern of length m, partitioned by X, and a text T containing all the strings of length at most m as substrings. Select a random order on the elements of A l , and let A l,i be the i-th element of A l . By the linearity of the expectation, #str(S, X, σ, n)

= l≥1 n l i=1 Pr T ∈Σ n [A l,i is a substring of T ],
where n l = |A l |. For any l and i, the string A l,i is a random string with uniform distribution over Σ l . Therefore, the probability that A l,i is a substring of T can be approximated by 1 -e -n/σ l using the Chen-Stein method [START_REF] Chen | Poisson approximation for dependent trials[END_REF]. Therefore, #str(S, X, σ, n) ≈ m l=1 n l (1 -e -n/σ l ).

(

) 1
In order to compute the values of n l , we give some definitions. Let n l,d be the number of strings in A l of length l with Hamming distance d to the prefix of P of length l. Let U ′ be a string obtained from U by replacing each character U [i] of U by a run of U [i] of length x π(i) , where x j is the length of the j-th part in the partition X. The string L ′ is defined analogously. In other words, the values

L ′ [i], U ′ [i]
give a lower and upper bounds on the number of allowed mismatches for an enumerated string of length i. The values of n l are given by the following recurrence.

n l = U ′ [l] d=L ′ [l] n l,d , n l,d =      n l-1,d + (σ -1)n l-1,d-1 if l ≥ 1 and L ′ [l] ≤ d ≤ U ′ [l] 1 if l = 0 and d = 0 0 otherwise (2)
For a specific search, a closed form formula can be given for n l . Consider equation [START_REF] Belazzougui | Versatile succinct representations of the bidirectional burrows-wheeler transform[END_REF]. The value of the term 1 -e -n/σ l is very close to 1 for l ≤ log σ n-O(1). When l ≥ log σ n, the value of this term decreases exponentially. Note that n l increases exponentially, but the base of the exponent of n l is σ -1 whereas the base of 1 -e -n/σ l is 1/σ. We can then approximate #str(S, X, σ, n) with function #str ′ (S, X, σ, n) defined by

#str ′ (S, X, σ, n) = ⌈log σ n⌉+cσ l=1 n l (1 -e -n/σ l ), (3)
where c σ is a constant chosen so that ((σ -1)/σ) cσ is sufficiently small. If a search scheme S contains two or more searches with the same π-strings, these searches can be merged in order to eliminate the enumeration of the same string twice or more. It is straightforward to modify the computation of #str(S, X, σ, n) to account for this optimization.

Computing an optimal partition

Let p be the number of parts. An optimal partition can be naively found by enumerating all m-1 p-1 possible partitions, and for each partition X, computing #str ′ (S, X, σ, n). We now describe a more efficient dynamic programming algorithm that computes an optimal partition for a given search scheme S.

The algorithm takes advantage of the fact that the value of #str ′ (S, X, σ, n) does not depend on the entire partition X, but only on the partition of a substring of P of length N = ⌈log σ n⌉ + c σ induced by X. We first give some definitions. Define S(i, j) to be the set containing every search S ∈ S such that i appears before j in the π-string of S (we assume S(i, p+1) = S). A partition of P is a partition of the set {1, 2, . . . , m} into disjoint sets of the form {i, i + 1, . . . , j} which will be called intervals. A partition of a substring P [i..j] is a partition of {i, i + 1, . . . , j} into disjoint intervals.

The algorithm builds partial partitions of prefixes of P , incrementally from left to right. That is, the partitions of P [1..m ′′ ] are built from partitions of P [1..m ′ ] for m ′ < m ′′ by extending them with interval {m ′ + 1, . . . , m ′′ }. Consider a partition X ′ of P [1..m ′ ], m ′ ≥ N , into p ′ parts, and let p ′ 2 be the left-toright rank of the interval in X ′ containing the position m ′ -N + 1. For a search S ∈ S(p ′ 2 , p ′ + 1), the value of #str ′ (S, X, σ, n) is the same for every extension of X ′ to a partition X of P , since #str ′ (S, X, σ, n) is determined by the partition of a substring of P of length N contained in P [1..m ′ ]. Thus, for every partial partition X ′ , the algorithm computes the value of S∈S(p ′ 2 ,p ′ +1) (S, X, σ, n) for some extension of X ′ to a partition X. This value will be denoted v(X ′ ). Note that for a partition X of P , v(X) = #str ′ (S, X, σ, n). When the algorithm extends the partial partition X ′ to a partial partition X ′′ = X ′ ∪ {m ′ + 1, . . . , m ′′ }, it has to compute v(X ′′ ). This is done by adding v(X ′′ ) -v(X ′ ) to the already computed v(X ′ ). By definition, v(X ′′ )-v(X ′ ) = S∈S(p ′′ 2 ,p ′′ +1)\S(p ′ 2 ,p ′ +1) #str ′ (S, X, σ, n), where p ′′ = p ′ + 1 and p ′′ 2 is the left-to-right rank of the interval in X ′ containing the position m ′′ -N + 1. If X ′ 1 and X ′ 2 are partitions of P [1..m ′ ] with the same number of parts and these partitions induce the same partition on

P [m ′ -N + 1..m ′ ], then for every partition X ′′ of P [m ′ + 1..m ′′ ] we have v(X ′ 1 ∪ X ′′ ) -v(X ′ 1 ) = v(X ′ 2 ∪ X ′′ ) -v(X ′ 2 )
. Thus, for every p ′ ≤ p and every partition X ′ of P ′ = P [m ′ -N + 1..m ′ ], we need to store only one partition X * of P [1..m ′ ] into p ′ parts that induces the partition X ′ on P ′ . The partition X * is chosen such that v(X * ) is minimum among all partitions of P [1..m ′ ] into p ′ parts that induce the partition X ′ on P ′ . We obtain an algorithm whose time complexity is

O(m 2 + (|S|N k + mp) p p ′ =1 N -1 p ′ -1
). Further details are omitted due to space limitations.

Properties of optimal search schemes

Designing an efficient search scheme for a given set of parameters consists of (1) choosing the number of parts, (2) choosing the searches, (3) choosing the partition of the pattern. While it is possible to enumerate all possible choices, and evaluate the efficiency of the resulting scheme using Section 3.1, this is generally infeasible due to a large number of possibilities. It is therefore desirable to have a combinatorial characterization of optimal search schemes.

The critical string of a search scheme S is the lexicographically maximal U -string of a search in S. A search of S is critical if its U -string is equal to the critical string of S. For example, the critical string of S LLTWWY is 022, and S f is the critical search. For typical parameters, critical searches of a search scheme constitute the bottleneck. Consider a search scheme S, and assume that the L-strings of all searches contain only zeros. Assume further that the pattern is partitioned into equal-size parts. Let ℓ be the maximum index such that for every search S ∈ S and every i ≤ ℓ, U [i] of S is no larger than the number in position i in the critical string of S. From Section 3, the number of strings enumerated by a search S ∈ S depends mostly on the prefix of the U -string of S of length ⌈⌈log σ n⌉/(m/p)⌉. Thus, if ⌈⌈log σ n⌉/(m/p)⌉ ≤ ℓ, a critical search enumerates an equal or greater number of strings than a non-critical search.

We now consider the problem of designing a search scheme whose critical string is minimal. Let α(k, p) denote the lexicographically minimal critical string of a k-mismatch search scheme that partitions the pattern into p parts. The next theorem gives the values of α(k, k + 2) and α(k, k + 1). We omit the proof due to space limitations. ) for 2-mismatch search schemes (recall that 4 is the size of the alphabet, and 4 16 is the length of the text). The second column gives #str values for the 3-part search scheme with equal-size parts. The other columns give #str values for different search schemes using an optimal partition of the pattern. For each search scheme, the optimal value of #str is shown in the first sub-column, and the optimal partition in the second sub-column. 5 Case studies

m
In this section, we examine the efficiency of several 2-mismatch and 3-mismatch search schemes. The search schemes were generated by a greedy algorithm. At each step, the algorithm considers the uncovered string M of weight k such that the lexicographically minimal U -strings of searches that covers M is maximal. Among the searches that cover M with minimal U -string, a search that covers the maximum number of uncovered strings of weight k is chosen. The L-string of the search is chosen to be lexicographically maximal among all possible L-string that do not decrease the number of uncovered strings. For each search scheme and each choice of parameters, we obtained an optimal partition and computed the efficiency of the scheme according to Section 3. Results for 2 mismatches are given in Table 1 and Table 2, for 4-letter and 30letter alphabets respectively, and results for 3 mismatches are given in Table 3.

Our theoretical analysis indicates that an increased number of parts combined with uneven partitioning improves over a scheme with k + 1 equal sized parts when m/(k + 1) is smaller than log σ n (details omitted due to lack of space). This is confirmed by the numerical results in Tables 1,2, and 3, which show significant improvement when m/(k + 1) is small.

For big alphabets (Table 2), we observe a larger gain in efficiency. The is due to the fact that the values of n l (see equation ( 2)) increase more rapidly when the alphabet is large, and thus a change in the size of parts can have a bigger influence on these values.

For 3 mismatches (Table 3), we observe smaller gain. This is partly explained by Theorem 1, as with 3 mismatches and 4 parts, the critical string starts with 01 (compared to 02 for 2 mismatches and 3 parts), therefore 4 parts provide already a competitive search. Another interesting observation is that with 4 parts, the optimal partition is an even one, as the U -strings in all searches in the 4-part scheme are the same.

We implemented our method using the 2BWT library, provided by [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF] (available at http://i.cs.hku.hk/2bwt-tools/) and experimentally compared different search schemes. The experiments were done on the sequence of human chromosome 14. The sequence is 88M long, with nucleotide distribution 29%, 21%, 21%, 29%. Searched patterns were obtained from Next-Generation Sequencing reads by cutting out strings of required length. Both the sequence and the reads were downloaded from http://gage.cbcb.umd.edu/data/. For every pattern length and every search scheme, 10 6 patterns were searched and the average number of enumerated strings was computed.

For the case of 2 mismatches, we implemented the 3-part and 4-part schemes from Section 5, as well as their equal part versions for comparison. For each pattern length, we computed an optimal partition, taking into account a nonidentical distribution of nucleotides. Results are presented in Table 4. Using unequal parts for 3-part schemes yields a notable time decrease (8-17%) for patterns of length between 21 and 33. Furthermore, using 4 parts leads to an even more important improvement, showing a significantly better results for all pattern length compared to the 3-equal-parts scheme of [START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF]. For pattern lengths from 21 to 42 we observe 16-36% improvement in running time and 12-35% improvement in the number of enumerated strings. For pattern lengths 24, 27, 39, 42, we observe that using unequal part lengths for 4-part schemes is beneficial. Overall, the experimental results are consistent with numerical estimations of Section 5. However, for pattern lengths 30-36, the 4-equal-parts scheme performs better than the 4-unequal-parts one, which illustrates that the optimal partition found for random texts may not be the best one for genomic sequences.

Conclusions

This paper can be seen as the first step towards an automated design of efficient search schemes for approximate string matching, based on bidirectional indexes.

More research has to be done in order to allow an automated design of optimal search schemes. It would be very interesting to study an approach when a search scheme is designed simultaneously with the partition, rather than independently as it was done in our work. The results of this paper can be extended to approximate string matching under edit distance. The estimation of n l (Section 3) becomes more complicated though.

Theorem 1 .

 1 α(k, k + 1) = 013355 • • • kk for every odd k, and α(k, k + 1) =

Table 1 .

 1 Values of #str(S, X, 4, 416

Table 2 .

 2 Values of #str(S, X, 30, 30 7 ) for 2mismatch search schemes.

	m 3 equal 3 unequal 4 unequal 5 unequal
	15 846 286 6,4,5 231 5,2,3,5 286 5,1,3,1,5
	18 112 111 7,6,5 81 6,2,4,6 111 6,1,4,1,6
	21 24	24 7,7,7 23 7,3,4,7 24 7,1,6,1,6

Table 3 .

 3 Values of #str(S, X, 4, 416 ) for 3mismatch search schemes.

	m	4 unequal	5 unequal
	24 11222 6,6,6,6 8039 4,6,5,1,8
	36 416	9,9,9,9	549 6,11,5,1,13
	48 185 12,12,12,12 213 11,11,11,1,14

Table 4 .

 4 Total time (in seconds) of searching for one million patterns in human chromosome 14, up to 2 mismatches.

	m 3 equal	3 unequal	4 equal	4 unequal
	18 230 230 (100%) 6,6,6 209 (91%) 203 (88%) 5,5,1,7
	21 175 161 (92%) 8,6,7 147 (84%) 150 (86%) 6,4,3,8
	24 142 120 (85%) 10,7,7 117 (82%) 107 (75%) 7,4,4,9
	27 119	99 (83%) 12,7,8 96 (81%) 82 (69%) 9,4,5,9
	30 101	84 (83%) 12,9,9 66 (65%) 68 (67%) 10,4,6,10
	33 83	70 (84%) 13,10,10 53 (64%) 58 (70%) 11,5,6,11
	36 68	66 (97%) 13,11,12 49 (72%) 50 (74%) 12,5,7,12
	39 56	56 (100%) 13,13,13 48 (86%) 45 (80%) 13,6,7,13
	42 45	45 (100%) 14,14,14 44 (98%) 38 (84%) 14,6,8,14

Acknowledgements. GK has been supported by the ABS2NGS grant of the French government (program Investissement d'Avenir ) as well as by a EU Marie-Curie Intra-European Fellowship for Carrier Development. KS has been supported by the co-tutelle PhD fellowship of the French government. DT has been supported by ISF grant 981/11.


Gregory Kucherov

Dekel Tsur

Improved filters for the approximate suffix-prefix overlap problem

Computing suffix-prefix overlaps for a large collection of strings is a fundamental building block for the analysis of genomic next-generation sequencing data. The approximate suffix-prefix overlap problem is to find all pairs of strings from a given set such that a prefix of one string is similar to a suffix of the other. Välimäki et al. (Information and Computation, 2012) gave a solution to this problem based on suffix filters. In this work, we propose two improvements to the method of Välimäki et al. that reduce the running time of the computation.

Introduction

Genomic sequences are deciphered by reading short overlapping fragments. Modern next-generation sequencing technologies, can produce, in a single run, tens or hundreds of millions of such fragments, called reads, each of the order of a hundred of letters. Dealing with these gigabytes of sequence data raises a number of algorithmic challenges.

A basic operation on a collection of genomic reads is the computation of overlaps: we need to be able to quickly retrieve reads which have a significant overlap with a given read. This operation is a prerequisite for many algorithms dealing with reads, and most prominently for genome assembly algorithms which follow the so-called overlap-layout-consensus paradigm [START_REF] Li | Comparison of the two major classes of assembly algorithms: overlap-layout-consensus and de-Bruijn-graph[END_REF]. These algorithms are based on overlap graphs (also called string graphs or assembly graphs) that represent all significant overlaps between reads. A recent example is provided by SGA assembler [START_REF] Simpson | Efficient de novo assembly of large genomes using compressed data structures[END_REF]. Earlier, this approach was taken by several "first-generation" methods of genome assembly, such as Celera assembler [START_REF] Myers | A whole-genome assembly of Drosophila[END_REF] used to assemble one of the first versions of the human genome.

The goal of this work is to propose an efficient way of computing significant approximate suffix-prefix overlaps of a set of strings. Previously, several solutions have been proposed to compute all exact suffix-prefix overlaps [START_REF] Gusfield | An efficient algorithm for the all pairs suffix-prefix problem[END_REF][START_REF] Ohlebusch | Efficient algorithms for the all-pairs suffix-prefix problem and the all-pairs substring-prefix problem[END_REF][START_REF] Simpson | Efficient de novo assembly of large genomes using compressed data structures[END_REF]. However, in practice, we are interested in the approximate case when strings can overlap within a certain number of errors.

Most practical methods for computing approximate string similarities are based on the filtering approach, when the search is done in two steps: at the first step, candidate regions are identified that potentially correspond to sought matches, and at the second step, those candidates are checked to actually verify the desired matching condition. Filtering algorithms usually do not yield interesting theoretical time bounds but are often very efficient in practice. As an example, spaced seeds [START_REF] Burkhardt | Better filtering with gapped q-grams[END_REF][START_REF] Ma | PatternHunter: Faster and more sensitive homology search[END_REF] constitute one of the filtering techniques that has been successfully used for DNA sequence comparison, e.g. [START_REF] Ma | PatternHunter: Faster and more sensitive homology search[END_REF][START_REF] Noé | YASS: enhancing the sensitivity of DNA similarity search[END_REF].

To compute approximate suffix-prefix overlap, the above-mentioned SGA assembler [START_REF] Simpson | Efficient de novo assembly of large genomes using compressed data structures[END_REF] uses a basic substring filtering. Välimäki et al. [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] proposed to apply a modified version of suffix filters earlier proposed by Kärkkäinen and Na [START_REF] Kärkkäinen | Faster filters for approximate string matching[END_REF]. Suffix filter provide a more selective filtering criterion and therefore a more efficient algorithm.

In this paper, we show how the method of [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] can be further improved. We propose two improvements that reduce the search space and therefore the running time of the algorithm: a new family of suffix filters and a new partitioning scheme. We report on estimations on random datasets that support the superiority of our schemes.

Throughout the paper, we present our method for the Hamming distance between strings, although it can be generalized to the edit distance, similar to [START_REF] Kärkkäinen | Faster filters for approximate string matching[END_REF][START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF]. This, however, would entail additional technical details and a more involved presentation that we wanted to avoid.

Preliminaries

Notation

For a sequence of integers A, let PrefixSum(A) be a sequence of integers of the same length as A in which PrefixSum(A)[i] = Let A be a string that has a partition

A = A 1 A 2 • • • A k into k disjoint parts.
Let B be a string with the same length as A. The partition of A induces a partition

B = B 1 B 2 • • • B k of B in which |B i | = |A i | for all i.
The partition distance between A and B, denoted pd(A, B), is a sequence of integers of length k in which pd(A, B)[i] = Ham(A i , B i ). The accumulated partition distance between A and B, denoted apd(A, B), is the sequence PrefixSum(pd(A, B)). That is, apd(A, B)[i] is the total number of mismatches between the first i parts of A and B.

Suffix filters

For the problem of approximate pattern matching with q errors, the most basic filtering method, called PEX in [START_REF] Navarro | Flexible Pattern Matching in Strings -Practical on-line search algorithms for texts and biological sequences[END_REF], consists in splitting the pattern into k = q +1 parts and searching for those parts independently. Once one of the parts is found, it provides a candidate location for the whole pattern. Kärkkäinen and Na [START_REF] Kärkkäinen | Faster filters for approximate string matching[END_REF] proposed an interesting extension of this principle called suffix filters. Suffix filters have been shown to generate many fewer candidates than substring filters, and therefore to be much more efficient. Since the present work builds on this idea, we briefly explain it here.

We still split pattern P into k = q + 1 parts P = P 1 P 2 • • • P k , but instead of searching for substrings P i , we will be searching for suffixes P i P i+1 • • • P k for i = 1, . . . , k allowing errors distributed according to a specific pattern

filter k,i = 01 • • • (k -i).
This means that for every i, we are searching for strings B such that apd(

P i • • • P k , B) ≤ filter k,i .
The key observation of [START_REF] Kärkkäinen | Faster filters for approximate string matching[END_REF] is that this scheme detects all possible occurrences of P within q errors. For example, let q = 2, namely, we want to find the substrings of the text that 2-match to P = P 1 P 2 P 3 . All such substrings are detected using three filters, denoted by sequences 012, 01, and 0. Filter 012 detects substrings B such that apd(P 1

P 2 P 3 , B) ≤ 012. That is, B = B 1 B 2 B 3 , such that |B i | = |P i | for all i, B 1 = P 1 , Ham(B 2 , P 2 ) ≤ 1,
and Ham(B 2 B 3 , P 2 P 3 ) ≤ 2. By a slight abuse of language, the set of all such strings B will be said to be enumerated by the filter. Similarly, filter 01 detects substrings B = B 2 B 3 such that B 2 = P 2 and B 3 is within one error from P 3 . Each such string is a suffix of a candidate approximate occurrence of P . Finally, filter 0 detects substrings B 3 with B 3 = P 3 , which provides again a suffix of a candidate occurrence of P .

Observe that there are 9 cases for the partition distance between P and B -011, 101, 110, 002, 020, 200, 100, 010, and 001 -which are all covered by suffix filters 012, 01 and 0. Indeed, filter 012 covers cases 011, 002, 010, and 001, filter 01 covers cases 101, 200 and 100, and filter 0 covers cases 020 and 110.

The set of strings enumerated by a filter can be naturally represented by a trie (see Figure 1), where branching nodes correspond to positions where the filter allows a possible mismatch to occur. The number of nodes in the tries of all the filters is a crucial parameter for the efficiency of a filtering scheme.

Suffix filters and full-text indexes

One of the advantages of suffix filters (as opposed e.g. to spaced seeds) is that they can be naturally implemented using full-text indexes that support incremental string matching. Those indexes include classical indexes such as suffix trees, but also succinct indexes such as FM-index and its variants [START_REF] Vyverman | Prospects and limitations of full-text index structures in genome analysis[END_REF]. These indexes allow reading a pattern left-to-right3 and quickly updating the index point after each letter, so that all occurrences of the prefix read so far can be retrieved efficiently.

Implementing suffix filters on a full-text index can be done simply by enumerating all "approximate suffixes" of the pattern detected by a filter, and reporting the occurrences of all these strings in the text, thus generating the candidate set. More precisely, for each of the tries that correspond to the filters, the algorithm traverses the trie in depth-first search order, and updates the index point after each descend in the trie. When reaching a node in the trie whose corresponding string does not appear in the index, the search does not continue to the descendants of this node in the trie. When reaching a leaf of the trie, all occurrences of the string corresponding to the leaf are retrieved from the index, and are added to the set of candidates.

Suffix filters applied to suffix-prefix overlap problem

Given a set S of strings, the approximate suffix-prefix overlap problem is to compute all significant approximate overlaps between pairs of strings of S. "Significance" is defined by a lower threshold l min on the overlap size. Since the overlap size is variable, imposing a fixed number of errors is not reasonable, and a relative error rate is specified instead. Formally, given an integer l min , and ǫ > 0, we have to find all pairs of strings S, S ′ ∈ S such that there is an integer l ≥ l min for which the prefix of S of length l and the suffix of S ′ of length l ⌈ǫl⌉-match.

In this Section, we explain how suffix filters can be used to solve the approximate suffix-prefix overlap problem. We proceed by enumerating all strings S ∈ S. For each S ∈ S and for each l ≥ l min , we want to identify all strings S ′ ∈ S whose suffix of length l ⌈ǫl⌉-match the prefix S[1..l]. We want to do it by applying suffix filters designed for patterns of length l and ⌈ǫl⌉ mismatches. If such a filter applies, then a candidate overlap is generated, which is a triplet (S, S ′ , l). At the verification stage, the actual Hamming distance between the prefix of S of length l and the suffix of S ′ of length l is computed, and the overlap is reported if this distance is no more than ⌈ǫl⌉.

Let us now explain how the filtering algorithm works with the filters of Section 2.2. Fix S and a partition S = S 1 S 2 • • • of S into disjoint parts. Fix l ≥ l min and let P = S[1..l]. The partition of S induces a partition P = P 1 • • • P k of P , where P i = S i for i < k and P k is a prefix of S k . Suppose that the partition of S was chosen a way to ensure that k ≥ ⌈ǫl⌉ + 1. This allows us to apply suffix filters of Section 2.2.

Consider the above filters filter k,1 , . . . , filter k,k , where filter k,i = 01 • • • (k -i). Each filter filter k,i enumerates all strings B such that apd(P i • • • P k , B) ≤ F , and B appears as a substring of some string in S. For each such string B, the algorithm further selects all the strings S ′ ∈ S that end with B and adds the triplets (S, S ′ , l) to the list of candidates.

The main difficulty in applying suffix filters to the approximate suffix-prefix overlap problem is that the length l and, consequently, the number of errors are not fixed. Once the partition of S is defined, our goal is to deal with all values of l in one left-to-right traversal of S. For each l ≥ l min , we should be able to efficiently identify appropriate filters that apply to the corresponding partition of S[1..l]. We now describe how it is done for the filters of Section 2.2.

A key point here is that the enumeration processes for different values of l are connected. Let k l denote the number of parts in the partition of S[1..l] that is induced by the partition of S. Let B S,i,l be the set of strings enumerated by the filter k l ,i when it is applied to S[1..l]. That is, B S,i,l are the strings enumerated by filter k l ,i when considering the filtering scheme for the fixed value of l. Let trie S,i be the trie representing the set B S,i,|S| . We have the following property: For every l, where l min ≤ l < |S|, and every i ≤ k l , the prefix of filter k |S| ,i of length k l + 1 -i is equal to filter k l ,i . It follows that the set B S,i,l is equal to the set of strings that correspond to the nodes of trie S,i at depth l -j<i |S j |. Therefore, generation of candidates can be done for all values of l by the following algorithm. For i = 1, 2, . . . , k S , traverse trie S,i . When the traversal is at a node that corresponds to a string B, if B ∈ B S,i,l for some l ≥ l min , find all the strings S ′ ∈ S that ends with B, and for each such S ′ add the triplet (S, S ′ , |B| + j<i |S j |) to the list of candidates.

Checking whether B ∈ B S,i,l is done with the following candidate generation condition.

Condition 1 B ∈ B S,i,l if and only if

|B| + j<i |S j | ≥ l min (C1)
2.5 The filtering scheme of Välimäki et al. [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] Välimäki et al. [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] observed that the filtering procedure of Section 2.4 is very inefficient. The inefficiency is caused by the filters filter k,k = 0 that has to be applied, during the search, to short strings including those consisting of a single letter. Formally, when traversing the trie trie S,i , a node of the trie at depth 1, whose corresponding single-letter string is B = S[1 + j<i |S j |], can generate candidates if Condition (C1) is satisfied. Since B has length 1, there can be many strings S ′ ∈ S that ends with B, generating many spurious candidates. More generally, assuming the strings of S are sampled randomly from an i.i.d. source, the expected number of candidates generated by a string B corresponding to some node in a trie is m/σ |B| , where m is the number of strings in S, and σ is the size of the alphabet. Therefore, the total number of candidates is dominated by the number of candidates generated by nodes of small depth. The solution of Välimäki et al. to this problem is to drop the filters filter k,k from the filtering scheme. In order to handle the cases of partition distances that were covered by this filter, filters filter k,1 = 01

• • • k are modified to filter k,1 = 12 • • • (k + 1)
. It is shown [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] that with this modification, all combinations of partition distance are still covered.

Due to the dropping of filters filter k,k , Condition 1 should be modified for the new filtering scheme. The modified candidate generating condition will now be as follows.

Condition 2 B ∈ B S,i,l if and only if Condition (C1) is satisfied and

|B| > |S i | (C2)
Clearly, the additional condition (C2) reduces the number of generated candidates.

3 New filtering scheme

In the filtering scheme of Välimäki et al., the filters filter k,1 start with 1 whereas the other filters start with 0. This has the consequence that the number of nodes in the trie trie S,1 is much larger than the number of nodes in the tries trie S,i for i > 1. We now present a new filtering scheme without this property. Our filtering scheme consists of filters filter k,1 , . . . , filter k,k-1 , where filter k,i is a sequence of length k -i + 1 whose first k -i elements are 0, 1, . . . , (k -i -1), and the last element is k -i -1. For example, for k = 4 the filters are filter 4,1 = 0122, filter 4,2 = 011, and filter 4,3 = 00. Our filtering scheme requires a difference of at least 2 between the number of parts in the partition of S[1.

.l] and ⌈ǫl⌉ (recall that in the scheme of Välimäki et al., the required difference is at least 1). We show the correctness of this scheme in the following lemma.

Lemma 1. For every k ≥ 2 and every sequence of integers M of length k whose sum is at most k -2, there is an integer i such that

PrefixSum(M [i..k]) ≤ filter k,i .
Proof. The proof is by induction of k. The base of the induction k = 2 is trivial since in this case M = 00 and therefore PrefixSum(M ) ≤ filter k,k-1 . Now consider k > 2. Let M be a sequence of length k whose sum is at most k -2. If PrefixSum(M ) ≤ filter k,1 we are done. Otherwise, there is an index i such that PrefixSum(M

)[i] > filter k,1 [i]. Since PrefixSum(M )[i] ≤ k -2 and filter k,1 [j] = k-2 for j > k-2, it follows that i ≤ k-2. Therefore filter k,1 [i] = i-1 and PrefixSum(M )[i] ≥ i. Let M ′ = M [i + 1..k]. The length of M ′ is k -i and the sum of M ′ is at most k -2 -PrefixSum(M )[i] ≤ k -2 -i. By the induction hypothesis, there is an index i ′ such that PrefixSum(M ′ [i ′ ..k -i]) ≤ filter k-i,i ′ . The lemma follows since M ′ [i ′ ..k -i] = M [i + i ′ ..k] and filter k-i,i ′ = filter k,i+i ′ . ⊓ ⊔
To illustrate Lemma 1, observe that the three above-mentioned filters filter 4,i , i = 1, 2, 3 cover all possible partition distances for the case of 4 parts and 2 errors. Indeed, filter 0122 cover cases 0020, 0002, 0101 and 0110, filter 011 cover cases 0011, 1001 and 1010, and finally filter 00 covers cases 2000, 0200 and 1100.

In addition to reducing the number of nodes in the tries, our filtering scheme also reduces the number of generated candidates. In the filtering scheme of Välimäki et al., filter k,i is a prefix of filter k ′ ,i for all k < k ′ and i. Our filtering scheme does not have this property. Therefore, we need a new condition for checking whether B ∈ B S,i,l .

Condition 3 B ∈ B S,i,l if and only if conditions (C1) and (C2) are satisfied and apd(S ′ , B)

≤ filter k,i (C3)
where S ′ is the prefix of S i S i+1 • • • of length |B| and k is the number of parts in the partition of S ′ induced by the partition of S.

Observe that in the filtering scheme of Välimäki et al., every node of trie S,i with depth |S i | + 1 generated candidates. However, in our new scheme, only the node whose corresponding string B is equal to S ′ generates candidates among the nodes of depth |S i | + 1 (for every other node, apd(S ′ , B) = 01 and therefore Condition (C3) is not satisfied). The same is true for depths

|S i | + 2, . . . , |S i | + |S i+1 |.
Our scheme can be generalized by introducing a parameter s ≥ 2. The filters for a given value of s are filter k,i = 01 • • • (k -i -s)(k -i -s + 1) s for i = 1, . . . , k -s + 1. This scheme requires a difference of at least s between the number of parts in the partition of S[1..l] and ⌈ǫl⌉.

Partition schemes

The efficiency of the algorithm of the previous section depends on the sizes of the parts in the partition of S: having larger parts reduces the number of trie nodes and the number of candidates. For correctness of the algorithm, the partitioning of a string S must satisfy the following property.

(P1) For every l ≥ l min , the number of parts in the partition of S[1..l] induced by the partition of S is at least ⌈ǫl⌉ + s, where s = 1 for the filtering scheme of Välimäki et al. and s = 2 for our new scheme (or some fixed value of s for our extended scheme).

Välimäki et al. used a partition of S into equal sized parts of size p, except for the last part whose size is at most p. The value of p is chosen to be the maximum integer for which Property (P1) is satisfied.

We propose a partitioning scheme in which most parts are larger than those of the equal sized parts partitioning. Since the efficiency of the filtering approach depends on the sizes of the parts, the new partitioning scheme gives better performance.

Let S be a string to be partitioned. Let l 0 < l 1 < • • • < l q be all the indices in the range l min , . . . , |S| for which ⌈ǫ(l -1)⌉ < ⌈ǫl⌉, and let l q+1 = |S| + 1. Let k = ⌈ǫl 0 ⌉ + s -1. We partition S as follows. The sizes of the first k parts are chosen in order to satisfy the following properties.

1. The total length of the first k parts is l 0 -1. 2. The length of the k-th part is at least l 0 -l min . We can set for example the length of the k-th part to be L = max(⌈(l 0 -1)/k⌉, l 0 -l min ), and set the lengths of the first k -1 parts to be p or p + 1, where p = ⌊(l 0 -1 -L)/(k -1)⌋. The lengths of the remaining parts in the partition are l 1 -l 0 , l 2 -l 1 , . . . , l q+1 -l q . We now show that this partition satisfies Property (P1). Moreover, the partitioning is optimal in the sense that the inequalities of Property (P1) are satisfied with equality.

Lemma 2. For every l ≥ l min , the number of parts in the partition of S[1..l] induced by the partition of S is ⌈ǫl⌉ + s.

Proof. For l min ≤ l < l 0 (assuming l 0 > l min ) we have by construction that the number of parts in the partition of S[1..l] is k = ⌈ǫl 0 ⌉ + s -1. By the definition of l 0 , ⌈ǫl⌉ = ⌈ǫl 0 ⌉ -1, so the equality of the lemma is satisfied. Similarly, if l i ≤ l < l i+1 then the number of parts in the partition of S[1..l] is k + 1 + i. Moreover, ⌈ǫl⌉ = ⌈ǫl 0 ⌉ + i. Therefore, the lemma follows.

⊓ ⊔

As an example, consider partitioning a string of length 200 with the parameters l min = 40 and ǫ = 0.1. If s = 1, the equal sized partition uses parts of size 8. In our new partitioning scheme, the first 5 parts have size 8, and the remaining parts have size 10.

Experimental results

In this Section, we compare the performance of 4 filtering schemes:

(1) the filtering scheme of Välimäki et al. [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] using partitioning into equal parts, (2) the filtering scheme of Välimäki et al. combined with our new partitioning scheme (Section 4), (3) our new filtering scheme (Sections 3-4), and (4) our extended filtering scheme (see end of Section 3), with s = 3.

For our filtering scheme (3), string partitioning is done with our new partitioning scheme. The comparisons have been done using the technique described in Kucherov et al. [START_REF] Kucherov | Approximate string matching using a bidirectional index[END_REF]. In this analysis, we assume that characters of the strings of S are randomly chosen uniformly and independently from the alphabet. Under this assumption, we analytically estimate the expected number of nodes in the tries Table 1. Expected performance of the filtering scheme of Välimäki et al. [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] and our filtering schemes. It is assumed that S contains m random strings of length 300 over an alphabet of size 4. For each scheme, the first column shows the expected number of nodes in the tries trieS,i for all i (for a single S ∈ S), and the second column is the expected number of candidates generated for S. m lmin ǫ method of [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] [13] with unequal parts our scheme our scheme extended for s = 3 10 6 20 0. Note that for different parameters, the bottleneck of the computation can be either the size of traversed tries, or the number of generated candidates. In all cases, we observe a significant decrease of both these measures compared to the original method of Välimäki et al. [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF]. When the threshold l min is small (in our experiments, 20 for the sequence length 300), the filters of [START_REF] Välimäki | Approximate all-pairs suffix/prefix overlaps[END_REF] combined with our partitioning scheme presents a trade-off with our filtering scheme: our scheme yields a smaller number of traversed trie nodes but a larger number of generated candidates. However, when l min is large enough (in our experiments, 40 for the sequence length 300), our scheme outperforms the one of Välimäki et al. in both the number of nodes in the tries and the number of generated candidates. The extended scheme with s = 3 yields a smaller number of candidates, but at the cost of increased number of traversed trie nodes.

Conclusions

In this paper, we proposed an improved filtering scheme for the approximate suffix-prefix overlap problem directly raised by bioinformatics applications. Two improvements are proposed: we provide a more efficient filtering scheme as well as new way of partitioning the query string. We show, through analytical estimations, the superiority of our scheme in terms of the size of the search space (size of traversed tries) as well as the selectivity (number of generated candidates).

Several directions for future work can be envisaged. We did not compare the actual performance of the different filtering schemes on real data. However, previous work [START_REF] Kucherov | Approximate string matching using a bidirectional index[END_REF] provides strong grounds to assume that the better performance will be supported by real data too. This, however, remains to be verified experimentally. Another direction, already mentioned in Introduction, concerns the generalization of our results to the case of edit distance. While we don't expect significant obstacles in this generalization, it does bring an additional technical difficulty.

  i j=1 A[j]. For two sequences of integers A and B of the same length, we define A ≤ B iff A[i] ≤ B[i] for all i. The Hamming distance between strings A and B of the same length, denoted Ham(A, B), is the number of indices i for which A[i] = B[i]. If Ham(A, B) ≤ k, we say that A and B k-match.

Fig. 1 .

 1 Fig. 1. The tries of the filters 012, 01, and 0 applied on the pattern P = 000110, for binary alphabet and q = 2 errors. The pattern P is partitioned into 3 parts P = P1P2P3 where P1 = 00, P2 = 01, and P3 = 10. The filter 012 enumerates the strings corresponding to the leaves of the tries in Figure (a), namely 000110, 000111, . . . , 001100. The filter 01 enumerates the strings 0110, 0111, and 0100, and the filter 0 enumerates the string 10.

  S,i and the expected number of generated candidates, following the method developed in[START_REF] Kucherov | Approximate string matching using a bidirectional index[END_REF]. The results are summarized in Table1. Columns 2 to 5 of the Tablecorrespondto schemes (1) to (4) above, respectively.

	1 31257 18950 11782	144 8582 341 20026	95
	10 7 20 0.1 64201 189506 22161 1449 13219 3412 46662	952
	10 6 40 0.1 10416	839 8260	65 6912	28 8868	0.5
	10 7 40 0.1 14391 8391 11138	651 8921 280 12916	4
	10 7 40 0.15 207504 857318 71271 82559 40671 18842 82164	116

trie

For some indexes, such as FM-index, matching is performed right-to-left, in which case we can just assume that the indexed sequences are reversed. On the other hand, FM-index can also be modified to perform matching left-to-right, see[START_REF] Lam | High throughput short read alignment via bi-directional BWT[END_REF].

Acknowledgements. GK has been supported by the ABS2NGS grant of the French government (program Investissement d'Avenir ) as well as by a EU Marie-Curie Intra-European Fellowship for Carrier Development. DT has been supported by ISF grant 981/11.