Learning to be attractive: probabilistic computation with dynamic attractor networks

Alexander Gepperth 1 and Mathieu Lefort 2

Abstract-In the context of sensory or higher-level cognitive processing, we present a recurrent neural network model, similar to the popular dynamic neural field (DNF) model, for performing approximate probabilistic computations. The model is biologically plausible, avoids impractical schemes such as log-encoding and noise assumptions, and is well-suited for working in stacked hierarchies. By Lyapunov analysis, we make it very plausible that the model computes the maximum a posteriori (MAP) estimate given a certain input that may be corrupted by noise. Key points of the model are its capability to learn the required posterior distributions and represent them in its lateral weights, the interpretation of stable neural activities as MAP estimates, and of latency as the probability associated with those estimates. We demonstrate for in simple experiments that learning of posterior distributions is feasible and results in correct MAP estimates. Furthermore, a preactivation of field sites can modify attractor states when the data model is ambiguous, effectively providing an approximate implementation of Bayesian inference.

I. INTRODUCTION

A. Context and scope

This contribution is conducted in the context of cognitive neural modeling, essentially trying to elucidate how neural populations can perform probabilistic computations which are speculated to play an important role in the considerable performance of biological brains. However, especially the currently influential Bayesian paradigm, stating that the basic quantities to be manipulated and passed around should be probability distributions, poses considerable difficulties. Considering the limited operations available to neurons (essentially just weighted summation), the following questions may be asked: How should such distributions be represented on a neural level? How can the multiplications required for manipulating them be performed, not to speak of the necessary renormalizations after each such step? And, lastly, how can neural populations take optimal decisions based on distributions? Please see Sec. V for a discussion of these questions in hindsight.

This article presents a network model that proposes answers to some of these questions, is biologically plausible (see Sec. V for a discussion of this claim) and yet computationally (rather) inexpensive. The basic properties of the model, some of which are already demonstrated in [START_REF] Gepperth | Processing and transmission of confidence in recurrent neural hierarchies[END_REF], some of which shall be demonstrated or at least made plausible here, are as follows:

1 UIIS lab, FLOWERS team, INRIA, ENSTA ParisTech, 828 Blvd des Marechaux, 91762 Palaiseau (France), alexander@gepperth.net 2

LIRIS

laboratory, university Lyon 1, Lyon (France), mathieu.lefort@liris.cnrs.fr

• Bayesian processing is approximate: we approximate full posterior probabilities in each layer by their argmax (the MAP estimate) and its associate posterior probability, see Fig. 1 for a visualization of this process. Fig. 1. Illustration of the basic approximation process proposed in this article: the approximation of a full probability distribution by its argmax (the MAP estimate), and its associated probability. We assume a 2D neural layer where each neuron represents a certain stimulus, and a probability distribution defined over stimuli (and thus over neuronal sites). The depicted Gaussian distribution can then be approximated by the three numbers: 0.2 and 0.5 for the "position" of the MAP state, as well as 0.75 for its probability under the given distribution. In our neural implementation, these quantities are expressed by the 2D position of the attractor state as well as its latency (time to convergence).

MAP estimation learned data model

Fig. 2. Functionality of the proposed neural model. Given an afferent synaptic input (left), supposed to be population-coded so that activity at a certain place codes for a certain quantity/concept, a neural layer converges to an attractor state (right) with activity corresponding to the MAP estimate of the underlying true stimulus. Input is supposed to be a noisy version of this stimulus, which is recovered from input by making use of learned statistical properties of the data, indicated on the left by synaptic connections whose strength (indicated by brightness) represents the probability of co-activation. The time until the attractor state is reached is a measure for the posterior probability of the MAP estimate. In this concrete example, the left stimulus is suppressed because it has zero probability of co-occurrence with the right stimulus, and a weaker probability of co-occurrence with the central stimulus.

stimulus. We have already shown in [START_REF] Gepperth | Processing and transmission of confidence in recurrent neural hierarchies[END_REF] that the time to convergence to this estimate, denoted latency, is a measure of its actual probability under the posterior distribution.

It is these two quantities, the MAP estimate and its probability, that are therefore passed to subsequent neural layers, thus replacing the actual posterior distribution by an unimodal approximation (see Fig. 1 for a visualization). This approximation allows a mapping of probabilistic processing to a neural model, as the competitive dynamics prevents the representation of arbitrary distributions in the form of neural activities. Lastly, this model can be explicitly fed with prior distributions by pre-activating certain field sites, thus approximating Bayesian inference in the course of the MAP calculation process.

B. Rigorous problem formulation

The basic representational unit we consider here is a neural layer X where neural activities (firing rates) z X (t) are arranged on a rectangular grid.

For a layer receiving an input I, we suppose the input to be a noisy version of the underlying "true" stimulus T . Input and underlying stimulus are related by the probabilistic input model p( I| T ). The basic task we wish to accomplish is the estimation of the most likely "true" stimulus T * given the current input I, as well as its associated probability p( T = T * | I). In other words, what we are looking for is the socalled maximum a posteriori (MAP) estimate given the input stimulus I. This task must involve some way of estimating the "inverse" distribution p( T = T | I) which is related to the input model by the law of Bayes:

p( T = T | I) = p( I| T = T )p( T = T ) p( I) (1)
.

Summing up, we wish to approximate the following quantities:

T * = argmax T p( T = T | I) (2)
p * = p( T = T * | I)
Furthermore, we wish to influence the outcome of this computation by specifying priors for the "true" stimulus T . This is totally consistent with Bayes' formula, which in our case states that the basic distribution p( T = T | I) can be re-stated as

p( T = T | I) ∼ p( I| T = T )p( T = T ), (3)
where the a priori probability for the true stimulus appears explicitly on the right-hand side of the equation.

C. Modeling approach: competitive neural dynamics and space-latency coding

We will demonstrate by Lyapunov analysis that a competitive neural dynamic approach very similar to the dynamic neural field (DNF) model [START_REF] Erlhagen | Dynamic field theory of movement preparation[END_REF] converges to an attractor state that is maximally compatible with afferent and lateral input, the latter defined by lateral connections. Learning and representation of the distribution p( T = T | I) required for the computation of the MAP estimate in eqn.( 2) is realized by the lateral connections weights by linear regression. If the lateral input thus obtained can be associated with a likelihood distribution, this would make it very plausible that the attractor state of the model indeed corresponds to, or at least approximates, the MAP estimate 1 .

In previous works on the subject [START_REF] Gepperth | Processing and transmission of confidence in recurrent neural hierarchies[END_REF], we already showed that this model can represent both the MAP estimate and its probability under the posterior distribution, T * and p * , by its attractor state and the associated latency (time-toconvergence). This naturally gives rise to a novel way of encoding information, which we term space-latency coding [START_REF] Gepperth | Processing and transmission of confidence in recurrent neural hierarchies[END_REF], into neural populations, which is exclusively due to the dynamic properties of the DNF model. We already showed in [START_REF] Geppertha | Latency-based probabilistic information processing in a learning feedback hierarchy[END_REF], [START_REF] Gepperth | Latency-based probabilistic information processing in recurrent neural hierarchies[END_REF] that this space-latency code is well suited for transmitting information in a neural hierarchy, and demonstrated its appropriateness for realistic information processing.

D. Related work

There exists a large body of literature [START_REF] Wj Ma | Bayesian inference with probabilistic population codes[END_REF], [START_REF] Zemel | Probabilistic interpretation of population codes[END_REF], [START_REF] Gold | Neural computations that underlie decisions about sensory stimuli[END_REF], [START_REF] David | The bayesian brain: the role of uncertainty in neural coding and computation[END_REF], [START_REF] Cuijpers | Implementing bayes' rule with neural fields[END_REF], [START_REF] Rajesh | Bayesian computation in recurrent neural circuits[END_REF] on probabilistic aspects of neural coding. Most authors explicitly assume that neural population activity is directly related to probability distributions [START_REF] Wj Ma | Bayesian inference with probabilistic population codes[END_REF], [START_REF] Zemel | Probabilistic interpretation of population codes[END_REF], [START_REF] David | The bayesian brain: the role of uncertainty in neural coding and computation[END_REF], [START_REF] Cuijpers | Implementing bayes' rule with neural fields[END_REF], [START_REF] Rajesh | Bayesian computation in recurrent neural circuits[END_REF]; posing the question of how to appropriately multiply and renormalize neural activities in a biologically plausible way. A very influential idea for addressing this issue posits that single-neuron activity is in fact related to log-probability [START_REF] Cuijpers | Implementing bayes' rule with neural fields[END_REF], [START_REF] Rajesh | Bayesian computation in recurrent neural circuits[END_REF], [START_REF] Gold | Neural computations that underlie decisions about sensory stimuli[END_REF], which would allow to perform multiplication by summation, which neurons do easily. Other authors have questioned the practicability of this scheme [START_REF] Wj Ma | Bayesian inference with probabilistic population codes[END_REF] as it would require all sorts of cut-offs due to the unbounded nature of logarithms as probabilities go to zero, and which would have to be performed at each hierarchy level which is deemed unfeasible and lossy. An alternative approach [START_REF] Wj Ma | Bayesian inference with probabilistic population codes[END_REF] is to consider single neuron's firing rates as the realizations of Poisson-like random variables whose mean is determined by the match of neural preferences to afferent input. Under certain conditions, sums of two such variables can be proven to come from a distribution whose mean corresponds to the product of individual means, thus realizing a multiplication by summation.

E. Novelty

When comparing our approach to related work, several differences are evident: First of all, and different from [START_REF] Gold | Neural computations that underlie decisions about sensory stimuli[END_REF], [START_REF] Rajesh | Bayesian computation in recurrent neural circuits[END_REF], [START_REF] Cuijpers | Implementing bayes' rule with neural fields[END_REF], our approach does not treat neural activities as log-probabilities. More generally, and in line with [START_REF] Wj Ma | Bayesian inference with probabilistic population codes[END_REF], we do not treat the set of input activations S as a probability distribution but as the expression of match to individual neuron's preferences, corrupted by noise. In contrast to [START_REF] Wj Ma | Bayesian inference with probabilistic population codes[END_REF], however, we do not require this noise to have any particular form in order for our approach to work. Conversely, our approach does not profit from noise but tries to remove it to recover the underlying stimulus. This will work best if the noise is not on the single-neuron level but consists, e.g., in the apparition or fluctuation of local activity bubbles. A further, related difference to [START_REF] Wj Ma | Bayesian inference with probabilistic population codes[END_REF], who treat each neuron in absolute isolation from is neighbours, is that we consider data models that are global in the sense that they consider the values of other, not necessarily adjacent neurons in the approximation of the MAP estimate. In contrast to all of the approaches listed here, our model alone exploits the fact of using a dynamic model for encoding information, in this case by latency. Furthermore, as we could show in previous work [START_REF] Gepperth | Processing and transmission of confidence in recurrent neural hierarchies[END_REF], [START_REF] Geppertha | Latency-based probabilistic information processing in a learning feedback hierarchy[END_REF], our model accounts for the encoding, the transmission and the decoding of this so-called spacelatency code in a deep hierarchy as a direct consequence of its dynamics. Lastly, the idea of using a DNF-like model with a plastic, learned kernel is, to our knowledge, unique w.r.t. to related work in the subject. Summarizing, this article suggests a new way of approximately representing and processing probabilistic information in neural hierarchies which is quite different from what has been proposed in previous works, although it will be validated only for very simple stimuli. To be fair, on the other hand, most previous work on the subject uses test stimuli of similar simplicity with the possible exception of [START_REF] Rajesh | Bayesian computation in recurrent neural circuits[END_REF].

F. Goals of the article

As this article extends our previous work mainly by the aspects of Lyapunov analysis as well as the learning of lateral connection weights, we wish to show that • convergence towards the MAP estimate is plausible • learned lateral connections (empirically) produce the correct MAP estimates in simple cases • the space-latency code produced by learned connections is consistent with data statistics • pre-activation of field sites acts as a prior distribution in ambiguous cases

II. METHODS

Throughout the paper, we will make use of an extension of the dynamic neural field (DNF) model [START_REF] Erlhagen | Dynamic field theory of movement preparation[END_REF], [START_REF] Amari | Mathematical foundations of neurocomputing[END_REF], which represents two-dimensional sheets of rate-coded model "neurons" who have an internal state termed the "membrane potential" that evolves according to a differential equation taking into account the membrane potential, afferent synaptic input and lateral input mediated via a symmetric weight matrix termed "convolution kernel" or "interaction kernel". In its original formulation, the DNF model assumes infinitesimal model neurons, which allows analytical treatments of simple situations which would otherwise be impossible. In this article, all insights are gained through simulation experiments for which naturally a neural sheet has to be discretized.

A. Notation

We will consider, for simplicity of notation, the case of discretized fields in two dimensions having a size of (N, M ) elements. The set of membrane potentials in the whole field at a given instant shall be denoted u(t), a particular activity value u i (t) being identified by its one-dimensional index i. Please note that this is done for simplicity of notation: the field is considered to be two-dimensional. The dimension of the field is in fact irrelevant except in the definition of the lateral interaction, as all the other operations updating the field are applied in an element-wise fashion. The elementwise application of a non-negative, monotonous transfer function f to all field potentials, giving the firing rate vector z(t) shall be denoted as z(t) = f [ u(t)]. In this article, we will always a transfer function of the form

f (x) = 1 1 + exp -x-θ ν , (4)
with two free parameters, the threshold θ and the slope ν.

The application of a symmetric 2D convolution kernel K of size (N , M ) < (N, M ), represented by a vector of coefficients (k j ), j = 0 . . . N × M , to a vector of firing rates z(t) (producing the transformed version z ), is denoted by

z i = (K * z) i . (5)
The afferent input vector to a neural field is denoted by I, whose entries we suppose to be always positive and normalized in the interval [0, 1]. Using this notation, the classical DNF model evolves according to:

τ ui = -u i + I i + (K * f [ u]) i + h, (6)
h denoting the so-called resting potential.

B. Lyapunov analysis of the DNF model

Lyapunov analysis denotes the analysis of the asymptotic behavior of a dynamical system by means of a Lyapunov functional L, more commonly termed "energy functional" by analogy to physics. This is a functional that maps the current system state to a real number (the current "energy"): L [χ(t)] → R, and which has the property that it is being minimized by the system's dynamics. In addition, the Lyapunov functional is bounded from below, which means that at least one stable state exists, the ground state where L reaches the lower bound and thus cannot decrease any further, making the state stable. All these properties make Lyapunov analysis interesting for studying non-linear dynamic systems such as the extended DNF model which is proposed here, given of course that a Lyapunov functional exists, which is not self-evident and has to be rigorously proven for each dynamical system under consideration. Fur a more in-depth discussion of Lyapunov analysis, see [START_REF] Derek P Atherton | An introduction to nonlinearity in control systems[END_REF].

In [START_REF] Kubota | Anayzing global dynamics of a neural field model[END_REF], a Lyapunov functional for the DNF model was shown to exist, and an explicit expression for it was presented. In our notation, it has the following form:

L [ u(t)] = - i z i (I i -h) - 1 2 i z i (K * z i ) i + i f -1 (z i ) = = -z, I -h - 1 2 z, K * z + i F (z i ), (7)
where F (z) denotes a monotonous function of z which is derived from the transfer function f . As we can see, the first two terms can be restated as scalar products, whereas the third acts as a bound on total field activity. In qualitative terms, the evolution of a neural field thus aims simultaneously to achieve the highest possible similarity between firing rate and input, between firing rate and lateral input, all of this while maintaining a bounded overall activity.

C. Extension of DNF model

We now propose the following formal generalizations:

• Position-dependent global kernel: in order to learn complex data models as shown in Fig. 2, separate lateral interaction kernels must be learned for all neuronal sites. Furthermore, these local kernels must take into account all other field activities and not just adjacent ones, as statistical relations will not necessarily be local in the 2D plane. We therefore replace the convolution operation in eqn.( 6) by a simple linear transformation L, which amounts to a normal matrix multiplication. • Separation into excitatory and inhibitory lateral interactions: We wish to split the lateral interaction into a part that always gives positive results and which is potentially plastic, and a part that always acts in an inhibitory fashion and is fixed, ensuring a global attractor dynamics of the model. Therefore the lateral interaction term in eqn. ( 6) is split into a negative part with a constant coefficient γ, and a positive part depending on the position-dependent global kernel L. • Learning of excitatory lateral interactions: A suitable learning algorithm adapts the excitatory part of lateral interactions (parametrized by L) to the statistics of input data All of this motivates us to reformulate the DNF equation ( 6) as follows, additionally introducing the coefficients α, β and γ before key terms:

τ ui = -u i + αI i + β(Lf [ u]) i -γ i z i + h, (8)
where we introduce a symmetric linear transformation L = (l ij ) to a vector of firing rates z, producing the transformed version z of same dimensionality (the matrix L is therefore square and has dimensions (N * M, N * M )):

z i = (L z) i = j l ij z j ,
which does not in any way change the Lyapunov function of the model given in eqn. [START_REF] Gold | Neural computations that underlie decisions about sensory stimuli[END_REF] as the lateral interaction terms are effectively still described by a symmetric convolution kernel (the matrix L which is supposed symmetric, and the constant γ which is position-independent and thus symmetric by definition. As long as the learning rule for L therefore guarantees symmetry of L, according to [START_REF] Kubota | Anayzing global dynamics of a neural field model[END_REF], the Lyapunov functional stays well defined.

D. Learning of kernel matrices

We suppose that the task of the kernel is make fields converge, starting with an empty field, to the underlying true stimulus that has the highest probability given the input. This will depend upon statistical properties of the data which are encoded in the kernel matrices L by a learning process. Assuming an unsupervised learning scenario where the true stimuli T are unavailable, we suppose that the noise on inputs is such as to cancel by averaging over longer time scales. Thus, we are in an "autoencoder" type of learning scenario, where kernel matrices are adapted to make the field state reproduce the (average) input, see Fig. 3. In the experiments presented here, we do not even go this far as all training stimuli are "clean", and just test stimuli are noisy. In this case, learning supports pure disambiguation, i.e., deciding which training stimulus is closest to the noisy test stimulus. We therefore propose a simple linear regression model that adapts kernel values guided by the difference between neural field state (i.e., firing rates) and inputs, which we suppose normalized in the same interval [0, 1]:

l ij (t + 1) = l ij (t) + 2 z j [(L z) i -I i ] , (9)
. This is stochastic gradient descent approximating a meansquared-error loss function, which has the advantage of online learning capability: eqn. ( 9) is applied after each iteration of the field using a suitable value of in order to learn at the correct time scale.

III. MAP CALCULATIONS: HOW AND WHY IT WORKS

When performing Bayesian inference, the terms appearing in the likelihood and in the a priori distribution are usually of different quality. An example would be p(rain=yes|season=summer) ∼ p(season=summer|rain=yes) p(rain=yes), where the quantities "rain" and "season" live in different spaces that are not directly related. In contrast to this, in our case the true stimulus T and the input I live essentially in the same space, which does not in the least invalidate our approach as we can still define meaningful likelihoods and priors. What a neural field does is really a sort of auto-encoder process, trying to derive a less noisy version of its own input according to a learned data model. This could be of immense practical value if learning could efficiently be conducted in an unsupervised manner, which is exactly what we propose in Sec. II-D.

A. How

Supposing a kernel matrix L has been learned on a set of input samples following a certain distribution, we draw samples from the same distribution and put them as input to a neural field with kernel L (and parameters to be specified in more detail in Sec. IV). This input will evoke a timevarying firing rate z(t) that will converge after a certain time. Time-to-convergence can be measured in various ways, most simply by demanding that the sum over absolute values of the left-hand side of the extended DNF equation ( 8) be smaller than a threshold: As will be shown in the rest of this section, the converged attractor state represents the maximum a posteriori (MAP) estimate of the true stimulus T .

τ i | ui | < Ω (10)

B. Why

As we supposed that the kernel matrix L has been learned in such a way as to transform the current field activity z(t) into a superposition of possible true stimuli, z (t) = L(t) z(t). This is, as explained in Sec. II-D, because we assume that noise influences are weak and cancel each other over the course of the learning process. From the Lyapunov analysis in Sec. II-B we know that the field will evolve such that its firing rate resembles as closely as possible the input and z while maintaining a bounded activity. This process acts as a disambiguation between the possible true stimuli, represented as a superposition, only one of which will be selected due to global competitive dynamics.

Due to the nature of linear regression, which is sensitive to the relative frequency of training samples, the sum in z will already be weighted by the prior probability p( T = T ) of each true stimulus T . By applying preshape in the form of certain true stimuli, or by adding an attenuated version of the stimulus we wish to enhance to the input, we can manipulate these prior probabilities directly, although the exact amount of probability that is added by a certain amount of preshape is difficult to quantify for the moment.

IV. EXPERIMENTS

All experiments are conducted using the exact same parameters and differ just by the nature of the input stimuli. The parameters, as defined in Sec. II-C, are the following: α = β = 2, γ = 0.15, h = -1, θ = 0, ν = 2.5,τ = 25. To ensure numerical stability, field potentials are clipped to the range [-2, 5] after each iteration. The neural field itself has a size of 100x100 elements, and so do the inputs. Input patterns are always presented for N = 400 iterations, and fields are initialized to h before each pattern presentation. The lateral connection learning rate is set to = 0.001.

Common to all experiments is a training phase, where the three training patterns (see Fig. 4) are repeatedly presented, one after the other, for a total duration of 86400 iterations. As each pattern presentation lasts 400 iterations, this amounts to 216 pattern presentations, or 72 presentations for each pattern. Each test pattern (see Figs. 5,6) is presented once for N = 400 iterations. Inputs are superpositions of unnormalized Gaussians having the same variance of 3 "pixels". For simplicity, each Gaussian can appear at one of 5 predetermined positions, see Fig. 5 for details. The test stimuli used to perform experiments are depicted in Fig. 6. Some Gaussians are multiplied by a constant factor, reducing their amplitude in order to test a particular point. Results are presented in two forms: first of all, we wish to show the overall behavior of the field at all of the sites A-E. For this purpose, we simply plot the maximal value in a small region around each of these sites over time, which is shown in Fig. 7. Results for each pattern can be stated as follows:

• Stimulus 1: since the stimulus at site A was never observed during training (compare with the training stimuli depicted in Fig. 4, it can not create any lateral input to support itself and is thus strongly suppressed by global inhibition. Fig. 8. High temporal resolution of the development of field activities at "winning" sites for stimuli 3 (strong competition), 5(weak competition) and 7 (strong cooperation). It is evident that the time to convergence (latency) depends on the amount of competition in the field. This is exactly what we need to find if we wish to interpret latency as the posterior probability of the converged state.

stimuli who have same amplitude, and who are both supported by lateral input as they were part of the training stimuli. Competition occurs as they never appeared together during training, and thus did not form mutually supporting lateral connections. Competition finally suppresses one of the two Gaussians in a random fashion. By putting a very slight Gaussian preshape, i.e., a pre-activation of the field potential with amplitude A = 0.05 to one of the two sites, the randomness disappears and the site with the preshape always wins the competition. Preshape can thus effectively play the role of a prior probability distribution. • Stimulus 4: same as for stimulus 3 except that the Gaussian at site B is weaker than 1.0, so it has a competitive disadvantage and is suppressed. • Stimuli 5, 6: here, we observe competition between two stimuli as they were never seen together during training. however the stimulus at site C was seen during training and thus supports itself via lateral input, whereas the stimulus at site D does not, and is therefore quickly suppressed depending on its amplitude. • stimuli 7, 8: interaction between stimuli that were always seen together during training, and therefore have mutually supporting lateral connections. Not only does the weaker stimulus "survive", it is even boosted depending on its original strength.

We therefore observe that the field always converges towards the state that is evidently the most probable one given the training data, see Fig. 4. A second result concerns the continued validity of the space-latency coding now that lateral interactions weights are learned from data statistics.

To this end, we consider the stimuli 3, 5 and 7 from Fig. 6.

Here, we observe the following interactions between peaks: strong competition for stimulus 3, as both Gaussians are selfsupporting via lateral input, weak competition for stimulus 5 as one Gaussian can support itself and the other cannot, and strong cooperation for stimulus 7 as both Gaussians have bilateral exciting connections between them. Stronger competition is intuitively equivalent to lower posterior probability, as competition means that the field and the lateral input (the likelihood) do not match. This is indeed what we observe when we plot the time evolution of the field in response to these three stimuli, magnified so that small latency differences may be perceived (see Fig. 8): stimulus 3 has the largest latency due to strong competition, stimulus 5 has lower latency as competition is weaker, and stimulus 7 has the lowest latency since there is collaboration going on. Summarizing, at least for these simple cases, the spacelatency code continues to work correctly by expressing a high level of conflict in a field, corresponding to a low posterior probability, by increased latency. The effect is rather weak here, but by correct parameter choices it can be magnified to visible time scales.

V. DISCUSSION, CONCLUSIONS AND FUTURE WORK

Continuing our previous work on the subject, we have presented a neural model that gives a complete account of how neural systems may process and transmit probabilistic information in a way that is applicable to real-world processing. In particular, the presented model is capable of learning from observations, creating its own internal models of posterior distributions that allow Bayesian inference. We have made it very plausible, by means of Lyapunov analysis, that, if the model learns a true posterior distribution from data, the final attractor state will closely resemble a MAP estimate. Lastly, we have verified for some simple cases that the field indeed converges towards the MAP estimate, and that the latency of the attractor state is a measure of competition in the field, or conflict to the learned posterior distribution.

In Sec. I, three fundamental questions are posed: 1) How can distributions be represented on a neural level? 2) How can the multiplications required for manipulating them be performed, not to speak of the necessary renormalizations after each such step? 3) How can neural populations take optimal decisions based on distributions? In the light of the results and methods presented here, the answers we give are the following: first of all, no distributions are represented at the neural level, just an approximation given by their dominant mode and its associated probability (represented by latency). Then, as only the dominant mode of the resulting "multiplied" distribution needs to be represented by the following layer, multiplications are indeed unnecessary and can be replaced by a simple competitive method of determining the dominant mode while taking into account the associated probabilities. And lastly, decision making occurs when competitive dynamics "selects" the dominant mode and suppressed the others.

It may be criticized that we give up the idea of a fully "Bayesian brain" in favor of an approximation that involves the MAP estimate and its posterior probability. However, in a behavioral context, it is often only the most likely interpretation T * of an input that is required to take an action. Furthermore, sub-leading interpretations are just suppressed but continue to exist, and can resurface if additional external input is provided. In addition, we also compute a measure of the posterior probability associated with the MAP estimate, which effectively realizes an unimodal approximation of the posterior distribution. If this probability is low, we know that the other modes of the distribution are non-negligible and we would do well not to take the MAP estimate too seriously. This manner of thresholding and, if need be, of discarding information is very relevant for intelligent agents, as it is in general better to wait for additional inputs than to take decisions based on inputs that are known to be inconclusive.

Although we do not present a spiking neural network model, we nevertheless claim high biological plausibility for the proposed model. Just as the original DNF model, it is formulated at a lower spatial "resolution" and thus not on the single-neuron level but rather on the population level (e.g., minicolumns or macrocolumns). In this context, introducing spiking models would be very inconsistent and is therefore avoided. Moreover, what we primarily aim to model is the dynamical interplay of lateral excitatory and inhibitory influences and its computational significance, especially since it is well-known that lateral connections are subject to learning as well [START_REF] Risto Miikkulainen | Computational maps in the visual cortex[END_REF]. In this sense, we feel it is well justified to claim biological plausibility for our model.

It is very evident that the next step must be a rigorous proof of the MAP property. Or rather: it must be shown that the lateral connection weights learn a posterior distribution (which will depend largely on the used learning rule). Subsequently, it must be established to what extent, and to what precision, the field converges to the maximum of this distribution. it would be very helpful if, at the same time, some statements could be made about the time to convergence (i.e., the latency) and find a rigorous link to probability here as well.

Fig. 3 .

 3 Fig. 3. Schematic sketch of the learning of kernel matrices, i.e., lateral connection weights. Basically, weights are adapted so that the difference between field activity (center,top) and input (left,bottom) is minimized. Modified kernel matrices act upon field activities in a non-linear way via lateral inputs (right,bottom).

Fig. 4 .

 4 Fig. 4. 100x100 pixel input patterns used for training the lateral connection weights of the model (axes are scaled differently).

Fig. 5 .

 5 Fig. 5. Discrete positions of stimuli in the test patterns used for subsequent experiments. At each discrete position, a Gaussian stimulus can appear, either in isolation or in combination with Gaussians at other positions.

Fig. 6 .Fig. 7 .

 67 Fig.6. The eight 100x100 pixel test stimuli used in the experiments, consisting of superpositions of Gaussians (axes are scaled differently). Gaussian amplitudes are 1.0 where not otherwise stated. In stimulus 4, the Gaussian at site B has amplitude 0.9. In stimulus 6, the Gaussian at site D has amplitude 0.8. In stimulus 8, the Gaussian at site D has amplitude 0.5.

•

  Stimulus 2: the stimulus at site B was observed during training, leading to a lateral response that sustains activity despite inhibition.• Stimulus 3: here, we observe competition between two

•

  Neural activities represent MAP estimates: given an afferent input, we assume that neural dynamics, guided by lateral connections, converge to a MAP estimate given the input. Neural activities are therefore not at all related to probability distributions. The basic functionality achieved by the model is depicted in Fig.2: following a period of training with clean stimuli, the model receives a noisy stimulus and attempts to converge to the noise-free ("true") underlying stimulus based on the results of the training phase. In this contribution, we shall make it very plausible (by considering the Lyapunov functional of the model) that the model in fact converges to the maximum a posteriori (MAP) estimate for the given

	0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8	0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

• Weights represent distributions: for ensuring convergence towards MAP estimates, the lateral connections must represent in some way a posterior distribution of the data • posterior distributions are learned: as a model with pre-specified posterior distributions would be useless in practice, we propose a simple learning scheme for adapting the lateral connections to data statistics

We would be delighted about rigorous mathematical proofs on this subject


Cem Karaoguz

Incremental

Incremental Learning for Bootstrapping Object Classifier Models

HAL is

Incremental Learning for Bootstrapping Object Classifier Models

Cem Karaoguz 1,2 and Alexander Gepperth 1,2 Accepted for 2016 IEEE 19th International Conference on Intelligent Transportation Systems (ITSC 2016) Abstract-Many state of the art object classification applications require many data samples, whose collection is usually a very costly process. Performing initial model training with synthetic samples (from virtual reality tools) has been proposed as a possible solution, although the resulting classification models need to be adapted (fine-tuned) to real-world data afterwards. In this paper, we propose to use an incremental learning from cognitive robotics, which is is particularly suited for perceptual problems, for this bootstrapping process. We apply it to a pedestrian detection problem where a synthetic dataset is used for initial training, and two different realworld datasets for fine-tuning and evaluation. The proposed scheme greatly reduces the number of real-world samples required while maintaining high classification accuracy. We also demonstrate an innovative incremental learning schemes for object detection which training object and background samples one after the other: this keeps models simple by representing only those background samples that can actually be confused with pedestrians.

I. INTRODUCTION

Applications in domains of driver assistance systems and autonomous driving demand visual recognition of traffic participants invariant to changes like illumination, viewpoint etc. To cope with this, modern object classification systems based on statistical learning require datasets with large number of annotated samples recorded under different conditions. Building up such large datasets is usually a tedious process and comes with high costs in resource and time. Furthermore, changes in various aspects of the application (e.g. hardware) might require construction of a modified database and retraining of models from scratch. This retraining "from scratch" is necessary since none of the methods usually used for classification have what is termed incremental learning capacity (see [START_REF] Gepperth | Incremental learning algorithms and applications[END_REF] for a discussion of the term) that would allow to update models with new samples without complete retraining, and without "damaging" already learnt knowledge.

In this context, a very popular setting for activities like domain adaptation [START_REF] Ben-David | A theory of learning from different domains[END_REF] is to consider model training on a source database with easy-to-obtain synthetic samples, e.g., from a virtual reality tool, and then to adapt (or fine-tune) the model to a target database with hard-to-obtain real-world samples. The beauty of this bootstrapping approach is that the number of hard-to-obtain samples required for fine-tuning is usually far inferior than the number required for training Talence, France cem.karaoguz@ensta.fr, alexander.gepperth@ensta.fr from scratch, and thus a great increase in efficiency can be achieved.

In this article, we show that dedicated incremental learning algorithms, as proposed in the domain of machine learning and developmental robotics, can be used as a ready-made tool to greatly facilitate the bootstrapping process for perceptual tasks in intelligent vehicles. In particular, we use the incremental learning approach presented in [START_REF] Gepperth | A bio-inspired incremental learning architecture for applied perceptual problems[END_REF] which is particularly suited for high-dimensional perceptual problems, for implementing the bootstrapping approach between a synthetic and several real-world databases for pedestrian classification. The bootstrapping is done in two phases: with positive (i.e. object) data first, and subsequently with negative data (i.e. background). Positive bootstrapping aims to form a coarse representation of an object using the source database and adapt it to a target database with only a limited number of examples. Negative bootstrapping aims at filtering out simple negative samples from the training process using already learned object representations so that only hard negative samples are represented. This scheme eliminates the necessity of multiple rounds of training employed by many state of the art object classifiers as well as leaves more resources in the model to represent object characteristics.

When performing experiments using synthetic data from [START_REF] Vazquez | Virtual and real world adaptation for pedestrian detection[END_REF] and real-world data from the KITTI and Daimler databases [START_REF] Geiger | Vision meets robotics: The kitti dataset[END_REF], [START_REF] Enzweiler | Monocular pedestrian detection: Survey and experiments[END_REF], we show that only a few annotated real samples are enough to sufficiently adapt models to the (slightly different) statistics of real samples. Hence, the proposed framework significantly reduces the number of real images and corresponding annotations required by the model training process and renders models reusable across different datasets, eliminating the necessity of model re-training from scratch.

A. Related Work

The related work can be addressed in two main groups: incremental learning and domain adaptation. A common strategy for incremental learning is to partition the input space and use local models for each partition. This avoids common problems of machine learning like catastrophic forgetting [START_REF] Goodfellow | An empirical investigation of catastrophic forgeting in gradient-based neural networks[END_REF] or concept drift [START_REF] Kulkarni | Incremental learning from unbalanced data with concept class, concept drift and missing features: a review[END_REF] since learning is always localized in the input space, in the sense that a change of statistics in one part of that space will not affect learning in other, distant parts. The manner of performing this partitioning is very diverse, ranging from kd-trees [START_REF] Cederborg | Incremental local online gaussian mixture regression for imitation learning of multiple tasks[END_REF] to genetic algorithms [START_REF] Butz | Computational complexity of the xcs classifier system[END_REF], adaptive Gaussian receptive fields [START_REF] Vijayakumar | Locally weighted projection regression: An o(n) algorithm for incremental real time learning in high-dimensional spaces[END_REF]. Equally, the choice of local models varies between linear models [START_REF] Vijayakumar | Locally weighted projection regression: An o(n) algorithm for incremental real time learning in high-dimensional spaces[END_REF], Gaussian mixture regression [START_REF] Cederborg | Incremental local online gaussian mixture regression for imitation learning of multiple tasks[END_REF] or Gaussian Processes [START_REF] Nguyen-Tuong | Local gaussian processes regression for real-time model-based robot control[END_REF]. The choice has to be made regarding the constraints on the computational complexity imposed by the application.

Bootstrapping object classifiers is closely related to domain adaptation problem where an object classifier trained on a source dataset needs to be operated reliably on a different target dataset. Two major approaches exist to face this problem: feature transformation and model adaptation. Feature transformation relies on projecting feature vectors to a space compatible with the classifier of the source domain (e.g. [START_REF] Hoffman | Continuous manifold based adaptation for evolving visual domains[END_REF], [START_REF] Kulis | What you saw is not what you get: Domain adaptation using asymmetric kernel transforms[END_REF]). On the other hand, model adaptation is based on adjusting the parameters of an already learned model or learning complementary models to cope with changing data statistics. This is also the approach taken in the current work. An incremental domain adaptation framework is presented in [START_REF] Xu | Incremental domain adaptation of deformable part-based models[END_REF] where two separate classifiers are used. A linear combination of domain and target classifiers gives the final classification result and the weight of each classifier is determined by their recorded performance. In [START_REF] Jain | Online domain adaptation of a pre-trained cascade of classifiers[END_REF] a Gaussian process regression model is constructed from confident outputs of a classifier and the scores of data instances with low prediction values are modified by this model. A-SVM is introduced in [START_REF] Yang | Cross-domain video concept detection using adaptive svms[END_REF] which enables domain adaptation for SVM based classifiers by learning a perturbation function between source and target classifiers. This idea is extended to cope with multiple target domains by hierarchically organizing these target domains in [START_REF] Xu | Hierarchical adaptive structural svm for domain adaptation[END_REF]. The majority of these works are based on discriminative models where direct adaptation of the model to changing data statistics is problematic. Hence, these methods often train new models (target models) on top of the existing ones (source models) or learn a residual function i.e. statistical difference between datasets. In contrast, in the presented work models are updated directly and continuously in the presence of new data. Hence, the approach is generic and works without prior knowledge about datasets.

II. METHODS

The proposed architecture is a three-layer neural network which is illustrated in Fig. 1. Adopting the common notation of neural networks, we utilize superscripts I , H and O for entities related to input, hidden and output layers, respectively. The input layer of the network is composed of a feature vector which is generated from the input data. The hidden layer of the network projects the input vector onto the prototype space based on a distance metric. Subspaces of the input space are coarsely approximated by hyperspheres whose centers are defined by the prototypes in the hidden layer. The output layer is composed of allto-all connected neurons that map local input space regions (i.e., sets of prototypes) to class memberships using simple linear regression learning.

A. Projection

The hidden layer of the network is composed of topologically organized prototypes represented as weight vectors w H m where w H ∈ IR N×M . Prototypes are distributed in a two dimensional grid (see Fig. 1), hence prototype locations are indicated as vectors m. However, we drop the vector notation for brevity and simply use m which can also be interpreted as prototype ID. The hidden layer acts similar to the well-known self-organizing map (SOM) algorithm [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF]: the projection of the input onto the hidden layer starts with computing the distance between the input vector and all prototypes:

zH (m) = ||w H m -z I || (1)
where z I is the input vector, || • || is the Euclidean norm.

The prototype m * with the smallest distance is called the best matching unit. In our model, the hidden layer re-encodes the input in a way that enables incremental learning while preserving information. Therefore, instead of reducing the output of the hidden layer to the best-matching unit (as it is usually done for SOMs), we calculate the (graded) activations of all hidden layer units:

zH = g κ zH (2)
where the activation function g κ is Gaussian with standard deviation κ. The activation function converts the distance measures into similarity and keep them in the [0, 1] interval. A transfer function is further applied to sparsify these activations:

z H = TF p zH (3)
where TF(•) represents a monotonous non-linear transfer function, TF : [0, 1] → [0, 1] which maintains the best matching unit value unchanged while non-linearly suppressing smaller values:

TF p (z H ) = zH p (z H (m * )) p-1 (4)

B. Prediction

Hidden layer is connected to the output layer in all-toall fashion with weights w P ∈ IR M×C . Generation of output layer activities is performed by a simple linear transformation of hidden layer activities z H :

z O (m) = w O m • z H (5)
The class associated with the unit that has the strongest activity in the output layer becomes simply the predicted class if the activity exceeds a threshold.

C. Learning model parameters

Prototype adaptation is performed online using the conventional SOM update step except that it takes into account a control signal λ coming from the output level of the hierarchy:

w H m ← w H m + λ H g σ (||m -m * ||)(z I -w H m ) (6)
where g σ (x) is a zero-mean Gaussian function with standard deviation σ. The control signal λ is a binary value that is set to 1 only if the current estimate of class membership, i.e., the output layer activities z P is either uncertain or wrong. The uncertainty is measured from the bounded difference between first and second maximum of the output layer The overview of the proposed architecture applied to pedestrian detection task. In this example, image intensity values are used directly as feature vector. The input vector is projected on the topologically organized prototype space that involves comparison with all prototypes. This renders an activation map of prototypes where prototypes similar to the input vector give higher values. At the prediction step the resulting activation map is mapped to a class-membership vector via regression.

activities. If the difference is below a threshold θ m the control signal is set to 1. In accordance with standard SOM training practices, the SOM learning rate and radius, H and σ, are maintained at 0 , σ 0 for t < T 1 iterations and are exponentially decreased afterwards in order to attain their long-term values ∞ , σ ∞ at t = T conv .

Since the output layer performs linear regression, the weights are modified via online gradient descent, optimizing the mapping of hidden layer activities z H to the target representation z T containing the "true" class of a sample:

w O m ← w O m + 2 O z H z O (m) -z T (m) (7)
In contrast to the hidden layer learning rate, the learning rate of linear regression, P remains constant at all times.

D. Incremental Learning for Bootstrapping

We employ incremental learning for bootstrapping of models in two phases: Positive bootstrapping aims that a basic prototype structure is learned with the source dataset. Prototypes are later fine-tuned with the target dataset where necessary. In order to realize this, the system is exposed to positive samples from the source dataset for T bsp+ iterations. After this, fine-tuning is performed by exposing positive samples from the target dataset to the system. The incremental learning scheme outlined in Sec.II-C ensures that weights are adapted only if the system cannot correctly classify a given sample. Hence, prototypes which describe the target dataset already sufficiently are not touched.

Negative bootstrapping is done subsequently and follows a scheme similar to the positive bootstrapping with additional heuristics. The prototype-based representation previously trained on positive samples enables early rejection of negative samples before actual classification due to the generative nature of the model: Eqn. (2) computes the activation of an input vector, negative samples yielding low activation values in this process if they are dissimilar to objects. Setting a threshold θ bsp can identify such samples in order to leave them out from further computations, i.e. classification and learning. This has two benefits: first, the model only represents negative samples which are very similar to the objects. With minimum allocation of the model's resources to negative samples, more representational resources (prototypes) can be allocated for representing the object class. Second, the system can use its own output directly during tests with annotated images without having to crop and prepare negative samples beforehand.

The training schedule employed in experiments (see Fig. 2) starts with positive bootstrapping for T bsp+ iterations, followed by fine tuning with positive samples for T f t+ iterations and finalized by negative bootstrapping for T bsp- iterations.

III. EXPERIMENTS A. Setup and Parameters

Three different datasets are used in the experiments: for synthetic dataset, the Virtual Pedestrian dataset presented in [START_REF] Vazquez | Virtual and real world adaptation for pedestrian detection[END_REF] is used whereas for real-world ones the KITTI Vision Benchmark Suite [START_REF] Geiger | Vision meets robotics: The kitti dataset[END_REF] and the Daimler Mono Pedestrian Detection Benchmark Dataset [START_REF] Enzweiler | Monocular pedestrian detection: Survey and experiments[END_REF] are used. Sample images from these datasets are shown in Fig. 3, and the number of samples used in the experiments are shown in Tab. I for each dataset. Experiments are conducted in the following settings:

• Baseline: measures the performance of the system on the source dataset. In this setting, both training and test are done in the same dataset. [START_REF] Dalal | Histograms of oriented gradients for human detection[END_REF] to be used as input vectors. This method is chosen due to the availability of the feature extraction framework however, any vectorized representation of images can be used with the system. The number of iterations is set to T bsp+ = 500000, T f t+ = 300000 and T bsp-= 700000.

For comparison to the state of the art, a discriminative approach is also implemented and evaluated to assess how close training on the synthetic dataset can get to training on real-world datasets, and whether it is feasible to be used as a basis for bootstrapping. We chose a state-of-the-art boosting algorithm for the purpose, which is widely used in discriminative object detection (see [START_REF] Dollar | Pedestrian detection: An evaluation of the state of the art[END_REF] for a survey). Our implementation follows [START_REF] Dollár | Fast feature pyramids for object detection[END_REF] with the number of weak classifiers set to 1000 which are selected from a pool generated using the aggregated channel features explained in the paper.

(a) Samples from the synthetic dataset [START_REF] Vazquez | Virtual and real world adaptation for pedestrian detection[END_REF] (b) Samples from KITTI dataset [START_REF] Geiger | Vision meets robotics: The kitti dataset[END_REF] (c) Samples from Daimler dataset [START_REF] Enzweiler | Monocular pedestrian detection: Survey and experiments[END_REF] Fig. 3: Representative samples from datasets used in experiments

B. Cross-Database Performance of Discriminative Models

Fig. 4 shows the performance of Adaboost classifiers trained and tested on various datasets. Boosting-RonR is the classifier trained and tested on KITTI dataset hence it serves as a baseline for the real-world dataset. Similarly, Boosting-SonS is trained and tested on synthetic dataset and serves as a baseline for this dataset. Boosting-SonR is the classifier trained on the synthetic dataset and tested on KITTI dataset. The results suggest that synthetic data has indeed statistical drift w.r.t. the real-world one. However, the extent of this drift is still below a margin that would allow us to use the synthetic dataset for positive bootstrapping purposes. The results are also in accordance with [START_REF] Marin | Learning appearance in virtual scenarios for pedestrian detection[END_REF] where a similar study is done for a comparable synthetic dataset and Daimler dataset using SVM classifiers.

C. Effects of Positive Bootstrapping

Firstly, we examine the performance of the system after positive bootstrapping only. This corresponds to the state of the models achieved after box 1 in Fig. 2 is processed. Results shown in Fig. 5 indicate that already at this stage, it is possible to filter out half of the negative samples and still achieve around 1% of miss rate on both KITTI and Daimler datasets. Adding the negative bootstrapping on top of this (i.e. after both box 1 and box 2 in Fig. 2 are processed) greately reduces the amount of false positives as shown in Fig. 6. Highest detection rates achieved after the full bootstrapping process are summarized in Tab. II. Compared to discriminative model (Fig. 4), false positive performance of the generative models does not vary extensively. Cross-Database performance of the model trained on the synthetic dataset varies. When applied to Daimler dataset, this model can achieve a performance close to the setting Baseline-Small whereas on KITTI, the performance is the lowest among all. Bootstrapping improves the models in both cases: the detection performance is increased from 62% to 88% on KITTI and from 92% to 95% on Daimler. On both datasets, the Bootstrapping setting achieves better results than Baseline-Small indicating that virtual dataset indeed provides the models with a useful knowledge basis.

D. Effects of Negative Bootstrapping

As shown earlier, a generative model trained with only positive samples can filter out "easy" negative samples: for a sample, if no similar prototype is found, the sample can be regarded as negative. Hard negatives on the other hand, are used in the training. This process is referred to as negative bootstrapping. In order to analyze the effects of negative bootstrapping, we keep track of the negative samples that are rejected during the training process. For both Daimler and KITTI datasets around 22% of the negative samples are filtered out by the models during the whole negative bootstrapping process. This is not to be confused with the performance of the system where models are trained with positive bootstrapping only (Fig. 5). During negative bootstrapping, as hard negative samples are accepted by the system the generative model starts building their explicit representations (i.e. prototypes). As the training process advances these prototypes accept more negative samples, which similar to them, reducing the total number of negative samples filtered out by the end of training. Fig. 7 shows some examples of the accepted and skipped negative samples from Daimler dataset. It can be noticed that the system not only eliminates easy samples such as sky, but also samples with texture. On the other hand, samples with high texture seem to be included in the training process. These may generate feature vectors that can confuse the models depending on texture patterns e.g. specific constellation of tree branches and leaves or shadows cast on the road surface.

IV. CONCLUSION

We present an object classification architecture that combines generative and discriminative models that facilitates incremental learning. The main contribution of this work is the utilization of the incremental learning approach for bootstrapping of object classification models. Bootstrapping is practised with both positive and negative samples. The essence of positive bootstrapping is training the generative model with positive samples from one dataset followed by a fine tuning process with positive samples from another dataset. The incremental learning capability of the system allows to make local changes in the model hence, the knowledge acquired from the first dataset is retained if it is useful for the second and statistical properties of the second dataset, which is not covered, is incorporated into the model gracefully. The feasibility of the system is demonstrated on pedestrian detection task. The positive bootstrapping is employed with synthetic data from the synthetic Pedestrian dataset and tested on a real-world dataset. It is possible to further improve the performance to the baseline level via the fine tuning step with only a few labeled examples from the real-world dataset. This is especially useful for applications where the amount of annotated data is limited. Cheaper and more convenient synthetic data can be used to bootstrap the models and state of the art performance can still be achieved.

In addition to positive bootstrapping, we also proposed negative bootstrapping where the system can reject a portion of negative samples using the internal object representations built by positive training. Incremental learning yields the benefit of updating only the parts of the model where the positive/negative discrimination stays weak. In this case, the representation of negative samples in the model is kept at minimum, allowing more resources for positive samples for a better object representation. This also eliminates the necessity of employing conventional bootstrapping methods where the models are initially trained with positive samples, false positives are collected as hard negative samples via the learned models and models are re-trained with positive and hard negative samples. In the proposed approach, the system is trained with positive and negative samples sequentially. Initial training only with positive samples already builds object representations and later at test stage only hard negatives are determined by the system and incorporated in the training. This property, combined with low computational complexity of the models and GPU paralellization renders relatively short training times. For the work presented, the whole training process takes less than 30 minutes.

  Fig.1: The overview of the proposed architecture applied to pedestrian detection task. In this example, image intensity values are used directly as feature vector. The input vector is projected on the topologically organized prototype space that involves comparison with all prototypes. This renders an activation map of prototypes where prototypes similar to the input vector give higher values. At the prediction step the resulting activation map is mapped to a class-membership vector via regression.

Fig. 2 :

 2 Fig. 2: Temporal sequence of training steps during the bootstrapping process from a source dataset (usually synthetic) to a target dataset (usually real-world, here: KITTI). The first step of positive bootstrapping includes training on positive samples from the source dataset, then training on positive samples from the target dataset. The subsequent negative bootstrapping is conducted in an analogous fashion except that it excludes "easy" negative samples.

•

  Baseline-Small: measures the performance of the system on the source dataset but training is done with a much smaller sub-set of positive samples chosen randomly. This setting is a control scenario to analyze if adapting models with few examples after bootstrapping does not converge to the same solution as learning models with few examples. • Cross-Dataset: performs training on source dataset using all data available while testing on a target dataset. • Bootstrapping: performs training on the source dataset and incremental model updated with limited positive data from target dataset as explained in Sec. II-D. For bootstrapping the same number of samples shown in Tab.I are used. For update smaller number of samples are used. The number of positive samples for the settings Baseline-Small and Bootstrapping is set to 50. We use the following fixed parameters for our system: the number of prototypes in the hidden layer: M = 20×20 = 400, O = 0.001, 0 = 0.1, σ 0 = 6, ∞ = 0.001, σ ∞ = 1, θ m = 0.7, T 1 = 50000, T conv = 150000, p = 20 and τ = 0.001. Both SOM and LR weight matrices are initialized to random uniform values between -0.001 and 0.001. Training examples are always randomly and uniformly drawn from the current training set. The Histograms of Oriented Gradients features extracted from samples following

Fig. 4 :

 4 Fig. 4: Results of experiments with boosting models on KITTI as the real-world dataset.

Fig. 5 :

 5 Fig. 5: Performance of the systems trained with only positive bootstrapping (synthetic and real-world) on KITTI and Daimler datasets.

Fig. 6 :

 6 Fig. 6: Performance of the systems trained under various settings.

Fig. 7 :

 7 Fig. 7: Some of the negative samples included and excluded by the system during training on Daimler dataset.

TABLE I :

 I Number of samples used in experiments

		Training	Test
		Positive Negative	Positive Negative
	Synthetic	1700	1000	716	2000
	KITTI	1000	10000	532	4968
	Daimler	5000	5000	10660	8488

TABLE II :

 II Detection performance obtained in experiments

		Detection Rate	FPR
	KITTI Baseline	0.94 0.025
	KITTI Baseline-Small	0.86 0.030
	Cross-Dataset (KITTI)	0.62 0.030
	Bootstrapping (KITTI)	0.88 0.035
	Daimler Baseline	0.97 0.015
	Daimler Baseline-Small	0.91 0.015
	Cross-Dataset (Daimler)	0.92 0.015
	Bootstrapping (Daimler)	0.95 0.015

ACKNOWLEDGMENT

We gratefully acknowledge the support of NVIDIA Corporation with GPU donation for this research.


Cyril Allignol

Nicolas Barnier

Nicolas Durand

Guido Manfredi

Éric Blond

Integration of UAS in Terminal Control Area

Keywords: UAS, self-separation, conflict resolution, geometrical algorithm

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

I. INTRODUCTION

In the last years, the demand to operate civilian UASs for very diverse missions (ranging from fire detection and river bed surveillance to small items delivery and archaeological exploration...) has rocketed [START_REF] Abiresearch | Small unmanned aerial systems market exceeds US$8.4 billion by 2019, dominated by the commercial sector and driven by commercial applications[END_REF]. Most of these UASs operate within the lower airspace (under FL180) and their trajectories may interfere with commercial aviation traffic in Terminal Control Areas (TCA). Developing new algorithms to help separating UASs from the rest of the traffic is becoming critical for safety reasons. In this context, several approaches could be investigated.

First, separation can be entirely managed by Air Navigation Service Providers. Experiments [START_REF] Dcf | ODREA demonstration report[END_REF] have shown that Air Traffic Controllers resolution process is hindered by UASs mixed with conventional traffic because they have unusual performance specifications and interact with different time responses.

Otherwise, separation could be delegated to both commercial aircraft and UASs which could autonomously maneuver to resolve potential conflicts. However, complex processes of coordination should be considered in such a context to keep Air Traffic Controllers aware of the resolution process and able to interfere in it.

Finally, conflict resolution could be taken care of by UASs only such that they do not disrupt the commercial traffic. This seems a more realistic approach, provided that the positions and speeds of surrounding aircraft are available (through ADS-B for example) and that the performances of UASs and the resolution anticipation are sufficient to solve all traffic situations. In this article, we investigate the latter by adapting a self-separation algorithm used in robotics to our context and experiment with various parameters and strategies on real traffic samples.

1) Detect & Avoid: This algorithm was designed by van den Berg et al. [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF] and tested with different speed constraint hypotheses by Durand et al. [START_REF] Durand | Does ATM need centralized coordination? Autonomous conflict resolution analysis in a constrained speed environment[END_REF] in the context of autonomous air conflict resolution. Compared to these first experiments, we here tailor the geometrical approach of van den Berg further to model the performances of UASs and consider specific fallback strategies to handle cases for which the first approach fails to maintain separation.

More specifically, adapting the algorithm to integrate UASs in commercial traffic leads to the following major differences from the preliminary work presented in [START_REF] Durand | Does ATM need centralized coordination? Autonomous conflict resolution analysis in a constrained speed environment[END_REF]:

• The UAS is here supposed to avoid other aircraft, which do not avoid the UAS. This means that the whole avoidance maneuver is endorsed by the UAS. • UASs used in a civilian context generally fly with low speeds compared to commercial aircraft. The ratio we used in this article can go from 1.5 to 5. We focus on the lower airspace where the aircraft speed is theoretically limited at 250 kn, but recorded data show that in practice some aircraft fly much faster (up to 400 kn). In this study, we consider two types of UASs: Fast UASs flying at 160 kn and Slow UASs flying at 80 kn. • Moreover, most civilian UASs have very poor speed up performances compared to conventional aircraft. We will therefore only consider maneuvers at constant speed for UASs, as this degree of freedom would have almost no effect on the resolution process with realistic traffic. • Commercial aircraft flying in the lower airspace are generally climbing or descending and their speeds are constantly changing, either increasing when climbing, or decreasing when descending and changing direction as well. This factor has a great influence on the detect and avoid strategy in order to ensure that a reasonable distance to the encountered traffic can be maintained. Using real traffic data is therefore essential to validate a resolution algorithm for such evolving and intricate traffic.

• Air traffic trajectory prediction, which is one of the main component of a conflict solver, is always tainted with uncertainties that must be taken into account to assess the efficiency of an algorithm. We show how our approach can handle uncertainties by providing robust resolution maneuvers. 2) Related Works: When the concept of Free-flight emerged in the 90s, one of the ideas was then to equip every aircraft with a detect and avoid algorithm able to ensure separation with the rest of the traffic.

The first effective approach used sliding forces to coordinate maneuvers between aircraft [START_REF]A comparison of different approaches based on force fields for coordination among multiple mobiles[END_REF]. Potential or vortex fields [START_REF] Kosecka | 2-1/2 D Conflict Resolution Maneuvers for ATMS[END_REF] as well as a model based on an analogy with electrical particle repulsion [START_REF] Eby | Free flight separation assurance using distributed algorithms[END_REF] were also used. In 2001, we proposed a token allocation strategy combined with an A * algorithm to solve conflicts with realistic maneuvers [START_REF] Granger | Token allocation strategy for free-flight conflict solving[END_REF], [START_REF] Granger | Optimal resolution of en route conflicts[END_REF]. Even if some maneuvers could be simultaneously decided, a complete ranking of aircraft was necessary and finding an optimal ranking has been shown to be problem-dependent [START_REF] Archambault | Scheduling heuristics for on-board sequential air conflict solving[END_REF]. We also tried artificial Neural Networks on the two-aircraft problem [START_REF] Durand | Neural nets trained by genetic algorithms for collision avoidance[END_REF] but they could not be generalized to handle more aircraft. All these approaches have been tested on en-route traffic, mainly with leveled aircraft.

Geometrical algorithms have also been widely studied in robotics [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF], [START_REF] Hwang | Protocol-based conflict resolution for air traffic control[END_REF]- [START_REF] Snape | Navigating multiple simple-airplanes in 3D workspace[END_REF]. The powerful technique developed by Van den Berg et al. [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF] can handle thousands of agents in a small space. It was applied to aircraft by Snape et al. [START_REF] Snape | Navigating multiple simple-airplanes in 3D workspace[END_REF], but the hypotheses of the algorithm require simultaneous vertical and horizontal speed changes. We also tested them [START_REF] Durand | Does ATM need centralized coordination? Autonomous conflict resolution analysis in a constrained speed environment[END_REF] in the horizontal plane with speed constraints and showed that this algorithm is unable to deal with high densities of traffic when the speed norm cannot be changed.

More generally, conflict resolution has been proven to be a highly combinatorial optimization problem [START_REF] Durand | Automatic aircraft conflict resolution using genetic algorithms[END_REF]. Most centralized approaches that have been proposed to solve conflicts can be broadly divided into two main categories. The first ones [START_REF] Krella | Arc 2000 scenario (version 4.3)[END_REF]- [START_REF] Hu | Optimal coordinated maneuvers for three dimensional aircraft conflict resolution[END_REF] use greedy sequential algorithms to optimize trajectories one by one after ranking the aircraft (ordering aircraft is however very challenging [START_REF] Archambault | Scheduling heuristics for on-board sequential air conflict solving[END_REF]). The others try to find the global optimum without the need to prioritize aircraft. Among this second category, many models define aircraft trajectories through simple analytic expressions that introduce strong limits on the type of situations that can be dealt with, as the ones described in [START_REF] Oh | Design and analysis of conflict resolution algorithms via positive semidefinite programming [aircraft conflict resolution[END_REF]- [START_REF] Gariel | 3D conflict avoidance under uncertainties[END_REF]. In [START_REF] Durand | Optimal resolution of en route conflicts[END_REF], [START_REF] Allignol | A New Framework for Solving En-Route Conflicts[END_REF], we proposed a model to solve multiple aircraft conflicts based on Metaheuristics (Genetic Algorithm and Tabu Search) using trajectory simulation with uncertainties. However, these works mainly targeted en-route traffic control and used simulated traffic only with the BADA model on real flight plans.

3) Outline: In this article, we come back to a simpler problem in a more realistic environment. We consider UASs flying in the lower airspace (under FL180) and design various conflict scenarios with real recorded commercial aircraft trajectories in TCAs. The aim of the study is to assess the performances of a "detect and avoid" strategy for UASs to maintain a reasonable horizontal separation with commercial traffic. We first consider conflicts involving one UAS and one commercial aircraft only, then we introduce multiple aircraft encountering the UAS.

A B d d/τ # » v A # » v B #» v r
In section II of this article, we detail the geometrical approach developed by Van den Berg et al. [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF] in the case of a single UAS avoiding non cooperative aircraft. We show how the algorithm can be adapted to take into account the constant speed constraint and propose several strategies to choose the maneuver when a future conflict is detected. We also introduce a fallback scheme that reduces the targeted separation distance when no solution is available to maintain separation. Moreover, we describe how our tailored model can be easily extended to take more than two aircraft into account. Finally, we show how some uncertainties on trajectory prediction can be simply taken into account in our model. In section III-A, we describe how we built the conflict scenarios from real traffic data and the hypotheses that were chosen for the UAS. In section III-B, we give some results obtained on the different scenarios and show the influence of the various parameters on the quality of the results in terms of separations achieved and maneuvers imposed on the UAS. The last section draws some conclusions about the results obtained with the simulations and highlights directions for future work.

II. DETECT AND AVOID MODEL

This section describes the algorithm developed in [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF] and its adaptation to the case where only one aircraft (the UAS) maneuvers. We first detail the model with a single aircraft and a minimizing strategy that requires to change the speed of the UAS, then show how the speed of the maneuvering UAS can be constrained to have a constant norm at the cost of greater heading deviations. We also describe a fallback strategy to minimize the separation violation when our model has no solution within the allowed turning range. Our model is then further extended to simultaneously take into account several aircraft while keeping the UAS speed constant.

A. Single Aircraft

Let d be the standard separation between aircraft and τ be a look ahead time. In figure 1, let us consider UAS A and When a potential conflict occurs, i.e. the endpoint of #» v r lies in the forbidden zone, #» v r must be changed to escape it. If the speed change is not constrained (by the performances of the aircraft for example), #» v r can be projected on the nearest bold line (see figure 2) to solve the conflict pairwise between UAS A and aircraft B while minimizing the speed change. In [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF], this change is shared in half by both aircraft. Here, we constrain the necessary speed change #» s c that moves #» v r out of the forbidden zone to be entirely attributed to UAS A.

A B authorized half-plane # » v A # » v B #» v r #» s c # » v A
#» v r = # » v A -# » v B lies
We can then define a half-plane bounded by the line going through point A + # » v A + #» s c parallel to the corresponding edge of the cone on the side determined by #» s c : if the new speed of the UAS # » v A is chosen such that its endpoint lies in this halfplane, then the new relative speed #» v r is out of the forbidden speed zone, as illustrated on figure 2.

However, this minimizing strategy generally results in changing the norm of the speed of the UAS, which is not realistic considering the average performances of commercial UASs. We show in section II-B how to take this additional constraint into account in our geometrical model.

B. Constant Speed

We consider here that a UAS can only change its heading and not its speed norm. In reality, a UAS can change its speed norm, but the change rate is very limited and almost negligible with respect to the speed of surrounding aircraft (typically airliners).

To ensure that the norm of the speed of the UAS remains constant throughout the conflict resolution process, we must

A B (A -# » v B ) # » v A # » v B #» v r # » v A
Fig. 3. Model with constant speed norm for the UAS: the end of #» v r must lie on the circle but outside of the forbidden zone (in light red) and in the allowed turn angle range of 30 • around # » v A , which leaves the tiny green safe zone as possible new heading for the UAS.

have:

# » v A = # » v A as #» v r = # » v A -# » v B
, we also have:

(x v r + x v B ) 2 + (y v r + y v B ) 2 = # » v A 2
which means that the possible endpoints of

#» v r belongs to a circle of radius # » v A centered at A-# » v B as shown on figure 3.

But

#» v r must also lies outside the forbidden zone defined in the previous section, which removes all angle ranges of the circle included in it (in red on the figure). The remaining angles must be further filtered by intersecting the allowed turn angle range θ corresponding to the performance of the UAS, i.e. an arc of ±θ around the current speed shown in blue on figure 3 (where θ = 30 • ).

The conflict-free heading change thus ranges over an arc of a circle (pictured in green on figure 3) which is the difference of the allowed (in blue) and forbidden (in red) arcs. Note that the resulting allowed headings may comprise two disjoint angle ranges, which is generally the case for facing aircraft for example (as shown on figure 4) and corresponds to the combinatorial decision to turn right or left to solve the future conflict.

C. Heading Change Strategies

As mentioned in the previous section, the conflict-free solution set is generally not a singleton, so a decision must be taken to choose the "best" angle within possibly discontinuous angle ranges.

1) Closest: The new speed for the UAS can be aggressively chosen with the available angle that minimizes the heading change, which will result in optimizing the "maneuvers quantity" (a measure defined in section III-B). This strategy will be referred to as Closest in the remaining sections and is pictured in figure 5.

A B (A -# » v B ) # » v A # » v B
Fig. 4. Disjoint angle ranges (in green) for facing aircraft: the conflict can be solved either by turning left or right.

As it is expected that the relative performances of airliners and UASs will incur many unsolvable scenarios, for which losses of separation will occur, it may be wiser to choose an angle that leaves some leeway for the next resolution steps, i.e. which does not saturate the separation constraint.

2) Safest: To maximize the expectation of escaping a future conflict in case of further maneuvering of aircraft B (without any knowledge of its intended trajectory), the new heading could instead be chosen as the median value of the largest allowed angle range (in case of several disjoint angle ranges), which will optimize the robustness of the maneuver.

But the largest angle range must be selected before filtering the conflict-free angle ranges with the allowed turn angle (of ±θ, cf. section II-B), because these two kinds of constraints does not have the same significance: the new heading should be as far as possible from aircraft constraints (therefore chosen as the median value of the largest range) but can be arbitrarily close to the allowed turn angle constraint (e.g. the best maneuver may be to turn at the maximal possible angle). Once the corresponding angle has been computed with aircraft constraints only, the allowed turn angle constraint can be applied: if the median angle is no longer available, the closest allowed heading is chosen instead. This strategy will be referred to as Safest in the following and its result is shown in figure 5.

However, the overestimate of the conflicting situation inherent to the Safest strategy can lead to unacceptable amounts of heading changes and totally disregards the UAS mission and flight plan.

3) Hybrid: Actually, the two strategies can be combined to benefit from their advantages while getting rid of their drawbacks: the Closest strategy could be favored while there is still time and space to solve potential future conflicts, resorting to the Safest one before the situation becomes dire. So we define the Hybrid strategy parameterized by γ ∈ R ≥0 for a set of (adverse) aircraft positions F at a given time step: So the Closest strategy is chosen if the minimal distance between the UAS and all other aircraft is greater than γd, and the Safest one otherwise.

Hybrid = Closest if min B∈F # » AB > γd Safest otherwise A # » v A # » v target # » v Safest # » v Closest
Therefore, this parameter allows to continuously slide from a pure Closest strategy for γ = 0 to a pure Safest strategy for γ = +∞. The influence of parameter γ, which should a priori be chosen greater than 1 to switch to the Safest strategy when there is still enough leeway to escape the conflict, is discussed in section III-B.

4) Fallback Strategy: If the permitted heading range is empty, it means that no turning angle can guarantee the separation distance d for the next τ minutes. However, even if a conflict occurs, a brief one at a distance just below the target separation distance is preferable to a lasting close encounter.

To take into account this criterion while choosing a maneuver whenever a loss of separation is inevitable, our algorithm resorts to a fallback strategy consisting in reducing the separation distance until a solution can be found anew. Our current implementation uses a linear scheme with a constant step (0.1 NM in our experiments) and returns the corresponding maneuver as soon as a solution is foundthough a dichotomous scheme could be used to speed up the search if more precision is deemed relevant.

The resulting maneuver therefore attempts to minimize the conflict, i.e. to maximize the minimal distance between the aircraft and the UAS, hopefully achieving a minimization of the conflict duration as well.

D. Multiple Aircraft

In real traffic situations, conflicts involving more than two aircraft often occur and cannot be solved pairwise in sequence, but must be globally handled. Fortunately, our geometrical model can be easily extended to take into account any number of aircraft.

For each aircraft, a forbidden angle range is determined exactly as in the pairwise case detailed in section II-B and removed from the possible turning range of the UAS. The resulting allowed headings may then comprise several (possibly more than two) disjoint angle ranges, because each new

A B (A -# » v B ) # » v A # » v B C (A -# » v C ) # » v C
Fig. 6. Several disjoint allowed heading ranges (in green) for UAS A in presence of multiple aircraft (B and C).

constraint may remove a new hole in the remaining angle ranges (as illustrated in figure 6).

III. EXPERIMENTS

Our geometrical model has been implemented and thoroughly tested for the single aircraft case on recorded real TCA traffic enhanced by injecting various conflicting UASs scenarios. The following section described the data and scenarios generation, and the next one reports the results and analysis of our experiments.

A. Test Description

Experiments were conducted on real traffic data, recorded on 2013/09/14 in the south-west of the French airspace. This gives a much more realistic picture of the types of trajectories that UASs may encounter when flying in real traffic. This is all the more necessary because UASs are essentially meant to fly in the lower airspace where most of the commercial traffic is climbing or descending.

UASs trajectories are built from recorded tracks of aircraft with similar performances for two UAS types (flying at 80 kn and 160 kn respectively), three vertical profiles (leveled flight, climb and descent) and two horizontal profiles (constant heading and circle around a fixed point), for a total of 6 different patterns for each UAS type.

We filtered trajectories that had at least six minutes of flight under FL195, leaving us with a set of 475 aircraft for a total of 2850 simulation scenarios. For each aircraft and each UAS pattern defining a scenario, we then isolate a plot (x AC , y AC , FL AC ) (respectively (x UAS , y UAS , FL UAS )) in the trajectory of the aircraft (respectively the UAS) that has at least 3 minutes of past positions. These points are used to create a conflict for this scenario, i.e. we adjust the UAS trajectory in space and time dimensions so that x AC = x UAS , y AC = y UAS and FL AC = FL UAS . If no maneuver were issued, there would be a collision between the aircraft and the UAS.

A fast-time simulator enables to play the trajectories (both recorded and built) and to modify them by sending maneuver messages consisting of a heading change and a turn rate. Those messages are send to UASs only, other aircraft are left unchanged.

In the experiments, we chose τ = 5 minutes and tried to achieve 3 NM of separation when possible. We do not take the vertical separation into account because we want to test the efficiency of the detect and avoid process in the horizontal plane only. Further research will be conducted in the vertical dimension. We chose to detect and avoid conflicts every step seconds (step = 10 s when not specified otherwise).

B. Results with Various Speed and Turn Rate Ranges

Fast-time simulations were run for all the 2850 scenarios on a 3.4 GHz Intel Xeon workstation, each run being executed within one tenth of a second. We measured various indicators of the efficiency of the conflict avoidance maneuvers and counted the occurrences of "close distance events", which we set up to be the simulations were the distance between the UAS and the aircraft went under 1 NM. In the following, we will call an airprox (aircraft proximity) this kind of event.

Figure 7 depicts the number of airprox events as a function of the turn rate for various combinations of speed and resolution method. The speed of the UAS is clearly decisive to avoid conflicts efficiently, as the 160 kn UAS experiences less than 10 airproxes whichever the turn rate and the strategy, whereas for slow UASs, escaping conflicts becomes easier with a greater turn rate. For the latter, the strategy is also influential and, as expected, the Safest and Hybrid ones provide better separation than the Closest strategy.

The histograms on figure 8 show the distribution of scenarios with respect to the closest distance of approach during the simulation for both Closest and Safest strategies, with a UAS flying at 160 kn and maneuvering at 3 • /s. We observe a peak at 3 NM (i.e. the target separation distance) for the Closest strategy, which attests the correctness of the algorithm. Some conflicts noticeably remain in the left part of the curve, which the Safest strategy (red curve) helps to solve, although it yields larger separation distances, so that the UAS loses more time and fuel during maneuvers. The same simulations were carried out with the Hybrid strategy for two values of parameter γ, as depicted on figure 9, exhibiting the expected trade-off between maintaining separation and limiting deviations from the original trajectories. With γ = 1.4, the results are similar to the Closest strategy, however with significantly fewer losses of separation less than 3 NM. When γ increases, the 3 NM peak is shifted and spread towards the right part of the graph, indicating less airproxes but more costly maneuvers.

The number of airprox events is the primary indicator for the efficiency of the detect and avoid strategy. Yet, the deviation of the UAS from its trajectory and the number of maneuvers matter as well: they can be seen as the "cost" of maintaining separation or as a secondary objective to discriminate solutions with equivalent amounts of airproxes. Figure 10 shows the mean maneuver quantity, defined as the sum of the absolute values of all heading changes given to the UAS, for our three strategies (with two values of γ for the Hybrid one) as a function of the turn rate. As expected, the Safest strategy (in red) incurs the greatest maneuver quantity, whereas the Closest strategy is much more sparing. Once again, the Hybrid strategy occupies the middle ground between the two others, occasionnally providing better results than the Closest one. Not so intuitively, the turn rate seems to have little (predictable) effect on the maneuver quantity (contrarily to the number of airproxes, cf. figure 7), except for the 80 kn UAS with the Safest strategy for which the number and range of deviation increase with the ability to take steeper turns.

IV. CONCLUSION AND FUTURE WORK

In the context of the integration of UASs in TCAs, we have implemented a geometrical "detect and avoid" algorithm based on data collected from ADS-B-equipped surrounding aircraft. This algorithm attempts to solve air traffic conflicts between a UAS and airliners in initial climb and approach phases by changing the UAS heading only. The geometrical resolution is based upon the model described in [START_REF] Van Den Berg | Reciprocal n-body Collision Avoidance[END_REF], modified with additional constraints enabling the UAS to take the responsibility of the entire maneuver at constant speed.

Once the constraints induced by the aircraft have been computed, we propose two basic strategies, one of them (Closest) leading to fewer deviations from the initial trajectory but with significant amounts of remaining airproxes, the other one (Safest) providing more robust trajectories at the cost of many more maneuvers. To overcome the drawbacks of these strategies while retaining their advantages we present a tradeoff strategy (Hybrid), switching from one to the other whenever the distance between the UAS and one of the surrounding aircraft falls below a given threshold. These methods were intensively tested on 2850 realistic traffic scenarios issued from data recorded in a French TCA.

In most cases, the aircraft was avoided as planned, with a separation greater than the target. Nevertheless, in particularly difficult circumstances where the aircraft turns out to change its heading towards the escape trajectory computed during the previous time steps, airprox events could not be avoided, especially with the Closest strategy. However, the Safest strategy was found very efficient at limiting the number of these occurrences, but at a high cost in terms of "maneuver quantity". Eventually, the Hybrid strategy was found the most promising one, with a threshold ranging from one to three separation norms. Further experiments need to be performed to identify the Pareto front w.r.t. this threshold parameter.

One of the pitfalls of our method is that it only takes into account the current state, so that any further change in the aircraft state could break the resolution, especially with the Closest strategy. In order to improve the robustness of the maneuvers, we plan to try and anticipate better, both on the aircraft intentions and the UAS capabilities.

Knowing the past positions of the aircraft, it is possible to build a short-term predicted trajectory, based on the analysis of the derivatives of its speed and turn angle. For example, the beginning or the end of a turn, a climb or a descent could be inferred. Particular care would have to be taken during the calibration phase, especially when choosing the number of past states to consider: too much states would create some latency in the prediction, whereas too few states would yield unreliable ones.

If the aircraft trajectory could be predicted this way, then it becomes particularly interesting to anticipate several maneuvers for the UAS. This could be planned optimally with an A * or Dijkstra algorithm, using the geometrical algorithm at each step to prune the search tree or validate the existence of partial solutions, at the cost of a significantly longer computation. It could also be performed geometrically by an approximation of a few maneuvers aggregated into a single one.

For the hardest scenarios, we plan to check the interaction of the current algorithm with the TCAS of airliners and, if necessary, to implement a vertical strategy for the UAS that could comply with airliner TCAS maneuvers.

Fig. 1 .

 1 Fig. 1. Conflicting aircraft model: a conflict will occur within time τ if and only if the relative speed #» vr lies in the forbidden zone in red.

Fig. 2 .

 2 Fig. 2. Conflict resolution model: the necessary speed change to move #» vr out of the forbidden zone is endorsed by the UAS.

Fig. 5 .

 5 Fig. 5. Closest and Safest speeds, computed from # » v A , given the constraints in red and turn angle range in blue.

γ = 3 Fig. 7 .

 37 Fig. 7. Number of airprox events w.r.t. UAS turn rate. The graphs where γ is defined correspond to the Hybrid strategy.

Fig. 8 . 3 Fig. 9 .

 839 Fig. 8. Comparison of Closest (blue) and Safest (red) strategies: distribution of scenarios w.r.t. the closest distance of approach.

γ = 3 Fig. 10 .

 310 Fig. 10. Maneuver quantity w.r.t. UAS turn rate.

  in this cone, a conflict will occur in the future. If we draw a circle of size d τ tangent to the two previous lines, we obtain a new zone (in light red) bounded by the bold line in figure1. It is then straightforward to understand that a conflict will occur within time τ if and only if #» v r lies in this zone.