Stéphane Bazeille

Vincent Lebastard

Sylvain Lanneau

Frédéric Boyer

Model based object localization and shape estimation using electric sense on underwater robots

Keywords: Object localization, Shape estimation, Electric sensing, Modelization, Bio-robotics

.

CONTEXT OF THE STUDY

Weakly electric fish

Electric sense is a bio-inspired sensorial ability. It has been observed almost exclusively in aquatic or amphibious animals [START_REF] Bullock | Electroreception[END_REF]. Several species of fish have this capacity to sense changes in electric fields in their vicinity. Among fishes, we can distinguish two typical modes of electro-perception: some fish passively sense changes in the nearby electric fields (passive electric sense [START_REF] Bullock | Electroreception[END_REF]), while others generate their own weak electric fields and sense its distortions with their skins (active electric sense) [START_REF] Caputi | The electric image in weakly electric fish: physical images of resistive objects in gnathonemus petersii[END_REF]].

Here, we will only consider the active modality of electric sense. The fish that use active electric sense are called weakly electric fish. These fish are principally nocturnal and live in confined turbid waters of the equatorial forests [START_REF] Lissmann | The mechanism of object location in Gymnarchus Niloticus and similar fish[END_REF]. In such environments, waters are generally rich in suspended particles (turbib) and cluttered with a lot of plant roots or trees. In these harsh conditions, these fish use electric sense to avoid obstacles while navigating, to communicate with conspecifics, or to hunt their preys. Technically, electric fish perceive their environment by self-generating an electric field thanks to an electric organ located at the base of the tail, and by measuring the distortions of this field through a dense array of electro-receptors distributed over their skin. Electric sense is an omnidirectional short range sense (typical range equal to 1 fish length) that provides to the fish a lot of information about their environment. Behavioral experiments have shown that these fish can discriminate the size and the shape of objects in the dark. In particular, in [von der Emde et al., 1998], biologists demonstrated that when perceiving an object a fish first identifies its electric nature, then localizes the object and finally recognizes the shape. This behavior is astonishing as it has been shown that using electric sense the separation between localization and object geometric properties is not obvious at all. Indeed, localization and shape are both contained in the electrical measurements but are mixed together in a non linear relationship as shown in [START_REF] Rasnow | The effects of simple objects on the electric field of apteronotus[END_REF]. As an example, a small sphere close to the fish can produce the same electric response as a bigger one located further.

Robot equipped with electric sense

Electric sense is a short range omnidirectional perceptual ability ideally suited for muddy/turbib waters and confined underwater environments, i.e. in conditions in which our today underwater perception devices such as vision and sonar are useless. The first, because turbidity decreases the visibility, the second because multiple echoes make the measured signals difficult to interpret. Few artificial electric sensors inspired by weakly electric fish exist today [START_REF] Bai | Sensing capacitance of underwater objects in bio-inspired electrosense[END_REF][START_REF] Servagent | Electrolocation sensors in conducting water bio-inspired by electric fish[END_REF][START_REF] Dimble | Electrolocation-based underwater obstacle avoidance using wide-field integration methods[END_REF][START_REF] Truong | Bioinspired underwater electrolocation through adaptive system identification[END_REF]. In the paper we will use the device presented in [START_REF] Servagent | Electrolocation sensors in conducting water bio-inspired by electric fish[END_REF]. The robot is designed as a long and thin plastic cylinder (length l F = 22cm and radius r F = 1cm), on which 4 pairs of conductive electrodes are arrayed. Among reasons, the slender shape of sensor has been chosen to provide a simple model [START_REF] Boyer | Model for a sensor bio-inspired from electric fish[END_REF] and ensures the electric field to exhibit a dipolar shape mimicking that of the fish. The electric field is generated by setting a voltage U on the 2 tail electrodes (considered as the emitter) while the 6 others electrodes (considered as receivers) are grounded. The voltage (U = 10V at 22kHz) is imposed through an off-board sine wave generator. Note here that a continuous voltage would generate an undesirable electrolysis. Then, the currents that flow across the 6 receivers are measured with an ampere-meter circuit.We note I k the measure on the electrode e k that belongs to the macro-electrode K (see Fig. 1). Here, we restrict the measure to the amplitude of the electric current, the phase is not considered.

Definition of the measurement vector

The 6 measured currents are not used directly. In order to facilitate the robot control and to obtain more information about the surroundings of the robot, we use an electric measurement vector denoted M and defined as:

M = I lat δI ax , (1)
with (3)

I lat = I lat,1 I lat,2 I lat,3 = I 1 -I 2 I 3 -I 4 I 5 -I 6 , (2)
I lat (or "lateral current") represents the differential part of left and right currents. The I 0 ax stand for the "basal current", that is to say the current measured without any perturbation and δI ax ("or axial current") represents the common part of the left and right currents flowing across each i . As explained in [START_REF] Boyer | Model for a sensor bio-inspired from electric fish[END_REF], the vector of axial currents δI ax models the variations of the total resistance of the scene while its lateral counterparts I lat is proportional to the lateral incident field. From these considerations, δI ax allows to determine if the object is conductive or insulating, and I lat allows to determine whether the object is on the left or on the right side of the sensor. As a particular case, when I lat = 0, the sensor axis is necessarily aligned along the incident field.

PROBLEM STATEMENT

Problem definition

In this paper, we concentrate our efforts on the cognitive tasks of localizing and estimating the shape of an object using active electric sense. The considered scene is composed of 1 static object and 1 robot moving along a trajectory (from T (0) to T (n)) and at a distance d with respect to the object (see Fig. 2).Our robot is controlled by its forward velocity V and its angular yaw velocity Ω. In the following, we consider that we have at our disposal an accurate measure of these control inputs. In [START_REF] Khairuddin | Characterization of objects by electrosensing fish based on the first order polarization tensor[END_REF]Lionheart, 2016, Ammari et al., 2014], authors show that at leading order, the electric response of any shaped object can be described as that of an ellipsoid. Therefore, objects are described by ellipsoids and their electrical response are modeled analytically by their first order generalized polarization tensor [START_REF] Ammari | Generalized polarization tensors for shape description[END_REF]. Due to experimental constraints which prevent us from measuring the object influence only (the aquarium walls were perceived by the robot because of its smallness), all experiments were performed twice with and without object to be able to remove the effects of the aquarium. For this reason, for all experiments, the robot trajectory has been deliberately chosen as a simple straight line alongside the object according to the simple control law:

V = C and Ω = 0, with C a constant.

(4) In the global reference frame, the robot pose T (k) will be defined as: x r , y r , θ r with k ∈ [1, n].

Mathematical definition

We search for the object O (defined by its electric conductivity, localization, and shape) that disturbs the electric field generated by the robot. Moreover, we will do the following assumptions:

• the water conductivity γ is perfectly known and uniform, • the displacement of the robot is known at each currents measurement, • the object can be modeled as a prolate ellipsoid (i.e. axisymmetric about its major axis) whose length is smaller than a half of the sensor length.

Under these experimental conditions, the object O (named "real object") is entirely described by 6 parameters: its localization (x, y) and its orientation (θ) in the global frame of reference, its semi-axis (a, b) and its electric conductivity σ. The third semi-axis of the ellipsoid is equal to the second one. All these parameters are constrained in the following intervals:

• x, y ∈ [-l F , l F ] (limited by the perception range),

• θ ∈ [0, π/2], • a, b ∈ [0, l F /4] with a ≥ b, • σ ∈ {1e 5 , 1e -5 }.
Giving the fact that the parameters are constrained in small intervals and that we dispose of an analytical model able to predict the measured currents imaging a given scene (see Section 3), we propose to address the inverse problem of object localization and estimation by using a greedy algorithm that tests all possible direct models while selecting the optimal solution. It has to be noted that to obtain an accurate approximation of the ellipse parameters a discriminative robot motion along the object is essential. In the rest of the paper, we suppose that we test m candidates among which O 0 = (x 0 , y 0 , θ 0 , a 0 , b 0 , σ 0 ) is the best one.

MODELIZATION OF ELECTRIC SENSE

The analytical model for our slender robot

Assuming that we used the robot presented in the previous section (l F = 22cm and r F = 1cm), we derived an analytical model of the electrical response of an object based on its leading order dipolar tensor [START_REF] Boyer | Model for a sensor bio-inspired from electric fish[END_REF]. The model considers that the object is small enough to assume the electric field as uniform on its domain.Such a model has given good recognition results in simulation with such small objects [START_REF] Lanneau | Object shape recognition using electric sense and ellipsoid's polarization tensor[END_REF]. For a given scene, the currents measured on the receiving electrodes is given by:

I lat (T (k), O x ) = 1 4π .P ⊥ .H.R θs .P.R t θs .G t .C 0 .U, (5)
δI ax (T (k), O x ) = 1 4πγ
.C 0 .G.R θs .P.R t θs .G t .C 0 .U, (6) where I lat , δI ax are vectors of currents 4 × 1, C 0 is a 4 × 4 matrix encoding the robot morphology and conductivity, P ⊥ is a 4×4 diagonal matrix depending on the polarization of the electrodes, R θs is a rotation matrix depending on the angle between the sensor and the object θ s = θ -θ r , P a 3 × 3 diagonal matrix modeling the ellipsoid object electric response (called tensor), G and H are 4×3 matrices encoding the distance between the object center and the electrodes, and U is the 4 × 1 polarization vector encoding the voltage imposed to the emitting electrodes with respect to the receiving ones. In [START_REF] Boyer | Model for a sensor bio-inspired from electric fish[END_REF], P ⊥ , C 0 , U, G, H are defined by:

P ⊥ =    5.29e -4 0 0 0 0 6.74e -4 0 0 0 0 6.95e -4 0 0 0 0 5.52e -4    , C 0 = γ.
I lat,1 2) I lat,2 , 3) I lat,3 4) δI ax,1 5) δI ax,2 , 6) δI ax,3 , 7) Legend 8) Scores 9) Result. U =    0 0 0 U 0    , G =       x c -x 1 r 3 1 y c r 3 1 0 . . . . . . . . . x c -x 4 r 3 4 y c r 3 4 0       , H =        3y c (x c -x 1 ) r 5 1 2y 2 c -(x c -x 1 ) 2 r 5 1 0 . . . . . . . . . 3y c (x c -x 4 ) r 5 4 2y 2 c -(x c -x 4 ) 2 r 5 4 0        , with r k=1...4 = (x c -x k ) 2 + (y c ) 2 ,
x k the axial coordinates of the electrode e k and x c , y c the object center position in the sensor frame (see Fig. 2). [START_REF] Ammari | Generalized polarization tensors for shape description[END_REF] As shown in [START_REF] Ammari | Generalized polarization tensors for shape description[END_REF], at the leading order the electric response of an object can be modeled by the following tensor:

Expression of the ellipse tensor

P = V. λ 1 (η) 0 0 0 λ 2 (η) 0 0 0 λ 2 (η)
, with V = 4/3.π.a.b 2 the object's volume, η = a/b its aspect ratio, and λ 1 , λ 2 defined as:

• (λ 1 , λ 2 ) = (1/A, 1/B) for conductive objects, • (λ 1 , λ 2 ) = (1/(A-1), 1/(B-1)) for insulating objects.
where we introduced the elliptic integrals A, B as:

         A(η) = η -2 +∞ 1 1 t 2 (t 2 -1 + η -2 ) dt, B(η) = η -2 +∞ 1 1 (t 2 -1 + η -2 ) 2 dt. (7)
For a sphere (i.e. a = b), which is a particular case of the ellipse we have λ 1 = λ 2 = χ.a 3 with a the radius of the sphere, and χ a contrast factor. χ = 1 for conducting objects and χ = -1/2 for insulating objects. Moreover, as the sphere is isotropic, R and R t can be removed in Eq. 5,6.

On Fig. 3 we present a first illustration of localization and shape estimation using this ellipse tensor for an ellipsoid of size 33 × 16mm with θ s = 0. This figure shows the good approximation of the 6 M i in comparison with the real experimental measurements. The score evolutions for the real object (black) compare to the estimated one (red) are shown on subplot 8 and the result is shown on subplot 9.

OBJECT LOCALIZATION AND ESTIMATION

Developed approach

Our method can be structured into 2 stages. First, the recognition of the electric conductivity of the object and its side with respect to the robot is performed. Second, the estimation of its localization and of its geometric properties is done.

Stage 1: Detect the object, find σ 0 and sensor side on which the object is. Since the sensor goes forward, it discovers its environment with its front head electrodes. While δI ax,1 = 0 the robot do not perceive the object so we do not use the measurements. When δI ax,1 = 0, depending on its sign we extract the electric conductivity of the object σ 0 . Then from σ 0 , the side of the object is found depending on the sign of I lat,1 as shown in [START_REF] Lebastard | Reactive underwater object inspection based on artificial electric sense[END_REF] where the following scenarios were listed:

• δI ax,1 , I lat,1 > 0 it is conducting on the left,

• δI ax,1 , I lat,1 < 0 it is insulating on the left,

• δI ax,1 > 0, I lat,1 < 0 it is conducting on the right,

• δI ax,1 < 0, I lat,1 < 0 it is insulating on the right,

• δI ax,1 > 0, I lat,1 = 0 it is conducting facing the robot,

• δI ax,1 < 0, I lat,1 = 0 it is insulating facing the robot.

Stage 2: Estimate the 5 unknown parameters: x 0 , y 0 , θ 0 , a 0 and b 0 . Giving the object electric conductivity and its side with respect to the robot, while we perceive the object we continue the localization and the shape estimation. This estimation requires from the user the 3 discretization parameters xy , θ , ab that divide the research space defined in Section 2.2. Then, for O m an object candidate, T (k) a pose of the robot (x r,k , y r,k , θ r,k ) on the trajectory T , our model computes a vector of 6 estimated currents through the function f (see Eq. 8). Supposing n to be the number of points along the trajectory, we have ∀k ∈

[1, n], f (T (k), O m ) = M (k) = I lat (T (k), O m ) δI ax (T (k), O m ) . (8)
Then, we defined an evaluation function g as:

g(M, T, O m ) = n k=1 6 i=1 |M i (k) -f i (T (k), O m )| |M i (k)| , (9)
which represents a scalar equivalent to the cumulative residual error between the measured currents M and the estimated currents f (for an given object O m ) along the whole trajectory. By estimating g(M, T, O m ) for all O m , the best object estimation disturbating the electric field is the object O 0 that minimizes g according to:

O 0 = argmin Om g(M, T, O m ). ( 10
)

From an exhaustive testing to an optimized approach

The algorithm introduced in Section 4 is described in Alg. 1. It consists in systematically evaluating all the solutions, input : n number of point of the trajectory input : T the robot trajectory input : M the measured currents along the trajectory input : xy , θ , ab the discretization parameters output: x 0 , y 0 , θ 0 , a 0 , b 0 , σ 0 characterizing an ellipsoid.

/* Stage 1: Detect object and find σ 0 , side */ t ← 0; while M (t, 4) = 0 do t+ = 1;

end if M t,4 > 0 then σ ← 1e5; if I t,1 > 0 then side = -l F ; else side = l F ; end else σ ← 1e -5; if I t,1 > 0 then side = l F ; else side = -l F ; end end /* Stage 2: Find x 0 , y 0 , θ 0 ,
| > 0)&&(t < n) do for k ← 1 to m do for i ← 1 to 6 do O k .score+ = |M (t, i) -f i (T t , O k .params))|; end end end /* Select the best candidate */ O m = sort least to greatest score(O m ); return O 0 .params;
Algorithm 1: Global algorithm: Ellipse localization and shape estimation from a moving underwater robot. and then selecting the set of parameters that minimizes the cumulative residual error between the model and the experimental measures. As the algorithm is designed, it gives in any case a solution after browsing all the parameter space and testing all parameters combinations (the evaluation of a single candidate includes the computation of a tensor followed by a multiplication of matrices with a maximum dimension of 4). However, the exhaustive testing of all candidates is time consuming. The complexity only depends on the space discretization parameters xy , θ , ab , as their intervals are fixed by the sensor length: [-l F , l F ] for x 0 , y 0 , [0, l F /4] for a 0 , b 0 , and [0; π/2] for θ s . To compute the complexity we will suppose θ = 10 • (i.e. a constant), xy = ab , and we will call X the number of samples in the interval [0, l F /4] (i.e. X = l F / xy ). Then, the complexity can be expressed as: nb = (8X).(4X).( 10).(C 2 X ) = (8X).(4X).( 10).(X!/(2!(X -2)!) = (8X).(4X).( 10).(X.(X -1).

(X -2)!)/(2. (X -2)!) = (160)X 4 -(160)X 3 where C 2 X is the binomial coefficient estimating all the combinations of the two semi-axis parameters without repetition since a 0 ≥ b 0 . For example, with our robot, l F = 22cm with xy = ab = 1mm (which seems to be a desirable accuracy) the number of candidates reaches 1.29e 9 . To reduce the number of operations and improve the performance of our algorithm, we developed an optimized approach that drastically reduces the number of operations while keeping the same performance. This optimized approach is an iterative procedure that consists of 3 stages. We begin by localizing and estimating a sphere (only 3 parameters: x 0 , y 0 and radius a 0 ). Then, from the first approximation of the sphere, we localize and roughly estimate an ellipse. Finally, from this approximated ellipse we reduce again all the intervals and increase the accuracy to the desired one. Assuming again θ = 10 • , and keeping the definition of X. The parameters intervals and accuracies at each of the steps are defined such as:

• x 0 ∈ [-l F , l F ],
with a grid of 4X, y 0 ∈ [0, l F ], with a grid of 2X and a 0 ∈ [0, l F /4], with a grid of X/4 (rough approximation of a sphere), • x 0 , y 0 ∈ [-l F /2, l F /2], with a grid 2X and a 0 , b 0 ∈ [-l F /8, l F /8], with a grid of X/2 (rough approximation of an ellipse), • x 0 , y 0 ∈ [-l F /8, l F /8], with a grid X and a 0 , b 0 ∈ [-l F /16, l F /16], with a grid of X/2 (accurate approximation of an ellipse).

The complexity at each stage is computed as:

• nb 1 = (4X).(2X).(X/4) -Rough sphere,

• nb 2 = (2X).(2X).( 4).(C 2 X/2 ) -Rough ellipse, • nb 3 = (X).(X).( 10).(C 2 X/2 ) -Accurate ellipse, Then, the overall complexity can be written as Eq. 4.2.

nb = nb 1 + nb 2 + nb 3 = 2X 4 + 8X.(X/2).(X/2 -1) + 5X 2 (X/2).(X/2 -1) = 2X 4 + 2X.(X).(X -2) + 5/4X 2 (X).(X -2) = (13/4)X 4 -(9/2)X 3
This strategy reduces dramatically the number of candidates evaluated as this number is reduced to 5.5e 7 for an accuracy of 1mm.

EXPERIMENTAL RESULTS

Presentation of experiments

To validate our object localization and estimation method we performed a large set of experiments in a 1m 3 tank with the slender shape robot presented in the Section 1. All experiments were performed under the same conditions: a straight line trajectory of about 40cm, and one object situated in the middle of the trajectory at a distance d (see Fig. 2). In these conditions, we experiment on 4 objects: 2 ellipsoids shapes (33×16mm and 27×18mm) with 2 different electric conductivity each (aluminum that is conductive and plastic that is insulating), 4 different angles with respect to the trajectory (0 • , 30 • , 60 • , 90 • ), and 7 different distances (d = 50, 60, 70, 80, 90, 100, 110mm). In total, we performed 108 experiments instead of 112, because for each object at 50mm and 90 • the robot was colliding the object. All results were obtained post-processing the raw experimental data with c Matlab on a desktop computer equipped with an c Intel Core i5 3.3 GHz CPU and 8 Go of RAM. Globally, the average processing time was around 10 min for an average experiment time of 3 min.

5.2

Results with the ellipsoid tensor [START_REF] Ammari | Generalized polarization tensors for shape description[END_REF] In Tab. 1, we present the mean results over the 108 experiments using the dipolar tensor of [START_REF] Ammari | Generalized polarization tensors for shape description[END_REF]. This table summarizes errors on both x and y-axis in mm, on the orientation in degrees, and on the shape in percentage of shape error. This shape error is defined in Eq. 11.

E s = (a -a 0 ) 2 + (b -b 0 ) 2 (a 2 + b 2 ) . ( 11
)
Moreover, we add the average processing time in minutes for a single experiment. Globally, errors are smaller on the localization than on the shape and the orientation. Moreover, errors are more important on insulating objects than on conductive objects. These errors are mainly due to the violation of an assumption on which is based our analytical model. In [START_REF] Boyer | Model for a sensor bio-inspired from electric fish[END_REF] we suppose that the object was small enough (a 0 ≤ r F with r F = 1cm) to consider the electric field uniform on its domain. But, in this paper, we experiment on objects bigger than 1cm as their length is about 6.r F (our biggest object is 6.6cm long). In fact, this unconsidered assumption is not problematic when the object and the sensor are aligned, i.e. θ s = 0 as the electric field is globally uniform over the object in this case (see Fig. 6.a). Unfortunately, when the angle θ s increases towards π/2, the non uniformity of the electric field is more visible (see Fig. 6.b) and causes an underestimation of the amplitude of the estimated I lat and δI ax . This can be seen on Fig. 4 which reproduces experiments shown Fig. 3 after increasing θ s to π/2. Because localization and shape are intricately linked in the electric measurements, this modelization error unavoidably leads to important errors on the object localization and on the shape estimation. It has to be noted that this phenomenon appears for both conducting and insulating object and that it increases with θ s . To overcome this underestimation of the amplitude, we proposed a corrected tensor suited to bigger objects. This tensor models objects by 3 dipoles instead of one (1 ellipse and 2 spheres, see Fig. 7) to take into account the non uniformity of the electric field on the object (Fig. 6. b).

A new tensor for our analytical model

This new tensor consists in computing a 2 spheres electric responses superimposed to the ellipse electric response (see Fig. 7). The 2 spheres parameters (localization

x 1 , y 1 , x 2 , y 2 and 2 radius r 1 , r 2 ) depends on the ellipse localization and size (x, y, a, b) and are obtained such as:

∀i ∈ [1, 2], x i = x + ((-1) i .0.8. (a 2 + b 2 )),

y i = y, r i = b + a 10 . (12)
The estimated currents are modeled as a weighted sum between the ellipse currents (M ell ) and the 2 spheres currents (M sph1 and M sph2 ) using Eq. 13.

M = (1 -µ).M ell + µ 2 .(M sph1 + M sph2 ). ( 13
)
The weights depend on the angle θ s and the electric conductivity of the object according to Eq. 14 and 15.

The weights are computed as:

µ      [0; π 2 ] → [0; 1] θ s → 1 - 2 π .|θ s | for conducting object. ( 14
) µ      [0; π 2 ] → [0; 1 2 ] θ s → 1 2 - 1 π .|θ s | for insulating object. ( 15
)
It has to be noted that the parameters of the corrected tensor have been identified empirically on few experiments in order to complement the underestimated signal due to big objects. Using this correction, in the conditions of Fig. 4, it can be seen on Fig. 5 that the estimated currents I lat and δI ax and the measurements are closer to each other. To clearly show the advantage of such a correction, we compare in Tab. 1 the mean results for the shape, the orientation and the localization obtained along our 108 experiments described at the beginning of this section. Tab. 1 shows that the results are improved, as we gain in average: 2.2% on the shape errors, 3 -4mm on the localization and 2 • on the orientation. This global improvement is even more visible on the shape error histogram presented on Fig. 9.a. Tab. 1 also shows that the results are a bit better for conducting objects than insulating objects but one can see that our corrected tensor is more beneficial for insulting object with an improvement of 3.5% on the shape instead of 1%. Now, to clearly illustrate the effectiveness of localization and shape estimation algorithm we show on Fig. 8 the results on 54 experiments that is to say half of the results: real objects are plotted in red, and the corresponding estimated objects are plotted in blue. On this figure, it can be seen that for both objects, results are getting worst (localization and position) as the distance increases since in this case the signal-to-noise ratio decreases. The errors that remain on the localization and on the shape are due to the fact that both are linked in the measurements. In fact, multiple couples (object/localization) produce exactly the same measurements, thus, the algorithm sometimes cannot distinguish between a small object situated close from a bigger situated further. As a consequence, to complement these results and evaluate our algorithm and model, we performed 2 additional experiments that consist in estimating the shape when the localization x 0 , y 0 is supposed known (resp. estimating the localization assuming the shape a 0 , b 0 is known). Experiments for the 2 particular cases are presented on Fig. 10.a and 10.b. Remarkably, results are really good in both cases even for long distances when the signal-to-noise ratio is very low. This is trivially explained by the fact that we estimate 3 parameters instead of 5 which reduces a lot the number of candidates. To see globally this result on the 108 experiments we show on Fig. 9.b the shape error histogram with the nominal model and the new one estimating the shape with the knowledge of the localization. This histogram shows again a global decrease of 3% of the shape errors which demonstrates the improvements brought by the new model.

CONCLUSION AND FUTURE WORK

We demonstrate in this paper that an underwater robot equipped with electric sense can estimate at the same time the pose and the geometric properties of an object while navigating in its surroundings. We built an algorithm that: first detects the object, then finds its electric nature and the side of the object with respect to the robot and, finally, estimates the pose and the object shape. Our algorithm has been tested on 4 different big objects with 2 different electric conductivity, 4 different orientations and at 7 distances. We obtain as average 18% of errors on the shape, 25 degrees on the orientation and less than 1cm over 108 experiments. The observed errors are mainly due to 2 factors that will be explored in future work. First inaccuracies of the model remain, in particular we still neglect the polarization of the sensor by the object and this phenomena has to be taken into account for the big objects we used. Secondly, the trajectory we chose is sometimes not discriminative enough to estimate with accuracy the object localization and geometric properties. As an example in nature, the electric fish approach the object and turn around it to estimate its pose and shape. In our future work we plan to enrich the sensing information by following more complex trajectories. For example, we will try to mimic the trajectory of the fish by navigating towards the object and turning around it. This will be achieved by using reactive behavior based on electric sense recently presented in [START_REF] Lebastard | Reactive underwater object inspection based on artificial electric sense[END_REF].

Fig. 1 .

 1 Fig. 1. Our robot equipped with 4 pairs of electrodes.

Fig. 2 .

 2 Fig. 2. Experimental set up with our notations.

Fig. 3 .

 3 Fig.3. Experiment using Ammari tensor with a conductive ellipsoid 33 × 16mm with θ s = 0 • , d = 70mm. For all plots, experimental measurements are in blue, model measurements with the real object are in dotted black, and model measurements with the estimated object are in red. 1) I lat,1 2) I lat,2 , 3) I lat,3 4) δI ax,1 5) δI ax,2 , 6) δI ax,3 , 7) Legend 8) Scores 9) Result.

  for x ← -l F to l F by xy do for y ← 0 to side by xy do for θ ← 0 to π by θ do for a, b ← 0 to l F 4 by ab do m+ = 1; O m .score = 0; O m .params = set params (x, y, θ, a, b, σ 0 ); O m .tensor = eval tensor (θ, a, b, σ 0 ); end end end end /* Evaluate all the candidates */ while (|M(t, 4)

Fig. 4 .Fig. 5 .

 45 Fig. 4. Experiment using Ammari tensor with a conductive ellipsoid 33 × 16mm with θ s = 90 • , d = 70mm. 1) I lat,1 2) I lat,2 , 3) I lat,3 4) δI ax,1 5) δI ax,2 , 6) δI ax,3 , 7) Legend 8) Scores 9) Result.

Fig. 6 .

 6 Fig. 6. Robot electric field applied on an object depending on its orientation. Up) θ s = 0 • , Down) θ s = 90 • . The black vectors show the orientation of the electric field in the 2 focus and in the center of an ellipse.

Fig. 7 .

 7 Fig. 7. The 3 objects defining the new tensor (right) compared to the nominal tensor of the ellipse (left).

Fig. 8 .

 8 Fig. 8. Object localization and estimation with the new tensor for 2 different object, 4 orientations, 7 distances.

  Fig. 9. Histograms of the shape error for the 108 experiments: Nominal model (blue), New model (red).

  Fig. 10. Shape estimation giving the localization (left) and object localization giving the shape (right). Both results are obtained using the new tensor.

Table 2 .

 2 Shape estimation giving the localization and object localization giving the shape.Mean on 108 experiments.

This project is supported by the European Unions Horizon 2020 research and innovation program under the grant agreement No 640967. Project subCULTron: http://www.subcultron.eu.


Bruno Quesson

Mathilde Merle

Max Köhler

Charles Mougenot

Sebastien Roujol

Baudouin Denis De Senneville

Chrit Moonen

A method for MRI guidance of intercostal high intensity focused ultrasound ablation in the liver

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

  Figure 1


Quentin Scanvion
email: quentin.scanvion@gmail.com

Thierry Perrez

François Cassim

Olivier Outteryck

Aurélia Lanteri

Pierre-Yves Hatron

Marc Lambert

Sandrine Morell-Dubois

Thierry Perez

Pierre- Yves Hatron

Pierre-Yves Hatron

Sandrine Morell

Pierre-Yves Hatron

Neuralgic amyotrophy triggered by hepatitis E virus: a particular phenotype

Keywords: Brachial neuritis, Parsonage-Turner syndrome, peripheral neuropathies, HEV, phrenic paresis

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

Introduction

The neuralgic amyotrophy syndrome (NA) is an acute multiple mononeuropathy predominantly affecting the brachial plexus described around 1880. An intensive neuropathic pain in the scapular area or upper arm, often exacerbated at night and when arm is at rest, is generally followed by multifocal paresis (stereotypically including the serratus or infrascapular muscles) and/or hypoesthesia, then an amyotrophy. However, several distinct subtypes can occur ─ as with phrenic nerves involvement ─ and generate the NA phenotypic variability [START_REF] Van Eijk | Neuralgic amyotrophy: An update on diagnosis, pathophysiology, and treatment[END_REF].

In fact, there are 2 forms of NA: (i) the hereditary neuralgic amyotrophy (HNA) associated with a mutation in the SEPT9 gene, on chromosome 17q25, in 55% of families [START_REF] Van Alfen | Clinical and pathophysiological concepts of neuralgic amyotrophy[END_REF]; and (ii) the idiopathic form (INA), also known as brachial neuritis or Parsonage-Turner syndrome (the eponymous names of the authors of the first case series, 1948). A 1/1000/year incidence, previously underestimated, has recently been reported [START_REF] Van Alfen | Incidence of neuralgic amyotrophy (Parsonage Turner syndrome) in a primary care setting -A prospective cohort study[END_REF].

The exact etiopathogenesis is unknown, but post infectious immune reactions are frequent (HIV, Cytomegalovirus, Epstein-Barr virus, parvovirus B19, Campylobacter, Yersinia, Borrelia) [START_REF] Van Alfen | Clinical and pathophysiological concepts of neuralgic amyotrophy[END_REF]. Recurrences happen in 30-40% patients with INA [4,5] mostly athletic, therefore mechanical factors and an individual genetic susceptibility are assumed [START_REF] Van Eijk | Neuralgic amyotrophy: An update on diagnosis, pathophysiology, and treatment[END_REF][START_REF] Dartevel | Hepatitis E and neuralgic amyotrophy: Five cases and review of literature[END_REF].

Hepatitis E virus (HEV) infection could trigger many NA. HEV has been found in 10% NA patients [START_REF] Van Eijk | Neuralgic amyotrophy and hepatitis E virus infection[END_REF] and it is the leading cause of hepatitis worldwide (50%) [START_REF] Dartevel | Hepatitis E and neuralgic amyotrophy: Five cases and review of literature[END_REF]. It also accounts for an autochthonous emerging disease, spreading through industrialized countries (genotype 3, 4) [START_REF] Dalton | Hepatitis E virus: time to change the textbooks[END_REF]. Animals account for the primary reservoir, for instance 65% French pigs are HEV+ [START_REF] Rose | High prevalence of hepatitis E virus in French domestic pigs[END_REF]. Although usually asymptomatic, genotype 3-HEV extrahepatic and neurological involvement has been reported more frequently since 2000 [START_REF] Dalton | Hepatitis E virus: time to change the textbooks[END_REF][START_REF] Kamar | Extrahepatic manifestations of hepatitis E virus[END_REF].

The study seeks to identify whether specific clinical features characterize patients bearing INA triggered by HEV, implying the arise of a specific phenotype.

Methods

The case of a man referred to the university hospital of Lille, France with INA induced by hepatitis E associated with a bilateral phrenic nerves involvement has been thoroughly described. No date or language restriction was used for A] as opposed to B] where search was restricted to English articles on Pubmed from 1999 to 2016. After deletion of duplicates, we grouped all case-reports of INA induced by hepatitis E (A) and all cases reports of INA with phrenic paresis (B), including our patient.

The following features (i.e. gender ratio, mean age, prevalence of phrenic paresis, bilateralism and average transaminases in group A; prevalence of bilateral phrenic paresis in group B; and pain duration and clinical improvement (with the median of follow-up) in both groups) were analyzed with respect to those of the largest cohort of NA available, described by Van Alfen and Van Engelen in 2006 [4]. Quantitative variables are expressed as mean ± standard deviation (SD). Qualitative variables are expressed as frequencies and percentages. Gender ratio and proportion of bilateral and asymmetric NA were compared using Chi-square test.

Fisher Exact test was used to compare proportion of painless and proportion of phrenic paresis. The Student's ttest was used to compare mean age. In particular for the historic cohort, SD was estimated from the range (range/4) as such data was unavailable. Statistical testing was done at the two-tailed  level of 0.05 and data were analyzed using the SAS software package, release 9.3 (SAS Institute, Cary, NC).

Results

Case report

A 42-year-old athletic (marathon, tennis) man rushed to emergency on October 1 st , 2015 at night, stating a brutal, non-traumatic and very intense (numeric rating scale 10/10) para-spinal cervicothoracic, predominantly right back pain, associated with vagal symptoms (nausea, vomiting). His medical records revealed a 20-year-old malaria (Plasmodium falciparum) episode and a penicillin allergy (skin rash) without any relevant neurological family history.

Chest emergencies, including dissecting aorta were excluded. Neither pleural or parenchymal injury nor osteolytic or sclerotic foci were revealed. Intravenous morphine was administered to ease the pain. Soon after (<24h), the patient described transient paresthesiae in the right arm, hypoesthesiae of the right thenar eminence and inside of the left knee (territory of the saphenous nerve). Decreased deep tendon reflexes of the upper limbs and moderate high dysphagia were reported with neither amyotrophy nor pyramidal irritation. Brain MRangiography and cerebrospinal fluid analysis excluded myelitis or subarachnoid hemorrhage. EMG did not reveal any signs of demyelination or neuropathy on any limb. However, the serratus anterior, infraspinatus, flexor pollicis longus and pronator quadratus (high-yield muscles for NA) had not been EMG-explored at that time.

Biological analysis depicted a sharp hepatic cytolysis and a moderate cholestasis, with no fulminant hepatitis.

Lactate dehydrogenase (350 IU/L) and creatine kinase (260 IU/L) were also increased. Blood and PCR analyses excluded HBV, HCV, HIV, VZV, CMV, HSV, EBV, parvovirus B19, HHV6, HHV8, Borrelia, syphilis, ANCA, antinuclear antibodies, antibodies associated with autoimmune hepatitis as well as an inflammatory syndrome, thrombocytopenia (125,000/mm 3 ), lymphopenia (500/mm 3 ) with 10% mononucleosis-like cells. IgM were increased to 4.06g/L, along with a monoclonal IgM- spike on the blood protein electrophoresis.

Transferred to the internal medicine department, the patient narrated he had spent three weeks in Corsica in August. Further analyses were run to test for HAV, Hantavirus, Chikungunya, Leptospirosis, Enterovirus and HEV. HEV serology (IgM and IgG), confirmed by positive blood RT-PCR (32.100 replicas/mL) confirmed acute hepatitis E. Genotyping attested a subtype 3f. The infection was recognized via the Corsican sausage (Figatelli) consumption 5 to 8 weeks prior to the symptoms.

The intermittent spinal pain lasted for a week with a subsequent spontaneous, progressive improvement of the pain along with the liver parameters. The MRI revealed isolated T2 STIR hypersignal of the right infraspinatus, with no signs of muscular atrophy, suggesting myositis.

In November, the patient had no pain, no amyotrophy, no more paresthesia and biology revealed the normalization of lymphocytes, IgM, electrophoresis, the inversion IgM/IgG on HEV serology and the decrease of RT-PCR (1000 replicas/mm 3 ).

In January, the patient referred orthopnea, paradoxical breathing and dyspnea on exertion with a gradual appearance since he had been dismissed in November. Respiratory function tests revealed forced vital capacity, 5.15L (i.e. 92% predicted), with 50% decrease in the supine position; total lung capacity, 7.42L (94% predicted); forced expiratory volume in 1 second, 4.3L (83% predicted). The global inspiratory muscle strength was low in voluntary moves, with a PImax (maximum inspiratory pressure) of 63 cmH 2 0 (predicted=105 cmH 2 0). His sniff nasal pressure was underperforming. The phrenic EMG confirmed the left paralysis (no motor potential after stimulation) and the very severe right paresis.

A final diagnosis of HEV-induced INA with severe bilateral phrenic paresis was established. The patient was offered physiotherapy and respiratory rehabilitation, with no treatment for hepatitis as RT-PCR was negative at that time.

At the 12-month-follow up visit, there were no motor deficits, no amyotrophy, and full recovery of the tendon reflexes. Nonetheless, hypoesthesia persisted, and the phrenic paresis had clinically moderately Neuralgic amyotrophy triggered by HEV. DOI: 10.1007/s00415-017-8433-z 6 improved. Supine fall of forced vital capacity improved at 25% and maximal inspiratory pressure was normalized at 123 cmH 2 0. The left phrenic paralysis decreased.

Review of the relevant literature

Thirty-eight cases (group A) relating HEV-induced INA and 41 cases (group B) reporting INA with phrenic paresis were found in the literature. Sources and characteristics are detailed in the table 1 [START_REF] Dartevel | Hepatitis E and neuralgic amyotrophy: Five cases and review of literature[END_REF][START_REF] Van Eijk | Neuralgic amyotrophy and hepatitis E virus infection[END_REF] and appendix table e-1, respectively. The features of each group are described in table 2.

The patients with HEV-INA were 46 years old (46±9.8), with statistically more males (87.2%) than INA of all etiologies (68.3%) and the syndrome is more frequently bilateral and asymmetric (86.8% vs 27.7%). A phrenic paresis is more frequently reported: 18.0% (7/39), vs 6.6%. The clinical improvement is poor: remission only occurs in 15.6% cases over a 10 month follow-up. The cytolysis predominates on the ALT with a 2/1 ratio.

The patients affected by INA with phrenic paresis were 53 years old (52.9±13.0). The absence of initial pain (23.8%) emerges as a peculiar feature of this group.

Discussion

Our case report illustrates the phenotypic variability of NA and underlines the difficulty of a prompt diagnosis. The two essential criteria are: a typical pain (acute onset of very severe and analgesic-resistant pain in the shoulder area) and a multifocal sensory and/or motor neuropathy [START_REF] Van Eijk | Neuralgic amyotrophy: An update on diagnosis, pathophysiology, and treatment[END_REF]. Other signs are optional. Amyotrophy depends on motor neuropathic impairment. Electrophysiological support findings are not fundamental because negative results are often due to sampling error [START_REF] Van Eijk | Neuralgic amyotrophy: An update on diagnosis, pathophysiology, and treatment[END_REF]: it is essentially a clinical diagnosis. In our case, the NA initial presentation combines several subtypes: a pure sensory form (patient did not show motor weakness or amyotrophy of limbs); a lumbosacral involvement (8.2% of cases, and only <6.7% for sensory signs [4]) with anesthesia of the territory of the left saphenous nerve; a distal reached (thenar hypoesthesia); an involvement of the cranial nerves (<2.6% of cases [4]) with a brief high dysphagia. In addition, when the phrenic involvement occurs, it exacerbates the severity of the syndrome.

The analysis of the literature review has yielded to a particular phenotype for the HEV-induced INA that is described by: a prevalence of men, moderately older, a bilateral and very asymmetric form, and a frequent phrenic nerve involvement.

Our findings are in agreement with a previous suggestion reported by Van Alfen et al. derived from their large cohort: "…there were the grossly elevated liver enzymes for which no cause could be identified in 10, usually male, patients, most of whom seemed to have a clinical subtype of NA with extensive bilateral brachial plexus and phrenic nerve involvement. This subtype may be triggered by a specific but unknown antecedent infection that also causes a transient hepatitis" [4].

At that time, HEV had not yet been identified (and the HEV-INA prevalence was unknown.) In the past decade, HEV has been reported to (i) trigger neurological disorders including Guillain-Barre syndrome, meningoencephalitis and INA and (ii) to be implicated in extra hepatic manifestations [START_REF] Kamar | Extrahepatic manifestations of hepatitis E virus[END_REF][START_REF] Perrin | Neurologic disorders in immunocompetent patients with autochthonous acute hepatitis E[END_REF]. In fact, research on HEV has recently evidenced the HEV neurotropism by proving its capability of replicating, assembling and entering specific nervous system-derived cells [START_REF] Drave | Extra-hepatic replication and infection of hepatitis E virus in neuronal-derived cells[END_REF]. Another hypothesis is that HEV may promote INA in predisposed people [START_REF] Dartevel | Hepatitis E and neuralgic amyotrophy: Five cases and review of literature[END_REF]. Recent guidelines recommend to systematically search HEV for all INA [START_REF] Kamar | Extrahepatic manifestations of hepatitis E virus[END_REF].

Overall, this literature review suggests that the phrenic paresis is an important symptom in the phenotype HEV-INA. This preponderance may be underestimated because some cases with phrenic paresis within the historic cohort could have been triggered by HEV. Ten patients were described with hepatitis [4] but HEV like etiology was unknown in 2006. In contrast, the increased prevalence of phrenic involvement in this subgroup must be moderated. Now Van Eijk et al. estimate 10% of phrenic reached in their current cohort (>1,200 patients), but they don't precise the proportion of HEV or HNA [START_REF] Van Eijk | Neuralgic amyotrophy: An update on diagnosis, pathophysiology, and treatment[END_REF]. This impairment not uncommon in HNA family (14%) [4], could overestimate the prevalence of our group A (cases 30 and 33 suggest a possible HNA, but no mutation of the SEPT9 gene were found [START_REF] Dartevel | Hepatitis E and neuralgic amyotrophy: Five cases and review of literature[END_REF]). The neuralgic amyotrophy physiopathology associates biomechanical factors (on the plexus brachial) to an immune trigger [START_REF] Van Eijk | Neuralgic amyotrophy: An update on diagnosis, pathophysiology, and treatment[END_REF][START_REF] Van Alfen | Clinical and pathophysiological concepts of neuralgic amyotrophy[END_REF]. The authors hypothesize that the strong neurotropism of HEV may induce a more diffuse form of neurologic disorders. The phrenic involvement is likely the most symptomatic, therefore the most often diagnosed and mentioned, but not the only one (like the transient high dysphagia of our patient). Furthermore, NA is a disease of young adults (35-40 years on average) [4,[START_REF] Cruz-Martínez | Neuralgic amyotrophy: variable expression in 40 patients[END_REF] although those with INA triggered by HEV tend to be on the older side of such young adulthood. On the other hand, hepatitis E more frequently affects older patients (median age 64 years for Woolson et al. [START_REF] Woolson | Extra-hepatic manifestations of autochthonous hepatitis E infection[END_REF]) however, those patients with neurological manifestations tend to be among the younger ones. This curiously overlap seems to correspond to a population that interesting matches most of the reported HEV-INA cases.

The general characteristics of HEV-INA phenotype appear independent of phrenic involvement. When all INA with phrenic involvement are considered, the hereby reviewed cases are older subjects, as previously reported in two small series [START_REF] Hughes | Long-term recovery of diaphragm strength in neuralgic amyotrophy[END_REF][START_REF] Podnar | Nosology of idiopathic phrenic neuropathies[END_REF] of 19 patients (mean age 59.5 and 61 years, respectively). As described in our case, a symptomatic phrenic involvement may appear much later than the hepatitis E. Moreover, pain was not reported in many cases of group B. We can suggest that some asymptomatic hepatitis E has gone misdiagnosed.

A better understanding of this NA phenotype seems crucial to early diagnose HEV-INA. In fact, in our case the clinical picture has been confused with HEV-induced myositis, as it has previously been described [START_REF] Mengel | Hepatitis E-induced severe myositis[END_REF].

Recovery of strength is rather slow (several years) and often incomplete as pain, disabilities, little residual paresis or intolerance to effort are still present in 3 patients out of 5 [START_REF] Van Alfen | Clinical and pathophysiological concepts of neuralgic amyotrophy[END_REF][START_REF] Cup | Residual complaints after neuralgic amyotrophy[END_REF]. A worse prognosis is linked to phrenic involvement [START_REF] Van Eijk | Neuralgic amyotrophy: An update on diagnosis, pathophysiology, and treatment[END_REF][START_REF] Hughes | Long-term recovery of diaphragm strength in neuralgic amyotrophy[END_REF][START_REF] Tsairis | Natural history of brachial plexus neuropathy. Report on 99 patients[END_REF], which accounts for a severity factor affecting sport and work abilities and leading to a worse quality of life and sleep. It is therefore necessary to be aware of and search for phrenic involvement, especially in cases of hepatitis E.

There is no current treatment for NA. A combination of nonsteroidal anti-inflammatory drugs and opiate is the best symptomatic option to treat initial neuropathic pain [4]. It is uncertain whether a short corticotherapy at NA onset decrease the symptom duration [START_REF] Van Eijk | Evaluation of prednisolone treatment in the acute phase of neuralgic amyotrophy: an observational study[END_REF] and it has been reported as non-effective on phrenic damage [START_REF] Tsairis | Natural history of brachial plexus neuropathy. Report on 99 patients[END_REF].

Acute hepatitis E doesn't usually require a specific treatment. In general, immunocompromised patients are to be treated because of the risk of chronic hepatitis. In particular, ribavirin is effective on acute or chronic HEV [START_REF] Pischke | Ribavirin treatment of acute and chronic hepatitis E: a single-centre experience[END_REF]. However, we are still debating whether patients with extra-hepatic manifestations should be treated for HEV. Strong arguments in favor of treating male patients are poor clinical improvement of INA and rather great proportion of phrenic nerve damage. As a continuous effort to (i) decrease neurological symptoms, (ii) improve quality of life and (iii) decrease overall sick leave, we feel that stronger therapeutic studies for INA triggered by HEV and guidelines are required. In fact, only empirical treatments are currently available in the literature.

The multiple limitations of a retrospective analysis of cases of independent publications must be considered (missing information and misclassification, publication bias and risk of non-representability and evaluation bias).

We suggest that our observational findings be supported by prospective studies with multivariate analysis.

At this time, it seems that a particular phenotype for the INA induced by HEV has emerged: male with bilateral and asymmetric type, and a frequent phrenic nerve involvement. The excess risk of phrenic paresis supports the HEV earlier and curative treatment in addition to the symptomatic therapy.

  Additionally, a review of the literature was performed in June 2016, by searching the following combinations of keywords: A] (((neuralgic amyotrophy) OR brachial neuritis) OR parsonage turner) AND ((hepatitis E) OR HEV). Neuralgic amyotrophy triggered by HEV. DOI: 10.1007/s00415-017-8433-z 4 B] (((neuralgic amyotrophy) OR brachial neuritis) OR parsonage turner) AND (phrenic OR diaphragmatic).

Table 1 :

 1 Cases of neuralgic amyotrophy presumed triggered by hepatitis E virus ALP=alkaline phosphatase; AST=aspartate aminotransferase; CSF=cerebrospinal fluid; g-GT=gamma glutamyl transferase; HEV g.=heptatitis E virus genotype; L=left; LDH=lactate dehydrogenase; Nv=normal value; PP=phrenic paresis; R=right; Y/N=yes/no; a Unknown=No source for HEV transmission, including travel, cases in the family, eating undercooked pork or drinking unsafe water

	No. First	Age	Sex Associated	H	Transmission Hepatic biology	Pain location	Pain	Associated neurologic	P	Follow-up and clinical
	authors	(y)		medical	E			(at the onset)		duration	features	P	improvement
	year			condition	V			IU/L
					g.
	1 Fong 2009	53	M	None		Travel in		ALT 2547,	Shoulder,	6wk	Amyotrophy, proximal paresis,	N 2y / remission
	[11]					Equator		bilirubin 20	radiating to		sensory loss, abnormal
						(faecal-oral	µmol/L, ALP	the chest		reflexes in the upper limbs.
						transmission)	247
	2 Rianthavorn	49	M		3f Travel in India ALT 795, AST	Shoulders	5wk	Amyotrophy, paresis, sensory	N 10mo / partial improvement
	2010 [12]							155, bilirubin 32			loss, abnormal reflexes in the
								µmol/L, ALP			upper limbs
								284,
	3 Kamar 2011	38	M	Diabetes	3e	_		ALT 1160,	Upper arms	_	Amyotrophy, paresis, sensory	N 18mo / partial improvement
	[13]			type 1				bilirubin			loss, abnormal reflexes in the
								70µmol/L			upper limbs. Left vocal cord
											paralysis, normal functions of
											the other cranial nerves
	4 Carli 2012	30	M	None		_		ALT 1518, AST	L shoulder	_	Amyotrophy, paresis inthe	N _ , treatment by bolus of
	[14]							517, γ-GT 59			upper limbs. No sensory loss,	corticosteroids / partial
											normal reflexes. Serpiginous	improvement
											choroiditis in the L eye.
	5 Inghilleri	28	M			Pork sausage	ALT 1007, AST	Shoulders	15d	Amyotrophy, paresis,	N	_
	2012 [15]					and ingestion	147, ALP 162	and upper		abnormal reflexes in the upper
						of raw meat		arms		limbs.
	6 Cheung	56	M			Pork sausage AST 234,	Upper limb	10mo Paresis in a upper limb	N 10mo / no improvement
	2012 [16]							bilirubin 30
	7 Peri 2013	53	M		3 Unknown	a	µmol/L ALT 1768, AST	Shoulders	6 wk Proximal paresis, sensory	N 3mo / remission
	[17]							696, bilirubin			loss, abnormal reflexes in the
								129µmol/L, g-			upper limbs.
								GT 404
	8 Deroux	38	M			_		ALT 1288, AST	_	_	_	N 6mo / remission
	2013 [18]							988
	9 Valence	62	M			Travel in		Cytolysis and	Shoulders	_	Paresis in the upper limbs.	Y	_
	2014 [19]					Africa		cholestasis			Dysphonia
								2xNV
	10 Chirio 2014	32	M		3	_		ALT 10xNV	Shoulders	_	Amyotrophy, proximal paresis,	N 2mo / partial improvement
	[20]										sensory loss, in the upper
											limbs.
	11 Chirio 2014	34	M		3	_		ALT 10xNV	Shoulders	_	Amyotrophy, proximal paresis,	N 2mo / partial improvement
	[20]										sensory loss, in the upper
	12 Motte 2014	52	M	Prostatic	3f Unknown	a	ALT 590, AST	Cervicodynia	limbs. 2mo Amyotrophy, proximal paresis,	Y	_
	[21]			adenoma				179, LDH 262,	to the upper		in the upper limbs. Paresis
								ALP 48, g-GT	limbs		and abnormal reflexes in the
								212			lower limbs. No sensory loss
	13 Moisset	36	M	None	3f	_		ALT 1707, AST	Shoulders	_	Proximal paresis in the upper	N 6mo, treatment by IVIg 5d
	2014 [22]							554, bilirubin 54			limbs. Hypophonia, dysphagia	and ribavirin 35d / partial
								µmol/L, ALP				improvement
								311, g-GT 737
	14 Deroux	38	M	None		_		ALT 1612, AST	L shoulder	_	Amyotrophy, proximal paresis	N 4mo / partial improvement
	2014 [23]							772, (bilirubin,			in the upper limbs.
								ALP, g-GT :
								normal)
	15 Van Eijk	37	M			_		ALT 313,	_	_	Patchy bilateral brachial	N 6mo / partial improvement
	2014 [7]							bilirubin 4			plexus involvement
								µmol/L, ALP
								113
	16 Van Eijk	35	M			_		_	_	_	Patchy bilateral brachial	N 6mo / partial improvement
	2014 [7]										plexus involvement
	17 Van Eijk	40	F		3	_		ALT 183,	_	_	Patchy bilateral brachial	N 6mo / partial improvement
	2014 [7]							bilirubin 3			plexus involvement
								µmol/L, ALP
								177
	18 Van Eijk	36	M			_		_	_	_	Patchy bilateral brachial	N 6mo / partial improvement
	2014 [7]										plexus involvement
	19 Van Eijk	34	M			_		ALT 224	_	_	Patchy bilateral brachial	N 6mo / partial improvement
	2014 [7]										plexus involvement
	20 Woolson	late	M		3	_		ALT 1160,	_	_	Bilateral NA	N 12mo / partial improvement
	2014 [24]	30th						bilirubin 70
								µmol/L
	21 Woolson	late	M		3	_		ALT 319,	_	_	Bilateral NA	N 12mo / partial improvement
	2014 [24]	30th						bilirubin 4
								µmol/L
	22 Woolson	late	M			_		_	_	_	Bilateral NA	N 12mo / partial improvement
	2014 [24]	30th
	23 Theochari	65	M			_		ALT 1368, ALP	_	_	Bilateral NA	N 10mo, treatment by
	2015 [25]							167				corticosteroid / remission
	24 Décard	47	M	Hypertensi		_		ALT 106, AST	Shoulders	Severa	Amyotrophy, proximal paresis	N 12mo, treatment by
	2015 [26]			on, type 2				221, g-GT 239	and upper	l days	in the upper limbs. No sensory	corticosteroids / partial
				non-insulin					arms		loss, normal reflexes.	improvement
				requiring
				diabetes

The authors thank Laura Ravasi who provided medical writing services and Hélène


Nicolas Champagnat

Denis Villemonais

Uniform convergence to the Q-process

Keywords: quasi-stationary distribution, Q-process, uniform exponential mixing property, conditional ergodic theorem 2010 Mathematics Subject Classification. 60J25, 37A25, 60B10

Introduction

Let (Ω, (F t ) t≥0 , (X t ) t≥0 , (P x ) x∈E∪{∂} ) be a time homogeneous Markov process with state space E ∪ {∂}, where E is a measurable space. We assume that ∂ ∈ E is an absorbing state for the process, which means that X s = ∂ implies X t = ∂ for all t ≥ s, P x -almost surely for all x ∈ E. In particular, τ ∂ := inf{t ≥ 0, X t = ∂} is a stopping time. We also assume that P x (τ ∂ < ∞) = 1 and P x (t < τ ∂ ) > 0 for all t ≥ 0 and ∀x ∈ E.

A probability measure α on E is called a quasi-stationary distribution if

P α (X t ∈ • | t < τ ∂ ) = α, ∀t ≥ 0.
We refer the reader to [START_REF] Méléard | Quasi-stationary distributions and population processes[END_REF][START_REF] Van Doorn | Quasi-stationary distributions for discrete-state models[END_REF][START_REF] Collet | Quasi-Stationary Distributions: Markov Chains, Diffusions and Dynamical Systems[END_REF] and references therein for extensive developments and several references on the subject. It is well known that a probability measure α is a quasi-stationary distribution if and only if there exists a probability measure µ on E such that

lim t→+∞ P µ (X t ∈ A | t < τ ∂ ) = α(A) (1.1)
for all measurable subsets A of E.

In [START_REF] Champagnat | Exponential convergence to quasistationary distribution and Q-process[END_REF], we provided a necessary and sufficient condition on X for the existence of a probability measure α on E and constants C, γ > 0 such that

P µ (X t ∈ • | t < τ ∂ ) -α T V ≤ Ce -γt , ∀µ ∈ P(E), t ≥ 0, (1.2)
where • T V is the total variation norm and P(E) is the set of probability measures on E. This immediately implies that α is the unique quasistationary distribution of X and that (1.1) holds for any initial probability measure µ.

The necessary and sufficient condition for (1.2) is given by the existence of a probability measure ν on E and of constants t 0 , c 1 , c 2 > 0 such that

P x (X t 0 ∈ • | t 0 < τ ∂ ) ≥ c 1 ν, ∀x ∈ E and P ν (t < τ ∂ ) ≥ c 2 P x (t < τ ∂ ), ∀t ≥ 0, x ∈ E.
The first condition implies that, in cases of unbounded state space E (like N or R + ), the process (X t , t ≥ 0) comes down from infinity in the sense that, there exists a compact set

K ⊂ E such that inf x∈E P x (X t 0 ∈ K | t 0 τ ∂ ) > 0.
This property is standard for biological population processes such as Lotka-Volterra birth and death or diffusion processes [START_REF] Cattiaux | Quasi-stationary distributions and diffusion models in population dynamics[END_REF][START_REF] Champagnat | Lyapunov criteria for uniform convergence of conditional distributions of absorbed Markov processes[END_REF]. However, this is not the case for some classical models, such as linear birth and death processes or Ornstein-Uhlenbeck processes.

Many properties can be deduced from (1.2). For instance, this implies the existence of a constant λ 0 > 0 such that

P α (t < τ ∂ ) = e -λ 0 t
and of a function η : E → (0, ∞) such that α(η) = 1 and

lim t→+∞ sup x∈E e λ 0 t P x (t < τ ∂ ) -η(x) = 0 (1.3)
as proved in [START_REF] Champagnat | Exponential convergence to quasistationary distribution and Q-process[END_REF]Prop. 2.3]. It also implies the existence and the exponential ergodicity of the associated Q-process, defined as the process X conditioned to never be extinct [START_REF] Champagnat | Exponential convergence to quasistationary distribution and Q-process[END_REF]Thm. 3.1]. More precisely, if (1.2) holds, then the family (Q x ) x∈E of probability measures on Ω defined by

Q x (Γ) = lim t→+∞ P x (Γ | t < τ ∂ ), ∀Γ ∈ F s , ∀s ≥ 0, (1.4)
is well defined and the process (Ω, (F t ) t≥0 , (X t ) t≥0 , (Q x ) x∈E ) is an E-valued homogeneous Markov process. In addition, this process admits the unique invariant probability measure (sometimes refered to as the doubly limiting quasi-stationary distribution [START_REF] Flaspohler | Quasi-stationary distributions for absorbing continuous-time denumerable Markov chains[END_REF])

β(dx) = η(x)α(dx)
and there exist constants C , γ > 0 such that, for any x ∈ E and all t ≥ 0,

Q x (X t ∈ •) -β T V ≤ C e -γ t . (1.5)
The measure β

The first aim of the present note is to refine some results of [START_REF] Champagnat | Exponential convergence to quasistationary distribution and Q-process[END_REF] in order to get sharper bounds on the convergence in (1.3) and to prove that the convergence (1.4) holds in total variation norm, with uniform bounds over the initial distribution (see Theorem 2.1). Using these new results, we obtain in Corollary 2.3 that the uniform exponential convergence (1.2) implies that, for all bounded measurable function f : E → R and all T > 0,

E x 1 T T 0 f (X t ) dt | T < τ ∂ - E f dβ ≤ a f ∞ T , (1.6)
for some positive constant a. This result improves the very recent result obtained independently by He, Zhang and Zu [6, Thm. 2.1] by providing the convergence estimate in 1/T . The interested reader might look into [START_REF] He | On the quasi-ergodic distribution of absorbing Markov processes[END_REF] for nice domination properties between the quasi-stationary distribution α and the probability β.

The second aim of this note is to prove that the existence of the Q-process with uniform bounds in (1.4) and its uniform exponential ergodicity (1.5) form in fact a necessary and sufficient condition for the uniform exponential convergence (1.2) toward a unique quasi-stationary distribution.

Main results

In this first result, we improve (1.3) and provide a uniform exponential bound for the convergence (1.4) of the conditioned process toward the Qprocess.

Theorem 2.1. Assume that (1.2) holds. Then there exists a positive constant a 1 such that

e λ 0 t P x (t < τ ∂ ) -η(x) ≤ a 1 e λ 0 t P x (t < τ ∂ )e -γt , (2.1)
where λ 0 and η are the constant and function appearing in (1.3) and where γ > 0 is the constant from (1.2).

Moreover, there exists a positive constant a 2 such that, for all t ≥ 0, for all Γ ∈ F t and all T ≥ t,

Q x (Γ) -P x (Γ | T < τ ∂ ) T V ≤ a 2 e -γ(T -t) , (2.2)
where (Q x ) x∈E is the Q-process defined in (1.4).

We emphasize that (2.1) is an improvement of (1.3), since the convergence is actually exponential and, in many interesting examples, inf x∈E P x (t < τ ∂ ) = 0. This is for example the case for elliptic diffusion processes absorbed at the boundaries of an interval, since the probability of absorption converges to 1 when the initial condition converges to the boundaries of the interval. The last theorem has a first corollary.

Corollary 2.2. Assume that (1.2) holds. Then there exists a positive constant a 3 such that, for all T > 0, all probability measure µ T on [0, T ] and all bounded measurable functions f :

E → R, E x T 0 f (X t )µ T (dt) | T < τ ∂ - E f dβ ≤ a 3 f ∞ T 0 e -γ t + e -γ(T -t) µ T (dt). (2.3)
This follows from (2.2), the exponential ergodicity of the Q-process stated in (1.5) and the inequality

E x T 0 f (X t )µ T (dt) | T < τ ∂ - E f dβ ≤ T 0 E x (f (X t ) | T < τ ∂ ) -E Qx (f (X t )) µ T (dt) + T 0 E Qx (f (X t )) - E f dβ µ T (dt),
where E Qx is the expectation with respect to Q x .

In particular, choosing µ T as the uniform distribution on [0, T ], we obtain a conditional ergodic theorem.

Corollary 2.3. Assume that (1.2) holds. Then there exists a positive constant a 4 such that, for all T > 0 and all bounded measurable functions f :

E → R, E x 1 T T 0 f (X t ) dt | T < τ ∂ - E f dβ ≤ a 4 f ∞ T .
Considering the problem of estimating β from N realizations of the unconditioned process X, one wishes to take T as small as possible in order to obtain the most samples such that T < τ ∂ (of order N T = N e -λ 0 T ). It is therefore important to minimize the error in (2.3) for a given T . It is easy to check that µ T = δ t 0 with t 0 = γT /(γ + γ ) is optimal with an error of the order of exp(-γ γT /(γ + γ )). Combining this with the Monte Carlo error of order 1/ √ N T , we obtain a global error of order

e λ 0 T /2 √ N + e -γγ T /(γ+γ ) .
In particular, for a fixed N , the optimal choice for T is T ≈

log N λ 0 +2γγ /(γ+γ )
and the error is of the order of N -ζ with ζ = γγ 2γγ +λ 0 (γ+γ ) . Conversely, for a fixed T , the best choice for N is N ≈ exp((λ 0 + 2γγ /(γ + γ ))T ) and the error is of the order of exp(-γγ T /(γ + γ )).

We conclude this section with a converse to Theorem 2.1. More precisely, we give a converse to the fact that (1.2) implies both (1.5) and (2.2).

Theorem 2.4. Assume that there exists a Markov process (Q x ) x∈E with state space E such that, for all t > 0, lim

T →+∞ sup x∈E Q x (X t ∈ •) -P x (X t ∈ • | T < τ ∂ ) T V = 0 (2.4)
and such that

lim t→+∞ sup x,y∈E Q x (X t ∈ •) -Q y (X t ∈ •) T V = 0. (2.5)
Then the process (P x ) x∈E admits a unique quasi-stationary distribution α and there exist positive constants γ, C such that (1.2) holds.

It is well known that the strong ergodicity (2.5) of a Markov process implies its exponential ergodicity [8, Thm. 16.0.2]. Similarly, we observe in our situation that, if (2.4) and (2.5) hold, then the combination of the above results implies that both convergences hold exponentially. For all x ∈ E, we set

η t (x) = P x (t < τ ∂ ) P α (t < τ ∂ ) = e λ 0 t P x (t < τ ∂ ),
and we recall from [2, Prop. 2.3] that η t (x) is uniformly bounded w.r.t. t ≥ 0 and x ∈ E. By Markov's property

η t+s (x) = e λ 0 (t+s) E x (1 t<τ ∂ P Xt (s < τ ∂ )) = η t (x)E x (η s (X t ) | t < τ ∂ ) .
By (1.2), there exists a constant C independent of s such that

E x (η s (X t ) | t < τ ∂ ) - E η s dα ≤ C e -γt .
Since η s dα = 1, there exists a constant a 1 > 0 such that, for all x ∈ E and s, t ≥ 0,

η t+s (x) η t (x) -1 ≤ a 1 e -γt .
Hence, multiplying on both side by η t (x) and letting s tend to infinity, we deduce from (1.3) that, for all x ∈ E, |η(x) -η t (x)| ≤ a 1 e -γt η t (x), ∀t ≥ 0, which is exactly (2.1). We also deduce that

1 -a 1 e -γt η t (x) ≤ η(x) ≤ 1 + a 1 e -γt η t (x) (3.1)
and hence, for t large enough,

η(x) 1 + a 1 e -γt ≤ η t (x) ≤ η(x) 1 -a 1 e -γt . (3.2)
Let us now prove the second part of Theorem 2.1. For any t ≥ 0, Γ ∈ F t and 0 ≤ t ≤ T ,

P x (Γ | T < τ ∂ ) = P x (Γ ∩ {T < τ ∂ }) P x (T < τ ∂ ) = e λ 0 T P x (Γ ∩ {T < τ ∂ }) η(x) η(x) e λ 0 T P x (T < τ ∂ )
.

We deduce from (2.1) that

η(x) e λ 0 T P x (T < τ ∂ ) -1 ≤ a 1 e -γT ,
while, for all T > log a 1 γ , (3.2) entails

e λ 0 T P x (Γ ∩ {T < τ ∂ }) η(x) ≤ η T (x) η(x) ≤ 1 1 -a 1 e -γT .
Hence, for all t ≥ 0 and all T > log a 1 γ ,

P x (Γ | T < τ ∂ ) - e λ 0 T P x (Γ ∩ {T < τ ∂ }) η(x) ≤ a 1 e -γT 1 -a 1 e -γT . (3.3)
Now, the Markov property implies that

P x (Γ ∩ {T < τ ∂ }) = E x (1 Γ P Xt (T -t < τ ∂ )) ,
and we deduce from (3.3) that, for all T > t + log a 1 γ ,

e λ 0 (T -t) P Xt (T -t < τ ∂ ) -η(X t ) ≤ a 1 e -γ(T -t) 1 -a 1 e -γ(T -t) η(X t ).
Thus we have

e λ 0 T P x (Γ ∩ {T < τ ∂ }) η(x) - e λ 0 t E x (1 Γ η(X t )) η(x) ≤ e λ 0 t η(x) E x 1 Γ e λ 0 (T -t) P Xt (T -t < τ ∂ ) -η(X t ) ≤ a 1 e -γ(T -t) 1 -a 1 e -γ(T -t) e λ 0 t E x (η(X t )) η(x) = a 1 e -γ(T -t) 1 -a 1 e -γ(T -t) ,
where we used the fact that E x η(X h ) = e -λ 0 h η(x) for all h > 0 (see [START_REF] Champagnat | Exponential convergence to quasistationary distribution and Q-process[END_REF]Prop. 2.3]). This and (3.3) allows us to conclude that, for all t ≥ 0 and all T > t + log a 1 γ ,

P x (Γ | T < τ ∂ ) - e λ 0 t E x (1 Γ η(X t )) η(x) ≤ 2a 1 e -γT 1 -a 1 e -γT . Since Q x (Γ) = e λ 0 t E x (1 Γ η(X t )) /η(x) (see [2, Thm. 3.1 (ii)]), we deduce that (2.2) holds true.
This concludes the proof of Theorem 2.1.

Proof of Theorem 2.4

We deduce from (2.4) and (2.5) that there exists t 1 > 0 and T 1 > 0 such that, for all T ≥ T 1 , sup

x,y∈E

P x (X t 1 ∈ • | T < τ ∂ ) -P y (X t 1 ∈ • | T < τ ∂ ) T V ≤ 1/2.
In particular, for all s ≥ 0 and all T ≥ s + T 1 , sup

x,y∈E

δ x R T s,s+t 1 -δ y R T s,s+t 1 T V ≤ 1/2, (3.4)
where, for all 0 ≤ s ≤ t ≤ T , R T s,t is the linear operator defined by

δ x R T s,t f = E x (f (X t-s ) | T -s < τ ∂ ) = E(f (X t ) | X s = x, T < τ ∂ ) = δ x R T -s 0,t-s f,
where we used the Markov property. Now, for any T > 0, the family (R T s,t ) 0≤s≤t≤T is a Markov semi-group. This semi-group property and the contraction (3.4) classically imply that, for all T ≥ T 1 , sup x,y∈E δ x R T 0,T -δ y R T 0,T T V ≤ (1/2) T -T 1 /t 1 .

Then, proceeding as in [2, Section 5.1], we deduce that (1.2) holds true. This concludes the proof of Theorem 2.4.

Université de Lorraine, IECL, UMR 7502, Campus Scientifique, B.P. 70239, Vandoeuvre-lès-Nancy Cedex, F-54506, France

CNRS, IECL, UMR 7502, Vandoeuvre-lès-Nancy, F-54506, France

Inria, TOSCA team, Villers-lès-Nancy, F-54600, France. E-mail: Nicolas.Champagnat@inria.fr, Denis.Villemonais@univ-lorraine.fr


Philippe Tigréat
email: philippe.tigreat@telecom-bretagne.eu

Pierre-Henri Horrein
email: ph.horrein@telecom-bretagne.eu

Vincent Gripon
email: vincent.gripon@telecom-bretagne.eu

Massively Parallel Implementation of Sparse Message Retrieval Algorithms in Clustered Clique Networks

Keywords: Clustered Clique Networks, Parallel processing, CUDA, GPGPU

Auto-associative memories are a family of algorithms designed for pattern completion. Many of them are based on neural networks, as is the case for Clustered Clique Networks which display competitive pattern retrieval abilities. A sparse variant of these networks was formerly introduced which enables further improved performances. Specific pattern retrieval algorithms have been proposed for this model, such as the Global-Winners-Take-All and the Global-Losers-Kicked-Out. We hereby present accelerated implementations of these strategies on graphical processing units (GPU). These schemes reach interesting factors of acceleration while preserving the retrieval performance.

I. INTRODUCTION

Associative memories are a type of computer memories that are accessed by content instead of address. Hetero-associative memory link different patterns together by pairs, and when requested with one pattern they return the associated one. Auto-associative memories link together the subparts of a pattern, and perform pattern completion when submitted a subset of a stored memory element.

Among the latter category are models such as Hopfield [START_REF] Hopfield | Neural networks and physical systems with emergent collective computational abilities[END_REF] and Willshaw [START_REF] Willshaw | Nonholographic associative memory[END_REF]. Both models are based on neural networks. Hopfield neural networks come with a fully-interconnected set of neurons, and the connection weights are adapted to the set of stored messsages so that each message becomes a local minimum of the network's energy function. Willshaw networks differ from this scheme in that their connectivity is very sparse, and only when two neurons are part of a same message do they get connected. This characteristic grants Willshaw networks with a greater storage capacity than Hopfield, meaning they allow a better message retrieval success rate on average for given numbers of neurons and stored messages.

Clustered Clique Networks are another model of associative memories originally developed by Gripon and Berrou [START_REF] Gripon | Sparse neural networks with large learning diversity[END_REF]- [START_REF]Nearly-optimal associative memories based on distributed constant weight codes[END_REF]. They leverage further benefits from sparseness than Willshaw, as a result of a clustering of the neurons set. Neurons are parted in equal sized clusters and intra-cluster connections are prohibited. In the first published version of these networks, a stored message is supported by a set of neurons with one picked in each cluster of the network. This fullyinterconnected pattern is called a clique. Aliabadi et al [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] has investigated the storage of sparse messages in these networks, meaning messages composed by neurons located in a subset of the network's clusters. This specific case can benefit from an adaptation of the retrieval procedure used in classic CCNs.

Aboudib et al [START_REF] Aboudib | A study of retrieval algorithms of sparse messages in networks of neural cliques[END_REF] introduces two such variants of the Winner-Take-All algorithm. These are the Global Winners-Take-All (GWsTA) and the Global Losers-Kicked-Out (GLsKO), and come with a significant improvement in retrieval success rate.

The present work focuses on the acceleration of these algorithm variants by means of Graphical Processing Units (GPUs). These chipsets allow the massively parallel computation of operations on data, and are currently widely used to accelerate neural network implementations. To achieve this, one must ensure a sufficient amount of operations to execute can actually be performed in parallel without altering the accuracy of the whole computation. The main requirement for this is the independence of the parallelly processed data. We aim at making clear here how we adapt the former serial implementation under this constraint. We use the Compute Unified Device Architecture (CUDA) which is a coding framework specifically designed to program GPUs, supported by C/C++. Yao [START_REF] Yao | A GPU-based associative memory using sparse Neural Networks[END_REF], [START_REF] Yao | A massively parallel associative memory based on sparse neural networks[END_REF] formerly achieved acceleration of the classic CCN algorithm using CUDA.

The paper is organized as follows. Section II summarizes the Clustered Clique Network model and algorithm, and its adaptation to sparse messages. Section III describes the strategies we use to adapt the specific retrieval algorithms for sparse messages on a massively parallel architecture. Section IV presents the results we achieve in terms of acceleration and message retrieval performance depending on different parameters of the network.

II. CCNS

A. Structure

Clustered Clique Networks (CCN) are composed of a set of N neurons, divided in clusters all containing the same number of neurons. A connection can be drawn between two of these neurons only if they belong to different clusters. A memory element is then supported by a choice of neurons all located in different clusters. The storage of such a message in the network consists in connecting all of its neurons together. The pattern of connections so created is called a clique.

By the separation of the network into clusters, it makes use of sparseness to a higher extent than former states-ofthe-art like Willshaw networks, thus lessening the amount of overlapping between stored messages and easing the search of the unknown elements of patterns to complete. Hence, CCN display a high storage capacity and the number of messages they are able to store and retrieve successfully grows quadratically as a function of the number of neurons in the network. During retrieval, the commonly used procedure consists in propagating activation from the neurons of a request through all their existing connections to other neurons, followed by a selection of neurons based on the resulting activity scores. A local Winner-Takes-All rule is generally used to perform the selection. It consists in keeping active only the neurons possessing the highest activity score in each cluster. This is an iterative process, and it is repeated as long as the stopping criterion is not reached. Activation is computed again with the selected neurons, and the rule is applied once again.

B. Sparse Messages

Sparse Clustered Clique Networks as introduced by Aliabadi et al [START_REF] Aliabadi | Storing sparse messages in networks of neural cliques[END_REF] have the same principle as classic CCNs. The main difference is that stored messages no longer need to possess a neuron in every cluster of the network, but only in a selected subset of the clusters. This characteristic raises the diversity of messages, that is the total number of different messages one can store in the network. However, it brings the added difficulty of not knowing which clusters support the searched message. The high level of sparsity usually makes up for this, and the retrieval error rates are much lower in general as compared to full messages CCNs. Figure 1 shows an example of such a network.

C. Specific Retrieval Algorithms for Sparse Messages

The retrieval of sparse messages can be performed more efficiently by using strategies specifically tailored for this case. These will typically make use of the set of activity scores observed after stimulation instead of solely focusing on local winners. Two such strategies can be found in [START_REF] Aboudib | A study of retrieval algorithms of sparse messages in networks of neural cliques[END_REF] namely the Global Winners-Take-All and the Global Losers-Kicked-Out.

In Global Winners-Take-All, one wants to retain a number of neurons equal or superior to a predefined minimum. The minimal threshold score for a fanal to be retained is adapted in order to fulfill this requirement. After the first activity computation step, we iteratively select a number of fanals which are then used to compute activity for next iteration. For instance, if after a stimulation of the network, there are 7 activated fanals with scores {2,3,1,4,3,4,4} and we know at least 4 neurons must be retained, a lower threshold of 3 will be applied on the scores of neurons to remain active. This way the 5 neurons with an activity of 3 or more are kept as output for the ongoing iteration, while the 2 neurons with scores 1 and 2 are discarded.

In Global Losers-Kicked-Out, focus is put on active neurons with relatively low scores. After the initial activity computation, the aim is to discard a number of fanals iteratively, instead of selecting them. Among the non-zero scores observed after stimulation, a fixed number of lower score values is set so that neurons with these activity levels are put to zero. In the previous example, if we know 3 lower score values must be discarded, this will designate the 4 neurons in the list with scores {2,1,3,3} and only the 3 neurons with an activity of 4 will remain active. Nonetheless, it has been shown that banning only the one minimal score in the list gives the best results.

III. ACCELERATIONS

A. GPU usage

Graphics Processing Units (GPU) are highly parallel processors. They are based on a high number of processing units, which can be used to process data using a Single Instruction Multiple Data (SIMD) paradigm. In this paradigm, the same instruction can be applied simultaneously on all elements of a vector. As a result, GPU can be very efficient when processing algorithm without data dependency. For example, loops in which each iteration can be performed regardless of previous iterations can be very efficiently accelerated. In order to fully benefit from the GPU, the selection algorithms have been adapted to follow this paradigm.

B. GWsTA

Global Winners-Take-All aims at retaining a number of active fanals superior or equal to a predefined minimum. In order to do this, we compute for each possible score value the number of fanals that reach it. In practice, scores counts are first computed in parallel for the different clusters, with one thread computing the counts for the different scores inside of one cluster, as shown on Figure 2. In a second step, the scores counts are accumulated by reduction to get their totals over the whole network. A number of threads is launched, equal to half the number of clusters. Each thread sums the vectors of scores counts of two clusters located in different halves of the network. The operation is repeated on the new set of counts lists and so forth, each time halving the number of computed values until one single list remains with the global sums. A final thread then runs a loop over the possible score values in descending order, accumulating scores counts until the sum equals or exceeds the minimum expected number of active neurons. As soon as this criterion is met, looping stops and the last considered score value is kept as the minimum score for nodes to remain active. A thresholding with binarization is then applied on nodes activities, with scores above the threshold being put to one while other neurons are all at zero. This operation is performed in parallel, where each thread computes the output values for one neuron as illustrated by Figure 3. The resulting set of active fanals is then returned as output for the current iteration.

C. GLsKO

For our implementation of the Global Losers-Kicked-Out we set to 1 the number of low scores to be discarded. This is known to give the best results, and it reduces the complexity of the algorithm, since we only need to find the minimum score and discard neurons with this score. Hence we focus on finding the global minimum among all non-zero scores present after stimulation. To this end the minimum positive scores inside of the different clusters are first computed in parallel. A reduction strategy is employed on the resulting set of local minimums to compute the global one. At every step of the reduction, each thread compares the minimums obtained for two clusters, keeping the lowest of the two values. The set of values is thus halved iteratively to result in a single minimum. Along with the strictly positive minimum, a global boolean value is computed indicating the presence of several different non-zero scores, as described in Figure 5. Here the fact that several threads modify the output boolean does not alter its accuracy, as it indicates whether two different non-zero scores are seen at least once. If this boolean is true at the end of the reduction procedure, the ThresholdBinarize operation is applied to compute the iteration output, setting off all nodes with the non-zero minimum score. Otherwise, only one strictly positive score is present and all its representatives are kept active. This case where all stimulated fanals have equal scores is also a good stopping criterion for the GLsKO.

IV. EXPERIMENTS

Our GPU experiments were performed with an NVIDIA GeForce 780 Ti, with a frequency of 900MHz and 3GB of memory. The CPU used is a 3.3GHz Intel Core i7. Each message is generated randomly by first picking a fixed size subset of the available clusters, and then choosing one fanal in each selected cluster. In the retrieval phase, requests are formed from stored messages by dropping out a fix number of their nodes.

Figure 7 shows the execution times in seconds for the retrieval of 500 stored messages as a function of the number of erasures per queries, in networks with 16 clusters and 32 fanals per cluster. Each message is made of 12 neurons and the corresponding query is a subsample from the message with a size going from 11 down to 1.

The execution times are very close at the beginning and differ when the number of erasures exceeds about half the message size, due to the different stopping criteria employed. For higher number of erasures the GLsKO takes more and more time, with a maximum when the request is made of a single neuron. The stopping criterion used is the equality of scores of all neurons active after stimulation. The GWsTA scheme takes less time for higher numbers of erasures. Its stopping criteria are convergence, i.e. equality between input and output at a given iteration, and a maximum number of iterations, hereby 8.

Figure 8 shows the execution time of the two schemes for networks of 64 clusters of 724 fanals each and 500 stored messages. Here messages are made of 56 nodes and the number of erasures applied to create the corresponding queries varies from 2 to 54. The two schemes make no retrieval error except for the GLsKO which fails to retrieve one pattern out of 500 when provided with 2-nodes requests. The number of iterations is impacted and alters the speed of retrieval.

With 2 erasures the GLsKO scheme takes 0.093573 seconds to complete retrieval of the whole set of messages, when a serial implementation of the same algorithm takes 4.74341 seconds running on CPU. This brings an acceleration factor of 50.69. Parallel and serial versions of the GWsTA use respectively 0.0246 seconds and 4.16881 seconds in the case with 2nodes requests, giving an acceleration of 169.319. To be fair, it has to be stated that our serial and parallel implementations are not tailored to leverage a high sparseness of the connectivity graph, which could bring further improvements. To this end, one would want to represent the set of connections as lists of indexes instead of a boolean matrix.

V. SUMMARY We present parallel implementations of two specific strategies for the retrieval of sparse messages in Clustered Clique Networks. The accelerations we obtain are interesting and bring prospects for the use of these schemes in real-world applications, such as machine learning. Future work may indeed focus on parallel implementations for the retrieval of sequences in Clustered Clique Networks, with practical trials on high-scale datasets. Using Clustered Clique Networks in conjunction with Deep Learning for image classification is also one of our main directions of research.

Figure 1 .

 1 Figure 1. Example of a sparse Clustered Clique Network with 12 clusters and messages of order 4.

Figure 2 .

 2 Figure 2. Example of the parallel computation of local scores counts in the different clusters. Each thread computes a vector of the numbers of fanals possessing the different possible score values in one cluster.

Figure 3 .

 3 Figure 3. Example of the threshold-binarize operation as performed on GPU. Non-zero scores below the threshold are put to 0 while those above it are set to 1. Input and output values are in one-to-one relationships, hence computation can be performed in parallel.

Figure 4 .

 4 Figure 4. Illustration of the GWsTA procedure after the first iteration of retrieval. All operations in rectangular boxes are performed on GPU.

Figure 5 .

 5 Figure 5. Illustration of the operation SeveralScoresPresent performed in the parallel implementation of GLsKO. This operation is performed in the same threads as the search of the non-zero minimum by reduction among local minimums. Each thread compares the minimum scores of two clusters to assess whether they are different and strictly positive. Several threads may access and modify the output boolean without affecting the result's accuracy.

Figure 6 .

 6 Figure 6. Illustration of the GLsKO procedure after the first iteration of retrieval. The global minimum non-zero score is computed by reduction, along with a boolean indicating if several non-zero scores are present. The output is computed using the parallel ThresholdBinarize operation, excluding the lowest non-zero score in the latter case.

Figure 7 .

 7 Figure 7. Execution time of the parallel versions of GWsTA and GLsKO for the retrieval of 500 sparse messages as a function of the number of erasures per query, in a network of 16 clusters and 32 fanals per cluster. Messages are made of 12 nodes and the size of sub-sampled queries goes from 11 nodes to 1.

Figure 8 .

 8 Figure 8. Parallel execution times of the two schemes for the retrieval of 500 messages of 56 nodes in a network with 64 clusters and 724 fanals per cluster. Queries size goes from 54 down to 2.

ACKNOWLEDGEMENTS

This work was supported by the European Research Council under Grant ERC-AdG2011 290901 NEUCOD.

The authors would like to thank NVIDIA for providing us with a free graphics card allowing to speed up computations for the experiments performed during this work.